Multi-Scale Modeling and Simulation of Transport Processes in an Elastically Deformable Perforated Medium

In this paper, we derive an effective model for transport processes in periodically perforated elastic media, taking into account, e.g., cyclic elastic deformations as they occur in lung tissue due to respiratory movement. The underlying microscopic problem couples the deformation of the domain with a diffusion process within a mixed Lagrangian/Eulerian formulation. After a transformation of the diffusion problem onto the fixed domain, we use the formal method of two-scale asymptotic expansion to derive the upscaled model, which is nonlinearly coupled through effective coefficients. The effective model is implemented and validated using an application-inspired model problem. Numerical solutions for both, cell problems and macroscopic equations, are investigated and interpreted. We use simulations to qualitatively determine the effect of the deformation on the transport process.


Introduction
Heat and mass transport in materials characterized by a complex microgeometry are actively researched topic in various fields such as material sciences or geosciences as well as in the biological/medical field. Here, mathematical models can help to identify underlying mechanisms as well as complement or replace experiments by numerical simulations.

3
This is of particular importance in the medical field, where experiments and animal studies or clinical trials are expensive if at all feasible. The main question of this work is, how transport processes are influenced by a complex microstructure of the carrier medium and its deformation. Applications are especially found in the context of diseases whose clinical picture includes a dysfunction of the lung, such as COVID-19, pneumonia or sepsis. Experimental evidence already indicates that in addition to the microstructure formed at the cellular level by pneumocytes, the periodic deformation caused by inhalation and exhalation also plays an important role in transport processes such as those for nutrients or respiratory gases, see e.g. Huh et al. (2010).
Unfortunately, the heterogeneous microstructure of such problems can only be resolved to a very limited extent for small domains in numerical simulations. For problems involving multiple scales, one has to rely on mathematical tools that are capable of deriving approximate problems, which are easier to solve but still incorporate the information on the microgeometry. To this end, we derive in this paper an effective model for transport processes in an elastically deformable, periodically perforated medium, where we restrict ourselves to the case of a purely elastic solid phase that is not coupled with a fluid phase. This model is suitable for applications involving low-density fluids, for example, and it is also a first step for dealing with more general problems. We approach the problem by formulating first the microscopic problem on the perforated domain in a mixed Lagrangian/Eulerian framework. The deformation process is described by a linear elasticity equation on the microscopic reference domain as opposed to the transport problem, which is initially formulated on the current, moving microscopic domain. In this form, the problem is not accessible for standard upscaling tools. Therefore, we use a transformation, which is defined by means of the solution to the elasticity problem, to reformulate the transport problem on the reference domain for the elasticity equation. Then we can exploit the formal method of two-scale asymptotic expansion in order to obtain a macroscopic model describing transport processes and deformation for a periodically perforated domain. The resulting model is formulated on a homogeneous domain and involves effective coefficient functions. These effective quantities are computed by means of solutions to so-called cell problems which are formulated on the reference cell and carry the information on the microscopic geometry. Additionally, we develop and study numerical methods for the derived effective model. A lot of work has been done concerning the upscaling of fluid-structure-interaction-models, pure diffusion or elasticity, or (linear) coupled models for heat transport and deformation. In contrast, there are hardly any papers considering modeling and numerical simulations for an elasticity-diffusion problem, where the diffusion process is formulated within the Eulerian framework and transformed into a common framework, leading to coupling through nonlinear coefficients. We discretize both, the cell problems as well as the macroscopic problem using the finite element method and the implicit Euler method for discretization in time. Afterwards, we use a application-inspired model problem with oscillating Dirichlet boundary conditions in the deformation for numerical convergence tests to validate the code. For the model problem, we investigate and interpret the numerical results. Finally, again within the scope of the model problem, we give a quantitative answer to the question how the deformation affects the transport of a diffusing substance. To this end, we perform a parameter study for key parameters involved in the deformation of the domain and analyze the sensitivity of the effective model with respect to these parameters. This is accompanied by the visualization of the homogenized coefficient functions.
Homogenization techniques for elastic heterogeneous media (modelling, e.g., deformable composite materials) have been extensively studied in Oleinik et al. (1992). Transport processes in porous media involving an elastically deformable solid phase interacting with a fluid occupying the pore space where firstly modelled under the assumption of small deformations. Hence the equations for the fluid flow as well as the transport equation were formulated with respect to Lagrangian coordinates. For such problems, a rigorous homogenization can be found in Clopeau et al. (2001); Gilbert and Mikelić (2000) for the pure fluid-structure interaction without transport. See also references therein for similar results, especially using formal upscaling. In Clopeau et al. (2001) in the limit the famous Biotlaw is obtained, see (Biot 1956(Biot , 1962. The rigorous homogenization for the transport of biochemical substances within a deformable porous medium and their interaction with the mechanical properties of the elastic solid phase was performed in Jäger et al. (2011). A similar framework involving a deformable porous medium coupled to heat transport is upscaled in Brun et al. (2018). Furthermore, the authors of Allaire et al. (2017) consider the transport of electrolytes through the pore space of a deformable, charged porous medium and rigorously derive an effective model. More recently, models involving a mixed Lagrangian/Eulerian framework were considered, where the equations for the fluid flow as well as the transport equation were formulated on the deformed domain, i.e., with respect to Eulerian coordinates. Here, we mention, e.g., Collis et al. (2017), where a macroscopic model for transport processes in porous media, involving fluid-structure-interaction and growth has been derived by formal upscaling, see also (Brown et al. 2014) for the case without transport. However, the explicit representation of the effective coefficients was not specified. This was overcome for the pure fluid-structure interaction problem for porous media (without transport) in Miller and Penta (2020) for linear elastic and in Miller and Penta (2021) for hyperelastic multi-composite media. Further results related to our investigations deal with processes in porous media with an evolving microstructure. In Peter (2007Peter ( , 2009) reaction-diffusion models were homogenized assuming a given evolution of the microstructure. For the homogenization of a model of thermoelasticity, assuming again a given evolution of the microstructure due to phase transformation, see (Eden and Muntean 2017). In the context of crystal precipitation and dissolution, where the evolution of the microstructure is dictated by the local dissolution/precipitation rate, an upscaled model has been derived in van Noorden (2009), and the macroscopic model for a similar micro-model was analyzed in Schulz et al. (2017) for specific geometries leading to a dependence of the effective diffusion on the porosity. We emphasize that rigorous homogenization results of problems including a free boundary are quite rare, see for example (Gahn and Pop 2023) for a microstructure including spherical grains, where the radii depend on the solute concentration at the surface.
The paper is organized as follows: In Sect. 2, we introduce the microscopic problem and rewrite it in non-dimensionalized form on a fixed, heterogeneous reference domain. This is the starting point for the derivation of the effective model, which we obtain by means of a two-scale asymptotic expansion. In Sect. 3, we numerically compute, evaluate and interpret the solutions to the macroscopic problem. Before we finish with some final remarks and an outlook, we use computer experiments in Sect. 4 to quantitatively investigate how the cyclic deformation of the domain affects the transport process in the macroscopic model.

The Mathematical Model
Let us consider a macroscopic domain Ω given by the hyper-rectangle Ω = (a, b) ⊂ ℝ n , n ≥ 2 (the physically relevant cases are obtained for n = 2, 3 ), with a, b ∈ ℤ n such that a i < b i for all i = 1, … , n . This macroscopic domain contains a subdomain Ω s which exhibits a periodic microstructure and can be interpreted as the reference configuration of the solid part of a porous medium. The small parameter with −1 ∈ ℕ describes the ratio between the length of one microscopic cell and the size of the whole domain Ω . We require −1 ∈ ℕ in order to avoid dissected microscopic cells near the boundary of Ω . The microscopic domain Ω s is constructed in the following way: Let Y = (0, 1) n be the n-dimensional open unit cube and Y s ⊊ Y be a connected domain in Y that intersects all faces of Y, i.e. we have Additionally, we require that opposite faces of Y s are equal, i.e. it holds that By Γ ∶= int( Y s ⧵ Y) we denote the internal boundary of Y s in Y. Then we define the subdomain Ω s as and the internal boundary of the heterogeneous subdomain in Ω We note that by our assumptions Ω s is connected. This allows us to make sense of elastic deformation and diffusion processes in Ω s for n ≥ 2 . For a sketch of the microscopic domain, see Fig. 1.

The Microscopic Problem
Due to elastic deformation, the physical domain changes with time. We denote the current solid domain at time t ∈ [0, T] , T > 0 , by Ω s (t) and the current internal interface by Γ (t) ∶= Ω s (t) ⧵ Ω(t) . The reference domain Ω s is given by the current solid domain at time t = 0 , i.e. Ω s = Ω s (0) , and, analogously, for the internal reference interface we have Γ = Γ (0) . In the following, quantities and operators defined on or associated with the current domain Ω s (t) are marked by a hat ⋅ whereas quantities and operators defined on or associated with the reference domain Ω s do not carry a distinct label. The outer boundary of Ω s is divided into a Dirichlet and a Neumann boundary for the elasticity problem, i.e. we have In similar fashion, the boundary of Ω is split into corresponding Dirichlet and Neumann parts for the elasticity problem ( Ω = Γ elast D ∪ Γ elast N ) and for the diffusion problem ( Ω = Γ diff D ∪ Γ diff N ). We assume that the displacement u and the concentration ĉ satisfy the system Here, s > 0 is the constant density of the solid phase, A is a constant fourth-order elasticity tensor and by f elast ∶ [0, T] × Ω → ℝ n we denote a body force acting on the solid phase. e(w) = 1 2 (∇w + (∇w) T ) denotes the symmetric gradient and on the Dirichlet boundary of the elasticity problem, the displacement is prescribed by a function h ∶ [0, T] × Γ elast {t} × Ω(t) → ℝ is a source/sink term and ĝ ∶ ⋃ t∈(0, T) {t} × Γ diff D → ℝ , ĉ 0 ∶ Ω → ℝ denote the Dirichlet boundary value and the initial value for the diffusion problem, respectively. By n and n , we denote the outer unit normals to the reference

3
domain Ω s and the current domain Ω s (t) , respectively. Further, v is the velocity field induced by the deformation of the domain.
The elasticity problem (2.1a)-(2.1d) is defined within the Lagrangian framework on the reference domain Ω s . In contrast, the natural setting for the diffusion problem (2.1e)-(2.1h) is the Eulerian framework, i.e. it is defined on the current deformed domain Ω s (t) . In the present context, the shape of Ω s (t) is solely dictated by the elasticity problem, whereas in general, other processes, e.g., growth, might also be conceivable as driver of the deformation of the domain. Thus, using the elastic deformation defined by the current deformed domain Ω s (t) is given by and its internal, microscopic boundary is given by The velocity field v in equation (2.1e) is given by As previously remarked, the advection term for the concentration is solely due to the elastic deformation of the domain. Additional contributions to the velocity field which take into account other transport processes are of course possible but not considered in the following.

Transformation of the Diffusion Equation to the Reference Domain
In order to make the elasticity-diffusion system accessible for standard homogenization techniques, we transform equations (2.1e)-(2.1h) to the reference domain Ω s by using the deformation S defined in (2.2). For given, sufficiently smooth, scalar-and vectorvalued functions on the current deformed space-time domain Q T we obtain their respective counterparts on the reference domain via Additionally, we define the deformation gradient and its determinant by The following well-known computation rules connect differential operators and integrals from Ω s and Ω s (t) via the transformation mapping x = S (t, x) . Similar transformations are used, e.g., in Peter (2007).
Accordingly, the equations governing the diffusion process on the reference domain are given by where we use the definitions We note that the advective term in Eq. (2.1e) (on the current deformed domain) now vanishes. The effect of the deformation of the domain is encoded in the coefficients J and D which depend on the displacement u in a nonlinear manner.

The Non-dimensional Model
We analyse system (2.1a)-(2.1e), (2.4a)-(2.4d) further by introducing dimensionless quantities. The relationship between these quantities gives an idea of the qualitative behaviour of the system and also might influence significantly the homogenization result. As stated in the introduction, we have in mind the experimental setup from (Huh et al. 2010), where the crosssectional size of the microchannels is 4 ⋅ 10 −4 m × 7 ⋅ 10 −5 m . To mimic physiological breathing motion, the domain was cyclically stretched on two opposing boundaries by approximately 10% to reach a maximal length of 4.4 ⋅ 10 −4 m . In order to maintain applicability of the model to a broader variety of settings, we choose as characteristic size of the domain x = 5 ⋅ 10 −3 and the characteristic cell size l = 5 ⋅ 10 −6 m . This results in the characteristic size of heterogeneities = l∕x = 10 −3 .
Next, let us replace the dependent and independent variables by scaled versions: Here, quantities marked by a dagger ⋅ † are dimensionless variables and they carry a characteristic reference value which is marked by a bar ⋅ . Note that, in the definition of D , the quantities J and F do not carry a unit and therefore it is sufficient to associate the non-dimensional diffusion tensor D † with a reference value D . Then, with the definition The introduction of the non-dimensional spatial variables implies that the equations are now posed on the non-dimensional domain Ω s † of unit size. Realistic values for the reference quantities according to the application we have in mind can be found in Table 1.
Using the relations in (2.3) and dropping the daggers for the sake of an easier notation we obtain According to Table 1, we compute (since = l∕x = 10 −3 ) Characteristic reaction rate of source/sink term and therefore obtain the non-dimensionalised elasticity-diffusion system on the non-dimensional reference domain:

Remark 1
The existence of a weak solution for the elasticity Eqs. (2.5a)-(2.5d) is quite standard. The crucial point is the existence for the reaction-diffusion equation (2.5e)-(2.5h), since the coefficients depend on the microscopic displacement u and may degenerate in space and time for vanishing J . To overcome this problem higher regularity for the displacement is necessary, to guarantee the positivity of J as well as the injectivity of u , at least locally in time. However, this local existence interval may depend on and could vanish for → 0 . Even for finite strong assumptions and compatibility conditions on the data (and especially the microscopic domain) are necessary to obtain enough regularity for the displacement. Hence, the rigorous analytical treatment of this problem is an open question.

Upscaling in the Reference Domain
Although the domain for the elasticity-diffusion system is now fixed in time, we still have to deal with its heterogeneous microstructure. This issue is addressed by application of the formal two-scale asymptotic expansion method with the aim of deriving effective equations on the spatially homogeneous macroscopic domain Ω . The effective parameters in these equations are computed by means of solutions to cell problems which account for the heterogeneous nature of the microscopic problem.

Upscaling of the Elasticity Problem
Following the standard approach of two-scale asymptotic expansion, we postulate for the microscopic displacement u and concentration c expansions of the form with functions u i (t, x, y) and c i (t, x, y) depending on the the macroscopic variable x ∈ Ω , and the microscopic variable y ∈ Y s and being Y-periodic with respect to y. For derivatives of such functions, we have by the chain rule and analogous rules hold for the divergence operator. We will see that the expansion for u induces a two-scale asymptotic expansion of type where Ψ and Ψ i , i = 1, 2, ..., can be scalar-or vector-valued, also for the coefficients J and D . We start with the derivation of a macroscopic equation for the displacement, and identify the zeroth order term u 0 and the first order term u 1 in the expansion (2.6) for the microscopic displacement u . We emphasize that these results are well known in the literature, even for rigorous homogenization, see for example (Oleinik et al. 1992). However, for the sake of completeness and since we need the explicit representations for u 1 , we give some details for the upscaling process. Hence, let us plug in the expansion for u from (2.6) into Eqs. (2.5a)-(2.5d). We obtain and for the boundary conditions where we denote by n Γ (y) the normal to the boundary Γ and by n Ω (x) the normal to the boundary Ω and as well as for the initial condition By comparing the coefficients in Eq. (2.8), we obtain from the terms of order −2 together with the boundary condition of order −1 the following cell problem for u 0 : Since A is positive on the space of symmetric matrices we obtain e y (u 0 ) = 0 and therefore, since Y s is connected, u 0 is a rigid displacement with respect to y. However, the only periodic rigid displacements are constant functions and therefore we set Since e y (u) = 0 , the term of order −1 in (2.8) yields the problem i.e. u 1 can be determined in terms of u and due to the linearity of the equation we have the representation (unique up to a constant depending on t and x, which we choose equal to zero by assuming mean value zero for u 1 ) where ij , for i, j = 1, ..., n are the solutions to the cell problems Here, e i is the i-th canonical unit basis vector in ℝ n and (a ⊗ b) ij = a i b j is the dyadic product of two vectors a, b ∈ ℝ n . Finally, the term of order 0 , together with boundary conditions, gives an equation for u 2 in terms of u and u 1 : Integration over Y s and integration by parts with respect to y gives the macroscopic problem: which is actually an equation for u . In fact, exploiting now the representation (2.10) of u 1 , we see that u solves the quasi-static problem where the effective elasticity tensor A * is given by means of the solutions to the cell problems:

Upscaling of the Diffusion Problem
For the upscaling of the transport problem (2.5e)-(2.5h), we first expand the coefficients J and D with respect to according to (2.7) by using the expansion of u and the Taylor expansion resp. Neumann series. We start with the expansion for S , which is given by and therefore, using also representation (2.10) for u 1 , we get A ijkl kr ls + e y ( rs ) kl dy for i, j, r, s = 1, ..., n.
with E n the unit matrix in ℝ n×n . This is also an expansion F (t, x) = F 0 (t, x, y) + O( ) of the form (2.7) with Since we use a formal upscaling approach, we can assume that F 0 (t, x, y) is invertible. However, we emphasize that for a rigorous homogenization process this would be a critical point. As the determinant and the inverse are nonlinear operations, we employ the following linearizations: For general matrices A, B ∈ ℝ n×n with A non-singular and > 0 small enough, one has with the Taylor expansion for the determinant and the Neumann series for the inverse Consequently, we obtain for the determinant J and the diffusion tensor D In the following we will see that for the macroscopic model for the transport equation it will be only necessary to consider in detail the leading order terms of the expansions for J , D , as the higher order terms will not play a role in the homogenized equation. Since J 0 is bounded from below by a positive constant, we obtain that D 0 is positive. Now we are equipped to plug the expansions for c (2.6), J and D into Eqs. (2.5e)-(2.5h). Analogously to the asymptotic expansion for the elasticity subproblem (2.5a)-(2.5d) before, we obtain (2.14a) J 0 (t, x, y) = det F 0 (t, x, y) , together with the boundary conditions and and as well as the initial condition In the same spirit as for the elasticity problem, we obtain for order −2 from equation (2.15) the problem The positivity of D 0 implies that Therefore the problem for c of order −1 reduces to Integrating again over Y s we obtain the compatibility condition Let us define the effective coefficients J * and D * via By using representation (2.16) for c 1 we obtain after an elemental calculation that c solves the following macroscopic problem: x, y) kj + y k j (t, x, y) dy for i, j = 1, ..., n.
(2.19a) t (J * c) − ∇ ⋅ (D * ∇c) = J * f diff in (0, T) × Ω, 1 3 Note that, even under the assumption that D in the initial problem (2.1e) is constant, we obtain now time-and space-dependent homogenized coefficients J * , D * for the diffusion problem, since they carry the information on the deformation of the domain after the transformation into the reference domain.

The Upscaled Equations on the Reference Domain
We now summarize the upscaled equations which we have derived previously. They consist of coupled, effective problems for the displacement and the concentration, and are complemented by cells problems for the elasticity problem and for the diffusion problem, respectively. Thus, the macroscopic problem is given by The elasticity cell solutions ij , for i, j = 1, ..., n are obtained by solving and the solutions to the diffusion cell problems i , i = 1, ..., n are obtained from From the cell problems, the following effective quantities are computed: x, y)dy = 0.

Numerical Simulation of the Upscaled Model
In the following, we focus on the simulation of the upscaled model summarized in Sect. 2.5. Due to the multi-scale character of the problem, computations are performed on two distinct domains: Ω for the macroscopic elasticity-diffusion system and Y s for the cell problems for elasticity and diffusion (see Fig. 2). We avoid having to solve a nonlinear problem by organizing the scheme according to the special structure of the system. In each time step, we solve first the quasi-static macroscopic elasticity problem independently from the diffusion problem due to the coupling being only one-sided. Next, the effective coefficients for the diffusion equation, J * and D * , are computed from the solution of the macroscopic elasticity problem and the diffusion cell solutions. Finally, with the new homogenized coefficients, the macroscopic solution in the current time step is computed.
The scheme is also organized in a way that takes advantage of the specific properties of the cell problems. To this end, the treatment of the elasticity cell problems differs from the treatment of the diffusion cell problems due to the latter depending on time and location of the cell within the macroscopic domain. It is in fact sufficient to only solve the elasticity cell problems once for each i, j = 1, 2 at the beginning of the simulation. The resulting effective elasticity tensor A * is constant in time and space during the subsequent simulation. In contrast, the diffusion cell problems depend on the solution of the macroscopic elasticity problem through the coefficient D 0 (see (2.14b) and (2.17a)).
A ijkl kr ls + e y ( rs (y)) kl dy for i, j, r, s = 1, ..., n, x, y) dy for i, j = 1, ..., n.  Fig. 2 Non-dimensionalized domains Ω = (− 1 2 , 1 2 ) 2 for the effective, macroscopic elasticity-diffusion problem and Y s = ((a, 1 − a) × (0, 1)) ∪ ((0, 1) × (b, 1 − b)) , with a, b ∈ (0, 0.5) for the elasticity and diffusion cell problems Therefore, different diffusion cell problems are solved for each quadrature point in the macroscopic grid and for each time step as the effective quantities J * and D * change with time and space in the macroscopic diffusion equation. Let us underline again that this effect is introduced when the initial diffusion problem on the moving domain is transformed into an equation on a reference domain in Sect. 2.2 and independent of the initial diffusion tensor D (cf. (2.1e)) being time-and/or space-dependent or not.
The coupled two-scale system was implemented within the C ++ finite element library deal.II (Arndt et al. 2020). All computations are performed on uniformly refined grids based on quadrilateral elements in two space dimensions. Spatial discretization for both, the macroscopic equations and the cell problems, is achieved using (bi-)linear Lagrange finite elements. For temporal discretization of the effective, macroscopic diffusion Eqs. (2.19a)-(2.19d), the Crank-Nicolson method is employed. The process of solving all the diffusion cell problems consumes most of the computation time during the simulation. Optimization of this step in terms of parallelization or some adaptive scheme is an interesting and promising task to save computation time but not within the scope of this paper.

