# Numerical Study on Crack Propagation by Using Softening Model under Blasting

Authors: Rong Hu; Zheming Zhu; Jun Xie; Dingjun Xiao
Journal: Advances in Materials Science and Engineering
 (2015)

---

## Abstract

A mixed failure criterion, which combined the modified maximum principal stress criterion with the damage model of tensile crack softening, was developed to simulate crack propagation of rock under blasting loads. In order to validate the proposed model, a set of blasting models with a crack and a borehole with different incident angles with the crack were established. By using this model, the property of crack propagation was investigated. The linear equation of state (EOS) was used for rock, and the JWL EOS was applied to the explosive. In order to validate the numerical simulation results, experiments by using PMMA (polymethyl methacrylate) with a crack and a borehole were carried out. The charge structure and incident angle of the blasting experimental model were the same as those in the numerical models. The experiment results agree with the numerical simulation results.


---

## Body

## 1. Introduction

Blasting is an economic and efficient excavation method and is widely used in engineering of quarrying, mining, and tunnel excavation. The propagation of blasting-induced stress wave in cracked rock mass would introduce crack expansion and failure, which would possibly induce destructively geologic hazards, such as rock burst and slope slump. The propagation property of stress wave in cracked rock mass and its dynamic mechanical response have drawn massive attention since plentiful joints and cracks exist within rock mass. It is an important topic in rock dynamics and relevant disciplines to study the way cracks propagate under dynamic loading so as to know how the wing cracks develop at crack tips [1] as the case in static compressive loading, which has been well investigated by using complex function method [2–5]. This research helps with predicting the dynamic strength and structural stability of cracked rock mass and provides theoretical basis for blasting design, improving blasting efficiency.The time taken by the explosive loading is extremely short, ranging, usually, from a couple of microseconds to dozens of microseconds. That renders the recording of the whole physical field impossible experimentally, whereas the state of the field can be described by means of numerical simulation, which is less expensive, and can be easily realized. However, the numerical simulation accuracy depends largely on the quality of the numerical model employed. Hence, a combined use of numerical simulation and experiment was carried out in this paper. In the numerical study, Preece and Thorne [6] used 3D finite element techniques and a damage constitutive model to study the detonation timing and fragmentation. Donzé et al. [7] applied a model based on the discrete element to investigate the importance of stress waves on the initiation and propagation of radial fractures during the dynamic loading phase of explosion. By using UEDC code, Chen and Zhao [8] have simulated blasting wave propagation in joint rock mass. Using AUTODYN code and the modified maximum principal stress, Zhu et al. [9, 10] established dynamic computation model for rock under dynamic loading, and by using this model, the crater blasting [11] and the zonal disintegration phenomenon [12] have been investigated. However, brittle materials, such as rock, act in inelastic way, which is caused by the nucleation, propagation, and connection of microcracks. Softening models are often used to describe the inelasticity of brittle material [13, 14]. Feenstra and De Borst [15, 16] used the softening model based on Rankine plasticity criterion to simulate the dynamic fracture process of concrete. In this paper, the softening model applied is a function of fracture energy, plastic strain, and the grid size, and it can limit the results of the calculation of grid sensitivity to a certain extent. In experimental study, Rossmanith et al. [17] proved that the failure mechanism of PMMA is similar to that of rock under dynamic loading. The PMMA possesses high strength and transmittance. Its fracture modes can be easily observed. The dynamic fracture mechanism of it has been widely studied by researchers [18–20].A numerical code, AUTODYN [21, 22], is applied in this study. It is an explicit finite difference code for solving a variety of nonlinear problems in solid, fluid, and gas dynamics. AUTODYN code has been successfully applied in the study of rock fracturing by Zhu et al. [9–12]. Based on the computational model proposed by Zhu et al. [9, 10] and the softening damage model for tensile cracking based on Rankine plasticity criterion, a crack propagation model of crack propagation under blasting loads was established to simulate the crack initiation and propagation and the impact of the orientation angle between the borehole and the preexisting crack. PMMA plates with an existing preset crack were used in blasting experiment, where the charge structures and the distance between the borehole and the crack tip were kept as constants, while the orientation angle between the borehole and the crack varied. The results from blasting experiment and those from the simulation were compared.


## 2. SDMTC (Softening Damage Model for Tensile Cracking)

Before the description of SDMTC (softening damage model for tensile cracking), we need to introduce JMFDFV (the joint method of finite difference and finite volume) for the SDMTC employed in this paper.


### 2.1. Brief Introduction of JMFDFV

