# Experiment and Applications of Dynamic Constitutive Model of Tensile and Compression Damage of Sandstones

Authors: Chunliang Dong; Xiaoyu Lu; Guangming Zhao; Xiangrui Meng; Yingming Li; Xiang Cheng
Journal: Advances in Materials Science and Engineering
 (2021)

---

## Abstract

A dynamic constitutive model of tensile and compressive damage was constructed on the basis of the ZWT and statistical damage models, particularly by introducing the maximum tension and maximum shear stress criteria to solve the failure problem of the surrounding rock mass caused by deep excavation unloading. A shock compression and splitting test of sandstone specimens under different strain rates were performed by using a split Hopkinson pressure bar (SHPB). The constitutive model was developed again by LS-DYNA for the secondary numerical impact compression and split test of sandstones. Results demonstrated that the constructed dynamic constitutive model of tensile and compressive damage could considerably simulate tensile and compressive stress-strain relations and failure features of sandstones well. Lastly, the constitutive model was applied to conduct a numerical study on damage distribution and failure laws of the surrounding rocks at Gaochou Roadway, Luling Mine under cyclic excavation unloading. Results showed that the unloading failure of surrounding rocks has significant accumulation effects, and the accumulated damage on the floor is larger than those on the roof and roadway walls. The maximum breaking and damage depths are 0.4 m and 5.31 m, respectively. Circumferential damage showed an “umbrella-shaped” distribution pattern. With respect to trend, the damage accumulation effect at the rear part of the excavation face is stronger than that at the front part and the maximum influence distance is 6.4 m. However, the influencing degree of the accumulation effect attenuates gradually as advancing into the excavation face. The reliability of the numerical simulation is verified by combining the test results of the field geological radar on the roadway roof.


---

## Body

## 1. Introduction

With the increase of coal resource mining depth in China, the coal rock advancing and mining environment in deep areas are more complicated than those in shallow areas [1]. In a deep environment, such operations as roadway excavation and face mining destroy the balance of the original stress field, and rock units on the excavation face are under a stress state of a three-directional compressive one-way unloading [2, 3]. The result is in the frequent occurrences of the surrounding rock cracking and peeling. The research on dynamic unloading has attracted considerable attention [4–8]. The essence of disasters caused by deep rock unloading is that the tensile stress produced by the instantaneous unloading of crustal stress exceeds the tensile strength of rocks. Therefore, the dynamic tensile mechanical properties of rocks are vital in exploring the unloading problems of deep rock masses. A constitutive model that can reflect tensile and compressive mechanical properties of rocks is significant to the stability of underground rock mass engineering and the safety construction of workers.To date, among the semiempirical dynamic material model established based on relevant test data, constitutive theory and strength criteria are the most universal. Liu et al. [9] considered the distortional dilation and damage softening effects of rocks and proposed an elastic plastic damage constitutive model that is beneficial for the numerical calculation of embedding explosive stress waves. A numerical model based on bond-based Peridynamic theory was proposed by Ai et al. [10] to simulate crack propagation and dynamic mechanical properties of coal under the SHPB test. However, this model failed to consider the rate effect of rock materials. Zhai et al. [11] proposed a dynamic model of viscoelastoplastic damage, and its applicability was verified by the SHPB experiment of two rock materials. Based on the dynamic mechanical behavior results of frozen soil at different temperatures and high strain rates, a dynamic constitutive model of frozen soil under a multiaxial state was proposed by Zhu et al. [12]. The dynamic constitutive model of secondary loading surface of rock materials constructed by Zhou et al. [13] could reflect the hysteretic loop characteristics and rate effect under seismic loads. Li and Ma [14] observed the high-stress occurrence conditions of deep underground rock mass engineering and proposed an elastoplastic damage model applicable to complicated high-stress states based on segmented damage and the Hoek Brown strength criterion. Jiang et al. [15] constructed a dynamic constitutive model of sandstones with damage by using the viscoelastoplastic module model. Via discrete element method (DEM), Du et al. [16] investigated the mechanical behavior of granitic specimens under different hydrostatic confinements using SHPB testing, and the results showed that as the confinement increases, the failure mode changes from a mixed failure of surface splitting and inner X-type shear to a failure in a single shear band. Ou et al. [17] explored the dynamic mechanical properties of different dip angles of layered slate and constructed a dynamic damage constitutive model by considering macroscopic layers. Zheng et al. [18] constructed a dynamic constitutive model of coal and rock damage based on the DP failure criterion, which showed high consistency with test data within the strain rate range of 20−100 s−1. However, all of the preceding constitutive models were constructed under compression conditions, and none of them has considered applicability under tensile conditions. Zhu and Tang [19] constructed a dynamic constitutive model of rock damage based on the maximum tension stress criterion and Mohr-Coulomb criterion, with consideration to the tensile properties of rock materials. However, this model had difficulty in reflecting the prepeak mechanical properties of rocks because it was constructed based on the elastic constitutive theory. Tang et al. proposed the ZWT model [20], which is mainly applied to organic materials. The ZWT model can describe the prepeak nonlinear properties of materials and can also reflect the strain rate effect of materials. Hence, this model is widely used by geotechnical engineering scholars in recent years [21–24]. Nevertheless, only a few dynamic constitutive models of rocks can reflect the rate effect of materials under dynamic loading and also consider the tensile and compressive damage of rocks.Thus, uniaxial compression and split tests of sandstone specimens under different impact speeds were carried out by using SHPB. A dynamic constitutive model of tensile and compressive damage of sandstones was constructed based on the ZWT model and statistical damage theory. The secondary development of the model was performed using the Fortran language compiling constitutive subprogram in LS-DYNA. Numerical and laboratory test results were compared to verify the reasonability of the constructed constitutive model. Lastly, this constitutive model was applied to a deep roadway excavation project, and the reliability of the research results was verified by field measurement. This study provides a constitutive theoretical basis for studying the dynamic mechanical properties of rock masses and presents beneficial references to study the dynamic excavation unloading of deep rock engineering.


## 2. Impact Test of Rocks and Analysis of Results

### 2.1. Test Program and Principle