The Model Problem
We fix now data for further investigations of the upscaled problem from Sect. 2.5. Let Ω = − 1 2 , 1 2 2 be a quadratic domain and Y s = 1 3 , 2 3 × (0, 1) ∪ (0, 1) × 1 3 , 2 3 a cross-shaped domain (cf. Fig. 2). For Γ elast D , we chose the lateral parts of Ω and set Γ elast N ∶= Ω⧵Γ elast D . We do not consider a body force acting on the material and therefore set f elast ≡ 0 . The deformation of the domain is solely induced by the time-periodic Dirichlet boundary condition with a = 0.25 , the maximal displacement at the boundary, and f = 1 , the frequency. For an intuition about this boundary condition, see Fig. 7b and c. We further assume that the material is isotropic, i.e. its elastic properties can be described by two Lamé-constants , and the entries of the elasticity tensor are given by For the non-dimensionalised elasticity tensor in the elasticity cell problems (2.11a), we set , = 1 . In Table 2, we list the non-zero entries of A opposed to their counterparts in the effective elasticity tensor A * and see that the isotropy of A is not conserved during the upscaling process, i.e. we cannot find * , * ∈ ℝ such that A * has a representation as in (3.2). For the diffusion problem, we set Γ diff D to be the upper part of the boundary of Ω and Γ diff N = Ω⧵Γ diff D . We do not consider a sink/source, i.e. we set f diff ≡ 0 and prescribe a constant concentration g ≡ 1 at the Dirichlet boundary as well as no-flux conditions on the Neumann boundary together with initial condition c 0 ≡ 0 . The non-dimensionalized and constant initial diffusion tensor in (2.14b) is set to be As mentioned before, the effective quantity D * is not constant in time or space in general, but depends on the displacement via the deformation gradient. Nonetheless, for t = 0 , it is constant across the whole domain Ω since it is initially in the non-deformed state. Therefore, the difference between D and D * (t = 0, ⋅) is only due to the heterogeneity of the domain and not due to deformation, when we compare the non-zero entries of these quantities in Table 2.
For a more comprehensive study of the properties of the effective coefficients in dependence of the micro-geometry, see Fig. 3. There, we have computed and visualized the nonzero entries of A * and of D * for t = 0 as a function of the widths of the intersecting bars of the cross Y s . One can directly deduce that increasing the widths increases also the entries of the effective coefficients. In particular, it can be observed that the effective coefficients approach the values of the coefficients A,D of the underlying microscopic problem (cf. Table 2) when the widths of the bars approach 1, i.e. when the cross is close to occupying the whole unit cube [0, 1] 2 . We note that, in general, although two distinct crosses might occupy the same volume fraction of Y, they do not lead to the same effective coefficients, as asymmetries are also reflected in A * and D * .

