Modeling wave propagation induced fracture in rock with correlated lattice bond cell

International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

Contents lists available at ScienceDirect

International Journal of Rock Mechanics & Mining Sciences journal homepage: www.elsevier.com/locate/ijrmms

Modeling wave propagation induced fracture in rock with correlated lattice bond cell Zhennan Zhang a,n, Yuan Yao a, Xianbiao Mao b a b

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221008, China

art ic l e i nf o

a b s t r a c t

Article history: Received 19 January 2015 Received in revised form 21 May 2015 Accepted 4 June 2015 Available online 31 July 2015

The correlated lattice bond cell model (CLBC) is used to account for the characteristics of the mesostructure of rock. Each mineral grain is modeled as a discrete unit bond cell. A linear modified Stillinger– Weber potential is adopted to characterize the unit bond cell. By this method, the lattice model can not only represent the variable Poisson ratio, but can also simulate more complex fracture behaviors than the two-body interaction-based lattice model. It can precisely simulate the stress wave propagation in rock with different Poisson ratios. With a simple bond rupture criterion, the spalling fracture induced by stress wave is well simulated. To represent the heterogeneity of rock, the moduli of bond cells are assumed to be random, following a Weibull distribution. The simulation results demonstrate that the distinct spalling fracture is likely to initiate in a hard rock subjected to impact load, but not in a soft rock. The heterogeneity has little effect on the spalling fracture. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Lattice model Stress wave Dynamic fracture Rock Spalling fracture Modified Stillinger–Weber potential

1. Introduction The stress wave propagation can induce fracture in rock, leading to many geological disasters. However, the existed mechanical models are still inadequate to effectively simulate this phenomenon. The continuum mechanics-based method is the most popular approach to this problem by far, but it has many limitations in modeling stress wave-induced fracture. The fracture behaviors, such as initiation, propagation, branching and arresting, cannot be simulated unless a suitable fracture criterion is given in prior. However, the fracture criterion, especially that of the unstable dynamic fracture, is very difficult to provide. The lattice method pioneered by Hrennikoff1 provides an efficient approach to the simulation of dynamic fracture. For the lattice method represents a continuum with a discrete bond system, it simulates the continuum fracture process by the successive bond ruptures with a simple bond rupture criterion. The bond rupture criterion is much easier to establish than the fracture criterion of 3D continuum. In the early version of this method, the discrete structure was composed of springs.2,3 Because the spring only accounts for the normal interaction of two particles, the Poisson ratio it presents is fixed, namely 0.25 for 3D case and 1/3 for 2D case. This restricts its application to a wider range of materials. To n

Corresponding author. Fax: þ86 21 34204346. E-mail address: [email protected] (Z. Zhang).

http://dx.doi.org/10.1016/j.ijrmms.2015.06.006 1365-1609/& 2015 Elsevier Ltd. All rights reserved.

address this problem, the beam was introduced into the particle system to replace the spring.4,5 However, in the beam system the particle rotation, an additional degree of freedom to the particle displacement, was introduced. To represent the variable Poisson ratio only in terms of particle displacement, Zhao et al.6 developed the distinct lattice spring model. In this method, the normal and the shear deformation of spring are separately accounted for. However, the shear deformation of spring directly computed by the particle displacement includes the rigid body rotation. To exclude the rigid body rotation, a local strain-based method was adopted in.6 In order to account for the bond rotation effect, Wang et al.7,8 introduced the angular spring to the lattice model. By this method, the variable Poisson ratio can be represented. These advances of lattice method make the fracture simulation more efficient. However, in the lattice method, the bond parameter calibration is case-dependent. Different topological structure of unit cell leads to different micro–macro parameter relationship. This makes it difficult to calibrate the bond parameters of 3D lattice model. Moreover, the conventional lattice model is incapable of presenting the characteristics of mesostructure of rock, which consists of mineral grains on mesoscale. To develop a lattice model that can represent the characteristics of mesostructure of rock, Zhang9 developed the discretized virtual internal bond method (VIB) based on the VIB theory.10 The discrete structure of discretized-VIB consists of unit bond cells. Each unit bond cell can have any geometry with any number of bonds. A universal bond parameter calibration method was

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