The rock core of the impact test was collected from the −693 m rock strata at the Gaochou Roadway of a coal mine in Huaibei. Rock samples were off-white and dull, with an elastic modulus of 15.6 GPa, Poisson's ratio of 0.22, and uniaxial compressive strength of 53 MPa. The rock core was cut using a ZS-100 vertical drilling machine and polished thereafter using an SHM-200 double-sided stone grinding machine. Processing unevenness of specimens was below 0.05 mm, and the end surface was perpendicular to the axis. The maximum error was below 0.3°. Cylinder standard specimens with a diameter of 50 mm and heights of 50 mm and 25 mm were prepared (Figure1(a)). Impact test was performed using the SHPB system in the National Key Laboratory of Mining Dynamics of Anhui University of Science and Technology (Figure 1(b)).Figure 1
Rock specimens and experiment systems. (a) Part of rock specimens (b) SHPB experiment system.
(a)(b)The SHPB test was based on two hypotheses of one-dimensional stress wave and stress homogeneity. The typical dynamic stress balance diagram of the SHPB test is shown in Figure2. As shown in Figure 2, the superposition curve of the incident and reflected waves in the platform stage of the reflected wave is consistent with the transmitted wave curve. This result indicates that a constant strain rate loading is maintained throughout the impact loading process, and the specimens are in a balanced stress state. Hence, data were processed using the two-wave method.Figure 2
Dynamic stress equilibrium diagram.At this moment, the mean stressσt, strain rate ε˙t, and strain εt of the specimens under a shock compression can be expressed as follows:(1)σt=−EAAsεtt,ε˙t=−2CLsεrt,εt=2CLs∫0tεrtdt,where εrt and εtt are the strains of the incident and transmitted waves, respectively; Ls and As are the height and cross-sectional area of the specimen; C, E, and A are the longitudinal wave velocity, elastic modulus, and cross-sectional area of the bars.The mean stressσt, strain rate ε˙t, and strain εt of the specimens under impact split conditions can be expressed as follows:(2)σt=−2EAπDsBεtt,ε˙t=2CDsεrt,εt=2CDs∫0τεrtdt,where Ds and B are the diameter and height, respectively, of the specimens.


### 2.2. Test Result Analysis

The stress-strain curves of the impact compression and split specimens under different strain rates, which were obtained from equations (1) and (2), are shown in Figure 3.Figure 3
Stress-strain relation curves. (a) Impact compression. (b) Impact split.
(a)(b)As shown in Figure3(a), the slope of the compressive stress-strain curve and compressive strength of sandstones increased with an increase in strain rate. With an increase in strain, the stress-strain curve develops a plastic platform and even a secondary peak, which is the plastic yielding feature. Analysis indicates that such a plastic deformation phenomenon of sandstones under a uniaxial impact compression may be related to existing fissures in specimens. Under impact loads, fissures are expanded into a fracture surface, and the material generally develops slippage and shearing along the fracture surface. Consequently, the plastic deformation phenomenon occurs.As shown in Figure3(b), the slope of the tensile stress-strain curve and tensile strength of sandstone increase with an increase in strain rate. However, sandstone has an evident prepeak plastic deformation under a low strain rate. With an increase in strain rate, plastic deformation disappears and the material mainly presents brittle failure characteristics.The elasticity modulus and dynamic strength of sandstone under impact compression and split tension increase with an increase in strain rate. The dynamic parameters are positively related to strain rate. Peak and maximum strain also increase with an increase in strain rate, and the area enclosed by the stress-strain curve also expands. That is, energies absorbed by rock damage failure increased continuously.


## 2.1. Test Program and Principle

The rock core of the impact test was collected from the −693 m rock strata at the Gaochou Roadway of a coal mine in Huaibei. Rock samples were off-white and dull, with an elastic modulus of 15.6 GPa, Poisson's ratio of 0.22, and uniaxial compressive strength of 53 MPa. The rock core was cut using a ZS-100 vertical drilling machine and polished thereafter using an SHM-200 double-sided stone grinding machine. Processing unevenness of specimens was below 0.05 mm, and the end surface was perpendicular to the axis. The maximum error was below 0.3°. Cylinder standard specimens with a diameter of 50 mm and heights of 50 mm and 25 mm were prepared (Figure1(a)). Impact test was performed using the SHPB system in the National Key Laboratory of Mining Dynamics of Anhui University of Science and Technology (Figure 1(b)).Figure 1
Rock specimens and experiment systems. (a) Part of rock specimens (b) SHPB experiment system.
(a)(b)The SHPB test was based on two hypotheses of one-dimensional stress wave and stress homogeneity. The typical dynamic stress balance diagram of the SHPB test is shown in Figure2. As shown in Figure 2, the superposition curve of the incident and reflected waves in the platform stage of the reflected wave is consistent with the transmitted wave curve. This result indicates that a constant strain rate loading is maintained throughout the impact loading process, and the specimens are in a balanced stress state. Hence, data were processed using the two-wave method.Figure 2
Dynamic stress equilibrium diagram.At this moment, the mean stressσt, strain rate ε˙t, and strain εt of the specimens under a shock compression can be expressed as follows:(1)σt=−EAAsεtt,ε˙t=−2CLsεrt,εt=2CLs∫0tεrtdt,where εrt and εtt are the strains of the incident and transmitted waves, respectively; Ls and As are the height and cross-sectional area of the specimen; C, E, and A are the longitudinal wave velocity, elastic modulus, and cross-sectional area of the bars.The mean stressσt, strain rate ε˙t, and strain εt of the specimens under impact split conditions can be expressed as follows:(2)σt=−2EAπDsBεtt,ε˙t=2CDsεrt,εt=2CDs∫0τεrtdt,where Ds and B are the diameter and height, respectively, of the specimens.


## 2.2. Test Result Analysis

The stress-strain curves of the impact compression and split specimens under different strain rates, which were obtained from equations (1) and (2), are shown in Figure 3.Figure 3
Stress-strain relation curves. (a) Impact compression. (b) Impact split.
(a)(b)As shown in Figure3(a), the slope of the compressive stress-strain curve and compressive strength of sandstones increased with an increase in strain rate. With an increase in strain, the stress-strain curve develops a plastic platform and even a secondary peak, which is the plastic yielding feature. Analysis indicates that such a plastic deformation phenomenon of sandstones under a uniaxial impact compression may be related to existing fissures in specimens. Under impact loads, fissures are expanded into a fracture surface, and the material generally develops slippage and shearing along the fracture surface. Consequently, the plastic deformation phenomenon occurs.As shown in Figure3(b), the slope of the tensile stress-strain curve and tensile strength of sandstone increase with an increase in strain rate. However, sandstone has an evident prepeak plastic deformation under a low strain rate. With an increase in strain rate, plastic deformation disappears and the material mainly presents brittle failure characteristics.The elasticity modulus and dynamic strength of sandstone under impact compression and split tension increase with an increase in strain rate. The dynamic parameters are positively related to strain rate. Peak and maximum strain also increase with an increase in strain rate, and the area enclosed by the stress-strain curve also expands. That is, energies absorbed by rock damage failure increased continuously.


## 3. Dynamic Constitutive Model of Rock Damage

### 3.1. Damage Model

Cracks and broken behaviors are inevitable during the impact process of rocks. This study applied a statistical damage model. The size of material infinitesimal can increase to cover microcracks and pores and can also be sufficiently small, thereby meeting the macroscopic homogeneity and isotropy of materials. The infinitesimal strength observes to the Weibull statistical distribution and its probability density function (PDF) is as follows:(3)pS=mFmFm−1exp−SFm,where S is infinitesimal strength and m and F are distribution parameters that can reflect the mechanical properties of materials.Damage variable is defined as follows:(4)D=NfN,where Nf is the number of broken infinitesimals and N is the total number of infinitesimals.The number of broken infinitesimals is as follows:(5)NfS=∫0SNpydy=N1−exp−SFm.The following statistical damage model can be obtained by combining equations (4) and (5):(6)D=1−exp−SFm.


### 3.2. Infinitesimal Strength

