EXAMPLES

Grechka and Kachanov (2006a,b) have generated two distinct numerical data sets for crack densities $\rho = 0.05$, 0.10, 0.15, and 0.20. (Crack density is defined as $\rho = na^3$, where n is the number density per unit volume. Radius of a typical crack is a. If there is a distribution of crack sizes, then na³ is replaced by an appropriate average value.) Cracks are all penny-shaped, but not necessarily flat. We assume that the aspect ratio is small (but finite, i.e., nonzero) and, therefore, we can ignore it here.

The two models considered have very different Poisson's ratios for the isotropic background media: (1) $\nu_0 = 0.00$ and (2) $\nu_0 = 0.4375$. We will call these two models, respectively, the first model and the second model. The first model corresponds to a rigid medium very stiff under shear, while the second model corresponds to a medium much more compliant in shear. The two models approximately bracket the expected range of behavior in earth systems, for which a typical value might be $\nu_0 = 0.4$ in reservoir rocks. The first model has background stiffness matrix values C₁₁ = C₂₂ = C₃₃ = 13.75 GPa, C₁₂ = C₁₃ = C₂₃ = 0.00 GPa, and C₄₄ = C₅₅ = C₆₆ = 6.875 GPa. Bulk modulus for this model is therefore K₀ = 4.583 GPa and shear modulus is G₀ = 6.875 GPa. The second model has stiffness matrix values C₁₁ = C₂₂ = C₃₃ = 19.80 GPa, C₁₂ = C₁₃ = C₂₃ = 15.40 GPa, and C₄₄ = C₅₅ = C₆₆ = 2.20 GPa. Bulk modulus for this model is therefore K₀ = 16.86 GPa and shear modulus is G₀ = 2.20 GPa. The second model also corresponds to a background material having compressional wave speed V_p = 3 km/s, shear wave speed V_s = 1 km/s, and mass density $\rho_m = 2200.0$ kg/m³.

The numerical modeling of Grechka and Kachanov (2006a,b) for the first model resulted in 110 examples: 23 for $\rho = 0.05$, 21 for $\rho = 0.10$,35 for $\rho = 0.15$, and 31 for $\rho = 0.20$. Of these 110 cases, we used 79 (selected randomly for the three smaller crack densities): 8 for $\rho = 0.05$, 16 for $\rho = 0.10$, 24 for $\rho = 0.15$, and all 31 for $\rho = 0.20$.

The numerical modeling of Grechka and Kachanov (2005,2006a,b) for the second model resulted in 100 examples: 25 each for all four cases $\rho = 0.05, 0.10, 0.15, 0.20$. Of these 100 cases, we used 73 (selected randomly for the three smaller crack densities): 8 for $\rho = 0.05$, 16 for $\rho = 0.10$, 24 for $\rho = 0.15$, and all 25 for $\rho = 0.20$.

For both numerical data sets, the output of the modeling for each example was a stiffness matrix. All these matrices were observed to be essentially orthotropic in character (i.e., depending significantly only on the nine stiffness constants C₁₁, C₂₂, C₃₃, C₁₂ = C₂₁, C₁₃ = C₃₁, C₂₃ = C₃₂, C₄₄, C₅₅, C₆₆ -- meaning that the remaining off-diagonal coefficients are negligible). In order to produce results useful for graphical comparisons, these orthotropic matrices were used to compute the Voigt-Reuss-Hill averages of bulk and shear moduli. So the Voigt and Reuss averages K_V and K_R were computed for each matrix, and arithmetic average K_VRH = (K_R + K_V)/2 was obtained as a measure of the bulk modulus for the stiffness matrix. Then, a similar calculation was done for the shear modulus G_VRH. Formulas for all these averages are given for orthotropic elastic stiffness matrices by Watt (1979). The results (VRH) are plotted as red $\times$'s in Figures 2-9, where the numerical results are also compared to results from the differential scheme (DS), the non-interacting approximation (NI), and various bounds and estimates based on random polycrystals of cracked grains (R,HS^-,SC,HS⁺,V).

We also show two other estimates: SMP for ``sample,'' and GR for ``grain.'' Since these computations were all done using displacement boundary conditions, we have made use of results from Huet (1990) on sampling -- the pertinent result being that effective overall stiffness C_ij^* of a random medium (one having a well-defined statistical ensemble associated with it) satisfies $C_{ij}^* \le C_{ij}^{SMP} \le C_{ij}^V$, where  

$$C_{ij}^{SMP} = \frac{1}{N}\sum_{n=1}^N C_{ij}^{n}$$ (1)

