Microcrack statistics, Weibull distribution and micromechanical modeling of compressive failure in rock

Mechanics of Materials 38 (2006) 664–681 www.elsevier.com/locate/mechmat

Microcrack statistics, Weibull distribution and micromechanical modeling of compressive failure in rock Teng-fong Wong a

a,*

, Robina H.C. Wong b, K.T. Chau b, C.A. Tang

c

Department of Geosciences and Department of Mechanical Engineering, State University of New York, Stony Brook, NY 11794-2100, USA b Department of Civil and Structural Engineering, The Hong Kong Polytechnic University, Hong Kong, China c Centre for Rock Instability and Seismicity Research, Northeastern University, Shenyang, China Received 16 March 2005; received in revised form 28 October 2005

Abstract Statistics of the mechanical and failure properties on the grain scale are often assumed to follow the Weibull distribution in numerical simulations of failure and damage development. To investigate the microstructural basis for such a statistical model of compressive failure in a brittle rock, we consider the development of instability in a wing crack model and establish a methodology whereby the Weibull parameters can be inferred from microstructural data on microcrack density and length statistics for input into finite element simulations. Application of this methodology to six Yuen Long marble samples provides important insights into how different attributes of the microstructure may influence the progressive development of rock failure. The finite element simulations underscore the significant influence of microcrack length statistics, which has not been emphasized in continuum damage mechanics models that usually emphasize the roles of average crack size and crack density. The microstructural data indicate that strength heterogeneity increases with increasing grain size, and this plays a key role in lowering the uniaxial compressive strength, which contributes to the overall decrease of strength with increasing grain size.  2005 Elsevier Ltd. All rights reserved. Keywords: Rock fracture; Grain size; Weibull distribution

1. Introduction Brittle fracture is a failure phenomenon that has been intensively investigated in rock physics and geotechnical engineering. Laboratory data on the fracture strength of brittle rock as a function of lithology, pressure, strain rate, pore fluid, sample size and shape are commonly anaylzed using empir* Corresponding author. Tel.: +1 631 632 8212; fax: +1 631 632 8240. E-mail address: [email protected] (T.-f. Wong).

ical theories (such as the Mohr–Coulomb and Hoek–Brown criteria) that provide much of basis for the interpretation of the stresses associated with faulting and fracturing in engineering, mining, geologic and seismological applications (Jaeger and Cook, 1979; Hoek and Brown, 1980; Lockner, 1995). Although such empirical models have often been cast in physical terms, they actually say very little about the physical mechanisms (Handin, 1969; Paterson and Wong, 2005), and there are still fundamental questions related to the mechanics of brittle fracture that remain unanswered.

0167-6636/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmat.2005.12.002

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

In recent years a broad range of techniques, including optical (e.g. Moore and Lockner, 1995), scanning electron (e.g. Brace et al., 1972) and laser scanning confocal (e.g. Fredrich et al., 1995) microscopy, have been used to image the microstructure and damage, which help to elucidate their critical roles in causing compressive failure in rocks. It is now recognized that damage initiates primarily in the form of extensile microcracks that emanate from local heterogeneities. An applied stress induces the microcracks to propagate and the damage to distribute throughout the sample, and finally the microcracks coalesce and the damage localizes to develop a shear band, leading to strain softening and macroscopic failure (Tapponier and Brace, 1976; Wong, 1982; Mene´ndez et al., 1996). Additional insights have also been gained from mapping out the spatiotemporal distribution of acoustic emission (AE) activity (Lockner et al., 1992). A physical theory of brittle fracture aims to explain both the phenomenological and microstructural observations. The objective of a micromechanical model is to capture the important attributes inferred from microstructural and AE measurements so as to arrive at constitutive relations that describe the macroscopic failure behavior. Since the development of shear faulting and strain localization is a manifestation of crack coalescence, the mechanics of which hinges on the complex interaction of stress fields of the numerous neighboring cracks, micromechanical modeling of the instability requires fairly detailed characterization of the microcrack statistics in conjunction with intensive computation or analytic approximation. Two alternative approaches have been adopted for such an analysis. In the first approach based on continuum damage mechanics (Kachanov, 1986; Krajcinovic, 1989), the constitutive behavior is considered to evolve with ‘‘damage’’, which in a brittle rock refers to one or more internal variables characterizing the density and geometry of the multiplicity of microcracks and pores. Fracture mechanics and analytic approximation of crack interaction are then employed to derive the evolution of the damage parameter with the history of loading. Important advances have been made using the sliding wing crack model (Horii and NematNasser, 1986; Ashby and Sammis, 1990; Kemeny and Cook, 1991). This class of damage mechanics model establishes the micromechanical basis for the development of dilatancy and brittle faulting in terms of the nucleation, propagation, and coales-

665

cence of stress-induced microcracks. It also underscores the fracture mechanics and crack density parameters that control these micromechanical processes. An advantage of the wing crack model is its mathematical convenience. However, it is probably unrealistic to expect the analytic approximations to fully capture the complexity of pore geometry and micromechanical processes in a brittle rock over a relatively broad range of stress conditions. The pre-existng and stress-induced microcrack populations of a brittle rock are characterized by a broad range of sizes (Hadley, 1976; Kranz, 1983; Wong, 1985), and in response to an applied stress field the spatial distribution of microcracking can be highly heterogeneous. To realistically analyze the progressive development of damage in such a solid one may alternatively consider a stochastic model using numerical simulation. In recent years several numerical simulation techniques have been employed. One technique is to use a discrete lattice with cracks represented by interconnected bonds (Lockner and Madden, 1991; Krajcinovic and Vujosevic, 1998; Katsman et al., 2005 ). Another technique based on the discrete element method (Cundall and Strack, 1979; Potyondy et al., 1996) has been particularly effective in modeling compressive failure in porous rocks (Hazzard et al., 2000; Li and Holt, 2002; Boutt and McPherson, 2002). A third technique considers the mechanical response of a stochastic continuum with heterogeneous distributions of mechanical and failure properties that mimic those of a cracked solid. The elastic interaction, damage evolution and stress perturbation are simulated by the finite element (Tang, 1997; Amitrano, 2003; Liu et al., 2004) or finite difference (Fang and Harrison, 2002) method. Typically the statistics of the mechanical and failure properties in such a numerical model are assumed to follow the Weibull distribution (Tang, 1997; Gupta and Bergstro¨m, 1998; Fang and Harrison, 2002; Van Mier et al., 2002; Xu et al., 2004), but there are also models that adopt a Gaussian (Hazzard et al., 2000; Katsman et al., 2005) or uniform (Krajcinovic and Vujosevic, 1998; Amitrano, 2003) distribution. While these stochastic models can simulate many key attributes of the macroscopic failure mode and spatial clustering of AE and damage, it remains unclear to what extent the statistical assumptions on local strength can be justified on rock microstructure. How can quantitative observations on

666

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

