Combining analysis of random elastic polycrystals with poroelasticity for granular composites having orthotropic porous grains and fluid-filled pores

International Journal of Engineering Science 72 (2013) 11–21

Contents lists available at SciVerse ScienceDirect

International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci

Combining analysis of random elastic polycrystals with poroelasticity for granular composites having orthotropic porous grains and fluid-filled pores James G. Berryman ⇑ University of California, Lawrence Berkeley National Laboratory, One Cyclotron Road MS74R316C, Berkeley, CA 94720, USA

a r t i c l e

i n f o

Article history: Received 30 April 2013 Received in revised form 21 June 2013 Accepted 21 June 2013 Available online 21 July 2013 Keywords: Poroelasticity Random polycrystals Orthotropic systems

a b s t r a c t Analysis of random polycrystals has typically been applied to solid grains of anisotropic elastic materials. Poroelastic analysis has similarly been applied to otherwise isotropic systems with pores having a variety of shapes and filled with fluids. Present effort is focused on combining these two types of geomechanical analyses by treating anisotropic (specifically orthotropic) poroelastic grains jumbled together to form an overall isotropic polycrystalline poroelastic material. The resulting problem is approximately twice as difficult to solve as the typical elastic polycrystal problem because the polycrystal analysis must be carried through twice: once for the drained (pore-fluid free to escape) poroelastic constants, and again for the undrained (pore-fluid trapped) poroelastic constants. As should be anticipated, poroelastic effects induced by trapped fluid are significantly stronger for the effective bulk modulus than for the effective shear modulus. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Poroelastic analysis (Berryman, 2010, 2011; Brown & Korringa, 1975; Biot & Willis, 1957; Gassmann, 1951; Skempton, 1954; Thigpen & Berryman, 1985) considers the range of responses possible in fluid-saturated porous earth materials when mechanical, seismic, acoustical, ultrasonic or other types of applied stresses are acting on these systems. The inherent symmetries of common poroelastic systems are often transversely isotropic if the entire sample or region to be studied is layered, or they might be more complicated due to various geological processes involving uplift and folding, in which case one of the most general cases often considered for analysis purposes is (among others) that of orthotropic symmetry (Berryman, 2011). For random polycrystals, the well-known Voigt (1928) and Reuss (1929) bounds, and the related Hill (1952) elastic constant estimators, are all fairly easily calculated once the elastic constants of individual anisotropic crystal grains are known. The more sophisticated Hashin–Shtrikman bounds (Hashin & Shtrikman, 1962) and similarly related self-consistent (Hill, 1965) estimators for these same constants are all comparatively more difficult to compute. Some recent work (Berryman, 2005, 2011) has shown how to simplify methods used to quantify estimators to some extent. The present effort first provides an overview to show how this extra work does produce added value, since – in particular – the well-known Voigt–Reuss–Hill estimators often do not fall within the Hashin–Shtrikman bounds, but the self-consistent estimators (for good reasons Willis, 1981) have always been found to do so. Various estimators termed ‘‘self-consistent’’ can take different forms, and not all of these forms are equivalent or equally valid. In particular, Hill’s version of self-consistency is not the same as the ones that are based on physical arguments and ⇑ Tel.: +1 510 486 5349; fax: +1 510 486 5686. E-mail address: [email protected] 0020-7225/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2013.06.007

12

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

scattering theory such as the ones that are sometimes called the ‘‘coherent potential approximation,’’ including Soven (1967), Taylor (1967), Gubernatis and Krumhansl (1975), Berryman (1980a), and Olson and Avellaneda (1992). This class of approximations is one considered specifically by Willis (1981), and is one that has been shown to give generally reliable results – by which statement we mean to imply (in part) that the results have also been found to be consistent with the rigorous bounds. A related class of approximations, also called ‘‘self-consistent,’’ is discussed by Kanaun and Levin (2008). Earlier work on Hashin–Shtrikman bounds or on self-consistent estimators has concentrated mostly on one or the other of these approaches, without making an effort to compare, contrast, and/or mutually validate them. Some recent work of the author (Berryman, 2005, 2011) has addressed some of these cross-validation issues, including the orthorhombic case that will be the main emphasis here. In outline, the paper first introduces poroelastic analysis and then the polycrystal problem. Section 2 treats the poroelastic analysis. Sections 3 and 4 emphasize the drained (fluid freely flowing) and undrained (trapped fluid) constants. Section 5 then introduces the orthotropic polycrystal problem within this same context. Section 6 presents the mathematical formulation needed to analyze the bounds and estimators for these polycrystal problems. Section 7 presents an overview of the polycrystal results including results for seven examples, and also some comparisons among both rigorous and approximate methods. Section 8 then makes use of the self-consistent polycrystal results by treating these bulk and shear modulus values as effective macroscopic grain constants. Using these values, it is then straightforward to estimate drained constants for specific materials when some assumption is made about shapes and density of cracks within these average grains. Then, it also becomes possible to estimate undrained constants. Results for seven examples are tabulated, and three of these are also plotted here. Section 9 provides an overview of the various results obtained, and finally summarizes our main conclusions. 2. Analysis of anisotropic poroelastic grains Consider that each poroelastic grain is anisotropic possibly due to some nonrandom alignment of the solid constituents and/or to the presence of oriented fractures or pores. We therefore consider the orthorhombic anisotropic version of the poroelastic equations:

0

e11

1

0

S11

S12

S13

Be C B S B 22 C B 12 B C¼B @ e33 A @ S13 f b1

S22

S23

S23 b2

S33 b3

r11 1 B C b2 C CB r22 C CB C: b3 A@ r33 A pf c b1

10

ð1Þ

The eii (no summation over repeated indices) are strains in the i = 1,2,3 directions. The rii are the corresponding stresses. The fluid pressure is pf. The increment of fluid content is f, which (like the strains) is dimensionless. There are two standard types of boundary conditions used in poroelastic studies: (a) drained versus undrained and (b) jacketed versus unjacketed. The drained versus undrained language arose in the soils context (Skempton, 1954), and implies the existence of natural barriers to flow when undrained, and freely flowing fluid in the drained scenario. The jacketed versus unjacketed language arose in discussions of laboratory data on porous samples (Biot & Willis, 1957), typically being immersed in a fluid bath. In order to arrive at boundary conditions similar to those undrained examples found in nature, it is often necessary to wrap the porous (but fluid-saturated) sample in some type of jacketing material (often rubber). Drained conditions are then achieved by removing the jacketing material and placing the bare porous sample into a fluid (often water) bath. Dry conditions are not necessarily desirable because cementing materials in dry samples may be harder than those in fluid-saturated systems, and therefore provide unrealistic results for comparisons to field data. The drained (fluid not trapped) compliances are Sdij  Sij . The drained Reuss average bulk modulus (Reuss, 1929) is defined via

1 K dR



X

Sdij

ð2Þ

ij¼1;2;3

a quantity which is what one commonly takes to be the definition of the bulk modulus of such a simple (non-heterogeneous) anisotropic system. The corresponding undrained compliances will be symbolized by Suij . For the Reuss average undrained bulk modulus K uR , we have drained compliances replaced by undrained compliances in a formula otherwise identical in form to (2). Off-diagonal coefficients bi ¼ Si1 þ Si2 þ Si3  1=3K gR , where K gR is again the Reuss average modulus of the grains – i.e., simply replace superscript d’s with g’s in (2) to determine K gR . The alternative Voigt (1928) average (also see Hill, 1952) of the P stiffnesses plays no direct role in the poroelastic discussion. Then, finally, coefficient c ¼ i¼13 bi =B in (1), where B is the second Skempton (1954) coefficient, which will be defined shortly. The shear terms due to twisting motions (i.e., strains e23, e31, e12 and stresses r23, r31, r12) are excluded from this part of the presentation because they typically do not couple to the modes of interest for anisotropic systems having orthotropic symmetry, or to more symmetric systems such as those having either transversely isotropic or isotropic symmetry. Summing the off-diagonal coefficients bi, we find

b1 þ b2 þ b3 ¼

1 K dR



1 aR ¼ ; K gR K dR

ð3Þ

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

13

where we have introduced (similar to the isotropic case) a Reuss effective stress coefficient:

aR  1  K dR =K gR :

ð4Þ

Furthermore, we have

c¼

  b1 þ b2 þ b3 aR 1 1 ¼ d þ/  g ; Kf KR B KR

ð5Þ

since a rigorous definition in this notation by Berryman (2011) for Skempton (1954) B coefficient is given by

B

1  K dR =K uR 1

K dR =K gR

¼

aR =K dR

d R =K R

a

 : þ / 1=K f  1=K gR

ð6Þ

Another related formula that will prove very useful to us later is:

K uR ¼

K dR 1  aR B

ð7Þ

for the pertinent undrained bulk modulus. More general versions (Brown & Korringa, 1975) of the B definition may include another bulk modulus for pore response that differs from the grain response if the medium consists of a more heterogeneous collection of grains and/or pores, not being included within our current scope. (But see Brown & Korringa, 1975 for a thorough discussion of this point.) We will also return briefly to this issue in the final technical section of the paper. With this one caveat, all the formulas presented are rigorous statements based on the anisotropic analysis. Appearances of Reuss average quantities K dR and aR are rigorous statements, not approximations (except for the already noted limitations). The specific choices of notation made will also help us to emphasize the similarity between rigorous anisotropic and isotropic formulas, such as those of Gassmann (1951). 3. Off-diagonal poroelastic coefficients bi Results follow next for the off-diagonal bi coefficients. Then, a general proof of their correctness is presented. The coefficients bi are determined by

bi ¼ Sdi1 þ Sdi2 þ Sdi3 

1 ; 3K gR

ð8Þ

where K gR is the Reuss average of the grain modulus. Eq. (8) holds for homogeneous isotropic grains, such that K gR ¼ K g . However, when the grains themselves are anisotropic, we also need to allow for this possibility by defining three directional grain bulk moduli determined using:

1  Sgi1 þ Sgi2 þ Sgi3 ¼ Sg1i þ Sg2i þ Sg3i 3K gi

ð9Þ

for i = 1,2,3. The second equality follows from symmetry of the compliance matrix. We call these quantities in (9) the partial grain-compliance sums, and the K gi are pertinent directional grain bulk moduli. Then, the formula in (8) may be appropriately replaced by

