Effective balance equations for poroelastic composites

We derive the quasi-static governing equations for the macroscale behaviour of a linear elastic porous composite comprising a matrix interacting with inclusions and/or fibres, and an incompressible Newtonian fluid flowing in the pores. We assume that the size of the pores (the microscale) is comparable with the distance between adjacent subphases and is much smaller than the size of the whole domain (the macroscale). We then decouple spatial scales embracing the asymptotic (periodic) homogenization technique to derive the new macroscale model by upscaling the fluid–structure interaction problem between the elastic constituents and the fluid phase. The resulting system of partial differential equations is of poroelastic type and encodes the properties of the microstructure in the coefficients of the model, which are to be computed by solving appropriate cell problems which reflect the complexity of the given microstructure. The model reduces to the limit case of simple composites when there are no pores, and standard Biot’s poroelasticity whenever only the matrix–fluid interaction is considered. We further prove rigorous properties of the coefficients, namely (a) major and minor symmetries of the effective elasticity tensor, (b) positive definiteness of the resulting Biot’s modulus, and (c) analytical identities which allow us to define an effective Biot’s coefficient. This model is applicable when the interactions between multiple solid phases occur at the porescale, as in the case of various systems such as biological aggregates, constructs, bone, tendons, as well as rocks and soil.


Introduction
The Theory of Poroelasticity [5][6][7][8] is a widely known modelling framework which is usually embraced to model the effective mechanical behaviour of a fluid-filled porous elastic structure. There exists a large variety of physical systems where porescale interactions between a deformable solid and a fluid phase take place, thus motivating a poroelastic modelling approach. Key examples include hard hierarchical tissues, such as the bone [18,51], the interstitial matrix in healthy and tumorous biological tissues (see, e.g. [9,23]), the human eye [12,28], artificial constructs and biomaterials [13,29], as well as rocks and soil [30,50].
Physical systems of this kind are typically multiscale in nature. In fact, they are characterized by a porous structure where the interactions between the solid and fluid phases occur on a scale (the porescale) which is much smaller than the average size of the whole system (the macroscale).
The equations governing a poroelastic medium can be derived by upscaling the coupled fluid-structure porescale balance equation by means of suitable homogenization techniques. The latter include effective medium theory, mixture theory, volume averaging and asymptotic homogenization, as thoroughly discussed, for example, in the review [4].
The two former techniques represent micromechanical approximations which can provide an estimate of the macroscale coefficients for specific geometries of the pores (such as ellipsoidal, penny-shaped or spherical), or for diluted pores (see, e.g. [15]), as well as closed form analytic formulas relating drained and undrained poroelastic coefficients, as those reported for example in [4,33,50]. Volume averaging approaches are suitable for deriving the functional form of the macroscale equations and are based on relationships between the porescale and macroscale energy of the system, see also [27]. However, the macroscale coefficients are usually not related to the underlying microstructure, and physical arguments and/or experimental data are to be supplemented in order to determine them. Examples of problems concerning specific inclusions (such as rigid discs) embedded in an elastic and poroelastic matrix and their analytic solutions by means of suitable Green's functions are also found in [1,2,20,21,24,35,48].
The asymptotic homogenization technique [3,25,34,40] exploits the sharp length scale separation which exists in the system to decouple spatial variations that occur on different scales. The relevant fields (e.g. velocities, displacements, and pressures) are expressed in terms of power series of the ratio between representative microscopic and macroscopic length scales. The resulting governing equations describe the behaviour of the system in terms of the leading (zeroth)-order fields in a homogenized macroscale domain where the microstructure is smoothed out. This technique is in general characterized by a higher algebraic complexity compared to average fields approaches; however, it provides a precise prescription for computing the coefficients of the model. These are based on the solution of microstructural differential problems where the constitutive behaviour and geometrical arrangement of the individual phases is clearly specified.
In [11,32], the authors applied the asymptotic homogenization technique to the fluid-structure interaction problem between a linear elastic matrix and a Newtonian fluid flowing through the pores obtaining the Biot's equations of poroelasticity. The formulation has been further extended in [22] for inviscid fluid, as well as in [37] to account for appositional finite growth between the solid and fluid phases.
The standard formulation of the equations of poroelasticity deduced via asymptotic homogenization is derived accounting for an elastic matrix interacting with a slowly flowing Newtonian fluid.
In this work, we aim to determine the effective behaviour of a material where the microstructure comprises both an elastic fluid-filled porous matrix, and a number of embedded elastic inclusions/fibres which can interplay with both the matrix and the fluid phase. Our chief motivation is the study of poroelastic composites, i.e. complex, multiscale physical systems where multiple adjacent elastic phases interacting with a fluid can be identified on the porescale. For example, this is the case in the biological tissues interstitial matrix, which can be considered a composite made of multiple constituents, such as cells and different types of collagen fibres [31], which are interacting with the fluid flowing in the interstitial space [23] on the porescale [47]. A similar scenario is encountered when dealing with hard hierarchical materials, such as bone and tendons, where water is interplaying with both collagen and mineral at the finest hierarchical levels of organization [51], as well as with the constituents of the osteonal structures [18]. In [46], the authors present an asymptotic homogenization analysis concerning the upscaling of a fluid-structure interaction problem motivated by experimental research on fractures. The analysis is conducted by considering two elastic phases and one fluid phase and restricted to a simplified scenario where one of the solid phase envelops the other, and the fluid is only in contact with this latter one. The upscaling is then carried out assuming that there is a non-welded interface between the two solids, and this is reflected by prescribing the solid stress across this interface (which is assumed to be continuous) as being proportional to the jump in the elastic displacements and velocities.
In the present work, we embrace the asymptotic homogenization technique to upscale the interaction between a number of linear elastic inclusions and/or fibres embedded in a porous matrix saturated with an incompressible Newtonian fluid. Both the elastic matrix and the various subphases are in general assumed to be interacting with the fluid flowing in the pores. We assume that the scale at which the various solid subphases are clearly resolved is comparable with the porescale and collectively denote it as the microscale, which is in turn assumed to be much smaller than the macroscale. The upscaling is then carried out by accounting for continuity of stresses and displacements across the interfaces between solid phases, and continuity of stresses and velocities across the fluid-solid interfaces. The resulting governing equations, which are presented in the quasi-static case, are of Biot's type and read as a generalization of both the standard formulation for linear elastic composites [39], and standard poroelasticity [11]. These two models are recovered when assuming that the matrix is not porous and that no elastic phase other than the matrix is present, respectively. The coefficients of the model encode the properties of the microstructure and are to be computed by solving local differential problems which combine the properties of those arising in the context of multiscale composites [38,39] and poroelasticity [19].
The paper is organized as follows. In Sect. 2, we introduce the fluid-structure interaction problem which describes the interplay between the elastic matrix, subphases, and fluid percolating in the pores. In Sect. 3, we perform a multiscale analysis of the system of partial differential equations illustrated in Sect. 2 and derive the new macroscale model which governs the homogenized mechanical behaviour of poroelastic composites. In Sect. 4, we discuss the macroscale results and prove rigorous properties of the arising effective elasticity tensor, Biot's modulus, and Biot's coefficient. In Sect. 5, we conclude our work by discussing the limitations of the model and further perspectives.