Numerical Convergence Studies
Before we proceed with a more detailed view on the solutions of the proposed model problem, let us first investigate the convergence of the finite element scheme. Due to the two-scale character of the problem, there are two distinct computational domains, representing Ω and Y s , respectively. In the following, we restrict ourselves to the analysis of the macroscopic solutions u and c with respect to the discretization parameter h of the grid for Ω , where h is the diameter of the largest quadrilateral in the triangulation. For the simulations, we use the model problem presented in the previous Sect. 3.1. To quantify its   convergence properties, we compute the error between solutions at the fixed time t = 1.5 on subsequently refined grids, using an uniform refinement strategy to obtain finer meshes. Additionally, we compute the estimated order of convergence (EOC) according to the following formulas: where i = 1, 2, ... indicates the refinement cycle, ◻ ∈ {u, c}(t = 1.5, ⋅) , and * ∈ {L 2 (Ω), H 1 (Ω)} . The so-obtained data is listed in the upper part of Table 3. Note that the EOCs for the most part lie below the convergence orders known from theory for the implemented finite element method of linear Lagrangian elements (order two for L 2 -norm and order one for H 1 -norm). This suboptimal behavior might be attributed to the mixed boundary conditions in our model problem, which are known to be a potential source of decreased regularity of solutions. In fact, we are able to recover the expected convergence rates, if we simulate instead a related problem with pure Dirichlet boundary condition and matching initial condition for the diffusion problem. Let Γ elast with, as above, a = 0.25 , the maximal displacement at the boundary and f = 1 , the frequency. This boundary condition for the elasticity problem results in a parabola-shaped extension of the domain at the lateral parts of the boundary while the upper and lower parts of the boundary is clamped, see Fig. 4. Analogously to the study of the problem with mixed boundary conditions, we gather convergence data for simulations of this problem in the lower part of Table 3. Comparison clearly indicates that the loss of convergence speed can be attributed to the mixed boundary conditions. However, the actual error between subsequent solutions appears to be smaller in the case of mixed boundary conditions, at least for courser meshes, when compared to the case of pure Dirichlet boundary conditions, see Fig. 5.

