# Explosion-Induced Stress Wave Propagation in Interacting Fault System: Numerical Modeling and Implications for Chaoyang Coal Mine

Authors: Xiaojun Feng; Qiming Zhang; Muhammad Ali
Journal: Shock and Vibration
 (2019)

---

## Abstract

Exploring the propagation of stress waves in rocks with preexisting discontinuities is of great importance to reveal rock and geological engineering problems, particularly dynamic disasters like earthquakes and rockbursts in underground coal mining. In this paper, six 3D models established with COMSOL Multiphysics are employed to explore the influence of two preexisting faults with different orientations on the propagation process of explosion-induced stress waves and the reflection effect. Considering the propagation process of stress waves, the interactive effect between two different size faults is discussed. The results show that the dip angles of the preexisting fault and the differences of the elastic modulus, density, and Poisson’s ratio between faults and rocks have great influence on the distribution of stresses and strain-energy density. Immediately after the stress wave induced by blasting arrived at preexisting fault A, a relatively high concentration of the strain-energy density was observed at the last wave before passing through fault A. The presence of faults leads to the reflection of most of the blast energy. When the stress wave propagates across fault A, the strain energy stored in the stress wave becomes attenuated; thus, most strain energy was absorbed by the fault’s domain. Finally, the modeling results were implicated in Chaoyang Coal Mine to account for the distribution of the observed seismic events. This study has guiding significance for the attenuation law of stress waves passing through joint/fissure zones in geological engineering, earthquake engineering, and underground mining engineering.


---

## Body

## 1. Introduction

Understanding stress wave propagation is of great significance to the constitutive relationship of rock mass, slope stability, earthquakes, and underground coal mining. There is a large variety of artificial or natural cracks, gaps, various joints, faults, folds, and interlayers inside rocks, along which shear is more likely to occur and deformation is usually nonlinear [1]. These structural planes have nonuniformities and discontinuities, and the cores of these structural planes are usually wet or filled with water [2], which can significantly reduce the strength of the rock mass and have an important influence on the propagation of stress waves in the rock mass [3–8]. When elastic waves propagate in these rock structures, reflections, refractions, and diffractions occur.The shock waves generated from the explosion are abruptly attenuated with the increase of propagation distance. In the process of rock fragmentation, the stress wave plays an extremely important role. Generally, the natural rock mass is not a homogenous body. There are many weak structural planes in rock mass, such as faults, joints, and fissures, which seriously hinder the propagation of stress waves and exacerbate the attenuation of stress wave energy.In recent years, many scholars have studied the dynamic response of structural planes under the action of explosive stress waves. Kang et al. [9] carried out ground stress measurement work on fault areas and found that the regional stress field is closely related to the faults’ distribution characteristics [10–16]. Sebastian and Sitharam [17] used FLAC 3D simulation to study the attenuation law of stress wave propagation at the frictional structural plane and discussed the influence of joint parameters on the energy attenuation of stress waves. Xie et al. [18–20] studied the propagation characteristics of explosive stress waves in the interstitial space of a medium. The experimental results show that the porosity of the medium significantly changes the propagation of the stress wave. As the porosity increases, the stress wave passes through the pore. The larger the amplitude is of the reflected wave being generated, the smaller the amplitude is of the transmitted wave. Needham [21] conducted a detailed study of the formation and propagation characteristics of explosive stress waves and pointed out the difference between the stress wave propagation of the explosion stress wave in the incomplete detonation state and the complete detonation. Moreover, he also emphasized the propagation characteristics of stress waves in explosive stress waves when they encounter different structures. Xie et al. [22] analyzed the influence of the angle between the dip direction of the joint and the stress wave on the propagation law of explosive stress waves. The results showed that the smaller the angle and average transmittance of the joint were, the greater the attenuation rate of the stress wave amplitude and the faster the energy dissipation were. Deng et al. [23] developed a coupling algorithm for UDEC and AUTODYN-2D to study the attenuation of stress waves through a single joint, two parallel joints, and multiple joint faces. However, the propagation mechanism of explosive stress waves is complex, and faults in different orientations have an important influence on the propagation of explosive stress waves. Therefore, it is necessary to further study the influence of the structural plane/discontinuity in rock bodies (e.g., granite) on the propagation process of stress waves.In this paper, a 3D model is established to explore the influence of two preexisting faults with different orientations on the propagation process of explosion-induced stress waves and the reflection effect. The interactive effect between two different size faults is also discussed by considering the transmissive process of the stress waves. Finally, the modeling result is implicated in Chaoyang Coal Mine to account for the spatial distribution of observed seismic events. This study may shed new light on the attenuation law of stress waves passing through joint/fissure zones in geological engineering, earthquakes, and underground mining engineering.


## 2. Method and Model Setup