A multiphase fluid-structure interaction problem
We begin by considering a set Ω ∈ R 3 , where Ω represents the union of a solid porous matrix Ω II , a connected fluid compartment Ω f , and a collection Ω I of N disjoint subphases (which could represent either inclusions or fibres) Ω α , such that andΩ =Ω I ∪Ω II ∪Ω f . A sketch of a cross section of the three-dimensional domain Ω is shown in Fig. 1. The balance equations in the solid domains Ω α and Ω II , by neglecting volume forces and inertia then read, and The symbols T α and T II that appear in relationships (2, 3) denote the solid stress tensors corresponding to each subphase Ω α and the one corresponding to the matrix Ω II , respectively. We then assume that both the matrix and each subphase are anisotropic linear elastic solids, so that the constitutive equations for T α and T II are given by where u α and u II are the elastic displacements in each subphase and the matrix, respectively.
The fourth rank tensors C α and C II are the elasticity tensors in each subphase and the matrix, respectively, with corresponding components C α i jkl and C II i jkl , for i, j, k, l = 1, 2, 3. We note that each C α and C II are equipped with right minor and major symmetries, namely and therefore also left minor symmetries follow by combining (6,7). In particular, by applying right minor symmetries we can equivalently rewrite constitutive equations (4,5) as where is the symmetric part of the gradient operator. The balance equation in the fluid compartment reads where T f is the fluid stress tensor. We assume that the fluid compartment is an incompressible Newtonian fluid, so that the constitutive equation for T f is given by where v denotes the fluid velocity, p the pressure, and μ the viscosity, together with the incompressibility constraint Substituting relationship (12) in (11) and using the the divergence-free condition (13) yields the Stokes' problem In order to close the fluid-structure interaction problem in the whole domain Ω, we also require interface conditions between the fluid and the solid phases. We first define the interface between the fluid phase and the α inclusion/fibre as Γ α := ∂Ω α ∩ ∂Ω f and the interface between the matrix and the fluid phase as Γ II := ∂Ω II ∩ ∂Ω f . We then impose continuity of velocities and stresses across each Γ α and Γ II , namelẏ ∀α = 1 . . . N , whereu α andu II are the solid velocities in each subphase Ω α and the matrix Ω II , respectively. The unit outward (i.e. pointing into the fluid domain Ω f ) vectors normal to the interfaces Γ α and Γ II are denoted by n α and n II , respectively. Finally, we require continuity of stresses and displacements across the interface between each elastic subphase and the matrix. We define this boundary as Γ αII := ∂Ω α ∩ ∂Ω II , so that T α n αII = T II n αII on Γ αII . (19) ∀α = 1 . . . N , where n αII is the unit vector normal to the interface Γ αII pointing into the fibre/inclusion Ω α .
In the next section, we perform a multiscale analysis by (a) non-dimensionalizing the partial differential equations (PDEs) described in this section and introducing two well-separated length scales, (b) applying the asymptotic homogenization technique to the resulting non-dimensional systems of PDEs, and (c) deriving the effective governing equations for the material as a whole.