The Solution of the Elasticity Subproblem
Let us now give a more detailed analysis of the numerical solution to the effective elasticity subproblem using the aforementioned data. As stated earlier, it suffices to solve the elasticity cell problems (2.11) once at the beginning of the simulation since the coefficient A is constant. Note that the symmetry properties of A imply the symmetry 12 = 21 . The resulting cell solutions ij , i, j = 1, 2 , contain the information needed for the computation of the effective elasticity tensor A * in (2.13). For a visualization of the elasticity cell solution, see Fig. 6. The non-zero entries of A * are given in Table 2. The cell solutions are also contained in the representation of the first order term u 1 , see (2.10), in the expansion (2.6) for u . In fact, we are now equipped to give an approximation of u in terms of the asymptotic expansion up to terms of order O( 2 ) and higher. To this end, let us investigate the mth component, m = 1, 2 , of u = (u 1 , u 2 ) T . For a fixed point x ∈ Ω s , we write x = x * + y   with x * = k , k ∈ ℤ 2 , y ∈ Y s , and obtain by using Taylor expansion and the symmetry  Table 3. We compare convergence properties of the model problem with mixed boundary conditions described in Sect. 3.1 and the model problem with pure Dirichlet boundary condition described in Sect. 3.2 with respect to the discretization parameter h of the grid for the macroscopic domain Ω In Fig. 7a, we give a visualization of the -term in the above computation for y ∈ Y s at t = 0.5 , i.e. when the material is maximally extended, to illustrate how the deformation acts on individual cells. Additionally, the effective displacement u(t = 0.5, ⋅) is displayed in Fig. 7b. There, the displacement u , defined on Ω , is represented by arrows with size proportional to their magnitude. We can also use u to visualize the deformed domain Ω 0 (t) ∶= {x ∈ ℝ 2 |x = x + u(t, x), x ∈ Ω} as can be seen in Fig. 7c for time t = 0.5 . This gives us also a clear intuition for the effect of the boundary conditions of the elasticity problem stated earlier: Within one unit of time, the domain is stretched in lateral direction to attain the shape depicted in 7c. Then, the displacement at the boundary is relaxed until the domain is again in the reference configuration.