The vibration in coal and rock mass often has the characteristics of short duration and strong destructive power. In a short period, the vibration causes much damage to the coal and rock mass near the source, forming a fracture zone. Because the propagation of the shock wave in the coal and rock mass is attenuated, the coal and rock mass far from the source forms a plastic zone and an elastic zone in turn. Considering that the source has little impact relative to the entire coal formation, when the gas stress and gas coupling laws are included under the simulated vibration load, the coal rock mass is regarded as linear elastic material, and the vibration wave in the coal rock body is regarded as an elastic wave.Based on this, the following basic assumptions are proposed:(1)
The coal and rock mass are considered an elastic medium, and the relationship between the stress and strain of the coal rock body obeys Hooke’s law(2)
The coal and rock mass are treated as a homogeneous medium, and thus each coal rock layer is composed of the same material(3)
The coal and rock mass are an isotropic medium; that is, the physical properties of the object are the same in all directions(4)
The deformation of the coal and rock mass under the action of force is small, far less than the original geometric size of the objectUnlike the stress balance equation of coal and rock under the action of static load, the coal rock deforms with the vibration duration when the vibration occurs. Therefore, the differential equation of motion describes the deformation of the microbody under the action of vibration load:(1)ρcu¨i=σij,j+fi,where ρc is the density, u¨i is the two differentials of displacement versus time, σij,j is the coal and rock mass stress, and fi is the coal and rock mass volume force.This 3D model is established with the mechanics module of COMSOL Multiphysics. The model geometry consists of two preexisting faults with different orientations (detailed dip angles of faults can be found in Table1). The dip angles of fault A were varied from 0–135°. The dip angles of fault B were set to 75° for five cases and 105° for one case. The modeling domain is 400 m long in the X direction, 300 m wide in the Y direction, and 150 m high in the Z direction.Table 1
Designed experiments with employing different dip angles for the two faults.
Case (°)
D0
D60
D75
D90
D135
D135R
α
0
60
75
90
135
135
β
75
75
75
75
75
105With regard to the blasting load simulation, a short duration load was applied on the left surface of the model (Figure1), which is a typical load/response during tunnel construction and other excavations that uses blasting in underground construction work and/or mining. An example of pressure function is shown in Figure 1 and originated from equation (2) [24]. The boundary conditions use the low-reflecting boundary conditions to truncate the computational domain to a reasonable size:(2)Pt=13.3⋅P0exp−γt2−exp−2γt,where P0 is the load magnitude of the blasting load and γ is the decay rate.Figure 1
Schematic geometry of the designed model. Two faults (A and B) with different dip angles (α and β) are installed in the model area. A blast loading function (purple curve) is applied at a rectangle area which is located on the left wall of the model (indicated by the red lines), with a range of 6 × 4 = 24 m2 (indicated by the brown area).The low-reflecting boundary condition takes the material data from the adjacent domain in an attempt to create a perfect impedance match for both pressure waves and shear waves:(3)σ⋅n=ρCp∂u∂t⋅nn−ρCs∂u∂t⋅tt,where n and t are the units normal and tangential vectors at the boundary, respectively, and Cp and Cs are the speeds of the pressure and shear waves in the material:(4)Cp=λ+2μρ,Cs=μρ,where λ is the first Lamé parameter, μ is the shear modulus, and ρ is the density of the material through which the wave propagates.Regarding the other parameters (Table2), they are calculated according to the previous founds [25, 26], and the load magnitude is calculated according to the amount of explosive with 1.4 × 108 × Q2/3 = 1.03 × 109 N; the loading duration is obtained with 0.81 × 10−3 × Q1/3 = 2.2 × 10−3 s; and the decay rate is calculated with loading duration (t0) by −(20.5 × log (0.5))/t0 = 445.841/s. All these values were input into COMSOL for the calculation of the stress wave propagation induced by the blast loading.Table 2
Detailed parameters used in numerical modeling.
Parameters
Symbol
Value-units
Rock (granite)
Density
ρ
r
2700 kg/m3
Young’s modulus
E
r
2.36 × 1010 Pa
Poisson’s ratio
Nu
r
0.25
Fault domain
Density
ρ
f
1350 kg/m3
Young’s modulus
E
f
1.97 × 109 Pa
Poisson’s ratio
Nu
f
0.35
Fault A dip angle
α
—
Fault B dip angle
β
—
Gravitational acceleration
G
9.81
Decay rate
γ
445.84 1 s
Loading duration
t
0
2.2 × 10−3s
Amount of explosive
Q
20 kg
Load magnitude
P
0
1.03 × 109 N


## 3. Result

### 3.1. Distribution of Stress Influenced by the Dip Angle

#### 3.1.1. Horizontal and Vertical Fault Planes

Horizontal delamination and vertical fault can be used as extreme geological conditions for exploring the role of preexisting discontinuity on the propagation process of the stress wave in civil engineering, earthquake engineering, and underground mining engineering. This paper first examines the difference between the two cases D0 and D90.For case D0 (Figure2), the results show that the blast loading spreads out in a wave and it propagated forward to fault A at about 0.15–0.02 s. At about 0.02 s, the front part of fault A was surrounded by the stress wave, with a large range of tensile stress and a relatively small range of compressional stress. The difference in the Sxx stress between fault A and the surrounding granite rock body was obvious. At about 0.03 s and 0.04 s, the maximum tensile and compressional stresses was decreasing, and a dislocated feature of the stress wave was observed in a triangular area along fault A, which is different from Case D90 (indicated with light red dotted lines in the first column in Figure 2). This is probably due to the effects of natural vibration and/or absorption during the propagation process of the stress wave. At about 0.09 s, the distribution of the Sxx stress at the left part of fault A dislocated, while at the right part of fault A, the Sxx stress showed highly zoned areas of tensile and compressional stresses. In particular, the small fault B presented a truncated effect on the propagation of the stress wave. At about 0.15 s, most of the modeling domain presented a dislocated distribution of the Sxx stress, only some highly concentrated stress was observed between the two faults and in the middle part of fault A, which was probably due to the reflection effect on the stress wave as a result of the presence of fault B.Figure 2
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 0° and 90°, respectively. The fault positions were indicated with black lines and purple front (A and B).For Case D90, the stress wave was propagated in the vicinity of fault A at about 0.03 s. At about 0.04 s, fault A started to effectively act as a significant barrier in the propagation of the stress wave, presenting a highly concentrated tensile stress (as indicated as red dotted ellipse area in Figure2). At about 0.09 s, the distribution of the Sxx stress near fault A showed an obvious reflection effect according to the symmetrical semicircular characteristic of the tensile and compressional stresses. Compared with Case D0, this case was less dislocated for the distribution of the Sxx stress in the modeling domain. At about 0.15 s, the reflection effect was decreased compared with its absorption effect as the blast loading had finished. Highly concentrated stress was mainly located at the left block of the modeling domain.


#### 3.1.2. Inclined Fault Planes