proposed based on the ideal unit cell conception.9 The discretizedVIB employs the two-body potential to characterize the bond energy. The two-body potential only accounts for the normal interaction between two particles, missing the effect of bond rotation. Therefore, the discretized-VIB presents a fixed Poisson ratio. Furthermore, from the standpoint of fracture mechanism, the twobody potential indicates that the fracture initiation and propagation only result from the normal separation of two particles, irrelevant to the shear effect between particles. Therefore, the twobody potential is too simple to simulate more complex fracture behaviors. Instead, the Stillinger–Weber (SW) potential11 indicates that the fracture initiation is not only related to the normal deformation of bond, but also related to the bond angles. Though SW potential is specially developed for the silicon material, the fracture mechanism it represents is applicable to other materials. Zhang et al.12 modified the SW potential and extended it to other materials than silicon. To address the fixed Poisson ratio problem of the discretized-VIB, Zhang and Chen13 adopted the modified SW potential to describe the unit bond cell. For the bonds in a unit cell are correlated together by the SW potential, the discretized-VIB employing the SW potential is also called the correlated lattice bond cell model (CLBC). Since a more elaborate fracture mechanism is incorporated, the CLBC can capture more characteristics of dynamic fracture. Rock is a kind of heterogeneous quasi-brittle material. To develop a suitable lattice model for rock based on CLBC, we use the linear modified SW potential to describe the unit cell of rock and introduce the heterogeneity on the mineral grain level. With this model, we simulate the stress wave propagation in rock with variable Poisson ratios and the induced spalling fracture. The aim of the present work is to provide a simple and efficient lattice method to simulate the stress wave-induced fracture in rock.

2. Modeling rock with correlated lattice bond cells 2.1. Brief introduction to Stillinger–Weber potential and its modified version The SW potential11 was originally proposed to simulate the silicon material in molecular dynamics methodology. It takes the form

Φ=

j
(1)

where Φ2, Φ3 stands for the two- and three-body interaction, respectively; rij the bond length of bond vector rij ; θjik the bond angle subtended by the bond vector rij , rik at vertex i . The specific expressions of Φ2 and Φ3 are.

⎡ ⎛ ⎞ p ⎛ ⎞q⎤ ⎛ ⎞ ℓ ℓ 1 ⎟⎟ Φ2 rij = A⎢B⎜⎜ 0 ⎟⎟ − ⎜⎜ 0 ⎟⎟ ⎥exp⎜⎜ * ⎢ ⎝ rij ⎠ ⎥ ⎝ rij ⎠ ⎦ ⎝ rij/ℓ0 − r0 ⎠ ⎣

( )

(

the ‘ideal’ tetrahedral angle as the reference value of bond angle in the current configuration. Such three-body interaction tends to arrange the randomly distributed bonds into an ‘ideal’ tetrahedral structure. Thus, the original SW potential is only applicable to the silicon material, not to others. In Ref. 12, the SW potential was modified in such manner that the bond angle in the reference configuration, not the ‘ideal’ tetrahedral angle, was taken as the reference value of bond angle in the current configuration. The modified SW potential12 is

⎡ Φ2(r˜I ) = A(Br˜I−p − r˜I−q)exp⎢⎣ r˜I − r0*

−1⎤

)

⎥⎦

i≠j

i
Φ3 rij, rik, θjik

Fig. 1. Modeling rock with correlated lattice bond cells (a) the mesostructure of rock (SEM image of marble) consisting of mineral grains; and (b) a continuous mineral grain is modeled as a discrete unit bond cell.

(

∑ Φ2(rij) + ∑ Φ3(rij, rik, θijk)

263

)

2 ⎛ ⎞ ⎛ ⎞ ⎛ 1⎞ γ γ ⎟⎟exp⎜ = λ⎜cos θjik + ⎟ exp⎜⎜ ⎟ * * ⎝ r r r r 3⎠ / / ℓ − ℓ − ⎝ ⎝ ij 0 0⎠ ik 0 0 ⎠

2 ⎡ Φ3 r˜I , r˜J , θIJ = λ cos θIJ − cos θIJ 0 exp⎢γ r˜I − r0* ⎣

(

)

(

)

(

−1

)

