The role of microscale solid matrix compressibility on the mechanical behaviour of poroelastic materials European of Mechanics / A Solids

We present the macroscale three-dimensional numerical solution of anisotropic Biot ’ s poroelasticity, with co- efficients derived from a micromechanical analysis as prescribed by the asymptotic homogenisation technique. The system of partial differential equations (PDEs) is discretised by finite elements, exploiting a formal analogy with the fully coupled thermal displacement systems of PDEs implemented in the commercial software Abaqus. The robustness of our computational framework is confirmed by comparison with the well-known analytical solution of the one-dimensional Therzaghi ’ s consolidation problem. We then perform three-dimensional numerical simulations of the model in a sphere (representing a biological tissue) by applying a given constant pressure in the cavity. We investigate how the macroscale radial displacements (as well as pressures) profiles are affected by the microscale solid matrix compressibility (MSMC). Our results suggest that the role of the MSMC on the macroscale displacements becomes more and more prominent by increasing the length of the time interval during which the constant pressure is applied. As such, we suggest that parameter estimation based on tech- niques such as poroelastography (which are commonly used in the context of biological tissues, such as the brain, as well as solid tumours) should allow for a sufficiently long time in order to give a more accurate estimation of the mechanical properties of tissues.

The asymptotic homogenisation technique (see, e.g., Auriault et al. (2009) ;Bakhvalov and Panasenko (1989); Holmes (2013); Mei and Vernescu (2010); Bensoussan et al. (1978); Palencia (1980)) exploits the sharp length scale separation between the pore-scale and the average size of the domain (the macroscale) to derive the macroscale equations of poroelasticity via providing a precise prescription of the coefficients of the model (as opposed to most of the other alternative approaches), as done in Burridge and Keller (1981) and extended in Penta et al. (2014) for growing materials. These coefficients encode information on the microstructure, as they are to be computed by solving appropriate cell problems at the pore-scale.
In Dehghani et al. (2018) the authors conduct a three-dimensional computational analysis via computing the relevant macroscale poroelastic coefficients, i.e. the stiffness tensor, the hydraulic conductivity, Biot's modulus, and Biot's coefficients, as prescribed via the asymptotic homogenisation technique via solving elastic-type and Stokes'-type periodic cell problems on a connected porous structure. Their analysis, framed in the context of tumour modelling, serves as a basis to quantify the response of poroelastic materials in terms of the microstructure, and primarily focuses on the role of the matrix's porosity and compressibility.
In fact, the mechanical properties of biological tissues, such as the brain and (healthy and malignant) cell aggregates, can be estimated in vivo by means of poroelastography (Perriñez et al. (2010) Weaver et al. (2012)), elastography (Islam et al. (2018)), and viscoelastography (Kazemirad et al. (2016)). Macroscale displacements are induced in the target tissue from different sources, e.g. Valsalva manoeuvre (Mousavi et al. (2014)), blood pressure pulses (Weaver et al. (2012)), mechanical waves (Weaver et al. (2001)), etc., and are subsequently measured. On the basis of the constitutive assumptions characterising a given tissue, the effective poroelastic, elastic, or viscoelastic properties of the tissue can then be estimated. For example, in Perriñez et al. (2010) the elastic matrix is assumed to behave as a linear elastic material.
Here, we extend the analysis carried out in Dehghani et al. (2018) by solving the actual macroscale Biot's system of PDEs to provide displacements and pore pressure maps via performing a parametric study in terms of the Microscale Solid Matrix Compressibility (MSMC). By doing so, we provide ways of improving grounds for the reliability of current poroelastography protocols by quantifying the importance of the time period during which such a technique should be applied, as well as by highlighting the role of the loading conditions (i.e. transient vs. steady-state analyses).
We solve the quasi-static anisotropic Biot's system of PDEs for poroelastic media via a finite element (FE) discretisation exploiting a formal analogy with the fully-coupled thermo-mechanical system of equations, which is well implemented in the commercial software Abaqus for different types of elements. We first benchmark our numerical approach against the analytical solution of the one-dimensional Terzaghi consolidation problem by using a given set of parameters. We then solve the model by accounting explicitly for the relationship between the pore-scale and the macroscale by exploiting the results reported in Dehghani et al. (2018). In particular, we follow our approach by solving the model in a sphere with a cavity which represents a compressible, poroelastic biological tissue. The analysis is performed by varying the MSMC and by fixing the porosity and the other micromechanical parameters (such as the Young's modulus of the matrix) and by applying a constant pressure in the cavity. Furthermore, we quantify the arising model anisotropy which is also present when accounting for a simplified microstructure, invariant under permutation of the three orthogonal axes, such as the one exploited in Dehghani et al. (2018).
The remainder of this work is organised as follows: In Section 2 the macroscale governing equations and their FE formulation are provided. Section 3 provides an illustration on the numerical implementation which relies on the analogy between Biot's equations and the fullycoupled thermo-mechanical system of equations (implemented in the FE software Abaqus). In Section 4 our framework is benchmarked and the results are shown and discussed. To close with, a summary of the presented research work, as well as a discussion concerning further developments of this work, are provided in Section 5.