To explore the effect of the orientation of a preexisting fault on the reflection and absorption effect of the stress wave, the inclined fault plane system was organized as shown in Figures3 and 4. At about 0.02 s, the stress wave in the cases did not show a significant difference because it was still propagated to fault A. At about 0.03 s, the effect of a barrier on the stress wave due to the presence of fault A started to become obvious. Because of their different orientations, case D75 had the maximum area of concentrated tensile stress (red ellipse in the right column in Figure 3). At about 0.04 s, the reflection effect becomes dominant. The reflection of the Sxx stress is highlighted with the red dotted line in Figures 3 and 4. For cases D60 and D135 (and D135R), the area of the reflection was mainly located at lower and upper intersection zones between fault A and the top and bottom boundaries of the modeling domain, respectively. For the case D75, which employed a steeper fault A, the areas of the refection of the stress wave were much wider than those observed in cases D60 and D135 (and D135R), presenting a flat ellipse in the vicinity of fault A. At about 0.09 s, the peak of the transmissive stress wave was propagated forward to fault B. The fault acted as a barrier in the stress wave propagation, but due to it not being very thick or long, the effect was much lower as compared with fault A. According to the relative spatial position of fault B and the stress wave, the propagation velocity of the stress wave in cases D60 and D75 was relatively faster than that in case D135. At about 0.15 s, the distribution of the Sxx stress became generally dislocated in most areas of the modeling domain (Figures 3 and 4). Specifically, the presence of fault A divided the box into a higher dislocation (left block) and a relatively placid area (right block).Figure 3
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 60° and 75°, respectively. The faults were indicated with black lines and purple front (A and B).Figure 4
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 135°, with different dip angles (β) of fault B (75° and 105°). The fault positions were indicated with black lines and purple front (A and B).


### 3.2. Distribution of Strain Energy Density

To explore the influence of a preexisting fault with different orientations acting on varying the distribution of the strain energy density induced by stress waves at different stages of the propagation process, three cases were selected to show the strain energy density along the horizontal middle line of the cross section in the middle of the modeling domain at 0.03 s, 0.04 s, and 0.05 s.At 0.03 s, when the stress wave induced from the blast loading did not reach the preexisting fault A, three major peaks of the strain energy density pulses were observed in the three cases (Figure5(a)). The distribution of the strain energy density shows that the maximum amplitudes are linearly increased from the left to right, generally fitting well with the inclined lines L1, L2, and L3. In case D60, the three maximum amplitudes (indicated with the purple points in Figure 5(a)) of the strain energy density were about 3052 J/m3, 3386 J/m3, and 4691 J/m3, suggesting an increasing rate of about 11% and 54% for the second and third peaks, respectively, compared with the amplitude of the first peak. In case D90, three similar peaks of the strain energy density yield amplitudes of 2568 J/m3, 3902 J/m3, and 4665 J/m3, respectively, showing an increasing rate of about 52% and 20% for the second and third peaks, respectively, as compared with amplitude of the first peak. In case D135, the strain energy density yields three maximum amplitudes peaks of about 5943 J/m3, 5705 J/m3, and 7332 J/m3, suggesting an increasing rate of about −4% and 29%, respectively, for the second and third peaks.Figure 5
The strain energy density along the horizontal middle line of the cross section in the middle of the model box (as shown in Figures2–4). The horizontal axis represents the distance from the blast loading area to the current position. The purple, blue, and black curves represent the distribution and evolution of the strain energy density at 0.03 s, 0.04 s, and 0.05 s, respectively. The preexisting fault is indicated with red solid line, the dip angle of fault A is set to α = 60°, 90°, and 135°, respectively. Dotted L1, L2, and L3 is used to indicate the evolution of the major pulses of the strain energy density with time. (a) Case D60. (b) Case D90. (c) Case D135.
(a)
(b)
(c)At about 0.04 s, the stress wave (indicated with blue solid line in Figure5) arrived at fault A. Three peaks of the strain energy density were observed and are identified with blue dots (Figure 5), generally showing a linear increase. For the three cases, the three peaks of strain energy density were observed on the left side of fault A, generally showing a flat, steeply increasing pathway for the yielded maximum strain energy density. In case D60, the three peaks (Figure 5(b)) yielded an amplitude of about 1076 J/m3, 1219 J/m3, and 3032 J/m3, and the increasing rates are about 13.2% and 149% for the second and third peaks, respectively. For case D90, the three peaks yielded an amplitude of about 1492 J/m3, 977 J/m3, and 3263 J/m3, and the maximum strain energy density reached a rate of about −35% and 234% for the second and third peaks, respectively, compared with its left peak value, representing a bathtub curve (blue lines in Figure 5(b)). In case D135, the three peaks of the strain energy density yielded amplitudes of about 1654 J/m3, 1911 J/m3, and 4178 J/m3, and the increasing rates were about 16% and 119% for the latter two peaks, respectively, compared with the first peak.At about 0.05 s, when the stress wave passed through the preexisting fault A, three peaks of the strain energy density were shown (indicated by the blue solid lines in Figure5). The trend line of these peaks presented a slow decrease (indicated by the black dotted lines in Figure 5). For case D60, the three peaks yielded amplitudes of about 1040 J/m3, 843 J/m3, and 593 J/m3, representing a decreasing rate of about −19% and −30% for the latter two peaks, respectively, compared with the first peak. In case D90, the maximum value of the strain energy density at points around the preexisting fault A yielded amplitudes of about 977 J/m3, 843 J/m3, and 623 J/m3, indicating a decreasing rate about −14% and −26%, in the second and third peaks, respectively, compared with the first one. In case D135, the amplitudes of the three peaks were about 1028 J/m3, 895 J/m3, and 825 J/m3 with decreasing rates about −13% and −5% in the second and third peaks, respectively, compared with the first peak.Comparing these three cases, it was found that when the stress wave induced by blast loading initially arrived at the preexisting fault A (at about 0.04 s), a relatively high concentration of the strain energy density was observed at the last wave before passing through fault A. The presence of a fault leads to the reflection of most of the blast energy [27]. This phenomenon indicates that the existence of fault A in the system decreases the propagation rate of the stress wave in granite rock. Thus, the strain energy carried by the stress wave is not smoothly transferred compared with the efficient transfer of energy before the stress wave reaches the fault areas. This results in a high concentration of strain energy before the stress wave reaches fault A, and fault A starts to absorb the energy. When the stress wave propagates across fault A at about 0.05 s, the strain energy stored in the stress wave attenuates; as a result, most of the strain energy is absorbed by the fault domain, as evidenced by yielding a decreasing trend line for the distribution of the strain energy density (indicated by the black dotted line in Figure 5).


### 3.3. Distribution of Strain Influenced by the Dip Angle