The maximum tensile stress and the maximum shear stress failure criteria were introduced and used as the infinitesimal strength of rocks. Therefore,(7)S=σ1,εv>0,σ1−σ3,εv<0,where σ1 and σ3 are the maximum and minimum principal stress, respectively.Hence, the following model of tensile and compressive damage was obtained:(8)D=1−exp−σ1Ftmt,εv>0,1−exp−σ1−σ3Fcmc,εv<0,where Ft, mt, Fc, and mc are the material parameters at the tensile and compressive damage of the rock materials.


### 3.3. Constitutive Model

The ZWT model is composed of one elastomerE0 and two Maxwell bodies in parallel connection. One Maxwell body E1,θ1 describes the material behavior under a low strain rate. Existing studies have proven that the time history was considerably smaller than θ1 in the high-speed impact test, which was substantially short for this Maxwell body to loose. Hence, it was disregarded in this study. The second Maxwell body E2,θ2 describes the material behavior under a high strain rate and conforms to the research scope of impact test in this study. Considering the parallel connection of the Maxwell body E2,θ2 and elastomer E0, the differential form of this constitutive model is as follows:(9)σ˙+σθ2=E0+E2ε˙+E0θ2ε,where E0 is the elasticity modulus of the elastic component and E2 and θ2 are the elasticity modulus and coefficient of viscosity, respectively, of the second Maxwell body.


## 3.1. Damage Model

Cracks and broken behaviors are inevitable during the impact process of rocks. This study applied a statistical damage model. The size of material infinitesimal can increase to cover microcracks and pores and can also be sufficiently small, thereby meeting the macroscopic homogeneity and isotropy of materials. The infinitesimal strength observes to the Weibull statistical distribution and its probability density function (PDF) is as follows:(3)pS=mFmFm−1exp−SFm,where S is infinitesimal strength and m and F are distribution parameters that can reflect the mechanical properties of materials.Damage variable is defined as follows:(4)D=NfN,where Nf is the number of broken infinitesimals and N is the total number of infinitesimals.The number of broken infinitesimals is as follows:(5)NfS=∫0SNpydy=N1−exp−SFm.The following statistical damage model can be obtained by combining equations (4) and (5):(6)D=1−exp−SFm.


## 3.2. Infinitesimal Strength

The maximum tensile stress and the maximum shear stress failure criteria were introduced and used as the infinitesimal strength of rocks. Therefore,(7)S=σ1,εv>0,σ1−σ3,εv<0,where σ1 and σ3 are the maximum and minimum principal stress, respectively.Hence, the following model of tensile and compressive damage was obtained:(8)D=1−exp−σ1Ftmt,εv>0,1−exp−σ1−σ3Fcmc,εv<0,where Ft, mt, Fc, and mc are the material parameters at the tensile and compressive damage of the rock materials.


## 3.3. Constitutive Model

The ZWT model is composed of one elastomerE0 and two Maxwell bodies in parallel connection. One Maxwell body E1,θ1 describes the material behavior under a low strain rate. Existing studies have proven that the time history was considerably smaller than θ1 in the high-speed impact test, which was substantially short for this Maxwell body to loose. Hence, it was disregarded in this study. The second Maxwell body E2,θ2 describes the material behavior under a high strain rate and conforms to the research scope of impact test in this study. Considering the parallel connection of the Maxwell body E2,θ2 and elastomer E0, the differential form of this constitutive model is as follows:(9)σ˙+σθ2=E0+E2ε˙+E0θ2ε,where E0 is the elasticity modulus of the elastic component and E2 and θ2 are the elasticity modulus and coefficient of viscosity, respectively, of the second Maxwell body.


## 4. Numerical Test

### 4.1. Numeralization of Constitutive Model

According to the method in [25], the differential constitutive of equation (9) can be expressed in the three-dimensional form as follows:(10)S˙ij+Sijθ2=2G0+G2e˙ij+2G0θ2eij,where Sij and eij are the stress and strain deviators, respectively; G0=E0/21+μ and G2=E2/21+μ are the shear modulus of the elastomer and second Maxwell body, and μ is the Poisson’s ratio.Given that the explicit dynamic algorithm of LS-DYNA is based on the incremented approach of the updated Lagrangian, the incremental form of the constitutive equation should be derived. The incremental form of equation (10) can be obtained as follows:(11)ΔSij+S¯ijθ2Δt=2G0+G2Δeij+2G0θ2e¯ijΔt,where ΔSij and Δeij are the stress and strain deviator increments; Δt is the time increment; and S¯ij and e¯ij are the mean deviatoric stress and mean deviatoric strain in one time incremental step.(12)S¯ij=12Sijt+1+Sijt,(13)e¯ij=12eijt+1+eijt,where Sijt are eijt are the stress and strain deviators at instant t; and Sijt+1 and eijt+1 are the stress and strain deviators at instant t+Δt.Combining equations (11)–(13), the incremental form of Sij can be obtained as follows:(14)ΔSij=4θ2G22θ2+Δt+2G0Δeij+2Δt2θ2+Δt2G0eijt−Sijt.It is considered to decomposing stress and strain tensors into spherical and deviatoric tensors as follows:(15)σij=Sij+δijσV,σV=13σkk,εij=eij+13δijεV,εV=εkk,where σij and εij are the stress and strain tensors; δij is the Kronecker symbol; and σV and εV are the spherical stress and strain tensors.Volume deformation is hypothesized linear elastic as follows:(16)σV=KεV,where K=E0+E2/31−2μ is the volume modulus of the material.Equations (15) and (16) were brought into Equation (14), resulting in the following equation:(17)Δσij=4θ2G22θ2+Δt+2G0Δεij+3K−2G0−2G23+3K−2G0Δt6θ22θ22θ2+ΔtδijΔεV+2θ22θ2+Δt2G0εijt+K−2G03δijεVt−σijt,where Δεij and ΔεV are the strain and spherical strain increments; εijt, σijt and εVt are the strain, spherical strain, and stress at instant t.Therefore, the nominal stress component att + 1 is as follows:(18)σijt+1=σijt+Δσij.According to the Lemaitre strain equivalent principle, effective stressesσ˜ij at t and t + 1 are as follows:(19)σ˜ijt=σijt1−Dt,(20)σ˜ijt+1=σijt+11−Dt+1.The difference between equations (19) and (20) was calculated, and the gained effective stress increment is as follows:(21)Δσ˜ij=Δσij−σijt+1Dt+1+σijtDt.Hence, effective stress can be expressed as follows:(22)σ˜ijt+1=σ˜ijt+Δσ˜ij.The preceding constitutive equation was compiled into a subprogram using Fortran language. The stress component and damage variable were updated at the end of each integral step. Meanwhile, the new stress component and damage value were transmitted to the main program of LS-DYNA in the next time step [26]. The calculation process of the constitutive subprogram is shown in Figure 4.Figure 4
Flowchart of the dynamic damage subprogram.


### 4.2. Finite Element Model

