Homogenization of coupled flow and deformation in a porous material

In this paper we present a framework for computational homogenization of the fluid-solid interaction that pertains to the coupled deformation and flow of pore fluid in a fluid-saturated porous material. Large deformations are considered and the resulting problem is established in the material setting. In order to ensure a proper FE-mesh in the fluid domain of the RVE, we introduce a fictitious elastic solid in the pores; however, the adopted variational setting ensures that the fictitious material does not obscure the motion of the (physical) solid skeleton. For the subsequent numerical evaluation of the RVE-response, hyperelastic properties are assigned to the solid material, whereas the fluid motion is modeled as incompressible Stokes' flow. Variationally consistent homogenization of the standard first order is adopted. The homogenization is selective in the sense that the resulting macroscale (upscaled) porous media model reminds about the classical one for a quasi-static problem with displacements and pore pressure as the unknown macroscale fields. Hence, the (relative) fluid velocity, i.e. seepage, "lives" only on the subscale and is part of the set of unknown fields in the RVE-problem. As to boundary conditions on the RVE, a mixture of Dirichlet and weakly periodic conditions is adopted. In the numerical examples, special attention is given to an evaluation of the Biot coefficient that occurs in classical phenomenological models for porous media.


Introduction
We consider the problem of fluid flow through deformable porous materials. Porous materials are present in a vast number of natural as well as engineered structures. Examples of natural structures include biological tissue and aquifers, while examples of engineered structures are foams and textiles. The microstructure of porous materials is generally very complex with characteristics at a lengthscale much smaller than the scale of the application; hence, it is computationally not feasible to solve the fully resolved problem. Consequently, macroscopic phenomenological material models, based on a priori homogenization, are commonly used. Starting from Biot [1], a large number of so-called "porous media theories" of various complexity have been developed. An important class is the so-called Porous Media Theory (PMT) [2]. Although computationally efficient, PMT models are less capable of representing the intrinsic physical properties of the material.
As a viable alternative, computational homogenization [3] may be used, whereby the material response is evaluated using a Representative Volume Element [4] (RVE) that contains a small subset of the fully resolved microstructure. To ensure that the response from the RVE is physically feasible, boundary conditions satisfying the Hill-Mandel condition [5,6] must be chosen.
This paper concerns the homogenization of a fluid-filled (saturated) porous material. The pore system is assumed to be open, and we restrict the analysis to two phases; one solid and one fluid phase. Since we consider 3D microstructures, the solid skeleton is assumed to be contiguous in order to be able to sustain loads. In particular, we account for the interaction between the deformable solid and fluid, which represents a Fluid-Structure-Interaction (FSI) problem. The homogenization of flow in deformable porous media is also addressed by Iliev et al. [7] where the special case of flow through a deformable channel using asymptotic expansion.
We restrict to the case of laminar and incompressible flow of the pore fluid. For the FSI-problem, we use a monolithic approach with a conforming interface mesh. In order to maintain a proper mesh in the fluid domain during deformation, we introduce a fictitious elastic material in such a way that the mesh in the fluid domain follows the deformation of the solid. Measures are taken in order to ensure that the presence of the fictitious elastic material does not contribute to the overall stiffness of the RVE.
The paper is outlined as follows: In Section 2, the fine scale FSI-problem is established (in the material format). In Section 3, we introduce the homogenization scheme in the setting of the Variational Multiscale method and identify the macroscale problem. In Section 4, we discuss the RVE-problem. In Section 5, two numerical examples are presented. Finally, in Section 6 conclusions and future work are discussed.

