Mesoscale Poroelasticity of Heterogeneous Media

Poromechanics is dedicated to the modeling and prediction of how porous materials deform in response to various external loadings. These loadings range from fluid–solid interactions by a variety of pressures at the liquid–solid interface to complex physical chemistry phenomena at the pore scale that produce a mechanical deformation (including fracture) of the solid. The classical backbone of poromechanics is based on continuum theories, ever since Maurice A. Biot defined the kinematics of deformation of the skeleton within the classical continuum mechanics framework as the reference for the description of the flow of the liquid phase through the pore space (Biot 1941), with the state equations for stress, $\mathbf{\Sigma }$, and porosity change, $\varphi -{\varphi }_0$, given in the linear poroelastic case bywhere $E_{\mathrm{pot}}$ = potential energy of the solid phase of the solid-pore composite of volume $V$ subjected to an average strain $\mathbf{E}={⟨\mathit{ϵ}⟩}_V$ at the boundary $\partial V$ and a pressure $p$ at the solid–pore interface; $\mathbb{C}$ = fourth-order elastic stiffness tensor; $\mathbf{b}$ = second-order tensor of Biot pore pressure coefficients; and $N$ = solid Biot modulus. This continuum framework also provided the backbone for the development of the close-to-equilibrium thermodynamics framework of irreversible deformation of porous media pioneered by Coussy (1995) and its extension to a large range of phase change and adsorption phenomena (Coussy 2010). In the same vein, microporomechanics theories can be viewed as refined extensions of the continuum framework to the micro scale, in that they adapt continuum micromechanics theory (Suquet 1987; Zaoui 2002) to the specific nature of porous materials viewed as solid–pore composite materials (Dormieux et al. 2002, 2006). Although continuum poromechanics theory has entered and transformed many engineering fields ranging from civil and environmental engineering and geophysics applications to biomechanics and the food industry (e.g., Hellmich et al. 2013), the intrinsic limitations of the theory relate to the very foundations of the continuum model, including scale separability and its impact on the relevance of the differential operators defining the momentum balance and displacement–strain operators. This is a serious limitation of the theory in its applicability to highly heterogeneous materials. For instance, such a continuum theory will fail for microstructure resolutions achieved by microcomputed and nanocomputed tomography (CT) imaging techniques of highly heterogeneous materials, in which the characteristic length scale of the heterogeneity is of a similar scale as the sample size, or for multiscale heterogeneous materials for which a single representative elementary volume (rev) cannot be defined. It is for such systems that a discrete form of poromechanics theory is proposed, in which physical interactions replace volume descriptors. This approach is much akin to molecular representations of material systems with interaction forces between mass points derived from potentials that define the out-of-equilibrium state of the system with regard to a relaxed equilibrium configuration.

Herein, the elements of such a discrete poromechanics approach are developed using statistical mechanics ensemble definitions within the context of the lattice element method (LEM) (Topin et al. 2007; Affes et al. 2012) using the framework of effective potentials (Laubie et al. 2017b). By way of validation, some pore–solid morphologies are revisited to determine poroelastic constants within and beyond the classical continuum limits of scale separability.

Consider a porous material composed of a solid (volume $V_s$) and pore space (volume $V_p$). Following the lattice element method (Topin et al. 2007; Affes et al. 2012; Laubie et al. 2017b), the two domains are discretized into a number of unit cells (or voxels), the center of which defines a mass point that interacts with a fixed number of neighboring mass points forming a regular or irregular lattice structure. The interaction forces and moments between two mass points $i$ and $j$ derive from an effective potential $U_{ij}$ as a function of the translational, ${\overrightarrow{\delta }}_i={\overrightarrow{x}}_i-{\overrightarrow{X}}_i$, and rotational, ${\overrightarrow{\vartheta }}_i$, degrees of freedom, where ${\overrightarrow{X}}_i$ and ${\overrightarrow{x}}_i$ denote the position vectors of mass point $i$ in the reference and the deformed configurations, respectively [see Laubie et al. (2017b) for a detailed derivation]where ${\overrightarrow{r}}_{ij}=l_{ij}^0{\overrightarrow{e}}_n^{ij}$ = vector connecting node $i$ to node $j$ of rest length $l_{ij}^0$ and oriented by the unit vector ${\overrightarrow{e}}_n^{ij}$ in a local orthonormal basis $({\overrightarrow{e}}_n,{\overrightarrow{e}}_b,{\overrightarrow{e}}_t)$. For such a discrete system, the stresses are modeled using the virial expression (Christoffersen et al. 1981) $\mathit{\sigma }={\rho }_c⟨\overrightarrow{r}\otimes \overrightarrow{F}⟩$, where ${\rho }_c$ represents the number of interaction bonds per unit volume, $⟨·⟩$ denotes the first moment of $\overrightarrow{r}\otimes \overrightarrow{F}$ distribution over interaction bonds, while neglecting the momentum term. In the LEM for mass point $i$, this virial expression can be written aswith $V_i$ denoting the volume of the unit cell and $N_i^b$ representing the number of node $i$’s neighboring mass points. The Virial expression provides a truly discrete description of the system as opposed to the continuum-based stress definition used in classical finite-element–based approaches. The stress in volume $V$ composed of a total of $N_t$ unit cells is simply the volume average of the local stresses; that is

