- Email: [email protected]
Process zone and size effect in fracture testing of rock
International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Process zone and size effect in fracture testing of rock Ali Fakhimi a,b,n, Ali Tarokh a a b
Department of Mineral Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA School of Civil and Environmental Engineering, Tarbiat Modares University, Tehran, Iran
a r t i c l e i n f o
abstract
Article history: Received 5 December 2011 Received in revised form 13 June 2012 Accepted 19 December 2012 Available online 30 January 2013
Fracture of rock involves formation of a localized region of damage or cohesive process zone, which controls size effects on strength and stability. Therefore, any attempt in predicting size of process zone is of prime importance in fracture study of quasi-brittle materials. In this paper, an approximate theoretical relationship between size of process zone and specimen size is proposed. The appropriateness of this relationship is examined by conducting a discrete element analysis of rock fracture. A softening contact bond model is used to study the process zone around a notch tip in three point bending tests. The numerical analysis is utilized to obtain the nominal tensile strength, apparent fracture toughness, and width of process zone. It is shown that the apparent fracture toughness is a function of the specimen size, and that the change in nominal tensile strength with specimen size can be captured by Bazant’s size effect law. In addition, both the theoretical arguments and the numerical results suggest that the inverse of width of the crack tip process zone has a linear relationship with the inverse of specimen size. The numerical results indicate a stronger relationship between width of process zone and specimen size for a material with a larger value of characteristic size. On the other hand, for a brittle material, specimen size has a small or no impact on the size of the process zone. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Process zone Specimen size Numerical modeling Fracture toughness
1. Introduction Damaged or process zone around a crack tip is an important subject in fracture study of quasi-brittle materials. This nonlinear region influences both the strength and stability of a rock structure, and is responsible for deviation of rock behavior from that dictated by linear elastic fracture mechanics. Even though there is a general agreement on the existence of a micro-fracture or process zone around a crack or notch tip, it is not clear if size of process zone is a material property. Li and Marasteanu [1] investigated the effect of aggregate size on fracture process zone of asphalt mixture at low temperature by means of bending tests and Acoustic Emission (AE) measurements. They claimed that the fracture process zone is a material property whose size may be dictated by many factors such as grain size and preexisting cracks or voids in the material. Zietlow and Labuz [2] claimed that size of process zone does not vary significantly for different size beams of the same material. Their evidence was the size of process zone in Berea sandstone and Sioux quartzite. On the other hand, for Rockville granite, sizes of the process zones for two beam sizes were significantly different. They stated that the issue of the effect of beam size on process zone dimensions for a given rock has not been resolved yet and needs further investigation.
n Corresponding author at: Department of Mineral Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA. Tel.: þ1 575 835 6577. E-mail address: [email protected] (A. Fakhimi).
1365-1609/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijrmms.2012.12.044
Zhang and Wu [3] reported that length of process zone is not a material parameter and is greatly influenced by specimen size. They conducted several three point bending tests on only one specimen size but with different notch depths for concrete. In their research, they assumed that the process zone size reaches a maximum value when the applied load approaches its peak value and that the length of the fracture process zone depends greatly on the length of the ligament; as the notch depth increases, the length of process zone decreases. Otsuka and Date [4] examined the influence of specimen size on process zone dimensions for concrete. They showed that for specimens of identical maximum aggregate size and initial notch length to depth ratio, the size (width and length) of the fracture process zone increases as the specimen size increases. Their results indicated that the increasing rate of the width of the process zone was smaller than that of the specimen size whereas the increasing rate of the length of the process zone was much larger than that of the specimen size. They also noted that for full development of process zone, an appropriate ligament length should be taken into consideration. Bazant and Kazemi [5] used Bazant’s size effect law to derive an equation for equivalent elastic crack length as a measure of length of process zone. However, there is a need to obtain a relationship between the actual size of process zone around a crack tip and the specimen size and to verify whether size of process zone can be considered as an intrinsic material property. The purpose of this paper is to study size of fracture process zone using both theoretical and numerical approaches. Through an approximate theoretical argument, it is shown that width of
96
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
process zone is in general a function of specimen size. A discrete element numerical model is used to verify the usefulness of the theoretical findings. Finally, published experimental data are utilized to support and justify the theoretical and numerical results.
2. Theoretical study Consider two rock specimens I and II in three point bending test (Fig. 1). The applied load in both specimens is assumed to be small enough so that both specimens are in linear elastic regime. The sizes of the specimens are proportional in the x1 x2 plane, i.e. for any point (A) with coordinates x1 and x2 in the smaller specimen, there is a corresponding point (A0 ) in the larger specimen that has the coordinates X1 and X2 such that X i ¼ mxi
ð1Þ
Hence, the size of specimen II is m times the size of specimen I. It is assumed that both specimens have the same thickness (b). The smaller specimen under load P is assumed to be in static equilibrium with the stress sij(x1,x2), strain eij(x1,x2) and displacement ui(x1,x2) fields at a typical point A. We seek for a load P0 applied to specimen II and a strain field (e0 ij) such that the stress (s0 ij) and displacement components (u0 i) at point A0 are related to the corresponding values at point A for specimen I through the following relationships: 0 ij
0 ij ðX 1 ,X 2 Þ ¼
s ¼s
sij ðx1 ,x2 Þ
ð2aÞ
u0i ¼ u0i ðX 1 ,X 2 Þ ¼ mui ðx1 ,x2 Þ
ð2bÞ
The above displacement field satisfies the displacement boundary conditions for specimen II. It is easy to see that the stress components s0 ij satisfy the equilibrium equations (assuming zero body forces): @s0ij @X j
¼
@sij @xk @sij 1 1 @sij ¼ d ¼ ¼0 @xk @X j @xk m kj m @xj
ð3Þ
x2
R
A
w
it is concluded that e0 ij ¼ eij. It is an easy task to show that P0 ¼mP by finding the integral of the vertical traction component along the corresponding semi-circles in Fig. 1 using the fact that sij ¼ s0 ij at the corresponding points along the circles. Above discussion suggests the well known fact that the compliances of beams I and II are identical: C0 ¼ C ¼
d0 P0
¼
md d ¼ mP P
ð5Þ
where d and d0 are the load-point displacements for specimens I and II, respectively. Now consider the scenario that damaged zones have been created around the notch tips in both specimens I and II (Fig. 1). In the remaining of this section we assume that for values of m41 but close to 1, the similitude discussed above still will be approximately valid. An immediate consequence of this assumption is that the width of the damaged zone in specimen II (W0 ) is m times that in specimen I (W0 ¼ mW). On the other hand, using Bazant’s size effect equation [6] and the nominal tensile strength of a beam (sN ¼3PmaxL/2bD2), the relationship between the peak load and the specimen height (D) can be obtained: Pmax ¼
bDBst qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi b 1 þ D=D0
ð6Þ
in which st is the tensile strength of the specimen, b ¼3L/2D, and B and D0 are constants that depend on the material and specimen geometry. To obtain the width of the fracture process zones in Fig. 1, two extreme situations are considered. For the first situation, assume that the specimen sizes are very small. In particular, for very small values of D and D0 ¼mD, using Eq. (6) it is easy to see that plastic behavior of material prevails which leads to no or negligible size effect on strength, i.e. P0max ffim Pmax
ð7Þ
Eq. (8) suggests a linear relationship between the width of process zone at peak load and specimen size, i.e.
D
WðDÞ ¼ lD
ao
x1
L mP
mR
Let the applied force P in Fig. 1a be equal to Pmax/Om. From Eq. (10) it follows that P0 ¼P0 max ¼PmaxOm¼ mP. These loading conditions for beams I and II result in the following approximate equation concerning the width of the process zone:
A' mD
w' mao
ð9Þ
where l is a constant for a given material and specimen shape. The second situation is when D is large. In this case, an approximate ratio of P0 max to Pmax is obtained from Eq. (6): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mDbBst =b 1 þ mD=D0 pffiffiffiffiffi P0max qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð10Þ ¼ ffi ffi m Pmax DbBst =b 1 þ D=D0
X2
(II )
ð4Þ
Therefore, for small specimens, the developed process zone widths at peak loads in beams I (W) and II (W0 ) are related through the following linear equation: ð8Þ W 0 P0 ¼ P 0max ,D0 ¼ mD ¼ mW ðP ¼ Pmax ,D ¼ DÞ
P
(I)
Considering the equation @u0i @ðmui Þ @xk @u ¼ ¼ i @xk @X j @X j @xj
W 0 9P0
X1
mL Fig. 1. Two geometrically similar beams I and II that are made of the same material.
max
¼ mW9Pmax =pffiffiffi ¼ mW 1 m
ð11Þ
where W1 is the width of the damaged zone in specimen I at a load equal to Pmax/Om. For large specimens, linear fracture mechanics gradually prevails and it is expected that the width of the process zone would change slowly with the specimen size. This fact together with the assumption that m is close to
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
1 suggests that the widths of the process zones for specimens I and II at peak load must be almost identical. Therefore, Eq. (11) leads to W9P ¼ Pmax ¼ mW9P ¼ Pmax =pffiffiffi m
ð12Þ
The solution for Eq. (12) is W ¼ zðD,b, st ÞP2
2
b D2 B2 s2t b ½1þ D=D0 2
ð14Þ
when specimen size approaches infinity, W is expected to approach a constant value. To fulfill this requirement, the function z(D,b,st) is assumed to be proportional to 1/( Ds2t b2) in which the tensile strength in the dominator has been added for convenience. After substitution in Eq. (14), the final relationship between W and D is obtained: W ¼C
DB2 W 1D ¼ b ½1 þ D=D0 D0 ½1 þ D=D0 2
The peak load data from numerical tests and Eq. (18) will be used to study the phenomenon of size effect on strength in fracture testing of three point bending specimens. Furthermore, the approximate validity of Eq. (15), and therefore the soundness of assumptions in deriving this equation will be shown by numerical analysis of rock fracture in the next sections.
ð13Þ
in which z(D,b,st) is in general a function of D and b for geometrically similar structures with a constant a0/D; z is a function of material strength (such as tensile strength) as well. Note that Eq. (13) is valid only for m values close to 1, i.e. P in this equation must be close to the peak load. Eq. (13) implies that for loads close to the peak load, the width of the process zone grows much faster than the increase in the applied load; majority of micro-cracks at the notch tip are developed when the applied load approaches its peak value. This approximate theoretical finding is consistent with the experimental observations [7,8]. The width of the process zone at the peak load can be obtained as a function of the specimen size if Eqs. (6) and (13) are combined: W ¼ zðD,b, st Þ
3. Numerical model A discrete element model was used to simulate the failure process [10,11]. The model consisted of many rigid circular cylinders or disks that interact through normal and shear springs to simulate elasticity. In order to withstand tensile and deviatoric stresses, cylinders or disks were bonded together at contact points. The micro-mechanical constants at a contact point in this model are normal and shear spring constants (Kn,Ks), normal and shear bonds (nb and sb) and friction coefficient (m). To capture softening behavior of a quasi-brittle material, a softening bond feature was implemented in the numerical model. In this model, normal bond at a contact point is assumed to reduce linearly after the peak tensile load (Fig. 2). This introduces a new microscopic constant in the model, the slope in the post-peak region of the normal forcenormal displacement between two cylinders in contact. This slope is indicated by Knp in Fig. 2a. No modification in the relationship of shear force (Fs)-relative shear displacement (us) of a contact is assumed in this simple model. The relationship between Fs and us is shown in Fig. 2b. The loading–unloading paths for both normal and shear contact forces are shown in Fig. 2 with arrows. Softening in shear was not included because this effect is pronounced for loading only under significant (410%qu) mean stress.
ð15Þ
Fn(Tension)
in which C is a dimensionless constant for a given material and geometry. The term WN in Eq. (15) is the width of process zone when the specimen size approaches infinity; it is equal to W1 ¼
CD0 B
b2
97
nb
C D
kn
2
knp
ð16Þ 1
B
As a consequence of Eq. (15), the width of the process zone for specimen II can be predicted if this width is specified for specimen I by using the following equation m½1 þ D=D0 ð17Þ W 0 ¼ W9mD ¼ W9D 1þ mD=D0 It is interesting to note that even though Eq. (15) was derived for large specimens, inserting a small D in this equation results in a linear relationship between width of process zone and specimen height, which is consistent with Eq. (9) that was obtained for small size specimens. Consequently, Eq. (15) shows an approximate relationship between W and D for both large and small specimens even though the prediction of the equation is expected to deviate to some extent from physical observations for intermediate beam sizes. In evaluation of the numerical fracture test data, stress intensity factor at peak load or apparent fracture toughness (KICA) defined as follows will be used: pffiffiffiffi a0 ð18Þ K ICA ¼ sN Df D in which f is a function of the ratio of initial notch length to beam height. The expression for f(a0/D) has been reported in [9]: pffiffiffi a½1:99að1aÞ 2:153:93a þ2:7a2 ð19Þ f ðaÞ ¼ 1:5 ð1 þ 2aÞð1aÞ
1
E
1
Kn
1
F
Un
kn
A Fn(Compression) Fs
Sb
B
μFn
C
A
D
Us
Fig. 2. Relationships in the softening contact bond model: (a) normal force and normal displacement and (b) shear force and shear displacement.
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
The discrete element analysis was performed using a computer program that was developed by Fakhimi [12,13]. CA2 is a hybrid discrete element-finite element code for large deformation analysis of solid materials. The details of sample preparation and calibration procedure have been described in the past [14]. The numerical model was calibrated to simulate the mechanical behavior of Berea sandstone. The grains in this sandstone range from 0.1 to 0.8 mm. The mechanical properties of the rock are the elastic modulus of E¼ 14 GPa, Poisson’s ratio of n ¼0.32, uniaxial compressive strength of qu ¼55–60 MPa, and bending tensile strength (80 240 mm2 beam) of st ¼8.6 MPa [2]. After calibration of the numerical model, uniaxial compressive (40 80 mm2 sample) and bending tensile tests (80 240 mm2 beam) were conducted to verify the accuracy of the model. An elastic modulus of 13.3 GPa, Poisson’s ratio of 0.19, uniaxial compressive strength of 60.5 MPa, and tensile strength of 8.7 MPa were obtained that are in close agreement with the corresponding rock properties. The micro-mechanical properties that were obtained through the calibration procedure are Kn ¼22.0 GPa, Ks ¼5.5 GPa, nb ¼2800 N/m, sb ¼12300 N/m, m ¼ 0.5, s0 (genesis pressure)¼2.2 GPa and Knp ¼1.83 GPa (Kn/Knp ¼12). The circular particles or disks radii (R) were assumed to have a uniform random distribution with a range of 0.27–0.33 mm (Rave ¼0.3 mm).
The peak load, nominal bending tensile strengths, apparent fracture toughness and apparent characteristic size (lchA ¼(KICA/ st)2) for numerical specimens with a0/D ¼0.375 are shown in Table 2. The apparent fracture toughness values in this table were obtained using Eqs. (18) and (19). The variation of apparent fracture toughness with specimen size for a0/D ¼0.375 is shown in Fig. 4a. Notice that as the material becomes more brittle, the value of apparent fracture toughness becomes less size dependent. The apparent characteristic size values in Table 2 and characteristic size (lch ¼(KIC/st)2) can be used to compare brittleness of specimens of the same size; specimens with larger lchA are associated with higher Kn/Knp values and are less brittle. To calculate the size effect parameters B and D0, two different linear regression analyses (methods I and II) have been suggested in [6]. In method I, Eq. (6) is written with respect to the nominal tensile strength (sN) and is linearized (y ¼Axþ C) by defining x¼D and y¼1/s2N. In the II approach, the apparent fracture toughness is written in the following form: pffiffiffiffi Bst K ICA ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Df a0=D 1 þ D=D0
ð20Þ
which can be modified to a linear equation if we choose x¼1/D and y¼(1/KICA)2:
3.1. Numerical tests Three point bending tests were conduced on the numerical samples 20 60, 40 120, 80 240, 160 480, and 320 960 mm2 in size. The applied vertical velocity at the beam center (at top surface) was 2.5 10 10 m/cycle to obtain a quasi-static solution. All specimens had a notch at their center. Three different notch length (a0) to beam height (D) ratios of 0.125, 0.25, and 0.375 were used. The slope of the softening line (Knp) was modified to obtain different quasi-brittle synthetic materials. Five different Kn/Knp values of 0, 10, 20, 50, and 100 were used while other micro-mechanical constants (Kn, Ks, nb, sb, s0, and R) were left unchanged. Uniaxial tensile and compressive tests on samples 40 80 mm2 in dimension resulted in strengths values reported in Table 1. Notice that as Kn/Knp increases, the ratio of the compressive to tensile strength reduces. This is expected because with increase in Kn/Knp, the model becomes less brittle. The numerical bending tests were quite time consuming considering the relative sizes of the beams and the average particle radius. For example, for the 320 960 mm2 beam, more than 1,200,000 particles were needed. The numerical analysis for this beam took more than 4 weeks to finish for Kn/Knp ¼100 on a desk top computer.
y¼
1 B2 st 2 f
2
xþ
1
ð21Þ
2
B2 st 2 f D0
50 Kn/Knp=100
45 40 35 Load (kN)
98
Kn/Knp=50
30 25 20
Kn/Knp=20
15 10
Kn/Knp=10 Kn/Knp=0
5 0 0.00
0.05
0.10 0.15 Displacement (mm)
0.20
0.25
300
3.2. Numerical bending test results
Kn/Knp=100 250 200 Load (kN)
The load-displacement curves for all beam sizes (five different sizes from 20 60 to 320 960 mm2), three different notch to beam height ratios (a0/D¼0.125, 0.25, 0.375), and five different synthetic materials (Kn/Knp ¼0, 10, 20, 50, 100) were obtained. Fig. 3 shows the typical load displacement curves for beams 20 60 and 160 480 mm2 in size and a0/D¼0.375. As expected, by increasing the Kn/Knp value, the material shows a more extensive non-linear hardening behavior as the peak load is gradually approached.
Kn/Knp=50
150 100 Kn/Knp=20
Table 1 Uniaxial tensile and compressive strengths of the numerical specimens.
50
Kn/Knp
st (MPa)
sc (MPa)
sc/st
0 10 20 50 100
4.2 5.6 7.0 8.6 10.7
57.9 63.1 64.8 66.3 73.1
13.8 11.4 9.2 7.7 6.8
Kn/Knp=10 Kn/Knp=0
0 0.0
0.1
0.2 0.3 0.4 Displacement (mm)
0.5
0.6
Fig. 3. Load vs. displacement of loading point for numerical specimens with Kn/Knp ¼ 0, 10, 20, 50, 100 and a0/D¼ 0.375: (a) 20 60 mm2 specimen and (b) 160 480 mm2 specimen.
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
Brittle
Table 2 Peak load, nominal tensile strength, apparent fracture toughness, and apparent characteristic size for numerical specimens with a0/D¼ 0.375.
Kn/Knp=10 Kn/Knp=20
4.0
20 60
40 120
80 240
160 480
320 960
Kn/Knp
Pmax (kN)
sN (MPa)
KICA (MPa m0.5)
lchA (mm)
Kn/Knp=50 3.5
0 10 20 50 100 0 10 20 50 100 0 10 20 50 100 0 10 20 50 100 0 10 20 50 100
10.80 13.44 19.36 30.79 45.90 13.94 21.89 30.33 53.73 79.98 21.36 33.00 47.16 96.51 149.21 34.63 47.56 74.82 153.87 242.69 45.45 63.19 110.63 253.37 387.25
2.43 3.02 4.36 6.93 10.33 1.57 2.46 3.41 6.04 9.00 1.20 1.86 2.65 5.43 8.39 0.97 1.34 2.10 4.33 6.83 0.64 0.89 1.56 3.56 5.45
0.42 0.53 0.76 1.21 1.80 0.39 0.61 0.84 1.49 2.22 0.42 0.65 0.93 1.89 2.93 0.48 0.66 1.04 2.14 3.37 0.45 0.62 1.09 2.49 3.80
10 9 12 20 28 9 12 14 30 43 10 13 18 48 75 13 14 22 62 99 11 12 24 84 126
Kn/Knp=100
3.0 KICA(MP a.m0.5)
Specimen size (mm2)
99
2.5 2.0 1.5 1.0 0.5 0.0 0
40
80
120
160 200 D (mm)
240
280
320
360
Strength criterion 1.0
The data from Table 2 together with fitting lines (Eq. (21)) have been plotted in Fig. 5. The numerical test data illustrate very good linear trends except for Kn/Knp ¼0 (brittle material) that shows greater scatter. The slopes and y-intercepts in Fig. 5 together with st in Table 1 are used to find B and D0 that are reported in Table 3. Table 3 includes the results for a0/D¼0.125 and 0.25 as well. It is important to realize that methods I and II regression analyses do not in general result in the same size effect parameters B and D0 as these regressions imply different weightings of the data points [6]. Notice that the II regression method has produced no unacceptable negative D0 values. Furthermore, the D0 values in method II show consistent increase as the material becomes less brittle which is expected; the brittleness of the material, i.e. D/D0 [6] must decrease with increase in the Kn/Knp value (for a fixed D). On the other hand, the results of method I regression analysis have produced some unacceptable data. The size effect parameters in Table 3 have been used to show how the nominal tensile strength changes with the specimen size (Fig. 4b). Note that the data points for more ductile materials (higher Kn/Knp values) are close to the line representing the strength of material theory. Conversely, more brittle materials follow more closely the LEFM theory. The main goal of this work was to study the width of the fracture process zone around the notch tip. To this end, the process zone in the numerical model is defined with the contact points between circular particles that are in the post-peak regime (e.g. point D in Fig. 2a). These damaged contact points for two specimen sizes 80 240 mm2 and 320 960 mm2 and Kn/Knp of 20 and 50 are shown in Fig. 6. The damaged contacts are shown in green while the actual sharp crack in the middle of the process zone has been shown in red in the figure. Fig. 6 suggests that the width of the process zone is greater for a more ductile material (higher Kn/Knp values). Furthermore, as the size of the specimen increases, the width of the process zone increases as well. Eq. (15) is reformulated to facilitate linear regression analysis of the width of process zone data as follows: 1 1 1 D0 ¼ þ W W1 W1 D
ð22Þ
σN/(Bσt)
LEFM
2 1 Brittle Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting 0.1 0.01
0.1
1 D/D0
10
100
Fig. 4. (a) Apparent fracture toughness vs. specimen size and (b) nominal tensile strength vs. specimen size. The ratio of notch length to specimen height (a0/D) is 0.375.
The variations of 1/W vs. 1/D for a0/D ¼0.375 and different Kn/Knp values are shown in Fig. 7. The linear trends of the data in Fig. 7 confirm the appropriateness of Eq. (15) and support the logic used in the theoretical section of the paper. Using the slopes and y-intercepts of the lines in Fig. 7 and Eq. (22), the D0 values and the widths of the process zones for very large specimens (WN) can be obtained. This information is reported in Table 4. Notice that similar to the situation with regression analysis in methods I and II that did not result in identical D0 values, the D0 values from this III regression method are not the same as those in Table 3. The values of WN and D0 in Table 4 have been used in Eq. (15) to obtain the variation of process zone width with specimen size (Fig. 8a). The good match between the fitting curves and the numerical data implies that Eq. (15) can accurately represent the dependence of the size of process zone on the specimen size. Note that as the material becomes more brittle, the specimen size imposes less impact on the size of the process zone. In Fig. 8b, the width of process zone together with the fitting curves (Eq. (15)) is shown using D0 values from size effect analysis (D0 values from Table 3, method II). Except for Kn/Knp ¼100, the
100
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
7
y = 0.0237x + 4.9795 R2= 0.2855
(1/KICA)2 MPa-2m-1
6
W=3.6 mm W=5 mm
5 y = 0.0244x + 2.2294 R2= 0.8454
4
Brittle Kn/Knp=10 Kn/Knp=20 Linear (Brittle) Linear (Kn/Knp=10) Linear (Kn/Knp=20)
3 y = 0.0187x + 0.8511 R2= 0.9545
2 1
Notch Notch
0 0
10
20
30
40
50
60
W=39 mm
W=22 mm
1/D (m-1) 0.8 y = 0.011x + 0.1447 R2= 0.992
(1/KICA)2 MPa-2m-1
0.7 0.6
Notch
Notch
0.5 y = 0.0052x + 0.0569 R2= 0.9906
0.4
Fig. 6. Width of the process zone around the notch tip for a0/D¼ 0.375: (a) 80 240 mm2 specimen with Kn/Knp ¼ 20, (b) 320 960 mm2 specimen with Kn/Knp ¼ 20, (c) 80 240 mm2 specimen with Kn/Knp ¼ 50 and (d) 320 960 mm2 specimen with Kn/Knp ¼50. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
0.3 Kn/Knp=50 Kn/Knp=100 Linear (Kn/Knp=50) Linear (Kn/Knp=100)
0.2 0.1 0.0 0
10
20
30
40
50
Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Linear (Kn/Knp=20) Linear (Kn/Knp=50) Linear (Kn/Knp=100)
0.4
60
1/D (m-1) Fig. 5. Linear regression analysis to obtain size effect parameters for beams with a0/D ¼0.375 and (a) Kn/Knp ¼0, 10, 20 and (b) Kn/Knp ¼50, 100.
1/W (mm-1)
Table 3 Size effect parameters for different a0/D values. a0/D
0.125
0.25
0.375
Kn/Knp
0 10 20 50 100 0 10 20 50 100 0 10 20 50 100
Method I
0.3
Method II
B
D0 (mm)
B
D0 (mm)
2.216 1.240 1.096 1.496 1.830 0.779 0.888 1.007 1.208 1.281 1.073 NA 0.729 0.839 1.024
8.9 26.4 56.7 150.0 260.0 39.1 21.6 28.7 92.0 176.7 6.8 0.4 32.0 95.5 103.8
3.283 1.873 2.025 1.960 2.183 2.147 1.181 1.289 1.384 1.463 1.256 0.934 0.846 0.895 1.048
3.8 9.4 12.5 50.2 75.9 3.8 11.4 16.8 51.2 81.5 4.8 10.9 22.0 76.0 91.4
prediction of the fitting equation is satisfactory suggesting that the D0 values from both methods II and III together with Eq. (15) can provide reasonable prediction of the evolution of process zone width as the specimen size is changed.
y = 3.3947x + 0.2181 R2 = 0.8769
0.2 y = 2.1531x + 0.0197 R2 = 0.9936 0.1 y = 1.3986x + 0.0099 R2 = 0.992
0.0 0
0.01
0.02
0.03
0.04
0.05
0.06
-1)
1/D (mm
Fig. 7. Linear regression analysis to obtain WN and D0. The ratio of notch length to specimen height (a0/D) is 0.375.
Table 4 Data from linear regression analysis of Eq. (15) and (22)). a0/D
Kn/Knp
Slope
y-intercept (mm 1)
WN (mm)
D0 (mm)
0.125
20 50 100 20 50 100 20 50 100
1.8318 1.5078 1.009 2.8469 1.8925 1.0949 3.3947 2.1531 1.3986
0.1936 0.0169 0.0062 0.1992 0.0179 0.0101 0.2181 0.0197 0.0099
5.2 59.2 161.3 5.0 55.9 99.0 4.6 50.8 101.0
9.5 89.2 162.7 14.3 105.7 108.4 15.6 109.3 141.3
0.25
0.375
4. Experimental data and discussion of the results Preceding discussions illustrate that the theoretical equation (Eq. (15)) for the size of fracture process zone around a notch tip can provide results that are consistent with those from discrete element analysis. To display further supports for findings of this study, some experimental observations regarding the width of process zone have been compiled in Fig. 9. The work of Otsuka
and Date [4] was concentrated on concrete. They performed tensile tests on different sample sizes. X-ray and three-dimensional Acoustic Emission (AE) techniques were used to investigate the process zone around the notch tip. It should be pointed out that in the work of Otsuka and Date, what is called damaged zone in Fig. 9 is in fact the Fracture Core Zone (FCZ) defined by these authors as a region
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting Kn/Knp=20 Fitting Kn/Knp=50 Fitting Kn/Knp=100
80 70
W (mm)
60
35
Berea sandtone, after Zietlow & Labuz (1998)
50 40 30 20
Charcoal granite, after Zietlow & Labuz (1998) Rockville granite, after Zietlow & Labuz (1998)
25 20 15 10 5
10
0
0 0
40
80
120
160 200 D (mm)
240
280
320
360
80 70
50 40 30 20 10 0 40
80
120
100
200
300
400
500
600
700
800
Fig. 9. Some published information regarding the width of process zone as a function of specimen size.
60
0
0
Specimen size (mm)
Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting Kn/Knp=20 Fitting Kn/Knp=50 Fitting Kn/Knp=100
90
W (mm)
Concrete, after Otsuka & Date (2000)
30 Width of damaged zone (mm)
90
101
160 200 D (mm)
240
280
320
360
Fig. 8. Width of process zone vs. specimen size. The fitting curves are based on D0 values from (a) method III regression analysis and (b) method II regression analysis. The ratio of notch length to beam height (a0/D) is 0.375.
around the notch tip and within the process zone that is subjected to more damage and releases a higher percentage of acoustic energy. Zietlow and Labuz [2] studied fracture process zone in different rock beams using AE measurement. Fig. 9 clearly shows the dependence of the width of process zone to specimen size for concrete, Charcoal granite, and Rockville granite. On the other hand, for Berea sandstone, no change in the size of process zone is observed by increasing the specimen size. This indicates that the contact points between grains of sandstone must have a relatively brittle behavior (small Kn/Knp values around 10–20, see Fig. 8). For granite specimens, the sizes of process zones are greater and in these cases, the specimen size has a greater impact on the size of the process zone. In summary, it appears that the numerical, theoretical, and experimental data all confirm that in a material with a larger process zone, the size of the process zone is expected to be more dependent on specimen size (compare the slopes of the lines passed through the data points in Fig. 9). This observation should resolve the controversy in the literature regarding whether the size of process zone is a material property or not. In general, for quasi-brittle materials, process zone like fracture toughness is dependent on size. If the material is highly brittle, the size of fracture process zone shows no noticeable change as the specimen size is modified; fracture process zone appears as a material property within the specimen sizes typically used in laboratory testing. On the contrary, fracture process
zone may not appear as a material property for laboratory size specimens if the material shows a less brittle behavior. From comparison of Figs. 4 and 8a, it appears that the width of process zone ceases changing with specimen size when apparent fracture toughness approaches a constant value (true fracture toughness). Since measurement of fracture process zone is a more expensive and laborious procedure, this observation indicates that measurement of apparent fracture toughness can be utilized as a tool to find out if the size of the specimen is large enough to simulate a mature and fully developed size independent fracture process zone (WN). Furthermore, a useful feature of Eq. (22) is its potential to provide an estimate of WN. By performing tests on small laboratory size specimens and measuring the widths of process zones, WN can be estimated from Eq. (22) without a need to perform tests on large specimens. The numerical results reveal that the width of process zone is not only a function of specimen size but also a function of a0/D and that as the material becomes more ductile, the influence of a0/D on the width of process zone becomes more important (Fig. 10). The numerical results suggest that the size of process zone increases as a0/D reduces. Therefore, for a shorter ligament length, a smaller fracture process zone is expected (Fig. 10). Experimental tests on quasi-brittle materials confirm this observation at least for the length of fracture process zone [1,15]. Nevertheless, for very large specimens, LEFM is expected to prevail and WN should be the same irrespective of the ratio of notch length to beam height. The data in Table 4 are approximately consistent with this expectation except for a0/D ¼0.125 and Kn/Knp ¼100 (see Fig. 10). It appears that for this relatively ductile synthetic material, even for large specimens, LEFM may not completely govern the behavior. This might be suggesting that for this synthetic material, the size of the process zone is much greater than the region of singular stress field dictated by stress intensity factor or fracture toughness leading to invalid asymptotic prediction of LEFM theory. Note that in discrete element analysis of the process zone, random location of the particles can have some influence on the size of the damaged zone. This issue was not investigated in this study and is not expected to affect the general conclusions of the paper as many particles exist within the damaged region that can smooth out the data scatter. Only for very small structures that are made of large grains with respect of the size of process zone, the random location of particles can cause noticeable scatter of the data. This issue needs further investigation. Finally, it is important to realize that even though the main focus of this
102
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
brittle, the dependence of the process zone on specimen size and notch length is reduced. The results of numerical modeling indicate that dependence of process zone on specimen size disappears when the size of the specimen is large enough such that the apparent fracture toughness is equal to the true material toughness. Attention to the peak loads in the bending tests shows that the discrete element analyses were able to accurately mimic the phenomena of size effect in quasi-brittle materials. Consequently, this numerical technique with the softening bond model can be used as a reliable tool to study crack initiation, crack development, and size effect in rock structures that are predominantly subjected to mode I loading. References
Fig. 10. Width of process zone vs. specimen size and the influence of the ratio of notch length to beam height (a0/D). WN values are in mm. WN values that are smaller than 60 mm are for Kn/Knp ¼ 50 and the rest are for Kn/Knp ¼100 (see Table 4).
study was evaluation of the width of fracture process zone, the same idea is applicable for prediction of the length of process zone. Interested readers are referred to [16] for further information.
5. Conclusion Based on Bazant’s size effect law, an approximate theoretical equation for size of process zone in quasi-brittle materials was proposed. It was shown that this equation results in a linear relationship between the inverse of width of process zone and the inverse of specimen size. A discrete element numerical model with tension softening was utilized to verify the appropriateness and accuracy of the theoretical equation. Synthetic materials with different degrees of ductility and sizes were produced and analyzed. The results obtained for the width of the process zone in the numerical simulations show consistency with the theoretical equation. In particular, it was found that the size of process zone is in general a function of specimen size and the ratio of notch length to specimen height. As the material becomes more
[1] Li X, Marasteanu M. The fracture process zone in asphalt mixture at low temperature. Eng Fract Mech 2010;77:1175–90. [2] Zietlow WK, Labuz JF. Measurement of the intrinsic process zone in rock using acoustic emission. Int J Rock Mech Min Sci 1998;35(3):291–9. [3] Zhang D, Wu K. Fracture process zone of notched three-point-bending concrete beams. Cem Concr Res 1999;29:1877–92. [4] Otsuka K, Date H. Fracture process zone in concrete tension specimen. Eng Fract Mech 2000;65:111–31. [5] Bazant ZP, Kazemi MT. Determination of fracture energy, process zone length and brittleness number from size effect, with application to rock and concrete. Int J Fract 1990;44:111–31. [6] Bazant ZP, Planas J. Fracture and size effect in concrete and other quasi-brittle materials. Boca Raton, FL: CRC Press; 1998. [7] Dai ST, Labuz JF. Damage and failure analysis of brittle materials by acoustic emission. J Mater Amer Soc Civ Eng 1997;9(4):200–5. [8] Jankowski LJ, Stys DJ. Formation of the fracture process zone in concrete. Eng Fract Mech 1990;36(2):245–53. [9] Strawley JE. Wide range stress intensity factor expressions for ASTM E399 standard fracture toughness specimens. Int J Fract 1976;12(2):475–6. [10] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329–64. [11] Lin Q, Fakhimi A, Haggerty M, Labuz JF. Initiation of tensile and mixed-mode fracture in sandstone. Int J Rock Mech Min Sci 2009;46:489–97. [12] Fakhimi A. Application of slightly overlapped circular particles assembly in numerical simulation of rocks with high friction angle. Eng Geol 2004;74: 129–38. [13] Fakhimi A. A hybrid discrete-finite element model for numerical simulation of geomaterials. Comput Geotech 2009;36:386–95. [14] Fakhimi A, Villegas T. Application of dimensional analysis in calibration of a discrete element model for rock deformation and fracture. Rock Mech Rock Eng 2007;40(2):193–211. [15] Wu Z, Rong H, Zheng J, Xu F, Dong W. An experimental investigation on the FPZ properties in concrete using digital image correlation technique. Eng Fract Mech 2011;78:2978–90. [16] Tarokh A. Discrete element simulation of fracture process zone in quasibrittle materials. MS thesis, Department of Mineral Engineering. Socorro, NM: New Mexico Institute of Mining and Technology, 2012.
Contents lists available at SciVerse ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Process zone and size effect in fracture testing of rock Ali Fakhimi a,b,n, Ali Tarokh a a b
Department of Mineral Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA School of Civil and Environmental Engineering, Tarbiat Modares University, Tehran, Iran
a r t i c l e i n f o
abstract
Article history: Received 5 December 2011 Received in revised form 13 June 2012 Accepted 19 December 2012 Available online 30 January 2013
Fracture of rock involves formation of a localized region of damage or cohesive process zone, which controls size effects on strength and stability. Therefore, any attempt in predicting size of process zone is of prime importance in fracture study of quasi-brittle materials. In this paper, an approximate theoretical relationship between size of process zone and specimen size is proposed. The appropriateness of this relationship is examined by conducting a discrete element analysis of rock fracture. A softening contact bond model is used to study the process zone around a notch tip in three point bending tests. The numerical analysis is utilized to obtain the nominal tensile strength, apparent fracture toughness, and width of process zone. It is shown that the apparent fracture toughness is a function of the specimen size, and that the change in nominal tensile strength with specimen size can be captured by Bazant’s size effect law. In addition, both the theoretical arguments and the numerical results suggest that the inverse of width of the crack tip process zone has a linear relationship with the inverse of specimen size. The numerical results indicate a stronger relationship between width of process zone and specimen size for a material with a larger value of characteristic size. On the other hand, for a brittle material, specimen size has a small or no impact on the size of the process zone. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Process zone Specimen size Numerical modeling Fracture toughness
1. Introduction Damaged or process zone around a crack tip is an important subject in fracture study of quasi-brittle materials. This nonlinear region influences both the strength and stability of a rock structure, and is responsible for deviation of rock behavior from that dictated by linear elastic fracture mechanics. Even though there is a general agreement on the existence of a micro-fracture or process zone around a crack or notch tip, it is not clear if size of process zone is a material property. Li and Marasteanu [1] investigated the effect of aggregate size on fracture process zone of asphalt mixture at low temperature by means of bending tests and Acoustic Emission (AE) measurements. They claimed that the fracture process zone is a material property whose size may be dictated by many factors such as grain size and preexisting cracks or voids in the material. Zietlow and Labuz [2] claimed that size of process zone does not vary significantly for different size beams of the same material. Their evidence was the size of process zone in Berea sandstone and Sioux quartzite. On the other hand, for Rockville granite, sizes of the process zones for two beam sizes were significantly different. They stated that the issue of the effect of beam size on process zone dimensions for a given rock has not been resolved yet and needs further investigation.
n Corresponding author at: Department of Mineral Engineering, New Mexico Institute of Mining and Technology, Socorro, NM 87801, USA. Tel.: þ1 575 835 6577. E-mail address: [email protected] (A. Fakhimi).
1365-1609/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijrmms.2012.12.044
Zhang and Wu [3] reported that length of process zone is not a material parameter and is greatly influenced by specimen size. They conducted several three point bending tests on only one specimen size but with different notch depths for concrete. In their research, they assumed that the process zone size reaches a maximum value when the applied load approaches its peak value and that the length of the fracture process zone depends greatly on the length of the ligament; as the notch depth increases, the length of process zone decreases. Otsuka and Date [4] examined the influence of specimen size on process zone dimensions for concrete. They showed that for specimens of identical maximum aggregate size and initial notch length to depth ratio, the size (width and length) of the fracture process zone increases as the specimen size increases. Their results indicated that the increasing rate of the width of the process zone was smaller than that of the specimen size whereas the increasing rate of the length of the process zone was much larger than that of the specimen size. They also noted that for full development of process zone, an appropriate ligament length should be taken into consideration. Bazant and Kazemi [5] used Bazant’s size effect law to derive an equation for equivalent elastic crack length as a measure of length of process zone. However, there is a need to obtain a relationship between the actual size of process zone around a crack tip and the specimen size and to verify whether size of process zone can be considered as an intrinsic material property. The purpose of this paper is to study size of fracture process zone using both theoretical and numerical approaches. Through an approximate theoretical argument, it is shown that width of
96
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
process zone is in general a function of specimen size. A discrete element numerical model is used to verify the usefulness of the theoretical findings. Finally, published experimental data are utilized to support and justify the theoretical and numerical results.
2. Theoretical study Consider two rock specimens I and II in three point bending test (Fig. 1). The applied load in both specimens is assumed to be small enough so that both specimens are in linear elastic regime. The sizes of the specimens are proportional in the x1 x2 plane, i.e. for any point (A) with coordinates x1 and x2 in the smaller specimen, there is a corresponding point (A0 ) in the larger specimen that has the coordinates X1 and X2 such that X i ¼ mxi
ð1Þ
Hence, the size of specimen II is m times the size of specimen I. It is assumed that both specimens have the same thickness (b). The smaller specimen under load P is assumed to be in static equilibrium with the stress sij(x1,x2), strain eij(x1,x2) and displacement ui(x1,x2) fields at a typical point A. We seek for a load P0 applied to specimen II and a strain field (e0 ij) such that the stress (s0 ij) and displacement components (u0 i) at point A0 are related to the corresponding values at point A for specimen I through the following relationships: 0 ij
0 ij ðX 1 ,X 2 Þ ¼
s ¼s
sij ðx1 ,x2 Þ
ð2aÞ
u0i ¼ u0i ðX 1 ,X 2 Þ ¼ mui ðx1 ,x2 Þ
ð2bÞ
The above displacement field satisfies the displacement boundary conditions for specimen II. It is easy to see that the stress components s0 ij satisfy the equilibrium equations (assuming zero body forces): @s0ij @X j
¼
@sij @xk @sij 1 1 @sij ¼ d ¼ ¼0 @xk @X j @xk m kj m @xj
ð3Þ
x2
R
A
w
it is concluded that e0 ij ¼ eij. It is an easy task to show that P0 ¼mP by finding the integral of the vertical traction component along the corresponding semi-circles in Fig. 1 using the fact that sij ¼ s0 ij at the corresponding points along the circles. Above discussion suggests the well known fact that the compliances of beams I and II are identical: C0 ¼ C ¼
d0 P0
¼
md d ¼ mP P
ð5Þ
where d and d0 are the load-point displacements for specimens I and II, respectively. Now consider the scenario that damaged zones have been created around the notch tips in both specimens I and II (Fig. 1). In the remaining of this section we assume that for values of m41 but close to 1, the similitude discussed above still will be approximately valid. An immediate consequence of this assumption is that the width of the damaged zone in specimen II (W0 ) is m times that in specimen I (W0 ¼ mW). On the other hand, using Bazant’s size effect equation [6] and the nominal tensile strength of a beam (sN ¼3PmaxL/2bD2), the relationship between the peak load and the specimen height (D) can be obtained: Pmax ¼
bDBst qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi b 1 þ D=D0
ð6Þ
in which st is the tensile strength of the specimen, b ¼3L/2D, and B and D0 are constants that depend on the material and specimen geometry. To obtain the width of the fracture process zones in Fig. 1, two extreme situations are considered. For the first situation, assume that the specimen sizes are very small. In particular, for very small values of D and D0 ¼mD, using Eq. (6) it is easy to see that plastic behavior of material prevails which leads to no or negligible size effect on strength, i.e. P0max ffim Pmax
ð7Þ
Eq. (8) suggests a linear relationship between the width of process zone at peak load and specimen size, i.e.
D
WðDÞ ¼ lD
ao
x1
L mP
mR
Let the applied force P in Fig. 1a be equal to Pmax/Om. From Eq. (10) it follows that P0 ¼P0 max ¼PmaxOm¼ mP. These loading conditions for beams I and II result in the following approximate equation concerning the width of the process zone:
A' mD
w' mao
ð9Þ
where l is a constant for a given material and specimen shape. The second situation is when D is large. In this case, an approximate ratio of P0 max to Pmax is obtained from Eq. (6): qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mDbBst =b 1 þ mD=D0 pffiffiffiffiffi P0max qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð10Þ ¼ ffi ffi m Pmax DbBst =b 1 þ D=D0
X2
(II )
ð4Þ
Therefore, for small specimens, the developed process zone widths at peak loads in beams I (W) and II (W0 ) are related through the following linear equation: ð8Þ W 0 P0 ¼ P 0max ,D0 ¼ mD ¼ mW ðP ¼ Pmax ,D ¼ DÞ
P
(I)
Considering the equation @u0i @ðmui Þ @xk @u ¼ ¼ i @xk @X j @X j @xj
W 0 9P0
X1
mL Fig. 1. Two geometrically similar beams I and II that are made of the same material.
max
¼ mW9Pmax =pffiffiffi ¼ mW 1 m
ð11Þ
where W1 is the width of the damaged zone in specimen I at a load equal to Pmax/Om. For large specimens, linear fracture mechanics gradually prevails and it is expected that the width of the process zone would change slowly with the specimen size. This fact together with the assumption that m is close to
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
1 suggests that the widths of the process zones for specimens I and II at peak load must be almost identical. Therefore, Eq. (11) leads to W9P ¼ Pmax ¼ mW9P ¼ Pmax =pffiffiffi m
ð12Þ
The solution for Eq. (12) is W ¼ zðD,b, st ÞP2
2
b D2 B2 s2t b ½1þ D=D0 2
ð14Þ
when specimen size approaches infinity, W is expected to approach a constant value. To fulfill this requirement, the function z(D,b,st) is assumed to be proportional to 1/( Ds2t b2) in which the tensile strength in the dominator has been added for convenience. After substitution in Eq. (14), the final relationship between W and D is obtained: W ¼C
DB2 W 1D ¼ b ½1 þ D=D0 D0 ½1 þ D=D0 2
The peak load data from numerical tests and Eq. (18) will be used to study the phenomenon of size effect on strength in fracture testing of three point bending specimens. Furthermore, the approximate validity of Eq. (15), and therefore the soundness of assumptions in deriving this equation will be shown by numerical analysis of rock fracture in the next sections.
ð13Þ
in which z(D,b,st) is in general a function of D and b for geometrically similar structures with a constant a0/D; z is a function of material strength (such as tensile strength) as well. Note that Eq. (13) is valid only for m values close to 1, i.e. P in this equation must be close to the peak load. Eq. (13) implies that for loads close to the peak load, the width of the process zone grows much faster than the increase in the applied load; majority of micro-cracks at the notch tip are developed when the applied load approaches its peak value. This approximate theoretical finding is consistent with the experimental observations [7,8]. The width of the process zone at the peak load can be obtained as a function of the specimen size if Eqs. (6) and (13) are combined: W ¼ zðD,b, st Þ
3. Numerical model A discrete element model was used to simulate the failure process [10,11]. The model consisted of many rigid circular cylinders or disks that interact through normal and shear springs to simulate elasticity. In order to withstand tensile and deviatoric stresses, cylinders or disks were bonded together at contact points. The micro-mechanical constants at a contact point in this model are normal and shear spring constants (Kn,Ks), normal and shear bonds (nb and sb) and friction coefficient (m). To capture softening behavior of a quasi-brittle material, a softening bond feature was implemented in the numerical model. In this model, normal bond at a contact point is assumed to reduce linearly after the peak tensile load (Fig. 2). This introduces a new microscopic constant in the model, the slope in the post-peak region of the normal forcenormal displacement between two cylinders in contact. This slope is indicated by Knp in Fig. 2a. No modification in the relationship of shear force (Fs)-relative shear displacement (us) of a contact is assumed in this simple model. The relationship between Fs and us is shown in Fig. 2b. The loading–unloading paths for both normal and shear contact forces are shown in Fig. 2 with arrows. Softening in shear was not included because this effect is pronounced for loading only under significant (410%qu) mean stress.
ð15Þ
Fn(Tension)
in which C is a dimensionless constant for a given material and geometry. The term WN in Eq. (15) is the width of process zone when the specimen size approaches infinity; it is equal to W1 ¼
CD0 B
b2
97
nb
C D
kn
2
knp
ð16Þ 1
B
As a consequence of Eq. (15), the width of the process zone for specimen II can be predicted if this width is specified for specimen I by using the following equation m½1 þ D=D0 ð17Þ W 0 ¼ W9mD ¼ W9D 1þ mD=D0 It is interesting to note that even though Eq. (15) was derived for large specimens, inserting a small D in this equation results in a linear relationship between width of process zone and specimen height, which is consistent with Eq. (9) that was obtained for small size specimens. Consequently, Eq. (15) shows an approximate relationship between W and D for both large and small specimens even though the prediction of the equation is expected to deviate to some extent from physical observations for intermediate beam sizes. In evaluation of the numerical fracture test data, stress intensity factor at peak load or apparent fracture toughness (KICA) defined as follows will be used: pffiffiffiffi a0 ð18Þ K ICA ¼ sN Df D in which f is a function of the ratio of initial notch length to beam height. The expression for f(a0/D) has been reported in [9]: pffiffiffi a½1:99að1aÞ 2:153:93a þ2:7a2 ð19Þ f ðaÞ ¼ 1:5 ð1 þ 2aÞð1aÞ
1
E
1
Kn
1
F
Un
kn
A Fn(Compression) Fs
Sb
B
μFn
C
A
D
Us
Fig. 2. Relationships in the softening contact bond model: (a) normal force and normal displacement and (b) shear force and shear displacement.
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
The discrete element analysis was performed using a computer program that was developed by Fakhimi [12,13]. CA2 is a hybrid discrete element-finite element code for large deformation analysis of solid materials. The details of sample preparation and calibration procedure have been described in the past [14]. The numerical model was calibrated to simulate the mechanical behavior of Berea sandstone. The grains in this sandstone range from 0.1 to 0.8 mm. The mechanical properties of the rock are the elastic modulus of E¼ 14 GPa, Poisson’s ratio of n ¼0.32, uniaxial compressive strength of qu ¼55–60 MPa, and bending tensile strength (80 240 mm2 beam) of st ¼8.6 MPa [2]. After calibration of the numerical model, uniaxial compressive (40 80 mm2 sample) and bending tensile tests (80 240 mm2 beam) were conducted to verify the accuracy of the model. An elastic modulus of 13.3 GPa, Poisson’s ratio of 0.19, uniaxial compressive strength of 60.5 MPa, and tensile strength of 8.7 MPa were obtained that are in close agreement with the corresponding rock properties. The micro-mechanical properties that were obtained through the calibration procedure are Kn ¼22.0 GPa, Ks ¼5.5 GPa, nb ¼2800 N/m, sb ¼12300 N/m, m ¼ 0.5, s0 (genesis pressure)¼2.2 GPa and Knp ¼1.83 GPa (Kn/Knp ¼12). The circular particles or disks radii (R) were assumed to have a uniform random distribution with a range of 0.27–0.33 mm (Rave ¼0.3 mm).
The peak load, nominal bending tensile strengths, apparent fracture toughness and apparent characteristic size (lchA ¼(KICA/ st)2) for numerical specimens with a0/D ¼0.375 are shown in Table 2. The apparent fracture toughness values in this table were obtained using Eqs. (18) and (19). The variation of apparent fracture toughness with specimen size for a0/D ¼0.375 is shown in Fig. 4a. Notice that as the material becomes more brittle, the value of apparent fracture toughness becomes less size dependent. The apparent characteristic size values in Table 2 and characteristic size (lch ¼(KIC/st)2) can be used to compare brittleness of specimens of the same size; specimens with larger lchA are associated with higher Kn/Knp values and are less brittle. To calculate the size effect parameters B and D0, two different linear regression analyses (methods I and II) have been suggested in [6]. In method I, Eq. (6) is written with respect to the nominal tensile strength (sN) and is linearized (y ¼Axþ C) by defining x¼D and y¼1/s2N. In the II approach, the apparent fracture toughness is written in the following form: pffiffiffiffi Bst K ICA ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Df a0=D 1 þ D=D0
ð20Þ
which can be modified to a linear equation if we choose x¼1/D and y¼(1/KICA)2:
3.1. Numerical tests Three point bending tests were conduced on the numerical samples 20 60, 40 120, 80 240, 160 480, and 320 960 mm2 in size. The applied vertical velocity at the beam center (at top surface) was 2.5 10 10 m/cycle to obtain a quasi-static solution. All specimens had a notch at their center. Three different notch length (a0) to beam height (D) ratios of 0.125, 0.25, and 0.375 were used. The slope of the softening line (Knp) was modified to obtain different quasi-brittle synthetic materials. Five different Kn/Knp values of 0, 10, 20, 50, and 100 were used while other micro-mechanical constants (Kn, Ks, nb, sb, s0, and R) were left unchanged. Uniaxial tensile and compressive tests on samples 40 80 mm2 in dimension resulted in strengths values reported in Table 1. Notice that as Kn/Knp increases, the ratio of the compressive to tensile strength reduces. This is expected because with increase in Kn/Knp, the model becomes less brittle. The numerical bending tests were quite time consuming considering the relative sizes of the beams and the average particle radius. For example, for the 320 960 mm2 beam, more than 1,200,000 particles were needed. The numerical analysis for this beam took more than 4 weeks to finish for Kn/Knp ¼100 on a desk top computer.
y¼
1 B2 st 2 f
2
xþ
1
ð21Þ
2
B2 st 2 f D0
50 Kn/Knp=100
45 40 35 Load (kN)
98
Kn/Knp=50
30 25 20
Kn/Knp=20
15 10
Kn/Knp=10 Kn/Knp=0
5 0 0.00
0.05
0.10 0.15 Displacement (mm)
0.20
0.25
300
3.2. Numerical bending test results
Kn/Knp=100 250 200 Load (kN)
The load-displacement curves for all beam sizes (five different sizes from 20 60 to 320 960 mm2), three different notch to beam height ratios (a0/D¼0.125, 0.25, 0.375), and five different synthetic materials (Kn/Knp ¼0, 10, 20, 50, 100) were obtained. Fig. 3 shows the typical load displacement curves for beams 20 60 and 160 480 mm2 in size and a0/D¼0.375. As expected, by increasing the Kn/Knp value, the material shows a more extensive non-linear hardening behavior as the peak load is gradually approached.
Kn/Knp=50
150 100 Kn/Knp=20
Table 1 Uniaxial tensile and compressive strengths of the numerical specimens.
50
Kn/Knp
st (MPa)
sc (MPa)
sc/st
0 10 20 50 100
4.2 5.6 7.0 8.6 10.7
57.9 63.1 64.8 66.3 73.1
13.8 11.4 9.2 7.7 6.8
Kn/Knp=10 Kn/Knp=0
0 0.0
0.1
0.2 0.3 0.4 Displacement (mm)
0.5
0.6
Fig. 3. Load vs. displacement of loading point for numerical specimens with Kn/Knp ¼ 0, 10, 20, 50, 100 and a0/D¼ 0.375: (a) 20 60 mm2 specimen and (b) 160 480 mm2 specimen.
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
Brittle
Table 2 Peak load, nominal tensile strength, apparent fracture toughness, and apparent characteristic size for numerical specimens with a0/D¼ 0.375.
Kn/Knp=10 Kn/Knp=20
4.0
20 60
40 120
80 240
160 480
320 960
Kn/Knp
Pmax (kN)
sN (MPa)
KICA (MPa m0.5)
lchA (mm)
Kn/Knp=50 3.5
0 10 20 50 100 0 10 20 50 100 0 10 20 50 100 0 10 20 50 100 0 10 20 50 100
10.80 13.44 19.36 30.79 45.90 13.94 21.89 30.33 53.73 79.98 21.36 33.00 47.16 96.51 149.21 34.63 47.56 74.82 153.87 242.69 45.45 63.19 110.63 253.37 387.25
2.43 3.02 4.36 6.93 10.33 1.57 2.46 3.41 6.04 9.00 1.20 1.86 2.65 5.43 8.39 0.97 1.34 2.10 4.33 6.83 0.64 0.89 1.56 3.56 5.45
0.42 0.53 0.76 1.21 1.80 0.39 0.61 0.84 1.49 2.22 0.42 0.65 0.93 1.89 2.93 0.48 0.66 1.04 2.14 3.37 0.45 0.62 1.09 2.49 3.80
10 9 12 20 28 9 12 14 30 43 10 13 18 48 75 13 14 22 62 99 11 12 24 84 126
Kn/Knp=100
3.0 KICA(MP a.m0.5)
Specimen size (mm2)
99
2.5 2.0 1.5 1.0 0.5 0.0 0
40
80
120
160 200 D (mm)
240
280
320
360
Strength criterion 1.0
The data from Table 2 together with fitting lines (Eq. (21)) have been plotted in Fig. 5. The numerical test data illustrate very good linear trends except for Kn/Knp ¼0 (brittle material) that shows greater scatter. The slopes and y-intercepts in Fig. 5 together with st in Table 1 are used to find B and D0 that are reported in Table 3. Table 3 includes the results for a0/D¼0.125 and 0.25 as well. It is important to realize that methods I and II regression analyses do not in general result in the same size effect parameters B and D0 as these regressions imply different weightings of the data points [6]. Notice that the II regression method has produced no unacceptable negative D0 values. Furthermore, the D0 values in method II show consistent increase as the material becomes less brittle which is expected; the brittleness of the material, i.e. D/D0 [6] must decrease with increase in the Kn/Knp value (for a fixed D). On the other hand, the results of method I regression analysis have produced some unacceptable data. The size effect parameters in Table 3 have been used to show how the nominal tensile strength changes with the specimen size (Fig. 4b). Note that the data points for more ductile materials (higher Kn/Knp values) are close to the line representing the strength of material theory. Conversely, more brittle materials follow more closely the LEFM theory. The main goal of this work was to study the width of the fracture process zone around the notch tip. To this end, the process zone in the numerical model is defined with the contact points between circular particles that are in the post-peak regime (e.g. point D in Fig. 2a). These damaged contact points for two specimen sizes 80 240 mm2 and 320 960 mm2 and Kn/Knp of 20 and 50 are shown in Fig. 6. The damaged contacts are shown in green while the actual sharp crack in the middle of the process zone has been shown in red in the figure. Fig. 6 suggests that the width of the process zone is greater for a more ductile material (higher Kn/Knp values). Furthermore, as the size of the specimen increases, the width of the process zone increases as well. Eq. (15) is reformulated to facilitate linear regression analysis of the width of process zone data as follows: 1 1 1 D0 ¼ þ W W1 W1 D
ð22Þ
σN/(Bσt)
LEFM
2 1 Brittle Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting 0.1 0.01
0.1
1 D/D0
10
100
Fig. 4. (a) Apparent fracture toughness vs. specimen size and (b) nominal tensile strength vs. specimen size. The ratio of notch length to specimen height (a0/D) is 0.375.
The variations of 1/W vs. 1/D for a0/D ¼0.375 and different Kn/Knp values are shown in Fig. 7. The linear trends of the data in Fig. 7 confirm the appropriateness of Eq. (15) and support the logic used in the theoretical section of the paper. Using the slopes and y-intercepts of the lines in Fig. 7 and Eq. (22), the D0 values and the widths of the process zones for very large specimens (WN) can be obtained. This information is reported in Table 4. Notice that similar to the situation with regression analysis in methods I and II that did not result in identical D0 values, the D0 values from this III regression method are not the same as those in Table 3. The values of WN and D0 in Table 4 have been used in Eq. (15) to obtain the variation of process zone width with specimen size (Fig. 8a). The good match between the fitting curves and the numerical data implies that Eq. (15) can accurately represent the dependence of the size of process zone on the specimen size. Note that as the material becomes more brittle, the specimen size imposes less impact on the size of the process zone. In Fig. 8b, the width of process zone together with the fitting curves (Eq. (15)) is shown using D0 values from size effect analysis (D0 values from Table 3, method II). Except for Kn/Knp ¼100, the
100
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
7
y = 0.0237x + 4.9795 R2= 0.2855
(1/KICA)2 MPa-2m-1
6
W=3.6 mm W=5 mm
5 y = 0.0244x + 2.2294 R2= 0.8454
4
Brittle Kn/Knp=10 Kn/Knp=20 Linear (Brittle) Linear (Kn/Knp=10) Linear (Kn/Knp=20)
3 y = 0.0187x + 0.8511 R2= 0.9545
2 1
Notch Notch
0 0
10
20
30
40
50
60
W=39 mm
W=22 mm
1/D (m-1) 0.8 y = 0.011x + 0.1447 R2= 0.992
(1/KICA)2 MPa-2m-1
0.7 0.6
Notch
Notch
0.5 y = 0.0052x + 0.0569 R2= 0.9906
0.4
Fig. 6. Width of the process zone around the notch tip for a0/D¼ 0.375: (a) 80 240 mm2 specimen with Kn/Knp ¼ 20, (b) 320 960 mm2 specimen with Kn/Knp ¼ 20, (c) 80 240 mm2 specimen with Kn/Knp ¼ 50 and (d) 320 960 mm2 specimen with Kn/Knp ¼50. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)
0.3 Kn/Knp=50 Kn/Knp=100 Linear (Kn/Knp=50) Linear (Kn/Knp=100)
0.2 0.1 0.0 0
10
20
30
40
50
Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Linear (Kn/Knp=20) Linear (Kn/Knp=50) Linear (Kn/Knp=100)
0.4
60
1/D (m-1) Fig. 5. Linear regression analysis to obtain size effect parameters for beams with a0/D ¼0.375 and (a) Kn/Knp ¼0, 10, 20 and (b) Kn/Knp ¼50, 100.
1/W (mm-1)
Table 3 Size effect parameters for different a0/D values. a0/D
0.125
0.25
0.375
Kn/Knp
0 10 20 50 100 0 10 20 50 100 0 10 20 50 100
Method I
0.3
Method II
B
D0 (mm)
B
D0 (mm)
2.216 1.240 1.096 1.496 1.830 0.779 0.888 1.007 1.208 1.281 1.073 NA 0.729 0.839 1.024
8.9 26.4 56.7 150.0 260.0 39.1 21.6 28.7 92.0 176.7 6.8 0.4 32.0 95.5 103.8
3.283 1.873 2.025 1.960 2.183 2.147 1.181 1.289 1.384 1.463 1.256 0.934 0.846 0.895 1.048
3.8 9.4 12.5 50.2 75.9 3.8 11.4 16.8 51.2 81.5 4.8 10.9 22.0 76.0 91.4
prediction of the fitting equation is satisfactory suggesting that the D0 values from both methods II and III together with Eq. (15) can provide reasonable prediction of the evolution of process zone width as the specimen size is changed.
y = 3.3947x + 0.2181 R2 = 0.8769
0.2 y = 2.1531x + 0.0197 R2 = 0.9936 0.1 y = 1.3986x + 0.0099 R2 = 0.992
0.0 0
0.01
0.02
0.03
0.04
0.05
0.06
-1)
1/D (mm
Fig. 7. Linear regression analysis to obtain WN and D0. The ratio of notch length to specimen height (a0/D) is 0.375.
Table 4 Data from linear regression analysis of Eq. (15) and (22)). a0/D
Kn/Knp
Slope
y-intercept (mm 1)
WN (mm)
D0 (mm)
0.125
20 50 100 20 50 100 20 50 100
1.8318 1.5078 1.009 2.8469 1.8925 1.0949 3.3947 2.1531 1.3986
0.1936 0.0169 0.0062 0.1992 0.0179 0.0101 0.2181 0.0197 0.0099
5.2 59.2 161.3 5.0 55.9 99.0 4.6 50.8 101.0
9.5 89.2 162.7 14.3 105.7 108.4 15.6 109.3 141.3
0.25
0.375
4. Experimental data and discussion of the results Preceding discussions illustrate that the theoretical equation (Eq. (15)) for the size of fracture process zone around a notch tip can provide results that are consistent with those from discrete element analysis. To display further supports for findings of this study, some experimental observations regarding the width of process zone have been compiled in Fig. 9. The work of Otsuka
and Date [4] was concentrated on concrete. They performed tensile tests on different sample sizes. X-ray and three-dimensional Acoustic Emission (AE) techniques were used to investigate the process zone around the notch tip. It should be pointed out that in the work of Otsuka and Date, what is called damaged zone in Fig. 9 is in fact the Fracture Core Zone (FCZ) defined by these authors as a region
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting Kn/Knp=20 Fitting Kn/Knp=50 Fitting Kn/Knp=100
80 70
W (mm)
60
35
Berea sandtone, after Zietlow & Labuz (1998)
50 40 30 20
Charcoal granite, after Zietlow & Labuz (1998) Rockville granite, after Zietlow & Labuz (1998)
25 20 15 10 5
10
0
0 0
40
80
120
160 200 D (mm)
240
280
320
360
80 70
50 40 30 20 10 0 40
80
120
100
200
300
400
500
600
700
800
Fig. 9. Some published information regarding the width of process zone as a function of specimen size.
60
0
0
Specimen size (mm)
Kn/Knp=10 Kn/Knp=20 Kn/Knp=50 Kn/Knp=100 Fitting Kn/Knp=20 Fitting Kn/Knp=50 Fitting Kn/Knp=100
90
W (mm)
Concrete, after Otsuka & Date (2000)
30 Width of damaged zone (mm)
90
101
160 200 D (mm)
240
280
320
360
Fig. 8. Width of process zone vs. specimen size. The fitting curves are based on D0 values from (a) method III regression analysis and (b) method II regression analysis. The ratio of notch length to beam height (a0/D) is 0.375.
around the notch tip and within the process zone that is subjected to more damage and releases a higher percentage of acoustic energy. Zietlow and Labuz [2] studied fracture process zone in different rock beams using AE measurement. Fig. 9 clearly shows the dependence of the width of process zone to specimen size for concrete, Charcoal granite, and Rockville granite. On the other hand, for Berea sandstone, no change in the size of process zone is observed by increasing the specimen size. This indicates that the contact points between grains of sandstone must have a relatively brittle behavior (small Kn/Knp values around 10–20, see Fig. 8). For granite specimens, the sizes of process zones are greater and in these cases, the specimen size has a greater impact on the size of the process zone. In summary, it appears that the numerical, theoretical, and experimental data all confirm that in a material with a larger process zone, the size of the process zone is expected to be more dependent on specimen size (compare the slopes of the lines passed through the data points in Fig. 9). This observation should resolve the controversy in the literature regarding whether the size of process zone is a material property or not. In general, for quasi-brittle materials, process zone like fracture toughness is dependent on size. If the material is highly brittle, the size of fracture process zone shows no noticeable change as the specimen size is modified; fracture process zone appears as a material property within the specimen sizes typically used in laboratory testing. On the contrary, fracture process
zone may not appear as a material property for laboratory size specimens if the material shows a less brittle behavior. From comparison of Figs. 4 and 8a, it appears that the width of process zone ceases changing with specimen size when apparent fracture toughness approaches a constant value (true fracture toughness). Since measurement of fracture process zone is a more expensive and laborious procedure, this observation indicates that measurement of apparent fracture toughness can be utilized as a tool to find out if the size of the specimen is large enough to simulate a mature and fully developed size independent fracture process zone (WN). Furthermore, a useful feature of Eq. (22) is its potential to provide an estimate of WN. By performing tests on small laboratory size specimens and measuring the widths of process zones, WN can be estimated from Eq. (22) without a need to perform tests on large specimens. The numerical results reveal that the width of process zone is not only a function of specimen size but also a function of a0/D and that as the material becomes more ductile, the influence of a0/D on the width of process zone becomes more important (Fig. 10). The numerical results suggest that the size of process zone increases as a0/D reduces. Therefore, for a shorter ligament length, a smaller fracture process zone is expected (Fig. 10). Experimental tests on quasi-brittle materials confirm this observation at least for the length of fracture process zone [1,15]. Nevertheless, for very large specimens, LEFM is expected to prevail and WN should be the same irrespective of the ratio of notch length to beam height. The data in Table 4 are approximately consistent with this expectation except for a0/D ¼0.125 and Kn/Knp ¼100 (see Fig. 10). It appears that for this relatively ductile synthetic material, even for large specimens, LEFM may not completely govern the behavior. This might be suggesting that for this synthetic material, the size of the process zone is much greater than the region of singular stress field dictated by stress intensity factor or fracture toughness leading to invalid asymptotic prediction of LEFM theory. Note that in discrete element analysis of the process zone, random location of the particles can have some influence on the size of the damaged zone. This issue was not investigated in this study and is not expected to affect the general conclusions of the paper as many particles exist within the damaged region that can smooth out the data scatter. Only for very small structures that are made of large grains with respect of the size of process zone, the random location of particles can cause noticeable scatter of the data. This issue needs further investigation. Finally, it is important to realize that even though the main focus of this
102
A. Fakhimi, A. Tarokh / International Journal of Rock Mechanics & Mining Sciences 60 (2013) 95–102
brittle, the dependence of the process zone on specimen size and notch length is reduced. The results of numerical modeling indicate that dependence of process zone on specimen size disappears when the size of the specimen is large enough such that the apparent fracture toughness is equal to the true material toughness. Attention to the peak loads in the bending tests shows that the discrete element analyses were able to accurately mimic the phenomena of size effect in quasi-brittle materials. Consequently, this numerical technique with the softening bond model can be used as a reliable tool to study crack initiation, crack development, and size effect in rock structures that are predominantly subjected to mode I loading. References
Fig. 10. Width of process zone vs. specimen size and the influence of the ratio of notch length to beam height (a0/D). WN values are in mm. WN values that are smaller than 60 mm are for Kn/Knp ¼ 50 and the rest are for Kn/Knp ¼100 (see Table 4).
study was evaluation of the width of fracture process zone, the same idea is applicable for prediction of the length of process zone. Interested readers are referred to [16] for further information.
5. Conclusion Based on Bazant’s size effect law, an approximate theoretical equation for size of process zone in quasi-brittle materials was proposed. It was shown that this equation results in a linear relationship between the inverse of width of process zone and the inverse of specimen size. A discrete element numerical model with tension softening was utilized to verify the appropriateness and accuracy of the theoretical equation. Synthetic materials with different degrees of ductility and sizes were produced and analyzed. The results obtained for the width of the process zone in the numerical simulations show consistency with the theoretical equation. In particular, it was found that the size of process zone is in general a function of specimen size and the ratio of notch length to specimen height. As the material becomes more
[1] Li X, Marasteanu M. The fracture process zone in asphalt mixture at low temperature. Eng Fract Mech 2010;77:1175–90. [2] Zietlow WK, Labuz JF. Measurement of the intrinsic process zone in rock using acoustic emission. Int J Rock Mech Min Sci 1998;35(3):291–9. [3] Zhang D, Wu K. Fracture process zone of notched three-point-bending concrete beams. Cem Concr Res 1999;29:1877–92. [4] Otsuka K, Date H. Fracture process zone in concrete tension specimen. Eng Fract Mech 2000;65:111–31. [5] Bazant ZP, Kazemi MT. Determination of fracture energy, process zone length and brittleness number from size effect, with application to rock and concrete. Int J Fract 1990;44:111–31. [6] Bazant ZP, Planas J. Fracture and size effect in concrete and other quasi-brittle materials. Boca Raton, FL: CRC Press; 1998. [7] Dai ST, Labuz JF. Damage and failure analysis of brittle materials by acoustic emission. J Mater Amer Soc Civ Eng 1997;9(4):200–5. [8] Jankowski LJ, Stys DJ. Formation of the fracture process zone in concrete. Eng Fract Mech 1990;36(2):245–53. [9] Strawley JE. Wide range stress intensity factor expressions for ASTM E399 standard fracture toughness specimens. Int J Fract 1976;12(2):475–6. [10] Potyondy DO, Cundall PA. A bonded-particle model for rock. Int J Rock Mech Min Sci 2004;41:1329–64. [11] Lin Q, Fakhimi A, Haggerty M, Labuz JF. Initiation of tensile and mixed-mode fracture in sandstone. Int J Rock Mech Min Sci 2009;46:489–97. [12] Fakhimi A. Application of slightly overlapped circular particles assembly in numerical simulation of rocks with high friction angle. Eng Geol 2004;74: 129–38. [13] Fakhimi A. A hybrid discrete-finite element model for numerical simulation of geomaterials. Comput Geotech 2009;36:386–95. [14] Fakhimi A, Villegas T. Application of dimensional analysis in calibration of a discrete element model for rock deformation and fracture. Rock Mech Rock Eng 2007;40(2):193–211. [15] Wu Z, Rong H, Zheng J, Xu F, Dong W. An experimental investigation on the FPZ properties in concrete using digital image correlation technique. Eng Fract Mech 2011;78:2978–90. [16] Tarokh A. Discrete element simulation of fracture process zone in quasibrittle materials. MS thesis, Department of Mineral Engineering. Socorro, NM: New Mexico Institute of Mining and Technology, 2012.