Computers and Structures 287 (2023) 107116

Contents lists available at ScienceDirect

Computers and Structures

journal homepage: www.elsevier.com/locate/compstruc

Quantitative assessment of the effect of in-situ stresses on blast-induced

damage to rock
Guodong Qiao a, Zegong Liu a, *, Changping Yi b, **, Kui Gao a, Gaoyuan Xuan c
a
School of Safety Science and Engineering, Anhui University of Science and Technology, Huainan, Anhui 232001, China
b
Division of Mining and Geotechnical Engineering, Luleå University of Technology, Luleå, Sweden
c
School of Artificial Intelligence, Anhui University of Science and Technology, Huainan, Anhui 232001, China

A R T I C L E I N F O A B S T R A C T

Keywords: In-situ stress significantly affects rock blast damage but there is a paucity of quantitative assessments of damage
Blasting evolution in rocks affected by confining pressure. The present paper analyses the effect of envelope pressure on
In-situ stress blast-induced rock damage through theoretical analysis and numerical simulations. Damage clouds obtained
Damage evolution
from numerical simulations are processed using image processing techniques. The concept of the damage vari­
Image processing
Quantitative assessment
able (η) is proposed to facilitate the presentation of the image processing results. The damage variable is found to
be negatively correlated with the hydrostatic pressure (Px ) at the same moment, in equiaxial in-situ stress fields.
In contrast, in anisotropic in-situ stress fields, η is not negatively correlated with Px due to the presence of hoop
tensile stresses in the rock. The mathematical relationship between η and Px in equiaxial and anisotropic stress
fields are established. An anisotropic damage variable (ηk ) is introduced to describe the effect of the anisotropy
ratio (K) on rock damage, which is found to increase with increasing values of K. The sharp increase in K equal to
4 and 5 is explained in terms of the state of the rock stress distribution under static loading. This study provides
insights into the effect of in situ stress on rock blast damage and presents new approaches for analyzing and
presenting the data.

1. Introduction Nicholls and Duvall [8] observed that the pre-cracked rocks under static
stress fields are most likely to fracture in the direction of the maximum
With the depletion of shallow mineral resources, the mining industry compressive principal stress. Yang et al. [9] found through photoelastic
is turning to deep mining to extract all kinds of minerals. Many mines experiments that the in-situ stress perpendicular to the crack would
around the world now exceed depths of 2000 m [1]. The TauTona gold reduce the stress intensity factor at the crack tip, thus hindering crack
mine in South Africa is mined to a depth of 3,900 m with an in-situ stress growth. Zhang et al. [10] investigated the impact of confining pressure
level of about 100 MPa [2]. Other deep mines such as the Kloof gold on the shape and volume of the crushed zone using model experiments.
mine, East Rand Proprietary, and Driefontein gold mine have depths This effect demonstrates that the crushed zone initially increases and
exceeding 3,500 m [3]. In response to the increasing demand for natural then decreases with increasing confining pressure. He et al. [11] used
resources, the mining depth is increasing yearly to obtain more minerals. Digital Image Correlation (DIC) and strain measurement techniques to
Take China’s coal mines, for example, increasing by 10 ~ 25 m annually analyze the dynamic response of rock under active confining pressure
[4]. The drill and blast (D & B) technique is a widely used method for and blasting dynamic load. The results revealed that the crack propa­
mining operations, but there is unique mechanical behavior of rocks gation direction is controlled by both, circumferential tensile stress and
under high in situ stress compared with shallow rocks [5,6]. High in-situ biaxial preloading ratio. Yue et al. [12] used a digital laser dynamic
stress is becoming one of the main challenges for deep mining blasting caustic experiment system to carry out experiments, indicating that the
excavation [7]. pressure in the vertical direction hinders the crack growth, while pres­
Laboratory experiments are crucial for studying the influence of in- sure parallel to the crack direction promotes it. Also, state that when the
situ stress on the dynamic response of the rock during blasting. principal pressure is at a certain angle to the slit, the main crack

* Corresponding author.
** Corresponding author.
E-mail addresses: [email protected] (Z. Liu), [email protected] (C. Yi).

https://doi.org/10.1016/j.compstruc.2023.107116
Received 12 December 2022; Accepted 21 July 2023
Available online 31 July 2023
0045-7949/© 2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
G. Qiao et al. Computers and Structures 287 (2023) 107116