microcrack statistics be used to constrain the statistical model of strength in a micromechanical model of compressive failure? Here we focus on the Weibull distribution which is widely used in recent models. Its connection to the fractal characteristics of microcrack size distribution is established within the framework of linear elastic fracture mechanics and wing crack model. Results from the present study provide the theoretical basis for a methodology for inferring the Weibull parameters from rock microstructure. A robust observation on fine- and coarse-grained rocks of a given lithology is that the strength decreases with increasing grain size, approximately following an inverse square root relation (Brace, 1961; Olsson, 1974; Fredrich et al., 1990; Wong et al., 1996). In the context of the wing crack model such a relation is plausible provided that the lengths of the pre-existing microcracks scale with the average grain size (Fredrich et al., 1990; Wong et al., 1996). However, the micromechanical process responsible for such a grain size dependence is expected to be quite complex since it actually involves the interplay of a multiplicity of cracks spanning a broad range of sizes. What are the coupled roles of grain size and microcrack statistics in controlling the strength of a brittle rock? In this study we selected the Yuen Long marble, the brittle strength of which shows appreciable dependence on grain size (Wong et al., 1996). Microcrack statistics were acquired so that the Weibull statistical parameters can be inferred and then input into the RFPA (Rock Failure Process Analysis) finite element code (Tang, 1997) to simulate the progressive development of damage and shear localization. Comparison of the numerical simulations with experimental observations provide useful insights into the micromechanics of compressive failure and its dependence on grain size. 2. Statistical theory of compressive failure 2.1. The Weibull approach Because of the grain-scale heterogeneity, the failure strength in a polycrystalline rock can vary significantly from a local volume to another. To analyze the statistical variation of the bulk failure strength in such a heterogeneous material, Weibull (1951) adopted the statistics of extreme (Gumbel, 1958) to characterize the local failure strength by the probability distribution function

f ðrÞ ¼

 m1   m  m r r exp  r0 r0 r0

and the cumulative probability function   m  r P ðrÞ ¼ 1  exp  r0

ð1aÞ

ð1bÞ

respectively. Here the assumptions are made that the overall failure is primarily controlled by the weaker elements and that the strength of the weakest element is vanishingly small. As a consequence the statistics of failure would involve only two parameters: r0 which is proportional to the mean value of the strength distribution, and a dimensionless parameter m which characterizes the degree of homogeneity in the structure. An infinitely high m value corresponds to a homogeneous structure with a uniform strength, whereas a heterogeneous structure with a broad distribution of local strength is associated with a relatively low m value. 2.2. Flaw size statistics and fracture mechanics In materials science a number of models have been developed to explicitly relate the two Weibull parameters to microstructural and fracture mechanics properties (Freudentahl, 1968; Argon, 1974; Evans and Langdon, 1976; Bazˇant and Planas, 1997), and our analysis here derives from some of these earlier studies. For an ensemble of microcracks of similar geometry, fracture mechanics predicts that a shorter crack would require a higher stress to propagate. Consequently the strength heterogeneity in a brittle rock is intimately related to the microcrack length statistics, which should be characterized with reference to an ‘‘elemental volume’’ V0, with a linear dimension not less than that of the longest crack in the rock mass. Since the lengths of most pre-existing microcracks in rocks are less than the grain size (Hadley, 1976; Kranz, 1983; Wong et al., 1996) one expects the linear dimension of the elemental volume to be comparable to the average grain size. In engineering materials it has often been observed that the length statistics of the longer cracks can be described by an extreme value distribution of the Cauchy type   z q gðaÞ ¼ ; ð2Þ a where a denotes the half-length of a microcrack (McClintock and Argon, 1966; Argon, 1974; Evans

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

and Langdon, 1976). Indeed in geologic materials it has also been suggested that the geometry of cracks and fractures are fractal (Turcotte, 1992; Steacy and Sammis, 1992; Gupta and Bergstro¨m, 1998), which necessarily implies that the size distribution also follows a power law as in (2). The function here is defined so that the cumulative probability for cracks to have lengths longer than a value a1 Rin an elemen1 tal volume is given by the integral V 0 a1 gðaÞ da. Since the distribution (2) is for characterizing the longer cracks, if one indiscriminately applies the power law to cracks of any length, this integral will blow up as a1 becomes infinitely small. Hence the extreme value distribution is associated with a length scale a 0 which defines the lower limit of crack length described by the power law, such that the cumulative probability for cracks longer than this threshold value a 0 in an elemental volume is equal to unity Z 1 V0 gðaÞ da ¼ 1. ð3aÞ a0

To satisfy this equality the lower integration limit (for z > 1) is necessarily given by a0 ¼



ðq Þz V 0 ðz  1Þ

1 z1

ð3bÞ

.

It should also be noted that the mean crack length for such a distribution (for z > 2) is given by Z 1 2ðz  1Þ 0 a. 2 a ¼ 2V 0 agðaÞ da ¼ ð3cÞ ðz  2Þ a0 The crack length statistics in the elemental volume V0 results in a distribution of mechanical strength, with a mean strength denoted by r0. An important feature of the Weibull model is that the strength of a large volume V containing many elemental volumes is size dependent. To scale up in volume it is assumed that the failure events in the V/V0 elemental volumes are non-interactive so that the failure in each elemental volume represents an independent event. Accordingly if the volume V is homogeneously stressed and P denotes the probability for failure to arise from cracks with lengths longer than a1 then  VV Z 1 0 gðaÞ da . ð4aÞ 1P ¼ 1V0 a1

Noting that lim ð1  x=jÞj ¼ ex the cumulative probj!1

ability for failure to arise from cracks with lengths

667

longer than a1 as V/V0 ! 1 is (McClintock and Argon, 1966; Evans and Langdon, 1976):  Z 1  P ¼ 1  exp V gðaÞ da "

a1 z

ðq Þ V ¼ 1  exp  ðz  1Þ ða1 Þz1 "  z1 # V a0 ¼ 1  exp  . V 0 a1

#

ð4bÞ

Concepts of linear elastic fracture mechanics can then be employed to relate the crack size distribution to the probability of fracture. The stress intensity factor of a mode-I p crack ffiffiffiffiffiffi of length 2a subjected to a stress r is K I ¼ yr pa, where y is a geometric constant which is equal to 1 for a two-dimensional crack embedded in an infinite medium (Lawn, 1993). If KI attains the value of the critical stress intensity factor KpIcffiffiffiffiffithen failure develops at the ffi stress r ¼ K =ðy pa Þ which also implies that Ic pffiffiffi pffiffiffi a ¼ K Ic =ðry pÞ. If the volume V is subjected to a uniform stress, then according to (4b) the cumulative probability that failure develops at stresses less than r is given by "  pffiffiffiffiffiffiffi2z2 # V yr pa0 P ðrÞ ¼ 1  exp  . ð4cÞ V0 K Ic It can be seen that the probability of failure increases with increasing sample volume and applied stress. The Weibull parameters can be related to the statistical distribution of crack lengths by comparing Eqs. (1b) and (4b); the parameter m is related to the exponent z in the power law (2) by m ¼ 2z  2.

ð5aÞ

