Micromechanical analysis of the effective stiffness of poroelastic composites

Within this work we investigate the role that the microstructure of a poroelastic material has on the resulting elastic parameters. We are considering the effect that multiple elastic and fluid phases at the same scale (LMRP model (L. Miller and R. Penta, 2020)) have on the estimation of the materials elastic parameters when compared with a standard poroelastic approach. We present a summary of both the LMRP model and the comparable standard poroelastic approach both derived via the asymptotic homogenization approach. We provide the 3D periodic cell problems with associated boundary loads that are required to be solved to obtain the effective elasticity tensor for both model setups. We then perform a 2D reduction of the cell problems, again presenting the 2D boundary loads that are required to solve the problems numerically. The results of our numerical simulations show that whenever investigating a poroelastic composite material with porosity exceeding 5% then the LMRP model should be considered more appropriate in incorporating the structural details in the Young’s moduli E1 and E3 and the shear C44. Whenever the porosity exceeds 20% it should also be used to investigate the shear C66. We find that for materials with less than 5% porosity that the voids are so small that a standard poroelastic approach or the LMRP model produce the same results.


Introduction
Poroelasticity, developed by Biot (1955Biot ( , 1956aBiot ( ,b, 1962 is an extremely useful modelling framework which can be used to determine the effective mechanical behaviour of structures that comprise an elastic matrix with a permeating fluid flow. The theory is applicable to physical systems where there are interactions between an elastic solid and a fluid at the scale where both phases are distinct. This theory has been applied to the following biophysical examples. Firstly to hard hierarchical tissues, an example of which is bones and tendons (Cowin, 1999;Weiner and Wagner, 1998). It can also be applied to soft biological tissues such as the interstitial matrix in healthy and tumorous tissues (Bottaro and Ansaldi, 2012), the heart (myocardium) and artery walls (Cookson et al., 2012;May-Newman and McCulloch, 1998;Bukac et al., 2015;Chapelle et al., 2010). There are also applications to artificial constructs and biomaterials (Flessner, 2001;Karageorgiou and Kaplan, 2005;Chalasani et al., 2007).
The key relevance of studying poroelastic materials is to obtain a comprehensive understanding of the elastic properties of the given material based on the properties of the individual constituents.
The elastic properties of the poroelastic material depend on the interactions between all the various constituents involved in the structure. These interactions can be very complicated and they generally occur on the materials porescale. The porescale is much smaller than * Corresponding author.
E-mail address: Raimondo.Penta@glasgow.ac.uk (R. Penta). the scale of the entire tissue that is being modelled and this scale is called the macroscale.
It is key to a thorough understanding of the effective macroscale properties of the material to relate them to the properties and interactions of the porescale constituents. However, computationally it would be potentially impossible to resolve all the porescale details and interactions. It is for this reason that a variety of homogenization techniques have been developed. These homogenization techniques (Mei and Vernescu, 2010;Auriault et al., 2010;Holmes, 2012), aim to provide a macroscale model where the porescale details have been incorporated in the effective model coefficients. These techniques include mixture theory, volume averaging and asymptotic homogenization (Hori and Nemat-Nasser, 1999;Davit et al., 2013).
The asymptotic homogenization technique has been used to model poroelastic materials (Burridge and Keller, 1981;Wang, 2017;Lévy, 1979;Penta et al., 2020). It has more recently been extended to include growth (Penta et al., 2014), vascularized poroelastic materials (Penta and Merodio, 2017), active poroelastic materials (Collis et al., 2017;Brown et al., 2014), poroelastic composites (linear and non-linear) Penta, 2020, 2021b) and double poroelastic materials (Miller and Penta, 2021a). Numerically the role of porosity and microscale solid matrix compressibility on the mechanical behaviour of poroelastic materials have been investigated in Dehghani et al. (2018) and the macroscale behaviour investigated in Dehghani et al. (2020). Homogenization has also been used in the context of double porosity in fluid-saturated elastic media (Rohan et al., 2015). In this case the authors considered the coupling between a poroelastic phase and a fluid phase and then provide an application to compact bone in Rohan et al. (2012). It has also been used to investigate tissue perfusion (Rohan and Cimrman, 2010;Rohan et al., 2021).
Within this work we will compare the resulting elastic parameters arising from solving the LMRP model  for poroelastic composites compared to the parameters that arise from solving a model for an elastic composite where the matrix is poroelastic. In other words we determine the effect of considering the interactions of three different phases (two elastic and one fluid) at the porescale compared with the interactions of two elastic phases, one of which results from a further homogenization problem at a finer scale (Burridge and Keller, 1981;Dehghani et al., 2018). This means that when a material's microstructure comprises a matrix, embedded elastic subphases and fluid filled pores then using current models (excluding the LMRP model) means that the assumption that the matrix is homogeneous (ignore the subphases) has to be made or we must carry out a two-step process. The two-step process involves first solving the porous matrix problem and then solving a composite that comprises the subphases and the results of the porous matrix simulations (the so-called standard poroelastic approach in this work). This second approach means that even when considering the three phases, these are not all at the same scale, which is what the intended application actually possesses as a microstructure. This means that estimations of the parameters cannot be fully reliable. We therefore developed the LMRP model to remove this issue. This analysis will highlight under which circumstances the LMRP model provides a more accurate description of the effective elastic parameters of a poroelastic composite material. We can describe our computational platform for the LMRP model as robust, in the sense that it is very applicable to a variety of situations. That is, the platform can be altered for a variety of geometries including short fibres, various directions of fluid flow, a variety of different shaped inclusions and a wide range of constitutive properties of the constituents.
The paper is organized as follows. In Section 2, we introduce the LMRP model for poroelastic composites, derived in Miller and Penta (2020), and also introduce a comparative setup that focuses on a standard poroelastic (Biot's poroelasticity) approach. In Section 3 we have a variety of subsections each individually aimed at introducing the computational setup that would be required to solve the 3D and also the reduced 2D cell problems that arise from both models introduced in Section 2. In Section 4, we provide the results of our 2D simulations. In Section 5 we carry out 3D simulations for a different geometrical setup, namely the case of short fibre elastic inclusions, and obtain the elastic parameters. In Section 6, we conclude our work by discussing the limitations of the current simulations and provide further perspectives of the types of problems that the model could investigate and the biological scenarios where it would be best applied. We also have an Appendix which contains a detailed 2D reduction of the cell problems for the LMRP model for poroelastic composites for orthotropic constituents. This cases is the most general and under simplifying assumptions it can also be used as a framework for the 2D cell problems for standard poroelasticity and elastic composites with elastic properties with any possible symmetries. Within the appendix we also provide a guide to the numerical simulations and meshing as well as giving an example of the applicability of the model to investigating the elastic parameters of the human heart.

Governing equations
In this section we describe the governing equations for a poroelastic composite, the LMRP model , and the governing equations for standard Biot's poroelastic materials with elastic inclusions. In Fig. 1 we can see a comparison of the microstructure of each model setup. We note that in order to exemplify the difference that being able to account for multiple elastic and fluid phases all at the same scale has in comparison to the existing computational frameworks we have chosen the simplest geometry that considers uniaxial flow and uniaxial fibre/subphase elongation. Also for the sake of symmetry we arrange the fluid flow as four cylinders in the corners of the cell. We also chose this microstructure as it allows for the reduction of the microstructural cell problems to 2D, which we present clearly in Appendix A.1 in a way that can be used as guidance to the reader who would also like to reduce to 2D their own different cell problems.
On the left of Fig. 1 we have the porescale microstructure (periodic cell) of the LMRP model for poroelastic composites. We can see that this structure comprises a porous matrix II , and elastic inclusion I and fluid filled pores f . The interface between the matrix and the inclusion is denoted III and the interface between the matrix and the fluid is denoted II .
In Fig. 1 we have presented the different choices for the microstructures showing both the 3D and a 2D cross section of the domain.
Here we will introduce the effective balance equations for a poroelastic composite derived by the asymptotic homogenization technique in Miller and Penta (2020). The model is derived by considering the fluid-structure interaction between a linear elastic porous matrix, II , with embedded linear elastic subphases, I , with a Newtonian fluid, f , flowing in the pores. The fluid-structure interaction problem consists of balance equations for each elastic domain and the fluid domain, as well as constitutive laws. We also have the incompressibility constraint for the fluid, and interface conditions such as continuity of traction stresses, elastic displacements, and velocities. We make the assumption that the size of the materials pores (the porescale) is comparable with the distance between the adjacent subphases. This length is then taken to be much smaller than the size of the whole domain (the macroscale). We call the ratio of the length scales our scale separation parameter. This allows us to decouple the spatial scales, then by embracing the asymptotic homogenization technique we derive the new macroscale model. The asymptotic homogenization technique applies the assumption that all fields in the fluid-structure interaction problem can be written as a power series of the scale separation parameter and then performing a multiple scale expansion we can derive the cell problems that determine the model coefficients. The system of partial differential equations that arises from applying the technique is of poroelastic-type. The coefficients of the model encode the properties of the microstructure, and can be computed by solving appropriate cell problems which reflect the complexity of the underlying material microstructure. The macroscale model comprises the balance of linear momentum and the conservation of mass equatioṅ where we have that ∇ is the macroscale gradient operator, (we note that with the subscript this would be the microscale gradient operator), LMRP Ef f is the stress tensor (the superscript LMRP is used to show that this is the stress that arises specifically from this model), (0) is the macroscale pressure, is the symmetric part of the macroscale gradient operator,̇( 0) is the leading order solid velocity, is the average fluid velocity, and are the resulting Biot's modulus and tensor associated with the system respectively. The conservation of mass equation relates changes in the fluid pressure to changes in the fluid and solid volumes. The macroscale model also comprises Darcy's law where ⟨ ⟩ f is the hydraulic conductivity tensor, and the constitutive law L. Miller and R. Penta where C for = I, II is the elasticity tensor for the inclusion and matrix respectively and the fourth rank tensors M I , M II are to be computed by solving the microscale cell problems that will be discussed in the next section. We can define the effective elasticity tensorC LMRP as We should note that our effective elasticity tensor possesses tetragonal symmetry, that is, possessing six distinct elastic entries. The reason for this is that, while if we began with the orthotropic elasticity tensor, we would have nine different elastic entries, as the geometry we selected is a cube with embedded cylinders with circular basis, then the 1 and 2 directions are equivalent, and hence the reduction to six independent elastic parameters. We therefore have that the behaviour of the poroelastic composite material can be fully described by the effective elasticity tensorC LMRP , the hydraulic conductivity ⟨ ⟩ f , the Biot's tensor of coefficients and the Biot's coefficient . We have that these macroscale coefficients read where the fourth rank tensors M I , M II and the second rank tensor II are to be computed by solving the microscale cell problems that will be discussed in the next section and is the porosity. We note that the notation ⟨ ⟩ is a cell average defined as where is a general field in our system and | | is the volume of the domain and the integration is taken over the porescale. Now we wish to consider the governing equations for our alternative comparison setup concerning standard poroelasticity with an elastic inclusion (SP). That is, the governing equations for a poroelastic material containing an elastic inclusion derived via the asymptotic homogenization technique. This standard poroelastic setup is merely created to act as a comparison highlighting how to approach a material that possesses a microstructure comprising a matrix, embedded elastic subphases and fluid filled pores using the computational models already available in the literature (assuming we had not yet created the LMRP model). This approach uses a combination of standard poroelasticity, together with elastic composites, however does not allow for the multiple elastic and fluid phases to be considered all at the same scale. Hence, this firstly justifies the introduction of the LMRP model, and also shows that calculating the elastic parameters via this standard poroelastic approach is not always entirely appropriate. We emphasize that this work focuses on the comparison between the drained elasticity tensors computed via the LMRP and SP approaches, respectively, see also Remark 1.
We can determine that this microstructure is a limit case of the double poroelasticity model (Miller and Penta, 2021a) when there is no fluid in the inclusion. This is also the geometry considered in Chen et al. (2020) and Royer et al. (2019). This means that we are considering a linear elastic problem where the interactions take place between a matrix (that is porous at a finer scale), PM , and an embedded linear elastic inclusion I . An alternative setup in the context of hierarchical models for fluid-saturated elastic media derived by asymptotic homogenization consists in considering the coupling between a poroelastic phase (arising from a fine scale coupling between an elastic phase and a fluid phase) and a fluid phase, as done in Rohan et al. (2015). In this latter work the authors are considering the coupling between a poroelastic phase and a fluid phase at the mesoscale, rather than the coupling between a poroelastic phase and an elastic inclusion.
We should note that the elastic inclusion I is the same in both model setups, that is, it possesses exactly the same elastic properties and volume fraction, and it is only the scale at which the matrix is porous that varies between this setup and the LMRP model. The comparison between the LMRP and SP structures is shown in Fig. 1. The interface between the inclusion and the porous matrix is denoted PM . We assume that the distance between the embedded subphases (the microscale) is small compared with the size of the whole domain (the macroscale). By enforcing this scale separation we can decouple the spatial scales and derive the effective governing equations for the poroelastic material with an elastic inclusion. The governing equations are those presented in the appendix of Miller and Penta (2021a) in the limit of only fluid in one phase. The stress balance is given by with the constitutive law wherê( 0) is the macroscale pressure,̂( 0) is the leading order elastic displacement and We have used the superscript SP to denote that this approach compares with Standard Poroelasticity. The second rank tensors SP I , SP II are to be computed by solving microscale cell problems that encode details of the materials microstructure and the second rank tensor M is the Biot's tensor of coefficients that arises from the homogenization of the porous matrix at the finer scale. We can define the effective elasticity tensorC SP as Where we have that C PM is the effective elasticity tensor that arises from carrying out the asymptotic homogenization technique on the porous matrix, C SP I is the elasticity tensor for the inclusion in this model setup, and is equal to C I from the LMRP model. The fourth rank tensors M SP I , M SP II are to be computed by solving the microscale cell problems that will be discussed in the next section. We should note that our effective elasticity tensor possesses tetragonal symmetry for the same reasons geometrical reasons as described forC LMRP previously.
The macroscale model also comprises the conservation of mass equation given bẏ̂( wherė̂( 0) is the leading order solid velocity,  and̃are the resulting Biot's modulus and tensor associated with the system respectively and are given bỹ and the ef f is given as the final macroscale equation (Darcy's law) where  PM is the Biot's modulus of the porous matrix, and the second rank tensor̃is a modified hydraulic conductivity tensor that accounts for the differences in hydraulic conductivities at different points in the microstructure.

Remark 1 (Undrained Effective Elasticity Tensors).
We have presented a brief summary of the LMRP model  and the comparative setup SP that combines standard poroelasticity with elastic composites. We have presented the whole model for each of these setups for the sake of completeness however, the simulations focus on only computing the effective drained elasticity tensor. However, it would be possible to also compute the equivalent tensors for the undrained case. To do this first for the LMRP model we would assume that we had a static fluid filling phase and the undrained elasticity tensor could be obtained using a static (2) in (4). This gives the undrained effective elasticity tensor In exactly the same way we would assume a static (12) and use in (9) to obtain the undrained elasticity tensor for the standard poroelastic with elastic inclusion approach.
By carrying out further simulations (which are beyond the scope of this particular work), we could also compute this undrained effective elasticity tensor for each of these setups. We also note that despite it not being a focus of this particular work that all the poroelastic coefficients of both model setups are able to be obtained. The additional problems which are not presented here and that would need to be solved to compute the Biot's modulus and tensor of coefficients and hydraulic conductivity can be found in Miller and Penta (2020).

Computational setup
Within this section we consider the numerical setups and describe the 3D cell problems required to compute the effective elasticity tensor for both the LMRP model and the Standard poroelasticity with an elastic inclusion. We then also provide the 2D reduction for each of the cell problems by simplifying the framework for a 2D reduction of a poroelastic composite with orthotropic elastic phases that is provided in detail in the Appendix A.1.
Here we summarize the specific equations for each of the 9 particular cell problems that we present in detail in the following 4 subsections. For the 3D LMRP model cell problem see (19)-(23) and for the 2D equivalent of the LMRP model see (140)-(144) (in-plane) and (50)-(54) (anti-plane). For the standard poroelastic with inclusion setup we have two steps, first the standard cell problem for poroelasticity in 3D is given by (73)-(74) and the 2D equivalent of this is given by (141) and (144) (in-plane) and (95) and (96) (anti-plane). We then have the second step where the 3D elastic composite problem is given by (77)-(80) and the 2D equivalent of this is given by (140)-(143) (in-plane) and (105)-(108) (anti-plane).

3D cell problems for LMRP model
We are able to compute the fourth rank effective elasticity tensor C LMRP for the LMRP model and by using its components calculate the two Young's moduli and two shear moduli corresponding to our model. The effective elasticity tensor is given bỹ We can see that this comprises the fourth rank tensor M , where = I, II, which can be defined as We can then write the cell problems for third rank tensors I and II , found in Miller and Penta (2020), with corresponding components I and II as The solutions to the problem (19)-(23), ( I ) and ( II ), are found by solving six elastic-type cell problems by fixing the couple of indices ( , ). By doing this the ( I ) and ( II ) that appear in (19)-(23) represent a strain. Then for every fixed couple (k, l) we have a linear elastic problem which has the following forces driving each of the six cell problems (19)-(23) which depend on the jump in the elastic constants between the matrix and the subphase and on the geometry of the subphase, that is encoded in the normal III to the interface between the matrix and inclusion III where III is the unit outward normal corresponding to the interface III . In order to solve (19)-(23) we also have interface conditions between the matrix and the fluid, II . We are still fixing every couple ( , ) to find the following forces that account for the difference between the elastic matrix and the void where the fluid has been removed since we are computing the drained coefficients. The normal II encodes the geometry of the voids as it is the normal to the interface II . The forces therefore are where II is the unit outward normal corresponding to the interface II . We assume that the fourth rank elasticity tensors I and II are isotropic at the porescale. That is, We can then use (36) and (37) in the interface loads to determine the forces III on III and the forces II on II . We note that for the forces the superscript given refers to the interface on which the force is applied, we will use this convention throughout this work. Firstly the forces on the matrix inclusion interface III = 11 III = ( II − I ) III + 2( II − I ) III where we have used II 1 , II 2 and II 3 to mean the components of the unit vector normal to the interface II and we have used the standard unit vectors in the Cartesian coordinate system 1 , 2 and 3 .
The cell problem (19)-(23) is a three dimensional problem. Our interface conditions (38)-(43) and (44)-(49) are also 3D. For each pair of boundary loads given in (38)-(43) and (44)-(49) we compute a corresponding numerical solution of the elastic-type problem (19)-(23). This can be done using the finite element software Comsol Multiphysics employing its Structural Mechanics Module. We wish to perform 2D simulations so we must reduce our cell problems.

2D cell problems for LMRP model
We now wish to consider the 2D reduction of the cell problems. The geometry of our periodic cell is such that we have a cube with a cylindrical elastic inclusion extending in the 3 direction from the top to the bottom of the cell as well as 4 cylindrical voids placed in each of the corners of the cube also extending from the top to the bottom of the cell in the 3 direction. This means that at every cross-section in the 3 direction gives a square with a circular inclusion and 4 circular voids. This geometry is shown in the 2D sketch Fig. 2.
We can consider the 2D reduction presented in the appendix and assume that both the inclusion and the matrix are isotropic materials. This means that we have that 1111 = 2222 = 3333 = + 2 , where the notation I , II are the corresponding 2D slices of I and II . The interfaces in 2D are represented as I ∩ II for the matrix inclusion interface and f for the fluid-matrix interface. We also note that our normals to the interfaces I ∩ II and f are the same normals (but with only two components) as in the 3D case as these were normals to cylindrical surfaces so are still the normals to the curved surfaces of the circular voids and inclusion. The solutions to the problem (50)-(54) are found by solving the two anti-plane problems by fixing the couple (k, l) = (1, 3) = (3, 1), (2, 3) = (3, 2). Then for every fixed couple we have the Poisson problem with the following two pairs of forces on the two different interfaces respectively. The first force we require to solve (50)-(54)

L. Miller and R. Penta
We now need to consider the in-plane problems. Using (140)-(144) from the appendix we can consider the forces required on both the matrix inclusion interface I ∩ II and on the matrix fluid interface and f . In the LMRP model we assume that both elastic phases are isotropic. Each interface load is a vector with two components, due to = 1, 2. We fix the couple (k,l) and obtain on the matrix inclusion interface To solve the problem we also need the forces on the fluid matrix interface f . We have that

3D cell problems for standard poroelasticity with elastic inclusion
Here we wish to compute the fourth rank effective elasticity tensor C SP for the poroelastic material with elastic inclusion and by using its components calculate the two Young's moduli and two shear moduli corresponding to this setup. This model requires two steps. We begin by finding the effective elasticity tensor for a poroelastic material and we use components of that tensor as the parameters for the matrix PM at the next scale. The problem we consider is for a porous matrix M with fluid flowing in the pores, shown in the zoomed in area of Fig. 1, and we wish to find the effective elasticity tensorC PM (porous matrix). That is where C M is the elasticity tensor for the elastic matrix and the fourth rank tensor M M is defined as We have that the third rank tensor solves the following cell problem. We have where is the unit normal pointing into the fluid. This problem can be solved as done in Dehghani et al. (2018). The solution of the problem, which is the fourth rank tensor M M , can be obtained by solving six elastic-type cell problems by fixing the couple of indices in the component wise representation of problem (73)-(74). This allows each component of M M to be interpreted as a strain and this means that for each couple of indices that are fixed we have a linear elastic problem with inhomogeneous Neumann interface conditions. The component wise cell problem and the corresponding interface conditions can be found in Dehghani et al. (2018). The second problem within this setup that we consider is for a composite comprising a matrix that is poroelastic at a finer scale, with parameters supplied from the effective elasticity tensorC PM , and an isotropic elastic inclusion. We begin by formulating the 3D problems.
Here we will be obtaining the effective elasticity tensor,C SP , for our standard poroelastic material with inclusion. That is We can see that this comprises the fourth rank tensor M SP , where = I, II, which can be defined as We can then write the cell problems for third rank tensors I and II with corresponding components I and II as where PM is the unit normal to the interface PM . We also have introduced the notation I SP , this is the component representation of C SP I and the notation PM which is the component wise representation ofC PM . The solutions to the problem (77)-(80), ( I ) and ( II ), are found by solving six elastic-type cell problems by fixing the couple of indices ( , ). By doing this the ( I ) and ( II ) that appear in (19) and (23) In this case we assume that the fourth rank elasticity tensor I SP is isotropic at the microscale and PM has components calculated by solving the poroelastic problem (73)-(74) above. That is and where we have used Voigt notation to represent the components of the tensor C PM . We should note that this effective elasticity tensor possesses tetragonal symmetry, that is, possessing six distinct elastic entries. The reason for the six entries is due to our choice of geometry (see Fig. 1 cube with cylindrical voids) which means that our 1 and 2 directions are equivalent hence the reduction from the nine entries in the orthotropic case to the six entries we have here. We can then use (87) and (88) where we have used PM 1 , PM 2 and PM 3 to mean the components of the unit vector normal to the interface PM and we have used the standard unit vectors in the Cartesian coordinate system 1 , 2 and 3 . This current setup is for solving the 3D problem. We are again able to perform a reduction so that we instead study the 2D problem.

2D cell problems for standard poroelasticity with elastic inclusion
We now wish to consider the 2D reduction of the cell problems that have been presented in the previous subsection. Within this model set up we have two different cell problems to solve. The first is the cell problem for a porous matrix. The geometry of our periodic cell in this case is such that we have a cube with 4 cylindrical voids placed in each of the corners of the cube extending from the top to the bottom of the cell in the 3 direction. This means that at every cross-section in the 3 direction gives a square with 4 circular voids (see Fig. 3).
We can consider the 2D reduction presented in the appendix where we assume that there is only the matrix M and the voids and assuming the matrix is isotropic. This means that we have that 1111 where the notation M is the corresponding 2D slice of M and the interface f ∩ M is the 2D projection of M . We also note that the normal to the interface is the same normal in 3D as in 2D since it is the outward normal to a cylinder in 3D and is the same to the circle in the 2D slices but with only components 1 and 2 .
The solutions to the problem (95) and (96) where the forces can be written using the assumption of isotropy of both phases as We now need to consider the in-plane problems. We have the problem (140)-(144) from the appendix, however for the case of a porous matrix we only require (140) and (144). We can consider the interface loads on the matrix fluid interface (146). We assume that the matrix material is isotropic. Each interface load is a vector with two components, due to = 1, 2. We fix the couple (k,l) and obtain The second cell problem is for an elastic composite that comprises the porous matrix and the elastic inclusion. The geometry of our periodic cell in this case is a cube with an embedded cylinder placed in the centre extending from the top to the bottom of the cell in the 3 direction. This means that at every cross-section in the 3 direction gives a square with an embedded circle in the centre (see Fig. 4).
We can consider the 2D reduction presented in the appendix in the case where there is no void. We assume that the inclusion is isotropic where the forces I ∩ PM on the interface I ∩ PM can be written using the assumption of isotropy of the inclusion and using the entries of C PM for the matrix We now need to consider the in-plane problems. Using (140)-(143) from the appendix where we assume that the inclusion is isotropic and we use the values of C PM . Each interface load is a vector with two components, due to = 1, 2. We fix the couple (k,l) and obtain on the matrix inclusion interface

Applications & results
Within this section we present the results of solving the 2D cell problems that were presented for each of the two different setups described in Section 2 using COMSOL Multiphysics 4.3. The solution is computed by following the procedures detailed in Section 3. We present the difference in the elastic parameters, Young's and Shear moduli, for the two different model setups. We interpret under which scenarios each model should be used.

2D simulation results
Within this section we focus on long fibres and therefore can make use of the 2D cell problems setup illustrated in the previous Sections 3.2 and 3.4. We set up our problems on a unit square cell with the circular elastic subphase accounting for 20% volume fraction and the porosity varying from 2%-30% divided among the four cylinders.
We solve the cell problems using the following parameters. For the LMRP model we have the matrix II with Poisson ratio 0.4 and Young's modulus 80 kPa, we have the elastic inclusion I with volume fraction 20% with Poisson ratio 0.49 and Young's modulus 35 kPa.
For the standard poroelastic material with elastic inclusion we have two steps the first is the porous matrix problem where we have the matrix M with Poisson ratio 0.4 and Young's modulus 80 kPa, and for the second step we have the problem between the inclusion and the porous matrix, where we have the matrix informed by the results from the porous matrix simulations and the elastic inclusion I with volume fraction 20% with Poisson ratio 0.49 and Young's modulus 35 kPa. This keeps the parameters consistent between the two model setups.
We begin by considering the comparison of the two Young's moduli 1 and 3 for the LMRP model and the standard poroelasticity with inclusion setup. Using the components of the effective elasticity tensors that we compute for each of the models at varying porosities we have the formulas for the Young's moduli given by where the superscript = LMRP, SP determines which model we are using. We have plotted the results of these Young's moduli with a range of porosities from 2%-30% in the Figs. 5(a) and 5(b). By considering Fig. 5(a) we can see that the Young's modulus 1 (transverse) decreases with increasing porosity. We can also see that the transverse Young's modulus is lower for the LMRP model. For this reason we have plotted the difference between the LMRP model and the standard poroelastic with inclusion setup in Fig. 5(b). From this plot we can see that at approx 5% porosity there is already a difference of 5% between the two different model setups. This difference between the models increases to approx 22% at 30% porosity. This means that the influence that the porosity has on the overall material stiffness, when directly being considered with the other phases, is considerably more prominent than when considering the porosity at the finer scale.
We can now consider Fig. 6(a) and we can see that the Young's modulus 3 (axial) also decreases with increasing porosity. Similarly to the transverse Young's modulus, we can see that the axial Young's modulus is also lower for the LMRP model. We have again plotted the absolute difference between the LMRP model and the standard poroelastic setup in Fig. 6(b). From this plot we can see that at approx 10% porosity there is already a difference of 2.5% between the two different model setups, rising to almost 10% at 30% porosity. This again means that the influence that the porosity has on the overall material stiffness when directly being considered with the other phases is considerably more prominent than when considering the porosity at the finer scale. However, the effect is slightly less prominent on the axial Young's modulus than the transverse one.
The other two elastic parameters that we compare for the two different model setups are the shear moduli 44 and 66 . These parameters are taken directly from the computed effective elasticity tensor for each of the models. We have plotted the comparison of the shear moduli over a range of porosities from 2%-30%. This is shown in Figs. 7(a), 7(b), 8(a) and 8(b).
In Fig. 7(a) we see that the shear 44 decreases with increasing porosity and that the LMRP model decreases more than the standard poroelastic type setup. For 44 the force is being applied in the axial direction which is the direction in which the inclusion and the voids elongate. This means that when the force is applied the material deforms and the voids flatten out, the voids just make it softer allowing for the decrease in shear with increasing porosity. The standard poroelastic type model has less of a decrease in shear due to the fact that is does not have the voids present at this scale but accounts for the increasing porosity of the matrix at a finer scale (where the voids are present). We also see that the difference between the two models is increasing with increasing porosity. For these reasons we wish to consider the absolute L. Miller and R. Penta    difference between the two model setups, and this has been done in Fig. 7(b). We can see that for 5% porosity we have approximately a 2% difference between the two model setups. This increases to a 16% discrepancy for a 30% porosity. This means that when the porosity exceeds 5% it is more useful and accurate to use the LMRP model to describe the behaviour of the material parameter 44 . For porosities below 5% then the standard poroelastic setup can be realistically used.
In Fig. 8(a) we can see that there is very little difference between the shear for both models when the porosity is less than 20%, however after this point the difference becomes more pronounced. This is confirmed by the absolute difference plot, see Fig. 8(b), where we can see that up to 20% porosity the discrepancy between the models does not exceed 2%, however after this point it reaches 13% when porosity is 30%. For 66 the force is being applied in the -direction (transverse). Therefore for the LMRP model the force is being applied taking a cross section of structure which contains the voids and inclusions. At higher porosities this makes the material weaker as the larger voids deform and hence the larger decrease in shear compared with 44 . Up until approx 20% porosity the difference in the two models is negligible (<2%), this is explained by a critical level of porosity where the scale at which the porosity is being considered becomes important with this direction of shear. When the porosity is low, the pores in the LMRP structure are small and do not influence the shear. When the porosity is higher the pores are much larger so the distance between the voids and the inclusion becomes less and then the difference this make to the shear value become apparent. This critical level of porosity where we begin to see the difference between the models is also influenced by the length of the embedded fibre. The shorter the fibre the more pronounced the difference between the models is at a lower porosity. We can see this is the case in the following section Section 5 where we consider a variety of fibre lengths.
We can now summarize the findings and explain how they should be interpreted. When a material possesses a microstructure comprising a matrix, embedded elastic subphases and fluid filled pores, then using the models that are currently available in the literature (excluding the LMRP model) offers two approaches. These are either make the assumption that the matrix is homogeneous (ignore the subphases) which is in general not true, or carry out a two-step process, such as the comparative standard poroelasticity approach (SP) illustrated in this manuscript. This second approach means that even though we are considering the three phases, we are not considering them all at the same scale, which is what the intended application actually possesses as a microstructure. This means that the estimations of the parameters Table 1 Threshold porosities for when model discrepancy exceeds 2% (long fibres). cannot be fully reliable. We therefore developed the LMRP model to remove this issue. With this novel model it is possible to account for multiple elastic phases all with different properties as well as the fluid all at the same scale. It is for this reason that the LMRP model results are more accurate for these types of applications as we are truly capturing the correct microstructure by using our model. These are the results that we have shown in the previous plots, where we assume 2% discrepancy as a threshold for determining when the LMRP model should be used. The simulations presented here are solely related to the drained elastic coefficients. While we are not including the explicit contribution of the fluid, it is completely possible to calculate the undrained model coefficients, see Remark 1, as well as the hydraulic conductivity tensor (this latter being the sole cell problem to be computed in the fluid cell portion), as done for instance in Dehghani et al. (2018).

Model
After considering all four of our elastic parameters we can now summarize our findings in Table 1, where we have set the critical model discrepancy percentage to be 2%. We find that at porosities exceeding 20% then the LMRP model is the more effective at determining the true elastic parameters for our given geometry. However, even at lower porosities (<20%) the LMRP model is more effective at determining 1 , 3 and 44 and is equally effective at determining 66 as the standard poroelastic type model. We should note that there may be different choices of geometry where the results of these simulations become more or less pronounced and of course the specific application and elastic properties of the phases will also play a role. We can however say that the LMRP model can efficiently incorporate multiple phases all at the same scale and should be considered a valid choice when modelling poroelastic materials with inhomogeneous microstructures.

3D Application -short fibres
Within this section we discuss the comparison of elastic parameters for short fibre models. This can be an interesting application since the difference in the elastic parameters varies depending on whether or not the subphases are connected between the cells. It can also be L. Miller and R. Penta useful for biological applications where specific cells are inclusions and not subphases. Here we must use a 3D framework since our embedded inclusion does not run from the top of our periodic cell to the bottom but is in fact fully embedded (i.e. a short fibre), therefore the fibres do not connect between cells. In this case we cannot use the 2D simulations since some of the slices in the 3 direction would not include the fibre. We therefore perform 3D simulations to solve the 3D cell problems (19)-(23) for LMRP, and for the comparative standard poroelastic model (73)-(74) and (77)-(80), presented in the previous Sections 3.1 and 3.3. We can justify that our 3D simulations are accurate by considering the absolute error plots between the 2D and 3D long fibres simulations presented in Appendix A.3, where we see that there is less than 1% error for all parameters considered.
We begin by considering the difference in the Young's Moduli between the two model setups for the short fibres. In this case we choose that the fibres are length 0.8 out of a length one cube and are placed in the centre. For further details about how the model is setup in COMSOL Multiphysics see Appendix A.2. In Fig. 9 we have carried out the 3D simulations to give a comparison between our two different computational setups for both the Young's moduli 1 and 3 for short fibres. We see that both 1 and 3 are lower for the LMRP model. This can be explained by the fact that using the LMRP model explicitly considers the fluid contribution at the porescale along with the matrix and the inclusion. In particular, the LMRP model fully considers the influence of the porosity (which appears as a void in the microscale geometry) thus leading to a lower value of the stiffness moduli.
To confirm our deductions we have plotted the absolute difference between the model setups for both of the Young's moduli. In Fig. 10 we can see that the discrepancy between the two models for 1 increases from 1.5% to 19% with increasing porosity. This means that for very low porosities then both LMRP and Standard poroelastic type model both produce similar results, but for materials with higher porosities the LMRP model will provide the most accurate representation of the parameters. Similarly for 3 the absolute difference increases from 0.5% to 9% with increasing porosity. This means that for low porosities (less than 10%) then both LMRP and Standard poroelastic type model both produce similar results. When the materials however, have higher porosities then again the LMRP model will provide the most accurate representation of the parameters.
We also wish to consider the difference in the shear moduli. In Figs. 11(a) and 12(a) we have carried out the 3D simulations to give a comparison between our two different computational setups for both the Shear moduli 44 and 66 for short fibres. Here we are again using the 3D framework since our embedded inclusion does not run from the top of our periodic cell to the bottom but is in fact fully embedded (i.e. a short fibre).
In Fig. 11(a) we see that for shear 44 that the LMRP model decreases more than the standard poroelastic type model and that the difference between the two models is increasing with increasing porosity. The LMRP model considers the voids and the two elastic phases at the same scale and this contributes to the greater decrease in shear. The standard poroelastic material with inclusion has less of a decrease in shear due to the fact that is does not have the voids present at this scale but accounts for the increasing porosity of the matrix at a finer scale. We also plot the absolute difference between the two models. In Fig. 11(b) we see that for very low porosities that the difference between the two models is less than 2%, however, for increasing porosities the discrepancy reaches 18%. This means that if we have a material with a porosity greater than 5% then the standard poroelastic approach will not capture the true elastic parameter and it would be much more appropriate to use the LMRP model. Yet still for materials with very low porosities both models will produce similar results.
For short fibre 66 the physical description of the deformation is the same as the long fibre case. In Fig. 12(a) we can see that up until approx 15% porosity the difference in the two models is negligible, this could be explained by a critical level of porosity where the scale at which the porosity is being considered becomes important with this direction of shear. This is because at low porosities the size of the voids is very small in the LMRP model and so therefore do not influence the shear more than the pores at a finer scale. However, once the porosity exceeds 15% then the pores are large enough to influence the shear. This is confirmed by Fig. 12(b) where we can see that for porosities up to 15% that the difference between the models is less than 2%. This means that up to this level of porosity the standard poroelastic approach or the LMRP model are capturing the behaviour similarly. However for porosities greater than 15% then the discrepancy between the models approaches 16% this highlights than in this case it would be more realistic to use the LMRP model. This critical level of porosity where we begin to see the difference between the models is influenced by the length of the embedded fibre. The shorter the fibre the difference between the models becomes more pronounced at a lower porosity. This is due to us keeping the volume of the fibre consistent even though the length is decreasing, so the fibre becomes shorter but thicker meaning the distance between the voids and the fibres is smaller in the shorter fibre models and this means that the voids can have an L. Miller and R. Penta Fig. 11. Results of short fibre Shear Modulus 44 simulations. influence on the shear at a lower porosity. Note that this idea can be further enforced by the plot of 66 Fig. 8(a) where we see that up until approx 20% porosity the difference in the two models is negligible. We further confirm this idea by performing the simulations for a range of fibres lengths from 0.6-1, results shown in Table 3 and Fig. 13.
From considering all four elastic parameters we can now summarize or findings for the short fibre simulations in Table 2. We find that at porosities exceeding 15% then the LMRP model can be more effective at determining the true elastic parameters when multiple phases interact at the same scale. However, even at low porosities (<15%) the LMRP model is more effective at determining 1 , 3 and 44 and is equally effective at determining 66 as the standard poroelastic type model. We also note that the range of the discrepancy between the models with increasing porosity is higher for the short fibre than the long fibre setting. This means we can conclude that the LMRP model is an efficient choice for modelling poroelastic materials with inhomogeneous microstructures.

Conclusions and future perspectives
Within this work we have created a robust computational platform that has allowed for a valid comparison between the LMRP model for  poroelastic composites and an approach that uses standard poroelasticity with elastic inclusions. We describe our platform as robust due to the wide range of situations where it can be used for computations. That is, the platform can be altered for a variety of geometries including short fibres, various directions of fluid flow, a variety of differently shaped inclusions and a wide range of constitutive properties of the constituents. We investigated a variety of elastic parameters obtained by solving the 2D cell problems, which were derived in this work, to determine under which circumstances the approach by the LMRP model is most appropriate to be used. We begin our analysis by providing a summary of the macroscale LMRP model for poroelastic composites that has been derived via asymptotic homogenization in Miller and Penta (2020). We then present the 3D cell problem (19)-(23) that is to be solved to obtain the model coefficients such as the effective elasticity tensor. From this problem we write down explicitly the boundary loads that are required to solve the problem numerically. In order to be computationally less expensive we carry out, and present, the reduction of the cell problems to 2D, where we again present the appropriate boundary loads. Since the aim of this work is to have a valid comparison that uses a standard poroelastic approach we then derive the comparison model setup using asymptotic homogenization. The comparison standard poroelastic with inclusion model has two steps, firstly the porous matrix and then a composite between the porous matrix and an elastic inclusion. This means that we have two 3D cell problems (73)-(74) and (77)-(80) that are to be solved to provide the model coefficients. We again write down explicitly the boundary loads that are required to solve these two problems numerically. Again in keeping with the comparison, we carry out and present the reduction of both of these cell problems to 2D, again presenting the appropriate boundary loads. We are then able to solve numerically the cell problems for a simplified geometry where we have unidirectional flow and only one fibre direction. We then plot and compare relevant elastic parameters, i.e. Young's and shear moduli.
We then present an example of our computational platform solving the 3D cell problems, justified by the error plots between the 3D and 2D simulations showing less than 1% error.
The results of our numerical simulations show that when investigating a poroelastic composite material with porosity exceeding 5%, then the discrepancy between the LMRP model (that truly considers multiple phases at the same scale) versus a standard poroelastic approach when computing the Young's moduli 1 and 3 and the shear 44 exceeds 2% and increases dramatically with increasing porosity. We also find that when the porosity exceeds 20% the discrepancy between the two model setups goes beyond the 2% threshold for the shear 66 . We find that for materials with less than 5% porosity a standard poroelastic approach or the LMRP model produce approximately the same results. These results mean that the LMRP model should be considered as a valid alternative instead of the models in the literature (i.e. standard poroelastic approach) when modelling a material with an inhomogeneous microstructure.
Our current modelling approach is open for improvement. Within this work we have only focused on the parameters of the elastic matrix. It is also of interest to investigate the fluid flow and solve the cell problem in order to obtain the hydraulic conductivity tensor for the material.
The simulations in this work have also only been carried out for a simplified geometry (unidirectional flow and fibre direction) so as to show the difference in the two models in the simplest possible case. It would however be possible, due to the robustness of the 3D computational platform highlighted by the 3D example, to consider a much more complex geometry consisting of additional fibres and also fibre angles and fluid flow in many directions.
We should also note that in this work we are focusing solely on solving the microscale cell problems that determine the model coefficients. It would however, be possible to solve the complete macroscale model that is presented in Section 2.
The derivation and numerical simulations have focused on a general set of parameters since we are aiming just to capture the effects of the LMRP model and the settings where it is most applicable. The next step would be to apply the LMRP model to a realistic set of parameters and geometry and this will allow for model validation by comparing with experimental data, as done for example by Chen et al. (2020). Another important aspect which is to be considered is the model validation in terms of how well it converges to the actual behaviour of the physical system whenever scales become more and more separated. This is indeed a problematic issue primarily due to the required computational cost, although one next natural step is indeed the development of direct numerical simulations (for example performed on reference heterogeneous geometries in two-dimensions). This approach (which is carried out for example in Cruz-González et al. (2022) in the context of three-scale asymptotic homogenization for a one-dimensional example) would increase the reliability of these results, and in general, of any model derived via homogenization techniques, and will become more and more realistic following advances in computational resources available. This specific validation will also better elucidate the role of the heterogeneities as a discriminant in determining which homogenized model better represents the actual physical system at hand.
Finally, the simulations could be extended to investigate changing the volume fraction and/or geometry of the inclusion. This could have many relevant biological applications, e.g. to myocardial infarction, where the cardiac myocytes die and become replaced by fibrous collagen matrix. In this case we could perform a parametric analysis where the inclusion volume fraction and geometry are changed to simulate the loss of myocytes and the replacement with fibrous collagen matrix scar tissue. We note that of course our model here uses linear elasticity and the heart is of course nonlinear. However, we can obtain results by using a piecewise linear approach as done in Hu et al. (2003a), Hu et al. (2003b). By doing this we can approximate the nonlinear behaviour using simple, computationally cheap simulations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
No data was used for the research described in the article.

A.2. Numerical simulations and meshing
Within this section we aim to give an overview of the steps carried out in COMSOL Multiphysics to compute the results presented in this work. To do this in the most clear and useful way we will begin with the set up for the 3D simulations platform, and then explain the simplified 2D simulations which have results presented in Figs. 5(a)-8(a).
We can now discuss the 3D simulations that are applicable to long fibres embedded in the matrix with the fluid cylinders. We can then mention the small modification for the case of short fibres, with results shown in Figs. 9-12(a). For the LMRP model the cell problem (19)-(23) is a three dimensional problem. Our interface conditions (38)-(43) and (44)-(49) are also 3D. For each pair of boundary loads given in (38)-(43) and (44)-(49) we compute a corresponding numerical solution of the elastic-type problem (19)-(23). This can be done using the finite element software COMSOL Multiphysics employing its Structural Mechanics Module.
We use this software to compute the 36 entries of the tensors M I and M II . Then once we have these results they can be used in (5) to obtain the entries ofC LMRP . Once we have the complete tensorC LMRP then we can use the components in the formulas for our elastic moduli 1 and 3 and take the shears directly from the tensor.
We will now give some details of how this process in COMSOL is carried out. The finite element software creates a mesh for the periodic cell . It does so in such a way that it creates a surface mesh for the interfaces between the phases and around the voids. It then extends the surface mesh into a three-dimensional one for the entire periodic cell . By using this method it allows for both interface conditions described by boundary pairs and interfaces on the drained fluid voids. This is beneficial since it allows for a particularly refined mesh on the interfaces (where the important physics takes place) and there surrounding area which gets gradually coarser the further away from the subphase and void interfaces that we are.
We know from the cell problem (19)-(23) that the stress-jump condition on the matrix-subphase interface as well as the condition on the interface between the matrix and fluid are the driving forces for the solution of the cell problems. This means that we require a sufficiently fine mesh on these interfaces to ensure we obtain an accurate numerical solution. Therefore, to capture these areas where the important physics is taking place, our mesh in is set to be much more refined around the boundary pairs and voids representing the interfaces than in the remainder of the domain further away from these interfaces. It is important we mention that we can use a sequence of increasingly refined meshes of . These meshes are predefined by COMSOL Multiphysics and range in refinement from extremely coarse to extremely fine.
In cell problem (19)-(23) the stress balance equations (19)-(20) are simple since we have zero volume forces and a constant, isotropic elasticity tensor for both the matrix II and the subphase I . The stress jump and the continuity of the auxiliary tensors I and II ((21)-(22)) across the interface between the two elastic phases are encoded by conditions on each boundary pair. The condition between the matrix and the void (23) is encoded on the interface. We then impose periodic boundary conditions on the outer boundary . By using this framework the solution of the elastic-type problem will be unique up to a constant. The constant here is unimportant since it will disappear when the partial derivatives of the solution are taken to allow us to determine M I and M II . That being said, computationally we do require a unique solution of the elastic-type problems in our periodic cell. To obtain this we can add an additional constraint in COMSOL that fixes the auxiliary displacement to zero in one corner point of . This fixes the constant that is obtained. COMSOL Multiphysics uses the principle of virtual work to implement the elastic-type problem described above in weak form.
As we do not have continuity of stresses the problem for the auxiliary variables I and II is solved in the geometrical setting in L. Miller and R. Penta Fig. 14. Error between 2D and 3D simulations. the Comsol feature assembly. If we were to use the union setting then the subdomains II and I would be merged to form a simple union with continuity which is not the case in our materials. By using the assembly feature we are able to retain the boundaries for each phase of the domain, which allows for the necessary flexibility in the application of the interface conditions.
The entries of the third rank tensors I and II are numerically approximated once the six elastic-type problems (19)-(23) corresponding to the six pairs of interface loads (38)- (43) and (44)-(49) have been solved. The derivatives of the entries of I and II are linear functions that can be evaluated without additional error and therefore so can all entries of the auxiliary fourth rank tensors M I and M II . The entries of the effective elasticity tensorC LMRP are then computed using (5) by calculating the averages without additional errors. All the steps carried out such as the finite element approximations for I and II to the computation of the effective elasticity tensorC LMRP is obtained in COMSOL Multiphysics by using its integral post-processing tools.
For the short fibre simulations the same setup is used however, in this case the boundary pairs are on the full cylindrical surface (curved walls and circular ends) of the embedded subphase rather than just on the curved surface of the cylinder since it is fully embedded in the matrix. The rest of the setup and post-processing remains the same procedure.
A very similar setup can be employed for the standard poroelastic setting. Here we have created a platform comprising 3 phases (and many more could be added) but by assuming limit cases of the details just provided we can split this platform into the two steps required by the standard poroelastic approach. We first have this setup assuming the subphase does not exist which would mean it would be a porous matrix set up with only the interface between matrix and void. Secondly the setup can be modified assuming the voids do not exist, for further details on this second step see the (Penta and Gerisch, 2015) where the numerical simulations for a composite have been carried out.
In terms of the 2D problems, the setup is almost identical although this time we are using a 2D domain and the interfaces are lines not surfaces. We again do not have continuity of stresses so the problem for the auxiliary variables I and II is still solved in the geometrical setting using the COMSOL feature assembly where the interface is just a circle. We have our voids on which the forces are placed on the interfaces. We cave the corresponding periodic conditions which are applied on the external edges (lines) making the boundary. For uniqueness of solution a constraint in COMSOL is placed that fixes the auxiliary displacement to zero in one corner of the square domain. In the same way all the steps carried out such as the finite element approximations for I and II to the computation of the effective elasticity tensorC LMRP are obtained in COMSOL Multiphysics by using its integral post-processing tools. This time using surface integration rather than volume integration. L. Miller and R. Penta A.3. Error plots Within this section of the appendix we show the plots of the error we obtained between carrying out the 3D and 2D simulations for both the LMRP model and the standard poroelasticity type model. All the plots have at most a 1% error between the 3D and 2D simulations which justifies the accuracy of the 3D simulations that we have carried out and justifies the results obtained for the short fibres in Section 5 (see Figs. 14 and 15).