# Study on the Fracture Behavior of Cracks Emanating from Tunnel Spandrel under Blasting Loads by Using TMCSC Specimens

Authors: Bang Liu; Zheming Zhu; Ruifeng Liu; Lei Zhou; Duanying Wan
Journal: Shock and Vibration
 (2019)

---

## Abstract

Radial cracks may exist around tunnel edge, and these cracks may propagate and weaken tunnel stability under nearby blasting operations. In order to study the blast-induced fracture behavior of radial cracks emanating from a tunnel spandrel, a tunnel model containing a spandrel crack (TMCSC) with different inclination angles was proposed in this paper. Crack propagation gauges (CPGs) and strain gauges were used in the experiments to measure crack initiation moment and propagation time. Finite difference models were established by using AUTODYN code to simulate crack propagation behavior and propagation path. ABAQUS code was used to calculate dynamic stress intensity factors (SIFs). The results show that (1) crack inclination angles affect crack initiation angles and crack propagation lengths significantly; (2) critical SIFs of both mode I and mode II decrease gradually with the increase of the crack propagation speed; (3) the dynamic energy release rates vary during crack propagation; and (4) there are “crack arrest points” on the crack propagation paths in which the crack propagation speed is very small.


---

## Body

## 1. Introduction

Blasting method is often used in rock engineering, such as in mining and tunneling. There are many cracks emanating from tunnel edge which could be either induced by the tunnel blasting excavation or naturally existing. These cracks under nearby blasting operations may propagate and coalesce, damaging tunnel strength and weakening tunnel stability, which could cause some catastrophes, such as rockburst and coal-gas outburst. So, it is necessary to carry out the corresponding numerical and experimental studies.Currently, many scholars have studied the dynamic propagation behavior of the cracks under blasting and have made significant achievements. Yang et al. [1] investigated the dynamic behavior of jointed rock mass under blasting at high-stress conditions. Yue et al. [2] studied crack propagation in PMMA material with double holes under the directional controlled blasting. Chen et al. [3, 4] investigated the propagation characteristics of vibration waves induced in surrounding rock by tunneling blasting. By using the blasting technique, Verma et al. [5] studied the blasting damage of a large-scale rock mass quality with a fault tunnel, and the results were verified by the ultrasonic testing of core samples obtained from a tunnel site experiment. By the model with precracks on the dome, Masi et al. [6] simulated the explosion process, and the simulation results were similar to cracks produced in the process of damage evolution of no crack dome. Ozcelik [7] used reverse analysis method to study the cause of the expansion of cracks in buildings, which was the blasting construction of limestone quarry. With the aid of extended finite element, Yan et al. [8] studied the cracking and failure characteristics of the segmented lining structure of underwater shield tunnel under the impact of the high-speed derailment train. The lining structure of impact section leads to continuous and perforated section v cracks, and a continuous strip or polygonal crack formed around the impact zone of the joint and on continuous surface of the outer surface near the joint. Keys and Clubley [9] studied the failure mode and debris distribution of the brick slab under the explosive load and pointed out that the explosion overpressure pulse was the key parameter of the fracture failure mode, the initial breakage, and the debris distribution.Using dynamic load generated by explosives on the supporting system, Shirzadegan et al. [10] carried out the fracture study in the LKAB Kiirunavaara mine in Sweden and pointed out that the level of induced damage was limited to the fracture zone behind the supporting system and the crack propagation in the shotcrete. Slavko and Veljko [11] used the excavation process of hard rock with two joints to study the formation of the rock wedge and pointed out that the rock wedge could not be formed under the finite persistent shock wave produced by the blasting load. Stolz and Ruiz-Ripoll [12] studied the behavior of the buried tunnel section under the explosive load and pointed out that the actual loading force would cause large area cracking of the tunnel system, and the soil settlement could reach 1 cm. Li et al. [13] researched the essential reason of two empty holes on crack fracture rule under blasting and pointed out that the empty holes brought a hindrance influence on the outgoing cracks, and the arrest depended on the spacing between two empty holes. The smaller the distance between the two holes was, the better the crack arrest was. Liu et al. [14] studied the dynamic fracture toughness and other fracture parameters of type I crack under the explosive load and pointed out that there was an inflection point on the propagation path and the obvious crack arrest at the inflection point. For rock dynamic fracture toughness of mode I crack under blasting loads, Zhu et al. [15] put forward a new measuring method. At the same time, in order to calculate the displacements near the crack tip, the 1/4 node element from ANSYS code is used.Although some experts and scholars have studied the propagation behavior of cracks under explosive loading and have made some achievements, they have focused more on mode I crack, and less on mode I-II crack. Of course, some unknown or partly unknown aspects have not been studied, such as the measuring method of critical dynamic SIFs and the fracture rule of mode I-II crack under blasting. For this paper, the aspects will be researched. In terms of the simulation study of rock dynamic fracture rule, finite difference AUTODYN code has been largely used, and moreover, its feasibility has been widely proved [16–26]. Therefore, in the simulation study of crack fracture rule, it will be used. In terms of the calculation of crack dynamic stress intensity factors (DSIFs) and combining with the measurement results of crack initiation time, ABAQUS code will be used, and the critical DSIFs will be calculated.


## 2. Experimental Study

To study fracture rule of spandrel cracks under blasting, a tunnel model containing a spandrel crack (TMCSC) with different inclination angles was proposed.