Since the crack length distribution falls off rapidly for large values of z, the strength distribution is relatively narrow for large values of m and accordingly the strength distribution in such a structure is expected to be relatively homogeneous. The other Weibull parameter r0 depends on the fracture toughness and volume V according to  1=m V K Ic pffiffiffiffiffiffiffi . r0 ¼ ð5bÞ V0 y pa0 Microstructural observations have revealed that the compressive failure of brittle rock initiates from predominately extensile cracking on the grain scale (Wong, 1982; Myer et al., 1992; Moore and

668

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

Lockner, 1995). It is therefore justified pffiffiffiffiffiffi to adopt the tensile fracture criterion K Ic ¼ yr pa for failure initiation in an elemental volume and indeed a similar criterion was used by Lockner and Madden (1991) in their network model. However, microstructural and AE (Lockner et al., 1992) observations also show that compressive failure in rock generally involves the progressive development of microcracks that first propagate stably and distribute throughout the volume before they coalesce to develop a shear band, leading to strain softening and macroscopic failure (Tapponier and Brace, 1976; Wong, 1982; Moore and Lockner, 1995). Hence the parameter y should be chosen so that the fracture criterion in an elemental volume can reflect the initiation, propagation and coalescence of microcracks that ultimately lead to damage localization. 2.3. The wing crack model: compressive strength and its dependence on crack density Following Gupta and Bergstro¨m’s (1998) suggestion we adopted the wing crack model (Horii and Nemat-Nasser, 1986; Kemeny and Cook, 1991; Ashby and Sammis, 1990) to relate the compressive failure stress to the critical stress intensity factor KIc. This damage mechanics model considers sources of tensile stress concentration that are located at the tips of inclined pre-existing cracks (with length 2a). The applied stresses induce a shear traction on the crack plane and, if this resolved shear traction is sufficiently high to overcome the frictional resistance along the closed crack, frictional slip results in tensile stress concentrations at the two tips of the sliding crack, which, in turn, may induce ‘‘wing cracks’’ to nucleate there. The driving force for nucleation is characterized by the stress intensity factor KI at the site of wing crack initiation, which is a function of the applied stresses, friction coefficient l of the sliding crack, and angle w between the sliding crack and maximum principal stress r1 (Fig. 1a). With increased loading, KI will attain the critical value KIc, at which point a wing crack nucleates and grows out of the initial plane of the sliding crack along orientations sub-parallel to r1. Also, as the applied stress is increased, shorter and less favorably oriented sliding cracks will be activated, thus nucleating an increasing number of wing cracks distributed throughout the sample. With the progressive increase in both number and dimensions of wing cracks, their stress fields will interact with one another in a complex manner

Fig. 1a. Schematic diagram of the wing crack model, showing the orientations of the sliding crack and wing cracks with respect to the principal stresses.

and ultimately leads to crack coalescence and shear localization. Analytic approximations of these complex interactions can be derived with reference to certain idealized crack geometry, and here we will adopt Ashby and Hallam’s (1986) model for a rectilinear array of wing cracks emanating from sliding cracks. For this configuration under uniaxial compression Wong et al. (1996) have derived an analytic expression p for ffiffiffiffiffiffi the peak strength in the form of rc ¼ K Ic =ðy paÞ with pffiffiffiffiffiffiffiffiffiffiffiffiffi " # 1 þ l2  l 1 y¼ C 3 ‘cr þ pffiffiffi 3=2 1=2 3ð1 þ ‘cr Þ ð1 þ ‘cr Þ þ

½2e0 ð‘cr þ cos wÞ1=2 pffiffiffi ; p

ð6aÞ

where the critical normalized length ‘cr denotes a positive root of the implicit equation " # C1 C2 C4 þ C3  1=2 3=2 3=2 2ð‘ þ cos wÞ ð1 þ ‘Þ 2ð1 þ ‘Þ " # 3C 2 C4  C3‘ þ ¼0 ð6bÞ 5=2 2ð1 þ ‘Þ ð1 þ ‘Þ1=2 pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi with C 1 ¼ 2e0 =p, C 2 ¼ 1 þ l2  l, C3 = 0.23 and C4 = 0.577. The damage parameter e0 is related to the number of cracks per unit area NA (pre-existing in an undeformed sample) by e0 = NAa2, corresponding to the two-dimensional equivalent of the ‘‘crack density’’ introduced by Budiansky and O’Connell (1976).

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

In using expression (6a) to evaluate the parameter y that enters into expression (5b) for the Weibull parameter r0, our approach is different from the conventional applications of the Weibull theory and wing crack model in several important aspects. First, the value of y according to (6a) is typically less than unity for realistic values of e0, l and w, whereas y > 1 for standard geometric configurations encountered in materials science, such as a surface crack or an internal crack loaded in tension. Second, the value of y depends on the initial crack density and may increase significantly with increasing e0. For a relatively high crack density the interaction among the stress fields of the many cracks is sufficiently strong that the wing cracks can readily coalescence at a relatively low stress corresponding to a relatively high y value. Third, in the damage mechanics model the sliding cracks are implicitly assumed to have a uniform length of 2a and as soon pffiffiffiffiffiffi as the stress r1 attains the peak value of K Ic =ðy paÞ, an instability would develop catastrophically throughout the whole volume. In contrast, in a stochastic model the pre-existing cracks and consequently the local compressive strengths are described by statistical distributions, and hence failure would develop only locally in a subset volume of the sample when the conditions (5b) and (6a) are satisfied. The local stress fields would subsequently be equilibrated and accordingly failure develops in a progressive manner. Lastly, although the failure arises from the interaction and coalescence of numerous wing cracks that have nucleated from sliding cracks with a range of w values, it is implicitly assumed in most applications of the damage mechanics model that the cooperative effects of the multiplicity of cracks can be approximated by choosing an ‘‘effective’’ value of w. For example Ashby and Sammis (1990) used a value of 45, and in the equations above we have followed Wong et al. (1996) who used w = (1/2) tan1(1/l), the optimal angle for wing crack nucleation. Given the idealizations intrinsic to the wing crack model, the parameter y cannot be uniquely determined; nevertheless, the use of Eq. (6a) for fixed values of l and w provides a relative measure of the dominant control of e0 over the stress intensity factor and critical stress at the onset of instability in an elemental volume.

669

Fig. 1b. Experimental data on uniaxial compressive strength of Yuen Long marble as a function of average grain size. (After Wong et al., 1996.) The filled circles correspond to the six samples selected in this study for more detailed microstructural characterization.

To test this approach we characterized quantitatively the microcrack statistics of marble cores from Yuen Long in the northwestern part of the New Territories of Hong Kong, China. Three varieties of the Yuen Long marble can be identified: the white marble consists primarily of crystalline calcite, while the dark and grey marbles have some silt and clay (Yuen, 1990). The marble is relatively compact, with effective porosity less than 0.5%. According to the petrophysical measurements of Wong et al. (1996), average grain sizes of Yuen Long marble vary by one order of magnitude. The dark marble has relatively fine grains ranging in size from 47 to 86 lm. In contrast the white marble has the coarsest grains with sizes ranging from 368 to 465 lm, while grain sizes in the grey marble fall in the range 163– 217 lm. They also determined the uniaxial compressive strengths of 13 samples, which approximately scale with the inverse square root of the grain size (Fig. 1b). The mechanical tests were conducted on cylindrical samples with diameters of about 60 mm and length-to-diameter ratio ranging from 2 to 2.2. We selected six of these samples for more detailed microstructural investigation. Some of the petrophysical and fracture mechanics data are summarized in Table 1.

3. Quantitative characterization of the microstructure of Yuen Long marble

3.1. Microcrack statistics

Eqs. (5) and (6) prescribe the microstructural data to be acquired for input into a Weibull model.

In Yuen Long marble most of the pre-existing cracks are intragranular, but there are also some

670

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

Table 1 Mechanical properties and microcrack statistics of Yuen Long marble Marble type

White marble

Grey marble

Dark marble

Mean properties Density, kg m3 Young’s modulus, GPa Fracture toughness KIc, MPa m1/2 Tensile strength, MPa Coefficient of friction

2727 58.57 0.878 8.65 0.6–0.8

2705 72.94 1.567 7.52 0.6–0.8