To examine the role of the preexisting fault with different orientations (Table1) in the evolution and distribution of strain, five cases with employing the same dip angle for fault B at different timings (0.02 s, 0.03 s, 0.04 s, 0.05 s, 0.09 s, 0.13 s, and 0.15 s) were plotted in Figure 6. The strain recorded at the monitoring lines showed that it generally decreased with time, and this was due to the linear elastic model being employed for the materials. At about 0.15 s, the strain was close to 0 with some small disturbances; for case D0, this can be regarded as a horizontal delamination. The presence of horizontal fault A showed a significant influence on the distribution of strain. The strain at the first half of fault A presented a steep increase (the tensile strain is maximum compared to other cases) followed by regular decrease. After that, the strain at the monitoring line was approaching 0. This was probably due to the reflection of the blast energy at the forebody and then a switch to the absorption stage for the energy. Regarding the interaction between the two faults during the prorogation of the stress wave, for cases D60, D75, and D90, a steep increase in compressional strain was observed in the vicinity of fault B, while for case D135, a tensional strain was observed. This was probably due to the dip orientations of the inclined faults controlling the reflectional degree of the blast energy when the stress wave arrived in the vicinity of fault A.Figure 6
The distribution of strain in theXX direction (Figure 1) at different time influenced by the dip angle of fault A. The vertical axis represents XX strain. The horizontal axis is the middle line (the green dotted line in Figure 1) of the cross section (plane γ, boundaries are indicated by the orange dotted lines) in the middle of the model box (as shown in Figures 2–4), representing the distance from the blast loading area to the current position. (a) Case D0. (b) Case D75. (c) Case D135. (d) Case D60. (e) Case D90.
(a)
(b)
(c)
(d)
(e)


## 3.1. Distribution of Stress Influenced by the Dip Angle

### 3.1.1. Horizontal and Vertical Fault Planes

Horizontal delamination and vertical fault can be used as extreme geological conditions for exploring the role of preexisting discontinuity on the propagation process of the stress wave in civil engineering, earthquake engineering, and underground mining engineering. This paper first examines the difference between the two cases D0 and D90.For case D0 (Figure2), the results show that the blast loading spreads out in a wave and it propagated forward to fault A at about 0.15–0.02 s. At about 0.02 s, the front part of fault A was surrounded by the stress wave, with a large range of tensile stress and a relatively small range of compressional stress. The difference in the Sxx stress between fault A and the surrounding granite rock body was obvious. At about 0.03 s and 0.04 s, the maximum tensile and compressional stresses was decreasing, and a dislocated feature of the stress wave was observed in a triangular area along fault A, which is different from Case D90 (indicated with light red dotted lines in the first column in Figure 2). This is probably due to the effects of natural vibration and/or absorption during the propagation process of the stress wave. At about 0.09 s, the distribution of the Sxx stress at the left part of fault A dislocated, while at the right part of fault A, the Sxx stress showed highly zoned areas of tensile and compressional stresses. In particular, the small fault B presented a truncated effect on the propagation of the stress wave. At about 0.15 s, most of the modeling domain presented a dislocated distribution of the Sxx stress, only some highly concentrated stress was observed between the two faults and in the middle part of fault A, which was probably due to the reflection effect on the stress wave as a result of the presence of fault B.Figure 2
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 0° and 90°, respectively. The fault positions were indicated with black lines and purple front (A and B).For Case D90, the stress wave was propagated in the vicinity of fault A at about 0.03 s. At about 0.04 s, fault A started to effectively act as a significant barrier in the propagation of the stress wave, presenting a highly concentrated tensile stress (as indicated as red dotted ellipse area in Figure2). At about 0.09 s, the distribution of the Sxx stress near fault A showed an obvious reflection effect according to the symmetrical semicircular characteristic of the tensile and compressional stresses. Compared with Case D0, this case was less dislocated for the distribution of the Sxx stress in the modeling domain. At about 0.15 s, the reflection effect was decreased compared with its absorption effect as the blast loading had finished. Highly concentrated stress was mainly located at the left block of the modeling domain.


### 3.1.2. Inclined Fault Planes

To explore the effect of the orientation of a preexisting fault on the reflection and absorption effect of the stress wave, the inclined fault plane system was organized as shown in Figures3 and 4. At about 0.02 s, the stress wave in the cases did not show a significant difference because it was still propagated to fault A. At about 0.03 s, the effect of a barrier on the stress wave due to the presence of fault A started to become obvious. Because of their different orientations, case D75 had the maximum area of concentrated tensile stress (red ellipse in the right column in Figure 3). At about 0.04 s, the reflection effect becomes dominant. The reflection of the Sxx stress is highlighted with the red dotted line in Figures 3 and 4. For cases D60 and D135 (and D135R), the area of the reflection was mainly located at lower and upper intersection zones between fault A and the top and bottom boundaries of the modeling domain, respectively. For the case D75, which employed a steeper fault A, the areas of the refection of the stress wave were much wider than those observed in cases D60 and D135 (and D135R), presenting a flat ellipse in the vicinity of fault A. At about 0.09 s, the peak of the transmissive stress wave was propagated forward to fault B. The fault acted as a barrier in the stress wave propagation, but due to it not being very thick or long, the effect was much lower as compared with fault A. According to the relative spatial position of fault B and the stress wave, the propagation velocity of the stress wave in cases D60 and D75 was relatively faster than that in case D135. At about 0.15 s, the distribution of the Sxx stress became generally dislocated in most areas of the modeling domain (Figures 3 and 4). Specifically, the presence of fault A divided the box into a higher dislocation (left block) and a relatively placid area (right block).Figure 3
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 60° and 75°, respectively. The faults were indicated with black lines and purple front (A and B).Figure 4
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 135°, with different dip angles (β) of fault B (75° and 105°). The fault positions were indicated with black lines and purple front (A and B).


## 3.1.1. Horizontal and Vertical Fault Planes