Preliminaries
As a starting point, we consider the fully resolved domain in the material configuration depicted in Fig. 1, where the gray area represents the solid and the white area the fluid. We consider two phases, a solid phase, which is located in Ω s and a fluid phase located in Ω f . We also define the total domain and boundary as Ω = Ω f ∪Ω s and Γ = Γ f ∪Γ s . Γ f is the part of Ω f where fluid can enter and exit the domain and Γ s is the part of Ω s on the outer boundary. The interface between the solid and fluid phases is introduced as Γ int = ∂Ω s ∩ ∂Ω f . In order to allow for subsequent boundary conditions on both phases, we split Γ f and Γ s into a Dirichlet part and a Neumann part. Thus, for the solid boundary, we have Γ s = Γ s u ∪ Γ s t where the displacement is prescribed on Γ s u (Dirichlet) and the traction is prescribed on Γ s t (Neumann). Likewise, on the fluid part of the boundary, we have where the velocity is prescribed on Γ f v (Dirichlet) and the pressure on Γ f p (Neumann). N is the outward pointing normal to Γ and N s is the normal on Γ int , pointing outwards from Ω s . Note that, due to the forthcoming discretization of the domain and the interaction between the solid and the fluid, the computational mesh in Ω f needs to follow the deformation of the solid in order to avoid the possibility of a distorted mesh. In a staggered approach to FSI problems, a mesh of poor quality in the fluid part of the domain should be updated after the solid deformation is computed in each iteration, using e.g. Laplace smoothing or by remeshing. However, the aim here is to establish a monolithic method and we propose to incorporate the update in the set of equations. To that end, a fictitious elastic material is introduced in Ω f . The material is then attached to the internal boundary Γ int . As the internal boundary moves, the fictitious elastic material is deformed, and the nodes inside Ω f follow the deformation of the solid. To clarify, consider Fig. 2 where Fig. 2(a) is the undeformed mesh of both phases. In Fig. 2(b), the solid is deformed and no means for correcting the mesh in Ω f is taken. Thus, the positions of nodes inside Ω f are fixed which can cause the elements along Γ int inside Ω f to become ill-shaped. However, in Fig. 2(c), the deformation of the fictitious elastic material ensures that the nodes in the fluid follow the fictitious elastic solid, thus producing a proper mesh.
In what follows we assume that the motion of both solid and fluid is quasi-static in the sense that inertia terms are left out. However, the motion of the fictitious elastic material in the pore space is still time-dependent, i.e. d t u f ̸ = 0, until the deformation becomes (possibly) stationary.
We now introduce the displacement u s on Ω s , total velocity v on Ω f and pressure p on Ω f , all defined in the material configuration. We also introduce the displacement u f on Ω f pertinent to the fictitious elastic material. The total velocity v of the fluid is given as where d t u f is the velocity of the fictitious elastic material, and w is the relative fluid-solid velocity. Obviously, d t u f tends to zero as the quasi-static solution is approached.

Remark 1.
As u f is not a "placement" of fluid particles, it lacks physical interpretation.

FSI-problem-Canonical format
Here, we restrict the analysis to that of a laminar and incompressible flow without convective acceleration, i.e., Stokes' flow. The FSI problem on strong form in the material configuration is given as where P s is the first Piola-Kirchhoff stress of the solid, while P f = P v (v ⊗ ∇, F) − p J F −T is the first Piola-Kirchhoff stress of the fluid where P v is the viscous stress and p is the hydrostatic pressure. Furthermore, F = I + u f ⊗ ∇ is the deformation gradient of the fictitious elastic material, and J is the Jacobian thereof. Eqs. (2a) and (2b) are the momentum balance equations for the solid and fluid phase, respectively. Eq. (2c) is the mass balance equation under the assumption of incompressibility. Here, we also give the pertinent boundary conditions as where t s and t f define the solid and fluid tractions.t s andt f are the values assigned to the respective tractions,û is the prescribed displacement along Γ s u ,p is the prescribed pressure on Γ f p in the current configuration andŵ n is the prescribed fluid velocity out from Ω f .
On the internal boundaries, Γ int , the interaction conditions are given as and we note that the no slip condition in Eq. (4b) is equivalent to v = d t u.
In the subsequent derivation, we first view the FSI problem and the problem pertinent to the fictitious elastic material as two separate problems. We will then combine them such that the "mesh update apparatus" is incorporated in the solver. Finally, we will reduce the number of unknowns by manipulating the solution and test spaces. The problem defined in Eqs. (2)-(4) is given on weak form as: where the abstract operators in Eq. (5) are defined as The Lagrange multiplier t in Eq. (5) constitutes the reaction force due to the no-slip condition in Eq. (4b). The solution spaces are defined as and the test spaces are defined as