Multiscale analysis
We now perform a multiscale analysis of the fluid-structure interaction problem introduced in the previous section, which is summarized below where, by means of the constitutive relationships (8), (9), and (12), together with the incompressibility constraint (24), the balance equations (21), (22), and (23) can also be rewritten as (21)(22)(23)(24)(25)(26)(27)(28)(29)(30) is then to be closed by prescribing suitable external boundary conditions on ∂Ω. We assume that there exist two typical length scales in the system. In particular, we denote the average size of the whole domain Ω by L (the macroscale), while d refers to the porescale (the microscale), which in this work is assumed to be comparable with the distance between adjacent subphases interacting with the matrix and the fluid domain. In order to emphasize the difference between such scales, it is helpful to perform a non-dimensional analysis of the system of PDEs (21-30).

Non-dimensional form of the equations
We carry out the non-dimensional analysis by assuming that the system is characterized by a reference pressure gradient C, and that the characteristic (reference) fluid velocity is given by the typical parabolic profile proportional to that of a Newtonian fluid slowly flowing in a cylinder of radius d. This is the appropriate scaling that captures the scale separation between the microscale d and the macroscale L in a porous domain, as also discussed in [37]. Although different scaling choices for the fluid velocity are formally possible, these would not account for the appropriate effective behaviour of a fluid flowing through a porous solid matrix. An example of alternative choices which lead to an effective viscoelastic-type behaviour is illustrated in [11]. Therefore, in our case we have The time is scaled by where is the reference parabolic fluid profile which is embraced to upscale a fluid-structure interaction problem to a poroelastic problem, as also explained in [37]. Equation (35) therefore represents the (macro) reference time scale for a fluid slowly flowing in the pores and is assumed to be analogous to the time scale for deformation of the elastic material, for the sake of consistency with continuity of velocities at the interface. We then exploit (34) and observe that to obtain the non-dimensional form of the system of PDEs (21)(22)(23)(24)(25)(26)(27)(28)(29)(30), namely where we have dropped the primes for the sake of simplicity of notation. The nondimensionalized counterparts of constitutive relationships (8), (9), and (12) are given by so that the balance equations (38-40) rewrite In the next section, we introduce the asymptotic homogenization technique which is used to upscale the nondimensional system of PDEs (38)(39)(40)(41)(42)(43)(44)(45)(46)(47) by formally assuming that the microscale and the macroscale are well separated.

The asymptotic homogenization technique
In this section, we introduce the two-scale asymptotic homogenization technique which is used to derive a macroscale model for Eqs. (38)(39)(40)(41)(42)(43)(44)(45)(46)(47). We first assume that the microscale (where the pores and individual inclusions/fibres are clearly resolved), denoted by d, is small compared to average size of the domain L, i.e.

(55)
We then introduce a local spatial variable to capture microscale variations of the field via setting The spatial variables x and y are to be considered formally independent and represent the macroscale and the microscale, respectively. The gradient operator then transforms as We further assume that all the fields u II , u α , v, p, T f , T α , T II , as well as the elasticity tensors C II and C α , ∀α = 1 . . . N , are functions of both x and y. We also assume that the fields u II , u α , v, p, T f , T α , T II can be represented in terms of a series expansion in powers of , i.e.
where ϕ collectively denotes each field involved in the present analysis.