(

−1⎤

)

+ γ r˜J − r0*

⎦⎥ (3)

in which r˜I , r˜J are the normalized bond length, r˜I = rij/ℓ0, r˜J = rik/ℓ0 ; θIJ0 stands for the bond angle subtended by the bond rij and rik in the reference configuration whereas θIJ its value in the current configuration. 2.2. Modeling mineral grain with correlated lattice bond cells

(2)

in which r0* is the cutoff radii; the 1/3 being the minus cosine of the’‘ideal’ tetrahedral angle; ℓ0 is the bond length when undeformed. The parameters A and λ implicitly contain the energy scaling parameter ε of the original two- and three-body interactions in Ref. 11. Eq. (2) indicates that the bond energy is not only related to the normal deformation of itself, but also related to the bond angles subtended with other bonds. This equation is especially applicable to the silicon materials because the three-body interaction takes

The mesostructure of rock is shown in Fig. 1a, which consists of mineral grains on the mesoscale. Originally, the mineral grain is a continuum. In the present method, it is modeled as a discrete system, i.e., a correlated lattice bond cell (Fig.1b). Considering that rock is a typical elastic-brittle material, the energy of a unit bond cell is characterized by the linear modified SW potential.

Φ2 =

1 2 2 Aℓ 0(r˜I − 1) 2

Φ3 =

1 λ θIJ − θIJ 0 2

(

2

)

(4)

264

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

A=

2V 3E ⋅ N (N − 1)ℓ20 (1 − 2ν )

λ=

V 9E(1 − 4ν ) ⋅ N (N − 1)(N − 2) 2(1 + ν )(1 − 2ν )

(8)

where V is the volume of a unit cell.

2.4. Bond rupture criterion

Fig. 2. Conjugate bond vectors subtended at a given particle i in a discrete system with N particles. The total particle number is N . The bond number of rij is N (N − 1) and its conjugate bond rik number is N − 2 with rij ≠ r ji .

Before a bond breaks, it is linear elastic. To simulate the failure process of rock, a bond rupture criterion has to be imposed. In the lattice model, either the bond force or the bond deformation rupture criterion is adopted, e.g. Ref. 7,14. In the present paper, the bond strain criterion is adopted. That means when the bond strain exceeds a critical value, the bond breaks. By the same way, we set the critical bond angle over which the bond angle will lose the mechanical resistance. Since a bond angle is subtended by two bonds, its mechanical resistance is also subjected to the status of the two bonds. As long as one of them breaks, the bond angle is failed. So, the bond potential with bond rupture criterion can be written as

Through the three-body interaction, the bonds in a unit cell are correlated together. Fig. 2 shows the relationship between a bond and its conjugated bonds.

⎧1 2 2 ˜ r˜I − 1 ≤ εr ⎪ ⎪ 2 Aℓ 0(rI − 1) if Φ2 = ⎨ ⎪1 2 2 Aℓ ε else ⎪ ⎩2 0 r

2.3. Derivation of force-displacement constitutive relation of a unit cell

⎧1 ⎪ λ θ − θIJ 0 Φ3 = ⎨ 2 IJ ⎪ ⎩ constant

Since a mineral grain is modeled as a discrete system, the constitutive relation of the mineral grain is directly characterized by the particle force-displacement relation, bypassing the socalled stress–strain relation. Let u denote the particle displacement vector of a unit cell, f the corresponding force vector and K the stiffness matrix. For a unit cell with N particles, by Eq.(1), the total energy of a unit cell is identified as

where εr is the critical bond length strain over which bond breaks; εθ is the critical bond angle strain over which bond angle loses mechanical resistance. Here, the ‘constant’ in the three-body interaction is equal to the potential value at which bond angle loses resistance. The bond force and bond moment derived from Eq.(9) are respectively

Φ=

1 2

N (N − 1)

∑

1 2

Φ2(r˜I ) +

I

N (N − 1) N − 2

∑ ∑ I

( )

Φ3 θIJ

(5)

J

in which N denotes the total particle number in a unit cell with rij ≠ rji . By Eq. (5), the particle force vector of a unit cell is derived as

Fi =

1 ∂Φ = 2 ∂ui

N (N − 1)

∑ I