### 2.1. TMCSC Specimens

The TMCSC specimen dimension used in this study is shown in Figure1. The radius of the borehole was 3.5 mm, and the 8# mine detonator (approximately 0.9 g of explosive charge) was placed in the borehole. The tunnel was made in accordance with the proportion of 1 : 100. The crack inclination angle θ was designed from −75° to 75°, the increment was 15°, and totally, there were eleven group specimens. The sandstone TMCSC specimens were cut by high-speed waterjet. The length of the tunnel spandrel cracks was 40 mm, and it was sharpened artificially by the 0.5 mm thickness steel saw blade. The thickness of the specimen was 10 mm.Figure 1
Dimension of the TMCSC specimen (unit: mm).An ultrasonic speed testing system, Sonic Viewer-SX, as shown in Figure2, which can be used to measure rock P-wave and S-wave speeds with high precision, was employed to measure the green sandstone dynamic parameters. The built-in software can calculate the parameters of dynamic Poisson’s ratio and dynamic elastic modulus. The dynamic mechanical parameters of green sandstone were obtained by the testing system, and the results are presented in Table 1.Figure 2
Ultrasonic speed testing system.Table 1
Dynamic mechanical parameters of green sandstone.
Density (kg/m3)
P-wave velocityCd (km/s)
S-wave velocityCs (km/s)
Rayleigh wave velocityCR (m/s)
Dynamic elastic modulusEd (GPa)
Tensile strengthσt (MPa)
Compressive strengthσc (MPa)
Dynamic Poisson ratioνd
2370
2.43
1.607
1.458
12.5
25
22.08
0.26


### 2.2. Rationality of TMCSC Specimen Dimension

Before blasting tests, it was necessary to verify if the reflection tensile wave from the free boundary affected the crack fracture behavior in the CPG zone concerned. Because the tensile strength of rock material is lower than its compressive strength, with the arrival of the reflected tensile stress wave, it will significantly influence dynamic behavior of the crack tip. The principle was that, as the wave reached the crack tip concerned, the initiation or propagation behavior of the crack should have completed.The shortest time of stress wave travelling along the path CDA was set asTCDA, the initiation time of crack tip A or the breaking time of the first filament of CPG was set as TA, the shortest time of stress wave travelling along the path CEB was set as TCEB, and the breaking time of the last filament of the CPG was set as TB, as shown in Figure 3. All the values of these parameters are listed in Table 2.Figure 3
Rationality analysis of TMCSC specimens (θ = −75° to 75°) dimensions.Table 2
Stress wave propagation time and crack propagation time.
Crack angle (°)
T
CDA (μs)
T
A (μs)
T
CEB (μs)
T
B (μs)
−75
321.81
64.6
326.75
168.2
−60
323.87
68.2
329.22
169.6
−45
327.57
71.8
332.92
188.6
−30
331.69
74.6
337.04
177.4
−15
335.39
75.8
341.15
181.4
0
331.69
77.4
337.45
187.0
15
327.57
76.2
333.33
181.0
30
323.05
75.6
328.81
175.4
45
319.75
72.4
325.05
164.6
60
316.05
70.6
321.40
168.8
75
313.58
67.4
318.93
166.8
T
CDA and TCEB were obtained from geometric calculation, and TA and TB were derived from monitoring data of CPG in subsequent blasting tests.From Table2, one can find all the time TA was much less than TCDA, that means the crack initiation behavior was much earlier than the arrive time of the reflected tensile wave from the upper boundary, and the initiation of the crack was not affected by the wave. In addition, time TCEB was much larger than the time TB; namely, the reflected tensile wave from the upper boundary had no effect on the crack propagation. Therefore, the TMCSC specimen dimension was quite large, the reflected tensile stress wave will not affect the dynamic behavior of the crack, and the data obtained from the experiments were effective and reliable.


### 2.3. Measuring System

In this experiment, the 8# mine detonator (approximately 0.9 g of explosive charge) was placed in the blast hole center, and the ultradynamic resistance strain gauge of type CS-1D and oscilloscope of type DS1104 were used. By a crack propagation gauge (CPG) pasted at crack tip, the initiation moment and the propagation time of the cracks were obtained. The explosion load was tested by the strain gauges near the blast hole, and the test system of the experiments is shown in Figure4.Figure 4
The test system of the blasting experiments.The strains and the associated stresses were obtained by the radial and tangential strain gauges pasted near the borehole. The CPG has 21 Kama copper wires, and the initial total resistance of the wires was 3.5 Ω. ResistorR1 (50 Ω) was connected in parallel with the CPG. In order to avoid large voltage, R2 (50 Ω) was connected in series with the CPG, and a 16 V constant power was supplied (accuracy of 1 mv). Before the strain gauges and CPGs were pasted, the specimen surface was polished by a sand paper.


### 2.4. Measurement of Blast-Induced Pressures