Remark 1 (Microscale periodicity)
We are now interested in obtaining a closed system of PDEs in terms of the leading (zeroth)-order pressure, velocity, and displacement fields. This is done by applying relationships (57) and (58) to system (38)(39)(40)(41)(42)(43)(44)(45)(46)(47) and constitutive equations (48)(49)(50), and in turn to relationships (51)(52)(53). This way, by equating coefficients of the same power of , we can obtain various sets of differential conditions which can be combined to obtain a system of PDEs which holds on the macroscale x ∈ Ω H . The domain Ω H represents the homogenized domain where microscale heterogeneities are smoothed out. The coefficients of the model are then typically expressed averaging the solution of appropriate microscale local problems. As each coefficient ϕ (l) that appears in (58) is required to be well-defined for arbitrary small values of , it is in principle necessary to assume that all the fields are locally bounded, i.e. finite with respect to the microscale variable y when → 0, see also [25,40]. This is the least restrictive assumption that is to be embraced to successfully perform upscaling of a given system of PDEs when dealing with formal asymptotic homogenization. In [11], the authors derive Biot's equations of poroelasticity by assuming local boundedness of the fields. However, this approach is appropriate when the main goal is the functional form of the macroscale model only. This is since the prescriptions of the coefficients obtained this way are related to microscale problems which are, in principle, to be solved on the whole microstructure. Therefore, they cannot be used in practice unless further geometrical restrictions, such as periodicity of the microstructure, are imposed, as indeed suggested in [11]. We therefore assume that every ϕ (l) , C α , and C II are y-periodic. This latter technical assumption allows us to restrict the analysis of the microstructure to a single periodic cell, which could in any case contain several different subphases characterized by different geometry, arrangement and elastic properties, as shown in Fig. 2.

Remark 2 (Macroscopic uniformity)
The microscale geometry can in principle vary with respect to the macroscale; however, this potential dependence is often (most of the time implicitly) neglected. Here, we assume that the medium is macroscopically uniform, i.e. the microscale geometry does not depend on the macroscale variable x. In particular, this assumptions leads to the straightforward differentiation under the integral sign Whenever Ω = Ω(x), Eq. (59) is not satisfied, and application of the generalized Reynolds' transport theorem may lead to additional macroscale contributions, see, e.g. [25,36,37].
Finally, for the sake of clarity of presentation and without loss of generality with respect to the properties of the model, we can restrict our analysis by assuming that only one subphase is contained in each periodic cell, as shown in Fig. 3. The model can be easily extended to a number of subphases within the periodic cell if necessary for a particular application, as done in the context of simple elastic composites in [38]. Therefore, the index α is no longer needed and we adjust the notation accordingly. We identify the domain Ω with its corresponding periodic cell, with fibre/inclusion, matrix, and fluid cell portions denoted by Ω I , Ω II , and Ω f , respectively. The interfaces between the different phases are then denoted by Γ I := ∂Ω I ∩ ∂Ω f , Γ II := ∂Ω II ∩ ∂Ω f , and Γ III := ∂Ω I ∩ ∂Ω II , with corresponding unit normal vectors n I , n II , and n III , where Γ I is the interface between the fluid and the fibre/inclusion, Γ II is the interface between the fluid and the matrix and Γ III is the interface between the matrix and the fibre/inclusion.

Derivation of the macroscale model
We apply the asymptotic homogenization assumptions (57) and (58) to Eqs. (38)(39)(40)(41)(42)(43)(44)(45)(46)(47)(48)(49)(50)(51)(52)(53) to obtain, accounting for periodicity, the following multiscale system of PDEs We have one subphase shown in blue that is in contact with the porous, solid, elastic matrix shown in red, and the fluid flowing in the pores is shown in white. We also highlight the interfaces Γ I which is shown in green between the subphase and the fluid, Γ II shown in black between the matrix and the fluid and Γ III shown in grey between the subphase and the matrix (colour figure online) equipped with multiscale constitutive equations for the fluid and solid stress tensors T f , T I , T II , given by while the balance equations in terms of the elastic displacement, fluid velocity, and pressure u I , u I , v , p read We can now substitute power series of type (58) into the relevant fields in (60-75). Then by equating the coefficients of l for l = 0, 1, . . ., we derive the macroscale model for the material in terms of the relevant leading (zeroth)-order fields. Whenever a component in the asymptotic expansion retains a dependence on the microscale, we can take the integral average, which we define as where ϕ is a field, and also where the integral average can be performed over one representative cell due to y-periodicity and |Ω| is the volume of the domain and the integration is performed over the microscale. We note that |Ω| = |Ω f | + |Ω I | + |Ω II |. Due to the assumption of y-periodicity, the integral average can be performed over one representative cell. Therefore, (76) represents a cell average. For the sake of brevity, we also introduce the notation for fields ϕ with components ϕ I and ϕ II defined in the solid cell portions Ω I or Ω II , respectively. Equating coefficients of 0 in (60-69) we obtain and the constitutive equations (70-72) for T f , T I , T II have coefficients of 0 and the balance equations (73-75) have coefficients of 0 Similarly, we now wish to equate the coefficients of 1 in Eqs. (60-69) which gives and the constitutive equations (70-72) for T f , T I , T II have coefficients of 1 and the balance equations (73-75) have coefficients of 1 We can now see from (80) and (88) that the leading order pressure p (0) does not depend on the microscale y. That is We also have from (89) and (90) that u II , which are the leading order solid displacements, are rigid body motions and therefore, by y-periodicity, do not depend on the microscale y. That is Since we have the boundary condition u (0) I = u (0) II on Γ III given by (87) we can define which we will use throughout the following sections.