The JMFDFV introduced in this paper is a numerical method based on the Lagrangian description and can be used to solve solid kinetic problems. Figure1 is a set of meshes for rectangular element cell under dynamic loading. 1~4 denote cell number; A
~
F denote node number. The quantities, pressure (P), density (ρ), stress (σ), strain (ε), strain rate (ε
˙), temperature (T), and mass (m), are defined within the element cell and remain unchanged in a time cycle. The quantities, coordinate (
x
,
y
), displacement (
u
,
y
), speed (u
˙
,
v
˙), acceleration (u
¨
,
v
¨), and node force (F
x
,
F
y), are defined on nodes. According to Newton’s second law of motion, the acceleration components of node A, u
¨ and v
¨, can be written as(1)
u
¨
=
F
x
m
,
v
¨
=
F
y
m
,where m denotes the total mass of domain B
D
E
F. The node speed at time n
+
1
/
2 can be written as(2)
u
˙
n
+
1
/
2
=
u
˙
n
-
1
/
2
+
u
¨
n
Δ
t
,
v
˙
n
+
1
/
2
=
v
˙
n
-
1
/
2
+
v
¨
n
Δ
t
,where n denotes the number of time step and Δ
t denotes time step. The node displacement at time n
+
1 can be written as(3)
u
n
+
1
=
u
n
+
u
˙
n
+
1
/
2
Δ
t
,
v
n
+
1
=
v
n
+
v
˙
n
+
1
/
2
Δ
t
.Figure 1
Mesh for rectangular element cell.By making use of Green’s integral formula and the four node speeds of elementA
B
C
D, the strain rates of element A
B
C
D, ε
˙
x, ε
˙
y, and γ
˙
x
y, can be written as(4)
ε
˙
x
=
∂
u
˙
∂
x
=
1
2
A
E
n
+
1
/
2
u
˙
A
-
u
˙
C
y
B
-
y
D
-
u
˙
B
-
u
˙
D
y
A
-
y
C
,
ε
˙
y
=
∂
v
˙
∂
y
=
1
2
A
E
n
+
1
/
2
v
˙
A
-
v
˙
C
x
B
-
x
D
-
v
˙
B
-
v
˙
D
x
A
-
x
C
,
γ
˙
x
y
=
∂
u
˙
∂
y
+
∂
v
˙
∂
x
=
1
2
A
E
n
+
1
/
2
v
˙
A
-
v
˙
C
y
B
-
y
D
-
v
˙
B
-
v
˙
D
y
A
-
y
C
+
u
˙
A
-
u
˙
C
x
B
-
x
D
-
u
˙
B
-
u
˙
D
y
A
-
y
C
,where A
E denotes the area of element A
B
C
D.Using the strain rates in (4), the stress deviators, S
x, S
y, and S
x
y, in step n
+
1 can be obtained as(5)
S
x
n
+
1
=
S
x
n
+
2
G
ε
˙
x
n
+
1
/
2
-
1
3
e
˙
n
+
1
/
2
Δ
t
,
S
y
n
+
1
=
S
y
n
+
2
G
ε
˙
y
n
+
1
/
2
-
1
3
e
˙
n
+
1
/
2
Δ
t
,
S
x
y
n
+
1
=
S
x
y
n
+
2
G
γ
˙
x
y
n
+
1
/
2
Δ
t
,where G denotes shear modulus and e
˙ volume strain rate. In 2D cases, e
˙
=
ε
˙
x
+
ε
˙
y.The stresses of elementA
B
C
D can be written as(6)
σ
x
=
p
+
S
x
,
σ
y
=
p
+
S
y
,
τ
x
y
=
S
x
y
,where p denotes mean pressure and is determined by the material EOS in the model.The node force is determined by integral on both ends of the equations of wave propagation:(7)
∬
S
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
d
x
d
y
=
∬
S
ρ
∂
2
u
∂
t
2
d
x
d
y
,
∬
S
∂
σ
y
∂
y
+
∂
τ
x
y
∂
x
d
x
d
y
=
∬
S
ρ
∂
2
v
∂
t
2
d
x
d
y
.According to Green’s integral formula, (7) can be rewritten as(8)
∬
S
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
d
x
d
y
=
∫
L
σ
x
d
y
-
τ
x
y
d
x
,
∬
S
∂
σ
y
∂
y
+
∂
τ
x
y
∂
x
d
x
d
y
=
∫
L
τ
x
y
d
x
-
σ
y
d
y
,where S denotes the area of domain B
D
E
F and L is the closed curve B
D
E
F. The following equations can be obtained using (7) and (8):(9)
∫
L
σ
x
d
y
-
τ
x
y
d
x
=
∬
S
ρ
∂
2
u
∂
t
2
d
x
d
y
=
2
m
ρ
∂
2
u
∂
t
2
=
2
F
x
,
∫
L
τ
x
y
d
x
-
σ
y
d
y
=
∬
S
ρ
∂
2
v
∂
t
2
d
x
d
y
=
2
m
ρ
∂
2
v
∂
t
2
=
2
F
y
.Taking into consideration (9) and the mesh in Figure 1, the following finite difference equations for 2D can be obtained as(10)
F
x
=
1
2
σ
x
1
y
F
-
y
E
-
τ
x
y
1
x
F
-
x
E
+
σ
x
2
y
B
-
y
F
-
τ
x
y
2
x
B
-
x
F
+
σ
x
3
y
D
-
y
B
-
τ
x
y
3
x
D
-
x
B
+
σ
x
4
y
E
-
y
D
-
τ
x
y
4
x
E
-
x
D
,
F
y
=
1
2
τ
x
y
1
y
F
-
y
E
-
σ
y
1
x
F
-
x
E
+
τ
x
y
2
y
B
-
y
F
-
σ
y
2
x
B
-
x
F
+
τ
x
y
3
y
D
-
y
B
-
σ
y
3
x
D
-
x
B
+
τ
x
y
4
y
E
-
y
D
-
σ
y
4
x
E
-
x
D
,where superscript 1, 2, 3, and 4 denote elements 1, 2, 3, and 4 in Figure 1, respectively.In practice, besides the type of forces in (10), there are other nodal forces that should be considered. They are pseudoviscous forces, hourglass damping, and antitangle forces, and if the node is located on the boundaries of the load, external forces should be incorporated. The computational method adopted in this paper used the numerical method for pseudoviscous force proposed by Wilkins [23]. Adopt the method proposed by Century Dynamics Inc. [22] and the hourglass damping is determined, namely, adding a set of corrective forces in the four nodes of the cell. The external forces are divided equally between the two nodes on the boundary, and the boundary forces due to external forces are therefore determined. Above all, the nodal force is a resultant force. The process starting from (1) to (10) is a cycle of the computation. The nodal acceleration is calculated after calculating the nodal force. By iteration of this calculation process, the quantities, node speed, displacement, acceleration, element stress, and strain rate of each node and element in different times can be obtained.