The thickness of the specimens was set to 15 mm, a detonator was placed in the borehole, the main charge was placed in the center of the borehole, and no filling and coupling were applied, as shown in Figure5.Figure 5
A TMCSC specimen with two strain gauges and a CPG and a detonator.In order to confirm the location of the strain gauges pasted near the borehole, the preliminary tests were conducted, and the results showed that the radius of the crushed zone was less than 25 mm. Therefore, to calculate the pressure near the borehole and to avoid the impact of crushed zone, two strain gauges were pasted at 30 mm distance to the center of the borehole, as shown in Figure5.The tangential strain (εθ) and the radial strain (εr) obtained by the tangential and radial strain gauges, respectively, are shown in Figure 6.Figure 6
The relation curves between strain and time.From the theory of elastic mechanics, the pressurePt can be calculated by(1)Pt=E1−ν2εrt+ν⋅εθt,where ν stands for Poisson’s ratio and E stands for the elastic modulus.By Figure6 and equation (1), the pressures Pt were calculated, and the results are shown in Figure 7. For the sandstone specimen, the maximum value of blast-induced pressure was 62.9 MPa from Figure 7, which was much less than the pressure on the borehole wall [27].Figure 7
The pressure curve.


### 2.5. Propagation Behavior of the Cracks

CPG was pasted on the propagation paths of the cracks to further monitor the dynamic fracture of the cracks, and the oscilloscope was used to display and collect the voltage data of the CPG. Eleven groups of specimens were tested. The crack inclination angles changed from −75° to 75° for different groups to analyze the dynamic fracture of the cracks under blasting.As the crack propagated, the CPG wires were broken consecutively. The relationship of voltage signal with time was a step-shaped curve. The measurement results of sandstone specimens with crack inclination angleθ = −45° and 45° were plotted and are shown in Figures 8(a) and 8(b). The crack speed was calculated, and the results are shown in Figure 8(c). Meanwhile, the average crack propagation speeds are shown in Figure 8(d).Figure 8
The voltage signal recorded by CPG and the propagation speed of the cracks. Voltage signal curve of specimen with inclination angle (a) −45° and (b) 45°. Relation between (c) crack speed and crack length and (d) average crack speed versus crack length.
(a)
(b)
(c)
(d)From Figures8(a) and 8(b), one can find that, between the 5th and 7th wire and between the 7th and 9th wire, the crack speeds spent a relatively long time to propagate, and the propagation speed of crack inclination angle −45° specimen shown in Figure 8(c) was only 80.65 m/s which was much lower than the average speed 342 m/s shown in Figure 8(d), and similarly, for other specimens of crack inclination angle 45° shown in Figure 8(c), the propagation speed was 128.21 m/s which was less than the average speed 434 m/s, as shown in Figure 8(d). When θ equaled 45°, the average crack speed was the maximum, and when θ equaled −45°, the average crack speed was the minimum.


## 2.1. TMCSC Specimens

The TMCSC specimen dimension used in this study is shown in Figure1. The radius of the borehole was 3.5 mm, and the 8# mine detonator (approximately 0.9 g of explosive charge) was placed in the borehole. The tunnel was made in accordance with the proportion of 1 : 100. The crack inclination angle θ was designed from −75° to 75°, the increment was 15°, and totally, there were eleven group specimens. The sandstone TMCSC specimens were cut by high-speed waterjet. The length of the tunnel spandrel cracks was 40 mm, and it was sharpened artificially by the 0.5 mm thickness steel saw blade. The thickness of the specimen was 10 mm.Figure 1
Dimension of the TMCSC specimen (unit: mm).An ultrasonic speed testing system, Sonic Viewer-SX, as shown in Figure2, which can be used to measure rock P-wave and S-wave speeds with high precision, was employed to measure the green sandstone dynamic parameters. The built-in software can calculate the parameters of dynamic Poisson’s ratio and dynamic elastic modulus. The dynamic mechanical parameters of green sandstone were obtained by the testing system, and the results are presented in Table 1.Figure 2
Ultrasonic speed testing system.Table 1
Dynamic mechanical parameters of green sandstone.
Density (kg/m3)
P-wave velocityCd (km/s)
S-wave velocityCs (km/s)
Rayleigh wave velocityCR (m/s)
Dynamic elastic modulusEd (GPa)
Tensile strengthσt (MPa)
Compressive strengthσc (MPa)
Dynamic Poisson ratioνd
2370
2.43
1.607
1.458
12.5
25
22.08
0.26


## 2.2. Rationality of TMCSC Specimen Dimension

Before blasting tests, it was necessary to verify if the reflection tensile wave from the free boundary affected the crack fracture behavior in the CPG zone concerned. Because the tensile strength of rock material is lower than its compressive strength, with the arrival of the reflected tensile stress wave, it will significantly influence dynamic behavior of the crack tip. The principle was that, as the wave reached the crack tip concerned, the initiation or propagation behavior of the crack should have completed.The shortest time of stress wave travelling along the path CDA was set asTCDA, the initiation time of crack tip A or the breaking time of the first filament of CPG was set as TA, the shortest time of stress wave travelling along the path CEB was set as TCEB, and the breaking time of the last filament of the CPG was set as TB, as shown in Figure 3. All the values of these parameters are listed in Table 2.Figure 3
Rationality analysis of TMCSC specimens (θ = −75° to 75°) dimensions.Table 2
Stress wave propagation time and crack propagation time.
Crack angle (°)
T
CDA (μs)
T
A (μs)
T
CEB (μs)
T
B (μs)
−75
321.81
64.6
326.75
168.2
−60
323.87
68.2
329.22
169.6
−45
327.57
71.8
332.92
188.6
−30
331.69
74.6
337.04
177.4
−15
335.39
75.8
341.15
181.4
0
331.69
77.4
337.45
187.0
15
327.57
76.2
333.33
181.0
30
323.05
75.6
328.81
175.4
45
319.75
72.4
325.05
164.6
60
316.05
70.6
321.40
168.8
75
313.58
67.4
318.93
166.8
T
CDA and TCEB were obtained from geometric calculation, and TA and TB were derived from monitoring data of CPG in subsequent blasting tests.From Table2, one can find all the time TA was much less than TCDA, that means the crack initiation behavior was much earlier than the arrive time of the reflected tensile wave from the upper boundary, and the initiation of the crack was not affected by the wave. In addition, time TCEB was much larger than the time TB; namely, the reflected tensile wave from the upper boundary had no effect on the crack propagation. Therefore, the TMCSC specimen dimension was quite large, the reflected tensile stress wave will not affect the dynamic behavior of the crack, and the data obtained from the experiments were effective and reliable.