Fluid flow on the macroscale
We now wish to investigate the leading order of the velocity which we denoted v (0) . We can define the relative fluid-solid displacement, w, by Using Eqs. (88), (82), (83), (96) and (104), exploiting notation (113), we have a Stokes'-type boundary value problem which is given by Now exploiting linearity and using (110), we can propose the following ansatz for the Stokes-type boundary value problem (115-117), where p (1) is defined up to an arbitrary y-constant function. Equation (119) is the solution to problem (115-117) provided that second rank tensor W and vector P satisfy the following cell problem where periodic conditions apply on the boundary ∂Ω f \Γ I ∪ Γ II and a further condition is to be placed on P for the solution to be unique (for example zero average on the fluid cell portion). Taking the integral average of (118) over the fluid domain leads to meaning that the fluid flow is described by Darcy's law in the macroscale. As expected, this is the same result as standard poroelasticity with only one solid phase.

Poroelasticity on the macroscale
We now require the macroscale equations to close the system for the elastic displacement u (0) and p (0) . Summing up the integral averages of Eqs. (94), (95) and (96) we have Applying the divergence theorem to the first three integrals and rearranging the last three integrals by means of macroscopic uniformity (59), we obtain where n I , n II , n III , n Ω I \Γ I ∪Γ III , n Ω II \Γ II ∪Γ III and n Ω f \Γ I ∪Γ II are the unit normals corresponding to Γ I , Γ II , Γ III , ∂Ω I \Γ I ∪ Γ III , ∂Ω II \Γ II ∪ Γ III and ∂Ω f \Γ I ∪ Γ II . Since the contributions over the external boundaries of Ω I , Ω II and Ω f cancel out due to y-periodicity (123) becomes The first six integrals in (124) cancel out due to (100), (101), and (102), and the final three terms can be written as where φ := |Ω f |/|Ω| is the porosity of the material.
Exploiting (110) and (113), we can write the following problem for u The solution to the problem given by (126-131), exploiting linearity is given as where A I and A II are third rank tensors and a I and a II are vectors. The auxiliary fields A I , A II , a I and a II solve the following cell problems.
and ∇ y · (C I ξ y (a I )) = 0 in Ω I (140) ∇ y · (C II ξ y (a II )) = 0 in Ω II (141) (C I ξ y (a I ))n III = (C II ξ y (a II ))n III on Γ III (142) a I = a II on Γ III (143) (C I ξ y (a I ))n I + n I = 0 on Γ I (144) (C II ξ y (a II ))n II + n II = 0 on Γ II We can also write the cell problems for A I , A II , a I and a II with corresponding components A I ikl , A II ikl , a I i and a II i as and where we have used the notation We should note that the problems in terms of a i and A i , where i = I, II, are to be solved on the solid cell portion Ω s , whereΩ s :=Ω I ∪Ω II .