2720 91.12 1.836 4.89 0.9–1.0

Samples

S19

S21

S6

S17

S16

S24

Mechanical properties Young’s modulus, GPa Poisson’s ratio Uniaxial compressive strength (UCS), MPa

34.65 0.25 55.85

82.49 0.35 56.22

107.18 0.29 87.61

37.68 0.22 62.79

83.13 0.31 143.67

99.11 0.27 166.09

Microstructural data Mean grain size, mm Mean crack length (2a), mm Crack density e0 Cauchy distribution exponent z (q*)V0 Correlation coefficient squared (r2)

0.445 0.141 0.243 2.017 0.007 0.69

transgranular cracks. Wong et al. (1996) observed that the lengths of microcracks in the fine-grained dark marble are comparable to the grain size, while cracks in the other two marbles with coarse grains have lengths only about 20% of the grain size. While it may be desirable to characterize the microstructural attributes using other microscopy techniques with finer resolutions, we consider it adequate to use optical microscopy for the present study since the Weibull model is based on the statistics of extreme, which in the present context corresponds to the subset of longer microcracks readily resolvable under an optical microscope at magnifications up to 400·. Standard thin sections (of thickness 30 lm) covering an area of 45 · 25 mm2 were prepared. On each thin section 10 areas (each equal to 3 · 3 mm2) were randomly selected for measuring the distribution of crack length. Typically each elemental area covers a minimum of six grains, and therefore for each thin section the crack distribution was based on a sample of more than 60 grains. The crack length data were typically grouped into bins corresponding to a length interval Da = 25 lm, except for the longest cracks for which larger length intervals might also be used for the more sparse population. To relate the two-dimensional data from thin section to a distribution defined with respect to the elemental volume, we will assume that the cavities can be approximated as ‘‘tunnel cracks’’ with uniform width t (in the third

0.368 0.120 0.056 2.195 0.012 0.88

0.187 0.094 0.0408 2.436 0.005 0.94

0.217 0.132 0.1352 2.452 0.008 0.78

0.051 0.079 0.035 2.646 0.002 0.89

0.047 0.084 0.022 2.548 0.003 0.91

dimension corresponding to the thickness of the thin section). We assume that the scanned area covers a sufficiently large number of grains to be representative of the crack statistics in the elemental volume, so that data acquired on the two-dimensional section can be related to the crack statistics in the volume by geometric probability concepts (Budiansky and O’Connell, 1976; Wong, 1985). If a total of N microcracks were observed over the elemental area A0 = V0/t and if n of these have half-length in the range (a  Da/2, a + Da/2), then we can calculate the relative frequency according to V 0 gðaÞ ¼ A0 tgðaÞ ¼

n . N Da

ð7aÞ

Our data so acquired for the relative frequency V0g(a) as a function of a (chosen to be the mean value in each bin) are shown in Fig. 2. By regression the microcrack data were fitted to a power law to determine the parameters z and (q*)zV0. These two parameters and the square of the correlation coefficient (r2) for the fit to a power law are compiled in Table 1. It should be noted that since we have used a unit of mm1 for V0g(a), the inferred values of (q*)zV0 are in unit of mmz2. The Weibull parameter m and the characteristic length a 0 can then be inferred by direct substitution into Eqs. (5a) and (3b), respectively (Table 2). Several features of the crack length distributions shown in Fig. 2 should be noted. First, the overall

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

671

Fig. 2. Microstructural data on relative frequency of crack length in Yuen Long marble samples (a) S19, (b) S21, (c) S6, (d) S17, (e) S16 and (f) S24. The data are fitted to a power law shown as the solid curve. The two regression parameters and correlation coefficients are compiled in Table 1.

420.5 94.18 (98.92) 354.7 79.66 (86.58) The bracketed values were from simulations with steel platens bonded to both sample ends. a

332.5 77.19 (78.72) Input parameter and simulated results for KIc = 0.90 MPa m1/2 r0, MPa 289.7 288.6 47.45 (48.98) 57.37 (59.14) Simulated UCSa, MPa

181.5 40.05 (45.01)

724.2 170.66 (181.87) 610.8 156.73 (158.07) 312.5 70.52 (70.67) 572.6 136.04 (138.6) Input parameter and simulated results for KIc = 1.55 MPa m1/2 r0, MPa 498.9 497.0 83.69 (89.17) 96.07 (102.92) Simulated UCSa, MPa

S24 S16

98 · 196 = 19 208 0.051 5 · 10 0.25 91.12 3.292 0.0188 0.330 92 · 184 = 16 928 0.217 20 · 40 0.25 72.94 2.904 0.0288 0.522

S17 S6 S21

81 · 162 = 13 122 0.368 30 · 60 0.25 58.57 2.39 0.0208 0.386

S19

79 · 158 = 12 482 0.445 35 · 70 0.25 58.57 2.034 0.0075 0.641 Number of square elements Linear dimension of element, mm Sample size, mm Poisson’s ratio Mean Young’s modulus (E0), GPa Weibull parameter m Characteristic length a 0 , mm y(e)

Dark marble Grey marble White marble

Table 2 Input parameters and predictions of the finite element model

106 · 212 = 22 472 0.047 5 · 10 0.25 91.12 3.096 0.0181 0.284

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

80 · 160 = 12 800 0.187 15 · 30 0.25 72.94 2.872 0.0193 0.347

672

trend of the data follows a power law, and especially for three samples (S6, S16 and S24) the fit were very good over the full range of crack lengths measured. Nevertheless, measurements for the shorter cracks (with half-lengths less than 50 lm) in the other three samples (S19, S21 and S17) are appreciably lower than predicted by a power law. Although this may indeed be an intrinsic feature of crack statistics in these samples, it is also possible that this arises from the limitation of the optical microscope to resolve the finer cracks. Second, the Weibull parameter m inferred from microcrack statistics falls in the range of 2.03 to 3.29 (Table 2 and Fig. 3). Such low m values imply that the three marbles are associated with highly heterogeneous distributions of strength on the grain scale. Third, there is an overall trend for m to increase with decreasing grain size, implying that microcrack lengths have narrower distributions and strength distributions are more homogeneous in the finegrained marbles. While such a trend is particularly pronounced when we compare the white Yuen Long marble with the grey marble, the contrast in m value between the grey and dark marbles are not as significant. The significant difference arises possibly because the dark Yuen Long marble has more transgranular cracks (Wong et al., 1996): its mean crack length is almost twice the average grain size (Table 1), and therefore the microcrack length scales with the grain size in a manner more complex that in the two coarse-grained marbles with predominately grain-boundary and intragranular cracks much shorter than the average grain size.

Fig. 3. The Weibull parameter m and characteristic length a 0 inferred from microcrack statistics as functions of average grain size in six Yuen Long marbles.

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

673

Fourth, the characteristic length a 0 is quite insensitive to the grain size (Table 2 and Fig. 3). This somewhat surprising result can be explained by the trade-off between crack size and its heterogeneity. The inferred values of the parameter z are larger than 2 and hence Eq. (3c) can be used to relate a 0 to the mean crack size  a and homogeneity index z. With increasing grain size the crack lengths and therefore  a increase, while the crack size distribution becomes broader and therefore z and (z  2)/(z  1) both decrease. Because of two opposite effects the value of a 0 is not sensitive to the change in grain size. 3.2. Crack density To evaluate the second Weibull parameter we need to first characterize the crack density to calculate the parameter y for the wing crack model. Since Eq. (6a) was derived with the implicit assumption of uniform crack size, there is not a unique method to determine the parameter e0 for an ensemble of cracks of non-uniform size and geometry (Budiansky and O’Connell, 1976). Following Hadley’s (1976) suggestion, Wong et al. (1996) adopted the following formula to evaluate the average crack density for N cracks observed in a two-dimensional section of area A0 e0 ¼