## 2.3. Measuring System

In this experiment, the 8# mine detonator (approximately 0.9 g of explosive charge) was placed in the blast hole center, and the ultradynamic resistance strain gauge of type CS-1D and oscilloscope of type DS1104 were used. By a crack propagation gauge (CPG) pasted at crack tip, the initiation moment and the propagation time of the cracks were obtained. The explosion load was tested by the strain gauges near the blast hole, and the test system of the experiments is shown in Figure4.Figure 4
The test system of the blasting experiments.The strains and the associated stresses were obtained by the radial and tangential strain gauges pasted near the borehole. The CPG has 21 Kama copper wires, and the initial total resistance of the wires was 3.5 Ω. ResistorR1 (50 Ω) was connected in parallel with the CPG. In order to avoid large voltage, R2 (50 Ω) was connected in series with the CPG, and a 16 V constant power was supplied (accuracy of 1 mv). Before the strain gauges and CPGs were pasted, the specimen surface was polished by a sand paper.


## 2.4. Measurement of Blast-Induced Pressures

The thickness of the specimens was set to 15 mm, a detonator was placed in the borehole, the main charge was placed in the center of the borehole, and no filling and coupling were applied, as shown in Figure5.Figure 5
A TMCSC specimen with two strain gauges and a CPG and a detonator.In order to confirm the location of the strain gauges pasted near the borehole, the preliminary tests were conducted, and the results showed that the radius of the crushed zone was less than 25 mm. Therefore, to calculate the pressure near the borehole and to avoid the impact of crushed zone, two strain gauges were pasted at 30 mm distance to the center of the borehole, as shown in Figure5.The tangential strain (εθ) and the radial strain (εr) obtained by the tangential and radial strain gauges, respectively, are shown in Figure 6.Figure 6
The relation curves between strain and time.From the theory of elastic mechanics, the pressurePt can be calculated by(1)Pt=E1−ν2εrt+ν⋅εθt,where ν stands for Poisson’s ratio and E stands for the elastic modulus.By Figure6 and equation (1), the pressures Pt were calculated, and the results are shown in Figure 7. For the sandstone specimen, the maximum value of blast-induced pressure was 62.9 MPa from Figure 7, which was much less than the pressure on the borehole wall [27].Figure 7
The pressure curve.


## 2.5. Propagation Behavior of the Cracks

CPG was pasted on the propagation paths of the cracks to further monitor the dynamic fracture of the cracks, and the oscilloscope was used to display and collect the voltage data of the CPG. Eleven groups of specimens were tested. The crack inclination angles changed from −75° to 75° for different groups to analyze the dynamic fracture of the cracks under blasting.As the crack propagated, the CPG wires were broken consecutively. The relationship of voltage signal with time was a step-shaped curve. The measurement results of sandstone specimens with crack inclination angleθ = −45° and 45° were plotted and are shown in Figures 8(a) and 8(b). The crack speed was calculated, and the results are shown in Figure 8(c). Meanwhile, the average crack propagation speeds are shown in Figure 8(d).Figure 8
The voltage signal recorded by CPG and the propagation speed of the cracks. Voltage signal curve of specimen with inclination angle (a) −45° and (b) 45°. Relation between (c) crack speed and crack length and (d) average crack speed versus crack length.
(a)
(b)
(c)
(d)From Figures8(a) and 8(b), one can find that, between the 5th and 7th wire and between the 7th and 9th wire, the crack speeds spent a relatively long time to propagate, and the propagation speed of crack inclination angle −45° specimen shown in Figure 8(c) was only 80.65 m/s which was much lower than the average speed 342 m/s shown in Figure 8(d), and similarly, for other specimens of crack inclination angle 45° shown in Figure 8(c), the propagation speed was 128.21 m/s which was less than the average speed 434 m/s, as shown in Figure 8(d). When θ equaled 45°, the average crack speed was the maximum, and when θ equaled −45°, the average crack speed was the minimum.


## 3. Numerical Simulation Analysis under Explosive Load

AUTODYN code was used in this paper to simulate crack fracture behavior under blasting which has been applied by more and more experts and scholars to solve crack dynamic fracture and proved to be feasible [16–26].


### 3.1. Numerical Models