Remark 3 (Compatibility condition for the cell problems)
We have the compatibility condition (also known as the solvability condition) for the cell problems, see, e.g. [16]. We first take the integral average of the sum of the left hand sides of (134) and (135) and apply the divergence theorem. That is, where n I , n II , n III , n Ω I \Γ I ∪Γ III and n Ω II \Γ II ∪Γ III are the unit normals corresponding to Γ I , Γ II , Γ III , ∂Ω I \Γ I ∪ Γ III and ∂Ω II \Γ II ∪ Γ III . Terms on the external boundaries of Ω I and Ω II cancel due to y-periodicity. So (160) becomes Where we have used the interface condition (136) and we can now also use (138) and (139) to write (161) as On the other hand we are able to rewrite (159) using (134) and (135) as We can then apply the divergence theorem to obtain Where the unit normals are as above. The contributions over the periodic boundaries cancel to give we can see that (162) and (165) are equal, and therefore, this proves the compatibility condition which is necessary for the problem to admit a solution, which can be made unique by imposing an additional condition, such as zero average of the auxiliary variables on the periodic cell.
We now consider the leading order solid stress tensors. Since from (105) and (106) we have that u (1) I and u (1) II are related to T and T (0) where we define Adding (166) and (167) and taking the integral average over the solid domain gives From (125), we have that where T Eff = T Equations (170) and (171) represent the average force balance equations for our poroelastic composite material. We now return to (97), the incompressibility condition and integrate to obtain Applying the divergence theorem twice to the first integral and using (98) and (99) and also rearranging the second integral, we obtain Therefore, we have Tr(ξ y (u (1) Using (132) and (133) with (168), we have that So using (176) then Eq. (175) becomes Now returning to (114), the expression for relative fluid-solid displacement, and taking the integral average over the fluid domain, we obtain where φ is the porosity of the material. Then rearranging, we obtain We then use this relation to rewrite (177) as We can expand the left hand side of (180) and then rearrange to obtain the following expression forṗ (0) . We note that we are able to express φ∇ x ·u (0) as φI : ξ x (u (0) ). Then, We can then defineM and then we can use (182) to write (181) aṡ Finally dividing through byM, we obtaiṅ We have now derived all the equations required to be able to state our macroscale model for a poroelastic composite. Within the next section, we will state our novel macroscale model for poroelastic composites and will then prove properties of the effective coefficients of this model.

The macroscale result and properties of the effective coefficients
The equations in the macroscale model describe the effective poroelastic behaviour of the material in terms of the pore pressure, the average fluid velocity, and the elastic displacement. The macroscale model is then given by where we have that p (0) is the macroscale pressure, u (0) is the solid displacement,u (0) is the solid velocity, and w is the leading order fluid velocity. The first equation of the macroscale model represents Darcy's law for w, where w is the relative fluid velocity. This is the same equation that we would obtain for standard poroelasticity. The second of the macroscale system of PDEs is the stress balance equation for the poroelastic composite material. The constitutive equation is of poroelastic type, with drained elasticity tensor given by The final equation in our macroscale model describes the conservation of mass for a poroelastic material and relates changes in the fluid pressure to changes in the fluid and solid volumes. We therefore have that the mechanical behaviour of the poroelastic composite material can be fully described by the effective elasticity tensor C, the hydraulic conductivity tensor W f , the tensorα, and the scalar coefficientM.
Our new model has a key difference from the model of classical poroelasticity. That is, our model is able to account for multiple elastic phases interacting at the porescale, whereas the model for classical poroelasticity is applicable when the matrix can be approximated as homogeneous at the porescale. The addition to the model of the extra interactions between multiple phases is particularly beneficial to physical applications. For example, in bones water is interplaying with both collagen and mineral at the finest hierarchical levels. It is useful to be able to account for the mineral and collagen fibres separately in the model, especially for numerical simulations, as both constituents have very different elastic and mechanical properties. The differences in elastic and mechanical properties are accounted for by the multiple elasticity tensors C I , C II and by M I , M II , Q I and Q II in the coefficients of the model. A similar argument is applicable to biological tissues such as the interstitial matrix, which is a composite consisting of different types of cells and collagen fibres with fluid flowing in the interstitial space; we again believe that we would obtain a much more representative description of the behaviour of the material by accounting for the multiple constituents and their varying properties in the model. Another key difference between the model in our work and classical poroelasticity is the model coefficients. Here, we propose novel cell problems from which the model coefficients are calculated. The cell problems given in our work are different from the cell problems in classical poroelasticity and also the cell problems for elastic composites. Overall, our new model reads as a comprehensive framework to describe materials where their constituents cannot be assumed to be homogeneous at the porescale.

Remark 4 (Limit cases for the macroscale model)
It is important to note that our macroscale model (185) reduces to previously obtained results when we consider the following limit cases. In the limit of only one elastic phase then this macroscale model reduces to the macroscale model for a standard poroelastic material (See the no growth limit in [37], as well as [11,34]). In this case, our model would retain all four equations presented in (185); however, the coefficient of p (0) in the third equation, the effective elasticity tensor C, the tensorα, and the scalar coefficientM would reduce to only one elastic phase. This is because the contributions on the interface between the different elastic phases are no longer present in this case. We also note that in the limit of zero fluid (no pores) then this macroscale model reduces to the macroscale model for a standard elastic composite (See [39]). In this case, the mechanical behaviour of the system is entirely described by the equations for the balance in the solid phase, the pressures, and fluid velocity reduces to zero, and the cell problem does not comprise any contribution related to the fluid phase, i.e. the only relevant cell problem (i.e. which results in non-trivial solutions) is the one involving the discontinuity in the elastic constants (134-139).
That is, the model would coincide with the standard one described in the literature for elastic composites, see, e.g. [39].
Next, we rigorously prove the following properties of the effective coefficients. We prove a) the symmetries of the effective elasticity tensor, b) an analytical identity that allows us to identify the tensorα with Biot's coefficient tensor, and c) the positive definiteness of theM, which can therefore be identified with the Biot's modulus for the whole material. These proofs generalize those proposed in [34] for standard poroelastic materials.
Theorem 1 (Symmetries of C) The fourth rank effective elasticity tensor C defined by is major and minor symmetric.
Proof We wish to show that the effective elasticity tensor C is major and minor symmetric. That is, We know that the elasticity tensor C is major and minor symmetric. We have that the first two equalities in (188) follow by definition. That is, C i M i for i = I, II, is both left and right minor symmetric since we have the left minor symmetry of C and the right minor symmetry of M. In order to prove major symmetry, we have to show that To show this, we begin with the cell problems (146) and (147) and multiplying by A I irs , A II irs , respectively, and then integrating these terms over Ω I and Ω II , respectively. We have Then integrating by parts, we obtain Applying the divergence theorem, we obtain where n I , n II , n III , n Ω I \Γ I ∪Γ III and n Ω II \Γ II ∪Γ III are the unit normals corresponding to Γ I , Γ II , Γ III , ∂Ω I \Γ I ∪ Γ III and ∂Ω II \Γ II ∪ Γ III . The terms on the boundaries ∂Ω I \Γ I ∪ Γ III and ∂Ω II \Γ II ∪ Γ III cancel due to periodicity, and we can rewrite (192) as The terms in the bracket cancel due to the cell problems (148) and (150, 151), and we obtain Which we can rearrange to obtain Rewriting the second equality in (195) interchanging r and s and k and l and using the symmetry C i j pq = C pqi j we obtain Since the left hand sides of (195) and (196) are the same then so to are the right hand sides. Taking i j as pq in (195) and (196), we have that as required. Therefore, C posessess major and minor symmetries.
The following two theorems relate to the resulting Biot's coefficient and Biot's modulus of the macroscale model. The coefficients of the macroscale model can be defined aŝ whereα is from (182) and β is the denominator of the Biot's modulusM also in (182). We obtain γ from (171), where γ is the coefficient of the leading order pressure in the effective stress, T Eff . We are able to give a physical interpretation of these coefficients following the descriptions given in [41]. The coefficient β (denominator ofM) can be thought of as a variation in the fluid volume in response to a variation in pore pressure.M is a poroelastic coefficient that depends on the porescale geometry and porosity and on the properties of the elastic matrix [11]. We can describeα as the ratio of change in interstitial fluid volume to changes in solid volume. We can now state the following theorem that will provide an analytical identity allowing us to define an effective Biot's coefficient.