Essential balance relations
We now introduce the homogenized Biot's poroelastic equations as derived via asymptotic homogenisation in Burridge and Keller (1981). We refer to the quasi-static case, ignore body forces and make use of a notation which is consistent with that used in Dehghani et al. (2018).
The homogenized model is derived by considering a fluid-structure interaction problem between a linear elastic matrix and a Newtonian fluid phase flowing in the pores. Relevant fields are then expressed in terms of power series expansions in terms of the ratio between the pore-scale and the macroscale (this latter typically representing the average size of the domain). The model is then derived by closing a system of PDEs for the leading (zeroth) order fields, namely the macroscale fluid pressure p ð0Þ and displacement u ð0Þ . The porous media problem is solved for small strains, such that the small strain measure is defined as usual as where u ð0Þ denotes the solid displacement. The essential balance equations at the macroscale are the balance of linear momentum and the continuity equation, respectively the conservation of mass. Their local forms allow representation as wherein the quasi-static case is assumed and where additional body forces are neglected. Equation (3) provides us with a continuity statement that relates the change in the volume of fluid and solid to the time derivative of pore pressure. Moreover, σ is the Cauchy stress tensor, p ð0Þ describes the present pore pressure, 〈v〉 rf represents the relative pore fluid velocity, M and α are the Biot modulus and the Biot tensor, respectively. The Biot tensor is defined as α ¼α I, where α represents the Biot coefficient and I denotes the second order identity tensor. The notation _ � is used to define the time derivative of the quantity �. In order to close the overall system of equations, the following constitutive relations shall be introduced for the stresses as well as for the relative pore fluid velocity, to be specific σ ¼Ẽ : ε À p ð0Þα ; (4) where Ẽ denotes the effective elasticity tensor and where 〈W〉 f symbolises the hydraulic conductivity.
The above system refers to the anisotropic quasi-static Biot's equations by using the notation that can be found in Penta et al. (2014), where the authors obtain a more general form of the above mentioned equations from the microstructure.