∂Φ2 ∂r˜I 1 ⋅ + 2 ∂r˜I ∂ui

N (N − 1) N − 2

∑ ∑ I

J

∂Φ3 ∂θIJ ⋅ ∂θIJ ∂ui

(6)

and the stiffness matrix is

Kij =

1 ∂ 2Φ = 2 ∂ui ∂uj +

N (N − 1)

∑ I

1 2

⎛ ∂ 2Φ ∂r˜ ∂r˜ ∂Φ2 ∂ 2r˜I ⎞ 2 I I ⎜⎜ ⎟⎟ ⋅ ⋅ + ⋅ 2 ∂r˜I ∂ui ∂uj ⎠ ⎝ ∂r˜I ∂uj ∂ui

N (N − 1) N − 2

∑ ∑ I

(

J

2 ⎛ ∂ 2Φ ∂θ ∂θ ∂Φ3 ∂ θIJ ⎞⎟ IJ IJ 3 ⎜ ⋅ ⋅ + ⋅ ⎜ 2 ∂u ∂u ∂θIJ ∂ui ∂uj ⎟⎠ j i ⎝ ∂θIJ

(7)

To calibrate the physical parameters (A, λ ) of Eq. (4), Zhang et al.12,13 took use of the concept of ideal unit cell and derived the relationship between micro parameters (A, λ ) and macro material constants (E, ν ) as

fn =

mIJ =

2

)

(

)

max r˜I , r˜J ≤ 1 + εr and θIJ − θIJ 0 ≤ εθ θIJ 0

if else

(9)

r˜ − 1 ≤ εr ∂Φ2 ⎧ Aℓ0(r˜ − 1) if =⎨ ∂r ⎩0 else ⎧λ θ − θ ∂Φ3 ⎪ IJ 0 ⎨ IJ =⎪ ∂θIJ ⎩0

(

)

if else

(

)

max r˜I , r˜J ≤ 1 + εr and θIJ − θIJ 0 ≤ εθ (10)

where fn means the bond force and mIJ the bond moment. The curves of the bond force-length and bond moment-angle are shown in Fig. 3. In Fig. 3, it is clearly shown that an abrupt drop in the bond force and bond moment occurs when the bond deformation exceeds a critical value. Such abrupt failure pattern is quite suitable to characterize the failure process of rock as suggested in Ref. 15. In Ref. 15, the rock failure process from the elastic-brittle to the residual plastic behavior was analyzed, but in the present paper, the plastic behavior of rock is not accounted. To check the effects of εr and εθ on the macro mechanical response of material, we simulate a uniaxial tensile specimen composed of tetrahedral unit cells (Fig. 4a) with the bond rupture criterion. The input material constants are E = 10 GPa and ν = 0.2. The check results are shown in Fig.4b, which suggests that εr has significant effect on the peak stress whereas εθ has only a little effect on the post-peak behavior.

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

265

3.2. Computation of bond angle and its derivatives The cosine of the bond angle between RI , RJ in the reference and current configuration can be respectively calculated as

(X¯ X)(Y¯ Y¯ )

¯ T Y¯ / cos θIJ 0 = X cos θIJ = x¯ T y¯ /

T

T

(x¯ x¯ )(y¯ y¯ ) T

T

(13)

The first derivative of θIJ with respect to the particle displacement vector is

θ′IJ (u) =

⎛ ∂Q 1 T ⎜ ⋅ 12 − 2⎜ Q Q u 2 Q ∂ui ∂ ⎝ i 11 1−T 11 22 ⎞ ∂Q T ∂Q 22 ⎟ ⋅ 11 − ⋅ 2Q 22 ∂ui ⎟⎠ ∂ui ∂θIJ (u)

−1

=

rITrI ,

(14)

rITrJ ,

rTJ rJ ,

where the scalars Q11 = Q12 = Q 22 = T = Q12/ Q11Q 22 , θIJ = arccos(T ); their derivatives with respect to node displacement vector are T ∂Q 12/∂u i = ⎡⎣ −rI1 − r J1, −rI 2 − r J 2, −rI 3 − r J 3, r J1, r J 2, r J 3, rI1, rI 2, rI 3 ⎤⎦ T ∂Q 11/∂u i = 2⎡⎣ −rI1, −rI 2, −rI 3, rI1, rI 2, rI 3, 0, 0, 0 ⎤⎦