FSI-problem-Auxiliary problem and monolithic formulation
When a spatially fixed FE-mesh is used in Ω f for fluid motion, the mesh quality will gradually deteriorate due to large deformation in the solid and corresponding large motion of the boundary Γ int , as illustrated in Fig. 2(b). Some sort of mesh update is obviously needed, whereby the classical method is "mesh-smoothing". Another option, which is adopted in this paper, is to introduce an auxiliary problem in Ω f that represents the motion of a (fictitious) solid material, whose displacements along Γ int are prescribed to those of the (real) solid material in Ω s . The constitutive properties assigned to this fictitious material are conveniently chosen as elastic. The mesh may now be continually updated (in time) from the resulting fictitious displacement field. On strong form, this auxiliary problem is given as where u s is the displacement solution to Eq. (5), and N s : Ω s → Γ f is a nonlocal boundary operator defining the external flow boundary in terms of the surrounding solid placement. The corresponding weak problem for the fictitious elastic material becomes: and P X is the first Piola-Kirchhoff stress tensor pertaining to the fictitious elastic material. The solution space is defined as and the test space as So far, we can solve the two problems in an iterative manner. However, we can combine them and thus incorporate the apparatus for maintaining a proper mesh in the system of equations. In addition, we choose to reformulate the combined equations using the relative velocity w rather than the total velocity v. We thus proceed as follows: Due to the fact that T is complete, the condition in Eq. (5d) is fulfilled identically. Thus, without altering the solutions u s , u f and p, we can redefine the velocity in terms of v = d t u f + w, whereby the pertinent solution space for the relative velocity w becomes Here, we have restated the Dirichlet boundary condition for the flow by introducingŵ :=v − d t N s u s . In order to impose (9b) explicitly, we can now construct a global displacement field u, defined as and introduce the function space U as By construction, we note that u is continuous across Γ int . The resulting problem can be posed as follows: Find As a final step, in order to eliminate the interfacial traction, t, we also merge the test spaces U s,0 and V 0 as follows: As a result, we can add Eqs. (17a) and (17b), whereby the terms with the traction t cancel. We thus pose the problem in its final weak form as: where we have introduced the new notation W * ≡ U f,0 and To conclude, we now test both the equilibrium balance for the solid and fluid phases with the same test function, δu, which exists on the whole Ω . Thus, force equilibrium between the two will be fulfilled at all times. Furthermore, due to the fact that δw = 0 on Γ int , the fictitious elastic material will not contribute to the overall stiffness, i.e. the coupling between the deformation in the solid and fluid goes only one way.

Selective VMS
In the spirit of the Variational Multiscale method, cf. [8,9], we can make the assumption that the fields u and p can be split uniquely into a smooth (slowly varying) macroscale part and a non-smooth (rapidly fluctuating) subscale part. We thus introduce the hierarchical decomposition P = P M ⊕P µ and U = U M ⊕U µ where P M , U M contain the slowly varying (macroscale) parts, whereas P µ , U µ contain the rapidly fluctuating (subscale) parts. Hence, p = p M + p µ and u = u M + u µ , where p M ∈ P M , p µ ∈ P µ , u M ∈ U M and u µ ∈ U µ . No such decomposition is made of W; hence, it is assumed that w ∈ W "lives" only on the subscale. Now, by solving local problems for p µ , u µ and w = w µ in terms of p M and u M (as the driving fields), we may formally obtain the functional dependence The curly brackets {•} in Eq. (21) indicate implicit dependence on the arguments. Furthermore, we split the pertinent test functions in Eq. (19) as δp = δp M + δp µ , and δu = δu M + δu µ whereby we are able to split the problem into a macroscale and a subscale problem. Consequently, we get for the subscale problem: and for the macroscale problem: Find (u M , p M ) ∈ U M × P M such that 1 We now reformulate Eq. (23b) upon integrating by parts To obtain the last term, we used the identity J F −T : Moreover, we used that w = 0 on Γ int to conclude that the surface integral does not get any contribution from Γ int . We now restate the macroscale problem in Eq. (23) as: where we introduced the additional forms and Γ f w is the part of Γ where w is prescribed.

