Variational estimates of the poroelastic coefficients

Saturated, isotropic, poroelastic materials are classically described by their elastic stiffness, one Biot coefficient and one Biot modulus (Coussy 2010; Dormieux, Molinari, and Kondo 2002). The situation becomes more complex for unsaturated, isotropic, poroelastic materials that require multiple Biot coefficients and moduli (Coussy and Brisard 2009). These poroelastic coefficients are dependent and linked by several linear relationships. Micromechanical estimates of these coefficients have been proposed by several authors (Ulm, Constantinides, and Heukamp 2004; Pichler and Hellmich 2010). However, these estimates may fail to fulfill the linear relationships that relate the exact poroelastic coefficients. This might be regarded as an undesirable inconsistency of the model. In this work, we propose new, consistent (in the sense that the above mentioned linear relationships are preserved) estimates of the poroelastic coefficients. Our point of departure is the principle of Hashin and Shtrikman, suitably extended to eigenstressed materials (Bornert et al. 2001). Adopting stress-polarization fields that are similar to the eigenstress-free case (Willis 1977; Ponte Castañeda and Willis 1995) allows us to derive variational estimates of the poroelastic coefficients, which can be shown to fulfill all known linear relationships required from the exact values. After outlining the derivation within the general framework of eigenstressed, heterogeneous materials, the results will be specialized to poroelasticity. This will lead to a variational justification of the ad-hoc pore isodeformation assumption (Coussy and Brisard 2009).

This is the post-print (ie final draft after peer-review) version of the conference paper Variational estimates of the poroelastic coefficients by Sébastien Brisard and Siavash Ghabezloo. The published version is available at http://ascelibrary.org/ (https://doi.org/10.1061/ 9780784480779.157. It is recalled that the ASCE Authorship Originality and Copyright Transfer Agreement (last retrieved 2017-07-13) states Authors may post the final draft of their work on open, unrestricted Internet sites or deposit it in an institutional repository when the draft contains a link to the published version at www.ascelibrary.org. "Final draft" means the version submitted to ASCE after peer review and prior to copyediting or other ASCE production activities; it does not include the copyedited version, the page proof, a PDF, or full-text HTML of the published version.

INTRODUCTION
The present paper is dedicated to the homogenization of eigenstressed, N-phase materials. We consider a representative volume element (RVE) Ω of such a material; phase α = 1, . . . , N occupies the domain Ω α ⊂ Ω. The eigenstressed, linear elastic constitutive law of phase α reads where C α denotes the elastic stiffness and α the eigenstress (σ: stress; ε: strain). 2 The effective behavior of such materials is governed by the following macroscopic equations (Dvorak and Benveniste 1992) where Σ = σ and E = ε denote the macroscopic stress and strain (volume averages over the whole RVE Ω); ε α denotes the average strain over phase α. Unless otherwise stated, all sums in this paper span all N phases. A 1 , . . . A N are the well-known strain-localization tensors, and are the effective stiffness and eigenstress, respectively. The fourth-rank tensors D αβ are the eigenstrain influence tensors introduced by Dvorak and Benveniste (1992). They enjoy the following properties where f α denotes the volume fraction of phase α = 1, . . . , N.
The macroscopic properties of the eigenstressed, heterogeneous material are classically obtained from the application of equations (2a), (3a) and (3b) to the solution of the following problem on the RVE Ω div σ( However, solving the above local problem is often quite difficult, and estimates of the A α and D αβ are of great practical value. Revisiting the work of Mori and Tanaka (1973), Benveniste (1987) proposed estimates of the strain localization tensors. His approach was then extended by Dvorak and Benveniste (1992) to the eigenstrain influence tensors. The estimates of Dvorak and Benveniste were however restricted to materials with aligned inclusions of identical (ellipsoidal) shape. This limitation was later overcome by Pichler and Hellmich (2010), who derived closed-form expressions of the estimates of the eigenstrain influence tensors for any distribution of ellipsoidal inclusions.
It is well-known that the Mori-Tanaka approach might deliver estimates of the effective stiffness that are not symmetric. Such unphysical estimates might occur in the case of non-aligned inclusions and/or inclusions of different shapes (Benveniste, Dvorak, and Chen 1991;Ferrari 1991;Schjødt-Thomsen and Pyrz 2001). Pichler and Hellmich (2010) have shown that in those cases where the classical Mori-Tanaka estimate of the effective stiffness is symmetric, their estimates of the eigenstrain influence tensors are consistent with properties (4a), (4b), (4c) and (4d). Conversely, if the Mori-Tanaka estimate of the effective stiffness is not symmetric, the Pichler-Hellmich estimates of the eigenstrain influence tensors are not physically acceptable.
Our aim is to propose alternative estimates of these tensors, which are always acceptable. Our estimates are derived from the variational framework introduced by Hashin and Shtrikman (1962b) and later extended to eigenstressed materials by Bornert et al. (2001). Our main contribution is to show that these estimates always verify identities (4a), (4b), (4c) and (4d)]. Such robustness comes with a price, though. Indeed, it requires additional information (namely, the so-called P-tensors of Ponte Castañeda and Willis 1995) on how the phases are distributed.