Fig. 3. Depiction of bond rupture criterion (a) normal bond force versus bond length; (b) bond moment mIJ versus bond angle. The bond moment mIJ = rIfIS = r Jf JS . As long as one bond length exceeds the critical value, the bond angle loses its mechanical resistance.

3. Implementation aspects 3.1. Computation of bond length derivatives Eqs. (6) and (7) involves the bond length, bond angle and their derivatives. To facilitate the implementation of the present method, we briefly list the computation of these quantities. For details, refer to Ref. 13. Take a pair of conjugate vectors as example (Fig.5). The bond vectors in the reference configuration are RI = X2 − X1 and RJ = X3 − X1 with X1, X2, X3 being the coordinate vectors of the particle 1,2,3. Their corresponding bond vectors in the current configuration are

r J = R J + u3 − u1

placement vector can be written as u = [u1, u2 , u3, u4 , u5, u6]T . The derivatives of the normalized bond length with respect to the node displacement vector are

∂r˜I 1 ⎡ = ⎣rI1(δ4i − δ1i ) + rI2(δ5i − δ2i ) + rI 3(δ6i − δ3i )⎤⎦ ∂ui ℓ0ℓ ∂ 2r˜I ℓ ∂r˜ ∂r˜ 1 p − 0⋅ I ⋅ I = ∂ui ∂uj ℓ0ℓ ij ℓ ∂ui ∂uj RTI RI ; ℓ =

(12) rITrI ; rIi is the ith component of the

(

)

(

)

bond vector rI ; pij = (δ4i − δ1i) δ4j − δ1j + (δ5i − δ2i) δ5j − δ2j + (δ6i

(

−δ3i) δ6j − δ3j

)

3.3. Constitutive relation with determined parameters By Eqs. (4) and (6–8), we have the node force vector of the unit cell

Fi =

∂Φ 1 = ∂ui 2

with δij being the Kronecker delta.

N (N − 1)

∑ I

∂r˜ 1 A˜ (r˜I − 1)⋅ I + 2 ∂ui

N (N − 1) N − 2

∑ ∑ I

(

J

∂θ

) ∂uIJ

λ θIJ − θIJ 0 ⋅

i

(16)

and the stiffness matrix

∂ 2Φ 1 = ∂ui ∂uj 2

(11)

where u1, u2, u3 being the displacement vectors of the particle 1,2,3. For the bond rI with local particle index 1,2, its particle dis-

(15)

with r Ji being the ith component of the bond vector rJ . Because the second derivative of bond angle with respect to particle displacement is extremely complicated, we omit it in the stiffness matrix, Eq. (7), for the sake of simplification. Such omission doesn’t change the convergence nature of computation according to the principle of the Newton–Raphson method.

Kij =

rI = R I + u2 − u1

in which ℓ0 =

T ∂Q 22/∂u i = 2⎡⎣ −r J1, −r J 2, −r J 3, 0, 0, 0, r J1, r J 2, r J 3 ⎤⎦

N (N − 1)

1 + 2

∑ I

⎛ ∂r˜ ∂r˜ ∂ 2r˜I ⎞ ⎟⎟ A˜ ⎜⎜ I ⋅ I + (r˜I − 1)⋅ ∂ui ∂uj ⎠ ⎝ ∂uj ∂ui

N (N − 1) N − 2

∑ ∑ I

λ⋅

J

∂θIJ ∂θIJ ⋅ ∂uj ∂ui

(17)

where the physical parameters are A˜ = 2V /[N (N − 1)]⋅3E /(1 − 2ν ) and λ = V /[N (N − 1)(N − 2)]⋅9E (1 − 4ν ) /[2(1 + ν )(1 − 2ν )]. 3.4. Consideration of heterogeneity A typical feature of rock material is its heterogeneity. The mechanical property varies with location in rock. To model the heterogeneity of rock, the Young’s modulus of a unit cell is assumed to be random number, which complies with a certain statistical distribution law. According to Refs. 16,17, the Weibull distribution law is suitable to characterize the Young's modulus distribution in rock:

f (x) =

⎛ ⎛ x ⎞m⎞ m ⎛⎜ x ⎞⎟m − 1 exp⎜−⎜ ⎟ ⎟ Ε¯ ⎝ Ε¯ ⎠ ⎝ ⎝ Ε¯ ⎠ ⎠

(18) in which Ε¯ denotes the mean value of Young's modulus; m the

266

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

Fig. 4. Influence of the critical bond strain on macro mechanical response: (a) the simulated uniaxial tensile specimen; and (b) the effects of critical bond strains on the uniaxial tensile stress–strain curve. ( E = 10GPa , ν = 0.2).

4. Stress wave simulation

Fig. 5. Displacement of a pair of conjugate bonds.

index of homogeneity level. The higher the m value, the stronger the homogeneity. Fig. 6 shows the gray images of rock sample with different homogeneity levels. Fig. 6d suggests that when m = 3, the rock has become strongly heterogeneous.

To simulate the stress wave propagation induced fracture, the mechanical model should be firstly capable of simulating stress wave propagation. To verify the present method in simulating stress wave propagation, we use it to simulate a 2D case reported in [7]. The schematic of computation domain setup is shown in Fig. 7. To represent the 2D case with 3D element, we let this plate have a thickness of an element. We equally discretized the plate into 100 × 10 × 1 cubes. Each cube is subdivided into 5 tetrahedral unit cells. The material parameters of plate are: Young’s modulus E = 50.0 GPa ; Poisson ratio ν = 0.1; material density ρ = 2400kg/m3. In the dynamic boundary condition case where a constant stress 0.1 MPa with duration 1.0μs is applied on the left end of plate, the simulated displacements at the designated points in Fig. 7 are shown in Figs. 8a and 8b. The profile of the simulated horizontal displacement-time curves (Fig. 8a) is consistent with that reported in Refs 7 and that observed in the usual Hopkinson

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

267

Fig. 6. Gray images of rock sample with different homogeneity levels (a) Weibull distribution density; (b) m = 50 ; (c) m = 10 ; (d) m = 3. The right gray column denotes the Young's modulus (GPa). The mean value is 38.4 GPa. The generated Young's modulus distribution density (k ¼38.4 GPa).

Fig. 7. Schematic of computation domain setup for 2D stress wave propagation with two designated points.

compressive bar experiment. The magnitude of oscillation in the transverse direction is much smaller than that in the longitudinal direction, shown as Fig. 8b.

In the kinematic boundary condition case where a constant velocity 0.01m /s is applied on the left end of plate, the simulated displacement at point X = L/4 and X = L/2 is shown in Fig.8c, which is consistent with that reported in [7]. To further check whether the present method can quantitatively simulate the wave speed in rock with variable Poisson ratios, we simulate the wave speed with different Poisson ratio inputs in the kinematic boundary condition case. The theoretical longitudinal stress wave speed is

Cp(2D) =

(1 − ν )E (1 + ν )(1 − 2ν )ρ

(19)

268

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

Table 1 Simulated stress wave speed and its relative error with different Poisson ratios. The → Kinetic Boundary Condition is V = 0.01m/s; Speed unit: m/s; Material parameters: E = 50GPa , ρ = 2400kg/m3. Poisson ratio

0.1

0.15

0.20

0.25

Simulated speed Theoretical speed Relative error (%)

4488.1 4615.8  2.766

4573.9 4690.1  2.478

4674.0 4811.4  2.856

4817.0 4999.9  3.660

Table 2 Simulated stress wave speed with different boundary conditions. Material parameters: E = 50 GPa ,ν = 0.15, ρ = 2400kg/m3; the theoretical wave speed: 4690.1 m/s. KBC

Simulated wave speed Relative error (%)

DBC

0.01 m/s

0.1 m/s

0.1 MPa

0.05 MPa

4573.9  2.48

4559.5  2.78

4545.2  3.09

4645.3  0.95

Fig. 9. Setup for the simulation of stress wave-induced spalling fracture (cm).