### 2.2. SDMTC

Basic physical quantities can be obtained by the computational method introduced above. To describe the failure process in reality, however, an effective failure criterion is needed. Thus MMTS (the modified maximum tensile stress) is adopted in this paper. The MMTS is an instantaneous failure law. The meaning can be expressed as follows: element loses its ability to bear tensile and shear stresses, when the principal stress of the calculating element reaches the tensile strength or when the maximum shear stress reaches the shear strength of the material, but its ability to bear compression is maintained. The law can be written as(11)
σ
1
≤
σ
T
;
τ
max
⁡
=
σ
1
-
σ
3
2
≤
τ
c
,where σ
1, σ
3 denote two principal stresses, τ
max
⁡ denotes the maximum shear stress, and σ
T, τ
c denote the tensile strength and shear strength of material, respectively.In this paper, stress, strain, and strain rate are defined and kept unchanged within element in a time cycle. The numerical method using MMTS can easily lead to the dependence of calculation results on grid. In fact, element loses its bearing capacity gradually as crack propagates through it. Figure2 sheds light on the gradual failure process of an element under static loading; that is, σ
max
⁡ increases to a certain value that is larger than or equal to the tensile strength σ
T before it decreases to a value smaller than σ
T and reaches 0. According to the SIF (stress intensity factor) of edge crack and the judgment formula for critical stress in LEFM (linear elastic fracture mechanics), the relationship between critical stress of element and crack length a can be listed as(12)
K
=
1.12
σ
max
⁡
π
a
=
K
IC
,
σ
max
⁡
=
K
IC
1.12
π
a
=
E
G
c
1.12
π
a
,where K
IC denotes fracture toughness, G
c fracture energy, E elastic modulus, and a crack length.Figure 2
The maximum tensile stress varies with crack propagation.
(a)
Intact
(b)
Fractured
(c)
FailureIn order to illustrate the process of gradual failure for an element, the author adopted the linear SDMTC algorithm and Figure3 shows the stress-strain curve of this method.Figure 3
Relationship of stress-strain for SDMTC.In this model, the MMTS determines the initiation stage of element failure, after which the maximum tensile principal stressσ
max
⁡, crack strain ε
cr, and the crack strain ε
u for the totally fractured element are recorded, respectively. ε
u is defined as a function of G
c, dynamic tensile strength σ
T, and the critical length L in the direction of maximum tensile principal stress:(13)
ε
u
=
2
G
c
σ
T
L
.The damage variable Dam is used to describe the damage level, defined as a function ofε
cr and ε
u, where ε
cr is obtained by the Backward-Euler method:(14)
Dam
=
ε
cr
ε
u
=
ε
cr
σ
T
L
2
G
c
.In damage process, the maximum tensile principal stress is defined as a function of dynamic tensile strengthσ
T and damage variable Dam:(15)
σ
max
⁡
=
σ
T
1
-
Dam
.


### 2.3. The Calculation of Crack Strain

The calculation of crack strain is the key to the realization of SDMTC. The crack strain calculated by Backward-Euler method must reflect and incorporate the essence of material failure. Brittle damage is tensile failure developed along the direction of principal stress. Hence, Rankine’s plastic yield criterion is chosen:(16)
f
1
=
σ
1
-
σ
T
~
κ
T
=
σ
1
-
σ
3
2
+
σ
1
+
σ
3
2
-
σ
T
~
κ
T
,where σ
1, σ
3 denote the principal stresses and σ
T
~
(
κ
T
) denotes the equivalent stress.This paper uses Backward-Euler method to obtain the crack strain. When the trial stress surpasses the yield surface, associated Backward-Euler method is employed to back-flow the stress to the yield surface. This method is similar to that described by Crisfield [24] showed in Figure 4. Associated flow is used in the meridian plane and π-plane.Figure 4
Sketch of associated flow.


## 2.1. Brief Introduction of JMFDFV