3N c2 ; 4pA0

ð7bÞ

where c2 denotes the mean value of the crack length (c) squared. For a simplest case of a two-dimensional tunnel crack we have c = 2a. We compile in Table 1 the values of e0 for our six Yuen Long marble samples determined by Wong et al. (1996). It should be noted that they also presented corresponding crack density values for subsets of ‘‘sub-vertical cracks’’ in these samples which were not considered in this study. The higher values adopted here are considered to be more consistent with our length statistics which were acquired from microcracks in not only the subvertical but all orientations as observed in the thin sections. The initial crack density values seem to be somewhat lower in the fine-grained dark marble samples, but as noted by Wong et al. (1996) there is not a clear trend for e0 to increase with grain size. The value of y as a function of e0 according to Eq. (6) is plotted in Fig. 4 for three values of the friction coefficient l ranging from 0.4 to 0.8. It can be seen that y is not very sensitive to such a variation in l. As discussed in the next section we used

Fig. 4. The parameter y as a function of the crack density e0 according to Ashby and Hallam’s (1986) wing crack model, for friction coefficient l=0.4, 0.6 and 0.8, respectively. (After Wong et al., 1996.)

in our numerical simulations the intermediate values of y corresponding to l = 0.6 (Table 1). The y values are all less than unity, ranging from 0.28 to 0.64 (Table 2). In our two-dimensional numerical simulations we take the cross-section to be a square, with width equal to the average grain size D as observed in the thin section. The volume V = D2t of this element in our finite-element grid is considered to be comparable to the elemental volume V0 = A0t, and according to (5b) the mean strength of an element pffiffiffiffiffiffiffi in our numerical model is simply r0 ¼ K Ic =ð pa0 yÞ. This second Weibull parameter were evaluated for KIc = 0.90 MPa m1/2 and 1.55 MPa m1/2, which correspond approximately to the lower end of fracture toughness values of these Yuen Long marble samples (Table 1, unpublished data of R.H. Wong) determined from Chervon Bend (V-notch) specimen in accordance with the standard testing method of the International Society for Rock Mechanics Commission on Testing Methods (1988). As shown in Table 2 the values of y seem to be somewhat lower and r0 seem to be somewhat higher in the finegrained dark marble samples, but there is not a clear trend for either parameters to vary systematically with grain size. 4. Numerical simulation of progressive failure under compressive loading Our microstructural data for the Yuen Long marble samples provide constraints on the Weibull parameters pertinent to the micromechanical

674

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

analysis of compressive failure and its dependence on grain size. In recent years finite element, finite difference and network simulations have been employed for such a micromechanical analysis. In this study the two-dimensional finite element code RFPA (Tang, 1997) was selected to simulate the failure behavior. 4.1. The RFPA code and choice of constitutive and statistical parameters In our finite element model the spatial distribution of failure stress is given by the two-parameter Weibull distribution (1). A state of plane strain is assumed and an element may fail in a shear or tensile mode. Specifically the Mohr–Coulomb failure criterion with a tension cut-off was adopted. With the convention that compression is positive, if we denote the maximum and minimum principal stresses by r1 and r3, respectively, then the criterion for shear failure to develop in an element is r1 

ð1 þ sin /Þ r3 P rc ; ð1  sin /Þ

ð8aÞ

where rc and / denote the uniaxial compressive strength and internal friction angle, respectively. In our finite element model the angle / is fixed while the uniaxial compressive strength rc in an element is assigned a value following the Weibull distribution (1), with the homogeneity index m and mean value r0 constrained by our microstructural data according to Eqs. (5a) and (5b), respectively. Failure may alternatively develop by a tensile mode if the minimum principal stress satisfies the criterion r3 6 rt .

ð8bÞ

This tension cut-off is assumed to be proportional to the uniaxial compressive strength, with the ratio g = rt/rc on the order of 0.1. Since the ratio g is fixed in our model, the tensile strength also follows a Weibull distribution (1) characterized by the two parameters m and gr0. The elements are assumed to be isotropic and elastic. While the Poisson’s ratio is assumed to be identical in all the elements throughout the deformation, the Young’s modulus is spatially distributed according to a Weibull parameter m identical to that for the shear and tensile failure stresses. In response to an externally applied stress, the elastically heterogeneous material develops a spatially heterogeneous stress field and whenever the stress

field in an element satisfies the failure criterion (8a) or (8b), it instantaneously drops to a residual strength level. A degradation factor k is introduced to characterize the ratio between the residual and peak strength in an element, and to simulate such a stress drop the initial Young’s modulus E is assumed to drop concomitantly to the value E 0 = kE in the elastic damaged material. Such an instantaneous degradation of stiffness is accommodated by corresponding redistribution of the local stresses. As the sample continues to deform, more elements will be damaged and a scalar parameter D is introduced to characterize the damage evolution. Specifically it is defined in terms of the relative change of the stiffness, such that the Young’s modulus in a damaged element is given by E0 ¼ ð1  DÞE.

ð9aÞ

If the residual stress level remains constant in a damaged element, then its Young’s modulus E 0 can be considered to be inversely proportional to the accumulated strain, so that E 0 /E will approach 0 and the damage parameter D approach unity as the strain becomes infinitely large. Furthermore the damage parameter should satisfy the initial condition 1  D = k immediately after the onset of damage, when the maximum principal strain e1 = rc/E for shear failure and the minimum principal strain e3 = rt/E for tensile failure. Accordingly the scalar damage parameters are given by D¼1

krc Ee1

ð9bÞ

for shear failure, and D¼1þ

krt Ee3

ð9cÞ

for tensile failure. (Note that e3 < 0 for tensile strain.) 4.2. Uniaxial compressive failure of Yuen Long marble Before we present our simulation results, it is important to first discuss the strengths and limitations of the code and our approach. A useful feature of the RFPA code is that it can simulate the spatial evolution of damage that mimics the AE and microstructural observations, thus providing useful insights into the progressive development of damage and shear localization through the post-peak stage

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

(Tang, 1997; Tang et al., 2000). A similar numerical approach was adopted by Amitrano et al. (1999) and Amitrano (2003) to simulate the effect of pressure on the brittle-ductile transition. To incorporate mineralogical heterogeneity Liu et al. (2004) developed a ‘‘rock and tool interaction code’’ on the basis of RFPA. Nevertheless, such a numerical approach has the limitation that for any micromechanical models to be mathematically tractable, it is necessary to make a number of idealized assumptions. In this study we have clarified the microstructural basis for the Weibull statistical parameters for strength, but our analysis cannot be immediately generalized to the statistics of the elastic moduli. Since the elastic moduli in a cracked solid depend sensitively on the crack density (Walsh, 1965; Budiansky and O’Connell, 1976) it is plausible that they can be described by Weibull statistics as for strength, but it is unclear whether the parameter m should be identical for statistical distributions of Young’s modulus and strength. For a similar reason it is unlikely that the Poisson’s ratio can be assumed to be homogenous and constant during deformation. Nevertheless the variation is probably not very large, and here in our simulations a constant value of t = 0.25 was assumed as in previous RFPA studies, and the mean values of Young’s modulus E0 were assumed to be identical to the experimental values (Tables 1 and 2). Three other parameters are kept constant in our model. In previous RFPA simulations a value of / = 30 was commonly used. Amitrano et al.