Homogenization
For the subsequent homogenization, we introduce the intrinsic averaging operators ⟨•⟩ α and ⟨⟨•⟩⟩ α , where α denotes either the solid phase, s, or the fluid phase, f , as where | • | denotes the size of •, the size being a measure of volume, area or length, depending on •. Here, we chose to split the fully resolved domain Ω into N adjacent subdomains as Ω = ∪ N i=1 Ω ,i , such that each Ω ,i is an RVE that is assumed to exhibit a periodic geometry. We also introduce the decomposition Ω = Ω f ∪ Ω s , as depicted in Fig. 3(a), where Ω f is the fluid part and Ω s is the solid part. Furthermore, we introduce the boundary Γ = Γ f ∪ Γ s . Finally, we introduce the union of all internal boundaries as Γ int . We note that the intrinsic averaging operator has the following important property for the split-up domain: where n α is the volume fraction of phase α, whereas n α Γ is the surface fraction of phase α. Note that n f is the porosity of the RVE. Without altering the problem, we now reformulate the integrals in Eqs. (22) and (25) as

Macroscale problem
We choose first order homogenization, i.e p M and u M are represented within each RVE (for X ∈ Ω ) in terms of macroscale fieldsū andp as follows: wherep =p(X) is the macroscale pore pressure,ḡ def = (∇p)(X) is the macroscale pore pressure gradient, H def = (ū⊗∇)(X) is the macroscale displacement gradient,X f andX s are the center points for Ω f and Ω s , respectively, andX is the macroscopic coordinate being the center of the whole Ω . 2 The macroscale momentum balance in Eq. (23a) thus becomes: where the macroscopic total stress tensorP{H,p,ḡ} is computed as Similarly, the macroscale mass balance equation (Eq. (24)) becomes  where we have introduced the material fluid velocity W = J w · F −T , the prescribed material fluid velocityŴ, and the storage terms which are the relative pore space and the first moment thereof, respectively.Φ (2) can be neglected for large scale separation, which is henceforth assumed. We may also introduce the macroscale seepage velocity (of the first Piola-Kirchhoff type) as follows: The complete macroscale problem is now given as: Find (ū,p) ∈Ū ×P such that is the prescribed seepage on Γ w . The macroscale boundary is decomposed as Γ = Γ u ∪ Γ t = Γ p ∪ Γ w . Furthermore, we introduce the solution and test spaces Note thatū andp are introduced on the whole of Ω and not just the solid or fluid part.

Preliminaries
Suppose that the macroscale solution fieldsū andp are known in a given (nested) iteration stage of the macrosubscale problem. For any given macroscale point, we may thus compute the values ofH,p andḡ, which are fed into the subscale problem formally defined by Eq. (22). This problem is localized to the RVE in what follows.
Along with the assumption of separation of scales, it is natural to assume that the transient of the fluctuating fields is very rapid within the RVE. As a consequence, we introduce the approximation of stationary subscale response and, therefore, neglect the term d t u, cf. [10]. This is also in complete analogy with the initial assumption of a stationary flow. As another preliminary, we choose Dirichlet boundary conditions on the displacement, i.e. we choose u µ = 0 on Γ .
As a result, the (still un-closed) problem reads as follows: Find (u, w, p) ∈ U × W × P such that where we introduce the RVE-spaces Moreover, t R is the reaction force due to the cutting up of the domain and is given as where the subscale fluid reaction term is defined as t R,sub

