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.
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 (Figure 1(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)).
(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 Figure 2. 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.
At this moment, the mean stress , strain rate , and strain of the specimens under a shock compression can be expressed as follows:where and are the strains of the incident and transmitted waves, respectively; and are the height and cross-sectional area of the specimen; , E, and A are the longitudinal wave velocity, elastic modulus, and cross-sectional area of the bars.
The mean stress , strain rate , and strain of the specimens under impact split conditions can be expressed as follows: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.
(a)
(b)
As shown in Figure 3(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 Figure 3(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: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:where Nf is the number of broken infinitesimals and N is the total number of infinitesimals.
The number of broken infinitesimals is as follows:
The following statistical damage model can be obtained by combining equations (4) and (5):
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,where and are the maximum and minimum principal stress, respectively.
Hence, the following model of tensile and compressive damage was obtained:where , , , and 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 elastomer and two Maxwell bodies in parallel connection. One Maxwell body describes the material behavior under a low strain rate. Existing studies have proven that the time history was considerably smaller than 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 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 and elastomer , the differential form of this constitutive model is as follows:where E0 is the elasticity modulus of the elastic component and and 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:where and are the stress and strain deviators, respectively; and 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:where and are the stress and strain deviator increments; is the time increment; and and are the mean deviatoric stress and mean deviatoric strain in one time incremental step.where are are the stress and strain deviators at instant t; and and are the stress and strain deviators at instant .
Combining equations (11)–(13), the incremental form of can be obtained as follows:
It is considered to decomposing stress and strain tensors into spherical and deviatoric tensors as follows:where and are the stress and strain tensors; is the Kronecker symbol; and and are the spherical stress and strain tensors.
Volume deformation is hypothesized linear elastic as follows:where is the volume modulus of the material.
Equations (15) and (16) were brought into Equation (14), resulting in the following equation:where and are the strain and spherical strain increments; , and are the strain, spherical strain, and stress at instant t.
Therefore, the nominal stress component at t + 1 is as follows:
According to the Lemaitre strain equivalent principle, effective stresses at t and t + 1 are as follows:
The difference between equations (19) and (20) was calculated, and the gained effective stress increment is as follows:
Hence, effective stress can be expressed as follows:
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.
4.2. Finite Element Model
The finite element models of the numerical impact compression and split test are shown in Figure 5. 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.
The mechanical parameters of rocks in the numerical impact test are shown in Table 1. In particular, is the initial elasticity modulus, which is gained from the static compression test of rocks; and are material parameters that reflect the degree of rock strength concentration (i.e., 2.5 and 4.35, respectively); and are the macroscopic strengths of rock mass; and , , and are adjusted according to the principle of consistency between the numerical simulation and test results.
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.
(a)
(b)
As shown in Figure 6, 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 Figures 7 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.
(a)
(b)
(c)
(d)
(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 Figure 9.
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 and the horizontal values were and . 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.
The 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 Figure 10.
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 Figure 11(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.
(a)
(b)
The depth of circumferential damage is further analyzed in Figure 12(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.
(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 Figure 11(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 Figure 12(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 Figure 13.
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.
(a)
(b)
(c)
(d)
As shown in Figure 14, 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.Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (Grant nos. 52074006, 51404011, 51974009, 51774012, 51874002, and 52004005) and Anhui University of Science and Technology Introduction of Talents Research Fund.