propagates in the pressure direction. Wang [13] conducted blasting basis for analyzing the distribution law of rock damage under dynamic
experiments with similar materials and showed that an increase in the blasting load. ANSYS/LS-DYNA numerical simulation software was then
maximum principal stress causes the blasting crack to deflect towards used to replicate the dynamic response of rock blasting under static
the direction of the maximum principal stress, but slows down the loads. Then, the numerical simulation images were processed using a
growth rate of the blasting crack. proposed method for quantitative assessment of blast damage areas
Numerical simulation to study the influence of in-situ stress on rock based on image processing techniques. Finally, some new useful quan­
blasting has become increasingly popular with the development of titative conclusions on blast damage patterns under the influence of in-
computer technology. Among the many numerical simulation software, situ stress magnitude and anisotropy were obtained through the analysis
ANSYS/LS-DYNA performs well-suited simulations for the nonlinear of image processing results.
dynamic response of structures [14] and has been widely used to study
the dynamic response of rock blasting under in-situ stress. Ma and An 2. Hoop stress distribution under static loads
[15] used the Johnson–Holmquist material model to simulate rock
material behavior during a blast and concluded that the blast cracks are 2.1. Hoop stress distribution at the wall of the blasthole
aligned in the direction of the prestress axis. This phenomenon becomes
more pronounced with increasing prestress. Lu et al. [16] used a coupled In engineering structural analysis, when the geometry and forces of a
SPH-FEM model and found that the radius of the crushed zone increases structure have certain characteristics, the space problem can be
with increasing in-situ stress. Xie et al. [7] reported that in-situ stresses simplified to a plane problem by appropriate simplification and me­
hinder tensile stresses but enhance compressive stresses in the blasthole chanical abstraction. Since the length of explosive is much larger than its
direction. Xie et al. [17] studied multi-hole blasting to cut rock under in- diameter, it can be simplified into a plane strain problem. Fig. 1 shows a
situ stress conditions and optimized the placement parameters for deep single blasthole in an infinite medium subjected to static stress load,
rock cutting and blasting. Yi et al. [18] showed that the early expansion assuming an elastic dynamic problem under plane strain conditions.
of blast cracks is determined by the dynamic blast load, and the later When the dynamic blast loads are not considered, the model can be
expansion of cracks is mainly influenced by in-situ stresses as the blast treated as an infinite plate containing a hollow hole. The plate is loaded
shock wave decays. Tao et al. [19] simulated the dynamic response of by biaxial stresses of magnitudes Px and Py in the horizontal and vertical
rocks under coupled static-dynamic loading and concluded that the in­ directions, respectively, and the radius of the hole is a.
crease in stress anisotropy concentrates cracks near the maximum The hoop stresses induced around the blasthole under biaxial stress
principal stress. Wang et al. [20] investigated the effect of in-situ stress loading were analyzed, and the complete solution for the stress distri­
on cyclic blasting damage extension and found that the contribution of bution was obtained [25]:
cyclic blasting load to rock fragmentation was mainly concentrated in ⎧ [ (
1 ( ) (a)2 ) ( )
( (a)2 (a)4 ) ]
the maximum principal stress direction, while the contribution to rock ⎪
⎪
⎪ σ = P + P 1 − + P − P 1 − 4 + 3 cos2θ
⎪
⎪ rr y x x y
fragmentation in the direction of the minimum principal stress was not ⎪
⎪
⎪
2 r r r
⎪ [ ( (a)2 ) ( ( (a)4 ) ]
significant. ⎨ 1 ( ) )
σ θθ = Py + Px 1 + − Px − Py 1 + 3 cos2θ
The study of the effect of in situ stress on the dynamic response of ⎪
⎪ 2 r r
⎪
rock blasting drew similar conclusions to the aforementioned study. ⎪
⎪
⎪ ( )
( (a)2 (a)4 )
⎪
⎪ 1
Specifically, blast crack expansion was enhanced in the direction of ⎪
⎩ τrθ = 2 Py − Px 1 + 2 r − 3 r sin2θ
maximum principal stress, and the crack expansion was suppressed in
other directions. However, most existing literature only reports the (1)
relationship between in-situ stress and blast cracking qualitatively, with
where σ rr is radial stress, σθθ is t hoop stress, τrθ is shear stress, r and θ are
few quantitative analyses of crack extension under the influence of in-
the distance and angle from the source of the blast, respectively.
situ stress. In recent years, some scholars have started to work on the
According to Eq. (1), when r = a, the magnitudes of σrr and τrθ are
quantitative description of blast cracking. Guo et al. [21] and Zhai et al.
zero, and the radial and shear stresses at the wall of the blasthole are
[22] quantified experimentally obtained blast cracks using fractal di­
free. The hoop stress distribution on the blasthole wall can represent the
mensions, while Yue et al. [23] and Zhang et al. [24] proposed the
stress distribution state. Fig. 2 illustrates the distribution of hoop stresses
concept of fracture degree based on CDEM simulation software. The
magnitude of fracture degree was used to quantitatively assess the
damage degree of concrete. Tao et al. [19] used the ImageJ software to
process the resultant images derived from ANSYS/LS-DYNA and reached
some quantitative useful conclusions. However, there are some disad­
vantages with this image-processing method. Specifically, ImageJ is
weak in processing color-rich images, and it cannot process images in
batches automatically. Rock fracture images obtained using the HJC
model simulations are available in few colors, which allow the ImageJ
software to be used to its full potential. For images of damage obtained
using the RHT model, the richness of the colors and the dispersion of the
damaged areas make it difficult to capture the damaged areas using
ImageJ software. Additionally, the existing studies have only processed
images of the final state of blasting, and there is a lack of quantitative
analysis of a large number of images throughout the whole rock blasting
process. This lack of analysis limits the in-depth understanding of the
evolution of the dynamic response of rocks under blast loads. Further­
more, most of the studies at this stage focus on the fracture phenomena
of rocks, and there is also a relative lack of research on the mathematical
relationship between rock damage and in-situ stress under the action of
dynamic blasting loads.
In this study, the stress distribution in rock mass under static loads
was determined using theoretical analysis, which forms the theoretical Fig. 1. Elasticity model under static load.

2
G. Qiao et al. Computers and Structures 287 (2023) 107116

Fig. 2. Hoop stress distribution on the blasthole wall: (a) in equiaxial in-situ stress fields, and (b) in anisotropic in-situ stress fields.

induced by 10 different in-situ stress conditions on the blasthole wall. the blasthole. The magnitude of the hoop stress in the x-axis direction
The biaxial stress loading conditions are presented in Table 1. shows an opposite trend. In the far field, the magnitude of the hoop
As shown in Fig. 2(a), the hoop stress on the blasthole wall exhibits stress in the Y-axis direction gradually tends to the value of Px , the
symmetric compressive stresses under hydrostatic pressure, and the magnitude of the hoop stress in the X-axis direction tends to gradually
magnitude of hoop stresses increases with the increase of hydrostatic reach the value of Py , and the magnitude of the hoop stress in the other
pressure. This indicates that the expansion of blast damage in all di­ directions is determined by the joint action of Px and Py . In addition, the
rections in the hydrostatic stress field will be hindered to the same magnitude of the hoop compressive stress in the far field in the Y-axis
extent, and the higher the in-situ stress, the stronger the hindrance. In direction is much larger than that in the X-axis direction. For this reason,
the anisotropic stress field, the loading stress in the vertical direction it can be expected that the blast damage will be mainly distributed in the
(Py ) is kept constant, and as the loading stress in the horizontal direction direction of maximum principal stress (the X-axis direction), and the
(Px ) gradually increases, the circumferential stress distribution at the development of damage in the direction of minimum principal stress
wall of the blasthole is no longer symmetric, as shown in Fig. 2(b). The (the Y-axis direction) will be more strongly suppressed.
hoop stress appears as compressive stress at θ = 90 and θ = 270 , in­
◦ ◦

creases sharply with increasing of Px . At θ = 0 and θ = 180 , the hoop 3. Numerical simulation of rock damage under in-situ stress
◦ ◦

compressive stress gradually decreases with the increase of the value of