The finite element models of the numerical impact compression and split test are shown in Figure5. The relevant size was consistent with the impact test. A total of 64,020 elements were divided, including 11,519 specimen elements. The stress wave loading mode was used to simulate the impact loads of bullets, and the ends of a transmission bar were set as nonreflecting boundary conditions. When the damage of the rock element reached 1, the element was viewed as failed and deleted by ∗MAT_ADD_EROSION.Figure 5
Impact finite element model.The mechanical parameters of rocks in the numerical impact test are shown in Table1. In particular, E0 is the initial elasticity modulus, which is gained from the static compression test of rocks; mt and mc are material parameters that reflect the degree of rock strength concentration (i.e., 2.5 and 4.35, respectively); Ft and Fc are the macroscopic strengths of rock mass; and Ft, Fc, E2 and θ2 are adjusted according to the principle of consistency between the numerical simulation and test results.Table 1
Mechanical parameters of rocks in the numerical simulation.
ε˙ (s−1)E0 (GPa)Ft (106)mtE2 (GPa)Fc (106)mcθ2 (10−6s)Compression271.422.25.12.52.0504.3577.3334.222.26.72.54.0624.3577.3412.522.27.02.510.01054.3577.3501.622.29.32.514.01304.3577.3Split143.042.220.12.5203504.357.3262.142.218.12.5255004.357.3303.742.220.12.5285404.357.3340.942.225.12.5296104.357.3


### 4.3. Stress-strain Relation

Similar to the data processing in the impact test, the axial stresses of elements at 1 m of the incidence bar and 0.6 m of the transmission bar were extracted and brought into equations (1) and (2) to obtain the stress-strain curves of the numerical impact compression and splitting under different strain rates. The results are shown in Figure 6.Figure 6
Comparison between the numerical simulation and test results. (a) Impact compression. (b) Impact split.
(a)(b)As shown in Figure6, the numerical simulation results from the differential form of the constructed three-dimensional dynamic constitutive model are consistent with the impact test results. Given that underground rock masses are mainly under a three-dimensional compression state, with reference to the relevant research results [27], it shows plastic yielding before the stress-strain curve peak of sandstones under confining pressure, which is known as the brittle-ductile transformation characteristics. However, it shows brittle failure under uniaxial impacts. Therefore, plastic yield characteristics before the peak under uniaxial impact compression were disregarded in the constructed constitutive model. The impact tensile stress-strain curve of sandstones fits better with the testing curves. This result proves that the constructed constitutive model could reflect the rate effect and tensile and compressive mechanical properties of rock materials.


### 4.4. Failure Characteristics of Rocks

The comparison of specimen failures in the numerical test and impact test is shown in Figures7 and 8. Evidently, specimens mainly develop compressive failures under impact compression, forming debris and irregular rock pieces. Moreover, specimens produce additional debris and small rock pieces and delete numerous elements when the strain rate is high. Owing to compressive dilatation, the peripheral elements of specimens relatively develop tensile damage, but the damage value is not high.Figure 7
Rock failures under impact compression. (a)ε˙=271.4 (b) ε˙=334.2 (c) ε˙=412.5 (d) ε˙=501.6.
(a)(b)(c)(d)Figure 8
Rock failure under impact split. (a)ε˙=143.0 (b) ε˙=262.1 (c) ε˙=303.7 (d) ε˙=340.9.
(a)(b)(c)(d)Two ends of specimens under impact splitting present compressive failures. Along the vertical loading direction, specimens are broken into two large rock pieces radially. With the increase in strain rate, the quantity of the produced debris increases. However, the quantity and size of rock pieces are consistent.


## 4.1. Numeralization of Constitutive Model

According to the method in [25], the differential constitutive of equation (9) can be expressed in the three-dimensional form as follows:(10)S˙ij+Sijθ2=2G0+G2e˙ij+2G0θ2eij,where Sij and eij are the stress and strain deviators, respectively; G0=E0/21+μ and G2=E2/21+μ are the shear modulus of the elastomer and second Maxwell body, and μ is the Poisson’s ratio.Given that the explicit dynamic algorithm of LS-DYNA is based on the incremented approach of the updated Lagrangian, the incremental form of the constitutive equation should be derived. The incremental form of equation (10) can be obtained as follows:(11)ΔSij+S¯ijθ2Δt=2G0+G2Δeij+2G0θ2e¯ijΔt,where ΔSij and Δeij are the stress and strain deviator increments; Δt is the time increment; and S¯ij and e¯ij are the mean deviatoric stress and mean deviatoric strain in one time incremental step.(12)S¯ij=12Sijt+1+Sijt,(13)e¯ij=12eijt+1+eijt,where Sijt are eijt are the stress and strain deviators at instant t; and Sijt+1 and eijt+1 are the stress and strain deviators at instant t+Δt.Combining equations (11)–(13), the incremental form of Sij can be obtained as follows:(14)ΔSij=4θ2G22θ2+Δt+2G0Δeij+2Δt2θ2+Δt2G0eijt−Sijt.It is considered to decomposing stress and strain tensors into spherical and deviatoric tensors as follows:(15)σij=Sij+δijσV,σV=13σkk,εij=eij+13δijεV,εV=εkk,where σij and εij are the stress and strain tensors; δij is the Kronecker symbol; and σV and εV are the spherical stress and strain tensors.Volume deformation is hypothesized linear elastic as follows:(16)σV=KεV,where K=E0+E2/31−2μ is the volume modulus of the material.Equations (15) and (16) were brought into Equation (14), resulting in the following equation:(17)Δσij=4θ2G22θ2+Δt+2G0Δεij+3K−2G0−2G23+3K−2G0Δt6θ22θ22θ2+ΔtδijΔεV+2θ22θ2+Δt2G0εijt+K−2G03δijεVt−σijt,where Δεij and ΔεV are the strain and spherical strain increments; εijt, σijt and εVt are the strain, spherical strain, and stress at instant t.Therefore, the nominal stress component att + 1 is as follows:(18)σijt+1=σijt+Δσij.According to the Lemaitre strain equivalent principle, effective stressesσ˜ij at t and t + 1 are as follows:(19)σ˜ijt=σijt1−Dt,(20)σ˜ijt+1=σijt+11−Dt+1.The difference between equations (19) and (20) was calculated, and the gained effective stress increment is as follows:(21)Δσ˜ij=Δσij−σijt+1Dt+1+σijtDt.Hence, effective stress can be expressed as follows:(22)σ˜ijt+1=σ˜ijt+Δσ˜ij.The preceding constitutive equation was compiled into a subprogram using Fortran language. The stress component and damage variable were updated at the end of each integral step. Meanwhile, the new stress component and damage value were transmitted to the main program of LS-DYNA in the next time step [26]. The calculation process of the constitutive subprogram is shown in Figure 4.Figure 4
Flowchart of the dynamic damage subprogram.


## 4.2. Finite Element Model