The Solution of the Diffusion Subproblem
Finally, we investigate the solution of the diffusion equation for the model problem from Sect. 3.1. As mentioned before, the diffusion problem depends on the elasticity problem through its coefficients, but not vice versa. Therefore, in each time step, we solve first for the displacement, compute from this the effective coefficients D * , J * for the diffusion problem and then solve the diffusion problem. The effective coefficients for the diffusion problem are defined by means of the diffusion cell solutions i and the coefficient D 0 for the diffusion cell problems is given in terms of F 0 , i.e. the lowest order term in the expansion of the deformation gradient F . Concerning the numerical scheme, this amounts to dim × number of time steps × number of elements × number of quadrature points per element diffusion cell problems that have to be solved numerically during the simulation. For illustration, numerical approximations of the functions i (t = 0.5, x = (0, 0) T , y) , i = 1, 2 , y ∈ Y s , are plotted in Fig. 8. Unsurprisingly, the computation of the diffusion cell solutions takes up most of the computation time, but there is some potential for optimization: The diffusion cell problems are independent from each other, so their treatment can be parallelized. Additionally, adaptive schemes such as clustering "similar" cell problems and solving only one representative cell problem per cluster are conceivable and have been successfully applied in the context of multi-scale models, see e.g. Bastidas et al. (2020); Gärttner et al. (2020). Due to the transformation of the diffusion problem onto a reference domain in Sect. 2.2, the concentration c, defined on Ω , does not have a direct physical interpretation, as concentrations in the real world would be associated with points in the current deformed domain. Nonetheless, the concentration ĉ , defined on the current deformed domain, can be obtained from the concentration c using the following transformation and Ω 0 (t) ∶= {x ∈ ℝ n ,x = S 0 (t, x), x ∈ Ω}.  Fig. 9. There, one can also spot an interesting effect of the timeperiodic deformation of the domain on the accumulation of particles in the domain: When the domain is expanded, the area increases and the concentration shrinks. Consequently, a larger concentration difference develops between the constant Dirichlet boundary condition at the top of the domain and the adjacent interior of the domain. This leads to an increased flux of particles into the domain. When the displacement is relaxed, the concentration rises again. After some time (see Fig. 9 at t = 25 ), the concentration can even increase to values that are higher than the Dirichlet boundary condition ĉ = 1 . This would not be possible in the case of a non-deforming domain.