Px . When Px reaches 80 MPa, the hoop stress starts to behave as tensile This study establishes a numerical model with a side length of 5 m,
stress. It is not reasonable to assume that the extension of blast damage and a blasthole with a diameter of 11 cm, and set up with non-reflecting
in the direction of the minimum principal stress will be hindered, while boundaries around it (Fig. 4). A plane strain model is used for numerical
blast damage extension will be promoted in the direction of the investigation in the study, which implies the assumption of infinite ve­
maximum principal stress when the hoop stress is expressed as tensile locity of detonation along the borehole and infinite charge length. Ten
stress. in-situ stress loading cases shown in Table 1 were numerically investi­
gated to replicate the rock blasting process. To take into consideration
the static loading effect of the in-situ stress, the rock is first pre-loaded to
2.2. Hoop stress distribution in rocks obtain an initial stress state. The rock in the initial stress state was then
introduced into the display analysis session as the initial condition.
Based on Eq. (1), the variation curve of the circumferential stress During the analysis phase, the explosive was detonated, and the rock is
with the distance to the blasthole wall was plotted (Fig. 3). Fig. 3(a) subjected to coupling dynamic-static load. Solid materials are generally
shows the stress concentration near the blasthole under equiaxial in-situ simulated using Lagrangian algorithms. However, the blasting process is
stress loading. As the distance from the blasthole increases, the magni­ considered a fluid–solid coupling problem, and therefore, fluid­
tude of the hoop stress gradually converges to the biaxially loaded stress –structure interaction methods are often applied. In this numerical
value for each case. Fig. 3(b) represents that the magnitude of the hoop simulation, the arbitrary Lagrangian-Eulerian (ALE) algorithm was
stress in the Y-axis direction gradually decreases with the distance from employed.

Table 1 3.1. Constitutive model and parameters of rocks

Biaxial stress loading conditions in this study.
Case Stress state Py /MPa Px /MPa K = Px /Py Several constitutive models have been developed to describe the
1 Equiaxial 0 0 1 response of rocks under dynamic loads, including the Holomquis-
2 20 20 1 Johnson-Cook (HJC) model [26], the Riedel-Hiermaier-Thoma (RHT)
3 40 40 1 model [27], and the JH series of models [28]. The RHT model is a
4 60 60 1
popular model for simulating the nonlinear response of rocks under
5 80 80 1
6 100 100 1 dynamic loading [29–31]. The pressure of the RHT model is described
7 Anisotropic 20 40 2 by the Mie-Gruneisen form, which includes a polynomial Hugoniot
8 20 60 3 curve and a p-α compaction relation [32]. The p-α compaction model is
9 20 80 4 shown schematically in Fig. 5(a). The model is in an elastic state when
10 20 100 5
the applied pressure is lower than the pore fragmentation pressure. As

3
G. Qiao et al. Computers and Structures 287 (2023) 107116

Fig. 3. Hoop stress variation with distance from the blasthole: (a) in equiaxial in-situ stress fields, and (b) in anisotropic in-situ stress fields.

FCAP(P) is the elastic yield surface cap function to limit the elastic bias
stress under elastic hydrostatic pressure.
Stage 2: Plastic hardening stage. The material is in the plastic
deformation phase when it is between the yield and failure surfaces. The
behavior of the material in this phase is described by strain-hardening
characteristics. The failure surface is defined as a function of the pres­
sure P, the Lode angle θ and the strain rate ε̇:
Yfail = YTXC(P) ⋅R3(θ) ⋅FRATE(ε̇) (3)
⃒ ⃒
⃒ N⃒
where, YTXC = fc ⃒⃒A(P* − P*spall FRATE ) ⃒⃒, fc is the uniaxial compressive

strength, A is the failure surface constant, N is the failure surface index,

P* is the normalized pressure according to fc , P*spall is defined as P* (ft /fc ),
and R3(θ) is the third invariant of the partial stress tensor.
Stage 3: Softening damage stage. Further, the material retains a
certain residual strength after failure, and the residual strength surface is
*
denoted as Yresid = B⋅P*M . B is the residual strength surface constant. M is
the residual strength surface index. The material behavior between the
failure surface and the residual surface is described by the damage
characteristics. Damage is accumulated by the equation:
Fig. 4. Configuration of the numerical model. ∑ ΔεPl
D= (4)
εfailure
P

the pores begin to collapse, the bulk stiffness and effective bulk modulus failure
of the material decrease. The variable α controls the porosity of the where, εP = D1 (P* − Pspall )D2 ⩾εmin
f , D1 and D2 are the damage con­
material, decreasing as the pressure increases, which results in an irre­ stants. εmin
f is the minimum failure strain. When D = 0, there is no
versible loading process. When the pressure reaches the compaction accumulated damage to the material. When D = 1, the residual surface
pressure, the material is completely compacted, at which point α equals is reached, and the material is completely damaged.
1. The strain rate effect is an important characteristic of dynamic load
The strength of the RHT model is described by three limit surfaces: that distinguishes it from the action of static load. The RHT model in­
the elastic yield surface, the failure surface, and the residual strength corporates the strain rate effect, and the dependence of strain rate on
surface, as shown in Fig. 5(b). The arrow in Fig. 5(b) depicts the typical strength can be expressed by the following equation:
loading process of the material, which is divided into three stages in ⎧ ( )β c
total. ⎪
⎪
⎪ c
⎪ ε̇ /ε̇ P⩾fc /3
Stage 1: Elastic deformation stage. The rock is in the elastic ⎪ p 0
⎪
⎪
⎪
⎪
deformation phase when the arrow does not reach the elastic yield ( ) ⎪ ⎨ ( )βc ( )β t
P + ft /3 P − fc /3
surface. the elastic yield surface Yelastic of the material is determined by Fr ε̇p =
⎪ fc /3 + ft /3
ε̇p /ε̇t0 − ε̇p /ε̇c0 − ft /3 < P < fc /3
⎪ fc /3 + ft /3
the failure surface Yfail , given by Eq. (2). ⎪
⎪
⎪ ( )β c
⎪
⎪
⎪
⎪
(2)
t
Yelastic = Yfail ⋅Felastic ⋅FCAP(P) ⎪
⎩ ε̇ p /ε̇ 0 P⩽ − ft /3

where Yelastic is the ratio of elastic strength to failure surface strength and (5)

4
G. Qiao et al. Computers and Structures 287 (2023) 107116

Fig. 5. RHT model: (a) Schematical description of the p-α equation of state, and (b) Stress limit surfaces and loading scenario in the RHT model.