The finite element models of the numerical impact compression and split test are shown in Figure5. The relevant size was consistent with the impact test. A total of 64,020 elements were divided, including 11,519 specimen elements. The stress wave loading mode was used to simulate the impact loads of bullets, and the ends of a transmission bar were set as nonreflecting boundary conditions. When the damage of the rock element reached 1, the element was viewed as failed and deleted by ∗MAT_ADD_EROSION.Figure 5
Impact finite element model.The mechanical parameters of rocks in the numerical impact test are shown in Table1. In particular, E0 is the initial elasticity modulus, which is gained from the static compression test of rocks; mt and mc are material parameters that reflect the degree of rock strength concentration (i.e., 2.5 and 4.35, respectively); Ft and Fc are the macroscopic strengths of rock mass; and Ft, Fc, E2 and θ2 are adjusted according to the principle of consistency between the numerical simulation and test results.Table 1
Mechanical parameters of rocks in the numerical simulation.
ε˙ (s−1)E0 (GPa)Ft (106)mtE2 (GPa)Fc (106)mcθ2 (10−6s)Compression271.422.25.12.52.0504.3577.3334.222.26.72.54.0624.3577.3412.522.27.02.510.01054.3577.3501.622.29.32.514.01304.3577.3Split143.042.220.12.5203504.357.3262.142.218.12.5255004.357.3303.742.220.12.5285404.357.3340.942.225.12.5296104.357.3


## 4.3. Stress-strain Relation

Similar to the data processing in the impact test, the axial stresses of elements at 1 m of the incidence bar and 0.6 m of the transmission bar were extracted and brought into equations (1) and (2) to obtain the stress-strain curves of the numerical impact compression and splitting under different strain rates. The results are shown in Figure 6.Figure 6
Comparison between the numerical simulation and test results. (a) Impact compression. (b) Impact split.
(a)(b)As shown in Figure6, the numerical simulation results from the differential form of the constructed three-dimensional dynamic constitutive model are consistent with the impact test results. Given that underground rock masses are mainly under a three-dimensional compression state, with reference to the relevant research results [27], it shows plastic yielding before the stress-strain curve peak of sandstones under confining pressure, which is known as the brittle-ductile transformation characteristics. However, it shows brittle failure under uniaxial impacts. Therefore, plastic yield characteristics before the peak under uniaxial impact compression were disregarded in the constructed constitutive model. The impact tensile stress-strain curve of sandstones fits better with the testing curves. This result proves that the constructed constitutive model could reflect the rate effect and tensile and compressive mechanical properties of rock materials.


## 4.4. Failure Characteristics of Rocks

The comparison of specimen failures in the numerical test and impact test is shown in Figures7 and 8. Evidently, specimens mainly develop compressive failures under impact compression, forming debris and irregular rock pieces. Moreover, specimens produce additional debris and small rock pieces and delete numerous elements when the strain rate is high. Owing to compressive dilatation, the peripheral elements of specimens relatively develop tensile damage, but the damage value is not high.Figure 7
Rock failures under impact compression. (a)ε˙=271.4 (b) ε˙=334.2 (c) ε˙=412.5 (d) ε˙=501.6.
(a)(b)(c)(d)Figure 8
Rock failure under impact split. (a)ε˙=143.0 (b) ε˙=262.1 (c) ε˙=303.7 (d) ε˙=340.9.
(a)(b)(c)(d)Two ends of specimens under impact splitting present compressive failures. Along the vertical loading direction, specimens are broken into two large rock pieces radially. With the increase in strain rate, the quantity of the produced debris increases. However, the quantity and size of rock pieces are consistent.


## 5. Engineering Applications

### 5.1. General Situation of the Project

The construction elevation of the Gaochou Roadway in Luling Coal III13 in Huaibei, Anhui Province, is from −715 m to −693 m. The roadway is next to the security coal column boundary of the 10# coal industrial square on the west and extends to the boundary of Mining Zone III and further to the mountains. There are the 10# coal mass in the II101 mining zone in the south and the III13# fully mechanized rock mining face in the north. High-angle and NNE-strike faults take the dominant role in the mining area, which are mainly compressive, torsional reverse faults. The construction section of the roadway is approximately 22 m (normal distance) away from the 8 and 9 floors and approximately 45 m away from the 10# roof. In the construction section, rock mass discloses the light gray argillaceous sandstone. The roadway profile and section are shown in Figure9.Figure 9
Profile and section of the Gaochou Roadway at the Luling coal mine.


### 5.2. Calculation Model and Parameters

Rock cores were collected from the wall of the excavation section of the III13 Gaochou Roadway. The laboratory three-dimensional geostress acoustic emission test method was used to calculate the geostress of rock masses; the vertical value wasσv=19.4MPa and the horizontal values were σH=23.7MPa and σZ=35.2MPa. The roadway section (4.6 m (width) × 3.5 m (height)) was a semicircular arch with straight walls. The mechanical properties of the rock mass are shown in Table 2.Table 2
Mechanical properties of the rock mass.
ρ/kg⋅m-3μE0(GPa)Ft (106)mtE2 (GPa)Fc (106)mcθ22547.50.2222.26.22.514.01024.3547.3The effect of faults was disregarded during modeling. For the convenience of applying excavation loads, the segmented millisecond blasting was simplified as the excavation mode of daily advancing by 2 m and full-section disposal excavation molding. It was advanced forward by 6 m according to the cyclic excavation of three segments (i.e., I, II, and III). A 5 m roadway was preserved before excavation. The model was divided into 541,000 elements and 562,176 nodes. The numerical calculation model is shown in Figure10.Figure 10
Numerical calculation model.The initial stress field (displacement field) was calculated using the implicit−explicit algorithm in LS-DYNA, and the dynamic excavation unloading-induced damage was calculated thereafter. Excavation loads were applied linearly. According to practical excavation situations in the Gaochou Roadway, the smooth blasting uses blast holes of 42 mm diameter, and the blast holes were set at an interval of 0.5 m. Blasting loads were calculated at 90 MPa according to the equivalent loading method [28]. The results in the following section are relative values because damage was disregarded during the calculation of geostress field. That is, the surrounding rocks are hypothesized to be free of initial damage before excavation.


### 5.3. Analysis of Numerical Results

#### 5.3.1. Circumferential Damage Laws of the Surrounding Rocks

During cyclic excavation, the maximum damage of surrounding rocks mainly concentrates in the middle section of the current excavation segment. In this section, the middle profiles (i.e., 6 m, 8 m, and 10 m profiles) of segments I, II, and III respectively, were chosen as representatives to analyze the circumferential evolutionary laws of damage in the surrounding rocks as in Figure11(a). Evidently, circumferential damage scope on the same profile of the surrounding rocks is expanded gradually as the excavation face advances forward, accompanied by the continuous growth of damage degree. The damage region presents an “umbrella-shaped” distribution circumferentially. Moreover, not all profiles surrounding rocks break during the excavation of the current segments. For example, failure of the 6 m profile occurs during the excavation of segment III. Breaking areas only concentrate on the floor. Moreover, the surrounding rock failures of the other two profiles first occur at the floor and transit to the roadway walls and roof thereafter. This situation indicates that the surrounding rock failures are closely related to the excavation of other segments, and excavation damage has an accumulation effect.Figure 11
Distributions of the surrounding rock damage. (a) Circumferential damage. (b) Striking damage.
(a)(b)The depth of circumferential damage is further analyzed in Figure12(a). As shown in Figure 12, the maximum damage depth of surrounding rocks on different profiles is at the floor, and damage depth on the wall is the second largest. However, the minimum damage depth is on the roof. Cyclic excavation can generally influence the damage of the 6 m and 8 m profiles. In particular, the maximum damage depth of the 6 m profile increases from 2.76 m to 5.31 m and the maximum damage depth of the 8 m profile increases from 0.39 m to 5.31 m. The maximum damage depths of the 6 m and 8 m profiles tend to be consistent during the excavation of segment III. The circumferential maximum damage depth of the 10 m profile is only 3.1 m, and the accumulation effect of the excavation damage is the smallest.Figure 12
Damage scope of the surrounding rocks. (a) Circumferential damage. (b) Striking damage.
(a)(b)