bi ¼ Sdi1 þ Sdi2 þ Sdi3 

1 3K gi

:

ð10Þ

The preceding results are for perfectly aligned grains. If the grains are instead perfectly randomly oriented, then it is clear that the formulas in (8) hold as before, but now the definition of the Reuss average grain bulk modulus K gR must be reformulated in analogy to (2). All of these statements about the bi’s are easily proven by considering a particular combination of the applied stresses, such that r11 = r22 = r33 = pc = pf. Then, from (1), we have:

    pf eii ¼ Sdi1 þ Sdi2 þ Sdi3 pc þ bi ðpf Þ ¼ Sgi1 þ Sgi2 þ Sgi3 pf  3K gi

ð11Þ

in the most general of the cases to be considered. This result holds true for each value of i = 1,2,3. Such statements about the strain eii (no sum over i) always hold in the situation discussed, as it must be the same if these anisotropic (or inhomogeneous) grains are immersed in fluid, while measurements are taken of strains observed in each of the three directions i = 1,2,3, during variations of this uniformly applied pressure pf. This argument is also entirely analogous to ones given by Biot and Willis (1957) for isotropic, homogeneous samples. The relationship of coefficient c to the other coefficients is easily established because the main issue involves determining the role of the various constants contained in the definition of Skempton’s coefficient B (Skempton, 1954). Again, from (1), we find that

14

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

f ¼ 0 ¼ ðb1 þ b2 þ b3 Þrc  cpf

ð12Þ

for undrained boundary conditions. Then we have

pf b þ b2 þ b3 B¼ 1 ; pc c

ð13Þ

where pc = rc is the external confining pressure. Thus, the scalar coefficient c is determined uniquely, and given by

c¼

  b1 þ b2 þ b3 aR =K dR 1 1 ¼ ¼ aR =K dR þ /  g : Kf KR B B

ð14Þ

This formula also provides an alternative (but nevertheless equivalent) definition of Skempton’s second coefficient:

B¼

aR ; cK dR

ð15Þ

although this special result is limited to systems having homogeneous (and isotropic) grains. 4. Undrained compliance matrix S uij The undrained compliance matrix Suij can now be determined easily. A general condition for undrained behavior is given by:

f ¼ 0 ¼ ðb1 r11 þ b2 r22 þ b3 r33 Þ  cpf ;

ð16Þ

which relates the undrained pf to all the values of the applied external stresses r11, r22, and r33. Then, (1) is replaced by

0

e11

1

0

Sd11

B C B d @ e22 A ¼ B @ S12 e33 Sd 13

Sd12

Sd13

Sd22

Sd23

Sd23

Sd33

10 r11 CB B r22 b2 C AB @ r33 b3 p b1

1 C C C; A

ð17Þ

f

where pf = (b1r11 + b2r22 + b3r33)/c from (16). Rewriting this expression in its most intuitive form gives:

0

1 0 u e11 S11 B C B u @ e22 A ¼ @ S12 Su13

e33

Su12 Su22 Su23

Su13

10

1

r11 C u CB S23 A@ r22 A; u r33 S33

ð18Þ

where

Suij ¼ Sdij 

bi bj

c

ð19Þ

:

To distinguish drained (d) and undrained (u) compliances, we have added superscripts accordingly. Compliances without superscripts are always assumed to be drained, so Sij ¼ Sdij . 5. Analysis of random polycrystals composed of orthotropic elastic materials Now we briefly review random polycrystal analysis. We use the common convention of reducing elastic tensors to 6  6 matrices using the Voigt prescription:

11 1 0 S11 B C BS B 22 C B 12 C B B B 33 C B S13 C B B B C ¼ B B 23 C B C B B @ 13 A @ 12 0

S12

S13

S22

S23

S23

S33

10

S44 S55

r11 1 CB r C CB 22 C CB C CB r33 C CB C CB r C: CB 23 C CB C A@ r13 A S66 r12

ð20Þ

Again, the rij’s are the stresses, and the ij’s are the strains, for i,j = 1,2,3, corresponding respectively to spatial axes x, y, z. Compliances are symmetric, so S21 = S12, etc. And we consider only orthotropic symmetry – which is not the most general case, but one case often treated in practice. The elastic compliances S44, S55, S66 are those for the twisting shear strains: 23, 13, 12, and their related stresses. For isotropic elastic materials, the stiffness (inverse of compliance) moduli satisfy C11 = C22 = C33 = k + 2l, C44 = C55 = C66 = l, and C12 = C13 = C23 = k, where k and l are the two Lamé constants, and the isotropic bulk and shear moduli are given (but only in the isotropic case) by K = k + 2l/3 and G = l, respectively.

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

15