Sensitivity of Transport Processes with Respect to Changes in the Deformation
Finally, we analyze the system's sensitivity with respect to experimental parameters which we can modulate by manipulating the boundary conditions of the elasticity problem. We adopt again the setup from the model problem in Sect. 3.1, i.e. the deformation of the domain is driven by the time-periodic Dirichlet boundary condition h (cf. Eq. (3.1)), while a diffusing substance spreads into the domain, originating from the Dirichlet boundary condition of the diffusion problem at the upper part of the boundary of Ω . Since we are concerned with the effect of the deformation on transport processes, we monitor the quantity over the course of each of the following computer experiments. M(t) is the approximation of the mass or the number of particles that resides at time t in the domain Ω . Different experimental scenarios arise from varying the frequency f or the amplitude a in (3.1). By comparing M(t) for different simulations we gain insights on the sensitivity of the transport of the diffusing substance with respect to the deformation parameters.
The macroscopic concentration ĉ , visualized on the current deformed domain Ω 0 (t) during simulations of the model problem from Sect. 3.1 Fig. 10a shows how the mass of the diffusing substance accumulates in the domain for different frequencies of the lateral displacement. The curves indicate no significant qualitative difference at first glance. Only a closer look reveals slight deviations in the small mass fluctuations. The question may arise why the mass in the domain can temporarily decrease during the advanced stage of the simulations. The answer is already partly given in Sect. 3.4: when the domain is released back into the relaxed configuration, the concentration inside the domain may reach values that are greater than the prescribed Dirichlet boundary condition, consequently leading to an outflow of particles across that very boundary.
The system clearly exhibits a higher sensitivity with respect to the amplitude of the lateral displacement, as can be seen from Fig. 10b. There, the mass M(t) is plotted for varying amplitudes of the lateral displacement. The percentage in the legend encodes the maximal length of the domain Ω 0 (t) in lateral direction during each experiment in relation to the edge length of Ω . Positive percentages indicate cyclic extension and relaxation into the initial state whereas negative percentages indicate cyclic compression and subsequent relaxation. From the numerical experiments it is obvious that deforming the domain with increasing amplitude leads to the accumulation of more mass in the domain. This effect can be attributed to the increased area of the extended configuration: During expansion, the concentration shrinks, hence leading to a greater concentration difference between the constant Dirichlet boundary condition and the concentration inside the domain, ultimately resulting in an increased flux into the domain across the respective boundary. Compression of the domain leads to the opposite effect.
Visualization of the homogenized quantities J * and D * in Ω for a fixed time point when the physical domain Ω 0 (t) is in the maximally deformed configuration, sheds more light into the properties of the diffusion process during the course of the simulations. In Fig. 11a-e, J * and the components of D * are plotted for the case where a = 0.25 , corresponding to 50% extension from Fig. 10b. Unsurprisingly, the homogenized determinant J * , which keeps track of volume changes of the domain, is increased everywhere (note that J * for the relaxed configuration is constant over space with J * = |Y| = 5∕9 due to the heterogeneity of the underlying microscopic problem). From Fig. 11b and e we deduce that lateral stretching of the domain facilitates diffusion in -50% -40% -30% -20% -10% 0% +10% +20% +30% +40% +50% (b) amplitude dependence  Table 2) in the non-deformed configuration. Additionally, we note that anisotropy effects are introduced due to non-zero off-diagonal entries D * 21 and D * 21 . This can be attributed to the deformation of the domain, as we have seen (cf. Table 2) that the isotropy property of the initial diffusion coefficient D (3.3) transfers to