The JMFDFV introduced in this paper is a numerical method based on the Lagrangian description and can be used to solve solid kinetic problems. Figure1 is a set of meshes for rectangular element cell under dynamic loading. 1~4 denote cell number; A
~
F denote node number. The quantities, pressure (P), density (ρ), stress (σ), strain (ε), strain rate (ε
˙), temperature (T), and mass (m), are defined within the element cell and remain unchanged in a time cycle. The quantities, coordinate (
x
,
y
), displacement (
u
,
y
), speed (u
˙
,
v
˙), acceleration (u
¨
,
v
¨), and node force (F
x
,
F
y), are defined on nodes. According to Newton’s second law of motion, the acceleration components of node A, u
¨ and v
¨, can be written as(1)
u
¨
=
F
x
m
,
v
¨
=
F
y
m
,where m denotes the total mass of domain B
D
E
F. The node speed at time n
+
1
/
2 can be written as(2)
u
˙
n
+
1
/
2
=
u
˙
n
-
1
/
2
+
u
¨
n
Δ
t
,
v
˙
n
+
1
/
2
=
v
˙
n
-
1
/
2
+
v
¨
n
Δ
t
,where n denotes the number of time step and Δ
t denotes time step. The node displacement at time n
+
1 can be written as(3)
u
n
+
1
=
u
n
+
u
˙
n
+
1
/
2
Δ
t
,
v
n
+
1
=
v
n
+
v
˙
n
+
1
/
2
Δ
t
.Figure 1
Mesh for rectangular element cell.By making use of Green’s integral formula and the four node speeds of elementA
B
C
D, the strain rates of element A
B
C
D, ε
˙
x, ε
˙
y, and γ
˙
x
y, can be written as(4)
ε
˙
x
=
∂
u
˙
∂
x
=
1
2
A
E
n
+
1
/
2
u
˙
A
-
u
˙
C
y
B
-
y
D
-
u
˙
B
-
u
˙
D
y
A
-
y
C
,
ε
˙
y
=
∂
v
˙
∂
y
=
1
2
A
E
n
+
1
/
2
v
˙
A
-
v
˙
C
x
B
-
x
D
-
v
˙
B
-
v
˙
D
x
A
-
x
C
,
γ
˙
x
y
=
∂
u
˙
∂
y
+
∂
v
˙
∂
x
=
1
2
A
E
n
+
1
/
2
v
˙
A
-
v
˙
C
y
B
-
y
D
-
v
˙
B
-
v
˙
D
y
A
-
y
C
+
u
˙
A
-
u
˙
C
x
B
-
x
D
-
u
˙
B
-
u
˙
D
y
A
-
y
C
,where A
E denotes the area of element A
B
C
D.Using the strain rates in (4), the stress deviators, S
x, S
y, and S
x
y, in step n
+
1 can be obtained as(5)
S
x
n
+
1
=
S
x
n
+
2
G
ε
˙
x
n
+
1
/
2
-
1
3
e
˙
n
+
1
/
2
Δ
t
,
S
y
n
+
1
=
S
y
n
+
2
G
ε
˙
y
n
+
1
/
2
-
1
3
e
˙
n
+
1
/
2
Δ
t
,
S
x
y
n
+
1
=
S
x
y
n
+
2
G
γ
˙
x
y
n
+
1
/
2
Δ
t
,where G denotes shear modulus and e
˙ volume strain rate. In 2D cases, e
˙
=
ε
˙
x
+
ε
˙
y.The stresses of elementA
B
C
D can be written as(6)
σ
x
=
p
+
S
x
,
σ
y
=
p
+
S
y
,
τ
x
y
=
S
x
y
,where p denotes mean pressure and is determined by the material EOS in the model.The node force is determined by integral on both ends of the equations of wave propagation:(7)
∬
S
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
d
x
d
y
=
∬
S
ρ
∂
2
u
∂
t
2
d
x
d
y
,
∬
S
∂
σ
y
∂
y
+
∂
τ
x
y
∂
x
d
x
d
y
=
∬
S
ρ
∂
2
v
∂
t
2
d
x
d
y
.According to Green’s integral formula, (7) can be rewritten as(8)
∬
S
∂
σ
x
∂
x
+
∂
τ
x
y
∂
y
d
x
d
y
=
∫
L
σ
x
d
y
-
τ
x
y
d
x
,
∬
S
∂
σ
y
∂
y
+
∂
τ
x
y
∂
x
d
x
d
y
=
∫
L
τ
x
y
d
x
-
σ
y
d
y
,where S denotes the area of domain B
D
E
F and L is the closed curve B
D
E
F. The following equations can be obtained using (7) and (8):(9)
∫
L
σ
x
d
y
-
τ
x
y
d
x
=
∬
S
ρ
∂
2
u
∂
t
2
d
x
d
y
=
2
m
ρ
∂
2
u
∂
t
2
=
2
F
x
,
∫
L
τ
x
y
d
x
-
σ
y
d
y
=
∬
S
ρ
∂
2
v
∂
t
2
d
x
d
y
=
2
m
ρ
∂
2
v
∂
t
2
=
2
F
y
.Taking into consideration (9) and the mesh in Figure 1, the following finite difference equations for 2D can be obtained as(10)
F
x
=
1
2
σ
x
1
y
F
-
y
E
-
τ
x
y
1
x
F
-
x
E
+
σ
x
2
y
B
-
y
F
-
τ
x
y
2
x
B
-
x
F
+
σ
x
3
y
D
-
y
B
-
τ
x
y
3
x
D
-
x
B
+
σ
x
4
y
E
-
y
D
-
τ
x
y
4
x
E
-
x
D
,
F
y
=
1
2
τ
x
y
1
y
F
-
y
E
-
σ
y
1
x
F
-
x
E
+
τ
x
y
2
y
B
-
y
F
-
σ
y
2
x
B
-
x
F
+
τ
x
y
3
y
D
-
y
B
-
σ
y
3
x
D
-
x
B
+
τ
x
y
4
y
E
-
y
D
-
σ
y
4
x
E
-
x
D
,where superscript 1, 2, 3, and 4 denote elements 1, 2, 3, and 4 in Figure 1, respectively.In practice, besides the type of forces in (10), there are other nodal forces that should be considered. They are pseudoviscous forces, hourglass damping, and antitangle forces, and if the node is located on the boundaries of the load, external forces should be incorporated. The computational method adopted in this paper used the numerical method for pseudoviscous force proposed by Wilkins [23]. Adopt the method proposed by Century Dynamics Inc. [22] and the hourglass damping is determined, namely, adding a set of corrective forces in the four nodes of the cell. The external forces are divided equally between the two nodes on the boundary, and the boundary forces due to external forces are therefore determined. Above all, the nodal force is a resultant force. The process starting from (1) to (10) is a cycle of the computation. The nodal acceleration is calculated after calculating the nodal force. By iteration of this calculation process, the quantities, node speed, displacement, acceleration, element stress, and strain rate of each node and element in different times can be obtained.