It is well-understood that, for such orthorhombic media, there are typically only three simple eigenvectors and eigenvalues, and these are the ones associated again with the twisting shear modes and the compliances, namely S44, S55, and S66, corresponding to rotations about spatial axes x, y, z, respectively. There are also three eigenmodes associated with the 3  3 submatrix in the upper lefthand corner of the full elastic matrix. However, these modes are not related in any simple way either to pure compression/extension or to pure shear modes. It follows that our analysis of moduli – such as effective bulk and shear modulus of polycrystals – requires understanding a rather complex relationship between two simpler concepts of a bulk modulus for pure compression or extension, and a shear modulus for one of the five (potentially) distinct shear modes of any elastic material. The resultant mixing of the modes is the main cause of all the difficulties we have when analyzing the average modal behavior of (by assumption) perfectly isotropic polycrystals (on average), and therefore makes necessary the use of the comparatively tedious methods now under discussion. Most studies of such systems are aimed at quantifying the behavior of random polycrystals – where the word ‘‘random’’ as used in this context usually implies that the polycrystals are composed of a large number of small crystallites typically oriented arbitrarily in space so the overall polycrystalline behavior is either isotropic or sufficiently close to isotropic for most purposes. The resulting effective isotropic constants can then be chosen to be effective bulk K and shear G moduli. The main goal of the rest of this paper will therefore be to find ways of estimating these effective bulk and shear moduli of such an overall isotropic average medium. The means we use to quantify these constants are first rigorous bounds, and then some approximate theoretical estimators. 6. Bounds and estimators

K SC

Watt (1979) for bounds on bulk modulus bounds K  HS and those of Middya and Basu (1986) for self-consistent estimates of bulk modulus can both be written in virtually the same form:

K HS ¼ K  þ

3B1 þ 2B2 ; 3 þ a ð3B1 þ 2B2 Þ

ð21Þ

where 

3B1 þ 2B2 ¼

9ðK V  K  Þ þ 2b ðd þ e  c Þ þ 3b2 X : 1  a b  9c ðK V  K  Þ þ D

ð22Þ

The coefficients a±, b±, c±, d±, e± are defined here (see the technical Appendix). The coefficients a±, b±, c± are defined below by (30)–(32), while X± is defined by (35) and (36). For self-consistency, we simply replace the subscripts or superscripts ± with either ⁄’s or SC symbols as desired, while also removing the HS (and/or PM) subscripts. The denominator of expression (22) is the same as the denominator of the first term in 

15B2 ¼





a  b þ b ð2d  2c  e Þ þ 3c ðd  c þ e Þ þ a b X 1  a b  9c ðK V  K  Þ þ D   3 1 1 1 þ 3ðG þ f Þ2    G þ f C 44 þ f C 55 þ f C 66 þ f

ð23Þ

and D± is defined in 

D ¼ b ðb þ 2c Þðc  d Þ  2e b c 

a b2 X 3

:

ð24Þ

The Voigt average for bulk modulus is

KV ¼

1 ½C 11 þ C 22 þ C 33 þ 2ðC 23 þ C 31 þ C 12 Þ: 9

ð25Þ

Similarly, the Voigt average for shear modulus is

GV ¼

1 ½C 11 þ C 22 þ C 33  C 23  C 31  C 12 þ 3ðC 44 þ C 55 þ C 66 Þ: 15

ð26Þ

For completeness, we also note that the corresponding Reuss averages (Skempton, 1954) for orthorhombic crystals are given by

1 ¼ ðS11 þ S22 þ S33 Þ þ 2ðS23 þ S31 þ S12 Þ KR

ð27Þ

15 ¼ 4ðS11 þ S22 þ S33 Þ  4ðS23 þ S31 þ S12 Þ þ 3ðS44 þ S55 þ S66 Þ; GR

ð28Þ

and

where the Sij’s are the compliance matrix elements, related to the stiffness matrix elements by the matrix equation S = C1.

16

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

The equation corresponding to bulk modulus K result (21) for the shear modulus G is given by

GPM ¼ G þ

B2 ; 1 þ 2b B2

ð29Þ

where PM indicates the contribution of Peselnick and Meister (1965), who were early evaluators of the pertinent HS bounds.    The Hashin–Shtrikman bounds themselves are then given exactly by K  HS  K PM and GHS  GPM . Again KV is the Voigt average of bulk modulus. Definitions of another useful shear factor Gveff depend on the specific crystal symmetry under consideration (see Berryman, 2005, 2011 for specifics). Factor B 2 itself was already defined in (23). Parameters a± and b± appear repeatedly above and can be related to Eshelby (1957) results by rewriting them in the forms:



1

a

¼ K  þ 4G =3

ð30Þ

and



1 ¼ G þ f : 2b

ð31Þ

Another combination of these two quantities frequently appearing in later formulas is

c 

a  3b 9

ð32Þ

:

Combining these definitions, we find that:

  1  B1 þ 2B2 =3 ðK  þ 4G =3Þ 1 ¼ ; K  þ 4G =3 K HS þ 4G =3

ð33Þ

which is valid for orthorhombic and some more symmetric crystal structures, such as hexagonal, tetragonal, and trigonal. Eq. (33) may then be compared to the analogous shear formula given by

GPM

1 1  B2 =ðG þ f Þ ¼ ; G þ f þ f

ð34Þ

which is valid for the same crystal symmetries, with (33) and (34) being analogous forms respectively for the bulk and shear moduli. As previously noted, these equations are for the upper and lower bounds K  HS on the bulk modulus. Resulting bounds are obtained when the constraints are optimal. Then the determinant of the matrix X± is given by



X  detðX  Þ ¼ X 11 X 22 X 33 þ 2X 12 X 23 X 13  X 11 X 23 ±

±

2

 2  2  X 22 X 13  X 33 X 12