#### 5.3.2. Striking Damage Laws of Surrounding Rocks

The evolutionary laws of the striking damage of the surrounding rocks are shown in Figure11(b). In excavating segment II, the surrounding rocks have evident damage, most of which mainly concentrated in the surrounding rocks of segment II. The floor of segment II is broken on a large scale. In excavating segment III, the floor is broken on a large scale, while the walls and roofs are broken on a small scale. Consequently, the failure of segment II is intensified, and the maximum damage depth reaches 0.4 m. The maximum damage depth of segment I is the second highest, which is 0.1 m. Evidently, the excavation damage has an accumulation effect.To clearly master the variation laws of such accumulation effects, the current excavation face was used as a basis and the damage scopes in the rear and front regions of the surface under the cyclic excavation were calculated. The results are shown in Figure12(b). As the excavation face advances forward, the damage accumulation effect in the rear region of the excavation face was evident, and the striking damage scope rapidly increased from 1 m to 6.4 m. However, the accumulation effect in the front regions of the excavation face was stable, and the damage scope increased from 2.9 m to 3.2 m and decreased thereafter to 2.9 m. This result reflected that cyclic excavation influences the front region more than in the rear regions. That is, the damage scope and damage degree of the surrounding rocks in the front region of the excavation face have an evident accumulation effect, but the accumulation degree decreases gradually.Given that the preceding calculated results did not consider the initial geologic conditions of rock mass, the results should be slightly lower than the measured value.


#### 5.3.3. Field Measurement

This study used the LTD–2200 geological radar made by the China Radio Wave Research Institute. Measuring lines were set at the floor and middle walls of the roadway along the striking direction shown in Figure13.Figure 13
Layout of the measuring lines.Point measuring method was used. The channel space was set at 0.5 m. There is 15 m in the front regions of the excavation face and 6 m in the rear regions of the cyclic excavation. The total length of the measuring line is 21 m. The antenna emitted electromagnetic waves, and the frequency was set at 100 MHz. Electromagnetic waves will generate reflection and transmission on the differential surface of the dielectric constant in the surrounding rock masses [29]. Reflected echo signals were collected by the receiving antenna. Signal data were processed by edition, signal gain, and filtering. Space constraints have limited this study to analyze only the geological radar profiles of floors at the Gaochou Roadway as in Figure 14.Figure 14
Floor test results. (a) Before excavation. (b) First excavation. (c) Second excavation. (d) Third excavation.
(a)(b)(c)(d)As shown in Figure14, the receiving quality of radar reflected waves is generally high, and the signal-to-noise ratio is high as well. Longitudinally, the reflected wave of the electromagnetic wave is strong, and the waveform is explicit in the 150 ns of two-way travel time. However, the reflected wave is weak after 150 ns. Horizontally, the shallow radar in-phase axis continuity is good, and the horizontal reflecting layer can be recognized. As shown in Figure 14(a), three large fracture zones have been formed on the roadway floor, owning to the coupling effect of the original geological conditions and excavation dynamic loads. In addition, the maximum depth is nearly 8.5 m. After the first excavation, the reflected wave of rock mass in the near excavation segment is strengthened, and the event scope is expanded. The original fractured zone is nearly connected to the event of the new fractured zone in the excavation segment, and the new fracture depth is approximately 5.7 m as shown in Figure 14(b). After the second excavation, the reflected wave of the fractured zone in the excavation segment is strengthened, and the new fracture depth is approximately 5.1 m. The amplitude of the reflected wave of rock masses in the rear region of the excavation face is relatively increased as in Figure 14(c). After the third excavation, the reflected wave in the adjacent fractured zone in the rear region of the excavation face is strengthened significantly. Fracture depth increases from 5.1 m to 6 m. Rock mass in the excavation segment is broken significantly, and depth is approximately 7 m.According to the geological radar detection results of the floor at the Gaochou Roadway, there are accumulation effects of the fracture zone and fracture degree in rock mass at the floor under cyclic excavation. Such accumulation effects mainly come from the influences of the adjacent segment excavation, and they weaken significantly during interval excavation. These results are consistent with those of previous numerical simulations. Influenced by initial conditions, the depth and scope of the measured fractured zone are higher than those in the numerical simulation.


## 5.1. General Situation of the Project

The construction elevation of the Gaochou Roadway in Luling Coal III13 in Huaibei, Anhui Province, is from −715 m to −693 m. The roadway is next to the security coal column boundary of the 10# coal industrial square on the west and extends to the boundary of Mining Zone III and further to the mountains. There are the 10# coal mass in the II101 mining zone in the south and the III13# fully mechanized rock mining face in the north. High-angle and NNE-strike faults take the dominant role in the mining area, which are mainly compressive, torsional reverse faults. The construction section of the roadway is approximately 22 m (normal distance) away from the 8 and 9 floors and approximately 45 m away from the 10# roof. In the construction section, rock mass discloses the light gray argillaceous sandstone. The roadway profile and section are shown in Figure9.Figure 9
Profile and section of the Gaochou Roadway at the Luling coal mine.


## 5.2. Calculation Model and Parameters

Rock cores were collected from the wall of the excavation section of the III13 Gaochou Roadway. The laboratory three-dimensional geostress acoustic emission test method was used to calculate the geostress of rock masses; the vertical value wasσv=19.4MPa and the horizontal values were σH=23.7MPa and σZ=35.2MPa. The roadway section (4.6 m (width) × 3.5 m (height)) was a semicircular arch with straight walls. The mechanical properties of the rock mass are shown in Table 2.Table 2
Mechanical properties of the rock mass.
ρ/kg⋅m-3μE0(GPa)Ft (106)mtE2 (GPa)Fc (106)mcθ22547.50.2222.26.22.514.01024.3547.3The effect of faults was disregarded during modeling. For the convenience of applying excavation loads, the segmented millisecond blasting was simplified as the excavation mode of daily advancing by 2 m and full-section disposal excavation molding. It was advanced forward by 6 m according to the cyclic excavation of three segments (i.e., I, II, and III). A 5 m roadway was preserved before excavation. The model was divided into 541,000 elements and 562,176 nodes. The numerical calculation model is shown in Figure10.Figure 10
Numerical calculation model.The initial stress field (displacement field) was calculated using the implicit−explicit algorithm in LS-DYNA, and the dynamic excavation unloading-induced damage was calculated thereafter. Excavation loads were applied linearly. According to practical excavation situations in the Gaochou Roadway, the smooth blasting uses blast holes of 42 mm diameter, and the blast holes were set at an interval of 0.5 m. Blasting loads were calculated at 90 MPa according to the equivalent loading method [28]. The results in the following section are relative values because damage was disregarded during the calculation of geostress field. That is, the surrounding rocks are hypothesized to be free of initial damage before excavation.