( ) in-situ stress effects. The rock parameters used for in the simulations are
where Fr ε̇p is the strain rate strength factor, ε̇p is the strain rate, ε̇c0 is presented in Table 2, while the explosives parameters are shown in
Table 3. To maintain consistency with Banadaki’s experiments [36], the
the reference strain rate under compression, ε̇c0 = 3.0 × 10− 5 s− 1 , ε̇t0 is
model boundary was not set as a reflection-free boundary. In the RHT
the reference strain rate under tension, ε̇t0 = 3.0 × 10− 6 s− 1 , ft is the
model, the LS-DYNA commercial software uses rock damage distribution
tensile strength, βc is the material constant in compression, βt is the
to simulate blast-induced crack extension. Fig. 6(b) illustrates the blast
tensile material constant.
shock wave strength was greater than the compressive strength of the
In this paper, parameters of the granite calibrated by Li [33] are used
rock, resulting in a crush zone near the blasthole. Fig. 6(c) shows that a
and listed in Table 2.
cracked area is created by the tensile component of the stress wave. As
the stress wave propagated towards the model boundary, the compres­
3.2. Jones–Wilkins–Lee (JWL) EOS sive stress wave is reflected at the free surface, producing a tensile wave.
When the rock is subjected to tensile stress waves, it incurred tensile
In LS-DYNA, the Jones-Wilkens-Lee (JWL) equation of state is used to damage, with caused circumferential cracking and spalling (Fig. 6(d)).
describe the relationship between pressure, volume, and energy of a The final state of the rock after blasting was divided into three zones, as
high explosive product. The equation of state (EOS) is given in the shown in Fig. 6(e). This division corresponds to the pattern of the crack
following equation [34]: distribution observed in the granite experiments conducted by Banadaki
ω ω ωE (Fig. 6(f)) [36], indicating that the simulation parameters are reliable.
P = A(1 − )e− R1 V
+ B(1 − )+ (6)
R1 V R2 V V
3.4. Results of numerical simulation
where, P is the detonation pressure, A, B, R1 , R2 , ω are the explosive
characteristics parameters, E is the internal energy, V is the relative The damage distribution of the rock at different moments under the
volume. PETN explosives were used in this calculation. The explosive combined effect of blast load and in-situ stress is illustrated in Fig. 7 and
parameters were confirmed by Ayman [35], and the parameters are Fig. 8. Fig. 7 shows the damage pattern of the rock in equiaxial in-situ
reproduced in Table 3. stress fields, while Fig. 8 displays the damage pattern of the rock in
the anisotropic in-situ stress fields.
3.3. Verification of the performance of simulation parameters without in- In Fig. 7(a), the variation of rock damage with time is presented in
situ stress the absence of in-situ stress. During the initial stage of the explosion, the
blast shock wave level is much larger than the dynamic compressive
Fig. 6(a) shows the numerical model used for rock blasting without strength of the rock, resulting in severe radial compression damage in
the near zone of the detonation of the rock (Fig. 7(a)). As time passes, the
Table 2 shock wave energy gradually dissipates and decays into a stress wave.
RHT model parameters for granite. Although, the strength of the stress wave is not sufficient to cause radial
crushing damage to the rock, the dynamic tensile strength of the rock is
Parameter Value Parameter Value Parameter Value
much less than its compressive strength, the tensile stress induced by the
ρ0 (kg⋅m - 3 ) 2700 A1 (GPa) 86.71 βt 0.0144
stress wave leads to tensile damage to the rock [37]. Radial tensile
G(GPa) 24.17 A2 (GPa) 145.67 G*c 0.4
damages continue to increase in the rock as the simulation proceeds to
fc (GPa) 0.119 A3 (GPa) 89.03 G*t 0.7
N1 0.56 Q0 0.64 XI 0.48
800 μs, and their extent as well with the propagation of the stress wave
βc 0.0106 Bq 0.0105 D1 0.042 until the stress wave is not powerful enough to cause damage to the rock
B0 1.68 A 1.6 D2 1.0 Fig. 7(a). The final state of the rock is shown in the 2000 μs image, where
B1 1.68 Np 4.0 Pcrush (GPa) 0.04 the damage has stopped expanding Fig. 7(a).
1.1 Pcomp (GPa) 5.5 Af 1.62
α0 As described in Section 2.1, the hoop stresses on the blasthole wall in
T1 (GPa) 86.71 ε̇c (ms - 1 )
0
3.0E-8 Nf 0.6
equiaxial in-situ stress fields are compressive, which impedes blast crack
T2 (GPa) 0 ε̇t0 (ms - 1 ) 3.0E-9 εm 0.012
p
extension (see Fig. 2(a)). Fig. 7(b), Fig. 7(c), Fig. 7(d), Fig. 7(e), and
F*t 0.1 ε̇c (ms - 1 ) 3.0E22
Fig. 7(f), demonstrate a significant decrease in the extent of radial
F*s 0.38 ε̇t (ms - 1 ) 3.0E22
damage in the rock with increasing in-situ stress. An increase in

5
G. Qiao et al. Computers and Structures 287 (2023) 107116

Table 3
Parameters of JWL EOS for the PETN.
ρ0 (kg⋅m - 3 ) D(m⋅s - 1 ) PCJ (GPa) A(GPa) B(GPa) R1 R2 ω E0 (GPa)

1500 7450 22 625 23.3 5.25 1.6 0.28 8.56

Fig. 6. Numerical simulation validation: (a) simulation model; (b) crush zone generation; (c) crack zone generation; (d) spalling failure zone generation; (e) final
failure mode of the model; (f) final pattern of granite blast cracking obtained by test [36].