and where the C_ijⁿ are subsamples of the collection of local stiffnesses in the ensemble. For present purposes, we can choose to think of the individual stiffness matrices found in the numerical experiments as being just such representatives selected randomly from the overall ensemble. Then, the SMP value we use is the one based on the values N = 8, 16, 24, etc., for the cases considered in the numerical experiments. Clearly, the value we obtain this way for C_ij^SMP is merely an estimate of the true average value, based on our rather limited statistics. But we find that the fluctuations in these averages are quite small already with the stated sample numbers used here, and so we conclude that the approach is in fact useful even for such apparently small sample statistics. Furthermore, the plotted values K_SMP and G_SMP are the Voigt averages for bulk and shear modulus obtained from the stiffness average C_ij^SMP. For comparison purposes, we also provide the values K_GR and G_GR, which are the Reuss averages associated again with stiffness average C_ij^SMP. These values have no special theoretical significance in terms of the analysis of Huet (1990), as the more appropriate bounding values from below should be obtained -- not from the stiffness obtained for displacement boundary conditions but -- instead for the compliance determined in a numerical experiment applying traction boundary conditions. But for the available data sets, we did not have this additional information. One further motivation, however, for presenting these values K_GR and G_GR is that the stiffness C_ij^SMP is in some fairly precise sense the stiffness of an ``average'' cracked grain in the overall polycrystal model. The bulk modulus of this ``average'' cracked grain is given precisely by the value K_GR, but there is no corresponding statement that can be made about the shear modulus of this same average grain. So G_GR is just a special estimate, or heuristic value, that can be easily computed; and its significance should therefore not be overinterpreted.

One class of approximations that has been found very useful for analysis of these cracked systems by Grechka and Kachanov (2006a,b) is the non-interaction approximation (NIA). Since our results will be presented in the form of plots of effective bulk and shear moduli, it is useful to consider the corresponding NIA formulas for bulk and shear moduli. Zimmerman (1991) gives such formulas, and in our present notation these formulas are:  

$$\frac{K_0}{K_{NI}} = 1 + \rho\frac{16(1-\nu_0^2)}{9(1-2\nu_0)}$$ (2)

for effective bulk modulus K_NI, and  

$$\frac{G_0}{G_{NI}} = 1 + \rho\frac{32(1-\nu_0)(5-\nu_0)}{45(2-\nu_0)}$$ (3)

for effective shear modulus G_NI. The host medium has bulk modulus K₀, shear modulus G₀, and Poisson's ratio $\nu_0 = (3K_0-2G_0)/2(3K_0+G_0)$. Again, the crack density $\rho = na^3$,where n is the number of cracks per unity volume, and a is the radius of the (assumed) penny-shaped cracks. In particular, we note that in a plot of inverse bulk modulus and/or inverse shear modulus versus crack density $\rho$, the NIA results are just straight lines. This general feature of NIA suggests that it is most fruitful to construct our plots in this way in order to distinguish easily whether the results are behaving according to NIA predictions -- or not. In fact, we find that, although the NIA gives good agreement for some of the numerical results, in general there are deviations from NIA, and that the polycrystal of cracked grains model gives a better representation of the numerical results.

Another method that we use here involves a crack-influence decomposition method of Sayers and Kachanov (1991); see Appendix B for details. When all cracks in the system have the same vertical (z-)axis of symmetry, then the cracked/fractured system is not isotropic, and we have the compliance correction matrix  

$$\Delta S_{ij} = \rho\left(\begin{array} {cccccc} 0 & 0 & \et... ... \cr & & & & 2\eta_2 & \cr & & & & & 0 \cr\end{array}\right).$$ (4)

(In fact, we will show later that $\eta_1$ corrections in (4) are usually negligible compared to the $\eta_2$ corrections. This also holds true in other formulas for compliance corrections, but we nevertheless carry $\eta_1$along in the formulas for completeness [see (15)].) Now it is also not difficult to see that, if the cracks were oriented instead so that all their normals were pointed horizontally along the x-axis, then we would have one permutation of this matrix and, if instead they were all pointed horizontally along the y-axis, then we would have a third permutation of the matrix. Averaging these three permutations by adding the three $\Delta S$'s together and then dividing by three, we obtain the isotropic compliance corrections matrix  