Finite element formulation
In this section, the weak form is introduced by considering the contributions from both, mechanical equilibrium in the absence of mechanical body forces for a porous medium and the continuity statement for the liquid phase in the porous medium. Based on this, the model equations are transformed into discrete forms. The FE implementation is briefly discussed in Section 3. Even though the general FE framework has been introduced in the literature before, the discussion of the FE model is included in this contribution in order to show in particular the similarities between the FE formulation of the poroelasticity and thermoelasticity.
The weak form is obtained by multiplying (2) with an arbitrary test function δu that is admissible regarding the continuity and the Dirichlet boundary conditions. The result is then integrated with respect to the volume of the underlying body B . Applying the divergence theorem and using Cauchy's theorem t ¼ σ⋅n, where t represents surface tractions per unit area and where n denotes the outward normal unit vector, results in In analogy, the same procedure is used for (3) with an admissible, continuous test function δp, resulting after integration over the domain Afterwards, the integral is split into three parts and divided by the Biot modulus M so that Finally, integration by parts and the divergence theorem can be applied such that As a next step, the discretisation of the problem in space and time has to be performed. Therefore, body B is decomposed into n el finite elements and the spacial discretisation of the weak forms is represented by In addition, the test functions are discretised elementwise for each element node n en in the following manner where N denotes the corresponding set of shape functions. The contributions δu eA and δp eC represent the respective values of the test function at the nodes. With these relations at hand and by additionally defining the assembly operator A, the discretised weak form can finally be given by The same approximation is applied analogously for the physical fields of u e and p e according to (12). For the time discretisation of the pressure rate _ p ð0Þ and strain rate _ ε an implicit Backward Euler integration scheme is adopted. For the respective variable � the definition _ � ¼ ½ nþ1 � À n ��=Δt is used, where n þ 1 denotes the current time step nþ1 t and Δt ¼ nþ1 t À n t.

Finite element implementation
In this section, we discuss the implementation of the macroscale system of equations into the commercial finite element program Abaqus by considering the mechanical equilibrium statement for the porous medium together with the continuity statement for its liquid phase. To solve this problem, a coupled pore fluid diffusion/stress analysis is usually used. Abaqus provides a built-in feature for such a problem. However, it is assumed that both solid and fluid phases are intrinsically incompressible which in turn means that Biot's modulus approaches infinity. Therefore, this assumption is not appropriate for the poroelastic problem at hand. In the following, the analogy to the heat equation shall be examined by briefly summerising basic equations related to the heat equation and its computational treatment. With help of Abaqus subroutines UMAT and UMATHT the porous media material model can then be realised. This allows the implementation of the framework at hand without writing a user element subroutine UEL in Abaqus.

Analogy to heat equation
The strong form of the energy balance can be expressed as where ρ is the mass density of the material, c represents the specific heat capacity, _ θ _ θ defines the temperature, q denotes the heat flux vector and where r symbolises heat sources. The heat sources in a thermomechanical problem can be additively decomposed into the sum of external heat supply r E and into the volumetric heat generation per unit time caused by mechanical working of the material r M defined as The heat conduction is represented by the Fourier law with the thermal conductivity tensor κ.
The weak form of (15) is derived by analogy with the procedure in Section 2.2. Multiplying (15) with an admissible test function δθ and integrating over the domain results in Applying integration by parts and the divergence theorem, the aforementioned equation can be reformulated to The similarity of (9) and (19) is already obvious. As a final step, the aforementioned equation is discretised in space and time in the same manner as in Section 2.2. By dividing body B into a finite number of elements n el and discretising the test functions for each element node n en , i.e.
equation (19) can be finally assembled as Comparing the latter equation with (14), it can be summarised that the heat equation as well as the continuity equation are both partially different equations in a similar form. With this property at hand the implementation into Abaqus is possible in a straightforward manner as defined in the next section.