equiaxial in-situ stress leads to a rise in the magnitude of the hoop between Px and Py is equal, the circumferential stress distribution curve
compressive stresses in the blasthole wall and the rock, further inhibit­ on the wall of the gun hole is symmetrical.
ing blast cracking. A comparison of the 2000 μs moment damage images
in Fig. 7(a) and Fig. 7(f) reveals a reduction in the number of radial 4. Quantitative analysis of numerical simulation results
major tensile damage cracks from 16 to 8 and inhibition of radial
branching development. 4.1. Image processing techniques
Fig. 8 shows the distribution of rock blast damage caused by aniso­
tropic in-situ stress fields. As the value of K (K = Px /Py ) increases, the The post-processing software LS-PrePost of ANSYS/LS-DYNA does
deflection of the major radial tensile damage towards the direction of not provide access to macroscopic damage ranges. The amount of
the maximum principal stress becomes more pronounced. damage to the rock in LS-PrePost is presented on an elemental basis,
The static analysis conducted in Section 2.2, explains that, under which poses difficulties for the quantitative assessment of the macro­
anisotropic in-situ stress conditions (Fig. 2(b)), maximum compressive scopic damage characteristics of the rock model. Quantitative analysis is
stresses occur in the direction of θ = 90 and θ = 270 on the blasthole essential for a deep understanding of the effects of in-situ stress on rock
◦ ◦

walls, inhibiting the crack expansion in this direction. Conversely, crack blasting damage. In this paper, a method for quantitative assessment of
extension is more likely to occur in the direction of θ = 0 and θ = 180 . blast damage areas based on image processing techniques is proposed.
◦ ◦

These theoretical findings explain the results of the numerical simula­ This method allows for the precise extraction of the target area (damage
tions in Fig. 8 very well. Furthermore, all the crush zones in Fig. 8 are area) from the numerically simulated image, which is then used as the
elliptical, with the long axis of the ellipse aligned with the direction of basis for rapid computational analysis of the target area. Previous
the maximum principal stress. As the value of K increases, the long axis studies used numerical simulations with CDEM software by Zhang et al.
of the ellipse grows, and the short axis shortens. These results are [24]. The fracture degree, defined as the ratio of the number of fractured
consistent with those reported in the literature [8,10]. contact elements to the total number of contact elements, was intro­
Note that by rotating the rock damage image in Fig. 8 by 90◦ in the duced to quantitatively describe the damage of rock materials [24]. A
instantaneous clockwise direction, a map of the rock damage under the new concept called the damage variable, which is similar to the fracture
influence of different magnitudes of vertical ground stress can be ob­ degree, to describe the characteristics of the ANSYS/LS-DYNA numerical
tained. This is because when the absolute value of the difference simulation software is proposed in this paper. The damage variable is

6
G. Qiao et al. Computers and Structures 287 (2023) 107116

Fig. 7. Rock damage distribution in equiaxial in-situ stress fields.

denoted by η, η is given by Eq. (7): simulation, S is the total area of the numerical model (total number of
pixel points in the picture). The magnitude of the η value reflects the
η = St × 100\% /S (7)
extent of the damage. The image processing technique consists of four
where, St is the damaged area (sum of damage pixels) in the numerical steps, as shown in Fig. 9.

7
G. Qiao et al. Computers and Structures 287 (2023) 107116

Fig. 8. Rock damage distribution in anisotropic in-situ stress fields.

Fig. 9. Steps in image processing technology.

Step1: Substract the white edges of the images from the numerical ⎡ ⎤
f (1, 1) ⋯ f (1, N)
simulations. The purpose of this step is to avoid the white blank areas
F(M, N) = ⎣ ⋮ ⋱ ⋮ ⎦ (8)
affecting the accuracy of the calculation results.
f (M, 1) ⋯ f (M, N)
Step 2: The original Red-Green-Blue (RGB) digital images are con­
verted into grey-scale images. The purpose of this step is to prepare for a where, F(M, N) is an M × N matrix of gray values (pixel). f(i, j) is the gray
binary image. value of a pixel, and f(i, j) = (R × 30 + G × 59 + B × 11 + 50)/100, R,
The digital image that has been converted to a greyscale format G, and B are RGB values.
contains an M × N matrix of grey values (pixels) of luminance infor­ Step 3: Based on the obtained grey-scale histograms, a threshold is
mation, expressed by Eq. (8) [38]. selected using the Otsu method (global thresholding algorithm) to

8
G. Qiao et al. Computers and Structures 287 (2023) 107116

binarise the grey-scale images.

In some simple images where the target regions and the background
are distinctly different, the selection of the threshold value is very
simple. The grey scale threshold of an image can be selected at the
bottom of the valley between the two peaks representing the target re­
gions and the background in the grey scale histogram [39]. However, for
the rock damage images covered in this paper, different levels of damage
are represented by different colors in the image, and each image con­
tains a large and complex range of colors. The grey scale thresholds
selected directly from the grey scale histogram are inaccurate, which
will cause significant errors in the calculation results. The Otsu method
offers a solution to the problem. The Otsu method is an automatic
threshold selection method that relies on the maximum between-class
variance between the background and target regions [40].
The total number of pixels in a greyscale image is given by equation
Eq. (9), where gmax is the maximum grey scale value in the image:
gmax
∑
N= ni (9)
i=1

n(g) is the total number of pixels in the image with a grey value of g. Pg
is the proportion of pixels in the image with a grey value g. The rela­ Fig. 10. Evolution of damage variables in equiaxial in-situ stress fields.
tionship between Pg and n(g) is expressed by Eq. (10).
n(g) evolution of the rock damage variable. That is, the higher the in-situ
Pg = (10)
N stress, the smaller the damage variable in all three phases.
In addition, the probability distribution of Pg is given by Eq. (11): When the blasting process is completed at the 2000th μs, the damage
distribution can be considered as the final damage state of the model. To
establish a mathematical relationship between the damage variables and
gmax
∑
Pg = 1 (11)
i=1
the magnitude of the in-situ stress in the final state of the model, the
damage variables for different equal biaxial loading are summarized in
According to the principles of the Otsu method for selecting Fig. 11. The damage variable is inversely proportional to the hydrostatic
thresholds, the threshold T0 is given by Eq. (12). pressure, as reflected in Fig. 11, and a good linear relationship was found
{
T0 = argmax σ2B (g)
}
(12) after fitting. In this paper, the damage variable related to the equiaxial
1⩽g⩽gmax pressure Px (Py ) can be expressed by Eq. (13):

where σ 2B (g) is the maximum between-class variance. η = 20.243 − 0.931Px (13)