$$\Delta S_{ij} = \rho\left(\begin{array} {cccccc} 2(\eta_1+\e... ... & 4\eta_2/3 & \cr & & & & & 4\eta_2/3 \cr\end{array}\right).$$ (5)

Then, since the unperturbed compliance matrix is related to Young's modulus E₀, Poisson's ratio $\nu_0$, and shear modulus G₀ by  

$$S_{ij}^{(0)} = \left(\begin{array} {cccccc} 1/E_0 & -\nu_0/E... ...cr & & & & 1/G_0 & \cr & & & & & 1/G_0 \cr\end{array}\right),$$ (6)

we find easily that  

$$\frac{1}{E^*} = \frac{1}{E_0} + \frac{2}{3}\rho(\eta_1 + \eta_2)$$ (7)

and  

$$\frac{-\nu^*}{E^*} = \frac{-\nu_0}{E_0} + \frac{2}{3}\rho\eta_1.$$ (8)

Solving these equations for the $\eta$'s, we have  

$$\eta_1 = - \frac{3}{2\rho}\left(\frac{\nu^*}{E^*} - \frac{\nu_0}{E_0}\right)$$ (9)

and  

$$\eta_2 = \frac{3}{2\rho}\left(\frac{(1+\nu^*)}{E^*} - \frac{(1+\nu_0)}{E_0}\right).$$ (10)

Since $1/G_0 = 2(1+\nu_0)/E_0$ for an isotropic system, we also have the consistency check that  

$$\eta_2 = \frac{3}{4\rho}\left(\frac{1}{G^*} - \frac{1}{G_0}\right).$$ (11)

Similarly, since $1/K_0 = 3(1-2\nu_0)/E_0$, we also have  

$$\frac{1}{K^*} - \frac{1}{K_0} = 2\rho(\eta_2 + 3\eta_1).$$ (12)

This construction shows in part both the power and the simplicity of the Sayers and Kachanov (1991) approach, also used by Bazant and Planas (1998).

The crack influence decomposition parameters $\eta_1$ and $\eta_2$ can be evaluated using the formulas presented for any convenient value of the crack density $\rho$. But this procedure is seen to be most useful if we evaluate the parameters at small $\rho$, since in that limit all the standard methods should give essentially the same results. This approach has been tested and found to be correct.

Once it is known that it makes little difference at low crack densities which theoretical methods we use to estimate the crack-influence parameters, we might as well consider the simplest one, which is surely the non-interaction approximation. Combining (2) and (3) with (11) and (12), we find easily that, within the NIA,  

$$\eta_2 = \frac{8(1-\nu_0)(5-\nu_0)}{15G_0(2-\nu_0)},$$ (13)

and  

$$\eta_1 = -\frac{4\nu_0(1-\nu_0)}{15G_0(2-\nu_0)}.$$ (14)

The ratio of these expressions is  

$$\eta_1/\eta_2 = -\frac{\nu_0}{2(5-\nu_0)}.$$ (15)

This shows that, when $0 \leq \nu_0 \leq 0.5$,  

$$\vert\eta_1/\eta_2\vert \leq 0.05.$$ (16)

So, $\vert\eta_1\vert$ is never larger than about $5\%$ of $\eta_2$,and, for small values of $\nu_0$, the ratio is substantially smaller. A typical value for tight sandstones is $\nu_0 \simeq 0.4$.This corresponds to a wave speed ratio $V_s/V_p = \left[\frac{2(1-\nu)}{(1-2\nu)}\right]^{1/2} \simeq 0.4$,which is typical of both shales and sandstones. It follows that the value of $\vert\eta_1\vert$ is about 4% of that for $\eta_2$ in many important geophysical applications. Thus, its value is suffficiently small so that we are often justified in neglecting $\eta_1$ in data analysis problems for real earth systems and rocks.

Once values of $\eta_1(0)$ and $\eta_2(0)$ are known in this way (using NIA as shown or some other method), we can use the Sayers and Kachanov (1991) method as one convenient way to study and evaluate anisotropic behavior in cracked systems -- thus, providing a simple method of extending the non-interaction approximation results, as quoted by Zimmerman (1991), to nonisotropic systems. We then introduce some crack-crack interactions here in a novel way by making use of the polycrystal of cracked-grains model.

Appendix B summarizes the main analysis using the crack-influence parameter approach of the Sayers and Kachanov (1991) method. We find that, in addition to the two parameters $\eta_1$ and $\eta_2$ that are easily found using effective medium theories such as NIA or DS (differential scheme), there are three more parameters that are expected to play a role in our results at higher crack densities: $\eta_3$, $\eta_4$, and $\eta_5$. These are the only significant quadratic corrections to (4). The analysis shows [see, for example, Eq. (40)] that to second order in the crack density $\rho$, we need to replace $\eta_1$ by $\eta_1 + \eta_4\rho$and $\eta_2$ by $\eta_2 + \eta_5\rho$. In addition, the parameter $\eta_3$ comes into play, but perturbs only the S₃₃ component of the compliance. The Reuss averages of shear modulus and bulk modulus are now given by  

$$\frac{1}{G_R} - \frac{1}{G_0} = \frac{4\rho}{3}\left[\eta_2 + (\eta_5 + 2\eta_3/5)\rho\right]$$ (17)

and  

$$\frac{1}{K_R} - \frac{1}{K_0} = 2\rho\left[\eta_2 + (\eta_3 + \eta_5)\rho + 3(\eta_1 + \eta_4\rho)\right],$$ (18)

respectively, instead of (11) and (12). Furthermore, since it has been found empirically that the shear modulus estimates given by the NIA are quite accurate for the polycrystals of cracked-grains model, we conclude that it is generally true that $\eta_5 + 2\eta_3/5 \simeq 0$, which approximately eliminates one degree of freedom in our three parameter ($\eta_3$, $\eta_4$,$\eta_5$) fitting method.

Also, recall that we are focusing here on penny-shaped (and therefore ellipsoidal) cracks. Phenomenology for other crack shapes may differ somewhat from the discussion presented here (Mavko and Nur, 1978). However, other shapes are beyond our present scope and so will necessarily be treated elsewhere. The numerical experiments considered here all used ellipsoidal cracks.