Conclusion
In this work, an effective model for transport processes in elastically deformable, heterogeneous media has been derived and investigated. Starting from a microscopic description of an elasticity-diffusion problem, which contained the underlying heterogeneity of the modelled process explicitly in a mixed Lagrangian/Eulerian framework, an upscaled model, formulated on a fixed, homogeneous, reference domain on the macroscale, has been derived. The underlying heterogeneity is addressed via cell problems for the deformation and the diffusion process. Although the initial problem is posed with constant coefficient functions, a transformation of the diffusion problem onto the reference domain introduces a diffusion coefficient for the cell problems which depends not only on the microscopic space-variable, but also on time and on the macroscopic space-variable. Solving individual cell problems and computing effective coefficients in each time step and for each quadrature point for the effective diffusion problem are the most time consuming tasks when simulating the effective model. There is some potential for optimization by computing cell solutions in parallel or by solving fewer, representative cell problems which are chosen by an adaptive algorithm.
Investigations associated to the numerical solutions of the upscaled problem such as convergence tests and interpretation of solutions within the context of realistic experimental setups, inspired by applications, are used to verify the implemented model and give insight on the qualitative and quantitative behavior of the model. The simulations emphasize, that the deformation can significantly influence the transport processes and therefore shows the necessity to consider e.g. the effect of respiratory movement when investigating transport processes on the cellular level in the lung or related problems.
The present model has to be understood as a starting point for further research in different directions. As a first generalization, it is reasonable to expand the model by including also an influence of the diffusing substance on the elastic properties of the solid, leading to a fully coupled system. In the medical context, this may be especially interesting when the effect of fibrosis on the functionality of the lung has to be considered. Furthermore, until now, the effects of the perforations of the microscopic domain are only considered in the sense that they introduce the heterogeneity of the problem. From a biological perspective, the space that is not occupied by cells can be interpreted as the extracellular matrix, which, as a first approximation, can be modelled by a fluid. Starting from this, a more comprehensive model could be derived, governing also the effects of fluid-structureinteraction in porous media. Here, due to the mixed Lagrangian/Eulerian framework, the effective model includes a nonlinear version of the Biot-equations for poro-elasticity. We remark that the homogenization of fluid-structure-interaction problems in porous media in Lagrangian/Lagrangian framework leads to linear Biot-equations where the effective elasticity coefficients are the same as the coefficients A * from (2.16) in our paper, see, e.g., (Jäger et al. 2011, (6.9)) or (Brun et al. 2018, (53)).
Finally, as we state in the introduction, the model is inspired by the lung and a related biomimetic microdevice that shares common features with the lung. Although the structures of these examples are based on thin layers (e.g. air-blood-barrier, consisting only of two cellular layers and the basal lamina), we chose to formulate the model for a domain that spreads in all space dimensions for the sake of simplicity.For the homogenization of a model governing an elasticity-diffusion process in a thin perforated elastic layer, an additional limit process has to be performed by letting the thickness of the layer tend to zero. There multi-scale methods for simultaneous homogenization and dimension reduction are needed. For the case of elastic perforated layers without transport such methods have been developed in , see also (Griso et al. 2020).