ð35Þ

±

and we must have X = det (X )  0. Here X is a 3  3 positive- or negative-semi-definite matrix, as defined by

4 X 11 ¼ C 11  K   G ; 3 4  X 22 ¼ C 22  K   G ; 3 4  X 33 ¼ C 33  K   G ; 3

2 X 12 ¼ C 12  K  þ G ; 3 2  X 13 ¼ C 13  K  þ G ; 3 2  X 23 ¼ C 23  K  þ G : 3

ð36Þ

Vanishing of det (X) is a necessary requirement because then, and only then, have we found either the greatest lower bound, or the smallest upper bound. Middya and Basu (1986) have already shown these same equations can be used to determine self-consistent estimates, as well as the bounds. Self-consistent values are determined specifically by considering two conditions: B2 = 0 and 3B1 + 2B2 = 0. Both conditions must apply simultaneously if the self-consistency conditions are to be satisfied. And so, it follows that B1 = 0; but we never need to impose this condition on B1 separately. The self-consistency conditions are therefore given directly by:

K SC ¼ K 

and GSC ¼ G ;

ð37Þ ⁄

⁄

where the conditions that determine the values of K and G are exactly the ones that guarantee the two factors B1 and B2 vanish simultaneously. Although this simultaneity condition might sound hard to satisfy, it has been found in practice to be very easy to achieve by applying an iterative process wherein some convenient initial K0 and G0 values are chosen and substituted into (33) and (29) for the K± and G± values. The results that are next obtained for the left-hand-sides of both these equations then become the new trial values for K0 and G0. Repeating this process has always been found to converge quickly as long as some reasonably intelligent choices are made for those initial values of K0 and G0. So this part of the overall search procedure is not difficult to implement in practice.

17

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

Determining the HS bounds from this same set of equations is harder by comparison, but some tricks were developed in previous work (Berryman, 2011) that made the computational process more efficient. This work will be only briefly reviewed. In particular, a useful ‘‘shooting method of optimization’’ (which makes use of the computed values of the self-consistent estimates in order to speed up the search for HS bounds) was developed previously to streamline the algorithm used to produce the Hashin–Shtrikman bounds (Berryman, 2005, 2011; Hashin & Shtrikman, 1962; Watt, 1979). We refer the interested reader especially to Berryman (2011) for more details of this approach. The basic outline is this: Having already computed the pertinent self-consistent values, and knowing the Reuss and Voigt bounds on both K and G, we scan from some crude bounding values towards the SC point. We then take care to observe where sign changes in these functionals occur. In this way, the optimal bounds (those closest) to the SC estimates can be quickly determined. 7. Discussion of results for bounds and estimators A variety of other bounds and estimators are to be found in the literature (Kanaun & Levin, 2008; Kröner, 1977; Pham, 2006, 2011) in addition to the ones that we have presented here in detail. Although we consider (in this section) only seven specific examples in our analysis of these bounds and estimates, we are nevertheless able to show that a number of rather plausible hypotheses about possible relationships between and among the various bounds and estimates can be quickly excluded via specific counterexamples. The numerical values of the various bounds and estimates for the cases considered are summarized in Table 1. First, we have correctly anticipated that the self-consistent (SC) estimates always lie inside the Hashin–Shtrikman (HS) bounds. This expectation follows from the form of the equations for the self-consistent approximation, being based as they are – here for orthorhombic materials only – on essentially the same equations as those for the HS bounds. However, a further hypothesis (which might have seemed reasonable) that the SC estimates could always lie at or near the center of the HS bounding region is quickly disproven. Of these seven orthorhombic examples considered, six produce results close to an edge of the HS bounding box. Only for Danburite among the examples computed do the SC estimates appear to be nearly centered in the HS bounding box. In all six of the remaining cases, the SC estimates lie either on or very near to a point on the boundary of the HS bounding box. Another general observation (Berryman, 2011) is that the relative differences between the SC estimates and the VRH estimates for five of these orthorhombic materials are quantitatively small, i.e., typically being fractions of 1%. This fact suggests that, for applications not requiring very high precision estimates, the VRH estimates will often continue to be of some value. 8. Drained and undrained constants for fractured poroelastic polycrystalline systems In order to complete our generalization of the polycrystalline analysis to poroelastic systems, we need to find ways to compute both drained and undrained constants for these systems. Similarly, Berryman and Grechka (2006) discuss related issues for fractured systems without fluids; this case is essentially the same as that for drained poroelastic systems. If Gg and Kg are the bulk and shear grain moduli for the nonporous (polycrystalline) system (in our examples these would correspond

Table 1 Comparisons of bounds and estimates for the bulk (K) and shear moduli (G) of polycrystals (composites of randomly oriented crystals of uniform type) for seven orthorhombic elastic materials: a-Uranium, Aragonite (CaCO3), Danburite (CaB2Si2O8), Enstatite (MgSiO3), Forsterite (Mg2SiO4), Sulfur (S), and Topaz (Al2(F,OH)2SiO4). Estimators shown are: Reuss lower bound (R), Hashin–Shtrikman lower bound (HS), Self-consistent estimate (SC⁄), Hashin–Shtrikman upper bound (HS+), and Voigt (V) upper bound, respectively, for both the bulk (K) and shear moduli (G). All effective constants are in units of GPa (gigapascal). (Musgrave, 2003; Simmons & Wang, 1971; Bass, 1995; Carmichael, 1989). K, G

R

HS

SC⁄

HS+

V

a-Uranium

K G

111.3 80.7

112.5 83.6

112.7 84.1

113.1 84.9

114.6 87.9

Aragonite

K G

44.71 36.62

45.56 37.56

46.36 38.31

46.41 38.35

49.04 40.41

Danburite

K G

90.52 62.47

91.37 63.50

91.89 64.27

92.40 65.06

92.89 65.87

Enstatite

K G

107.29 75.18

107.65 75.52

107.83 75.70

107.83 75.71

108.29 76.15

Forsterite

K G

127.27 79.54

128.49 80.35

128.49 80.35

128.49 80.89

131.67 82.58

Sulfur

K G

17.56 6.17

18.76 6.61

18.85 6.64

18.87 6.66

20.60 7.22

Topaz

K G

166.19 113.54

167.37 114.73

167.46 115.06

167.73 115.09

168.67 116.00

18

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

to the self-consistent values of G and K Berryman, 1980a, 1980b; Mavko, Mukerji, & Dvorkin, 2009; Middya & Basu, 1986; Walsh, 1965), and Gd and Kd are the effective isotropic drained moduli when cracks or pores are present. Then we have respectively