Horizontal delamination and vertical fault can be used as extreme geological conditions for exploring the role of preexisting discontinuity on the propagation process of the stress wave in civil engineering, earthquake engineering, and underground mining engineering. This paper first examines the difference between the two cases D0 and D90.For case D0 (Figure2), the results show that the blast loading spreads out in a wave and it propagated forward to fault A at about 0.15–0.02 s. At about 0.02 s, the front part of fault A was surrounded by the stress wave, with a large range of tensile stress and a relatively small range of compressional stress. The difference in the Sxx stress between fault A and the surrounding granite rock body was obvious. At about 0.03 s and 0.04 s, the maximum tensile and compressional stresses was decreasing, and a dislocated feature of the stress wave was observed in a triangular area along fault A, which is different from Case D90 (indicated with light red dotted lines in the first column in Figure 2). This is probably due to the effects of natural vibration and/or absorption during the propagation process of the stress wave. At about 0.09 s, the distribution of the Sxx stress at the left part of fault A dislocated, while at the right part of fault A, the Sxx stress showed highly zoned areas of tensile and compressional stresses. In particular, the small fault B presented a truncated effect on the propagation of the stress wave. At about 0.15 s, most of the modeling domain presented a dislocated distribution of the Sxx stress, only some highly concentrated stress was observed between the two faults and in the middle part of fault A, which was probably due to the reflection effect on the stress wave as a result of the presence of fault B.Figure 2
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 0° and 90°, respectively. The fault positions were indicated with black lines and purple front (A and B).For Case D90, the stress wave was propagated in the vicinity of fault A at about 0.03 s. At about 0.04 s, fault A started to effectively act as a significant barrier in the propagation of the stress wave, presenting a highly concentrated tensile stress (as indicated as red dotted ellipse area in Figure2). At about 0.09 s, the distribution of the Sxx stress near fault A showed an obvious reflection effect according to the symmetrical semicircular characteristic of the tensile and compressional stresses. Compared with Case D0, this case was less dislocated for the distribution of the Sxx stress in the modeling domain. At about 0.15 s, the reflection effect was decreased compared with its absorption effect as the blast loading had finished. Highly concentrated stress was mainly located at the left block of the modeling domain.


## 3.1.2. Inclined Fault Planes

To explore the effect of the orientation of a preexisting fault on the reflection and absorption effect of the stress wave, the inclined fault plane system was organized as shown in Figures3 and 4. At about 0.02 s, the stress wave in the cases did not show a significant difference because it was still propagated to fault A. At about 0.03 s, the effect of a barrier on the stress wave due to the presence of fault A started to become obvious. Because of their different orientations, case D75 had the maximum area of concentrated tensile stress (red ellipse in the right column in Figure 3). At about 0.04 s, the reflection effect becomes dominant. The reflection of the Sxx stress is highlighted with the red dotted line in Figures 3 and 4. For cases D60 and D135 (and D135R), the area of the reflection was mainly located at lower and upper intersection zones between fault A and the top and bottom boundaries of the modeling domain, respectively. For the case D75, which employed a steeper fault A, the areas of the refection of the stress wave were much wider than those observed in cases D60 and D135 (and D135R), presenting a flat ellipse in the vicinity of fault A. At about 0.09 s, the peak of the transmissive stress wave was propagated forward to fault B. The fault acted as a barrier in the stress wave propagation, but due to it not being very thick or long, the effect was much lower as compared with fault A. According to the relative spatial position of fault B and the stress wave, the propagation velocity of the stress wave in cases D60 and D75 was relatively faster than that in case D135. At about 0.15 s, the distribution of the Sxx stress became generally dislocated in most areas of the modeling domain (Figures 3 and 4). Specifically, the presence of fault A divided the box into a higher dislocation (left block) and a relatively placid area (right block).Figure 3
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 60° and 75°, respectively. The faults were indicated with black lines and purple front (A and B).Figure 4
The cross section in the middle of the model box (Figure1) showing the evolution and distribution of the stress in the XX direction when the dip angle (α) of fault A was set to 135°, with different dip angles (β) of fault B (75° and 105°). The fault positions were indicated with black lines and purple front (A and B).


## 3.2. Distribution of Strain Energy Density