675

(1999) observed that the choice of internal friction angle has strong influence on the development of failure mode and brittle-ductile transition. In particular shear localization would develop in a more diffuse manner with increasing /. Since we here focus on brittle failure, the angle / was fixed at 30. In previous RFPA simulations the ratio g = rt/rc was commonly assigned a value of 0.1, which is reasonable in light of laboratory data on the enhancement of uniaxial strength in compression (Jaeger and Cook, 1979). The degradation factor k presumably controls the macroscopic stress drop during strain softening, but we are not aware of any direct microstructural or rock mechanics justification for specific choice of this factor. While these issues clearly warrant more systematic investigations in the future, we follow earlier studies to assign a value of 0.1 to both parameters g and k. The numerical tests were conducted under constant displacement rate, and the mesh sizes are listed in Table 2. As discussed before, the linear dimension of each element is set to be comparable to the grain size. Fig. 5a shows the simulated stress–strain curves for marble S19, S21 and S6 assuming KIc = 0.90 MPa m1/2, and Fig. 5b shows the curves for marble S17, S16 and S24 assuming KIc = 0.15 MPa m1/2. The simulated uniaxial compressive strengths for all six samples are compiled in Table 2 and plotted in Fig. 6. It can be seen that the experimental data are bracketed by the simulated values for KIc = 0.90 MPa m1/2 and 1.55 MPa m1/2. The laboratory measurements indicate that the fine-grained marbles have higher toughness (Table 1), and indeed the

Fig. 5. Simulated axial stress–strain curves for marble samples (a) S19, S21 and S6, and (b) S17, S16 and S24. The assumed values of fracture toughness are as indicated.

676

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

The most conservative estimate of this effect is predicted by the convention model in which failure in an elemental volume is viewed as an independent event, in which case the expected strength can be calculated as the first moment of the Weibull distribution (McClintock and Argon, 1966; Argon, 1974), and the compressive strength is expected to depend on the Weibull parameters according to   rpeak 1 1 ¼C 1þ . ð10aÞ m r0

Fig. 6. Comparison of experimental data and simulated values of the uniaxial compressive strength for two different values of fracture toughness.

laboratory measurements of uniaxial compressive strength of the two dark marbles seem to be in better agreement with simulated values for the higher value of 1.55 MPa m1/2, whereas the laboratory measurements of the white and grey marbles seem to be in better agreement with simulated values for 0.90 MPa m1/2. To assess the sensitivity of the model to the parameters /, g and k, we conducted a number of simulations using the Weibull parameters for sample S19 assuming KIc = 0.90 MPa m1/2. Unlike the simulations in Table 2 (for which all three parameters were kept constant at / = 30 and g = k = 0.1) we allowed the value of one of these three parameters to vary while keeping the others constant. Our simulations show that the RFPA model results for uniaxial compressive strength are not very sensitive to changes in either / or k. For example, the strengths calculated for k = 0.067 and 0.20 were 50.39 MPa and 52.20 MPa, respectively and the strengths for / = 25 and 35 were 47.80 MPa and 48.59 MPa, respectively. These are within ±10% of the strength of 47.45 MPa for S19 in Table 2. However, our simulations also show that the model is very sensitive to the choice of g. The peak stress was observed to increase with increasing g, with calculated values of 35.76 MPa and 66.36 MPa for g = 0.067 and 0.2, respectively. Our simulations indicate that the uniaxial compressive strength scales as the Weibull parameter r0 and can increase significantly with increasing m.

Since pffiffiffi the gamma function increases from a value of p=2ð¼ 0:886Þ to 1 as m increases from 2 to 1, the variation in strength is expected to be relatively small (Fig. 7). A somewhat stronger dependence on the Weibull parameter m is predicted by the effective medium model of Xu et al. (2004), who attempted to account for damage evolution and its effect on the stress– strain curve to arrive at rpeak 1 1 ¼ ðemÞ m . r0

ð10bÞ

In a highly heterogeneous material the strength ratio is predicted to drop to <0.4 (Fig. 7). In our finite element model the interaction among elemental volumes and damage localization are accounted for in a more explicit manner than the effective medium model. Accordingly the RFPA model shows the strongest dependence on material heterogeneity. On the basis of their results for m ranging from 1.1 to 200 (with / = 30) Liu et al. (2004) obtained an empirical fit to their RFPA simulations

Fig. 7. Normalized peak stress under uniaxial compression according to predictions of the non-interactive (Eq. (10a)), effective medium (Eq. (10b)) and RFPA (Eq. (10c)) model. Simulations from this study are shown by the dark symbols.

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

  rpeak m 1 . ¼ 0:85928  0:80668 exp  10:68877 r0 ð10cÞ The strength ratio may drop precipitously to as low as 0.1 in a highly heterogeneous material. Our simulated results for the six marble samples basically agree with this empirical relation (Fig. 7). While our simulations do show an overall trend for compressive strength to decrease with increasing grain size, this phenomenon did not arise simply from the scaling of crack length with grain size (e.g. Fredrich et al., 1990) or the variation in crack density (e.g. Wong et al., 1996). In addition the statistics of microcrack length as characterized by the Weibull parameter m have dominant influence, as expressed by Eq. (10c). A coarse-grained marble has a more heterogeneous distribution of crack length and a lower m value (Fig. 3), and consequently its compressive strength may be appreciably lower than that of a fine-grained marble. The other Weibull parameter r0 is related to the characteristic length a 0 , geometric factor y and frac-

677

ture toughness KIc. Its value reflects the complex interplay of grain size, crack length statistics and crack density through Eq. (5b). The parameter y is very sensitive to the crack density. As discussed above, with increasing grain size the crack lengths increase while the crack size distribution becomes broader, bring in two opposite effects which cause the value of a 0 to be somewhat insensitive to grain size. 4.3. Damage evolution and shear localization To mimic the laboratory configuration we also conducted simulations with steel spacers bonded to both ends of the rock sample. Since the steel is assigned a Young’s modulus of 300 GPa which is at least three times that of a marble sample, the elastic mismatch induces a slight enhancement of confinement near the sample ends which results in compressive failure at a peak stress slightly higher than the unrestrained case (Table 2). The damage evolution for the grey marble S6 for such a loading configuration is illustrated in Fig. 8. Damage was

Fig. 8. Spatial distribution of damage according to RFPA simulation of the grey marble sample S6. The seven snapshots were from the pre-peak stage (0.6 and 0.8 of the peak stress), the peak stress, and post-peak stage (with stress drops corresponding to 0.2, 0.4, 0.6 and 0.8 of the peak stress). White circles represent earlier damage events, and size of the circles provides relative magnitude of strain energy release of a damage event.

678

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

not appreciable until the stress reached 0.6 of the peak stress. While the damage continued to accumulate as the stress was increased to 0.8 and finally the peak value, it was diffusely distributed in the pre-peak stage, and did not localize to develop a tabular cluster until well into the post-peak stage. Fig. 9 shows the localization patterns in three samples (S19, S6 and S16) at two different post-peak stages. The grey scale measures the maximum shear stress (equal to half the principal stress difference) at an element. While the relatively light zones represent ligaments subjected to enhanced shear stresses, the very dark zones can be considered to be proxies for shear bands with relatively low shear resistance. The top panel shows samples that had undergone stress drops equal to 0.6 of the peak stresses. While the curvature of shear band in the white marble (with more coarse grains and stronger mechanical heterogeneity) seems to change progressively from near the top to the bottom of the sample, the band