1 Gd



1 4q ’ g 3 2 Gg

ð38Þ



1 ’ 2qg2 : Kg

ð39Þ

and

1 Kd

Here we have kept only the contributions from the second crack parameter g2 in the pertinent expansions, because it is known that g1 normally contributes less than a 5% correction, while g3 and higher order terms produce still smaller corrections. The most significant factor g2 itself can be determined independently within the noninteraction approximation (NIA) by

g2 ¼

8ð1  m0 Þð5  m0 Þ ; 15G0 ð2  m0 Þ

ð40Þ

where m0 = (3K0  2G0)/2(3K0 + G0) is the pertinent value of Poisson’s ratio. For these estimates, the nonporous values obtained with the self-consistent method are the preferred numerical values that should be used. Results for drained, undrained, and nonporous values are shown in Table 2. For easy reference, the values of g2 used in producing Table 2 are also displayed in Table 3. The final step in the process of finding all the pertinent parameters for granular poroelastic systems is the one where we obtain the undrained constants from the drained constants. The formula for doing so is well-known to be:

K uR ¼

K dR ; 1  aR B

ð41Þ

where aR is the Biot–Willis coefficient (5), B is Skempton’s second coefficient (6), and K dR is the drained bulk modulus, which we have calculated above for various materials in the preceding section. For a crack density q = na3, we need to multiply by (4p/3)  (aspect ratio) to determine the pertinent porosity. For q = 0.1, we find porosity / ’ 0.04 for an aspect ratio of 0.1. We assume here that the pore-fluid is water, so Kf = 2.25 GPa. In Table 2, we also display the ratio of the differences: drained to undrained over those from drained to nonporous. Results show that there are significant differences in these ratios, both when comparing results for shear and bulk moduli for the same material, and also when comparing the results across a set of different elastic materials. Finally, we present three Figures in order to emphasize the relationships between results of the nonporous elastic polycrystal analysis and the poroelastic drained and undrained results. Examples displayed here in Figs. 1–3, respectively, are for Aragonite, Danburite, and Topaz. For all three we have the polycrystal results for Hashin–Shtrikman bounds and the related self-consistent estimates in the upper righthand corner of each example. The corresponding drained (fluid not trapped) and undrained (fluid trapped) results appear in the lower left part of each figure. For the drained results, we see that the fluid has very little effect on the moduli, but would nevertheless have an effect on overall density, and therefore on the wave speeds in these cases. For the undrained cases, the fluid has only a very small effect on shear modulus, but a much larger effect on the bulk modulus. Table 2 Seven sets of polycrystal constants for drained, undrained, and nonporous examples in units of GPa (gigapascal). The final column provides a measure of the relative importance (via ratio) of the changes in values from drained to undrained (U  D) in comparison to those from drained to nonporous (N  D). K, G

Drained

Undrained

Nonporous (SC⁄)

UD ND

112.7 84.1

.227 .036

a-Uranium

K G

85.67 72.69

91.8 73.1

Aragonite

K G

36.27 33.22

38.4 33.4

46.36 38.31

.211 .035

Danburite

K G

70.0 56.1

79.1 58.3

91.89 64.27

.416 .269

Enstatite

K G

81.5 65.7

88.9 66.7

107.83 75.70

.281 .100

Forsterite

K G

94.5 69.9

104.7 71.3

128.49 80.35

.300 .134

Sulfur

K G

12.07 5.87

18.85 6.64

.593 .039

Topaz

K G

125.45 99.76

167.46 115.06

.377 .264

16.09 5.90 141.3 103.8

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

19

Table 3 Values of the second crack parameter g2 used in computations for Table 2.

g2 (GPa1) a-Uranium

0.014 0.030 0.017 0.015 0.014 0.149 0.010

Aragonite Danburite Enstatite Forsterite Sulfur Topaz

39

Shear Modulus G (GPa)