Closing the RVE-problem-Auxiliary constraints
To close the RVE-problem, a number of extra conditions (constraints) must be introduced as complements to Eq. (36). Firstly, in order to enforce the proper mean pressure we add the condition corresponding to the Lagrange multiplier λ. Secondly, we adopt weakly periodic boundary conditions on w and p µ .
[Note that we have already assumed homogeneous Dirichlet condition on u µ .] For a background discussion on weakly periodic boundary conditions, we refer to [11]. We formulate the condition on p µ as The jump operator , where X is a coordinate on the positive (mirror) part on Γ , denoted Γ + henceforth, whereas X − (X) is the corresponding point on the negative (image) part of Γ . For the seepage velocity w, we pose the weakly periodic boundary condition as where we recall that J F −1 · w = W is the pull-back of w to the material configuration. Altogether, the conditions on weak periodicity give two Lagrange multipliers: β ∈ B := [L 2  Γ f+  ] 3 for the seepage velocity condition and for the pressure condition.

Operational format of the RVE-problem
The system of equations that defines the RVE-problem in its operational format becomes: Find (u, w, p, λ, γ , β) ∈  Here we note that, due to the choice of boundary conditions, the subscale part of t R vanishes on both Γ s and Γ f .

RVE-model
Numerical examples are presented subsequently for an RVE with the simple subscale geometry shown in Fig. 4. The RVE is of cubic shape with the side-length 1 m. All faces of the RVE are identical with a quadratic hole located in the center. The hole measures 0.33 × 0.33 m 2 ; consequently, the porosity of the undeformed RVE is n f = 0.26. The following constitutive equations were adopted for the two phases: Neo-Hooke hyperelasticity was assumed to represent the intrinsic deformation characteristics of the solid phase, whereby the parameter G s , K s are recovered in the small strain limit for the deviatoric and volumetric response, respectively. The used values for G s and K s vary in the numerical examples (as given below). The fictitious elastic material is also represented by a Neo-Hooke hyperelastic material. Here, we chose both the bulk and shear modulus as 50 Pa.
Standard Newtonian (linear) viscous flow was assumed for the fluid flow. The corresponding first Piola-Kirchhoff stress is given as Due to the assumption d t u = 0, we have v = w in practice. The viscosity µ is set to 1 Pa s. The Lagrange multipliers γ and β associated with the weakly periodic boundary conditions are approximated using global polynomials of order 2. FE-computations were carried out on a 3D-mesh of tetrahedrons, whereby the following element approximations were used: For the solid skeleton and the fictitious elastic material in the pore space, u is quadratic. For the fluid phase, we adopt the Taylor-Hood element with quadratic w and linear p. The mesh were generated using a criteria on the maximum edge length of 0.09 on the elements.
All computations are performed using the open source Finite Element solver OOFEM [12].

Assessment of the Biot coefficient
From [1], the total macroscopic stress in the case of small deformations and linear elastic solid skeleton can be given as where α is the (constant) Biot coefficient, andĒ is the homogenized overall stiffness. Thus, α measures how the pressure in the fluid affects the mean stress in the solid. Theoretically, α varies from n f (no interaction) to 1 (full interaction). Here, the term "no interaction" can be thought of as the case where a rigid interface boundary Γ int prevents the pressure in the fluid to contribute to the stress in the solid. Full interaction implies that the fluid pressure is "extrapolated" onto the solid. The effect of α should appear in the present definition ofP ifP{H,p,ḡ} is properly linearized and the result is evaluated for small values ofp. However, to simplify calculations, we rather redefine α as the secant of the total mean stress in the following sense Hence, α corresponds to the Biot coefficient for small values ofp in the light of the ansatz in Eq. (45). Whereas the Biot coefficient classically is defined only for linear problems, α thus defines the interaction also for non-linear formulations.
In the following example, we evaluate α for the special case whenH = 0 andḡ = 0. SinceP m (0, 0, 0) = 0, we thus obtain As to parameter values, we choose the shear modulus G s = 50 Pa and let the bulk modulus K s vary. Fig. 5 shows how α changes withp for different choices of the bulk modulus K s . For very low values of K s , it turns out that α approaches the porosity, n f , while for large values, it is close to 1. For low values of K s , α depends significantly onp while for large values it is quite insensitive, i.e. the ansatz in Eq. (45) is validated. Fig. 5(b) gives the same information in a different way. Fig. 6 shows the distribution of stress and deformation in the RVE due to the prescribed macroscale pore pressure.