Theorem 2 (Biot's coefficient)
We have the following analytical identity which allows us to define an effective poroelastic Biot's coefficient tensor.
Proof In index notation, we can write γ andα as We use (146) and (147) from the cell problems and multiply by a I i , a II i , respectively. We then integrate these expressions over Ω I and Ω II , respectively, to obtain Then integrating by parts, we obtain Applying the divergence theorem, and using symmetry properties (6-7), we have where n I , n II , n III , n Ω I \Γ I ∪Γ III and n Ω II \Γ II ∪Γ III are the unit normals corresponding to Γ I , Γ II , Γ III , ∂Ω I \Γ I ∪ Γ III and ∂Ω II \Γ II ∪ Γ III , and cancelling terms on the periodic boundaries due to y-periodicity and accounting for the interface conditions (148) and (150, 151) we obtain Hence we have We now wish to multiply (140) and (141) from the cell problems by A I ikl , A II ikl , respectively, and then integrate over Ω I and Ω II , respectively. Integrating by parts we obtain Applying the divergence theorem and using (144) and (145) we obtain where terms on the boundaries have cancelled due to periodicity. Then reversing the divergence theorem, we have Hence, Using (208) and (212), we have that Tr (M I + M II ) s = C I Q I + C II Q II s . Therefore using this in the definitions ofα i j and γ i j , we have that γ i j = −α i j as required.
We should note here that the coefficientsα and γ are defined by different cell problems for A I , A II and for a I , a II , respectively. By having an analytical identity as in Theorem 2, we can reduce computations as the numerics then do not have to be carried out for the cell problems involving a I , a II . Now that we have proved this analytical identity, we can use it to restate the macroscale model for a poroelastic composite material. We have We will now state our final theorem relating to the Biot's modulus of our system. Proof In order to prove this, we only need to prove that the denominator ofM which we defined as β is less than zero. We can begin by recalling β from (200) and then in index notation We can then multiply Eqs. (140) and (141) from the cell problems by a I i , a II i and integrate over Ω I , Ω II , respectively. That is, Integrating by parts, we obtain Then applying the divergence theorem we obtain where n I , n II , n III , n Ω I \Γ I ∪Γ III and n Ω II \Γ II ∪Γ III are the unit normals corresponding to Γ I , Γ II , Γ III , ∂Ω I \Γ I ∪ Γ III and ∂Ω II \Γ II ∪ Γ III . The terms on the periodic boundaries cancel due to periodicity, and then using Eqs. (142), (144), and (145) from the cell problems we obtain Accounting for y-periodicity and relationship (143), the sum of the first and third integrals above is equal to the sum of the corresponding integrals (with corresponding normals) on the boundaries ∂Ω I and ∂Ω II , so that we can apply the divergence theorem in reverse to obtain We can write this as Since the two terms on the right hand side are positive, we therefore have that and so using the integral average notation we have Since β < 0 we therefore have thatM > 0. That is, the Biot's modulus is positive definite.

Concluding remarks
We have presented a novel system of PDEs which describes the effective behaviour of poroelastic composites. In Sect. 2, we have begun by considering the quasi-static multiphase fluid-structure interaction problem which describes the mechanics of a number of linear elastic inclusions/fibres embedded in a porous, linear elastic matrix, filled by a slowly flowing incompressible Newtonian fluid. In Sect. 3, we have then enforced the length scale separation between the microscale and the macroscale to upscale the non-dimensionalized system of PDEs via asymptotic homogenization. In particular, we have assumed that both the pores and the elastic subphases (i.e. inclusions or fibres) are clearly resolved on the microscale, while the macroscale represents the average size of the macroscale domain. In Sect. 4, we show that the new model is both formally and substantially of poroelastic type, by proving minor and major symmetries of the effective elasticity tensor, positive definiteness of the macroscale Biot's modulus, and existence of a global macroscale Biot's coefficient tensor. The new model encodes the properties of the microstructure in its coefficients, i.e. the effective hydraulic conductivity tensor, elasticity tensor, Biot's modulus, and Biot's tensor of coefficients, which are to be computed by solving appropriate periodic cell problems. The latter comprise both stress jump conditions on the solidsolid interface (as in the cell problems for elastic composites), and inhomogeneous Neumann-type conditions on the fluid-solid interface (this is typical of the cell problems arising in poroelasticity). The results are derived by assuming periodicity of the microstructure and are presented by assuming that only one elastic subphase is contained in the representative periodic cell for the sake of simplicity. The new model is a natural generalization of the standard formulations for poroelastic media and composite materials derived via asymptotic homogenization, which are both recovered as particular cases.
Our model is relevant to the description of physical scenarios where the interactions between multiple elastic constituents take place at the porescale. The standard formulation of poroelasticity is appropriate when the interactions between the individual constituents of the solid phase and the fluid can be ignored, i.e. when the solid phase can be geometrically approximated as a homogeneous matrix at the porescale.
Our model also has some limitations and is open to a number of improvements that could enhance its range of applicability. The present formulation provides the effective governing equations in a quasi-static, linearized setting, and accounting for incompressibility of the fluid phase. Generalization of our model to linearized inertia and compressibility of fluid is in principle straightforward, as it could be carried out as in [11], thus resulting in the corresponding changes on the macroscale. Leading order linearized inertia would appear in the effective balance equations for the poroelastic stress. Furthermore, the definition of the effective Biot's modulus would comprise the fluid bulk modulus, while the former depends only on the properties of the microstructure when the fluid phase is incompressible, see also [37]. However, the functional form of the elastic cell problems would not be affected by such changes, while the fluid cell problem is actually not affected by considering multiple elastic phases, and is simply the same as in standard poroelasticity in both cases. A relevant system where such an extension could provide a more accurate poroelastic modelling framework is the lung, where the investigation of the acoustic properties can be used in the context of noninvasive diagnosis for pulmonary diseases, see, e.g. [49].
Mathematical modelling of nonlinear constitutive behaviour of the individual constituents is challenging in the context of asymptotic homogenization. However, there have been recent theoretical developments concerning both multiphase elastoplastic composites [42], and hyperelastic growing porous media, see [17]. Combining these approaches to extend our formulation would provide a comprehensive multiscale modelling framework for nonlinear poroelastic composites. The latter could then be used to formulate realistic predictions when large strains are relevant, as in the case of soft tissues such as arteries, see, e.g. [10,26].
The results that we have illustrated here can also serve as a basis towards a more realistic modelling of hierarchical materials characterized by multiple separated scales (see, e.g. [43,44] in the context of elastic composites). In this case, as the coefficients are to be computed at the porescale, our results could be exploited as a starting point to model the interaction between a poroelastic matrix and another fluid or elastic compartment, as done in the recent works [14,41,45,52] for vascularized and fibre/inclusion reinforced poroelastic materials, respectively.
The next natural step is to obtain solutions of the model on the basis of a given microstructure, ideally parameterized by real-world images, for example related to biological tissues or artificial constructs, as described in the Introduction. For instance, three-dimensional numerical simulations of the asymptotic homogenization cell problems for elastic composites and poroelastic materials have been recently performed in [19,38], respectively. As such, strategies developed therein could be adapted to compute the poroelastic coefficients presented here, which are obtained by solving cell problems which generalize those solved in [19,38]. This way, predictions of the model could be validated against experimental data and/or used to optimize the design of poroelastic artificial constructs.