Implementation in Abaqus
Due to the similarity of the continuity equation and the heat equation, it is possible to use the coupled thermal-displacement analysis within the commercial FEM software Abaqus. The mechanical and thermal material behaviour with the respective subroutines UMAT and UMATHT can be adapted by the user. This way, the heat equation is "manipulated", such that the temperature degree of freedom is equal to the pore pressure and that the heat flux corresponds to the relative pore fluid velocity. Then, Abaqus can be used to implement the discretised macroscale system of equations as defined in (14) and (15). Contrary to the formulation of a user element in Abaqus via the UEL subroutine, the time for development and for testing the model is less when using a coupled temperature-displacement analysis. This also enables the usage of all standard features within Abaqus. The heat equation with the subroutine UMATHT has been manipulated due to these advantages before, see for example Ostwald et al. (2019) for a finite deformation gradient-enhanced damage model. To gain a better understanding into the implementation into Abaqus, all return values for both subroutines shall be analysed in the following paragraphs.
For a thermo-mechanically coupled user material, the general return quantities STRESS, STATEV, DDSDDE and the thermo-mechanical quantities RPL, DRPLDT, DRPLDE and DDSDDT must be evaluated within the UMAT user subroutine for the mechanical material model. These quantities are further defined in Table 1.
Within the subroutine, not only the strain tensor ε is provided as an input variable, but also the temperature θ, here respectively the pore pressure p ð0Þ , as the heat transfer equation is "abused" to solve the partial differential equation of the pore pressure (3). Thus it is possible to update the stresses σ according to the constitutive equation (4) within the routine UMAT, where STRESS is used as the return variable in Abaqus.
In addition, volumetric heat generation caused by mechanical working r M is incorporated into RPL. For a thermo-mechanically coupled problem, the dissipation and the Gough-Joule effect would be expected here.
In view of (15), the only contribution dependent on the displacement field, i.e.
is incorporated into the subroutine of the mechanical material model. Furthermore, the respective derivatives, i.e. the contributions to the Jacobian within the Newton-Raphson-based solution scheme, are necessary, as summarised in Table 1.
In addition, the thermal material model can be adapted within the user subroutine UMATHT for heat transfer, where the quantities U, DUDT, DUDG for the internal energy and FLUX, DFDT, DFDG for the heat flux have to be specified. The temperature degree of freedom θ is substituted by the pore pressure p ð0Þ as mentioned before. For the internal energy U, the expression U ¼ p ð0Þ =M is valid when comparing the respective equations (14) and (21). The heat flux q is replaced with the relative pore fluid velocity 〈v〉 rf , where the constitutive relation (5) is used for the return variable FLUX. The respective derivatives are shown in Table 2.

Applications and results
In this section we present the results obtained by solving the macroscale system of poroelastic PDEs introduced in Section 2. The solution is computed by finite elements by following the methodology depicted in Section 3. Next, in Section 4.1, we present the numerical solution of the one-dimensional Terzaghi's problem, which is successfully benchmarked by the well-known analytical solution of the problem reported, for example, in Ferronato et al. (2010). In Section 4.2 we solve the macroscale Biot's equations in three dimensions by accounting explicitly for the role of the microstructure in the coefficients of the model Dehghani et al. (2018), and highlight the importance of the proposed framework in the context of poroelastography for biological tissues. In Section 4.3 we quantify the anisotropic behaviour of the macroscale fields of interest.

Terzaghi's problem
Terzaghi's problem is a special case of quasi-static Biot's consolidation problem in which both phases are assumed to be intrinsically incompressible.
In order to verify the reliability of our computational framework, we solve the problem for a saturated column of sand under an applied load p p ¼ 1000 Paat height z ¼ 0, and compare the results with the analytical solution given in Ferronato et al. (2010), which reads where, Here, M; c; p 0 ¼ 4662:4 Pa; u 0max ¼ 6:6720 � 10 À 4 m; and L ¼ 15 m are Biot's modulus, the consolidation coefficient, the instant maximum pore pressure, instant settlement (i.e. the maximum displacement), and the height of the undrained column of sand, respectively. The consolidation coefficient c is related to the sand hydraulic conductivity K, fluid specific weight ρg, Lame's constants λ and μ, and to Biot's coefficient α and modulus M as per the work Ferronato et al. (2010).
The expressions given by (23) and (24) represent the analytical solutions for the quasi-static Biot's consolidation problem when the latter can be reduced to one dimension due to the given symmetry, as in this case, see Fig. 1. However, here we focus on the Terzaghi's limit, that is α→1and M→∞. In practice, a sufficiently high Biot's modulus is chosen, and the parameters which in this case are prescribed and do not encode microstructural properties, are given in Table 3.
In this example, drainage is only allowed at z ¼ 0 where the load is assigned, zero displacements are imposed in all directions at z ¼ 15 m, and free axial displacements on the sides, so that the consolidation/ settlement is allowed. As expected, it is necessary to refine the mesh only along the axial direction where displacements are allowed, according to the given symmetry. We have used 100 elements with a quadratic interpolation function for the displacement, and linear ones for the pore pressure in the height of the column as shown in Fig. 1.
In a consolidation analysis, the time increment is of great importance as it is a time-dependent process. In principle, the higher the number of time increments, the closer the numerical and analytical results should be. For this specific setting, 1000 time increments are sufficient for the difference between these two solutions to be negligible. Fig. 2 shows the almost perfect match between the numerical and analytical solutions at different times.