## 2.2. SDMTC

Basic physical quantities can be obtained by the computational method introduced above. To describe the failure process in reality, however, an effective failure criterion is needed. Thus MMTS (the modified maximum tensile stress) is adopted in this paper. The MMTS is an instantaneous failure law. The meaning can be expressed as follows: element loses its ability to bear tensile and shear stresses, when the principal stress of the calculating element reaches the tensile strength or when the maximum shear stress reaches the shear strength of the material, but its ability to bear compression is maintained. The law can be written as(11)
σ
1
≤
σ
T
;
τ
max
⁡
=
σ
1
-
σ
3
2
≤
τ
c
,where σ
1, σ
3 denote two principal stresses, τ
max
⁡ denotes the maximum shear stress, and σ
T, τ
c denote the tensile strength and shear strength of material, respectively.In this paper, stress, strain, and strain rate are defined and kept unchanged within element in a time cycle. The numerical method using MMTS can easily lead to the dependence of calculation results on grid. In fact, element loses its bearing capacity gradually as crack propagates through it. Figure2 sheds light on the gradual failure process of an element under static loading; that is, σ
max
⁡ increases to a certain value that is larger than or equal to the tensile strength σ
T before it decreases to a value smaller than σ
T and reaches 0. According to the SIF (stress intensity factor) of edge crack and the judgment formula for critical stress in LEFM (linear elastic fracture mechanics), the relationship between critical stress of element and crack length a can be listed as(12)
K
=
1.12
σ
max
⁡
π
a
=
K
IC
,
σ
max
⁡
=
K
IC
1.12
π
a
=
E
G
c
1.12
π
a
,where K
IC denotes fracture toughness, G
c fracture energy, E elastic modulus, and a crack length.Figure 2
The maximum tensile stress varies with crack propagation.
(a)
Intact
(b)
Fractured
(c)
FailureIn order to illustrate the process of gradual failure for an element, the author adopted the linear SDMTC algorithm and Figure3 shows the stress-strain curve of this method.Figure 3
Relationship of stress-strain for SDMTC.In this model, the MMTS determines the initiation stage of element failure, after which the maximum tensile principal stressσ
max
⁡, crack strain ε
cr, and the crack strain ε
u for the totally fractured element are recorded, respectively. ε
u is defined as a function of G
c, dynamic tensile strength σ
T, and the critical length L in the direction of maximum tensile principal stress:(13)
ε
u
=
2
G
c
σ
T
L
.The damage variable Dam is used to describe the damage level, defined as a function ofε
cr and ε
u, where ε
cr is obtained by the Backward-Euler method:(14)
Dam
=
ε
cr
ε
u
=
ε
cr
σ
T
L
2
G
c
.In damage process, the maximum tensile principal stress is defined as a function of dynamic tensile strengthσ
T and damage variable Dam:(15)
σ
max
⁡
=
σ
T
1
-
Dam
.


## 2.3. The Calculation of Crack Strain