THE PRINCIPLE OF HASHIN AND SHTRIKMAN WITH EIGENSTRESSES
The variational principle of Hashin and Shtrikman was introduced in Hashin and Shtrikman (1962b) (see also Willis 1977) and later extended to eigenstressed materials by Bornert et al. (2001). It is stated below in a form amenable to the derivations proposed in the next section.
Introducing a reference (linear, elastic) material of stiffness C 0 , the functional HS of Hashin and Shtrikman is defined as follows where C(x) [resp. (x)] is the phase-wise constant local stiffness (resp. eigenstress): C(x) = C α [resp. (x) = α ] for x ∈ Ω α . Γ 0 denotes the Green operator for strains of the large (but finite) domain Ω (Korringa 1973;Zeller and Dederichs 1973;Kröner 1974;Willis 1977). The secondrank, symmetric tensor field τ that appears in the definition (6) of the functional HS of Hashin and Shtrikman is called stress-polarization. For arbitrary reference materials, the principle of Hashin and Shtrikman states that i. the functional HS is stationary at (C − C 0 ) : ε, where ε is the strain that solves problem (5) and ii. the value of HS at this critical point is the macroscopic strain energy Furthermore, if the reference material C 0 is stiffer (resp. softer) than all phases, meaning that C 0 ≥ C α (resp. C 0 ≤ C α ) in the sense of quadratic forms for all α = 1, . . . , N, then this critical 4 point is in fact the unique minimum (resp. maximum) of HS where the second inequality holds for any choice of the trial stress-polarization τ.

DERIVATION OF THE VARIATIONAL ESTIMATES
From the principle of Hashin and Shtrikman, finding the critical point of HS over the whole space of stress-polarizations delivers the exact values of the effective stiffness C eff , strain localization tensors A α and eigenstrain influence tensors D αβ [see equations (7a) and (7b)]. Likewise, finding the critical point of HS over a (usually, finite dimension) subspace of stress-polarizations delivers estimates.
The obvious choice for this subspace is the space of phasewise constant stress polarizations: τ(x) = τ α for all x ∈ Ω α (α = 1, . . . , N), where τ 1 , . . . , τ N are constant, second-rank, symmetric tensors. In the eigenstress-free case, this results in the classical estimates of Hashin and Shtrikman (1962a) of the effective stiffness; these estimates turn out to be bounds provided that either condition in equation (8) is satisfied. The same approach is extended here to eigenstressed materials.
It can be shown that for large RVEs Ω, the functional of Hashin and Shtrikman HS defined in equation (6) evaluates to HS(τ; E, 1 , . . . , N ) = 1 2 E : In the above equation, the so-called P-tensors P αβ are weighted averages of the Green operator for strains Γ ∞ 0 of the whole space where S αβ denotes the two-point probability function of phases α and β (see e.g. Torquato 2002).
It should be noted that in Ponte Castañeda and Willis (1995), the P-tensors were denoted A (rs) (where r and s are phase indices in this publication). For weakly isotropic materials, S αβ (r) depends on the norm r of r only, and (Willis 1977) Stationarity of HS with respect to the τ α leads to the following linear system of equations (indexed by α = 1, . . . , N) In the above system, the stiffness C 0 of the reference material, the stiffness C α of the constituants (α = 1, . . . , N) and the microstructural parameters P αβ [defined by equation (10), α, β = 1, . . . , N] are given. The macroscopic strain E and phase eigenstresses 1 , . . . , N are loading parameters. Finally, the stress-polarization τ 1 , . . . , τ N are unknown.
Equation (12) can in general not be solved analytically. However, owing to linearity, it is formally possible to express the solution as follows (Brisard and Ghabezloo 2017) where the tensors A var α and D var αβ are found (numerically) from the inversion of equation (12). Plugging equation (13) into equation (9) and introducing it is found that Comparing equations (7b) and (15) finally shows that C var , A var α and D var αβ ought to be considered as variational estimates of C eff , A α and D αβ , respectively. For isotropic microstructures, these estimates coincide with that of Dvorak and Benveniste (1992) (which are identical to the estimates of Pichler and Hellmich 2010, in this case). Outside this particular case, they generally differ from the estimates of Pichler and Hellmich (2010), as will be discussed in the next section.