To explore the influence of a preexisting fault with different orientations acting on varying the distribution of the strain energy density induced by stress waves at different stages of the propagation process, three cases were selected to show the strain energy density along the horizontal middle line of the cross section in the middle of the modeling domain at 0.03 s, 0.04 s, and 0.05 s.At 0.03 s, when the stress wave induced from the blast loading did not reach the preexisting fault A, three major peaks of the strain energy density pulses were observed in the three cases (Figure5(a)). The distribution of the strain energy density shows that the maximum amplitudes are linearly increased from the left to right, generally fitting well with the inclined lines L1, L2, and L3. In case D60, the three maximum amplitudes (indicated with the purple points in Figure 5(a)) of the strain energy density were about 3052 J/m3, 3386 J/m3, and 4691 J/m3, suggesting an increasing rate of about 11% and 54% for the second and third peaks, respectively, compared with the amplitude of the first peak. In case D90, three similar peaks of the strain energy density yield amplitudes of 2568 J/m3, 3902 J/m3, and 4665 J/m3, respectively, showing an increasing rate of about 52% and 20% for the second and third peaks, respectively, as compared with amplitude of the first peak. In case D135, the strain energy density yields three maximum amplitudes peaks of about 5943 J/m3, 5705 J/m3, and 7332 J/m3, suggesting an increasing rate of about −4% and 29%, respectively, for the second and third peaks.Figure 5
The strain energy density along the horizontal middle line of the cross section in the middle of the model box (as shown in Figures2–4). The horizontal axis represents the distance from the blast loading area to the current position. The purple, blue, and black curves represent the distribution and evolution of the strain energy density at 0.03 s, 0.04 s, and 0.05 s, respectively. The preexisting fault is indicated with red solid line, the dip angle of fault A is set to α = 60°, 90°, and 135°, respectively. Dotted L1, L2, and L3 is used to indicate the evolution of the major pulses of the strain energy density with time. (a) Case D60. (b) Case D90. (c) Case D135.
(a)
(b)
(c)At about 0.04 s, the stress wave (indicated with blue solid line in Figure5) arrived at fault A. Three peaks of the strain energy density were observed and are identified with blue dots (Figure 5), generally showing a linear increase. For the three cases, the three peaks of strain energy density were observed on the left side of fault A, generally showing a flat, steeply increasing pathway for the yielded maximum strain energy density. In case D60, the three peaks (Figure 5(b)) yielded an amplitude of about 1076 J/m3, 1219 J/m3, and 3032 J/m3, and the increasing rates are about 13.2% and 149% for the second and third peaks, respectively. For case D90, the three peaks yielded an amplitude of about 1492 J/m3, 977 J/m3, and 3263 J/m3, and the maximum strain energy density reached a rate of about −35% and 234% for the second and third peaks, respectively, compared with its left peak value, representing a bathtub curve (blue lines in Figure 5(b)). In case D135, the three peaks of the strain energy density yielded amplitudes of about 1654 J/m3, 1911 J/m3, and 4178 J/m3, and the increasing rates were about 16% and 119% for the latter two peaks, respectively, compared with the first peak.At about 0.05 s, when the stress wave passed through the preexisting fault A, three peaks of the strain energy density were shown (indicated by the blue solid lines in Figure5). The trend line of these peaks presented a slow decrease (indicated by the black dotted lines in Figure 5). For case D60, the three peaks yielded amplitudes of about 1040 J/m3, 843 J/m3, and 593 J/m3, representing a decreasing rate of about −19% and −30% for the latter two peaks, respectively, compared with the first peak. In case D90, the maximum value of the strain energy density at points around the preexisting fault A yielded amplitudes of about 977 J/m3, 843 J/m3, and 623 J/m3, indicating a decreasing rate about −14% and −26%, in the second and third peaks, respectively, compared with the first one. In case D135, the amplitudes of the three peaks were about 1028 J/m3, 895 J/m3, and 825 J/m3 with decreasing rates about −13% and −5% in the second and third peaks, respectively, compared with the first peak.Comparing these three cases, it was found that when the stress wave induced by blast loading initially arrived at the preexisting fault A (at about 0.04 s), a relatively high concentration of the strain energy density was observed at the last wave before passing through fault A. The presence of a fault leads to the reflection of most of the blast energy [27]. This phenomenon indicates that the existence of fault A in the system decreases the propagation rate of the stress wave in granite rock. Thus, the strain energy carried by the stress wave is not smoothly transferred compared with the efficient transfer of energy before the stress wave reaches the fault areas. This results in a high concentration of strain energy before the stress wave reaches fault A, and fault A starts to absorb the energy. When the stress wave propagates across fault A at about 0.05 s, the strain energy stored in the stress wave attenuates; as a result, most of the strain energy is absorbed by the fault domain, as evidenced by yielding a decreasing trend line for the distribution of the strain energy density (indicated by the black dotted line in Figure 5).


## 3.3. Distribution of Strain Influenced by the Dip Angle

To examine the role of the preexisting fault with different orientations (Table1) in the evolution and distribution of strain, five cases with employing the same dip angle for fault B at different timings (0.02 s, 0.03 s, 0.04 s, 0.05 s, 0.09 s, 0.13 s, and 0.15 s) were plotted in Figure 6. The strain recorded at the monitoring lines showed that it generally decreased with time, and this was due to the linear elastic model being employed for the materials. At about 0.15 s, the strain was close to 0 with some small disturbances; for case D0, this can be regarded as a horizontal delamination. The presence of horizontal fault A showed a significant influence on the distribution of strain. The strain at the first half of fault A presented a steep increase (the tensile strain is maximum compared to other cases) followed by regular decrease. After that, the strain at the monitoring line was approaching 0. This was probably due to the reflection of the blast energy at the forebody and then a switch to the absorption stage for the energy. Regarding the interaction between the two faults during the prorogation of the stress wave, for cases D60, D75, and D90, a steep increase in compressional strain was observed in the vicinity of fault B, while for case D135, a tensional strain was observed. This was probably due to the dip orientations of the inclined faults controlling the reflectional degree of the blast energy when the stress wave arrived in the vicinity of fault A.Figure 6
The distribution of strain in theXX direction (Figure 1) at different time influenced by the dip angle of fault A. The vertical axis represents XX strain. The horizontal axis is the middle line (the green dotted line in Figure 1) of the cross section (plane γ, boundaries are indicated by the orange dotted lines) in the middle of the model box (as shown in Figures 2–4), representing the distance from the blast loading area to the current position. (a) Case D0. (b) Case D75. (c) Case D135. (d) Case D60. (e) Case D90.
(a)
(b)
(c)
(d)
(e)


## 4. Discussion

### 4.1. Stress Wave Propagation with Time