The calculation of crack strain is the key to the realization of SDMTC. The crack strain calculated by Backward-Euler method must reflect and incorporate the essence of material failure. Brittle damage is tensile failure developed along the direction of principal stress. Hence, Rankine’s plastic yield criterion is chosen:(16)
f
1
=
σ
1
-
σ
T
~
κ
T
=
σ
1
-
σ
3
2
+
σ
1
+
σ
3
2
-
σ
T
~
κ
T
,where σ
1, σ
3 denote the principal stresses and σ
T
~
(
κ
T
) denotes the equivalent stress.This paper uses Backward-Euler method to obtain the crack strain. When the trial stress surpasses the yield surface, associated Backward-Euler method is employed to back-flow the stress to the yield surface. This method is similar to that described by Crisfield [24] showed in Figure 4. Associated flow is used in the meridian plane and π-plane.Figure 4
Sketch of associated flow.


## 3. Numerical Blasting Model

Using the numerical method introduced above and PMMA parameters obtained by experiment, the 2D numerical blasting model shown in Figure5 is established.Figure 5
Schematic of blasting models, wherea is the crack length, D is borehole diameter, d is the ANFO diameter, r is the distance from borehole center to crack tip, and β ranges from 0° to 75°, the units being mm.


### 3.1. Blasting Model

In the numerical blasting model shown in Figure5, the decoupling charge is adopted, decoupling coefficient being set as 2. Dimension of the model is 400 mm × 400 mm, with one borehole and one thorough centric crack. To delve into the impact of the orientation angle between crack and borehole on crack propagation, a set of models, with orientation angle being 0°, 15°, 30°, 45°, 60°, or 75°, are established. Free boundary is chosen as the boundary condition. The charge and PMMA are set with different subgrids, and the mesh near the borehole is shown in Figure 6. The mesh size determines the computational efficiency and the accuracy of result. However, it is difficult to estimate the sensitivity of our approach. Usually, we determine the best mesh size by trying different mesh sizes, like the method used by Y. Hao and H. Hao [25]. The total number of elements is more than 160000 and the maximum mesh size is less than 2 mm. These meshes are seen to be fine enough for the required precision of calculation. With additional details on the interactions between subgrids, such as gap, one type of Lagrange-Lagrange coupling can be found in [21, 22]. No mesh is set for air for the inaccessibility of computation caused by the large blasting-induced deformation.Figure 6
Mesh of the explosive and the PMMA.


### 3.2. Material Model

PMMA falls into the catalog of isotropic material, has been widely studied, shows similar characteristics to rock under dynamic loading, machines easily, and has excellent transparency for crack observation. In accordance with the charge amount, the linear EOS is adopted for PMMA [22]:(17)
p
=
k
ρ
ρ
0
-
1
,where p denotes pressure, k bulk modulus, and ρ
/
ρ
0 the ratio of present density and initial density. The JMFDFV and MMTS are used as the failure models for PMMA with the parameters listed in Table 1.Table 1
Dynamic parameters of PMMA [26].
ρ (g·cm−3)
E
d (Gpa)
K
d (Gpa)
G
d (Gpa)
C
p (m·s−1)
C
s (m·s−1)
ν
1.187
6.1
5.35
2.328
2320
1260
0.31The dynamic strength of PMMA is difficult to determine, because it changes largely at different environment temperature and strain rate [27]. The shear strength is determined from the charge amount and experiment results. The static values [28] of tensile strength and GC were adopted and the tensile strength was modified experientially because later experiment was taken at temperature 313 K. The dynamic strength parameters of PMMA are listed in Table 2.Table 2
Dynamic strength parameters of PMMA.
σ
T (Mpa)
τ
c (Mpa)
G
c (J/m2)
30
58
133The JWL (Johns-Wilkins-Lee) EOS (equation of state) has been adopted for the ANFO charge:(18)
p
=
A
1
-
ω
R
1
V
e
-
R
1
V
+
B
1
-
ω
R
2
V
e
-
R
2
V
+
ω
e
V
,where p denotes pressure, e denotes the initial internal energy of denotation product, V denotes relative volume, and A, B, R
1, R
2, and ω denote material constants. With strength model and failure criterion ignored, the charge will turn into ideal gas after blasting. The material parameters for ANFO (ammonium nitrate/fuel oil) are listed in Table 3.Table 3
Parameters of ANFO [29].
ρ (g·cm−3)
D (m·s−1)
p
c
j (Gpa)
A (Gpa)
B (Gpa)
R
1
R
2
ω
0.931
4160
5.15
49.46
1.891
3.907
1.18
0.333


## 3.1. Blasting Model

In the numerical blasting model shown in Figure5, the decoupling charge is adopted, decoupling coefficient being set as 2. Dimension of the model is 400 mm × 400 mm, with one borehole and one thorough centric crack. To delve into the impact of the orientation angle between crack and borehole on crack propagation, a set of models, with orientation angle being 0°, 15°, 30°, 45°, 60°, or 75°, are established. Free boundary is chosen as the boundary condition. The charge and PMMA are set with different subgrids, and the mesh near the borehole is shown in Figure 6. The mesh size determines the computational efficiency and the accuracy of result. However, it is difficult to estimate the sensitivity of our approach. Usually, we determine the best mesh size by trying different mesh sizes, like the method used by Y. Hao and H. Hao [25]. The total number of elements is more than 160000 and the maximum mesh size is less than 2 mm. These meshes are seen to be fine enough for the required precision of calculation. With additional details on the interactions between subgrids, such as gap, one type of Lagrange-Lagrange coupling can be found in [21, 22]. No mesh is set for air for the inaccessibility of computation caused by the large blasting-induced deformation.Figure 6
Mesh of the explosive and the PMMA.