Multiscale Biot's poroelasticity for biological tissues
Biological tissues, such as brain white and grey matter and cell aggregates, are typically characterised by a poroelastic mechanical behaviour, see, e.g. Basser (1992); Franceschini et al. (2006). Here, we provide displacements, pressure, and volumetric strains radial maps by solving Biot's equations in three dimensions. We perform a multiscale analysis by explicitly accounting for the role of the microstructure in the coefficients of the model. In particular, we carry out a parametric analysis in terms of the MSMC (i.e. the microscale Poisson's ratio of the elastic matrix), which in turn affects the macroscale coefficients as reported in Table 4.
The effect of the MSMC, as opposed to porosity and elastic stiffness (see, e.g., Netti et al. (2000)) on the poroelastic behaviour of biological tissues, including the brain, has been mostly overlooked in the biophysical literature. This is possibly because the elastic phase in biological tissues is often assumed to be intrinsically incompressible, see, e.g. Netti et al. (1995) and Franceschini et al. (2006), where the latter is related explicitly to brain tissues. However, there is a remarkable variability in the reported compressibility of biological tissues (and especially the brain), so that Poisson's ratios values down to 0.35 have been reported for brain tumours, as remarked in Stewart et al. (2017). This motivates us to investigate the poroelastic macroscale response of such tissues by varying the MSMC (i.e. the matrix Poisson's ratio) in the range mentioned in Tables 4 and i.e. from 0.35 to 0.49 (the latter approximately representing incompressibility), by keeping the remaining input parameters (i.e. the solid matrix Young's modulus, fluid viscosity, and porosity) fixed.
The coefficients are computed by means of the asymptotic homogenisation closed formulas derived in Burridge and Keller (1981), as revisited for example in Penta et al. (2020) under the assumption of a periodic microstructure. These formulas comprise auxiliary variables which are to be computed by performing three-dimensional simulations at the pore-scale, as illustrated in Dehghani et al. (2018). This analysis is aimed at providing useful hints to improve the accuracy of poroelastic parameters deduced via poroelastography techniques when the latter are applied to tissues that are expected to behave as compressible materials. These include for example various types of benign and malignant cell aggregates Dehghani et al. (2018), as well as the unopened skull, see, e.g. Mousavi et al. (2014).
In Dehghani et al. (2018), the authors show that the MSMC does not affect the effective Young modulus. As such, we prescribe the Young modulus of the elastic matrix (by fixing its value to 2:57 � 10 4 Pa) in  Table 3 List of parameters used to compute both the analytical and the numerical solution of the model. Hydraulic  order to obtain, by means of the approach described in Dehghani et al. (2018), a macroscale Young modulus E p ¼ 1:6 � 10 4 Pa, which is representative for the brain tissue, see, e.g. Rashid et al. (2012), Basser (1992), Weaver et al. (2012), Perriñez et al. (2010) and Soza et al. (2005). We solve the anisotropic, quasi-static Biot's equations on a sphere of radius R ¼ 2 cm (resembling a tumor size), from which a very small cavity of radius a ¼ 0:03 cm is cored. The latter mimics the tip of a needle related to an injection at a constant pressure, as described in Basser (1992), where these two values for the radii R and a are used.
We prescribe the pressure p ð0Þ to be equal to 666:4 Pa at t � 0on the internal boundary, while free drainage is applied on the external boundary. Free stress conditions are considered on the internal boundary of the sphere, while we consider free radial stress and zero polar and azimuthal displacements on the external boundary of the sphere.
Furthermore, it is worth noting that the model is anisotropic due to Table 4 The effective poroelastic properties derived from the underlying microstructural properties. the cubic symmetry of the effective elasticity tensor, which is in turn induced by a microstructure organised in cubic cells, Penta and Gerisch (2015); Dehghani et al. (2018). Therefore, the solution is genuinely three-dimensional, i.e., even though the results are presented in terms of the (non-dimensionalised) radius, the solution does, in principle, depend on the azimuthal and polar angles. The latter are fixed to zero in this section, while the role of anisotropy is highlighted in the following section. This problem is solved in Abaqus (using 4128 C3D20T elements, i.e. three dimensional Serendipity elements with quadratic approximation for the displacement field and linear approximation for the temperature, and, accordingly, pressure field) via transient analyses in which we choose Δt ¼ 5e À 02 s as a sufficiently accurate time step. The pressure in the core is applied at the beginning of the analysis and we present the results for the cases reported in Table 4 at various time points. Fig. 3a to d shows the displacement of case 1 to case 4 at different times, respectively. According to our results, time plays a crucial role in determining displacement maps. The latter are normally measured and interpreted by means of poroelastography (Perriñez et al. (2010); Weaver et al. (2012)).
The profiles of the displacements of all the cases together are shown in Fig. 4a to d at time points t ¼ 50; 100; 500; 1000 s. The effect of solid matrix compressibility becomes more evident at longer times. We can therefore observe that a sufficiently long time should be considered to obtain an accurate estimation of poroelastic properties via poroelastography techniques. In fact, modelling and simulation strategies such as the ones presented here can be used together with poroelastography to help quantify the minimum time point after which the results become less and less affected by the passage of time. Due to the fact that reaching the steady-state sometimes takes a long time, we set an indicator for the minimum time at which the effect of solid matrix Poisson's ratio becomes distinguishable. In general, as one important indicator of the evolution of the problem (fluid injection) is the distribution of volumetric strain, an appropriate criterion for identifying a time which is sufficiently long for reliable experimental measurements can be identified from its profile. The volumetric strain ratio in Fig. 5d is the ratio of volumetric strain at the dimensionless radius r=a ¼ 5 divided by its maximum value located at r=a ¼ 1. In the context of the present analysis, we propose that the time should be considered sufficiently long when the volumetric strain ratio at r=a ¼ 5 reaches 0.1 which happens at around t ¼ 500 s. Volumetric strain ε v ¼ trðεÞ ¼ ε 11 þ ε 22 þ ε 33 vs. the dimensionless radius at φ ¼ 0 is shown in Fig. 5a. The latter profile is, as expected, in a good agreement with the pore pressure distribution in radial direction demonstrated in Fig. 5b. The mentioned profiles are important for the determination of the effects of the macroscale poroelastic response on the microscale properties which will be explained afterwards. Fig. 5c exemplifies the pore pressure gradient which imposes pore fluid velocity in this framework. The times t ¼ 2000 s and t ¼ 3000 s are included in Fig. 5a, b, and 5c in order to show that the model response does not change considerably after t ¼ 1000 s.
Fluid injection, percolation, and drainage are very important factors  Table 4.
in the response of the model as well as the determination of the poroelastic material properties. Immediately after application of the cavity pressure, the tissue behaves as almost incompressible (which is also reported in other studies such as Islam et al. (2018)), which is shown in Fig. 5a (see e.g. graph for t ¼ 5 s). By comparing Fig. 5a and b, we can conclude that the volumetric strain is mostly driven by the hydrostatic pore pressure, which is negligible when evaluated immediately after imposing the cavity pressure. However, as time passes, the pressure which is continuously applied in the cavity drives the corresponding changes in the effective solid and fluid volume portions from which the compressible behaviour of the poroelastic material can be deduced. Finally, in order to further elucidate this behaviour and to demonstrate the qualitative reliability of our framework, we focus on Equation (3). The latter can be either derived from the microstructure (Burridge and Keller (1981)) or from standard mass-balance considerations and relates pressure variations to solid and fluid volume changes. We can rewrite (3) as The integration in time of Equation (28) can be approximated as where 〈u〉 rf is the fluid displacement and where Δr x ⋅〈u〉 rf ¼ Δε vf is related to the change in the volume of the fluid in one time increment.
We can also write r x ⋅u ð0Þ ¼ ε v which is the volume change of the solid phase. In the case of α being proportional to the identity, i.e. α ¼α I, we can write α : Δr x u ð0Þ ¼α Δr x ⋅u ð0Þ ¼α Δε v then, We can also write where m is the number of the current increment, ε m vf and ε m v are the total volume change of fluid and solid at that increment, and Δε n vf and Δε n v are the fluid and solid volume change of the increment n, respectively.
From Fig. 6 one can see that, as expected, the mean pore pressure between the mentioned cases in Table 4 is approximately constant, i.e. it is not affected by the MSMC. Consequently, only the macroscale coefficients M and α as well as ε v cause the different volumes of the injected fluid in different cases, which is shown together with the solid volume change in Fig. 7. The mean values are determined by taking the weighted average value (with respect to the element dimensions) of the variables of the whole discretised domain.
As predicted by Biot's theory of poroelasticity, Fig. 7 reveals that the higher the compressibility, the wider the difference between the change in volume of the solid and the fluid, which is in turn directly related to changes in pressure. This type of analysis can also be used a posteriori to evaluate the change in fluid volume fraction (i.e. porosity).

Anisotropy and steady-state analysis
In order to observe the effect of microscale solid matrix compressibility on the macroscale anisotropy, we consider the mechanical properties of the cases in Table 4. We calculate the shear modulus for an isotropic material with the same Young's modulus and Poisson's ratio of the mentioned cases as well as the deviation of each effective shear modulus from the isotropic one for each case via  where i indicates the case number, μ i iso is the shear modulus for an isotropic material with the same Young's modulus and Poisson's ratio of the case i, and where S i is the deviation coefficient of the effective shear modulus from μ i iso . The values of the mentioned parameters are provided in Table 5.
The parameter "Deviation from isotropy" in Table 5 highlights the fact that the less the value of the solid matrix Poisson's ratio ν, the more the deviation from the isotropy S i . The effect of the cubic anisotropy which arises as a direct consequence of a cubic cell microstructural arrangement (see also Penta and Gerisch (2015)) is a heterogeneous displacement distribution in all the directions in spherical coordinates (which resembles a cubic shape) and is shown in Fig. 8a and b. Moreove, from the latter figures, we observe that the displacements at different azimuthal locations start to deviate at a specific radius at which the volumetric strain is very small (see Fig. 5a). Poroelastography measurements taken at different angles could be compared to estimate the mechanical properties, in particular the degree of anisotropy, of poroelastic tissues. A comparison between all the cases at hand is made at r=a ¼ 20 and t ¼ 1000 s in Fig. 8c along the azimuthal direction (which has equivalent profile along the polar one) demonstrating that the more the matrix is compressible, the more the radial displacement at φ ¼ 45 deviates from the displacements at zero azimuthal angle, and that the impact of the cubic symmetry decreases with time. Fig. 8d highlights that the heterogeneous displacement profile is more pronounced at short times, such as t ¼ 50 s. However, the absolute values of the displacements are much smaller for longer times, such as t ¼ 1000 s (for example, u r ðφ ¼ 0Þj t¼50 s =u r ðφ ¼ 0Þj t¼1000 s ¼ 0:086). It is noteworthy that, in the previous section, all the values against the dimensionless radius have been extracted at φ ¼ ϕ ¼ 0 which are the azimuthal and polar directions in spherical coordinates. Fig. 8 indicates that although the load and geometry are totally symmetric and independent of azimuthal and polar directions, the anisotropic effective elasticity tensor (with cubic symmetry) imposes the dependency of the response on the mentioned spatial directions.

Suddenly vs. gradually applied loading
In this section, we clarify the effects of sudden versus gradual loading conditions. In the former, we apply the entire pressure instantaneously at the beginning of the first time increment (i.e. transient analysis in ABAQUS), whereas in the latter, in each time increment, we apply a portion of the total cavity fluid pressure which is equal to the portion of the total time. For gradual loading, an incremental loading of 0:001 p 0 is applied every second, for the total time considered being t ¼ 1000 s (steady-state analysis in ABAQUS). Fig. 9 shows the radial displacement and volumetric strain versus dimensionless radius of case 4 at φ ¼ 0 for the gradual loading considered.
The difference between Figs. 9a and 4d shows the fact that the loading condition also plays an important role. It can be seen from Fig. 9b that the difference between the volumetric strain for the different cases are comparatively small under gradual loading.
Comparing the results from gradual and instant cavity pressure application, the effects of solid matrix compressibility are more pronounced for sudden loading.

Summary and future works
In this work we have presented the homogenized solution of the three-dimensional Biot's system of equations as derived in Burridge and Keller (1981). We have obtained the results via finite element numerical simulations, as illustrated in Section 3. In Dehghani et al. (2018), the authors have presented a systematic computational approach which can be used to evaluate all the relevant macroscale poroelastic coefficients as prescribed by the asymptotic homogenisation technique. Here we have exploited the findings reported in Dehghani et al. (2018) to compute the corresponding macroscale solution, thus clearly highlighting the link between relevant poroelastic maps, such as displacements and pressures, and the pore-scale properties, in particular focusing on the microscale solid matrix compressibility.
We have performed a parametric analysis by varying the Poisson ratio of the elastic matrix (which is considered isotropic consistent with Dehghani et al. (2018)), and have obtained displacements and pressures profiles at various time points. The analysis shows that the role of the solid matrix compressibility in determining displacements and volumetric strains becomes more and more pronounced as time passes (see Fig. 4), thus suggesting that current poroelastographic measurements should account for a sufficiently long time interval to provide an accurate estimate of the poroelastic parameters.
Furthermore, the present results also serve as a basis to evaluate the mechanical anisotropy which arises from the microstructure. In particular, this is reflected in the macroscale displacements and pressure maps ε m v f dv, respectively) over time for four cases. C1, C2, C3, and C4 represent case 1, case 2, case 3, and case 4, respectively.