In the process of simulation, the samples are divided into as many as possible tetrahedral elements. To make the calculation more accurate, around the crack, borehole and tunnel, the meshes are condensed, as shown in Figure9.Figure 9
Meshes of the specimen.Due to small deformation and pressure, in this study, linear EOS was used which can be expressed as the following equation:(2)Pρ=kρρ0−1,where k stands for the bulk modulus, ρ0 stands for the density of the initial state, and ρ stands for that of the current state.The JWL state equation was used to the explosive hexogen of the detonator which can be written as the following equation:(3)P=A⋅1−ωR1V⋅e−R1V+B⋅1−ωR2V⋅e−R2V+ωE0V,where V stands for the specific volume of detonation products, E0 stands for the total initial energy, P stands for the pressure, and A, B, R1, R2, and ω are fixed values, which are as follows: A=778.3GPa, B=7.071GPa, R1=4.2, R2=1.0, and ω=0.03.The failure criterions of maximum shear stress and principal stress were used in the rock specimens; i.e., for an element, if its principal stressσ1 reached its dynamic value σT of tensile strength, then it will be destroyed. Similarly, if its maximum τmax of shear stress reached its dynamic value τc of shear strength, then it will be also destroyed, that is to say:(4)σ1≤σT,orτmax≤τc.


### 3.2. Numerical Simulation Results

To study fracture rules of the cracks under blasting, the pressure curve obtained in the test was used in the AUTODYN3D code simulation. The simulation parameters are derived from Table 1. From the simulation results, the crack propagation paths could be obtained and were compared with experiment results, as shown in Figure 10.Figure 10
Comparison between (A) numerical results (AUTODYN simulation) and (B) experimental results (unit of crack length: mm). (a)θ = −75°, (b) θ = −60°, (c) θ = −45°, (d) θ = −30°, (e) θ = −15°, (f) θ = 0°, (g) θ = 15°, (h) θ = 30°, (i) θ = 45°, (j) θ = 60°, and (k) θ = 75°.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)From Figure10, for the crack propagation path, numerical simulation and experimental results were all not a straight line which may be due to the fact that dynamic energy release rate of the crack tip was not a fixed value. On the other hand, the experimental results were incompletely consistent with the numerical simulation, which may be due to the inhomogeneity of the experimental materials.


### 3.3. Crack Initiation Angle

Initiation angle of the crack was the angle between the crack initiation direction and the initial crack direction. To study the variation rule of the initiation angle, by PicPick software, the initiation angles of numerical simulation and test results were measured from Figure10, as shown in Table 3.Table 3
Crack initiation angle.
Crack inclination angle (°)
−75
−60
−45
−30
−15
0
15
30
45
60
75
—
Initiation angle (°)
Simulation value
63.28
57.57
51.02
39.49
34.28
21.27
33.40
56.37
68.35
85.12
97.21
—
Experimental value
63.68
58.81
51.14
39.93
34.49
21.76
33.88
56.69
68.81
85.66
97.70
—To understand the rule of crack initiation angle intuitively, the data in Table3 were plotted, as shown in Figure 11.Figure 11
The curves of crack initiation angle versus crack inclination angle from the blasting tests and the numerical simulation.From Figure11, the simulated and experimental values had similar change trends which all decreased first and then increased and were all minimum when θ equaled 0°, and the experimental value was slightly greater than the simulated value, which may be caused by the heterogeneity of experiment materials.


### 3.4. Crack Propagation Length

To study variation rule of crack propagation length, by PicPick software, numerical simulation and test results were measured from Figure11, as shown in Table 4.Table 4
Crack propagation length (L).
Crack inclination angle (°)
−75
−60
−45
−30
−15
0
15
30
45
60
75
—
L (mm)
Simulation value
43.54
52.63
63.75
77.35
86.78
92.69
100.77
111.56
136.44
108.2
96.82
—
Experimental value
43.15
52.26
62.46
75.89
84.46
90.83
98.68
109.32
134.85
106.88
94.24
—To understand the rule more intuitively, the data in Table4 were plotted, as shown in Figure 12.Figure 12
The curves of crack propagation length versus crack inclination angle from the numerical simulation and blasting tests.From Figure12, the simulated and experimental values had similar change trends which all increased first and then decreased and were all maximum when θ equaled 45°, and the experimental value was slightly smaller than the simulated value, which may be caused by the heterogeneity of experiment materials.


## 3.1. Numerical Models

In the process of simulation, the samples are divided into as many as possible tetrahedral elements. To make the calculation more accurate, around the crack, borehole and tunnel, the meshes are condensed, as shown in Figure9.Figure 9
Meshes of the specimen.Due to small deformation and pressure, in this study, linear EOS was used which can be expressed as the following equation:(2)Pρ=kρρ0−1,where k stands for the bulk modulus, ρ0 stands for the density of the initial state, and ρ stands for that of the current state.The JWL state equation was used to the explosive hexogen of the detonator which can be written as the following equation:(3)P=A⋅1−ωR1V⋅e−R1V+B⋅1−ωR2V⋅e−R2V+ωE0V,where V stands for the specific volume of detonation products, E0 stands for the total initial energy, P stands for the pressure, and A, B, R1, R2, and ω are fixed values, which are as follows: A=778.3GPa, B=7.071GPa, R1=4.2, R2=1.0, and ω=0.03.The failure criterions of maximum shear stress and principal stress were used in the rock specimens; i.e., for an element, if its principal stressσ1 reached its dynamic value σT of tensile strength, then it will be destroyed. Similarly, if its maximum τmax of shear stress reached its dynamic value τc of shear strength, then it will be also destroyed, that is to say:(4)σ1≤σT,orτmax≤τc.


## 3.2. Numerical Simulation Results