The simulated stress wave speeds and their theoretical values are listed in Table 1. From Table 1, it is seen that the simulated wave speeds are very close to their theoretical values. The relative errors are about 3.0%. The stress wave speed is only related to the material properties, independent of the boundary conditions. To check whether the simulated stress wave speed is boundary condition-independent, we simulate the stress wave speed with different boundary conditions. The check results (Table 2) indicate that simulated wave speed is indeed independent of the boundary conditions.

5. Stress wave-induced spalling fracture 5.1. Simulation of spalling fracture process

Fig. 8. Simulated displacement at the designated point X ¼ L/4 and X ¼L/2 (a) horizontal displacement vs. time with Dynamic Boundary Condition; (b) vertical displacement versus time with DBC; (c) horizontal displacement versus time with Kinematic Boundary Condition.

To examine whether the present method can simulate the stress wave-induced spalling fracture, we simulate a rock beam (target) subjected to an impact load, shown in Fig. 9. A flyer im→ pacts the target at the speed V = 260.0 m/s. The flyer is linear elastic with material parameters: E = 60 GPa; ν = 0.2; ρ = 2400 kg/m3. The target is elastic-brittle with material parameters: E = 60 GPa; ν = 0.2; ρ = 2400kg/m3; εr = εθ = 0.1 × 10−3. After the flyer impacts on the target, a spalling fracture will initiate near the free boundary of target due to the stress wave propagation and reflection. The simulated spalling fracture process is

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

269

Fig. 10. Simulated spalling fracture process (a)t = 0.2 μs ; (b) t = 3.0 μs ; (c)t = 4.0 μs ; (d) t = 5.0 μs ; (e) t = 6.0 μs . (Node displacements are magnified 5 times.).

Fig. 11. Spalling fracture pattern simulated by MD [18,19].

shown in Fig. 10. From Fig. 10 it is seen that a spalling fracture is clearly observed near the free boundary of target. The pattern of the spalling fracture is similar to that simulated by molecular dynamics in Refs. 18,19, shown as Fig. 11. The reason that the spalling fracture initiates near the free boundary of target can be indicated in Fig. 12a. This figure shows the simulation results in which the target is modeled as a linear elastic continuum. According to Fig. 12a, the peak tensile stress appears near the free boundary, say about 0.2 cm from the free boundary, after the stress wave was reflected. Hence, the spalling fracture should initiate at there. The extraordinary large tensile strain in Fig. 12b indicates the spalling fracture indeed initiates where the peak stress tensile stress occurs. Fig. 10 also confirms that a spalling fracture clearly initiates at that location. 5.2. Effect of material properties on spalling fracture To examine the effect of Young’s modulus and homogeneity level of rock on the spalling fracture, we change the modulus and homogeneity level of target and re-simulate the example above. The simulation results are shown in Fig. 13, which suggests that the Young’s modulus has significant influence on the spalling fracture. When the modulus is high, e.g., E = 60.0 GPa , a distinct spalling fracture is observed, shown as Fig. 13a-b. However, when the modulus is low, e.g. E = 10.0 GPa , no distinct spalling fracture is observed. Instead, the granular loosening is observed at the free

Fig. 12. Simulated vertical normal stress and strain distribution along the middle line of target: (a) normal stress distribution where the target is modeled as a linear elastic continuum with parameters: E = 60 GPa , ν = 0.2, ρ = 2400kg/m3; (b) normal strain distribution where the target is modeled as the elastic-brittle unit bond cell with parameters: E = 60 GPa , ν = 0.2, ρ = 2400kg/m3 , εr = εθ = 0.1 × 10−3.

boundary in the soft rock cases, shown as Fig. 13c-d. Compared with the modulus, the homogeneity level has little influence on the spalling fracture, shown as Fig. 13b and d.

6. Conclusion remarks The correlated lattice bond cell method with linear modified Stillinger–Weber potential is an effective approach to modeling rock material. This lattice bond model can precisely simulate the stress wave speed in rock with variable Poisson ratios. With the bond strain rupture criterion, the present method can simulate the

270

Z. Zhang et al. / International Journal of Rock Mechanics & Mining Sciences 78 (2015) 262–270

China University of Mining & Technology (Grant SKLGDUEK1203). They are gratefully acknowledged.

no.

Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ijrmms.2015.06. 006.

References