## 5.3. Analysis of Numerical Results

### 5.3.1. Circumferential Damage Laws of the Surrounding Rocks

During cyclic excavation, the maximum damage of surrounding rocks mainly concentrates in the middle section of the current excavation segment. In this section, the middle profiles (i.e., 6 m, 8 m, and 10 m profiles) of segments I, II, and III respectively, were chosen as representatives to analyze the circumferential evolutionary laws of damage in the surrounding rocks as in Figure11(a). Evidently, circumferential damage scope on the same profile of the surrounding rocks is expanded gradually as the excavation face advances forward, accompanied by the continuous growth of damage degree. The damage region presents an “umbrella-shaped” distribution circumferentially. Moreover, not all profiles surrounding rocks break during the excavation of the current segments. For example, failure of the 6 m profile occurs during the excavation of segment III. Breaking areas only concentrate on the floor. Moreover, the surrounding rock failures of the other two profiles first occur at the floor and transit to the roadway walls and roof thereafter. This situation indicates that the surrounding rock failures are closely related to the excavation of other segments, and excavation damage has an accumulation effect.Figure 11
Distributions of the surrounding rock damage. (a) Circumferential damage. (b) Striking damage.
(a)(b)The depth of circumferential damage is further analyzed in Figure12(a). As shown in Figure 12, the maximum damage depth of surrounding rocks on different profiles is at the floor, and damage depth on the wall is the second largest. However, the minimum damage depth is on the roof. Cyclic excavation can generally influence the damage of the 6 m and 8 m profiles. In particular, the maximum damage depth of the 6 m profile increases from 2.76 m to 5.31 m and the maximum damage depth of the 8 m profile increases from 0.39 m to 5.31 m. The maximum damage depths of the 6 m and 8 m profiles tend to be consistent during the excavation of segment III. The circumferential maximum damage depth of the 10 m profile is only 3.1 m, and the accumulation effect of the excavation damage is the smallest.Figure 12
Damage scope of the surrounding rocks. (a) Circumferential damage. (b) Striking damage.
(a)(b)


### 5.3.2. Striking Damage Laws of Surrounding Rocks

The evolutionary laws of the striking damage of the surrounding rocks are shown in Figure11(b). In excavating segment II, the surrounding rocks have evident damage, most of which mainly concentrated in the surrounding rocks of segment II. The floor of segment II is broken on a large scale. In excavating segment III, the floor is broken on a large scale, while the walls and roofs are broken on a small scale. Consequently, the failure of segment II is intensified, and the maximum damage depth reaches 0.4 m. The maximum damage depth of segment I is the second highest, which is 0.1 m. Evidently, the excavation damage has an accumulation effect.To clearly master the variation laws of such accumulation effects, the current excavation face was used as a basis and the damage scopes in the rear and front regions of the surface under the cyclic excavation were calculated. The results are shown in Figure12(b). As the excavation face advances forward, the damage accumulation effect in the rear region of the excavation face was evident, and the striking damage scope rapidly increased from 1 m to 6.4 m. However, the accumulation effect in the front regions of the excavation face was stable, and the damage scope increased from 2.9 m to 3.2 m and decreased thereafter to 2.9 m. This result reflected that cyclic excavation influences the front region more than in the rear regions. That is, the damage scope and damage degree of the surrounding rocks in the front region of the excavation face have an evident accumulation effect, but the accumulation degree decreases gradually.Given that the preceding calculated results did not consider the initial geologic conditions of rock mass, the results should be slightly lower than the measured value.


### 5.3.3. Field Measurement

This study used the LTD–2200 geological radar made by the China Radio Wave Research Institute. Measuring lines were set at the floor and middle walls of the roadway along the striking direction shown in Figure13.Figure 13
Layout of the measuring lines.Point measuring method was used. The channel space was set at 0.5 m. There is 15 m in the front regions of the excavation face and 6 m in the rear regions of the cyclic excavation. The total length of the measuring line is 21 m. The antenna emitted electromagnetic waves, and the frequency was set at 100 MHz. Electromagnetic waves will generate reflection and transmission on the differential surface of the dielectric constant in the surrounding rock masses [29]. Reflected echo signals were collected by the receiving antenna. Signal data were processed by edition, signal gain, and filtering. Space constraints have limited this study to analyze only the geological radar profiles of floors at the Gaochou Roadway as in Figure 14.Figure 14
Floor test results. (a) Before excavation. (b) First excavation. (c) Second excavation. (d) Third excavation.
(a)(b)(c)(d)As shown in Figure14, the receiving quality of radar reflected waves is generally high, and the signal-to-noise ratio is high as well. Longitudinally, the reflected wave of the electromagnetic wave is strong, and the waveform is explicit in the 150 ns of two-way travel time. However, the reflected wave is weak after 150 ns. Horizontally, the shallow radar in-phase axis continuity is good, and the horizontal reflecting layer can be recognized. As shown in Figure 14(a), three large fracture zones have been formed on the roadway floor, owning to the coupling effect of the original geological conditions and excavation dynamic loads. In addition, the maximum depth is nearly 8.5 m. After the first excavation, the reflected wave of rock mass in the near excavation segment is strengthened, and the event scope is expanded. The original fractured zone is nearly connected to the event of the new fractured zone in the excavation segment, and the new fracture depth is approximately 5.7 m as shown in Figure 14(b). After the second excavation, the reflected wave of the fractured zone in the excavation segment is strengthened, and the new fracture depth is approximately 5.1 m. The amplitude of the reflected wave of rock masses in the rear region of the excavation face is relatively increased as in Figure 14(c). After the third excavation, the reflected wave in the adjacent fractured zone in the rear region of the excavation face is strengthened significantly. Fracture depth increases from 5.1 m to 6 m. Rock mass in the excavation segment is broken significantly, and depth is approximately 7 m.According to the geological radar detection results of the floor at the Gaochou Roadway, there are accumulation effects of the fracture zone and fracture degree in rock mass at the floor under cyclic excavation. Such accumulation effects mainly come from the influences of the adjacent segment excavation, and they weaken significantly during interval excavation. These results are consistent with those of previous numerical simulations. Influenced by initial conditions, the depth and scope of the measured fractured zone are higher than those in the numerical simulation.


## 5.3.1. Circumferential Damage Laws of the Surrounding Rocks