Step 4: The pixels of the binarised image only have grey values of
0 and 255. After the binarization process, the image is traversed, and the 4.2.2. Variation of damage variables in anisotropic in-situ stress fields
pixels with a grey value of 255 are added up to obtain the total number Fig. 12 illustrates the time course of rock damage evolution in
of pixels in the target region. The total number of pixels is defined as the anisotropic stress fields. The time-course curve trends for the rock
area of the target region St . Finally, by inputting St into the Eq. (7), the damage variables shown in Fig. 12 are consistent with those in Fig. 11
final value of η can be obtained. and can be divided into three phases. The phases correspond to the
compressional damage phase caused by the dominant blast shock wave
(Zone 1), the rapid expansion of radial cracks caused by the dominant
4.2. Quantitative analysis of rock damage

One image of damage was captured every 100 μs, starting from 0 μs
mark. A total of 21 images were captured for each calculation. The
damage variables in each of the numerical simulation plots are calcu­
lated using image processing techniques, and the variation curves of
damage variables are plotted with increasing time.

4.2.1. Variation of damage variables in equiaxial in-situ stress fields

Fig. 10 illustrates the dynamics of the damage varies over time,
which shows the three main stages in the evolution of damage variables
under equiaxial loading. The first stage (Zone 1) ranges from 0 μs to 400
μs it is dominated by the damage caused by the blast shock wave, as the
rock is compressed during this phase and the damage variables increase
rapidly. The second stage, between 400–μs––1400 μs (Zone 2), is
dominated by blast stress waves, as radial cracking continues to expand
under the tangential component of the blast stress wave. The growth rate
of the damage variable is slower than in the first stage. The third stage
ranges from 1400 μs to 2000 μs (Zone 3), at which point the blast stress
wave can no longer cause extensive damage to the rock. With the end of
the blast, the damage zone in the rock no longer increases, and the
growth of the damage variable tends to stagnate. In addition, there is an
observable suppressive effect of in-situ stress on all three stages of the Fig. 11. Relationship between damage variable and equiaxial static stress.

9
G. Qiao et al. Computers and Structures 287 (2023) 107116

the damage variable η can be described by the following equation:

η = 18.336 + 0.16271Px - 0.09929P2x (14)

4.2.3. Effect of anisotropy ratio on damage variables

The effect of anisotropy ratios on rock damage distribution is sig­
nificant. Blasting damage is approximately symmetrically distributed in
equiaxial in-situ stress fields (Fig. 7). In contrast, blast damage in
anisotropic stress fields is mainly distributed in the direction of the
maximum principal stress, with less damage distribution in the direction
of the minimum principal stress (Fig. 8). This conclusion has been
confirmed in previous studies [17–20,41,42]. However, the quantitative
description of the difference between blast damage in the direction of
maximum principal stress and minimum principal stress is rarely re­
ported. In addition, the quantitative assessment of the change in rock
damage with increasing in-situ stress anisotropy has also been reported
only in a few previous studies. Yue et al. [23] set up 10 regions in the
model to describe the degree of fracture in the loading direction and in
the vertical direction of loading. However, the 10 zones set did not cover
the entirety of the model, and no conclusions were drawn on the global
nature of in-situ stress damage to the rock. Thanks to the efficiency of
Fig. 12. Evolution of damage variables in anisotropic in-situ stress fields. image processing techniques, a global assessment of rock damage was
achieved in this paper. As shown in Fig. 14, the numerical model is cut
blast stress wave (Zone 2), and the final stage of the blast (Zone 3). An into four triangular regions. To quantify the effect of changes in the
interesting phenomenon emerges from the zoomed-in view in Fig. 12. value of the anisotropy ratio K(K = Px /Py ) on the damage distribution,
The magnitude of the rock damage variables η in Cases 9 and 10 exceeds the anisotropic damage variables ηk is presented, ηk defined by Eq. (15):
those in Cases 7 and 8, ranging from 400 μs to 800 μs. After 800 μs, the ηk = Sx × 100\% /Sy (15)
rock damage variables η in Cases 9 and 10 become progressively smaller
than in Cases 7 and 8. where Sx is the sum of the areas of damage appearing in zones II and IV
This phenomenon occurs because the maximum hoop tensile stresses and Sy is the sum of the areas of damage appearing in zones I and III. The
in the blasthole wall of Cases 9 and 10 occur in the horizontal direction four zones in Fig. 14 cover the entire model, so the results of the study
(see Fig. 2(b)), which contributes to blast crack extension and results in a can reflect how the anisotropy ratio affects the model as a whole.
rapid increase in blast damage variables. As one moves further away Step 4 of Section 4.1 was improved to apply to the new image pro­
from the hole, the annular stress gradually changes from tensile to cessing. This is done by giving traversal boundaries of each image ac­
compressive (see Fig. 3(b)), which changes the effect on blast damage cording to the range shown in Fig. 14 and thus obtaining the damaged
from a facilitator to a hindrance, resulting in a progressively smaller area of Sx and Sy . Finally, the value of ηk for the anisotropic damage
damage. variable is calculated according to Eq. (15). The trends over time for
Different degrees of anisotropy were stablished by fixing Py equal to each of the four anisotropic damage variables, obtained using image
20 MPa and making Px equal to 20 MPa, 40 MPa, 60 MPa, 80 MPa, and processing techniques, are shown in Fig. 15.
100 MPa, respectively. The damage variables for the final damage state It is observed in Fig. 15 that the value of ηk rises with K at the same
of the rock in different anisotropic in-situ stress fields were calculated moment in time. In addition, the curves in Case 7 (K = 2) and Case 8
and plotted (Fig. 13). Meanwhile, Fig. 13 suggests that as the horizontal (K = 3) are flat and have smaller values. The curves for Case 9 (K = 4)
stress Px (or anisotropy ratio) is raised, the damage variable is reduced. and Case 10 (K = 5) increase sharply at around 400 μs and remain at
By fitting the data points, it was found that when Py = 20MPa, Px and large values thereafter. The results of the theoretical analysis in Section
1. can explain the dramatic increase in the ηk -curves of Case 9 (K = 4)
and Case 10 (K = 5). When K equals 4 or 5, the hoop stress in the

Fig. 13. Relationship between damage variable and anisotropic static stress. Fig. 14. Division in the numerical model.

10
G. Qiao et al. Computers and Structures 287 (2023) 107116

MPa, the relationship between the damage variable η and Px in

the final state can be represented by
η = 18.336 + 0.16271Px - 0.09929P2x .
(4) To overcome the challenge of quantitatively describing the effect
of the anisotropy ratio K on damage distribution, the anisotropic
damage variable ηk was introduced. The analysis demonstrates
that ηk is positively correlated with K at the same moment. When
the hoop stress in the direction of the maximum principal stress is
tensile, rock damage is promoted, leading to a sharp increase in
ηk .