The interaction between the stress wave and the preexisting faults play an important role in estimating the mechanical conditions and timing for reactivation of the faults [28–31]. During the propagation of a dynamic/blast wave, part of the stress wave energy will be absorbed by the preexisting faults, resulting in an increase of energy before it arrives at the faults and the attenuation of energy after crossing the faults [32–34]. The stress wave probably originates from ruptures of rock and coal masses, breakage, and/or collapse of the roof or blast loading in the stope. The energy disturbance along and through the preexisting faults has a large influence on the evolution and distribution of dynamic disasters such as rockbursts and the sudden reactivation of faults in coal mining [33].To explore the evolutional process and attenuation effect on stress influenced by the existence of a fault while imposing the stress wave in the modeling domain (Figure1), six passive particles were placed along faults A and B (Figure 7). Regarding the location of the six passive particles (panel (d) in Figure 7), four were grouped into two groups of two and placed at the left and right blocks divided by fault A, respectively. With respect to constraining the related parameters with time around fault B, another two particles were placed at its two blocks, one per block. The result shows that, generally, the existence of faults A and B effectively attenuates the stress waves. We employ ΔSxx (orange dotted line in Figure 7) to show the stress difference of the amplitudes of the first and/or second waves, indicating the attenuation degree due to the absorption or resistance effect of weak faults.Figure 7
The evolution and attenuation of stress in theXX direction (Sxx, MPa) recorded at six passive particles (indicated by the black and red points, faults A and B are shown by the purple lines) influenced by the stress wave. Panel (d) is a cross section extracted in the middle of the model box (Figure 1). The Y-axis represents Sxx (MPa), and the X-axis is the duration of running model. ΔSxx is used to show the difference of the first and/or second waves’ amplitude (orange dotted line), which could help indicate the attenuation degree due to the absorption or resistance effect of weak faults.For panel (a), the two passive particles first experienced a pressure wave, followed by a relatively weaker tensile wave. Regarding the first crest of the pressure wave, recorded at the two passive particles, the fault showed a significant effect of absorption. Sxx, obtained at the right side of fault A, is decreased by 91.4% compared with the state of the passive particles located at the left side of fault A. With respect to the first trough of the stress wave, the effect of absorption of fault A was relatively weak; Sxx was decreased by 69% as a result of fault A.Regarding the passive particles in panel (b), they have a shortest distance to the blast loading areas compared with other passive particles analyzed in panels (a) and (c). For the first trough of the stress wave obtained at two particles, the particles yielded a maximum and minimum Sxx of −26.3 MPa and −4.5 MPa, respectively, with the second values decreasing by 82.9% because of the existence of fault A. The first crest of the stress wave, yielded about 5.5 MPa and 4.9 MPa with a ΔSxx of 0.6 MPa.In panel (c), the two particles can be used to explore the role of fault B in the absorption effect over the propagation process of the stress wave. At about 0.75 s, the black particle recorded the trough, which yielded values of about −1.2 MPa and −1.1 MPa. After that, a relatively high crest of the stress wave obtained by the black particle at about 0.85 s yielded amplitudes of about 3.9 MPa and 1.8 MPa with a ΔSxx of 2.1 MPa.By comparing the evolution and the effect of absorption of the stress wave, it can be found that a preexisting fault generally has a significant role in varying the amplitude/energy of the stress wave [27, 35, 36] owing to the difference in elastic modulus, density, and Poisson’s ratio [37–40] between the two materials. This indicates that a thick fault plane (fault A) has a greater effect on reducing the maximum amplitudes (Figure 6) compared with a thin fault plane (fault B).


### 4.2. Implications for Chaoyang Coal Mine

Chaoyang Coal Mine is located in Shandong Province, China. The mine is about 12 km long from east to west and 1.5 km wide from north to south. The mining depth is about −550 m to −1200 m. The 3206 working face is located at the boundary of the mine and the working depth is about −1000 m. The dip angle of the coal seam is about 16°. The thickness of the coal seam is 2.0 m to 7.5 m. Regarding the working face, the average length of the strike is 303 m, and the tendency length is about 111 m. Based on the Chinese standards GB/T 25217.1 and GB/T 25217.2 for assessing the risk of rockbursts of coal and/or rock mass, the coal seam has a weak impact tendency for the occurrence of rockbursts. The side of the working face is 3205 goaf, and the 3205 and 3206 working faces of the longwall panel are adjacent to the Fushan fault (300 m drop, 30°–50° inclination) (Figure8). A fault coal pillar with a width of about 100 m is left between the 3206 working face and the Laoshan fault, and a small fault has developed in the middle of the 3206-working face (F3206-1 fault, 0–10 m drop, and 70° inclination). Because of fault F3206-1, severe mine pressure appeared during the initial mining of the working face, and strong mine earthquakes appear frequently. The closer the mining distance is to fault F3206-1, the more serious the phenomenon is.Figure 8
Seismic events between February 2, 2013, and September 16, 2013 (modified from [33]). Two different color dots are employed to show the energy of seismic events.In Figure8, F3206-1 fault divides longwall panel 3206 into two blocks. Regarding the distribution of the seismic events that occurred between February 2, 2013, and September 16, 2013, the two grades of seismic events (>10000 J and 1000∼10000 J) appeared along fault F3206-1. In particular, most of them were concentrated at the footwall of the fault. This regular distribution of seismic events indicates that the strain energy carried by stress waves induced by blast loading, mining activities, and fracture of roofs (which happened at the left blocks of fault F3206-1) is not transferred smoothly through fault F3206-1, resulting in a high concentration of strain energy (Figure 5) before the stress wave reaches fault F3206-1. This high concentration of strain energy results in fractures of rock and coal masses and/or mining earthquakes at the footwall of fault F3206-1, evidenced by the large amounts of seismic events at the footwall. When the stress wave fully passes through fault F3206-1, most of the strain energy is absorbed by fault F3206-1, suggesting that less strain energy is concentrated at the hanging wall compared with that at the footwall. This inference about the uneven distribution of the strain energy carried by the strain waves can be made evident by the distribution of the seismic events as indicated by red and blue dots in Figure 8.


## 4.1. Stress Wave Propagation with Time