angles in the grey and dark marble samples are relatively constant. The lower panel shows samples that had undergone stress drops equal to 0.8 of the peak stresses. There seems to be an overall trend for the bands to be wider in the fine grained samples. The loading machine used by Wong et al. (1996) was not sufficiently stiff for stabilizing the post-peak behavior, and therefore they were not able to retrieve Yuen Long marble samples at different stages of strain softening for microstructural characterization and comparison with the numerical simulation. Nevertheless such characterizations have been conducted on Westerly granite, Indiana limestone, Berea sandstone and Darley Dale sandstone (Wong, 1982; Myer et al., 1992; Moore and Lockner, 1995; Mene´ndez et al., 1996; Wu et al., 2000), which show damage localization qualitatively similar to the numerical simulations. While the spatial clustering patterns reveal differences in post-peak failure modes among the three

Fig. 9. The spatial distribution of maximum shear stress (half of the principal stress difference) in three marble samples at two post-peak stages (with stress drops corresponding to 0.6 and 0.8 of the peak stress). Dark zones represent shear bands with relatively low shear resistance.

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

Fig. 10. Evolution of the areal ratio of damage with axial strain in three Yuen Long samples from RFPA simulations.

marbles, they provide little insights into the damage evolution before the peak stress was attained. To investigate this issue we consider the spatial extent of damage as a function of strain. If we denote the total number of elements in the grid by NFE and if at given stage of deformation we denote the cumulative number of damage events by NDE then the ratio NDE/NFE represents an areal ratio of damage in the sample. This ratio is plotted as a function of the axial strain for the samples S19, S6 and S16 in Fig. 10. The peak stresses are indicated by the arrows, and at these points there is an overall trend for the damage ratio NDE/NFE to increase with increasing grain size and decreasing m value. This implies that in preparation for the ultimate failure, relatively more damage is distributed in a sample with a low m value and highly heterogeneous strength distribution. On the other hand the upsurge in damage during the post-peak stage seems to be more pronounced in a sample with finer grains and higher m value, with the implication that the failure is more brittle. 5. Discussion The first objective of this study was to establish the microstructural basis for Weibull statistical parameters associated with the development of instability in a wing crack model. A methodology was also formulated whereby microcrack density and length statistics can be characterized quantita-

679

tively from microscopy observations for input into numerical simulations of the progressive development of compressive failure and damage localization. Application of this methodology to six Yuen Long marble samples provides important insights into how different attributes of the microstructure may influence rock failure. In particular our finite element simulations underscore the significant influence of microcrack length statistics, which has not been emphasized in continuum damage mechanics models that usually emphasize the roles of average crack size and crack density. In this study we have concentrated on the failure of marble under uniaxial compression. In principle our analysis can readily be extended to a triaxial stress field (Batdorf, 1977; Evans, 1978), although it is likely that the geometric parameter y and degradation factor k need to be calibrated accordingly. Certain aspects of the microcrack statistics in other rock types such as granite and sandstone have been investigated in previous studies, and indeed it will be of interest to explore to what extent their influences on compressive failure behavior are similar to our observations in the Yuen Long marble. Our analysis and microstructural observations require the Weibull parameter m to be greater than 2, and while our data fall in the typical range of m = 2–6 reported for engineering materials (McClintock and Argon, 1966), a value as low as 1 had been adopted in some previous simulations of rock failure without adequate microstructural or rock mechanics justification. These are important issues in rock physics that warrant more systematic studies in the future. The RFPA code was selected in the present study. While it provides a computationally efficient tool to simulate the micromechanics of compressive failure, our results also underscore several aspects of the finite element code that can be improved. The microstructural basis for the statistical distribution of elastic moduli has not been firmly established. More parametric studies should be conducted to gain a deep understanding of the roles of the degradation parameter k and the ratio g. Lastly, it will represent an important advance if a three-dimensional version of such a code can be developed in the future. Acknowledgments Preliminary microstructural measurements and finite element simulations were conducted by

680

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681

Ming-Ruo Jiao and the wing crack model parameters were calculated by Wei Zhu. The first author would like to thank Joanne Fredrich, who first pointed out to him the possible links between Weibull statistics and wing crack model some years ago during an earlier study of the grain size effect. We have benefited from comments of an anonymous reviewer. The research at Stony Brook was partially supported by the Office of Basic Energy Sciences, US Department of Energy under grant DE-FG02-99ER14996. The research at Hong Kong Polytechnic University was funded by RGC Internal Competitive Research Grant (DA) A-PD49 and A-PE31. The research at Shenyang was funded by the NSFC Grant #50374020. References Amitrano, D., 2003. Brittle-ductile transition and associated seismicity, Experimental and numerical studies and relationship with the b value. J. Geophys. Res. 108 (B1), 2044. doi:10.1029/2001JB000680. Amitrano, D., Grasso, J.-R., Hantz, D., 1999. From diffuse to localised damage through elastic interaction. Geophys. Res. Lett. 26, 2109–2112. Argon, A.S., 1974. Statistical aspects of fracture. Compos. Mater. 5, 153–189. Ashby, M.F., Hallam, S.D., 1986. The failure of brittle solids containing small cracks under compressive stress states. Acta metall. 34, 497–510. Ashby, M.F., Sammis, C.G., 1990. The damage mechanics of brittle solids in compression. Pure Appl. Geophys. 133, 489–521. Batdorf, S.B., 1977. Some approximate treatments of fracture statistics for polyaxial tension. Int. J. Fract. 13, 5–10. Bazˇant, Z.P., Planas, J., 1997. In: Fracture and Size Effect in Concrete and Other Quasi-Brittle Materials. CRC Press, Boca Raton, p. 640. Boutt, D.F., McPherson, B.J.O.L., 2002. Simulation of sedimentary rock deformation: Lab-scale model calibration and parameterization. Geophys. Res. Lett. 29 (4), 1054. doi:10.1029/ 2001GL013987. Brace, W.F., 1961. Dependence of fracture strength of rocks on grain size. In: Proc. 4th Symp. Rock Mech., pp. 99–103. Brace, W.F., Silver, E., Hadley, K., Goetze, C., 1972. Cracks and pores—a close look. Science 178, 163–165. Budiansky, B., O’Connell, R.J., 1976. Elastic moduli of a cracked solid. Int. J. Solids Struct. 12, 81–97. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Evans, A.G., 1978. A general approach for the statistical analysis of multiaxial fracture. J. Am. Ceram. Soc. 61, 302–308. Evans, A.G., Langdon, T.G., 1976. Structural ceramics. Prog. Mater. Sci. 21, 174–441. Fang, Z., Harrison, J.P., 2002. Development of a local degradation approach to the modelling of the brittle fracture in heterogeneous rocks. Int. J. Rock Mech. Min. Sci. 39, 443– 457. Fredrich, J., Evans, B., Wong, T.-f., 1990. Effects of grain size on brittle and semi-brittle strength, implications for micro-