To study fracture rules of the cracks under blasting, the pressure curve obtained in the test was used in the AUTODYN3D code simulation. The simulation parameters are derived from Table 1. From the simulation results, the crack propagation paths could be obtained and were compared with experiment results, as shown in Figure 10.Figure 10
Comparison between (A) numerical results (AUTODYN simulation) and (B) experimental results (unit of crack length: mm). (a)θ = −75°, (b) θ = −60°, (c) θ = −45°, (d) θ = −30°, (e) θ = −15°, (f) θ = 0°, (g) θ = 15°, (h) θ = 30°, (i) θ = 45°, (j) θ = 60°, and (k) θ = 75°.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)From Figure10, for the crack propagation path, numerical simulation and experimental results were all not a straight line which may be due to the fact that dynamic energy release rate of the crack tip was not a fixed value. On the other hand, the experimental results were incompletely consistent with the numerical simulation, which may be due to the inhomogeneity of the experimental materials.


## 3.3. Crack Initiation Angle

Initiation angle of the crack was the angle between the crack initiation direction and the initial crack direction. To study the variation rule of the initiation angle, by PicPick software, the initiation angles of numerical simulation and test results were measured from Figure10, as shown in Table 3.Table 3
Crack initiation angle.
Crack inclination angle (°)
−75
−60
−45
−30
−15
0
15
30
45
60
75
—
Initiation angle (°)
Simulation value
63.28
57.57
51.02
39.49
34.28
21.27
33.40
56.37
68.35
85.12
97.21
—
Experimental value
63.68
58.81
51.14
39.93
34.49
21.76
33.88
56.69
68.81
85.66
97.70
—To understand the rule of crack initiation angle intuitively, the data in Table3 were plotted, as shown in Figure 11.Figure 11
The curves of crack initiation angle versus crack inclination angle from the blasting tests and the numerical simulation.From Figure11, the simulated and experimental values had similar change trends which all decreased first and then increased and were all minimum when θ equaled 0°, and the experimental value was slightly greater than the simulated value, which may be caused by the heterogeneity of experiment materials.


## 3.4. Crack Propagation Length

To study variation rule of crack propagation length, by PicPick software, numerical simulation and test results were measured from Figure11, as shown in Table 4.Table 4
Crack propagation length (L).
Crack inclination angle (°)
−75
−60
−45
−30
−15
0
15
30
45
60
75
—
L (mm)
Simulation value
43.54
52.63
63.75
77.35
86.78
92.69
100.77
111.56
136.44
108.2
96.82
—
Experimental value
43.15
52.26
62.46
75.89
84.46
90.83
98.68
109.32
134.85
106.88
94.24
—To understand the rule more intuitively, the data in Table4 were plotted, as shown in Figure 12.Figure 12
The curves of crack propagation length versus crack inclination angle from the numerical simulation and blasting tests.From Figure12, the simulated and experimental values had similar change trends which all increased first and then decreased and were all maximum when θ equaled 45°, and the experimental value was slightly smaller than the simulated value, which may be caused by the heterogeneity of experiment materials.


## 4. Dynamic SIFs (DSIFs) in the Process of Crack Propagation

DSIFs stand for propagation status of the crack. For rock materials, the fracture rules under impact loading have been widely researched by many experts and scholars [28–31]; however, the rules under blasting have been less studied.


### 4.1. Calculation of Static SIFs with the Help of ABAQUS Code

AUTODYN code has not been widely recognized in computing SIFs, but ABAQUS code has been widely recognized by many experts and scholars. So, in this study, for TMCSC specimens, finite element models of ABAQUS code under blasting were established. Triangular element CPS6 meshes were set near the crack tip, and quadrilateral elements CPS8 were set in other regions, as shown in Figure13.Figure 13
Mesh of a TMCSC specimen.In order to avoid the influence of crushed zone, when numerical simulation was carried out, the borehole radius was set to 30 mm, and the borehole wall was exerted by the loading from Figure7. In order to obtain accurate static SIFs under blasting, more meshes were added around the crack. Finally, SIFs kI0t and kII0t can be extracted from numerical simulation results.


### 4.2. Calculation of DSIFs

