JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 116, B03207, doi:10.1029/2010JB007500, 2011
Pore effect on macroscopic physical properties: Composite elasticity determined using a two‐dimensional buffer layer finite element method model Akira Yoneda1 and Ferdous Hasan Sohag1 Received 25 February 2010; revised 24 December 2010; accepted 12 January 2011; published 29 March 2011.
[1] A porous system has distinct macroscopic proprieties that are very different from those of a nonporous matrix. There has been much interest in clarifying such characteristics. In the present study, a two‐dimensional buffer layer finite element method (FEM) model was developed to investigate the porosity effect on macroscopic elasticity without introducing assumptions or approximations. Buffer layers, introduced at the periphery of porous objects, were found to be effective in neutralizing irregularities caused by the boundary discontinuity. It is noted that the background concept of the buffer layer is similar to that of the window method used in many FEM analyses. In terms of the buffer layer FEM model, the porosity effect was systematically analyzed while considering anisotropy by changing the degree of porosity, aspect ratio of the pore, and elasticity of the matrix. Consequently, various systematic relations were found for two‐dimensional porosity effects. The relations are helpful in characterizing two‐dimensional cracks in a laboratory rock compression test. The two‐dimensional relations were further extended to three‐dimensional relationships to compare the present results with the results from previous models. Citation: Yoneda, A., and F. H. Sohag (2011), Pore effect on macroscopic physical properties: Composite elasticity determined using a two‐dimensional buffer layer finite element method model, J. Geophys. Res., 116, B03207, doi:10.1029/2010JB007500.
1. Introduction [2] Porous systems are found universally. Sandstone and pumice are porous rocks that accommodate underground water and hydrocarbons. Karst terrain is representative of porous structures on Earth [e.g., Hiltunen et al., 2007], and most meteorites and asteroids also have significant porosity [Britt et al., 2002]. Besides materials associated with Earth and planetary sciences, it is well known that some biological tissues, such as trunks and bone, are porous. Furthermore, there are many porous industrial goods, such as rubber sponges, metal foams, and semisintered ceramics. Therefore, evaluation of the macroscopic physical properties of porous objects has been studied intensively not only in Earth sciences but also in other disciplines, such as material engineering. [3] Among the various physical properties, macroscopic elasticity, or composite elasticity, has attracted special interest in the context of geophysical exploration for natural resources. Furthermore, the composite elasticity of porous objects is itself of fundamental interest because of nontrivial interactions among pores. Therefore, the porosity effect on composite elasticity has been investigated by employing
1 Institute for Study of the Earth’s Interior, Okayama University, Misasa, Japan.
Copyright 2011 by the American Geophysical Union. 0148‐0227/11/2010JB007500
various approaches including finite element methods (FEMs) [e.g., Mavko et al., 1998; Garboczi and Berryman, 2001]. [4] We note here the difference between the topic of the present work and poroelasticity [e.g., Mavko, 1980; Schmeling, 1985; Takei, 1998]. The study of poroelasticity focuses on complicated solid–liquid interactions; an important assumption made is that the solid frame and liquid frame are continuums intruding on each other. Needless to say, this setting is essential for underground water and oil. On the other hand, the pores in the present work are assumed to be isolated. Therefore, strictly speaking, the present work does not relate to poroelasticity in this context. However, we believe that the findings of this work are useful in the study of poroelasticity, in that they describe the macroscopic elasticity of the dry solid framework. [5] Although a single spherical pore in an isotropic elastic body can be analyzed explicitly by employing exact analytical solutions, interactions among multiple pores cannot be solved analytically. Therefore, many practical approaches, such as average or bound approaches and self‐consistent approaches, have been proposed and applied. The Voigt– Reuss bound [Hill, 1952] and the Hashin–Shtrikman bound [Hashin and Shtrikman, 1963] have been used as a first choice for constraining composite elasticity, because of their relatively simple expressions. However, it is well known that these bounds cannot constrain the composite elasticity of porous materials effectively; the lower sides of those bounds
B03207
1 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
B03207
[9] An important advantage of the FEM approach is that it can deal with anisotropy without difficulty. This means that we can extend the present FEM analysis for anisotropic cases in terms of not only physical properties of the constituent material but also the geometrical configuration of the pores and inclusions, even in three‐dimensional (3‐D) space. However, we restrict the present discussion to only the two‐dimensional (2‐D) porosity effect to attain a fundamental understanding of the porous system. [10] Finally, it is noted that we follow the Voigt notation for elasticity (stress s, strain ", stiffness elastic constant C, and elastic compliance constant S) throughout this study [e.g., Nye, 1985].
2. Procedures for 2‐D FEM Analysis: Introduction of the Buffer Layer Model Figure 1. Example of a multibuffer layer model for N = 5, or a 5 × 5 model in a square 1 m in length. Note that x and y directions are horizontal and vertical. In this case, the uniform compressive displacement (∼10−6 relative to the edge length) was applied vertically (in the y direction) from each end, although the present plot based on Eulerian coordinates does not show displacement of the object. The intensity of the von Mises stress is shown according to the side color bar (red, maximum; blue, positive minimum). The edges of the central black square (with 0.2 m length) are the observation lines for stress and displacement under the given condition. are uselessly too small owing to the infinitely small elasticity of pores [e.g., Watt et al., 1976]. [6] Therefore, the differential effective medium (DEM) method, a kind of self‐consistent approach, has been applied to the problem of porous materials [e.g., Budiansky, 1965; Hill, 1965; Wu, 1966; Walpole, 1969]. The reasonability of the DEM has been checked using experimental data for specimens of sintered hard material with suitable porosity [e.g., Berryman, 1995]. However, it has been difficult to prepare porous specimens with sufficient quality for evaluating the performance of these practical approaches. Thus, the FEM was introduced to provide unambiguous references [e.g., Garboczi and Berryman, 2001]. [7] The recent progress in using the FEM allows us to evaluate the composite elasticity of porous material numerically, without introducing any specific assumption or simplification. This superiority inspired us to conduct the present FEM analysis on the composite elasticity of porous material. As shown by the geometry model in Figure 1, the present study employed regular geometry free from randomness; this is important for searching a useful simple formula based on deterministic FEM results. [8] There have been extensive studies on “digital rock,” which is analogous of natural rock synthesized in computer by using FEM and other simulation techniques. One of the important concepts in current digital rock studies is how to generate random heterogeneous texture [e.g., Torquato, 2002; Chen et al., 2007; Nouy and Clement, 2010]. Therefore, the interpretations of digital rock simulation have been frequently stochastic rather than deterministic. In this sense, we claim that the approach taken in the present study is an alternative to current digital rock studies.
[11] While the wavelength of a teleseismic wave (∼1 Hz frequency) is of a few kilometers, the size of a crystal grain constituting typical rock is of a few millimeters. The discrepancy is as large as ∼106. Therefore, the size of the model geometry is a critical problem in conducting FEM analysis on rock texture [e.g., Garboczi and Berryman, 2001]. Using a quasi‐infinitely large model seems to be an ideal solution. However, handling too large a model is practically impossible owing to the limitations of computer capacity. [12] A breakthrough for the above problem came from the finding presented in Figure 1. Figure 1 shows an example of aligned pores in an isotropic elastic body. The example is referred to as a 5 × 5 buffer layer FEM model, in which uniform displacements were applied to the outermost edges as a boundary condition. We can recognize excellent periodicity within the geometry; the anomaly caused by the discontinuous boundary rapidly relaxes in peripheral layers. In other words, the peripheral layer is an extremely effective buffer for generating macroscopic periodicity of the stress–
Figure 2. An example of a shear deformation test. A quarter part of the 5 × 5 model is sufficient for the FEM analysis by taking account of symmetry. Unlike Figure 1, this also shows macroscopic deformation at a magnified scale for presentation. Colors indicate the von Mises stress with the same color bar as used in Figure 1. Blue and orange arrows specify off‐porosity and on‐porosity layers, respectively, in discussion of Figure 5.
2 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Figure 3. Two types of 7 × 7 analysis model with mesh plots. Note that only “3.5 × 3.5” pores are plotted because of the quarter‐part presentation. The difference between models is the position of the observation lines shown by red lines. The numbers of elements and degrees of freedom are ∼4700 and ∼20,000, respectively. It is confirmed that the present density of the mesh gives results exactly consistent with those determined using a rougher or finer mesh. The black line in Figure 3b is the observation line for stress shown in Figure 4. Horizontal and vertical axes are x and y coordinates, respectively, which are not shown because of the exact similarities with those of Figure 2.
strain state within the model geometry. This finding is the basic idea of the present FEM analysis of the composite elasticity of porous material. The edges of the central square are the observation lines of stress and displacement in the present FEM analysis. It is noted that the inside region of the observation lines can be considered as the counterpart of the window in the window method [e.g., Hain and Wriggers, 2008]. Although the outermost edge length is set as 1 m in the present work for simplicity, the absolute size of the geometry does not matter for the present purpose. This means that the ratio between the edge length and the diameter of the pore is essential. In this work, we conducted FEM analysis on COMSOL®, which is commercial FEM software. [13] Considering the symmetry along the central lines, the FEM analysis can be conducted only on a quarter part as shown in Figure 2; note that the horizontal and vertical central lines form a symmetric boundary in compressional deformation tests and an antisymmetric boundary in shear deformation tests. [14] Obviously, using a model with a larger number of multiple buffer layers would give us greater confidence in the present FEM analysis. Figure 3 shows examples of mesh plots for 7 × 7 models. The left‐hand and right‐hand models only differ in the position of the observation lines for stress and displacement. Note that the two 7 × 7 models and the 5 × 5 model give satisfactorily consistent results throughout. In this work, we decided to use the right‐hand 7 × 7 model for analysis. Using the model, it is confirmed that the present FEM analysis yields systematically consistent results for various mesh sizes of the FEM model between the porosities of 0.001 and 0.77. Note that the porosity of the present configuration is limited as p/4 (=0.785), at which point neighboring circles come into contact. [15] Three concepts must be taken into account in the FEM analysis, as pointed out by Garboczi and Berryman
B03207
[2001]: the finite size effect, spatial resolution, and statistical variation. The finite size effect was successfully overcome by introducing the buffer layer model. Figure 4 shows that sufficient periodicity was realized inside the buffer layer model. The spatial resolution was good enough because we carefully checked that the present FEM results were stable and consistent for different element sizes of the model (see Figure 3). The statistical variation was obviously satisfied because the present FEM model employed a regular lattice configuration of pores to find systematic relations. Therefore, the results we extracted from the FEM analysis are quite systematic and allow us to consider composite elasticity of porous and cracked objects. [16] Plane strain and plane stress are two typical cases for 2‐D analysis. In plane strain analysis, strains "3, "4, and "5, involving displacement in the z direction, are commonly assumed to be zero, while in plane stress analysis, stresses s3, s4, and s5 are commonly assumed to be zero. The strains and stresses are considered to occur in thick and thin 2‐D plates, respectively. Note that the result for the 3‐D problem lies between results for the 2‐D analyses. Therefore, the present 2‐D results are helpful in constraining the porosity problem in three dimensions. [17] After setting the model geometry as described above, uniform displacements were applied to the outer edges of the model geometry to generate a stress–strain relation inside the model. In this work, the matrix was first assumed to be elastically isotropic with a Young modulus Y =
Figure 4. Periodicity in the present FEM analysis. Solid line, s1 in the x direction compression test; dashed line, s2 in the y direction compression test; dotted line, s1 in the y direction compression test. Note that the horizontal axis corresponds to the x coordinate position on the observation line shown in Figure 3b. Geometry for the FEM analysis is exactly same as shown in Figure 3. The boundary conditions and the elastic properties of the model are consistent with the other FEM analyses conducted in the present work.
3 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
100 GPa; the Poisson ratio n was treated as a variable between 0 and 0.45. As an alternative boundary constraint, we can apply uniform stress to the edges instead of uniform displacement. It is confirmed that both constraints (uniform displacement and uniform stress) give consistent results throughout. It is noted that the displacement boundary condition was selected in the present work because of its simplicity. [18] Let us define u and v as displacements in x and y directions, respectively. In the 2‐D analysis we can classify three types of forced displacement on the outermost edges: (1) uniform forced displacement in the x direction, or u = −1 × 10−6 m on the x edge, while the model is free in the y direction, or there is no constraint on the y edge; (2) uniform forced displacement in the y direction, or v = −1 × 10−6 m on the y edge, while the model is free in the x direction, or there is no constraint on the x edge; and (3) coupled uniform forced displacement in the y direction, or v = 1 × 10−6 m on the x edge, and in the x direction, or u = 1 × 10−6 m on the y edge. It is noted that the ratio between the edge length and displacement is essential in an elastic problem. [19] We easily obtain the average strains from the resulting displacements Du and Dv along the observation lines shown in Figure 3b as Du Dv "1 ¼ ; "2 ¼ ; "6 ¼ x0 y0
Du Dv ; þ y0 x0
ð1Þ
ð2Þ
Note that we can make two sets of simultaneous equations according to the uniaxial compressions in the x and y directions. Consequently, we have four redundant equations for three unknown parameters, C11, C12, and C22; the three parameters were solved using the least squares method with an extremely small residual. For shear deformation, we have a single equation for C66: 6 ¼ C66 "6 :
ð3Þ
It is noted that other stiffness elastic constants, C13, C23, C44, and C55, cannot be evaluated in the plane strain model because they are related to strains "3, "4, "5 involving displacement in the z direction. [20] For plane stress analysis, we have similar relations: "1 ¼ S11 1 þ S12 2 "2 ¼ S12 1 þ S22 2 : "6 ¼ S66 6
[21] It is well known that the elastic wave velocity V is related to the elastic stiffness constant as according to V ¼
pffiffiffiffiffiffiffiffiffi C=;
ð5Þ
where r is density and C is the elastic stiffness constant relating to the propagation and polarization directions of the elastic wave. For this reason, plane strain analysis is more useful than plane stress analysis in geophysics. Therefore, in this article, we first focused on plane strain analysis and then conducted selected plane stress analysis for comparison. [22] Finally, we compare the present buffer layer model and the window method [e.g., Hain and Wriggers, 2008], which has been used frequently in recent FEM modeling. The concept of the window method is quite similar to that of the present buffer layer model in terms of neutralizing boundary discontinuity. However, the present buffer layer model differs in that (1) the observation line is set far more within the geometry compared with cases for models based on the window method, (2) we can recognize the outer peripheral region as a layer of pores owing to the systematically aligned geometry, and (3) the window method has been used to improve the robustness of average values obtained for a random geometry model.
3. Results of FEM Analysis
where the over bar for Du and Dv specifies their average and x0 and y0 are the positions of observation lines in x and y coordinates, respectively. We then formulate the simultaneous equations between the averaged stresses over the observed line and those averaged strains. For plane strain analysis, we have 1 ¼ C11 "1 þ C12 "2 : 2 ¼ C12 "1 þ C22 "2
B03207
ð4Þ
Thus, we can obtain the relations for the elastic compliance constants from the plane stress analysis. Similarly, other elastic compliance constants, S13, S23, S44, and S55, cannot be evaluated for the plane stress model.
3.1. Aligned Circular Pores [23] We first examined the case of aligned circular pores as shown in Figure 3. Figure 5a shows a selected result of the elastic moduli as a function of porosity at n = 0.3. Because of the alignment of the pores in 2‐D space, the symmetry of the model is equivalent to tetragonal symmetry in 3‐D space. Thus, the definition of the elastic moduli of C11 and C66 is obvious. CP (=(C11 + C12 + 2C66)/2) and CS (=(C11 − C12)/2) are the corresponding elastic moduli for longitudinal and shear waves, respectively, in the direction bisecting the x and y directions. Note that the direction is mentioned as [110] direction in this work from crystallographic analogy. We clearly see that each modulus converges to zero at the critical porosity of C = 0.785, or p/4, and the macroscopic elasticity is nearly isotropic below 5% porosity. It is noted that consistency between the two models shown in Figure 3 is satisfactory within the plotted line width. [24] Figure 5b shows the same results for normalized moduli C* (=C/C0), where C0 is that of the matrix itself or that at zero porosity. Above 10% porosity, significant anisotropy is recognized for both longitudinal and shear * and CS* is moduli. Furthermore, the difference between C66 * and C P* over a approximately three times that between C11 * * is significantly less than C *S, C11 wide range of While C66 is slightly greater than C P* in the intermediate porosity range. In the x and y directions, the alternation of off‐ porosity and on‐porosity layers (see blue and orange arrows in Figure 2) is configured in series for shear deformation and the alternation of off‐porosity and on‐porosity columns is configured in parallel for longitudinal deformation, while alternation in the [110] direction is less contrastive between off‐porosity and on‐porosity layers. It is understandable that * is controlled by the weak on‐porosity layer the smaller C66 * is controlled by the rigid off‐porosity and the larger C11
4 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
B03207
column. It is also understandable that the weak on‐porosity layer in the series configuration can control the phenomenon more effectively than the rigid off‐porosity column in the parallel configuration can. Although these concepts are qualitatively obvious even before conducting the FEM analysis, one of the advantages of the FEM analysis is that we can evaluate them quantitatively. [25] The relation between C* and was empirically fitted with sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 8 C*ð; 8Þ ¼ expfDð Þ8gf1 þ Eð Þ8 g 1 ; 8C 2
ð6Þ
where D(n) and E(n) are fitting parameters for the Poisson ratio n. D(n) is identical to the initial slope of C*, or Dij ð; Þ ¼ lim
!0
@Cij*ð; ; Þ : @
ð7Þ
While D(n) has clear physical meaning as the derivative of C* with respect to porosity, E(n) is introduced here just as a fitting parameter. The square root term in equation (6) obviously causes C* to converge to zero at = C. The results calculated with equation (6) are shown in Figure 5b with dashed lines; it is confirmed that equation (6) works well in the range n < 0.4. Figure 6a and 6b show D(n) and E(n), respectively, as functions of n. Note that their subscript suffixes correspond to those of C*. From the information in Figure 6b, we can easily reconstruct the relation between C* and . Although equation (6) seems satisfactory for practical use in the entire range of porosity, we would like to propose a more sophisticated functional form after introducing the ellipsoidal pore. [26] Figure 6a illustrates several remarkable findings, which are further discussed in section 3.2 with consideration of the aspect ratio of pores: (1) both D11 (=DP) and D66 (=DS) equal −3 at n = 0.25, (2) both D11 and D66 have integer values at n = 0, and (3) D66 increases linearly with the Poisson ratio n. Finding 1 may be due to the equality of the two Lame elastic constants at n = 0.25 for an isotropic elastic body. [27] From finding 1, we extract an interesting result for the relationship between the elastic wave velocity and porosity. Let us consider the normalized wave velocity V* = V/V0, where V is the wave velocity of porous material and V0 the original velocity of nonporous material. V* is also expressed as V* ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi C*=*;
ð8Þ
where r* is the normalized density defined as * ¼ =0 ¼ 1 :
ð9Þ
Thus, we derive dV * 1 dC* d* 1 þ dC*=d 1 þ Dð Þ ¼ ¼ ¼ d 2 d d 2 2
ð10Þ
and obtain dV * ¼ 1 or d8
V * / *
ð11Þ
Figure 5. (a) Elastic moduli versus porosity for a Poisson ratio of 0.3. Longitudinal and shear moduli are shown with solid and dashed lines, respectively. (b) Normalized plot for elastic moduli. Solid and dashed lines are equivalent to those in Figure 5a. Dashed lines are fitting results for each modulus calculated with equation (6). at around n = 0.25, regardless of whether we are considering a P or S wave. This finding is useful in geophysics because the typical Poisson ratio of rocks and minerals is about n = 0.25. Note that the present relation is only useful for a porous object with aligned cylindrical tubes. For a different type of porous object, such as a disk or sphere, a different relation is expected. For vertically aligned cracks, clear systematic relations are discussed in section 4. [28] It is also important to examine the case of a pore saturated with fluid; “dry” and “wet” pores are terms that specify the existence of fluid inside pores. Although the present FEM analysis is of a dry pore, knowledge of the composite elasticity of a wet pore is important in determining the acoustic properties of a reservoir of a natural resource such as water or hydrocarbons. [29] Considering the importance of this matter, further comprehensive study is needed to establish the composite
5 of 14
B03207
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
case, we expect 15% deviation from the dry value. Thus, equation (13) can be used to determine the acoustic properties of a not fully saturated wet pore system. 3.2. Aligned Ellipsoidal Pores [30] The case for an ellipsoidal pore was analyzed to evaluate the anisotropy caused by the shape of a pore. The aspect ratio of a pore is defined as a = a/b (see Figure 7), where a and b are the long and short axes of the ellipsoid in the x and y directions, respectively. In the present ellipsoidal pore analysis, special attention was paid to the initial slopes of D functions because findings 2 and 3 in section 3.1 are of sufficient interest to be investigated more precisely. [31] Figure 8 shows linear trends for D66 even at larger a. The trends are completely reproduced by the simple equation D66 ð; Þ ¼ Að1 Þ;
ð14Þ
where A ¼ 2 þ þ 1 ¼ 2 þ þ
1 :
ð15Þ
The parameter A, symmetric in terms of a and 1/a, is an important parameter characterizing the present observation. It is an unambiguous relation found employing an a posteriori approach on the present FEM results. In this sense, each FEM result is interpreted as a kind of experimental datum point to find a new systematic relation. [32] The simple functional form of equation (14) motivated us to find a similar formula for D11 and D22. After a little trial and error, we found 2A ðA1 Þ ; D11 ð; Þ ¼ 1 þ 21 þ ð1 2 Þ
Figure 6. (a) Results for D(n) or relation between D and Poisson ratio n. Note that D11 = DP and D66 = DS owing to the isotropy in an infinitely small porosity range. (b) Results for E(n). elasticity of the wet pore. However, we would like to propose a simple extension of equation (10) at the present time. Let us define the filling factor of fluid inside a pore as h and the effective density ratio between the matrix (r0) and fluid inside the pore (rf) as = h rf /r0. The normalized density is then * ¼ 1 ð1 Þ:
ð12Þ
D22 ð; Þ ¼ ð1 þ 2Þ þ
2A ðA1 Þ : ð1 2 Þ
ð16Þ
ð17Þ
Equations (16) and (17) were confirmed as being completely consistent with the present FEM results up to a = 30, as shown in Figures 8 and 9. Although these equations are not derived deductively, their elegant functional forms convince us of their legitimacy even at infinitely large a. The FEM analysis can be used not only to numerically simulate a complicated system but also to obtain a posteriori relations such as equations (14)–(17). Pores with larger a are interpreted as cracks. Therefore, the above equations can be applied to evaluate the effect of cracks in an elastic body as discussed in section 4.
If the pore is not fully saturated, or h < 1, we can ignore the change in the bulk modulus relating to there being fluid inside pore. In this case, we expect D not to change greatly. Therefore, substituting this expression into equation (8), we have dV * 1 þ Dð Þ 1 þ Dð Þ ¼ ¼ : d 2 2 2
ð13Þ
Thus, we expect a negatively larger dV*/d for a wet pore than that for a dry pore. If ≈ 0.3 as a representative
Figure 7. Geometry of an ellipsoidal pore. Axis lengths are a and b.
6 of 14
B03207
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Figure 8. FEM results for D66 (shown by asterisks) plotted on the solid line drawn using equation (14). The aspect ratio of the pore a is given at top of each plot. [33] It is also interesting to observe the behavior of D12. In the range of smaller porosity values, we can assume the 2‐D isotropic elasticity for porous objects as C11 ¼ C12 þ 2C66
ð18Þ
from the analogy of 3‐D isotropic elasticity using Lame elastic constants of l (=C12), m (=C66), and l + 2m (=C11). Therefore, in this porosity range, we expect @C12 @C11 @C66 ¼ 2 : @ @ @
ð19Þ
The relation can be modified by substituting equations (14)–(17) to give
1 2A ðA1 Þ 0 0 þ 2C66 D12 ð; Þ ¼ 0 C11 A þ A ð1 Þ : ð1 2 Þ C12 ð20Þ
This formula corresponds completely to the FEM results as shown in Figure 10. [34] After finding these characteristic relations regarding the initial slope of C*, we tried again to interpret the global trend of C* and found the simple functional forms * ¼ expðD11 Þ 1 D11 2 ; C11
ð21Þ
* ¼ expðD22 Þ 1 D22 2 ; C22
ð22Þ
* ¼ expðD66 Þ 1 þ D66 2 : C66
ð23Þ
The performance of these formulas is shown in Figure 11. At a glance, calculations using the above equations are generally consistent with the FEM results within a few percent, with some exceptions. In particular, the agreement appears perfect at n = 0, where the initial slopes of D11, D22, and D66 solely depend on the geometric factor a. [35] Although the above functions cannot cover the terminal range of porosity near C (=p a−1/2/4) unlike equation (6), it is beneficial to find practical functional forms approximating the FEM results better than equations (21)–(23) do in the intermediate range of . We see that equations (21)–(23) become inconsistent with the FEM results with increasing n. Through trial and error, we found the following functions as a practical approximation of the FEM results. * ¼ expðD11 Þ 1 1 1 C11 * C22 ¼ expðD22 Þ 1 1 1 * ¼ expðD22 Þ 1 þ 1 1 C66
2ð1 2 Þ
D11 2 ; ð24Þ
2ð1 2 Þ 2ð1 2 Þ
D22
2
; ð25Þ
D66 2 : ð26Þ
These functional forms are confirmed as being consistent with the FEM results within 10% at least for the plotted range in Figure 12. [36] After finding the relations on the basis of the plane strain FEM analysis, plane stress FEM analysis was also conducted to compare the two representative settings in 2‐D
7 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Figure 9. FEM results for D11 and D22 (shown by asterisks) plotted on the solid line drawn using equations (16) and (17). The aspect ratio of the pore a is given at top of each plot. (a) Plots for D11. (b) Plots for D22.
8 of 14
B03207
B03207
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Figure 10. FEM results for D12 (shown by asterisks) plotted on the solid line drawn using equation (20). The aspect ratio of the pore a is given at top of each plot.
analysis. In plane stress analysis, we define the initial slope of the function for elastic compliance constants as @Sij*ð; ; Þ ; !0 @
dij ð; Þ ¼ lim
ð27Þ
where the definition of S* (=S/S0) is similar to that of C*. From analogy with plane strain analysis, we find d66 ð; Þ ¼
A ; 1þ
3.3. The 3‐D Extension Based on the Present 2‐D Analysis [38] For comparison between the present 2‐D results and the existing 3‐D estimations, let us consider the simplest extension of the present 2‐D FEM results to an isotropic 3‐D structure. Owing to the symmetry of the present model, we can assume orthotropy for elastic constants. Thus, we can write a 6 × 6 stiffness matrix as 0
ð28Þ
d11 ð; Þ ¼ 1 þ 21 ;
ð29Þ
d22 ð; Þ ¼ ð1 þ 2Þ:
ð30Þ
[37] The relations are completely consistent with the FEM analysis as shown in Figure 13 up to a = 20. The similarity between D and d values again convinces us of the reasonability and accuracy of the present FEM analysis. We see that d11 and d22 are independent of the Poisson ratio. It is understandable that the plane stress condition does not constrain "3 determined by the Poisson ratio. Although the present plane stress FEM analysis became unstable over a = 20, we believe that equations (28)–(30) can be extrapolated for a range of much higher a values. The various findings based on the plane stress FEM analysis are useful in constraining the mechanical properties of perforated thin plates. This setting seems to be more important in material engineering than in Earth and planetary sciences.
C11 B C12 B B C13 C¼B B 0 B @ 0 0
C12 C22 C23 0 0 0
C13 C23 C33 0 0 0
0 0 0 C44 0 0
0 0 0 0 C55 0
1 0 0 C C 0 C C: 0 C C 0 A C66
ð31Þ
It is noted that C11, C12, C22, and C66 are determined through the present 2‐D FEM analysis. Other constants (C33, C13, C23, C44, and C55) are assumed simply as 0 Cmn ¼ ð1 ÞCmn ;
ð32Þ
0 where Cmn is the corresponding stiffness for the isotropic nonporous solid matrix material. Equation (32) must be the exact relation for C33, while it is the simplest assumption for the other constants. [39] We suppose an isotropic body is composed of porous objects characterized by equations (31) and (32). We then invoke the Voigt and Reuss averaging procedures for iso-
9 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
B03207
tropic elasticity (see Appendix A). The results for bulk modulus K and rigidity G are shown in Figure 14, together with those obtained using the reciprocal theory for a 3‐D honeycomb structure of dry circular tubes [e.g., Walsh, 1965; Mavko, 1980] and the DEM theory for the isotropically distributed needle‐type pore [e.g., Berryman, 1995]. Note that the comparison is only shown between n = 0 and n = 0.3 because the Reuss averages of K and G are divergent at n = 0.4. The difference between the reciprocal theory and the DEM may be attributed to the difference in assumed geometry. Basically, the present Voigt and Reuss averages are closer to the results of the reciprocal theory and the DEM, respectively. Especially, the Reuss average of rigidity is surprisingly consistent with that obtained using the DEM. Rigidities among the four models are consistent up to 10% porosity. On the other hand, the consistency for bulk modulus is only up to porosity of a few percent for a Poisson ratio of 0.2–0.3, corresponding to natural rocks. This contrast may be due to the appropriate assumption of equation (32) for C44 and C55 and an inappropriate assumption for C13 and C23. The functional form for C13 and C23 may be much more convex downward than specified as equation (32), which would decrease the Voigt and Reuss average values of the bulk modulus. [40] The DEM results have been compared with experimental data, although experimental data tend to have significant error [e.g., Berryman, 1995]. It is noted that such errors are mainly due to the difficulty in preparing well‐ controlled porous specimens. Therefore, the FEM analysis can provide an alternative measure of checking the performance of the DEM. The present FEM analysis suggests the DEM is more compatible with the Reuss average than with the Voigt average. [41] Because of the 3‐D isotropic nature of the DEM, the present 3‐D extension may be a good tentative measure for
Figure 11. FEM results (solid lines) versus simplified “experimental” functions (dotted lines), given as equations (21)–(23), for three aspect ratios a = 1, 4, 8; as shown in Figure 11a (top left), a is easily distinguished using the difference in maximum porosities (horizontal axis). Results for (a) D11, (b) D22, and (c) D66.
Figure 12. FEM results (solid lines) versus simplified experimental functions (dotted lines), given as equations (24)–(26), for three aspect ratios a = 1, 4, 8; as shown in Figure 11a (top left). Only the results for n = 0.4 are shown because the consistency between solid and dotted lines is better in cases of n < 0.4.
10 of 14
B03207
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Figure 13. FEM results for d values (shown by asterisks) are plotted on the solid line drawn using equations (28)–(30). The aspect ratio of the pore, a, is 20. evaluating the DEM. Obviously, future 3‐D FEM analysis is required for an exact comparison with the DEM.
oriented vertical cracks (Figure 15c) and linear effect of each crack in the low porosity range. Thus, we introduce a mean initial slope
4. Implication for the Acoustic Velocity in a Crack System
Dm ¼
[42] It is well known that the acoustic velocity decreases with the growth of a crack in a rock compression test [e.g., Schubnel et al., 2003]; this concept has been frequently analyzed theoretically and experimentally. The present results regarding D values are quite useful for interpreting the crack effect on acoustic properties. In a laboratory uniaxial compression test, the initial crack system is basically tensile or vertically aligned like the present 2‐D FEM model (Figure 15a). [43] First, let us examine the anisotropy of longitudinal elastic constants due to vertically aligned cracks. Figure 15b shows the concept of the vertically aligned cracks in a laboratory compression test. Using equations (16) and (17), we easily obtain the ratio between D11 and D22, or R1 = D22/D11 as shown in Figure 16a. There are several points of interest we note from the plots: (1) R1 increases monotonously with increasing a, and converges to a finite value except for n = 0. (2) R1 decreases monotonously with increasing n. (3) There is no further increase in R1 above a > 100 if the Poisson ratio n > 0.2. (4) There is the possibility of specifying the crack aspect ratio a from R1 and n. [44] From points 1 and 2, we conclude that a high Poisson ratio can neutralize geometric anisotropy due to an ellipsoidal pore. Related to point 3, we see R1 ≈ 9 for cracks in rock if a > 100 and n ≈ 0.25. [45] Furthermore, let us compare longitudinal and shear elastic constants. It seems reasonable to assume the randomly
D11 þ D22 2
ð33Þ
for a mean longitudinal elastic constant in the lateral direction of the rock specimen in a uniaxial compression test. Figure 16b shows the ratio R2 = Dm/D66 as a function of a and n. We note several points of interest from the plots: (1) R2 increases monotonously with increasing a and n. (2) There is no further increase in R2 above a > 100. (3) Practically, R2 > 1 for cracks in rock unless simultaneously a < 10 and n < 0.1. (4) It is possible to specify the crack aspect ratio a from R2 and n. These concepts are applicable for laboratory experimental data. [46] Schubnel et al. [2003] measured the elastic wave velocity in a rock compression test. From their Figure 3, we can evaluate the aspect ratio of cracks in their experiments. First, they did not observe significant horizontal anisotropy. This suggests a randomly oriented crack distribution as shown in Figure 15c. Therefore, we can calculate R2 from experimental data as R2 ¼
ðDVP =VP0 Þ2 ðDVS =VS0 Þ2
2:2;
ð34Þ
where VP and VS are P and S wave velocities, respectively. VP0 and VS0 specify the reference velocities just before cracks grow in the rock sample. The results suggest that the crack aspect ratio a should be higher than 100 if the Poisson ratio n ≈ 0.3. This preliminary comparison with their
11 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
B03207
[48] Furthermore, we presented the potential utility of the present results for clarifying the geometrical properties of cracks, which has been difficult to clarify using any other method. Therefore, this work may open new research methodology in rock mechanics, which is important not only in civil engineering but also in geophysics. [49] The outcome of the present study is recognition of the importance of porous systems in Earth and planetary sciences. Meteorites and asteroids have particularly high porosity; in this case, the porosity effect on shock compression must be evaluated precisely using elastic–plastic FEM analysis. Therefore, extensive FEM analysis of porous systems is not only fascinating but also indispensable both in Earth and planetary sciences.
Appendix A: Voigt, Reuss, and Hill Averages for Polycrystalline Aggregate [50] Note that the elastic stiffness constant matrix, Cmn in (31), is for a single crystal with orthorhombic symmetry. The corresponding matrix of elastic compliance constants, Smn, is given as the inverse matrix of Cmn. The averages are then given as C11 þ C22 þ C33 ; 3 C44 þ C55 þ C66 ; C¼ 3 A¼
C11 þ C22 þ C33 ; 3 C44 þ C55 þ C66 ; c¼ 3 a¼
B¼
C12 þ C23 þ C31 ; 3 ðA1Þ
b¼
C12 þ C23 þ C31 ; 3 ðA2Þ
Figure 14. A 3‐D extension of the present 2‐D composite elasticity compared with results for other theoretical models. (a) The results for isotropic bulk modulus. Solid lines are for the present scheme based on the FEM results. Among them, the upper and lower curves correspond to Voigt and Reuss averages, respectively. Dashed lines are the results from reciprocal theory (upper curve) and the DEM (lower curve) for reference. (b) The results for isotropic rigidity. Legend is the same as that in Figure 14a.
experimental data convinces us of the reasonability of our present equations derived from the FEM analysis.
5. Conclusion [47] The present FEM study is based on the multibuffer layer model newly founded in the preliminary stage of the present work. The reliability of FEM analysis of the multibuffer layer model is unequivocally supported by the present analytical results. The present FEM analysis is proven to have potential in suggesting various exact relations in an a posteriori approach. Therefore, the present 2‐D models must be extended to various 3‐D models for more direct comparison with real systems.
Figure 15. Schematic view of the rock compression test and the two representative cross sections with cracks. A couple of anvils (dark gray) compress a rock sample (light gray). (a) Lateral view of the rock compression test. Vertical thin lines show tensile cracks vertically distributed. The dashed arrow specifies the propagation direction of the elastic wave discussed here. (b) Schematic cross section of the aligned cracks. Note that the cracks are shown consistently with a lateral view. (c) Case of randomly oriented cracks.
12 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
B03207
Figure 16. Ratio between initial slopes of D functions. The curvatures are calculated as a function of a and n using equations (14)–(17) and equation (33). Poisson ratios n are given on the right axis. (a) Variation in R1 (=D22/D11) as an indicator of anisotropy in longitudinal elastic constants. (b) Variation in R2 (=Dm/D66) as the ratio between Dm for the longitudinal elastic constant and D66 for the shear elastic constant. Thus, Voigt and Reuss averages for bulk modulus K and rigidity G are A þ 2B ; 3
GV ¼
A B þ 3C ; 5
ðA3Þ
1 ; 3ða þ 2bÞ
GR ¼
5 : 4a 4b þ 3c
ðA4Þ
KV ¼
KR ¼
Obviously, the subscripts “V” and “R” refer to Voigt and Reuss averages, respectively. The Hill average frequently used is the average between the Voigt and Reuss averages [Hill, 1952]. Meister and Peselnick [1966] derived useful equations for other crystallographic symmetries.
[51] Acknowledgments. This work was supported by a Grant‐in Aid for Scientific Research (19204044) from the Japan Society for the Promotion of Science and by a grant from the COE‐21 program. The authors would like to express their gratitude to the Associate Editor and two anonymous reviewers for critical reading of the manuscript and numerous constructive comments.
References Berryman, J. P. (1995), Mixture theories for rock properties, in Rock Physics and Phase Relation: A Handbook of Physical Constants, AGU Ref. Shelf, vol. 3, edited by T. J. Ahrens, pp. 205–228, AGU, Washington, D. C. Britt, D. T., D. Yeomans, K. Housen, and G. Consolmagno (2002), Asteroid density, porosity, and structure, in Asteroids III, edited by W. F. Bottke Jr. et al., pp. 485–500, Univ. of Ariz. Press, Tucson.
13 of 14
B03207
YONEDA AND SOHAG: COMPOSITE ELASTICITY USING FEM MODEL
Budiansky, B. (1965), On the elastic moduli of some heterogeneous materials, J. Mech. Phys. Solids, 13, 223–227, doi:10.1016/0022-5096(65) 90011-6. Chen, S., Z. Q. Yue, and L. G. Tham (2007), Digital image based approach for three‐dimensional mechanical analysis of heterogeneous rocks, Rock Mech. Rock Eng., 40, 145–168, doi:10.1007/s00603-006-0105-8. Garboczi, E. J., and J. G. Berryman (2001), Elastic moduli of a material containing composite inclusions: Effective medium theory and finite element computations, Mech. Mater., 33, 455–470, doi:10.1016/S0167-6636 (01)00067-9. Hain, M., and P. Wriggers (2008), Numerical homogenization hardened cement paste, Comput. Mech., 42, 197–212, doi:10.1007/s00466-0070211-9. Hashin, Z., and S. Shtrikman (1963), A variational approach to the elastic behavior of multiple minerals, J. Mech. Phys. Solids, 11, 127–140, doi:10.1016/0022-5096(63)90060-7. Hill, R. (1952), The elastic behavior of crystalline aggregate, Proc. Phys. Soc. London, Ser. A, 65, 349–354. Hill, R. (1965), A self‐consistent mechanics of composite materials, J. Mech. Phys. Solids, 13, 213–222, doi:10.1016/0022-5096(65)90010-4. Hiltunen, D. R., N. Hudyma, T. P. Quigley, and C. Samakur (2007), Ground proving three seismic refraction tomography programs, Transp. Res. Rec., 2016, 110–120, doi:10.3141/2016-12. Mavko, G. M. (1980), Velocity and attenuation in partially molten rocks, J. Geophys. Res., 85, 5173–5189, doi:10.1029/JB085iB10p05173. Mavko, G., T. Mukerji, and J. Dvorkin (1998), The Rock Physics Handbook: Tools for Seismic Analysis in Porous Media, Cambridge Univ. Press, Cambridge, U. K. Meister, R., and L. Peselnick (1966), Variational method of determining effective moduli of polycrystals with tetragonal symmetry, J. Appl. Phys., 37, 4121–4125, doi:10.1063/1.1707986. Nouy, A., and A. Clement (2010), eXtended Stochastic Finite Element Method for the numerical simulation of heterogeneous materials with
B03207
random material interfaces, Int. J. Numer. Methods Eng., 83, 1312–1344, doi:10.1002/nme.2865. Nye, J. F. (1985), Physical Properties of Crystals: Their Representation by Tensors and Matrices, Oxford Univ. Press, New York. Schmeling, H. (1985), Numerical models on the influence of partial melt on elastic, anelastic and electric properties of rocks. Part I: Elasticity and anelasticity, Phys. Earth Planet. Inter., 41, 34–57, doi:10.1016/0031-9201(85) 90100-1. Schubnel, A., O. Nishizawa, K. Masuda, X. J. Lei, Z. Xue, and Y. Gueguen (2003), Velocity measurements and crack density determination during wet triaxial experiments on Oshima and Toki granites, Pure Appl. Geophys., 160, 869–887, doi:10.1007/PL00012570. Takei, Y. (1998), Constitutive mechanical relations of solid‐liquid composites in terms of grain‐boundary contiguity, J. Geophys. Res., 103, 18,183–18,203, doi:10.1029/98JB01489. Torquato, S. (2002), Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, New York. Walpole, L. J. (1969), On the overall elastic moduli of composite materials, J. Mech. Phys. Solids, 17, 235–251, doi:10.1016/0022-5096(69)90014-3. Walsh, J. B. (1965), The effect of cracks on the compressibility of rock, J. Geophys. Res., 70, 381–389, doi:10.1029/JZ070i002p00381. Watt, J. P., G. F. Davies, and R. O’Connell (1976), The elastic properties of composite materials, Rev. Geophys. Space Phys., 14, 541–563, doi:10.1029/ RG014i004p00541. Wu, T. T. (1966), The effect of inclusion shape on the elastic moduli of a two‐phase material, Int. J. Solids Struct., 2, 1–8, doi:10.1016/0020-7683 (66)90002-3. F. H. Sohag and A. Yoneda, Institute for Study of the Earth’s Interior, Okayama University, 827 Yamada, Misasa, Tottori 682‐0193, Japan. (
[email protected]‐u.ac.jp)
14 of 14