38

37

HS SC Undrained Drained

36

35

34

33 36

38

40

42

44

46

48

Bulk Modulus K (GPa) Fig. 1. Hashin–Shtrikman (HS) bounding box, and self-consistent (SC) estimates for effective elastic constants of polycrystalline Aragonite, plotted together with (poroelastic) undrained and drained results when the crack density is assumed to be q = 0.1. Influence of the pore-fluid is strongest in the undrained bulk modulus.

Shear Mudulus G (GPa)

70

65

HS SC Undrained Drained

60

55 65

70

75

80

85

90

95

Bulk Modulus K (GPa) Fig. 2. Hashin–Shtrikman (HS) bounding box, and self-consistent (SC) estimates for effective elastic constants of polycrystalline Danburite, plotted together with (poroelastic) undrained and drained results when the crack density is assumed to be q = 0.1. Poroelastic effects on shear modulus are small.

9. Overview and conclusions To summarize: Once the full set of elastic constants for an orthorhombic material constituting the random elements of the composite is known, the easiest isotropic estimators to compute are always the Voigt and Reuss averages for both bulk and shear moduli. Among the typical estimators most workers might consider computing, the next easiest and most useful quantities to compute are the self-consistent estimates. Self-consistent elastic constant estimation for polycrystals is robust,

20

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

116

Shear Modulus G (GPa)

114 112 110

HS SC Undrained Drained

108 106 104 102 100 98 120

130

140

150

160

170

Bulk Modulus K (GPa) Fig. 3. Hashin–Shtrikman (HS) bounding box, and self-consistent (SC) estimates for effective elastic constants of polycrystalline Topaz, plotted together with (poroelastic) undrained and drained results when the crack density is assumed to be q = 0.1. Poroelastic effects on shear modulus are again comparatively small.

by which statement we mean that the results follow easily from the pertinent equations. However, the optimal Hashin– Shtrikman bounds can be tedious to obtain directly, because careful checking for optimality is always a key issue. However, optimality can be more easily assured by first computing the self-consistent (CPA) estimates. These SC values can then be used to speed up the computation of the Hashin–Shtrikman bounds as shown previously (Berryman, 2011), and this approach was used to advantage here. There is also potential to redo some of the analysis presented using the popular differential scheme (Norris, 1985; Zimmerman, 1991). However, for the polycrystal analysis, it seems that Hashin–Shtrikman bounds and self-consistent estimates are very closely related estimators, and also somewhat easier to compute than similar results for the differential scheme. Nevertheless, one good reason to switch to a differential-scheme-based approach is likely to be for possible situations where it might be known that the microstructure (of actual earth materials to be studied) is better represented by the implicit microstructure of the differential scheme itself. With modern microtomography methods now being applied to rocks (Arns, Knackstedt, Pinczewski, & Garboczi, 2002), circumstances in which such additional structural information might become available for such analyses and will be more likely to arise in the near future. Acknowledgments Work performed under the auspices of the US Department of Energy, at the Lawrence Berkeley National Laboratory under Contract No. DE-AC02–05CH11231. Support was provided specifically by the Geosciences Research Program of the DOE Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences. Appendix A. Coefficients a±  e± The remaining constants appearing in (23) are given by

a ¼ X 11 þ X 22 þ X 33 ;



b ¼ X 12 þ X 23 þ X 13 ;  2  2  2  c ¼ X 11 X 22 þ X 22 X 33 þ X 33 X 11 ; d ¼ X 12 þ X 23 þ X 13 ;

ð42Þ

e ¼ X 12 X 13 þ X 13 X 23 þ X 23 X 13  X 11 X 23  X 22 X 13  X 33 X 12 ; where the X± matrix elements were defined in (36). (Note: The symbol ± always appears as a subscript for scalar quantities, except for the scalar Hashin–Shtrikman bounds themselves, where the bound label is used as a subscript. The symbol ± appears as a superscript for all quantities that are themselves matrix elements (therefore having additional subscripts), and for quantities that are combinations only of such matrix elements. For scalar quantities that are themselves combinations of scalars and also quantities derived from matrix elements, the subscript version is again used – except as already noted for the scalar bounds themselves.) References Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., & Garboczi, E. J. (2002). Computation of linear elastic properties from microtomographic images: Methodology and agreement between theory and experiment. Geophysics, 67, 1396–1405.

J.G. Berryman / International Journal of Engineering Science 72 (2013) 11–21

21