What thus differs between different material domains is the interaction potential from which forces and moments are derived.

The poroelastic properties of materials form much of the backbone of application of poromechanics theory. This includes the elasticity tensor $\mathbb{C}$, the tensor of Biot coefficients $\mathbf{b}$, and the solid Biot modulus $N$. From the composite structure of porous materials, it is readily understood that these macroscopic properties call for averages. Such averages are best defined, in statistical mechanics, within the context of specific statistical ensembles, which—at least theoretically—include every possible microscopic state of the system. The advantage of using statistical ensembles for the determination of the poroelastic properties is that each ensemble is associated with a characteristic state function or thermodynamic potential that uniquely define—upon minimization—the equilibrium state of the system in function of a few observable parameters; much akin to the classical minima theorems of elasticity used in continuum mechanics, e.g., for the derivation of the state equations of poroelasticity Eqs. (1) and (2) (Dormieux et al. 2002, 2006). It is thus shown that making the link between statistical ensembles and such boundary conditions is quite helpful for the determination of the poroelastic constants from discrete simulations.

By way of application, the proposed discrete model of poroelasticity is implemented for simple cubic (SC) lattice systems. The LEM approach here used follows the approach developed by Laubie et al. (2017b) for solids, where specific details about calibration and numerical implementation of the method can be found. In short, for all the cases considered herein, a cubic simulation box of side length $L=a_0(n-1)$ composed of ${(n-1)}^3$ unit cells of size $a_0$ is used, where $n$ stands for the number of mass points in any given direction. The mass points form a regular lattice with their interactions encapsulated by a network of links that connects a mass point to its 26 neighboring mass points. Thus, mechanical information is propagated through this lattice network in 13 directions. This forms the so-called D3Q26 lattice structure consisting of 6 box-links of rest-length $l^0=a_0$, 8 cross-diagonal links of length $\sqrt{3}{a}_0$, and 12 in-plane-diagonal links of length $\sqrt{2}{a}_0$ (Fig. 1).

The idealized structures considered in this study represent a microtexture best captured by the Mori-Tanaka homogenization scheme. The Mori-Tanaka scheme is often associated with a matrix-inclusion microtexture where the matrix phase overshadows the mechanical response of the inclusion phase(s) while considering interactions between inclusions (in contrast to the dilute scheme; see, for instance, Dormieux et al. 2006). Furthermore, for a two-phase composite with spherical inclusions, the Mori-Tanaka scheme corresponds to the Hashin-Shtrikman bounds (Weng 1984) and, specifically for spherical voids, the upper Hashin-Shtrikman bound. However, the presented methodology to estimate poroelastic properties of heterogeneous media is independent of the microtextures being considered.

Although the discrete simulation results compare well against their continuum poroelastic counterparts for both the isotropic and transversely isotropic cases, a deviation is observed at higher porosity values that merits further discussion. Specifically, for small porosities, ${\varphi }_0<5\times {10}^{-3}$ (or $R/L\le 0.1$), the two approaches provide similar results. This is not surprising because—within this limit—scale separability, delineating the domain of application of the continuum models (here the Mori-Tanaka model), strictly applies. Beyond that limit, however, the results obtained from the discrete and continuum approaches begin to differ. One possible reason for the observed deviations is related to finite size effects associated with the finite size of the simulation box, noting that the elementary voxel size ($a_0$) remains much smaller than the size of the elementary heterogeneity at high porosities. To explore this further, two quantities, ${\delta }_{\mathrm{iso.}}$ and ${\delta }_{ti.}$, are defined to capture any deviations from the imposed elastic solid symmetry for the isotropic and transversely isotropic cases; that is, for the isotropic caseand for the transversely isotropic case