PROPERTIES OF THE VARIATIONAL ESTIMATES
Careful analysis of the linear system (13) shows that the variational estimates C var , A var α and D var αβ enjoy the same properties (4a)-(4d) as their exact counterparts (see Brisard and Ghabezloo 2017).
In particular, the variational estimate C var of the effective stiffness is always symmetric, unlike the Mori-Tanaka estimate. This suffices to prove that the proposed variational estimates of the eigenstrain influence tensors differ in general from their Pichler-Hellmich estimates. While the latter might take unphysical values (when the Mori-Tanaka estimate of the effective stiffness is not symmetric), the former are always consistent with identities (4a)-(4d).
The proposed variational estimates are therefore appealing alternatives to the Pichler-Hellmich estimates in cases where the latter do not satisfy identities (4a)-(4d).

APPLICATION TO POROELASTICITY
The above general framework can be specialized to poroelastic materials: the solid skeleton is made of N phases α = 1, . . . , N, while the porous network is decomposed in P families α = N + 1, . . . , N + P. Pores of the familiy α are filled with a fluid at pressure p. The effective behavior of the porous material is deduced from equations (2) and (3b) with α = 0 (α = 1, . . . , N), C α → 0 (α = N + 1, . . . , N + P) and α = −p α I (α = N + 1, . . . , N + P) where The second-rank tensor B α is the tensor of Biot coefficients; the scalars N αβ are the Biot moduli. For isotropic materials, the tensors of Biot coefficients degenerate to scalars (B α = b α I), and the classical equations of unsaturated poroelasticity (Coussy and Brisard 2009) are retrieved.
Variational estimates of the poroelastic coefficients can therefore be deduced, for a given microstructure, from the above variational estimates of the eigenstrain influence tensors. For isotropic microstructures, the variational estimates coincide with the estimates of Pichler and Hellmich (2010) when the solid matrix is selected as a reference material.
We further assume that the solid skeleton is homogeneous (κ S : bulk modulus; µ S : shear modulus) and that the pore space is filled with a liquid (S L : saturation; p L : pressure) and a gas (S G = 1 − S L : saturation; p G : pressure). Since the variational estimates coincide with the Pichler and Hellmich (2010) estimates in that case, the closed-form expressions proposed in Coussy and Brisard (2009) can be used. In particular and it is observed that the variational estimate of the Bishop coefficient b L /b is the saturation in liquid. This approximation is usually called pore isodeformation assumption (Coussy and Brisard 2009). It can in fact be shown that this assumption is exact for microstructures that realize the Hashin-Shtrikman upper-bound on the bulk modulus. Indeed, inequality (8) can be specialized for isotropic elasticity where E denotes the macroscopic volumetric strain (E = tr E). The above inequality must hold for all E, p L and p G . Therefore, the following matrix must be positive. In particular, κ eff ≤ κ var (upper-bound on the bulk modulus) and If the microstructure is such that κ eff = κ var , then we must have b L = b var L and the pore isodeformation assumption is rigorous.

CONCLUSION
Using the framework of Hashin and Shtrikman (1962b), we have presented a derivation of variational estimates of the eigenstrain influence tensors of a heterogeneous microstructure. It was shown that these estimates are always physically acceptable, contrary to Mori-Tanaka estimates. This robustness has a price: the P-tensors defined by equation (10) must be evaluated.
For porous materials, estimates of the poroelastic coefficients are readily deduced from the estimates of the eigenstrain influence tensors. Our analysis shows that for isotropic porous media with homogeneous solid skeleton, the pore isodeformation assumption is rigorous provided that the drained bulk modulus coincides with the Hashin-Shtrikman upper-bound.
It is known that for a wide range of microstructures, the drained bulk modulus is actually very close to its upper-bound. For such microstructures, our results indicate that the pore isodeformation assumption should be an excellent approximation.