CRediT authorship contribution statement

Guodong Qiao: Conceptualization, Investigation, Software, Writing

– original draft, Visualization. Zegong Liu: Resources, Supervision,
Funding acquisition. Changping Yi: Writing – review & editing, Meth­
odology. Kui Gao: Software, Validation. Gaoyuan Xuan: Formal
analysis.

Declaration of Competing Interest

Fig. 15. Time course curves for each anisotropic damage variable.

The authors declare that they have no known competing financial

horizontal direction of the blasthole wall is tensile, and the initial blast
interests or personal relationships that could have appeared to influence
damage is created first in the direction of the maximum principal stress. the work reported in this paper.
The existing cracks have a direct effect on the further extension of blast
cracks, which provides good conditions for the rapid development of Data availability
blast damage in that direction. When K equals 2 or 3, there is no such
situation. It is speculated that this is consistent with the cause of the No data was used for the research described in the article.
anomalous change in the damage variable η at the 400 μs time in the
enlarged plot of Fig. 12 in Section 4.2.2.
Acknowledgments
5. Conclusions
This work was supported by the Graduate Research Project of Anhui
Education Department (YJS20210396), National Natural Science
This paper investigates the effects of in-situ stress on rock blast
Foundation of China (52074013), and Natural Science Foundation of
damage using theoretical analysis and numerical simulation methods.
Anhui Province (2208085ME125), which are gratefully appreciated.
Although numerical simulations have been widely used in the study of
Thanks to Carlota Rodriguez San Miguel for proofreading the article.
the dynamic response of deep rock, previous studies have lacked a
quantitative description of model damage phenomena. Therefore, in this
paper, a processing of images from ANSYS/LS-DYNA numerical simu­ References
lation results were processed using image processing techniques to
[1] Zhao X. Basic theory and development Trends of ultra-deep shaft construction.
quantitatively assess the distribution of rock blast damage. Metal Mine 2018;04:1–10. https://doi.org/10.19614/j.cnki.jsks.201804001.
[2] Lucier AM, Zoback MD, Heesakkers V, Ze R, Murphy SK. Constraining the far-field
(1) The analysis of rock stresses under static loading was carried out, in situ stress state near a deep South African gold mine. Int J Rock Mech Min Sci
2009;46(3):555–67. https://doi.org/10.1016/j.ijrmms.2008.09.005.
and the distribution of stresses on the borehole wall as well as in [3] Schweitzer JK, Johnson RA. Geotechnical classification of deep and ultra-deep
the far field was determined. The analysis demonstrates that Witwatersrand mining areas. South Africa Mineralium Deposita 1997;32(4):
under isotropic loading, the hoop compressive stresses are sym­ 335–48. https://doi.org/10.1007/s001260050100.
[4] Yuan L. Strategic thinking of simultaneous exploitation of coal and gas in deep
metrically distributed on the borehole wall, and the magnitude of mining. J China Coal Soc 2016;41(1):1–6. https://doi.org/10.13225/j.cnki.
the far-field stresses converges to the magnitude of the hydro­ jccs.2015.9027.
static pressure. Under anisotropic loading, the maximum [5] Malan DF. Manuel Rocha medal recipient simulating the time-dependent
behaviour of excavations in hard rock. Rock Mech Rock Eng 2002;35(4):225–54.
compressive stress on the blasthole wall occurs in the direction of https://doi.org/10.1007/s00603-002-0026-0.
the minimum principal stress, with the far-field stress tending [6] Kulatilake PHSW, Wu Q, Yu Z, Jiang F. Investigation of stability of a tunnel in a
towards Py in the x-axis direction and Px in the y-axis direction. deep coal mine in China. Int J Min Sci Technol 2013;23(4):579–89. https://doi.
org/10.1016/j.ijmst.2013.07.018.
(2) The damage variable was introduced, and the time-course curve
[7] Xie L, Lu W, Zhang Q, Jiang Q, Wang G, Zhao J. Damage evolution mechanisms of
of the damage variable was plotted using image processing rock in deep tunnels induced by cut blasting. Tunn Undergr Space Technol 2016;
techniques to analyze the damage cloud at different moments. 58:257–70. https://doi.org/10.1016/j.tust.2016.06.004.
[8] Nicholls, HR, Duvall WI. Presplitting rock in the presence of a static stress field. US
The dynamic evolution of rock blasting damage under in-situ
Department of the Interior, Bureau of Mines. Washington; 1966.
stress was categorized into three stages on the curves trend. In [9] Yang L, Yang R, Xu P, Song Y. Experimental study on the effect of initial
equiaxial in-situ stress fields, the evolution of blast damage is compression stress field on blast-induced crack behaviors. J China Coal Soc 2013;
suppressed throughout due to the effect of in-situ stress. In 38(3):404–10. https://doi.org/10.13225/j.cnki.jccs.2013.03.001.
[10] Zhang F, Peng J, Qiu Z, Chen Q, Li Y, Liu J. Rock-like brittle material fragmentation
anisotropic stress fields, the tensile stress in the direction of the under coupled static stress and spherical charge explosion. Eng Geol 2017;220:
maximum principal stress causes the damage variable deviate 266–73. https://doi.org/10.1016/j.enggeo.2017.02.016.
from the being smaller at a given time with higher in-situ stresses. [11] He C, Yang J, Yu Q. Laboratory study on the dynamic response of rock under blast
loading with active confining pressure. Int J Rock Mech Min Sci 2018;102:101–8.
(3) In the final state, the damage variable η and the hydrostatic https://doi.org/10.1016/j.ijrmms.2018.01.011.
pressure Px (Py ) can be fitted using the equation: η = [12] Yue Z, Tian S, Zhang S, Zong L, Xu S. Expanding law of cracks formed by slotted
20.243 − 0.931Px . In an anisotropic stress field, when Py is 20 cartridge blast under unidirectional confining pressure. J Vib Shock 2019;38(23):
186–95. https://doi.org/10.13465/j.cnki.jvs.2019.23.027.

11
G. Qiao et al. Computers and Structures 287 (2023) 107116