## 3.2. Material Model

PMMA falls into the catalog of isotropic material, has been widely studied, shows similar characteristics to rock under dynamic loading, machines easily, and has excellent transparency for crack observation. In accordance with the charge amount, the linear EOS is adopted for PMMA [22]:(17)
p
=
k
ρ
ρ
0
-
1
,where p denotes pressure, k bulk modulus, and ρ
/
ρ
0 the ratio of present density and initial density. The JMFDFV and MMTS are used as the failure models for PMMA with the parameters listed in Table 1.Table 1
Dynamic parameters of PMMA [26].
ρ (g·cm−3)
E
d (Gpa)
K
d (Gpa)
G
d (Gpa)
C
p (m·s−1)
C
s (m·s−1)
ν
1.187
6.1
5.35
2.328
2320
1260
0.31The dynamic strength of PMMA is difficult to determine, because it changes largely at different environment temperature and strain rate [27]. The shear strength is determined from the charge amount and experiment results. The static values [28] of tensile strength and GC were adopted and the tensile strength was modified experientially because later experiment was taken at temperature 313 K. The dynamic strength parameters of PMMA are listed in Table 2.Table 2
Dynamic strength parameters of PMMA.
σ
T (Mpa)
τ
c (Mpa)
G
c (J/m2)
30
58
133The JWL (Johns-Wilkins-Lee) EOS (equation of state) has been adopted for the ANFO charge:(18)
p
=
A
1
-
ω
R
1
V
e
-
R
1
V
+
B
1
-
ω
R
2
V
e
-
R
2
V
+
ω
e
V
,where p denotes pressure, e denotes the initial internal energy of denotation product, V denotes relative volume, and A, B, R
1, R
2, and ω denote material constants. With strength model and failure criterion ignored, the charge will turn into ideal gas after blasting. The material parameters for ANFO (ammonium nitrate/fuel oil) are listed in Table 3.Table 3
Parameters of ANFO [29].
ρ (g·cm−3)
D (m·s−1)
p
c
j (Gpa)
A (Gpa)
B (Gpa)
R
1
R
2
ω
0.931
4160
5.15
49.46
1.891
3.907
1.18
0.333


## 4. Simulation Results and Analysis

It is well known that the blasting load in the borehole is p wave in 2D case and p wave scattering near the crack is shown in Figure7. Crack blocks the spread of stress wave, diffracted wave generates at the crack tip, and reflected wave generates on the crack surface. As you can image, when the incident wave impinges the crack with a certain angle, the crack tip nearer to the borehole cannot completely block the incident wave. At the crack tip on the far end, there exists only the diffraction wave on the side of the crack surface due to the crack surface blocking.Figure 7
p wave scattering near the crack.


### 4.1. Crack Initiation at the Far End for Incident Angular Degree 0°

Only a single gauge point is set in the vicinity of crack tip, gauge #3 shown in Figure5, to analyze the cause for the initiation at far end, because the loading is symmetrical for the case when the incident angle is 0°. Figure 8(a) shows the time history of the major principal stress. It can be seen that the major principal stress turns from compression to tensile stress, and as it reaches the tensile strength, material is damaged and no tensile stress is sustained. This is because when the incident angle is 0°, neither reflection nor diffraction is aroused by stress wave. The driving force for crack initiation and propagation is the same as the case without a crack, that is, the circular tensile stress, which is generated by the radial movement of materials, result of compression pressure on the borehole. The circular tensile stress gives rise to crack initiation and propagation.Figure 8
Time history curves of major principal stress at gauges for incident angle of (a) 0° and (b) 45°.
(a)
(b)


### 4.2. Crack Initiation at the Far End for Incident Angle 45°

Two gauge points are set in the vicinity of crack tip, gauge #1 and gauge #2 shown in Figure5, to analyze the cause for the initiation at far end, one being in the domain where the crack surface faces the borehole and the other on the other side of the crack tip. Figure 8(b) shows the time history of major principal stress for the two gauge points. In this figure, gauge #1 starts with compression stress, while gauge #2 starts with tensile stress. When the detonation-induced compression stress encounters the crack tip, diffraction wave generates. The general wave field consists of both incident wave and diffraction wave. The domain in the crack back is in a stress-free state; the only acting wave is diffraction wave composed of tensile wave and shear wave. Hence, the crack initiation at the far end begins in the direction away from the borehole.


### 4.3. The Impact of Crack Orientation on Crack Initiation at Far End

Based on the initial analysis of the crack tip at far end for incident angular degrees 0° and 45°, two conclusions are reached. First, when the incident angular degree is relatively small, the cause for the crack initiation at far end is mainly because of the circular tensile stress. While the incident angular degree grows, the diffraction of stress wave plays the major role. Figure9 compares the damage status of the models, in angular degree range 0–75°, 80 μs after detonation. As the incident angle increases, the initiation angle at far end increases. When the angular degree reaches 75°, crack initiation takes form at a distance from the preset crack tip.Figure 9
Material status of rock samples as a function of degreeβ att = 80 us after detonation in the borehole.