Using the elasticity constants $C_{ij}$ obtained from the simulations, Figs. 2(e and h) plot ${\delta }_{\mathrm{iso.}}$ and ${\delta }_{ti.}$ versus $R/L$, showing that for $R/L>0.1$ the effective (i.e., composite) elasticity content captured by the simulations departs from the material symmetries of the solid phase. Within the range of considered values, ${\delta }_{\mathrm{iso.}}\le 4$ in the isotropic case and ${\delta }_{ti.}\le 8$ in the transversely isotropic case. On the other hand, the simulation results deviate from the continuum solution for the poroelastic constants [Figs. 2(a and b)], for which the continuum solutions [i.e., Eqs. (41) and (42) in the isotropic case and Eqs. (43)–(45) for the transversely isotropic case] hold irrespective of elastic homogenization scheme. Thus, the observed deviation between discrete simulations and continuum calculations in the high porosity limit for elastic and poroelastic properties seems to be rooted in the finite size of the system as it challenges both the application of Eshelby’s solution for an ellipsoidal inclusion in an infinite medium (Eshelby 1957) and the Mori-Tanaka homogenization scheme’s subjection of inclusions to the first moment (mean) of matrix stresses (Mori and Tanaka 1973; Beneviste 1987). The same deviation is observed for highly disordered systems (Laubie et al. 2017a) but attributed to the high stress concentrations between pore walls. In this vein, the probability density function (PDF) of normalized solid stresses of the considered idealized pore-matrix structures in the $\mu VT$ ensemble is plotted in Fig. 5 for three different $R/L$ ratios. In violation of scale separability, for $R/L=0.157$ and $R/L=0.229$, normalized stresses follow Gaussian distributions. However, the long tails for $R/L=0.443$ indicate areas of high stress concentration, a feature not captured by mean-field–based theories of micromechanics. This is intimately related to the requirement of scale separability in homogenization theory. A key property of scale separability exploited in the theory of homogenization is that the local problem cannot see the boundaries (Pavliotis and Stuart 2007), which is clearly violated in cases of high $R/L$ ratios studied here.

Surface energy effects are incorporated in poromechanics by making a distinction between the free energy stored elastically into the solid matrix ${\psi }_s$ and energy $u$ stored at the solid–fluid interface, such that ${\mathrm{\Psi }}_s={\psi }_s(\mathbf{E},p)+u(\mathbf{E},p)$, with the energy balance for the interface at equilibrium expressed as (Vandamme et al. 2010; Brochard et al. 2012)where ${\tilde{\sigma }}^s$ denotes surface stress and $s$ represents the actual area of the pore walls per unit volume of porous material in its reference configuration. Furthermore, for example, in the case of adsorption in a linearly elastic isotropic porous material, one can obtain material parameters ${\alpha }_{\epsilon }$ and ${\alpha }_{\phi }$ to quantify strain and porosity changes due to surface stresses, respectively (Vandamme et al. 2010). In this vein, the proposed method can be extended to capture adsorption-induced structural phase transitions in a porous material using an osmotic ensemble (Snurr et al. 1993; Mehta and Kofke 1994; Coudert et al. 2011)where $T$, $P$, ${\mathrm{\Psi }}_s$, $V$, $N_{ads.}$, and $V_m$ = temperature, pressure, the free energy of the solid in the absence of adsorbed molecules, the volume of the porous host, the number of adsorbed molecules inside the host, and the molar volume of the adsorbing species in its bulk state, respectively. Then, one seeks the structure that minimizes ${\mathrm{\Omega }}_{os.}$. Once this structure is obtained, $N_{ads.}(T,P)$ can be predicted with standard grand canonical Monte Carlo (GCMC) simulations. Classically, the main challenge of using Eq. (52) is access to ${\mathrm{\Psi }}_s$, which would be readily available via the LEM.

As the resolution of microtexture and heterogeneity of porous materials is progressing rapidly thanks to advancements in, e.g., CT imaging techniques (Hubler et al. 2017), there is a need to adapt the tools of poromechanics to model and predict the deformation of porous materials in response to various external loadings. The discrete poromechanics approach proposed and implemented in the lattice element method aims at contributing to this effort well beyond the classical mean-field–based theories of continuum microporomechanics, which do not capture microtextural information beyond one-point correlation functions and are confined in their application by the scale separability condition. Specifically, the discrete nature of the approach provides access to local stresses and displacements as well as force flow in a heterogeneous system, which can illuminate the path for understanding stress and strain localization in a multiphase porous composite and form a basis for subsequent refinements to include irreversible deformation (including fracture), deformation during flow, and so on. The following points of observation deserve attention:

This work was funded by X-Shale Hub: the Science and Engineering of Gas Shale, a collaboration between Shell, Schlumberger, and the Massachusetts Institute of Technology, enabled through MIT’s Energy Initiative. The authors would like to acknowledge Vincent Richefeu (Universite Joseph Fourier), Jean-Yves Delenne (Montpellier SupAgro), and Saeid Nezamabadi (Universite de Montpellier), who provided the backbone of the LEM code used for simulations. FR would like to acknowledge the support of the ICoME2 Labex (ANR-11-LABX-0053) and the A*MIDEX projects (ANR-11-IDEX-0001-02) cofunded by the French program Investissements d’Avenir, managed by the French National Research Agency (ANR). The authors would also like to thank the reviewers for their insightful comments.