[13] Wang M. Experimental study on explosion stress field and damage mechanism of [27] Riedel W, Thoma K, Hiermaier S, Schmolinske E. Penetration of reinforced
rock under confining pressure. Anhui University of Science and Technology; 2019. concrete by BETA-B-500 numerical analysis using a new macroscopic concrete
[14] John O. LS—DYNA keyword user’s manual. version 970. Livermore Software model for hydrocodes. In: Proceedings of the 9th International Symposium on the
Technology Corporation. California; 2003. Effects of Munitions with Structures; 1999.
[15] Ma G, An X. Numerical simulation of blasting-induced rock fractures. Int J Rock [28] Johnson GR, Holmquist TJ. An improved computational constitutive model for
Mech Min Sci 2008;45(6):966–75. https://doi.org/10.1016/j.ijrmms.2007.12.002. brittle materials. In: AIP conference proceedings; 1994.
[16] Lu W, Leng Z, Chen M, Yan P, Hu Y. A modified model to calculate the size of the [29] Jayasinghe LB, Shang J, Zhao Z, Goh ATC. Numerical investigation into the
crushed zone around a blast-hole. J Southern African Inst Mining Metallurgy 2016; blasting-induced damage characteristics of rocks considering the role of in-situ
116(5):412–22. https://doi.org/10.17159/2411-9717/2016/v116n5a7. stresses and discontinuity persistence. Comput Geotech 2019;116:103207. https://
[17] Xie L, Lu W, Zhang Q, Jiang Q, Chen M, Zhao J. Analysis of damage mechanisms doi.org/10.1016/j.compgeo.2019.103207.
and optimization of cut blasting design under high in-situ stresses. Tunn Undergr [30] Liu K, Li Q, Wu C, Li X, Li J. A study of cut blasting for one-step raise excavation
Space Technol 2017;66:19–33. https://doi.org/10.1016/j.tust.2017.03.009. based on numerical simulation and field blast tests. Int J Rock Mech Min Sci 2018;
[18] Yi C, Johansson D, Greberg J. Effects of in-situ stresses on the fracturing of rock by 109:91–104. https://doi.org/10.1016/j.ijrmms.2018.06.019.
blasting. Comput Geotech 2018;104:321–30. https://doi.org/10.1016/j. [31] Yi C, Sjöberg J, Johansson D. Numerical modelling for blast-induced fragmentation
compgeo.2017.12.004. in sublevel caving mines. Tunn Undergr Space Technol 2017;68:167–73. https://
[19] Tao J, Yang X-G, Li H-T, Zhou J-W, Fan G, Lu G-D. Effects of in-situ stresses on doi.org/10.1016/j.tust.2017.05.030.
dynamic rock responses under blast loading. Mech Mater 2020;145:103374. [32] Borrvall T, Riedel W. The RHT concrete model in LS-DYNA. In: Proceedings of the
https://doi.org/10.1016/j.mechmat.2020.103374. 8th European LS-DYNA Users Conference; 2011.
[20] Wang H, Wang Z, Wang J, Wang S, Wang H, Yin Y, et al. Effect of confining [33] Li H. The study of the rock RHT model and to determine the values of main
pressure on damage accumulation of rock under repeated blast loading. Int J parameters. China Univ Mining Technol 2016.
Impact Eng 2021;156:103961. https://doi.org/10.1016/j.ijimpeng.2021.103961. [34] Lee E, Hornig H, Kury J. Adiabatic expansion of high explosive detonation
[21] Guo T, Zhang S, Ge H, Wang X, Lei X, Xiao B. A new method for evaluation of products. TID4500-UCRL 50422. USA: Lawrence Livermore National Laboratory,
fracture network formation capacity of rock. Fuel 2015;140:778–87. https://doi. University of California. Livermore; 1968.
org/10.1016/j.fuel.2014.10.017. [35] Ayman T. Hard rocks under high strain-rate loading. Queen’s university; 2010.
[22] Zhai C, Xu J, Liu S, Qin L. Fracturing mechanism of coal-like rock specimens under [36] Banadaki MMD. Stress-Wave Induced Fracture in Rock Due to Explosive Action.
the effect of non-explosive expansion. Int J Rock Mech Min Sci 2018;103:145–54. University of Toronto; 2010.
https://doi.org/10.1016/j.ijrmms.2018.01.037. [37] Hustrulid W. Blasting principles for open pit mining. Rotterdam: A.A. Balkema;
[23] Yue Z, Zhou J, Feng C, Wang X, Peng L, Cong J. Coupling of material point and 1999.
continuum discontinuum element methods for simulating blast-induced fractures [38] Zhang P, Zhao Q, Tannant DD, Ji T, Zhu H. 3D mapping of discontinuity traces
in rock. Comput Geotech 2022;144:104629. https://doi.org/10.1016/j. using fusion of point cloud and image data. Bull Eng Geol Environ 2019;78(4):
compgeo.2021.104629. 2789–801. https://doi.org/10.1007/s10064-018-1280-z.
[24] Zhang Q, Zhi Z, Feng C, Cai Y, Pang G, Yue J. Investigation of concrete pavement [39] Prewitt JM, Mendelsohn ML. The analysis of cell images. Ann N Y Acad Sci 1966;
cracking under multi-head impact loading via the continuum-discontinuum 128(3):1035–53. https://doi.org/10.1111/j.1749-6632.1965.tb11715.x.
element method. Int J Impact Eng 2020;135:103410. https://doi.org/10.1016/j. [40] Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst
ijimpeng.2019.103410. Man Cybern 1979;9(1):62–6. https://doi.org/10.1109/TSMC.1979.4310076.
[25] Brady BH, Brown ET. Rock mechanics: for underground mining. London: Kluwer [41] Donzé FV, Bouchez J, Magnier SA. Modeling fractures in rock blasting. Int J Rock
Academic Publishers; 1985. Mech Min Sci 1997;34(8):1153–63. https://doi.org/10.1016/S1365-1609(97)
[26] Holmquist TJ, Johnson GR. A computational constitutive model for glass subjected 80068-8.
to large strains, high strain rates and high pressures. J Appl Mech 2011;78(5). [42] Panchadhara R, Gordon PA, Parks ML. Modeling propellant-based stimulation of a
https://doi.org/10.1115/1.4004326. borehole with peridynamics. Int J Rock Mech Min Sci 2017;93:330–43. https://
doi.org/10.1016/j.ijrmms.2017.02.006.

12