mechanical modelling of failure in compression. J. Geophys. Res. 95, 10907–10920. Fredrich, J.T., Menendez, B., Wong, T.-f., 1995. Imaging the pore structure of geomaterials. Science 268, 276–279. Freudentahl, A.M., 1968. Statistical approach to brittle fracture. In: Liebowitz, H. (Ed.), Fracture, vol. II. Academic Press, NY, Chapter 6. Gumbel, E., 1958. Statistics of Extremes. Columbia University Press, New York, p. 375. Gupta, V., Bergstro¨m, J., 1998. Compressive failure of rocks by shear faulting. J. Geophys. Res. 103 (B10), 23875–23895. Hadley, K., 1976. Comparison of calculated and observed crack densities and seismic velocities in Westerly granite. J. Geophys. Res. 81, 3484–3494. Handin, J., 1969. On the Coulomb–Mohr failure criterion. J. Geophys. Res. 74, 5343–5348. Hazzard, J.F., Young, R.P., Maxwell, S.C., 2000. Micromechanical modeling of cracking and failure in brittle rocks. J. Geophys. Res. 105, 16683–16697. Hoek, E., Brown, E.T., 1980. Underground Excavation Engineering. Institute of Mining and Metallurgy, London. Horii, H., Nemat-Nasser, S., 1986. Brittle failure in compression, splitting, faulting and brittle-ductile transition. Philos. Trans. Roy. Soc. London 319, 337–374. International Society for Rock Mechanics Commission on Testing Methods, 1988. Suggested method for determining fracture toughness of rock (F. Ouchterlony, Working Group Coordinator). Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 25, 71–96. Jaeger, J.C., Cook, N.G.W., 1979. Fundamentals of Rock Mechanics, third ed. Chapman and Hall, p. 593. Kachanov, L.M., 1986. Introduction to Continuum Damage Mechanics. Martinus Nijhoff, Dordrecht, p. 135. Katsman, R., Aharonov, E., Scher, H., 2005. Numerical simulation of compaction bands in high-porosity sedimentary rock. Mech. Mater. 37, 143–162. Kemeny, J.M., Cook, N.G.W., 1991. Micromechanics of deformation in rocks. In: Shah, S.P. (Ed.), Toughening Mechanisms in Quasi-Brittle Materials. Kluwer Academic, pp. 155–188. Krajcinovic, D., 1989. Damage mechanics. Mech. Mater. 8, 117–197. Krajcinovic, D., Vujosevic, M., 1998. Strain localization—short to long correlation length transition. Int. J. Solids Struct. 35, 4147–4166. Kranz, R.L., 1983. Microcracks in rocks, a review. Tectonophysics 100, 449–480. Lawn, B., 1993. Fracture of Brittle Solids, second ed. Cambridge University Press, Cambridge, p. 378. Li, L., Holt, R.M., 2002. Particle scale reservoir mechanics. Oil Gas Sci. Technol.-Rev. IFP 57, 525–538. Liu, H.Y., Roquete, M., Kou, S.Q., Lindqvist, P.-A., 2004. Characterization of rock heterogeneity and numerical verification. Eng. Geol. 72, 89–119. Lockner, D.A., 1995. Rock failure. In: Ahrens, T.J. (Ed.), Rock Physics and Phase Relations. A Handbook of Physical Constants. American Geophysical Union, Washington, DC, pp. 127–147. Lockner, D.A., Madden, T.R., 1991. A multiple crack model of brittle fracture. 1. Non-time-dependent simulations. J. Geophys. Res. 96, 19623–19642. Lockner, D.A., Byerlee, J.D., Kuksenko, V., Ponomarev, A., Sidorin, A., 1992. Observations of quasistatic fault growth from acoustic emissions. In: Evans, B., Wong, T.-f. (Eds.),

T.-f. Wong et al. / Mechanics of Materials 38 (2006) 664–681 Fault Mechanics and Transport Properties of Rocks. Academic, San Diego, CA, pp. 3–32. McClintock, F.A., Argon, A.S., 1966. Mechanical Behavior of Materials. Addison-Wesley, Reading, MA, p. 770. Mene´ndez, B., Zhu, W., Wong, T.-f., 1996. Micromechanics of brittle faulting and cataclastic flow in Berea sandstone. J. Struct. Geol. 18, 1–16. Moore, D.E., Lockner, D.A., 1995. The role of microcracking in shear-fracture propagation in granite. J. Struct. Geol. 17, 95– 114. Myer, L.R., Kemeny, J., Cook, N.G.W., Ewy, R., Suarez, R., Zheng, Z., 1992. Extensile cracking in porous rock under differential compressive stress. Appl. Mech. Rev. 45, 263– 280. Olsson, W.A., 1974. Grain size dependence of yield stress in marble. J. Geophys. Res. 32, 4859–4862. Paterson, M.S., Wong, T.-f., 2005. Experimental Rock Deformation—The Brittle Field, second ed. Springer-Verlag, New York, p. 348. Potyondy, D.O., Cundall, P.A., Lee, C., 1996. Modeling rock using bonded assemblies of circular particles. Proc. North Am. Rock Mech. Symp. 2, 1937–1944. Steacy, S.J., Sammis, C.G., 1992. A damage mechanics model for fault zone friction. J. Geophys. Res. 97, 587–594. Tang, C.A., 1997. Numerical simulation on progressive failure leading to collapse and associated seismicity. Int. J. Rock Mech. Min. Sci. 34, 249–262. Tang, C.A., Liu, H., Lee, P.K.K., Tsui, Y., Tham, L.G., 2000. Numerical studies of the influence of microstructure on rock failure in uniaxial compression—Part I: effect of heterogeneity. Int. J. Rock Mech. Min. Sci. 37, 555–569.

681

Tapponier, P., Brace, W.F., 1976. Development of stress-induced microcracks in Westerly granite. Int. J. Rock Mech. Min. Sci. 13, 103–112. Turcotte, D.L., 1992. Fractals and Chaos in Geology and Geophysics. Cambridge U.P., Cambridge, p. 221. Van Mier, J.G.M., Van Vliet, M.R.A., Wang, T.K., 2002. Fracture mechanisms in particle composites, statistical aspects in lattice type analysis. Mech. Mater. 34, 705–724. Walsh, J.B., 1965. The effect of cracks on the uniaxial elastic compression of rocks. J. Geophys. Res. 70, 399–411. Weibull, W., 1951. A statistical distribution function of wide applicability. J. Appl. Mech. 18, 293–297. Wong, T.-f., 1982. Micromechanics of faulting in Westerly granite. Int. J. Rock Mech. Min. Sci. 19, 49–64. Wong, T.-f., 1985. Geometric probability approach to the characterization and analysis of microcracking in rocks. Mech. Mater. 4, 261–276. Wong, R.H.C., Chau, K.T., Wang, P., 1996. Microcracking and grain size effect in Yuen Long marbles. Int. J. Rock Mech. Min. Sci. 33, 479–485. Wu, X.Y., Baud, P., Wong, T.-f., 2000. Micromechanics of compressive failure and spatial evolution of anisotropic damage in Darley Dale sandstone. Int. J. Rock Mech. Min. Sci. 37, 143–160. Xu, X.H., Ma, S.P., Xia, M.F., Ke, F.J., Bai, Y.L., 2004. Damage evaluation and damage localization of rock. Theor. Appl. Fract. Mech. 42, 131–138. Yuen, K.M., 1990. Characteristics and petrography of the Carboniferous Yuen Long marble. Karst Geol. H. K., Geol. Soc. H. K. Bull. 4, 73–87.