Fig. 13. Simulated spalling fractures with different Young's modulus and homogeneity levels (a) E = 60.0 GPa, m = ∞; (b) E = 60.0 GPa, m = 3; (c) E = 10.0 GPa, m = ∞; (d) E = 10.0 GPa, m = 3. ( t = 20.0 μs ).

spalling fracture process induced by stress wave. It can capture the characteristic of spalling fracture. The modulus of rock has significant impact on the spalling fracture. In the hard rock, the spalling fracture distinctly initiates while in the soft rock no distinct spalling fracture is observed. The heterogeneity of rock has little impact on the spalling fracture. The present paper provides an efficient lattice method to simulate the stress wave-induced fracture in rock.

Acknowledgments The present work is supported by the National Natural Science Foundation of China (no. 11172172), the National Basic Research Program of China (Grant no.2011CB013505) and the State Key Laboratory for GeoMechanics and Deep Underground Engineering,

[1] Hrennikoff A. Solution of problems of elasticity by the framework method. J Appl Mech. 1941:A169–A175. [2] Greenspan D. Supercomputer simulation of cracks and fractures by quasimolecular dynamics. J Phys Chem Solids. 1989;50:1245–1249. [3] Wang G, Ostoja-Starzewski M. Particle modeling of dynamic fragmentation – I: theoretical considerations. Comput Mater Sci. 2005;33:429–442. [4] Griffiths DV, Mustoe GGW. Modeling of elastic continua using a grillage of structural elements based on discrete element concepts. Int J Numer Methods Eng. 2001;50:1759–1775. [5] Lilliu G. Van Mier JGM. 3D lattice type of fracture model for concrete. Eng Fract Mech. 2003;70:927–941. [6] Zhao GF, Fang J, Zhao J. A 3D distinct lattice spring model for elasticity and dynamic failure. Int J Numer Anal Methods Geomech. 2011;35:859–885. [7] Wang G, Al-Ostaz A, Cheng AHD, Mantena PR. Hybrid lattice particle modeling of wave propagation induced fracture of solids. Comput Methods Appl Mech Eng. 2009;199:197–209. [8] Wang G, Al-Ostaz A, Cheng AHD, Mantena PR. Hybrid lattice particle modeling: theoretical considerations for a 2D elastic spring network for dynamic fracture simulations. Comput Mater Sci. 2009;44:1126–1134. [9] Zhang Z. Discretized virtual internal bond model for nonlinear elasticity. Int J Solids Struct. 2013;50:3618–3625. [10] Gao HJ, Klein PA. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bonds. J Mech Phys Solids. 1998;46:187– 218. [11] Stillinger FH, Weber TA. Computer simulation of local order in condensed phases of silicon. Phys Rev B. 1985;31:5262–5271. [12] Zhang Z, Chen Y, Zheng H. A modified Stillinger–Weber potential-based hyperelastic constitutive model for nonlinear elasticity. Int J Solids Struct. 2014;51:1542–1554. [13] Zhang Z, Chen Y. Modeling nonlinear elastic solid with correlated lattice bond cell for dynamic fracture simulation. Comput Methods Appl Mech Eng. 2014;279:325–347. [14] Ostoja-Starzewski M, Sheng P, Alzebdeh K. Spring network models in elasticity and fracture of composites and polycrystals. Comput Mater Sci. 1996;7:82–93. [15] Zheng H, Liu DF, Lee CF, Ge XR. Principle of analysis of brittle-plastic rock mass. Int J Solids Struct. 2005;42:139–158. [16] Tang CA, Liu H, Lee PKK, Tsui Y, Tham LG. Numerical studies of the influence of microstructure on rock failure in uniaxial compression – Part I: effect of heterogeneity. Int J Rock Mech Min Sci. 2000;37:555–569. [17] Tang CA, Wong RHC, Chau KT, Lin P. Modeling of compression-induced splitting failure in heterogeneous brittle porous solids. Eng Fract Mech. 2005;72:597–615. [18] Krivtsov AM. Relation between spall strength and mesoparticle velocity dispersion. Int J Impact Eng. 1999;23:477–487. [19] Krivtsov AM. Structural Dynamics – EURODYN’99. Rotterdam: Balkema; 1999. p. 475–477.