Table 5
The effective mechanical properties derived from the underlying microstructural properties compared with the shear moduli calculated for isotropic cases with the same effective Young's modulus and Poisson's ratio.  Fig. 8), and is supposed to be present even when considering a geometry which is invariant under permutation of the three orthogonal axes, see Dehghani et al. (2018). Measurements taken at different angles may provide possibilities towards a better understanding of the anisotropic response of poroelastic materials, as well as suggestions on the pore-scale structural arrangement when the latter is not accessible via proper imaging techniques. Finally, results in terms of mean solid and volumetric strains (see Fig. 7) provide a successful qualitative benchmark of the implementation of Biot's theory of poroelasticity and could quantitatively be used to estimate the variation in porosity that takes place depending on the loading condition.
We have presented our results in the context of biological tissues, which are often characterised by a poroelastic mechanical response. The next natural step is the application of the presented framework to support specific experimental measurements, for example related to cell aggregates (as highlighted in Dehghani et al. (2018)) and organs, such as the brain Mousavi et al. (2014). Furthermore, the present analysis could be generalized to incorporate inertial effects (see also the analysis carried out in Gao et al. (2015) in the context of representative element volumes), so that this framework could be applied also in scenarios where ultrasound techniques are involved, such as the bone structure discussed in Fellah et al. (2013).
These results are to be considered as a first step towards a better understanding of the mechanics of poroelastic tissues. The framework presented here could also be used as a basis for a suitable implementation of the recent advances in the literature concerning multiscale poroelastic tissues, see, e.g. Penta et al. (2014); Penta and Merodio (2017). Predictions for the latter models would account for growth and vascularisation and therefore could be applicable to more realistic biological scenarios of interest.

Declaration of competing interest
The authors declare that they have no known competing financial  Table 5.