Assessment of the apparent permeability
The apparent permeability is investigated in terms of the seepage velocityw for unit pressure gradient, i.e. for g = −e 1 . Clearly the result will be affected by the state (H,p). More specifically, Fig. 7 shows howw 1 varies with the macroscopic displacement gradientH =H 11 e 1 ⊗ e 1 orH =H 12 e 1 ⊗ e 2 (representing uniaxial stretch and simple shear) for increasing values ofH 11 andH 12 , respectively. We choose the bulk modulus K s = 50 Pa and let the shear modulus G s vary. The behavior is clearly nonlinear, as shown in Fig. 7. Figs. 8 and 9 show snapshots of the deformed pore space when the RVE is subjected to different modes of macroscale deformation.

Influence of RVE-size
Next, we consider RVEs consisting of several unit cells, positioned in an array as illustrated in Fig. 10. In order to achieve an objective comparison, we choose constant values of the Lagrange multipliers γ and β, i.e. 0th order. Both G s and K s are chosen as 50 Pa. The prescribed mean pressurep is 20 Pa. As in the previous computations, a unit pressure gradient in the e 1 -direction is imposed on the RVE. Due to the high computational cost, the mesh used in this numerical example is somewhat coarser than the mesh used in the previous computations. The maximum length of an element edge was chosen as 0.1. Fig. 11(a) shows how the seepage velocity changes with strain for the three RVE-sizes. The difference between the results presented here and the corresponding results in Fig. 7 is due to the choice of lower order weakly periodic boundary conditions and a coarser mesh. In absolute terms, however, the permeability differs only a little between the RVEs. Fig. 11(b) shows how the current (deformed) porosity φ changes with strain. The observed trend for increasing RVE-size can be explained by the fact that the overall flexibility of the solid skeleton increases as the constraining boundary effect is reduced.

Conclusions and outlook
A novel approach to modeling the fluid-structure interaction in porous media, based on computational homogenization, is presented. A monolithic FSI problem, representing finite strains, is devised. By introducing a  fictitious elastic material in the pore system, we prevent the computational mesh from being overly distorted. Special care is taken such that the fictitious material does not contribute to the overall stiffness. By homogenizing the FSI problem, we arrived at a coupled macroscale equation that is recognized in the porous media literature.
In the numerical examples, the classic Biot coefficient is evaluated, and it is shown to remain within its theoretical bounds. Furthermore, the dependence of deformation and pressure on the macroscale permeability is assessed. The numerical examples also verify that the mesh-smoothing method prevents distorted elements.
A list of future extensions of the model framework can be made long due to the large number of complex physical phenomena that take place in porous media. For example, we could choose to incorporate contact models to handle pore collapse, models for erosion of the internal surfaces and adding more phases. However, as many of these processes are very complex and, generally, computationally expensive to simulate, a natural first step is to reduce the computational cost. To that end, assessing non-linear solvers and investigating Dirichlet boundary conditions on the fluid phase are pressing matters.
As of now, the model has only been evaluated using an elastic solid material. Thus, the model should be extended to other types of subscale material models, such as visco-plasticity, in order to be able to simulate the behavior of e.g. soft clay, gels, biological tissues etc. Hence, the theoretical framework must be extended to include history effects.