Because SIFs simulated by ABAQUS code was static under blasting, to get DSIFs, it needed to be corrected by the following equation:(5)kdt=kv⋅k0t,where kdt stands for DSIFs, v stands for crack propagation speed, k0t stands for static SIFs, and kv stands for universal function, as shown in the following equation:(6)kv=1−v/cR1−hv,where h=2/cdcs/cR21−cs/cd2, cd stands for P-wave speed, cR stands for Rayleigh wave speed, and cs stands for S-wave speed. If v=0, then kv=1, that stands for initiation state of the crack. If v=CR, then kv=0, that stands for propagation state of the crack with the Rayleigh wave speed, DSIF=0.The crack with an inclination angle whose value equaled −45° was used to illustrate the calculation process of DSIFs.To obtain SIFsk0t at the crack initiation moment, the crack initiation was simulated by the ABAQUS code, and SIFs k0t was obtained according to crack initiation moment measured by CPG, as shown in Figure 14. When the stress wave just reached the crack tip, SIFs kI0t was negative, and then, the absolute value of SIFs increased gradually, and SIFs kII0t was positive and then increased gradually. From the experiment, crack initiation moment was 72.4 μs, and thus, kI0 and kII0 were obtained whose absolute values were, respectively, 0.71 MPa m1/2 and 0.36 MPa m1/2.Figure 14
Calculation of DSIFs.When the crack tip extended 8 mm, the fifth wire of the CPG was broken, and the SIFsK0t curve at the crack tip could be calculated by the ABAQUS code, and the results are shown in Figure 15 (the solid line). At this time, crack propagation speed was 121.95 m/s, kv could be obtained by equation (6), and DSIFs Kdt could be obtained by equation (5). The curve of DSIFs Kdt versus time is shown in Figure 15 (dashed curve).Figure 15
DSIFs in the process of crack propagation.From this experiment, it could be seen that the fracture moment of the fifth wire was 99.4μs, and absolute values of DSIFs KIdt and KIIdt were obtained whose values were, respectively, 0.74 MPa·m1/2 and 0.15 MPa·m1/2.Similarly, according to the fracture moment measured by CPG, the other DSIFsKIdt and KIIdt were obtained, and the polynomial fitting curves of KIdt and KIIdt versus crack propagation speed were made by origin code which are shown in Figure 16. Their equation and the constants in the equation are shown in Table 5.Figure 16
Relationship between crack propagation speed and DSIFsKId and KIId.Table 5
Polynomial fitting equation of dynamic stress intensity factorKdv.
SIF
Fitting equation
b
0
b
1
b
2
b
3
b
4
b
5
b
6
b
7
K
I
d
v
K
I
d
v
=
∑
n
=
0
7
b
n
⋅
v
n
n
1.583
−0.005
1.094e−5
6.653e−8
2e−10
3e−13
2e−16
2e−20
K
II
d
v
K
II
d
v
=
∑
n
=
0
6
b
n
⋅
v
n
n
0.5792
0.0028
2.745e−5
1.083e−7
2e−10
2e−13
6e−17
—From Figure16, DSIFs KId and KIId decreased gradually, and there existed a critical propagation speed which made DSIFs convert from KId<KIId to KId>KIId. When the crack propagation speed was less than 285 m/s, KIId was greater than KId, which indicated that the crack propagation was mainly caused by shear action of stress wave. When the speed was greater than 285 m/s, KIId was less than KId, which indicated that the crack propagation was mainly caused by tensile action of stress wave.


### 4.3. Dynamic Energy Release Rate (G)

G is defined as the rate of mechanical energy flow out of the body and into the crack tip per unit crack advance. Freund derived the relation between G and DSIFs, for plane stress, which can be written by the following equation [32]:(7)G=1EdAIv⋅KI2+AIIv⋅KII2,where AIv and AIIv are the universal functions which are related to the crack propagation speed v, and when v=0, AIv=AIIv=1; when v≠0, AIv=v2⋅ad/1−νd⋅Cs2⋅D, AIIv=v2⋅as/1−νd⋅Cs2⋅D, D=4as⋅ad−1+as22, ad=1−v2/Cd2, as=1−v2/Cs2, Ed stands for dynamic elastic modulus, Cp stands for P-wave speed, Cs stands for S-wave speed, and νd stands for the dynamic Poisson ratio.The relation betweenG at the crack tip and crack propagation length L can be calculated by equation (7), as shown in Figure 17 (take the model whose crack inclination angle equals 45°as an example).Figure 17
The relation betweenG and L.From Figure17, it can be seen that G decreased first and then increased with the increase of L, and as L was about 21 mm for TMCSC, G was the minimum. G was not a constant, and it can be able to explain curvilinear propagation paths shown in Figure 10.


## 4.1. Calculation of Static SIFs with the Help of ABAQUS Code

AUTODYN code has not been widely recognized in computing SIFs, but ABAQUS code has been widely recognized by many experts and scholars. So, in this study, for TMCSC specimens, finite element models of ABAQUS code under blasting were established. Triangular element CPS6 meshes were set near the crack tip, and quadrilateral elements CPS8 were set in other regions, as shown in Figure13.Figure 13
Mesh of a TMCSC specimen.In order to avoid the influence of crushed zone, when numerical simulation was carried out, the borehole radius was set to 30 mm, and the borehole wall was exerted by the loading from Figure7. In order to obtain accurate static SIFs under blasting, more meshes were added around the crack. Finally, SIFs kI0t and kII0t can be extracted from numerical simulation results.


## 4.2. Calculation of DSIFs