The interaction between the stress wave and the preexisting faults play an important role in estimating the mechanical conditions and timing for reactivation of the faults [28–31]. During the propagation of a dynamic/blast wave, part of the stress wave energy will be absorbed by the preexisting faults, resulting in an increase of energy before it arrives at the faults and the attenuation of energy after crossing the faults [32–34]. The stress wave probably originates from ruptures of rock and coal masses, breakage, and/or collapse of the roof or blast loading in the stope. The energy disturbance along and through the preexisting faults has a large influence on the evolution and distribution of dynamic disasters such as rockbursts and the sudden reactivation of faults in coal mining [33].To explore the evolutional process and attenuation effect on stress influenced by the existence of a fault while imposing the stress wave in the modeling domain (Figure1), six passive particles were placed along faults A and B (Figure 7). Regarding the location of the six passive particles (panel (d) in Figure 7), four were grouped into two groups of two and placed at the left and right blocks divided by fault A, respectively. With respect to constraining the related parameters with time around fault B, another two particles were placed at its two blocks, one per block. The result shows that, generally, the existence of faults A and B effectively attenuates the stress waves. We employ ΔSxx (orange dotted line in Figure 7) to show the stress difference of the amplitudes of the first and/or second waves, indicating the attenuation degree due to the absorption or resistance effect of weak faults.Figure 7
The evolution and attenuation of stress in theXX direction (Sxx, MPa) recorded at six passive particles (indicated by the black and red points, faults A and B are shown by the purple lines) influenced by the stress wave. Panel (d) is a cross section extracted in the middle of the model box (Figure 1). The Y-axis represents Sxx (MPa), and the X-axis is the duration of running model. ΔSxx is used to show the difference of the first and/or second waves’ amplitude (orange dotted line), which could help indicate the attenuation degree due to the absorption or resistance effect of weak faults.For panel (a), the two passive particles first experienced a pressure wave, followed by a relatively weaker tensile wave. Regarding the first crest of the pressure wave, recorded at the two passive particles, the fault showed a significant effect of absorption. Sxx, obtained at the right side of fault A, is decreased by 91.4% compared with the state of the passive particles located at the left side of fault A. With respect to the first trough of the stress wave, the effect of absorption of fault A was relatively weak; Sxx was decreased by 69% as a result of fault A.Regarding the passive particles in panel (b), they have a shortest distance to the blast loading areas compared with other passive particles analyzed in panels (a) and (c). For the first trough of the stress wave obtained at two particles, the particles yielded a maximum and minimum Sxx of −26.3 MPa and −4.5 MPa, respectively, with the second values decreasing by 82.9% because of the existence of fault A. The first crest of the stress wave, yielded about 5.5 MPa and 4.9 MPa with a ΔSxx of 0.6 MPa.In panel (c), the two particles can be used to explore the role of fault B in the absorption effect over the propagation process of the stress wave. At about 0.75 s, the black particle recorded the trough, which yielded values of about −1.2 MPa and −1.1 MPa. After that, a relatively high crest of the stress wave obtained by the black particle at about 0.85 s yielded amplitudes of about 3.9 MPa and 1.8 MPa with a ΔSxx of 2.1 MPa.By comparing the evolution and the effect of absorption of the stress wave, it can be found that a preexisting fault generally has a significant role in varying the amplitude/energy of the stress wave [27, 35, 36] owing to the difference in elastic modulus, density, and Poisson’s ratio [37–40] between the two materials. This indicates that a thick fault plane (fault A) has a greater effect on reducing the maximum amplitudes (Figure 6) compared with a thin fault plane (fault B).


## 4.2. Implications for Chaoyang Coal Mine

Chaoyang Coal Mine is located in Shandong Province, China. The mine is about 12 km long from east to west and 1.5 km wide from north to south. The mining depth is about −550 m to −1200 m. The 3206 working face is located at the boundary of the mine and the working depth is about −1000 m. The dip angle of the coal seam is about 16°. The thickness of the coal seam is 2.0 m to 7.5 m. Regarding the working face, the average length of the strike is 303 m, and the tendency length is about 111 m. Based on the Chinese standards GB/T 25217.1 and GB/T 25217.2 for assessing the risk of rockbursts of coal and/or rock mass, the coal seam has a weak impact tendency for the occurrence of rockbursts. The side of the working face is 3205 goaf, and the 3205 and 3206 working faces of the longwall panel are adjacent to the Fushan fault (300 m drop, 30°–50° inclination) (Figure8). A fault coal pillar with a width of about 100 m is left between the 3206 working face and the Laoshan fault, and a small fault has developed in the middle of the 3206-working face (F3206-1 fault, 0–10 m drop, and 70° inclination). Because of fault F3206-1, severe mine pressure appeared during the initial mining of the working face, and strong mine earthquakes appear frequently. The closer the mining distance is to fault F3206-1, the more serious the phenomenon is.Figure 8
Seismic events between February 2, 2013, and September 16, 2013 (modified from [33]). Two different color dots are employed to show the energy of seismic events.In Figure8, F3206-1 fault divides longwall panel 3206 into two blocks. Regarding the distribution of the seismic events that occurred between February 2, 2013, and September 16, 2013, the two grades of seismic events (>10000 J and 1000∼10000 J) appeared along fault F3206-1. In particular, most of them were concentrated at the footwall of the fault. This regular distribution of seismic events indicates that the strain energy carried by stress waves induced by blast loading, mining activities, and fracture of roofs (which happened at the left blocks of fault F3206-1) is not transferred smoothly through fault F3206-1, resulting in a high concentration of strain energy (Figure 5) before the stress wave reaches fault F3206-1. This high concentration of strain energy results in fractures of rock and coal masses and/or mining earthquakes at the footwall of fault F3206-1, evidenced by the large amounts of seismic events at the footwall. When the stress wave fully passes through fault F3206-1, most of the strain energy is absorbed by fault F3206-1, suggesting that less strain energy is concentrated at the hanging wall compared with that at the footwall. This inference about the uneven distribution of the strain energy carried by the strain waves can be made evident by the distribution of the seismic events as indicated by red and blue dots in Figure 8.


## 5. Conclusion

It is widely accepted that exploring the propagation process of the stress wave in a rock mass with preexisting discontinuities is important to reveal rock and geological engineering problems, particularly in understanding the mechanisms of landslides, earthquakes, and dynamic disasters in underground coal mining. In this study, we explored the influence of two preexisting faults with different orientations on the propagation process of the explosion-induced stress wave and the reflection effect by building six 3D models using COMSOL Multiphysics. Some results were achieved which are as follows:(1)
For a horizontal delamination, the maximum tensile and compressional stresses decreased with time, and a dislocated feature of the stress wave was observed in a triangular area along fault A due to the effects of natural vibration and/or absorption during the propagation process of the stress wave. At about 0.15 s, most of the modeling domain presented a dislocated distribution of the Sxx stress.(2)
For an inclined fault system, it was found that when the stress wave induced by blast loading initially arrived at the preexisting fault A (at about 0.04 s), and a relatively high concentration of the strain energy density was observed at the last wave before passing through fault A. The presence of a fault leads to the reflection of most of the blast energy. This phenomenon indicates that the existence of fault A in the system decreased the propagation rate of the stress wave in granite rock; the strain energy carried by the stress wave was not smoothly transferred compared with the efficient energy transfer before the stress wave reached the fault areas. When the stress wave propagated across fault A at about 0.05 s, the strain energy stored in the stress wave attenuated. As a result, most of the strain energy was absorbed by the fault domain.(3)
The preexisting fault generally had a significant role in varying the amplitude/energy of the stress wave owing to the difference in the elastic module, density, and Poisson’s ratio between the two materials (fault and surrounding rocks).(4)
Finally, the modeling results were implicated in Chaoyang Coal Mine to account for the distribution of observed seismic events.



---

*Source: 5856080-2019-11-11.xml*