Bass, J. D. (1995). Elasticity of minerals, glasses, and melts. In A. J. Ahrens (Ed.), Mineral physics and crystallography (pp. 45–63). Washington, DC: American Geophysical Union. Berryman, J. G. (1980a). Long-wavelength propagation in composite elastic media. Journal of the Acoustical Society of America, 68, 1820–1831. Berryman, J. G. (1980b). Long-wavelength propagation in composite elastic media. II. Ellipsoidal inclusions. Journal of the Acoustical Society of America, 68, 1809–1819. Berryman, J. G. (2005). Bounds and self-consistent estimates for elastic constants of random polycrystals with hexagonal, trigonal, and tetragonal symmetries. Journal of the Mechanics and Physics of Solids, 53, 2141–2173. Berryman, J. G. (2010). Poroelastic measurement schemes resulting in complete data sets for granular and other anisotropic porous media. International Journal of Engineering Science, 48, 446–459. Berryman, J. G. (2011). Bounds and self-consistent estimates for elastic constants of polycrystals composed of orthorhombics or crystals with higher symmetries. Physical Review E, 83, 046130. Berryman, J. G. (2011). Poroelastic response of orthotropic fractured porous media. Transport in Porous Media, 93, 293–307. Berryman, J. G., & Grechka, V. (2006). Random polycrystals of grains containing cracks: Model of quasistatic elastic behavior for fractured systems. Journal of Applied Physics, 100, 113527. Biot, M. A., & Willis, D. G. (1957). The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics, 24, 594–601. Brown, R. J. S., & Korringa, J. (1975). On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics, 40, 608–616. Carmichael, R. S. (1989). CRC practical handbook of physical properties of rocks and minerals. Boca Raton, FL: CRC Press [p. 66]. Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of the Royal Society of London A, 241, 376–396. Gassmann, F. (1951). Über die Elastizität poröser Medien. Vierteljahrsschrift der Naturfoschenden Gesellschaft in Zurich, 96, 1–23. Gubernatis, J. E., & Krumhansl, J. A. (1975). Macroscopic engineering properties of polycrystalline materials: Elastic properties. Journal of Applied Physics, 46, 1875–1883. Hashin, Z., & Shtrikman, S. (1962). A variational approach to the theory of elastic behaviour of polycrystals. Journal of the Mechanics and Physics of Solids, 10, 343–352. Hill, R. (1952). The elastic behaviour of crystalline aggregate. Proceedings of the Physical Society of London, A65, 349–354. Hill, R. (1965). Theory of mechanical properties of fiber-strengthened materials. III. Self-consistent method. Journal of the Mechanics and Physics of Solids, 13, 189–198. Kanaun, S. K., & Levin, V. M. (2008). Self-consistent methods for composites: Volume 1 – static problems. Dordrecht, The Netherlands: Springer. Kanaun, S. K., & Levin, V. M. (2008). Self-consistent methods for composites volume 2 – wave propagation in heterogeneous materials. Dordrecht, The Netherlands: Springer. Kröner, E. (1977). Bounds for effective elastic moduli of disordered materials. Journal of the Mechanics and Physics of Solids, 25, 137–155. Mavko, G., Mukerji, T., & Dvorkin, J. (2009). The rock physics handbook: Tools for seismic analysis of porous media (2nd ed.). Cambridge, UK: Cambridge University Press [pp. 185–188]. Middya, T. R., & Basu, A. N. (1986). Self-consistent T-matrix solution for the effective elastic properties of noncubic polycrystals. Journal of Applied Physics, 59, 2368–2375. Musgrave, M. J. P. (2003). Crystal acoustics: Introduction to the study of elastic waves and vibrations in crystals. New York: Acoustical Society of America [pp. 278–281]. Norris, A. N. (1985). A differential scheme for the effective moduli of composites. Mechanics of Materials, 4, 1–16. Olson, T., & Avellaneda, M. (1992). Effective dielectric and elastic constants of piezoelectric polycrystals. Journal of Applied Physics, 71, 4455–4464. Peselnick, L., & Meister, R. (1965). Variational method of determining effective moduli of polycrystals: (A) hexagonal symmetry and (B) trigonal symmetry. Journal of Applied Physics, 36, 2879–2884. Pham, D. C. (2006). New estimates for macroscopic elastic moduli of random polycrystalline aggregates. Philosophical Magazine, 86, 205–226. Pham, D. C. (2011). On the scatter ranges for the elastic moduli of random aggregates of general anisotropic crystals. Philosophical Magazine, 91, 609–627. Reuss, A. (1929). Berechung der Fliessgrenze von Mischkristallen. Zeitschrift für Angewandte Mathematik und Mechanik, 9, 55. Simmons, G., & Wang, H. (1971). Crystal elastic constants and calculated aggregate properties: A handbook. Cambridge, MA: The MIT Press. Skempton, A. W. (1954). The pore-pressure coefficients A and B. Géotechnique, 4, 143–147. Soven, P. (1967). Coherent-potential model of substitutional alloys. Physical Review, 156, 809–813. Taylor, D. W. (1967). Vibrational properties of imperfect crystals with large defect concentrations. Physical Review, 156, 1017–1029. Thigpen, L., & Berryman, J. G. (1985). Mechanics of porous elastic materials containing multiphase fluid. International Journal of Engineering Science, 23, 1203–1214. Voigt, W. (1928). Lehrbuch der Kristallphysik. Leipzig: Teubner [p. 962]. Walsh, J. B. (1965). The effect of cracks on the compressibility of rock. Journal of Geophysical Research, 70, 381–389. Watt, J. P. (1979). Hashin–Shtrikman bounds on the effective elastic moduli of polycrystals with orthorhombic symmetry. Journal of Applied Physics, 50, 6290–6295. Willis, J. R. (1981). Variational and related methods for the overall properties of composites. In C.-S. Yih (Ed.), Advances in applied mechanics (pp. 1–78). New York: Academic Press. Zimmerman, R. W. (1991). Compressibility of sandstones. New York: Elsevier [pp. 110–127].