- Email: [email protected]
Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks
International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks Xiao-Ping Zhou a,b,c,n, Xin-Bao Gu b, Yun-Teng Wang b a
State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, China School of Civil Engineering, Chongqing University, Chongqing, China c Engineering Institute of Engineering Crops, PLA University of Science and Technology, Nanjing, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 12 August 2014 Received in revised form 24 April 2015 Accepted 8 September 2015
Propagation process of pre-existing flaws in brittle rock materials are simulated using peridynamic theory, in which it is assumed that particles in a continuum interact with each other across a finite distance. Since damage is incorporated at the level of these two particle interactions in peridynamic theory, localization and fracture of materials occur as a natural outgrowth of equations of motion and constitutive models. Numerical simulations of notched semi-circular bend (NSCB) tests and mix mode fracture in a tension–shear rock sample are performed. The effects of array of multiple pre-existing flaws on the propagation process in rock-like materials subjected to uniaxial tensile loads are investigated. The propagation process of macroflaws and microflaws in rock-like materials subjected to uniaxial tensile loads are simulated.The failure modes of rock sample containing the different array of pre-existing flaws are studied. It is found that the present numerical results obtained from peridynamic theory are in good agreement with the previous experimental and numerical results. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Peridynamic theory Propagation process of pre-existing flaws Numerical simulation Brittle rock.
1. Introduction Rock masses contain various types of pre-existing flaws. The propagation process of these pre-existing flaws under compressive and tensile loads are significant in the study of rock engineering. Extensive study has been done on crack propagation in the different rock-like materials under uniaxial compression in 2D experimental studies1–10 and numerical studies11-16. However, the previous experimental and numerical studies mainly focus on the propagation process of the pre-existing flaws in rock-like materials under uniaxial compressive loads, while the experimental and numerical studies on the propagation process of the pre-existing flaws in rock-like materials under tensile loads is very limited. The main reason is that it is very difficult to investigate crack propagation process in rock-like materials under tensile loads in experimental studies. In fact, the failure of rock masses is not only due to compressive loads, but also due to tensile loads. For example, the excavation of tunnels may lead to the occurrence of tensile stress in the surrounding rock masses. The failure of the surrounding rock masses around tunnels occurs when the value of tensile stress is more than that of the uniaxial tensile strength of rock masses. Moreover, since the uniaxial tensile n Corresponding author at: State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, China. Fax: þ86 2365123511 E-mail address: [email protected] (X.-P. Zhou).
http://dx.doi.org/10.1016/j.ijrmms.2015.09.006 1365-1609/& 2015 Elsevier Ltd. All rights reserved.
strength of rock masses is smaller than the uniaxial compressive strength of rock masses, the failure of rock masses subjected to tensile loads is easier than that of rock masses subjected to compressive loads. Therefore, it is important to numerically investigate the propagation process of the pre-existing flaws in rock-like materials under tensile loads. Many numerical methods are put forward to model the propagation process of flaws by researchers. The finite element method (FEM) has been applied to investigate crack growth and coalescence. In the finite element-based method, singular crack-tip elements are frequently encountered.17 Because of the crack tip stress singularity, external fracture criterion must be introduced to determine the crack coalescence and bifurcation, and the problem of the crack nucleation is still not solved.18 In order to overcome the above difficulties, the extended finite element theory19 is proposed to simulate the propagation of cracks. Although many problems of cracks are solved by virtue of the extended finite element method, bifurcation criterion must still be introduced when displacement is discontinuous, and when multiple crack interaction and bifurcation are involved. Besides, a series of difficulties are encountered for the problem of the three-dimensional cracks by the XFEM. In order to solve the problems of interactions among cracks and bifurcation phenomenon of multiple cracks, meshless methods are developed.20 For example, the propagation of cracks can be simulated by Smooth Particle Hydrodynamics (SPH). However, tensile instability problems are encountered in SPH.21 The molecular dynamics can be applied to model the propagation of
242
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
cracks. However, the molecular dynamics has some shortcomings, such as the longer computational time and the lower computational efficiency.22 The phase-field theory can also be applied to model the fracture of the brittle material, in which any external criteria does not need to be adopted.23 In order to overcome the aforementioned shortcomings, peridynamic theory (PD), which is a meshless numerical method based on the nonlocal concept, is introduced to model bifurcation and coalescence of cracks. It is assumed that particles in a continuum interact with each other across a finite distance, its formulation is by integral equations rather than partial differential equations.24 Therefore, the peridynamic method (PD) can be applied to model the problems of continuous or discontinuous displacements.25–27 The peridynamic method (PD) needs not to make use of any external criteria. Moreover, it breaks through the singularity problem of crack tips, it has the virtue of meshless methods and molecular dynamics method, and it avoids the limitations of calculation dimension of molecular dynamics method. Based on the peridynamic method, the bifurcation and coalescence phenomenon of cracks can be taken place spontaneously. Growth velocity and branching angle of the cracks can be modeled correctly, especially for the brittle material.28 Both the peridynamic theory and the phase-field theory can be applied to model the fracture of the brittle material, and they do not need any exterior criterion. However, there are some differences between the peridynamic theory and the phase-field theory. The differences are summarized as follows: (1) The peridynamic method is a meshless particle method, which is based on a nonlocal theory. In peridynamic method, materials are discretized into some particles with mass. While the phase-field theory is a mesh element method, in which materials are discretized into some meshes. (2) For the peridynamic theory, the governing equations are in integral form, whereas for the phase-field theory, the governing equations are of partial differential form. (3) The physical meaning of s is different. In peridynamic method, the parameter s is the bond stretch, whereas in the phase-field theory, the parameter s is an additional continuous field.
where f (u(j ) − u(k ) , x(j ) − x(k ), t ) and f (u(k ) − u(j ) , x(k ) − x(j ), t ) are respectively the pairwise force function between the material point x(k ) and x(j ), α is the coefficient of thermal expansion of the material, ΔT is uniform temperature change; s(k ) (j ) and s(j ) (k ) are respectively the bond stretch between material point x(k ) and x(j ), it is defined by:
x (j ) + u ( j ) − x ( k ) − u ( k ) − x ( j ) − x ( k )
s(k)(j) = s(j)(k) =
x (j ) − x ( k )
(4)
where μ (x(j ) − x(k ), t ) and μ (x(k ) − x(j ), t ) are respectively a history-dependent scalar-valued function between the material point x(k ) and x(j ). The function μ it is given by
⎧ 1 s ( t ′, x − x ) < s (j ) (k ) 0 μ ( x (j ) − x ( k ) , t ) = ⎨ ⎩0 otherwise ⎪
⎪
(5)
where s0 is the critical stretch for bond failure; for plane stress problem, when critical stretch value is obtained, it is defined as
s0 =
4π G 0 9Eδ
(6)
where E is the elastic modulus;δ is a positive number, which is called the horizontal radius; the work G0 required to break all bonds per unit fracture area is then found from
G02
δ
δ
∫0 ∫x ∫0
−1 z ξ
cos
⎛ 1 2 2⎞ 1 ⎜ cξ s0 ⎟dϕdξdz = cs02δ 4 ⎝2 ⎠ 4
(7)
where c(k ) and c(j ) are respectively a “spring constant” of material point x(k ) and x(j ), which are described respectively as follows:
2. Basic theory The peridynamic theory (PD) was proposed by Silling.29 The peridynamic equation of motion is given by29: ..
ρ (x ) u (x , t ) =
πδ 4
(8)
(1)
where ρ (x ) is the mass density, H is a neighborhood of the material point δ /δs = 1.0, b (x, t ) is a prescribed body force density field of the material point x(k ) , which represents the external force per unit reference volume square, t(k ) (j ) (u(j ) − u(k ) , x(j ) − x(k ), t ) and t(j ) (k ) (u(k ) − u(j ) , x(k ) − x(j ), t ) are respectively the force density vector of the material point x(k ) and x(j ). For an isotropic and elastic material, the following formulation can be expressed as29:
)
1 f u (j ) − u (k ) , x (j ) − x (k ) , t 2 1 = μ x (j) − x (k) , t c(k) s(k)(j) − α ΔT 2
t (k)(j) u (j) − u (k) , x (j) − x (k) , t =
(
)
) (
)
1 f u (k ) − u (j ) , x (k ) − x (j ) , t 2 1 = − μ x (k) − x (j) , t c(j) s(j)(k) − α ΔT 2
t (j)(k) u (k) − u (j) , x (k) − x (j) , t =
c (j ) =
9 E (j ) πδ 4
(
) (
)
∫H μ ( x (k) , t , ϵ) dVϵ
(2)
φ ( x (k ) , t ) = 1 −
(3)
where φ (x(k ), t ) is local damage value of any point x(k ) , its ranges are 0 ≤ φ (x(k ), t ) ≤ 1,with 0 representing virgin material, 1 representing the complete disconnection of a point from all of the points with which it initially interacted. Substituting Eqs. (8) into (2), the following expression can be
)
(
(9)
where E(k ) and E(j ) are the elastic moduli of the material points x(k ) and x(j ). In order to model the fracture, a notion of local damage30–32 at a point is introduced, which is defined as
)
(
(
9 E (k )
∫H ( t (k)(j) ( u(j) − u(k), x (j) − x (k), t ) − t (j)(k) ( u(k) − u(j), x (k) − x (j), t ) ) dH + b (x , t )
(
c (k ) =
∫H dVϵ
(10)
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
243
Fig.1. Initiation and propagation of the notch in NSCB tests.
given by
In order to quantitatively describe different influence factors on the crack propagation velocity, the crack propagation velocity at the time ti is computed by33
t(k)(j) ( u(j) − u(k) , x (j) − x (k) , t ) =
9 E (k ) 1 μ ( x (j) − x (k) , t ) 4 ( s(k)(j) − αΔT ) πδ 2
(11)
Combining Eqs. (9) with (3), the following expression can be written:
9 E (j ) 1 μ ( x (k) − x (j) , t ) 4 ( s(j)(k) − αΔT ) πδ 2
(12)
Substituting Eqs. (11) and (12) into Eq. (1), the equation of motion can be rewritten as ..
1 ( μ ( x (k ) H 2
ρ (x ) u(x, t ) = ∫
+ b (x , t )
− x(j ), t )(
‖xi − xi − 1‖ ti − ti − 1
(14)
where xi and xi − 1 denote the crack-tip positions at the current time ti and the previous data-dump time ti − 1, respectively.
3. Numerical method
t(j)(k) ( u(k) − u(j) , x (k) − x (j) , t ) =
νi =
9E(k ) πδ 4
+
9E(j ) πδ 4
)(s(j ) (k ) − αΔT )) dH (13)
The region is discretized into nodes with a known volume and quality in the reference configuration, the nodes form grid, so Eq. (1) can be discretized into the following form: n
ρu" i = ∑p f (upn − uin , xp − xi ) Vp + bin
(15)
where superscript n represents time step, subscript p represents node number, uin represents the displacement of node xi at the nth
244
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
time step; Vp is the volume of the node p; the expression of is written as follows:
⎧ ϵ+η f |ϵ + η| ≠ 0 ⎪ f = ⎨ |ϵ + η| ⎪ 0 |ϵ + η| = 0 ⎩
(16)
In Eq (15), an explicit central difference formula is used for the acceleration as follows: ..
uin =
uin + 1 − 2uin + uin − 1 Δt 2
(17)
where Δt is constant
4. The numerical results 4.1. Numerical simulation of notched semi-circular bend (NSCB) tests Fig. 3. fraction initiation toughness as a function of the loading rate.
The numerical case selected for the PD simulations involves a semi-circle Barre granite specimen with a notch under three-point bending conditions. The geometric configuration of Barre granite specimen is plotted in Fig. 1(a). In this paper, we refer to the preexisting fracture as a flaw, and the growing crack is fracture. The PD simulations are performed by 20000 collocation points with a spacing of Δx=400 μm . The horizon size is specified as δ = 3.015Δx . As in Ref. 34, a material behavior of the numerical model that is consistent with Barre granite is currently implemented in the PD formulation, so a brittle material model is used for Barre granite. P1 and P2 are the forces on the incident bar-sample interface and the transmitted bar-sample interface, respectively. The mechanical parameters of Barre granite are from Ref. 34. Because dynamic analysis is used, the recommended time step 5.33 10 8 s is specified. Based on the selected displacement increment, a total 2000 time steps are performed. The peridynamic theory predicts the crack propagation behavior, as shown in Fig. 1(b)–(e). When the time reaches to 58.63 μs, the fracture first emanates from the tip of the notch. When the time reaches to 71.955 μs, the propagation of the fracture continues. When the time reaches to 106.6 μs, the failure of a semicircle Barre granite specimen with a notch occurs. The growth path of the fracture obtained from PD theory is consistent with experimental results in Fig. 2.,34 The comparison of the fracture toughness between the present numerical results and the suggested method described in Ref. 35 is made in Fig. 3. It can be seen from Fig. 3 that the error of the value of KIC between the numerical results and the suggested method is little. 4.2. Numerical simulation of mix mode failure in a tension–shear
rock sample The geometric configuration of the double-edge-notched rock specimen is illustrated in Fig. 4(a). The dimensions of the specimen are taken to be 200 mm in both length and height. Two notches at edges are 25 mm in length, 5 mm in width. The PD simulations are performed by 2500 collocation points with a spacing of Δx = 4 mm . The horizontal size is specified as δ = 3.015Δx ; Young’s modulus is E¼ 30 GPa, the mass density ρ is 2265 kg/m3, Possion’s ratio v is 1/3, s0 ¼ 5.527 10 4.36 A horizontal shear force of 10 kN is first applied, and then vertical displacement V is applied on the top and bottom of the rock sample by imposing a constant velocity of 0.01 m/s as quasi-static loading. For the numerical calculation, time step is set to dt ¼1 10 7 s, which is less than the critical time step. Initiation and propagation process of mix mode cracks are plotted in Fig. 4(b)–(e). When the time reach to 2000 μs, the initiation of fracture occurs at tips of the left and right notches, as shown in Fig. 4(b). When the time reaches 2035 μs, fracture bifurcation occurs, as shown in Fig. 4(c). When the time reaches 2050 μs, an enclose area is gradually formed between two curvilinear fractures, as plotted in Fig. 4(d). The growth path of fracture obtained from PD theory is consistent with the previous numerical results in Fig. 5.36 At the same time, the comparison of shear deformation at point x ¼ 35 mm, y¼0 between experiments and the numerical results is made in Fig. 6. It can be found from Fig. 6 that the error of the value of shear deformation between the present numerical results and the previous experimental results is little. 4.3. Numerical simulation of propagation process of the multiple echelon flaws
Fig. 2. The experimental results of Barre granite in NSCB test34.
There is a rectangular rock sample with both the length and the width of 0.05 m. The rectangular rock sample is divided into 200 200¼ 40000 discrete nodes. The grid space between adjacent two points is Δx ¼2.5 10 4 m. The calculation parameters are listed as follows: δ = 3Δx , s0 ¼0.04472, Young’s modulus is E¼ 192 GPa, the mass density ρ is 2000 kg/m3, Possion’s ratio v is 1/3, temperature in the entire system keep constant T ¼100 K, temperature coefficient is α = 2.3 × 10−5. For convenience of description and discussion, the pre-existing flaws are marked with ①, ② and ②, as shown in Fig. 7(a). Fig. 7(b)–(g) shows the
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
245
Fig. 4. Numerical simulation of mix mode fracture in a tension–shear rock sample.
propagation process of three echelon flaws with the same length and the same spacing and the same inclination angle of 45°. A uniform time step of 13.3 μs is used and this is a stable time step for the finest model among all tests performed in this paper. The propagation and coalescence processes of three echelon flaws in the rectangular rock sample subjected to the applied velocity 40 m/s are plotted in Fig. 7(b)–(g). When the time reaches 5.985 μs, the tensile fractures first emanate from the left and right tips of the flaws ③ and the right tip of the flaw ② which propagate along the direction vertical to the applied velocity. When the time reaches 6.251 μs, the tensile fracture emanates from the left tip of
the flaw ② which grows along the direction vertical to the applied velocity, the tensile fractures initiating from the left and right tips of the flaw ③ and the right tip of the flaw ② continue to propagate. When the time reaches 6.65 μs, the secondary fractures coalesce, the flaw ③ links with the flaw ②, meanwhile the tensile fracture initiates from the right tip of flaw ①. When the time reaches 7.98 μs, the secondary fractures coalesce, the flaw ① coalesces with the flaw ② at the same time. When the time reaches 10.64 μs, the tensile fractures initiating from the left tip of the flaw ① and the right tip of the flaw ③ respectively propagate to the left and right sides of the rock sample, and the tensile failure
246
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
the flaw ① coalesces with the flaw ②, and the tensile fracture emanates from the right tip of the flaw ③. When the time reaches 6.55 μs, the tensile fracture initiating from the left tip of flaw ③ coalesces with the flaw ②. When the time reaches 7.98 μs, the tensile fractures continue to propagate. When the time reaches 13.3 μs, the tensile fracture emanating from the right tip of the flaw ② links with the flaw ③, the tensile fractures initiating from the left tip of the flaw ① and the right tip of the flaw ③ respectively propagate to the left and right side of the plate, the failure of rock material occurs at the same time. It is observed from Fig. 8 that the flaws ④–⑨ do not initiate and propagate. The main reason is that the interaction among flaws leads to the stress shielding and amplification effects, which lead to bifurcation of flaw growth pattern. Based on the interaction between cracks, Zhou et al.41,42 found the five modes of crack growth, as shown in Fig. 9. This present numerical results are in good agreement with the second mode of crack growth.41,42 4.4. Numerical simulation of propagation and coalescence process of twenty parallel flaws
Fig. 5. The numerical prediction of growth paths of flaws using the VL-coupling scheme.36
Fig. 6. The comparison of shear deformation at point x¼ 35 mm, y ¼0 between experiments and the present numerical results when δ /δs = 1.0 .
of the rock material occurs. It is found from Fig. 7 that the flaw ③ first propagates, the flaw ① finally grows. The main reason is that the stress amplification and stress shielding effects lead to the bifurcation of flaw growth pattern. This numerical results are in good agreement with the theoretical ones.37–41 Fig. 8 shows the propagation process of nine echelon flaws with the same inclination angle of 45° in rectangular rock sample subjected to the applied velocity 40 m/s. The length of flaws and spacing between flaws are plotted in Fig. 8(a). When the time reaches 3.99 μs, the tensile fractures first initiate from the left and right tips of the flaw ② and the left tip of the flaw ③ which propagate along the direction vertical to the applied velocity. When the time reaches 4.655 μs, the secondary fractures coalesce,
There is a rectangular rock sample with the length and width of 0.05 m. The rectangular rock sample contains twenty parallel flaws with the same length, the same spacing and the same inclination angle 0°. The length of flaws and spacing between flaws are depicted in Fig. 10(a). In order to analyze the propagation process of twenty parallel flaws in rectangular rock sample subjected to the applied velocity 150 m/s, the calculation parameters in the geometric configuration 5 and 6 are the same as those in the geometric configurations 3 and 4. Fig. 10(b)–(g) shows the propagation process of doubly periodic rectangular array of twenty flaws. When the time reaches 1.33 μs, the tensile fractures first emanate from the left the flaws ①–④ and the flaws ⑰–⑳, which propagate along the direction vertical to the applied velocity, while the other flaws do not initiate. When the time reaches 2.328 μs, coalescences of the flaws ①–④ and coalescences of the flaws ⑰–⑳ ocur simultaneously, while the other flaws do not initiate. When the time reaches 2.66 μs, the tensile fractures emanate from the left and right tips of the flaws ⑨–⑫, which propagate along the direction vertical to the applied velocity. When the time reaches 3.325 μs, coalescences of the flaws ⑨–⑫ occurs, meanwhile the tensile fractures emanating from the left tips of the flaws ① and ⑰ and the right tips of the flaws ④ and ⑳ continue to propagate. When the time reaches 4.256 μs, the tensile fractures emanating from the left tips of the flaws ① and ⑨ and ⑰, and from the right tips of the flaws ④ and ⑫ and ⑳ respectively propagate to the left and right sides of the rock sample, then the failure of the material occurs. Based on the “group theoretical bifurcation analysis”, Oguni et al.39 found that the alternative pattern of adjacent long and short crack is expected. Based on the interaction among flaws, Chau and Wang40 also found the alternative pattern of adjacent long and short cracks. Zhou et al.41,42 found the five modes of crack growth including the alternative pattern of adjacent long and short cracks. This numerical results are in good agreement with the previous theoretical analysis of the alternative pattern of adjacent long and short cracks.39–42 The length of flaws and spacing between flaws are depicted in Fig. 11(a). Fig. 11(b)–(g) shows propagation process of the diamond-shaped array of twenty flaws with the same inclination angle 0° in the rectangular rock sample subjected to the applied velocity 150 m/s. When the time reaches 1.33 μs, the tensile fractures first emanate from the left and right tips of the flaws ①–④ and the flaws ⑰–⑳, which propagate along the direction vertical to the applied velocity, while the other flaws do not initiate. When
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 7. The propagation process of three echelon flaws with the same length and the same spacing.
247
248
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 8. The propagation process of nine echelon flaws with different length and spacing.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
249
Fig. 9. A typical sketch for the bifurcations in growth pattern corresponding to the first five eigenvalue of homogeneous system of crack increments at the onset of bifurcation.41
the time reaches 2.328 μs, coalescences of the flaws ①–④ and coalescences of the flaws ⑰–⑳ occur simultaneously, while the other flaws do not initiate. When the time reaches 2.66 μs, the tensile fractures emanate from the left and right tips of the flaws ⑨–⑫, which propagate along the direction vertical to the applied velocity. When the time reaches 3.325 μs, coalescences of the flaws ⑨–⑫ occurs, the tensile fractures emanating from the left tips of the flaws ① and ⑰ and the right tips of the flaws ④ and ⑳ continue to propagate, the flaws ⑤ and ⑬ initiate at the same time. When the time reaches 4.256 μs, the tensile fractures emanating from the left tips of the fractures ①, ⑨ and ⑰, and from the right tips of the flaws ④, ⑫ and ⑳ respectively, propagate to the left and right sides of the plate, the tensile fractures emanating from the flaws ⑤ and ⑬ continue to grow, and failure of the rock material occurs. This numerical results are in good agreement with the previous theoretical analysis of the alternative pattern of adjacent long and short cracks.39–42 The difference between the geometric configuration Fig. 10(a) and Fig. 11(a) is that initiation and growth of the flaws ⑤ and ⑬ occur in the rectangular rock sample containing the diamond-shaped array of twenty flaws, while initiation and propagation of the flaws ⑤ and ⑬ do not occur in the rectangular rock sample containing the doubly
periodic rectangular of twenty flaws. It is implied that arrays of the flaws affect the initiation, propagation process of fractures. 4.5. Numerical simulation of the propagation process of the macroflaws and microflaws There is a rectangular rock sample with the length and width of 6 m. The rectangular rock sample contains two parallel macroflaws with the same length of 1.20 m and the same inclination angle of 45° and the seven random microflaws with the different length and different inclination angle, as shown in Fig. 12(a). The length of the microflaw ③ is 0.023 m, and its inclination angle is 60°. The length of the microflaw ④ is 0.029 m, and its inclination angle is 0°. The length of the microflaw ⑤ is 0.032 m, and its inclination angle is 50°. The length of the microflaw ⑥ is 0.047 m, and its inclination angle is 20°. The length of the microflaw ⑦ is 0.028 m, and its inclination angle is 80°. The length of the microflaw ⑧ is 0.02 m, and its inclination angle is 140°. The length of the microflaw ⑨ is 0.025 m, and its inclination angle is 110°. The rectangular rock sample is divided into 200 200 ¼ 40000 discrete nodes. The grid space is Δx ¼3 10 2 m. The calculation parameters in the geometric configuration Fig. 12(a) are the same
250
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 10. The propagation process of doubly periodic rectangular array of 20 flaws.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 11. The propagation process of the diamond-shaped array of 20y flaws.
251
252
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 12. The propagation process of the macro and microflaws.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
as those in the geometric configuration Fig. 7(a). The propagation process of the macroflaws and microflaws in the rectangular rock sample subjected to the applied velocity 40 m/s is plotted in Fig. 12 (b)–(g). When the time reaches 11.7 μs, the tensile fracture first emanates from the left and right tips of the macroflaw ② which propagates along the direction vertical to the applied velocity. When the time reaches 13.3 μs, the tensile fracture emanates from the left and right tips of the macroflaw ① which propagates along the direction vertical to the applied velocity, the tensile fractures emanating from the left and right tips of the macroflaw ② continue to grow at the same time. When the time reaches 15.96 μs, the macroflaw ② coalesces with the microflaw ③, the macroflaw ① links with the microflaw ④, and bifurcations of tensile fractures initiating from the left tip of the macroflaw ① and the right tip of the macroflaw ② occurs at the same time. When the time reaches 18.62 μs, the microflaw ③ coalesces with the microflaw ④, the macroflaw ① links with the microflaw ⑧, the macroflaw ② coalesces with microflaw ⑧. When the time reaches 23.94 μs, the tensile fractures emanating from the tips of the macroflaw ① and the macroflaw ② propagate to the top and bottom of the rock sample, and the failure of the rock sample occurs at the same time. In addition, the other microflaws do not propagate. It is implied that the interaction among macroflaws and microflaws significantly affects the growth paths of cracks.
6. Discussion From the aboved numerical examples, it can be found that for peridynamic numerical model, the initiation and propagation of cracks is autonomous, it does not need to set the propagation path in advance. Growth path, growth velocity, branching, mutual interaction are only predicted by the constitutive model and equations of motion, it need not control these things for any externally supplied relation, these are its advantages over other numerical methods. However, the problem of the failure of materials subjected to compressive loads is rarely studied. Therefore, the peridynamic theory will be developed to model the failure of materials subjected to compressive loads.
7. Conclusions Peridynamic theory is applied to study the propagation process of multiple pre-existing flaws in brittle rock materials. Numerical simulations of notched semi-circular bend (NSCB) tests and mix mode fracture in a tension–shear rock sample are first performed. The present numerical results obtained from peridynamic theory are in good agreement with the previous experimental and numerical results. Then, the effects of array of pre-existing flaws on the propagation process are investigated. The numerical results show alternative pattern of adjacent long and short cracks, which are in good agreement with the previous theoretical results. Finally, the effects of interaction among pre-existing macroflaws and microflaws on the propagation process of macroflaws and microflaws are studied. The numerical results show that the macroflaws first initiate and propagate, and the interaction among macroflaws and microflaws significantly affects the growth paths of cracks.
Acknowledgment This work was supported by Project 973 (Grant no. 2014CB046903), the National Natural Science Foundation of China
253
(Nos. 51325903 and 51279218). Natural Science Foundation Project of CQ CSTC (Nos CSTC, cstc2013kjrcljrccj0001 and cstc2013jcyjys30002).
References 1. Bobet A, Einstein HH. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci. 1998;35:863–888. 2. Bobet A, Einstein HH. Numerical modeling of fracture coalescence in a model rock material. Int J Fract. 1998;92:221–252. 3. Sagong M, Bobet A. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int J Rock Mech Min Sci. 2002;39:229–241. 4. Wong RHC, Chau KT, Tang CA, Lin P. Analysis of crack coalescence in rock-like materials containing three flaws – part I: experimental approach. Int J Rock Mech Min Sci. 2001;38:909–924. 5. Wong LNY, Einstein HH. Crack coalescence in molded gypsum and Carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech Rock Eng. 2009;42:475–511. 6. Wong LNY, Einstein HH. Crack coalescence in molded gypsum and Carrara marble: Part 2. microscopic observations and interpretation. Rock Mech Rock Eng. 2009;42:513–545. 7. Park CH, Bobet A. Crack coalescence in specimens with open and closed flaws: a comparison. Int J Rock Mech Min Sci. 2009;46:819–829. 8. Park CH, Bobet A. Crack initiation, propagation and coalescence from frictional flaws in uniaxial compression. Eng Fract Mech. 2010;77:2727–2748. 9. Lee H, Jeon S. An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int J Solids Struct. 2011;48:979–999. 10. Zhou XP, Cheng H, Feng YF. An experimental study of crack coalescence behaviour in rock-like materials containing multiple flaws under uniaxial compression. Rock Mech Rock Eng. 2014;47:1961–1986. 11. Tang CA, Lin P, Wong RHC, Chau KT. Analysis of crack coalescence in rock-like materials containing three flaws—Part II: numerical approach. Int J Rock Mech Min Sci. 2001;38:925–939. 12. Liu HY, Kou SQ, Lindqvist PA, Tang CA. Numerical simulation of shear fracture (mode II) in heterogeneous brittle rock. Int J Rock Mech Min Sci. 2004;41:14–19. 13. Ning YJ, Yang J, An XM, Ma GW. Modeling rock fracturing and blast-induced rock mass failure via advanced discretization within the discontinuous deformation analysis framework. Comput Geotech. 2011;38:40–49. 14. Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. Int J Rock Mech Min Sci. 2011;48:964–975. 15. Wu ZJ, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech. 2012;39:38–53. 16. Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock Eng. 2015;48:1097–1114. 17. Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng. 1999;46:131–150. 18. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng. 1999;45:601–620. 19. Wagner GJ, Moes N, Liu WK, Belytschko T. The extended finite element method for rigid particles in Stokes flow. Int J Numer Methods Eng. 2001;51:293–313. 20. Belytschko T, Krongauz Y, Organ D, Flemming M, Krysl. P. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng. 1996;139:3–47. 21. Libersky LD, Petsehek AG. Smoothed particle hydrodynamics with strength of materials. Adv Free Lagrange Method Lect Notes Phys. 1990;26:248–257. 22. Izumi S, Katake S. Molecular dynamics study of solid deformation. Trans Jpn Soc Mech Eng A. 1993;59:263–267. 23. Jin Y, Wang Y, Khachaturyan A. Three-dimensional phase field microelasticity theory and modeling of multiple cracks and voids. Appl Phys Lett. 2001;79:3071–3073. 24. Silling SA, Bobaru F. Peridynamic modeling of membranes and fibers. Int J Nonlin Mech. 2005;40:395–409. 25. Silling SA. Dynamic fracture modeling with a meshfree peridynamic code. In: Bathe KJ, editor. Computational fluid and solid mechanics. Amsterdam: Elsevier; 2003. p. 641–644. 26. Silling SA, Askari E. Peridynamic modeling of impact damage. In: Moody FJ, editor. PVP, vol. 489. New York: American Society of Mechanical Engineers; 2004. p. 197–205. 27. Silling SA, Zimmermann M, Abeyaratne R. Deformation of a peridynamic bar. J Mech Phys Solids. 2003;73:173–190. 28. Huang D, Zhang Q, Qiao PZ, Shen F. Peridynamic methods and application. Adv Mech. 2010;4:448–459. 29. Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids. 2000;48:175–209. 30. Gunzburger M, Lehoucq RB. A nonlocal vector calculus with application to nonlocal boundary value problems. Multiscale Model Simul. 2010;8:1581–1589. 31. Ha YD, Bobaru F. Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech. 2011;78:1156–1168. 32. Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int J Numer Methods Eng. 2003;58:1873–1905.
254
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
33. Emmrich E, Weckner O. The peridynamic equation and its spatial discretization. Math Model Anal. 2007;12:1–10. 34. Gao G, Huang S, Xia K, Li Z. Application of digital image correlation (DIC) in dynamic notched semi-circular bend (NSCB) tests. Exp Mech. 2015;55:95–104. 35. Zhou. YX, Xia. K, Li XB, Li HB, Ma GW, Zhao J, Zhou ZL, Dai F. Suggested methods for determining the dynamic strength parameters and mode-I fracture toughness of rock materials. Int J Rock Mech Min Sci. 2012;49:105–112. 36. Liu WF. Discretized bond-based peridynamics for solid mechanics [Ph.D. thesis]. Michigan State University; 2012. 37. Dyskin AV. On the possibility of bifurcation in linear periodic arrays of 2-D cracks. Int. J. Fract.. 1994;67:R31–R36.
38. Muhlhaus HB, Chau KT, Ord A. Bifurcation of crack pattern in arrays of twodimensional cracks. Int. J. Fract.. 1996;77:1–14. 39. Oguni K, Hori M, Ikeda K. Analysis on evolution pattern of periodically distributed defects. Int J Solids Struct. 1997;34:3259–3272. 40. Chau KT, Wang GS. Condition for the onset of bifurcation in en echelon crack arrays. Int J Numer Anal Methods Geomech. 2001;25:289–306. 41. Zhou XP, Xie WT, Qian QH. Bifurcation of collinear crack system under dynamic compression. Theor Appl Fract Mech. 2010;54:166–171. 42. Zhou XP, Qian QH, Yang HQ. Bifurcation condition of crack pattern in the periodic rectangular array. Theor Appl Fract Mech. 2008;49:206–212.
Contents lists available at ScienceDirect
International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms
Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks Xiao-Ping Zhou a,b,c,n, Xin-Bao Gu b, Yun-Teng Wang b a
State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, China School of Civil Engineering, Chongqing University, Chongqing, China c Engineering Institute of Engineering Crops, PLA University of Science and Technology, Nanjing, China b
art ic l e i nf o
a b s t r a c t
Article history: Received 12 August 2014 Received in revised form 24 April 2015 Accepted 8 September 2015
Propagation process of pre-existing flaws in brittle rock materials are simulated using peridynamic theory, in which it is assumed that particles in a continuum interact with each other across a finite distance. Since damage is incorporated at the level of these two particle interactions in peridynamic theory, localization and fracture of materials occur as a natural outgrowth of equations of motion and constitutive models. Numerical simulations of notched semi-circular bend (NSCB) tests and mix mode fracture in a tension–shear rock sample are performed. The effects of array of multiple pre-existing flaws on the propagation process in rock-like materials subjected to uniaxial tensile loads are investigated. The propagation process of macroflaws and microflaws in rock-like materials subjected to uniaxial tensile loads are simulated.The failure modes of rock sample containing the different array of pre-existing flaws are studied. It is found that the present numerical results obtained from peridynamic theory are in good agreement with the previous experimental and numerical results. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Peridynamic theory Propagation process of pre-existing flaws Numerical simulation Brittle rock.
1. Introduction Rock masses contain various types of pre-existing flaws. The propagation process of these pre-existing flaws under compressive and tensile loads are significant in the study of rock engineering. Extensive study has been done on crack propagation in the different rock-like materials under uniaxial compression in 2D experimental studies1–10 and numerical studies11-16. However, the previous experimental and numerical studies mainly focus on the propagation process of the pre-existing flaws in rock-like materials under uniaxial compressive loads, while the experimental and numerical studies on the propagation process of the pre-existing flaws in rock-like materials under tensile loads is very limited. The main reason is that it is very difficult to investigate crack propagation process in rock-like materials under tensile loads in experimental studies. In fact, the failure of rock masses is not only due to compressive loads, but also due to tensile loads. For example, the excavation of tunnels may lead to the occurrence of tensile stress in the surrounding rock masses. The failure of the surrounding rock masses around tunnels occurs when the value of tensile stress is more than that of the uniaxial tensile strength of rock masses. Moreover, since the uniaxial tensile n Corresponding author at: State Key Laboratory of Coal Mine Disaster Dynamics and Control, Chongqing University, Chongqing, China. Fax: þ86 2365123511 E-mail address: [email protected] (X.-P. Zhou).
http://dx.doi.org/10.1016/j.ijrmms.2015.09.006 1365-1609/& 2015 Elsevier Ltd. All rights reserved.
strength of rock masses is smaller than the uniaxial compressive strength of rock masses, the failure of rock masses subjected to tensile loads is easier than that of rock masses subjected to compressive loads. Therefore, it is important to numerically investigate the propagation process of the pre-existing flaws in rock-like materials under tensile loads. Many numerical methods are put forward to model the propagation process of flaws by researchers. The finite element method (FEM) has been applied to investigate crack growth and coalescence. In the finite element-based method, singular crack-tip elements are frequently encountered.17 Because of the crack tip stress singularity, external fracture criterion must be introduced to determine the crack coalescence and bifurcation, and the problem of the crack nucleation is still not solved.18 In order to overcome the above difficulties, the extended finite element theory19 is proposed to simulate the propagation of cracks. Although many problems of cracks are solved by virtue of the extended finite element method, bifurcation criterion must still be introduced when displacement is discontinuous, and when multiple crack interaction and bifurcation are involved. Besides, a series of difficulties are encountered for the problem of the three-dimensional cracks by the XFEM. In order to solve the problems of interactions among cracks and bifurcation phenomenon of multiple cracks, meshless methods are developed.20 For example, the propagation of cracks can be simulated by Smooth Particle Hydrodynamics (SPH). However, tensile instability problems are encountered in SPH.21 The molecular dynamics can be applied to model the propagation of
242
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
cracks. However, the molecular dynamics has some shortcomings, such as the longer computational time and the lower computational efficiency.22 The phase-field theory can also be applied to model the fracture of the brittle material, in which any external criteria does not need to be adopted.23 In order to overcome the aforementioned shortcomings, peridynamic theory (PD), which is a meshless numerical method based on the nonlocal concept, is introduced to model bifurcation and coalescence of cracks. It is assumed that particles in a continuum interact with each other across a finite distance, its formulation is by integral equations rather than partial differential equations.24 Therefore, the peridynamic method (PD) can be applied to model the problems of continuous or discontinuous displacements.25–27 The peridynamic method (PD) needs not to make use of any external criteria. Moreover, it breaks through the singularity problem of crack tips, it has the virtue of meshless methods and molecular dynamics method, and it avoids the limitations of calculation dimension of molecular dynamics method. Based on the peridynamic method, the bifurcation and coalescence phenomenon of cracks can be taken place spontaneously. Growth velocity and branching angle of the cracks can be modeled correctly, especially for the brittle material.28 Both the peridynamic theory and the phase-field theory can be applied to model the fracture of the brittle material, and they do not need any exterior criterion. However, there are some differences between the peridynamic theory and the phase-field theory. The differences are summarized as follows: (1) The peridynamic method is a meshless particle method, which is based on a nonlocal theory. In peridynamic method, materials are discretized into some particles with mass. While the phase-field theory is a mesh element method, in which materials are discretized into some meshes. (2) For the peridynamic theory, the governing equations are in integral form, whereas for the phase-field theory, the governing equations are of partial differential form. (3) The physical meaning of s is different. In peridynamic method, the parameter s is the bond stretch, whereas in the phase-field theory, the parameter s is an additional continuous field.
where f (u(j ) − u(k ) , x(j ) − x(k ), t ) and f (u(k ) − u(j ) , x(k ) − x(j ), t ) are respectively the pairwise force function between the material point x(k ) and x(j ), α is the coefficient of thermal expansion of the material, ΔT is uniform temperature change; s(k ) (j ) and s(j ) (k ) are respectively the bond stretch between material point x(k ) and x(j ), it is defined by:
x (j ) + u ( j ) − x ( k ) − u ( k ) − x ( j ) − x ( k )
s(k)(j) = s(j)(k) =
x (j ) − x ( k )
(4)
where μ (x(j ) − x(k ), t ) and μ (x(k ) − x(j ), t ) are respectively a history-dependent scalar-valued function between the material point x(k ) and x(j ). The function μ it is given by
⎧ 1 s ( t ′, x − x ) < s (j ) (k ) 0 μ ( x (j ) − x ( k ) , t ) = ⎨ ⎩0 otherwise ⎪
⎪
(5)
where s0 is the critical stretch for bond failure; for plane stress problem, when critical stretch value is obtained, it is defined as
s0 =
4π G 0 9Eδ
(6)
where E is the elastic modulus;δ is a positive number, which is called the horizontal radius; the work G0 required to break all bonds per unit fracture area is then found from
G02
δ
δ
∫0 ∫x ∫0
−1 z ξ
cos
⎛ 1 2 2⎞ 1 ⎜ cξ s0 ⎟dϕdξdz = cs02δ 4 ⎝2 ⎠ 4
(7)
where c(k ) and c(j ) are respectively a “spring constant” of material point x(k ) and x(j ), which are described respectively as follows:
2. Basic theory The peridynamic theory (PD) was proposed by Silling.29 The peridynamic equation of motion is given by29: ..
ρ (x ) u (x , t ) =
πδ 4
(8)
(1)
where ρ (x ) is the mass density, H is a neighborhood of the material point δ /δs = 1.0, b (x, t ) is a prescribed body force density field of the material point x(k ) , which represents the external force per unit reference volume square, t(k ) (j ) (u(j ) − u(k ) , x(j ) − x(k ), t ) and t(j ) (k ) (u(k ) − u(j ) , x(k ) − x(j ), t ) are respectively the force density vector of the material point x(k ) and x(j ). For an isotropic and elastic material, the following formulation can be expressed as29:
)
1 f u (j ) − u (k ) , x (j ) − x (k ) , t 2 1 = μ x (j) − x (k) , t c(k) s(k)(j) − α ΔT 2
t (k)(j) u (j) − u (k) , x (j) − x (k) , t =
(
)
) (
)
1 f u (k ) − u (j ) , x (k ) − x (j ) , t 2 1 = − μ x (k) − x (j) , t c(j) s(j)(k) − α ΔT 2
t (j)(k) u (k) − u (j) , x (k) − x (j) , t =
c (j ) =
9 E (j ) πδ 4
(
) (
)
∫H μ ( x (k) , t , ϵ) dVϵ
(2)
φ ( x (k ) , t ) = 1 −
(3)
where φ (x(k ), t ) is local damage value of any point x(k ) , its ranges are 0 ≤ φ (x(k ), t ) ≤ 1,with 0 representing virgin material, 1 representing the complete disconnection of a point from all of the points with which it initially interacted. Substituting Eqs. (8) into (2), the following expression can be
)
(
(9)
where E(k ) and E(j ) are the elastic moduli of the material points x(k ) and x(j ). In order to model the fracture, a notion of local damage30–32 at a point is introduced, which is defined as
)
(
(
9 E (k )
∫H ( t (k)(j) ( u(j) − u(k), x (j) − x (k), t ) − t (j)(k) ( u(k) − u(j), x (k) − x (j), t ) ) dH + b (x , t )
(
c (k ) =
∫H dVϵ
(10)
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
243
Fig.1. Initiation and propagation of the notch in NSCB tests.
given by
In order to quantitatively describe different influence factors on the crack propagation velocity, the crack propagation velocity at the time ti is computed by33
t(k)(j) ( u(j) − u(k) , x (j) − x (k) , t ) =
9 E (k ) 1 μ ( x (j) − x (k) , t ) 4 ( s(k)(j) − αΔT ) πδ 2
(11)
Combining Eqs. (9) with (3), the following expression can be written:
9 E (j ) 1 μ ( x (k) − x (j) , t ) 4 ( s(j)(k) − αΔT ) πδ 2
(12)
Substituting Eqs. (11) and (12) into Eq. (1), the equation of motion can be rewritten as ..
1 ( μ ( x (k ) H 2
ρ (x ) u(x, t ) = ∫
+ b (x , t )
− x(j ), t )(
‖xi − xi − 1‖ ti − ti − 1
(14)
where xi and xi − 1 denote the crack-tip positions at the current time ti and the previous data-dump time ti − 1, respectively.
3. Numerical method
t(j)(k) ( u(k) − u(j) , x (k) − x (j) , t ) =
νi =
9E(k ) πδ 4
+
9E(j ) πδ 4
)(s(j ) (k ) − αΔT )) dH (13)
The region is discretized into nodes with a known volume and quality in the reference configuration, the nodes form grid, so Eq. (1) can be discretized into the following form: n
ρu" i = ∑p f (upn − uin , xp − xi ) Vp + bin
(15)
where superscript n represents time step, subscript p represents node number, uin represents the displacement of node xi at the nth
244
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
time step; Vp is the volume of the node p; the expression of is written as follows:
⎧ ϵ+η f |ϵ + η| ≠ 0 ⎪ f = ⎨ |ϵ + η| ⎪ 0 |ϵ + η| = 0 ⎩
(16)
In Eq (15), an explicit central difference formula is used for the acceleration as follows: ..
uin =
uin + 1 − 2uin + uin − 1 Δt 2
(17)
where Δt is constant
4. The numerical results 4.1. Numerical simulation of notched semi-circular bend (NSCB) tests Fig. 3. fraction initiation toughness as a function of the loading rate.
The numerical case selected for the PD simulations involves a semi-circle Barre granite specimen with a notch under three-point bending conditions. The geometric configuration of Barre granite specimen is plotted in Fig. 1(a). In this paper, we refer to the preexisting fracture as a flaw, and the growing crack is fracture. The PD simulations are performed by 20000 collocation points with a spacing of Δx=400 μm . The horizon size is specified as δ = 3.015Δx . As in Ref. 34, a material behavior of the numerical model that is consistent with Barre granite is currently implemented in the PD formulation, so a brittle material model is used for Barre granite. P1 and P2 are the forces on the incident bar-sample interface and the transmitted bar-sample interface, respectively. The mechanical parameters of Barre granite are from Ref. 34. Because dynamic analysis is used, the recommended time step 5.33 10 8 s is specified. Based on the selected displacement increment, a total 2000 time steps are performed. The peridynamic theory predicts the crack propagation behavior, as shown in Fig. 1(b)–(e). When the time reaches to 58.63 μs, the fracture first emanates from the tip of the notch. When the time reaches to 71.955 μs, the propagation of the fracture continues. When the time reaches to 106.6 μs, the failure of a semicircle Barre granite specimen with a notch occurs. The growth path of the fracture obtained from PD theory is consistent with experimental results in Fig. 2.,34 The comparison of the fracture toughness between the present numerical results and the suggested method described in Ref. 35 is made in Fig. 3. It can be seen from Fig. 3 that the error of the value of KIC between the numerical results and the suggested method is little. 4.2. Numerical simulation of mix mode failure in a tension–shear
rock sample The geometric configuration of the double-edge-notched rock specimen is illustrated in Fig. 4(a). The dimensions of the specimen are taken to be 200 mm in both length and height. Two notches at edges are 25 mm in length, 5 mm in width. The PD simulations are performed by 2500 collocation points with a spacing of Δx = 4 mm . The horizontal size is specified as δ = 3.015Δx ; Young’s modulus is E¼ 30 GPa, the mass density ρ is 2265 kg/m3, Possion’s ratio v is 1/3, s0 ¼ 5.527 10 4.36 A horizontal shear force of 10 kN is first applied, and then vertical displacement V is applied on the top and bottom of the rock sample by imposing a constant velocity of 0.01 m/s as quasi-static loading. For the numerical calculation, time step is set to dt ¼1 10 7 s, which is less than the critical time step. Initiation and propagation process of mix mode cracks are plotted in Fig. 4(b)–(e). When the time reach to 2000 μs, the initiation of fracture occurs at tips of the left and right notches, as shown in Fig. 4(b). When the time reaches 2035 μs, fracture bifurcation occurs, as shown in Fig. 4(c). When the time reaches 2050 μs, an enclose area is gradually formed between two curvilinear fractures, as plotted in Fig. 4(d). The growth path of fracture obtained from PD theory is consistent with the previous numerical results in Fig. 5.36 At the same time, the comparison of shear deformation at point x ¼ 35 mm, y¼0 between experiments and the numerical results is made in Fig. 6. It can be found from Fig. 6 that the error of the value of shear deformation between the present numerical results and the previous experimental results is little. 4.3. Numerical simulation of propagation process of the multiple echelon flaws
Fig. 2. The experimental results of Barre granite in NSCB test34.
There is a rectangular rock sample with both the length and the width of 0.05 m. The rectangular rock sample is divided into 200 200¼ 40000 discrete nodes. The grid space between adjacent two points is Δx ¼2.5 10 4 m. The calculation parameters are listed as follows: δ = 3Δx , s0 ¼0.04472, Young’s modulus is E¼ 192 GPa, the mass density ρ is 2000 kg/m3, Possion’s ratio v is 1/3, temperature in the entire system keep constant T ¼100 K, temperature coefficient is α = 2.3 × 10−5. For convenience of description and discussion, the pre-existing flaws are marked with ①, ② and ②, as shown in Fig. 7(a). Fig. 7(b)–(g) shows the
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
245
Fig. 4. Numerical simulation of mix mode fracture in a tension–shear rock sample.
propagation process of three echelon flaws with the same length and the same spacing and the same inclination angle of 45°. A uniform time step of 13.3 μs is used and this is a stable time step for the finest model among all tests performed in this paper. The propagation and coalescence processes of three echelon flaws in the rectangular rock sample subjected to the applied velocity 40 m/s are plotted in Fig. 7(b)–(g). When the time reaches 5.985 μs, the tensile fractures first emanate from the left and right tips of the flaws ③ and the right tip of the flaw ② which propagate along the direction vertical to the applied velocity. When the time reaches 6.251 μs, the tensile fracture emanates from the left tip of
the flaw ② which grows along the direction vertical to the applied velocity, the tensile fractures initiating from the left and right tips of the flaw ③ and the right tip of the flaw ② continue to propagate. When the time reaches 6.65 μs, the secondary fractures coalesce, the flaw ③ links with the flaw ②, meanwhile the tensile fracture initiates from the right tip of flaw ①. When the time reaches 7.98 μs, the secondary fractures coalesce, the flaw ① coalesces with the flaw ② at the same time. When the time reaches 10.64 μs, the tensile fractures initiating from the left tip of the flaw ① and the right tip of the flaw ③ respectively propagate to the left and right sides of the rock sample, and the tensile failure
246
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
the flaw ① coalesces with the flaw ②, and the tensile fracture emanates from the right tip of the flaw ③. When the time reaches 6.55 μs, the tensile fracture initiating from the left tip of flaw ③ coalesces with the flaw ②. When the time reaches 7.98 μs, the tensile fractures continue to propagate. When the time reaches 13.3 μs, the tensile fracture emanating from the right tip of the flaw ② links with the flaw ③, the tensile fractures initiating from the left tip of the flaw ① and the right tip of the flaw ③ respectively propagate to the left and right side of the plate, the failure of rock material occurs at the same time. It is observed from Fig. 8 that the flaws ④–⑨ do not initiate and propagate. The main reason is that the interaction among flaws leads to the stress shielding and amplification effects, which lead to bifurcation of flaw growth pattern. Based on the interaction between cracks, Zhou et al.41,42 found the five modes of crack growth, as shown in Fig. 9. This present numerical results are in good agreement with the second mode of crack growth.41,42 4.4. Numerical simulation of propagation and coalescence process of twenty parallel flaws
Fig. 5. The numerical prediction of growth paths of flaws using the VL-coupling scheme.36
Fig. 6. The comparison of shear deformation at point x¼ 35 mm, y ¼0 between experiments and the present numerical results when δ /δs = 1.0 .
of the rock material occurs. It is found from Fig. 7 that the flaw ③ first propagates, the flaw ① finally grows. The main reason is that the stress amplification and stress shielding effects lead to the bifurcation of flaw growth pattern. This numerical results are in good agreement with the theoretical ones.37–41 Fig. 8 shows the propagation process of nine echelon flaws with the same inclination angle of 45° in rectangular rock sample subjected to the applied velocity 40 m/s. The length of flaws and spacing between flaws are plotted in Fig. 8(a). When the time reaches 3.99 μs, the tensile fractures first initiate from the left and right tips of the flaw ② and the left tip of the flaw ③ which propagate along the direction vertical to the applied velocity. When the time reaches 4.655 μs, the secondary fractures coalesce,
There is a rectangular rock sample with the length and width of 0.05 m. The rectangular rock sample contains twenty parallel flaws with the same length, the same spacing and the same inclination angle 0°. The length of flaws and spacing between flaws are depicted in Fig. 10(a). In order to analyze the propagation process of twenty parallel flaws in rectangular rock sample subjected to the applied velocity 150 m/s, the calculation parameters in the geometric configuration 5 and 6 are the same as those in the geometric configurations 3 and 4. Fig. 10(b)–(g) shows the propagation process of doubly periodic rectangular array of twenty flaws. When the time reaches 1.33 μs, the tensile fractures first emanate from the left the flaws ①–④ and the flaws ⑰–⑳, which propagate along the direction vertical to the applied velocity, while the other flaws do not initiate. When the time reaches 2.328 μs, coalescences of the flaws ①–④ and coalescences of the flaws ⑰–⑳ ocur simultaneously, while the other flaws do not initiate. When the time reaches 2.66 μs, the tensile fractures emanate from the left and right tips of the flaws ⑨–⑫, which propagate along the direction vertical to the applied velocity. When the time reaches 3.325 μs, coalescences of the flaws ⑨–⑫ occurs, meanwhile the tensile fractures emanating from the left tips of the flaws ① and ⑰ and the right tips of the flaws ④ and ⑳ continue to propagate. When the time reaches 4.256 μs, the tensile fractures emanating from the left tips of the flaws ① and ⑨ and ⑰, and from the right tips of the flaws ④ and ⑫ and ⑳ respectively propagate to the left and right sides of the rock sample, then the failure of the material occurs. Based on the “group theoretical bifurcation analysis”, Oguni et al.39 found that the alternative pattern of adjacent long and short crack is expected. Based on the interaction among flaws, Chau and Wang40 also found the alternative pattern of adjacent long and short cracks. Zhou et al.41,42 found the five modes of crack growth including the alternative pattern of adjacent long and short cracks. This numerical results are in good agreement with the previous theoretical analysis of the alternative pattern of adjacent long and short cracks.39–42 The length of flaws and spacing between flaws are depicted in Fig. 11(a). Fig. 11(b)–(g) shows propagation process of the diamond-shaped array of twenty flaws with the same inclination angle 0° in the rectangular rock sample subjected to the applied velocity 150 m/s. When the time reaches 1.33 μs, the tensile fractures first emanate from the left and right tips of the flaws ①–④ and the flaws ⑰–⑳, which propagate along the direction vertical to the applied velocity, while the other flaws do not initiate. When
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 7. The propagation process of three echelon flaws with the same length and the same spacing.
247
248
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 8. The propagation process of nine echelon flaws with different length and spacing.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
249
Fig. 9. A typical sketch for the bifurcations in growth pattern corresponding to the first five eigenvalue of homogeneous system of crack increments at the onset of bifurcation.41
the time reaches 2.328 μs, coalescences of the flaws ①–④ and coalescences of the flaws ⑰–⑳ occur simultaneously, while the other flaws do not initiate. When the time reaches 2.66 μs, the tensile fractures emanate from the left and right tips of the flaws ⑨–⑫, which propagate along the direction vertical to the applied velocity. When the time reaches 3.325 μs, coalescences of the flaws ⑨–⑫ occurs, the tensile fractures emanating from the left tips of the flaws ① and ⑰ and the right tips of the flaws ④ and ⑳ continue to propagate, the flaws ⑤ and ⑬ initiate at the same time. When the time reaches 4.256 μs, the tensile fractures emanating from the left tips of the fractures ①, ⑨ and ⑰, and from the right tips of the flaws ④, ⑫ and ⑳ respectively, propagate to the left and right sides of the plate, the tensile fractures emanating from the flaws ⑤ and ⑬ continue to grow, and failure of the rock material occurs. This numerical results are in good agreement with the previous theoretical analysis of the alternative pattern of adjacent long and short cracks.39–42 The difference between the geometric configuration Fig. 10(a) and Fig. 11(a) is that initiation and growth of the flaws ⑤ and ⑬ occur in the rectangular rock sample containing the diamond-shaped array of twenty flaws, while initiation and propagation of the flaws ⑤ and ⑬ do not occur in the rectangular rock sample containing the doubly
periodic rectangular of twenty flaws. It is implied that arrays of the flaws affect the initiation, propagation process of fractures. 4.5. Numerical simulation of the propagation process of the macroflaws and microflaws There is a rectangular rock sample with the length and width of 6 m. The rectangular rock sample contains two parallel macroflaws with the same length of 1.20 m and the same inclination angle of 45° and the seven random microflaws with the different length and different inclination angle, as shown in Fig. 12(a). The length of the microflaw ③ is 0.023 m, and its inclination angle is 60°. The length of the microflaw ④ is 0.029 m, and its inclination angle is 0°. The length of the microflaw ⑤ is 0.032 m, and its inclination angle is 50°. The length of the microflaw ⑥ is 0.047 m, and its inclination angle is 20°. The length of the microflaw ⑦ is 0.028 m, and its inclination angle is 80°. The length of the microflaw ⑧ is 0.02 m, and its inclination angle is 140°. The length of the microflaw ⑨ is 0.025 m, and its inclination angle is 110°. The rectangular rock sample is divided into 200 200 ¼ 40000 discrete nodes. The grid space is Δx ¼3 10 2 m. The calculation parameters in the geometric configuration Fig. 12(a) are the same
250
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 10. The propagation process of doubly periodic rectangular array of 20 flaws.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 11. The propagation process of the diamond-shaped array of 20y flaws.
251
252
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
Fig. 12. The propagation process of the macro and microflaws.
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
as those in the geometric configuration Fig. 7(a). The propagation process of the macroflaws and microflaws in the rectangular rock sample subjected to the applied velocity 40 m/s is plotted in Fig. 12 (b)–(g). When the time reaches 11.7 μs, the tensile fracture first emanates from the left and right tips of the macroflaw ② which propagates along the direction vertical to the applied velocity. When the time reaches 13.3 μs, the tensile fracture emanates from the left and right tips of the macroflaw ① which propagates along the direction vertical to the applied velocity, the tensile fractures emanating from the left and right tips of the macroflaw ② continue to grow at the same time. When the time reaches 15.96 μs, the macroflaw ② coalesces with the microflaw ③, the macroflaw ① links with the microflaw ④, and bifurcations of tensile fractures initiating from the left tip of the macroflaw ① and the right tip of the macroflaw ② occurs at the same time. When the time reaches 18.62 μs, the microflaw ③ coalesces with the microflaw ④, the macroflaw ① links with the microflaw ⑧, the macroflaw ② coalesces with microflaw ⑧. When the time reaches 23.94 μs, the tensile fractures emanating from the tips of the macroflaw ① and the macroflaw ② propagate to the top and bottom of the rock sample, and the failure of the rock sample occurs at the same time. In addition, the other microflaws do not propagate. It is implied that the interaction among macroflaws and microflaws significantly affects the growth paths of cracks.
6. Discussion From the aboved numerical examples, it can be found that for peridynamic numerical model, the initiation and propagation of cracks is autonomous, it does not need to set the propagation path in advance. Growth path, growth velocity, branching, mutual interaction are only predicted by the constitutive model and equations of motion, it need not control these things for any externally supplied relation, these are its advantages over other numerical methods. However, the problem of the failure of materials subjected to compressive loads is rarely studied. Therefore, the peridynamic theory will be developed to model the failure of materials subjected to compressive loads.
7. Conclusions Peridynamic theory is applied to study the propagation process of multiple pre-existing flaws in brittle rock materials. Numerical simulations of notched semi-circular bend (NSCB) tests and mix mode fracture in a tension–shear rock sample are first performed. The present numerical results obtained from peridynamic theory are in good agreement with the previous experimental and numerical results. Then, the effects of array of pre-existing flaws on the propagation process are investigated. The numerical results show alternative pattern of adjacent long and short cracks, which are in good agreement with the previous theoretical results. Finally, the effects of interaction among pre-existing macroflaws and microflaws on the propagation process of macroflaws and microflaws are studied. The numerical results show that the macroflaws first initiate and propagate, and the interaction among macroflaws and microflaws significantly affects the growth paths of cracks.
Acknowledgment This work was supported by Project 973 (Grant no. 2014CB046903), the National Natural Science Foundation of China
253
(Nos. 51325903 and 51279218). Natural Science Foundation Project of CQ CSTC (Nos CSTC, cstc2013kjrcljrccj0001 and cstc2013jcyjys30002).
References 1. Bobet A, Einstein HH. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int J Rock Mech Min Sci. 1998;35:863–888. 2. Bobet A, Einstein HH. Numerical modeling of fracture coalescence in a model rock material. Int J Fract. 1998;92:221–252. 3. Sagong M, Bobet A. Coalescence of multiple flaws in a rock-model material in uniaxial compression. Int J Rock Mech Min Sci. 2002;39:229–241. 4. Wong RHC, Chau KT, Tang CA, Lin P. Analysis of crack coalescence in rock-like materials containing three flaws – part I: experimental approach. Int J Rock Mech Min Sci. 2001;38:909–924. 5. Wong LNY, Einstein HH. Crack coalescence in molded gypsum and Carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech Rock Eng. 2009;42:475–511. 6. Wong LNY, Einstein HH. Crack coalescence in molded gypsum and Carrara marble: Part 2. microscopic observations and interpretation. Rock Mech Rock Eng. 2009;42:513–545. 7. Park CH, Bobet A. Crack coalescence in specimens with open and closed flaws: a comparison. Int J Rock Mech Min Sci. 2009;46:819–829. 8. Park CH, Bobet A. Crack initiation, propagation and coalescence from frictional flaws in uniaxial compression. Eng Fract Mech. 2010;77:2727–2748. 9. Lee H, Jeon S. An experimental and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression. Int J Solids Struct. 2011;48:979–999. 10. Zhou XP, Cheng H, Feng YF. An experimental study of crack coalescence behaviour in rock-like materials containing multiple flaws under uniaxial compression. Rock Mech Rock Eng. 2014;47:1961–1986. 11. Tang CA, Lin P, Wong RHC, Chau KT. Analysis of crack coalescence in rock-like materials containing three flaws—Part II: numerical approach. Int J Rock Mech Min Sci. 2001;38:925–939. 12. Liu HY, Kou SQ, Lindqvist PA, Tang CA. Numerical simulation of shear fracture (mode II) in heterogeneous brittle rock. Int J Rock Mech Min Sci. 2004;41:14–19. 13. Ning YJ, Yang J, An XM, Ma GW. Modeling rock fracturing and blast-induced rock mass failure via advanced discretization within the discontinuous deformation analysis framework. Comput Geotech. 2011;38:40–49. 14. Ning YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. Int J Rock Mech Min Sci. 2011;48:964–975. 15. Wu ZJ, Wong LNY. Frictional crack initiation and propagation analysis using the numerical manifold method. Comput Geotech. 2012;39:38–53. 16. Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock Eng. 2015;48:1097–1114. 17. Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Methods Eng. 1999;46:131–150. 18. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Methods Eng. 1999;45:601–620. 19. Wagner GJ, Moes N, Liu WK, Belytschko T. The extended finite element method for rigid particles in Stokes flow. Int J Numer Methods Eng. 2001;51:293–313. 20. Belytschko T, Krongauz Y, Organ D, Flemming M, Krysl. P. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng. 1996;139:3–47. 21. Libersky LD, Petsehek AG. Smoothed particle hydrodynamics with strength of materials. Adv Free Lagrange Method Lect Notes Phys. 1990;26:248–257. 22. Izumi S, Katake S. Molecular dynamics study of solid deformation. Trans Jpn Soc Mech Eng A. 1993;59:263–267. 23. Jin Y, Wang Y, Khachaturyan A. Three-dimensional phase field microelasticity theory and modeling of multiple cracks and voids. Appl Phys Lett. 2001;79:3071–3073. 24. Silling SA, Bobaru F. Peridynamic modeling of membranes and fibers. Int J Nonlin Mech. 2005;40:395–409. 25. Silling SA. Dynamic fracture modeling with a meshfree peridynamic code. In: Bathe KJ, editor. Computational fluid and solid mechanics. Amsterdam: Elsevier; 2003. p. 641–644. 26. Silling SA, Askari E. Peridynamic modeling of impact damage. In: Moody FJ, editor. PVP, vol. 489. New York: American Society of Mechanical Engineers; 2004. p. 197–205. 27. Silling SA, Zimmermann M, Abeyaratne R. Deformation of a peridynamic bar. J Mech Phys Solids. 2003;73:173–190. 28. Huang D, Zhang Q, Qiao PZ, Shen F. Peridynamic methods and application. Adv Mech. 2010;4:448–459. 29. Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids. 2000;48:175–209. 30. Gunzburger M, Lehoucq RB. A nonlocal vector calculus with application to nonlocal boundary value problems. Multiscale Model Simul. 2010;8:1581–1589. 31. Ha YD, Bobaru F. Characteristics of dynamic brittle fracture captured with peridynamics. Eng Fract Mech. 2011;78:1156–1168. 32. Belytschko T, Chen H, Xu J, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int J Numer Methods Eng. 2003;58:1873–1905.
254
X.-P. Zhou et al. / International Journal of Rock Mechanics & Mining Sciences 80 (2015) 241–254
33. Emmrich E, Weckner O. The peridynamic equation and its spatial discretization. Math Model Anal. 2007;12:1–10. 34. Gao G, Huang S, Xia K, Li Z. Application of digital image correlation (DIC) in dynamic notched semi-circular bend (NSCB) tests. Exp Mech. 2015;55:95–104. 35. Zhou. YX, Xia. K, Li XB, Li HB, Ma GW, Zhao J, Zhou ZL, Dai F. Suggested methods for determining the dynamic strength parameters and mode-I fracture toughness of rock materials. Int J Rock Mech Min Sci. 2012;49:105–112. 36. Liu WF. Discretized bond-based peridynamics for solid mechanics [Ph.D. thesis]. Michigan State University; 2012. 37. Dyskin AV. On the possibility of bifurcation in linear periodic arrays of 2-D cracks. Int. J. Fract.. 1994;67:R31–R36.
38. Muhlhaus HB, Chau KT, Ord A. Bifurcation of crack pattern in arrays of twodimensional cracks. Int. J. Fract.. 1996;77:1–14. 39. Oguni K, Hori M, Ikeda K. Analysis on evolution pattern of periodically distributed defects. Int J Solids Struct. 1997;34:3259–3272. 40. Chau KT, Wang GS. Condition for the onset of bifurcation in en echelon crack arrays. Int J Numer Anal Methods Geomech. 2001;25:289–306. 41. Zhou XP, Xie WT, Qian QH. Bifurcation of collinear crack system under dynamic compression. Theor Appl Fract Mech. 2010;54:166–171. 42. Zhou XP, Qian QH, Yang HQ. Bifurcation condition of crack pattern in the periodic rectangular array. Theor Appl Fract Mech. 2008;49:206–212.