Because SIFs simulated by ABAQUS code was static under blasting, to get DSIFs, it needed to be corrected by the following equation:(5)kdt=kv⋅k0t,where kdt stands for DSIFs, v stands for crack propagation speed, k0t stands for static SIFs, and kv stands for universal function, as shown in the following equation:(6)kv=1−v/cR1−hv,where h=2/cdcs/cR21−cs/cd2, cd stands for P-wave speed, cR stands for Rayleigh wave speed, and cs stands for S-wave speed. If v=0, then kv=1, that stands for initiation state of the crack. If v=CR, then kv=0, that stands for propagation state of the crack with the Rayleigh wave speed, DSIF=0.The crack with an inclination angle whose value equaled −45° was used to illustrate the calculation process of DSIFs.To obtain SIFsk0t at the crack initiation moment, the crack initiation was simulated by the ABAQUS code, and SIFs k0t was obtained according to crack initiation moment measured by CPG, as shown in Figure 14. When the stress wave just reached the crack tip, SIFs kI0t was negative, and then, the absolute value of SIFs increased gradually, and SIFs kII0t was positive and then increased gradually. From the experiment, crack initiation moment was 72.4 μs, and thus, kI0 and kII0 were obtained whose absolute values were, respectively, 0.71 MPa m1/2 and 0.36 MPa m1/2.Figure 14
Calculation of DSIFs.When the crack tip extended 8 mm, the fifth wire of the CPG was broken, and the SIFsK0t curve at the crack tip could be calculated by the ABAQUS code, and the results are shown in Figure 15 (the solid line). At this time, crack propagation speed was 121.95 m/s, kv could be obtained by equation (6), and DSIFs Kdt could be obtained by equation (5). The curve of DSIFs Kdt versus time is shown in Figure 15 (dashed curve).Figure 15
DSIFs in the process of crack propagation.From this experiment, it could be seen that the fracture moment of the fifth wire was 99.4μs, and absolute values of DSIFs KIdt and KIIdt were obtained whose values were, respectively, 0.74 MPa·m1/2 and 0.15 MPa·m1/2.Similarly, according to the fracture moment measured by CPG, the other DSIFsKIdt and KIIdt were obtained, and the polynomial fitting curves of KIdt and KIIdt versus crack propagation speed were made by origin code which are shown in Figure 16. Their equation and the constants in the equation are shown in Table 5.Figure 16
Relationship between crack propagation speed and DSIFsKId and KIId.Table 5
Polynomial fitting equation of dynamic stress intensity factorKdv.
SIF
Fitting equation
b
0
b
1
b
2
b
3
b
4
b
5
b
6
b
7
K
I
d
v
K
I
d
v
=
∑
n
=
0
7
b
n
⋅
v
n
n
1.583
−0.005
1.094e−5
6.653e−8
2e−10
3e−13
2e−16
2e−20
K
II
d
v
K
II
d
v
=
∑
n
=
0
6
b
n
⋅
v
n
n
0.5792
0.0028
2.745e−5
1.083e−7
2e−10
2e−13
6e−17
—From Figure16, DSIFs KId and KIId decreased gradually, and there existed a critical propagation speed which made DSIFs convert from KId<KIId to KId>KIId. When the crack propagation speed was less than 285 m/s, KIId was greater than KId, which indicated that the crack propagation was mainly caused by shear action of stress wave. When the speed was greater than 285 m/s, KIId was less than KId, which indicated that the crack propagation was mainly caused by tensile action of stress wave.


## 4.3. Dynamic Energy Release Rate (G)

G is defined as the rate of mechanical energy flow out of the body and into the crack tip per unit crack advance. Freund derived the relation between G and DSIFs, for plane stress, which can be written by the following equation [32]:(7)G=1EdAIv⋅KI2+AIIv⋅KII2,where AIv and AIIv are the universal functions which are related to the crack propagation speed v, and when v=0, AIv=AIIv=1; when v≠0, AIv=v2⋅ad/1−νd⋅Cs2⋅D, AIIv=v2⋅as/1−νd⋅Cs2⋅D, D=4as⋅ad−1+as22, ad=1−v2/Cd2, as=1−v2/Cs2, Ed stands for dynamic elastic modulus, Cp stands for P-wave speed, Cs stands for S-wave speed, and νd stands for the dynamic Poisson ratio.The relation betweenG at the crack tip and crack propagation length L can be calculated by equation (7), as shown in Figure 17 (take the model whose crack inclination angle equals 45°as an example).Figure 17
The relation betweenG and L.From Figure17, it can be seen that G decreased first and then increased with the increase of L, and as L was about 21 mm for TMCSC, G was the minimum. G was not a constant, and it can be able to explain curvilinear propagation paths shown in Figure 10.


## 5. Conclusion

(1)
From explosion experiments of physical models and AUTODYN numerical simulation, it could be seen that the crack propagation path was not a straight line which may be caused by the heterogeneous property of sandstone materials and the dynamic energy release rate may not be a constant. There were inflection points on the crack propagation paths which may be due to the fact that the stress wave encountered the hard grains in the rock and bypassed them, and at inflection points, there existed a minimum propagation speed which was called “crack arrest speed.” At the same time, it was also found that whenθ equaled 45°, average crack propagation speed was maximum. The experimental results were incompletely consistent with the numerical simulation, which may be due to the inhomogeneity of the experimental materials.
(2)
The simulated and experimental values of crack initiation angle had similar change trends which decreased first and then increased and were minimum whenθ equaled 0°. The experimental value was slightly greater than the simulated value, which may be caused by the heterogeneity of sandstone materials in the experiment.
(3)
When crack propagation speed increased, DSIFsKId and KIId decreased gradually. Besides, there existed a critical crack propagation speed which made DSIFs convert from KId<KIId to KId>KIId. When the crack propagation speed was less than 285 m/s, KIId was greater than KId, which indicated that the crack propagation in CPG monitoring range was mainly caused by shear action of stress wave. When the crack propagation speed was greater than 285 m/s, KIId was less than KId, which indicated that the crack propagation was mainly caused by tensile action of stress wave.
(4)
The dynamic energy release rateG decreased first and then increased with the increase of crack propagation length L, and as L was about 21 mm for the TMCSC specimen, G was the minimum. G was not a constant, and it can be able to explain curvilinear propagation paths.



---

*Source: 2308218-2019-05-20.xml*