During cyclic excavation, the maximum damage of surrounding rocks mainly concentrates in the middle section of the current excavation segment. In this section, the middle profiles (i.e., 6 m, 8 m, and 10 m profiles) of segments I, II, and III respectively, were chosen as representatives to analyze the circumferential evolutionary laws of damage in the surrounding rocks as in Figure11(a). Evidently, circumferential damage scope on the same profile of the surrounding rocks is expanded gradually as the excavation face advances forward, accompanied by the continuous growth of damage degree. The damage region presents an “umbrella-shaped” distribution circumferentially. Moreover, not all profiles surrounding rocks break during the excavation of the current segments. For example, failure of the 6 m profile occurs during the excavation of segment III. Breaking areas only concentrate on the floor. Moreover, the surrounding rock failures of the other two profiles first occur at the floor and transit to the roadway walls and roof thereafter. This situation indicates that the surrounding rock failures are closely related to the excavation of other segments, and excavation damage has an accumulation effect.Figure 11
Distributions of the surrounding rock damage. (a) Circumferential damage. (b) Striking damage.
(a)(b)The depth of circumferential damage is further analyzed in Figure12(a). As shown in Figure 12, the maximum damage depth of surrounding rocks on different profiles is at the floor, and damage depth on the wall is the second largest. However, the minimum damage depth is on the roof. Cyclic excavation can generally influence the damage of the 6 m and 8 m profiles. In particular, the maximum damage depth of the 6 m profile increases from 2.76 m to 5.31 m and the maximum damage depth of the 8 m profile increases from 0.39 m to 5.31 m. The maximum damage depths of the 6 m and 8 m profiles tend to be consistent during the excavation of segment III. The circumferential maximum damage depth of the 10 m profile is only 3.1 m, and the accumulation effect of the excavation damage is the smallest.Figure 12
Damage scope of the surrounding rocks. (a) Circumferential damage. (b) Striking damage.
(a)(b)


## 5.3.2. Striking Damage Laws of Surrounding Rocks

The evolutionary laws of the striking damage of the surrounding rocks are shown in Figure11(b). In excavating segment II, the surrounding rocks have evident damage, most of which mainly concentrated in the surrounding rocks of segment II. The floor of segment II is broken on a large scale. In excavating segment III, the floor is broken on a large scale, while the walls and roofs are broken on a small scale. Consequently, the failure of segment II is intensified, and the maximum damage depth reaches 0.4 m. The maximum damage depth of segment I is the second highest, which is 0.1 m. Evidently, the excavation damage has an accumulation effect.To clearly master the variation laws of such accumulation effects, the current excavation face was used as a basis and the damage scopes in the rear and front regions of the surface under the cyclic excavation were calculated. The results are shown in Figure12(b). As the excavation face advances forward, the damage accumulation effect in the rear region of the excavation face was evident, and the striking damage scope rapidly increased from 1 m to 6.4 m. However, the accumulation effect in the front regions of the excavation face was stable, and the damage scope increased from 2.9 m to 3.2 m and decreased thereafter to 2.9 m. This result reflected that cyclic excavation influences the front region more than in the rear regions. That is, the damage scope and damage degree of the surrounding rocks in the front region of the excavation face have an evident accumulation effect, but the accumulation degree decreases gradually.Given that the preceding calculated results did not consider the initial geologic conditions of rock mass, the results should be slightly lower than the measured value.


## 5.3.3. Field Measurement

This study used the LTD–2200 geological radar made by the China Radio Wave Research Institute. Measuring lines were set at the floor and middle walls of the roadway along the striking direction shown in Figure13.Figure 13
Layout of the measuring lines.Point measuring method was used. The channel space was set at 0.5 m. There is 15 m in the front regions of the excavation face and 6 m in the rear regions of the cyclic excavation. The total length of the measuring line is 21 m. The antenna emitted electromagnetic waves, and the frequency was set at 100 MHz. Electromagnetic waves will generate reflection and transmission on the differential surface of the dielectric constant in the surrounding rock masses [29]. Reflected echo signals were collected by the receiving antenna. Signal data were processed by edition, signal gain, and filtering. Space constraints have limited this study to analyze only the geological radar profiles of floors at the Gaochou Roadway as in Figure 14.Figure 14
Floor test results. (a) Before excavation. (b) First excavation. (c) Second excavation. (d) Third excavation.
(a)(b)(c)(d)As shown in Figure14, the receiving quality of radar reflected waves is generally high, and the signal-to-noise ratio is high as well. Longitudinally, the reflected wave of the electromagnetic wave is strong, and the waveform is explicit in the 150 ns of two-way travel time. However, the reflected wave is weak after 150 ns. Horizontally, the shallow radar in-phase axis continuity is good, and the horizontal reflecting layer can be recognized. As shown in Figure 14(a), three large fracture zones have been formed on the roadway floor, owning to the coupling effect of the original geological conditions and excavation dynamic loads. In addition, the maximum depth is nearly 8.5 m. After the first excavation, the reflected wave of rock mass in the near excavation segment is strengthened, and the event scope is expanded. The original fractured zone is nearly connected to the event of the new fractured zone in the excavation segment, and the new fracture depth is approximately 5.7 m as shown in Figure 14(b). After the second excavation, the reflected wave of the fractured zone in the excavation segment is strengthened, and the new fracture depth is approximately 5.1 m. The amplitude of the reflected wave of rock masses in the rear region of the excavation face is relatively increased as in Figure 14(c). After the third excavation, the reflected wave in the adjacent fractured zone in the rear region of the excavation face is strengthened significantly. Fracture depth increases from 5.1 m to 6 m. Rock mass in the excavation segment is broken significantly, and depth is approximately 7 m.According to the geological radar detection results of the floor at the Gaochou Roadway, there are accumulation effects of the fracture zone and fracture degree in rock mass at the floor under cyclic excavation. Such accumulation effects mainly come from the influences of the adjacent segment excavation, and they weaken significantly during interval excavation. These results are consistent with those of previous numerical simulations. Influenced by initial conditions, the depth and scope of the measured fractured zone are higher than those in the numerical simulation.


## 6. Conclusions

(1)
The uniaxial impact compression and split experiments of sandstones under different strain rates demonstrate a rate correlation between the dynamic strength of sandstones and failure characteristics under impact loading conditions. Given an increase in strain rate, dynamic tensile and compressive strengths of the specimens increase, and the degree of crushing also increases accordingly. The stress-strain curves of sandstone under impact compression show evident prepeak plastic features but show only a low strain rate under impact split. Given an increase in strain rate, rocks show evident brittle mechanical properties.(2)
The dynamic constitutive model of the tensile and compressive damage of sandstones is constructed based on the ZWT viscoelastic model and statistical damage theory. According to the verification of numerical test, the proposed model can considerably reflect the mechanical behaviors and failure features of sandstones under different strain rates, thereby verifying the applicability of the proposed model.(3)
A case study based on the segmented excavation projects at the Gaochou Roadway at the Luling coal mine is carried out. Damage distribution and failure laws of deep surrounding rocks under cyclic excavation are analyzed through numerical simulation. The results indicate that surrounding rock damage has significant accumulation effects. Moreover, the floor experiences the most serious damage, followed by the walls and roof. The accumulation effect influences rocks in the rear regions of the excavation face compared with that in the front regions. As the excavation face advances, the damage accumulation of surrounding rocks in the rear regions attenuates gradually.(4)
According to a comparison between the numerical simulation and field measurement results, the measured fracture depth at the floor of Gaochou Roadway in the excavation segment also has an accumulation effect. Influenced by geological conditions, the measured results are relatively high. This primary reason is that the existing initial conditions of the roadway are disregarded in the numerical analysis.



---

*Source: 2492742-2021-05-25.xml*