### 4.4. Comparison of Experiment and Simulation Results

In order to verify the validation of the results from numerical simulation, corresponding blasting experiment was carried out in the blasting laboratory, Southwest University of Science and Technology, using PMMA. The environment temperature is 313 K. Specimen size is400
∗
400
∗
8 mm and ANFO is employed as the explosive. The experiment results were listed in Figure 10. This figure is part of the specimen and it was taken after explosion. Dark paper was used as the background paper for photographing (we put the specimen on a dark paper, and then we take it using camera, Canon power shot A3300IS, vertical to specimen.). It can be seen that the angular degree for crack initiation at far end is less than 90° when crack angle falls in 0°–30° range. When crack angle stays in range 45°–75°, the angular degree for crack initiation is a bit larger than 90°. The experiment results illustrated in Figure 9 are analogous to the simulation results in Figure 10, meaning the simulation is reliable.Figure 10
Experiment results of different incident angles in PMMA.


## 4.1. Crack Initiation at the Far End for Incident Angular Degree 0°

Only a single gauge point is set in the vicinity of crack tip, gauge #3 shown in Figure5, to analyze the cause for the initiation at far end, because the loading is symmetrical for the case when the incident angle is 0°. Figure 8(a) shows the time history of the major principal stress. It can be seen that the major principal stress turns from compression to tensile stress, and as it reaches the tensile strength, material is damaged and no tensile stress is sustained. This is because when the incident angle is 0°, neither reflection nor diffraction is aroused by stress wave. The driving force for crack initiation and propagation is the same as the case without a crack, that is, the circular tensile stress, which is generated by the radial movement of materials, result of compression pressure on the borehole. The circular tensile stress gives rise to crack initiation and propagation.Figure 8
Time history curves of major principal stress at gauges for incident angle of (a) 0° and (b) 45°.
(a)
(b)


## 4.2. Crack Initiation at the Far End for Incident Angle 45°

Two gauge points are set in the vicinity of crack tip, gauge #1 and gauge #2 shown in Figure5, to analyze the cause for the initiation at far end, one being in the domain where the crack surface faces the borehole and the other on the other side of the crack tip. Figure 8(b) shows the time history of major principal stress for the two gauge points. In this figure, gauge #1 starts with compression stress, while gauge #2 starts with tensile stress. When the detonation-induced compression stress encounters the crack tip, diffraction wave generates. The general wave field consists of both incident wave and diffraction wave. The domain in the crack back is in a stress-free state; the only acting wave is diffraction wave composed of tensile wave and shear wave. Hence, the crack initiation at the far end begins in the direction away from the borehole.


## 4.3. The Impact of Crack Orientation on Crack Initiation at Far End

Based on the initial analysis of the crack tip at far end for incident angular degrees 0° and 45°, two conclusions are reached. First, when the incident angular degree is relatively small, the cause for the crack initiation at far end is mainly because of the circular tensile stress. While the incident angular degree grows, the diffraction of stress wave plays the major role. Figure9 compares the damage status of the models, in angular degree range 0–75°, 80 μs after detonation. As the incident angle increases, the initiation angle at far end increases. When the angular degree reaches 75°, crack initiation takes form at a distance from the preset crack tip.Figure 9
Material status of rock samples as a function of degreeβ att = 80 us after detonation in the borehole.


## 4.4. Comparison of Experiment and Simulation Results

In order to verify the validation of the results from numerical simulation, corresponding blasting experiment was carried out in the blasting laboratory, Southwest University of Science and Technology, using PMMA. The environment temperature is 313 K. Specimen size is400
∗
400
∗
8 mm and ANFO is employed as the explosive. The experiment results were listed in Figure 10. This figure is part of the specimen and it was taken after explosion. Dark paper was used as the background paper for photographing (we put the specimen on a dark paper, and then we take it using camera, Canon power shot A3300IS, vertical to specimen.). It can be seen that the angular degree for crack initiation at far end is less than 90° when crack angle falls in 0°–30° range. When crack angle stays in range 45°–75°, the angular degree for crack initiation is a bit larger than 90°. The experiment results illustrated in Figure 9 are analogous to the simulation results in Figure 10, meaning the simulation is reliable.Figure 10
Experiment results of different incident angles in PMMA.


## 5. Conclusions

The crack propagation of rock mass with a preexisting crack under blasting load was studied by using numerical models. The corresponding confirmatory experiment was carried out using PMMA and the following conclusions can be obtained:(1)
When the incident angle is 0°, the cracks initiate and propagate at both tips, and at the far tip, they are driven by the circular tensile stress.
(2)
When the incident angle is 45°, the crack initiation and propagation at the far end are driven by the diffractive tensile stress in the crack back.
(3)
As the incident angle increases, the driving force for crack initiation and propagation at the far end transform from circular tensile stress to diffractive tensile stress. The experiment results agreed with the simulation results.



---

*Source: 108580-2015-08-02.xml*