A Strain-Softening Constitutive Model of Heterogeneous Rock Mass Considering Statistical Damage and Its Application in Numerical Modeling of Deep Roadways

During the construction of underground caverns, the stability of deep underground cavern excavation, which affects the safety and sustainable development of such projects, is a hot issue. First, based on the mechanical properties of surrounding rock in deep tunnels, the strain-softening behavior, damage, and heterogeneity of rock masses are analyzed. Then, a strain-softening model of heterogeneous jointed rock mass that considers statistical damage (SSD) is developed and implemented through FLAC3D simulation software. Finally, the SSD is applied to a deep roadway in the Jinchuan mining area, and a comparative analysis of the computation results of the Mohr–Coulomb (MC) model and the strain-softening (SS) model are carried out. The numerical results are compared with the field-monitoring results, which show that the SSD model simulated the behavior of the surrounding rocks well. The results show that the deformations of the roof and floor are larger, which may serve as a reference for the support pattern of deep roadways.


Introduction
Rock mass is a kind of complex geological body composed of structural planes and rock blocks.It is encountered constantly in water conservation and hydropower projects, tunnel engineering, and mining engineering.With the construction of a series of large-scale projects around the world, many countries have carried out experimental research on the rock mechanics, thermodynamics, seepage, and long-term stability of rock mass in deep underground engineering.In the process of underground engineering construction, a deep understanding of the physical and mechanical properties of the rock mass is required, especially as the buried depth increases.Because of the high ground stress, temperature and osmotic pressure in the deep rock mass, the rock burst risk will be further increased, the time effect of rock mass deformation will be further revealed, and the effect of stress redistribution caused by excavation on the stability of deep underground engineering is more prominent [1][2][3].
Researchers from various countries have carried out systematic research work on deep underground projects, such as the South-to-North Water Transfer Project and the Jinping II Hydropower Station in China [4,5], in situ tests at Yucca Mountain in the United States [6] and in the Aspö Hard Rock Laboratory (HRL) in Sweden, and the DECOVALEX project, which is jointly carried out by European countries [7,8].These large projects promote the development of rock mass mechanics.Current research methods mainly involve laboratory tests, on-site monitoring, numerical simulations, and theoretical analysis.As the main research method, laboratory tests have many disadvantages, such as difficult sample preparation, poor controllability of loading conditions, and poor visualization in the test process.So, it is difficult to obtain realistic results with this method.
As a supplement to the physical tests, numerical tests can overcome most of the difficulties mentioned above.J.C. Simo found that the softening law must be re-interpreted in a distributional sense for continuum solutions, and it provides a precise physical interpretation of the softening modulus and a new class of anisotropic rate-independent damage models for brittle materials [9].M. Nikolic proposed a model which is suitable for assessing the failure of heterogeneous rock with cracks or defects based on the discontinuity approach, and he presented two models (I and II) to represent the intact part and the microcracks, respectively [10,11].C. W. Boon put forward a systematic approach to provide useful insight for tunnel support in jointed rock masses based on the distinct element method (DEM) and analyzed the effects of support structures with different parameters [12].O.K. Mahabadi presented a new numerical code based on the combined finite-discrete element method (FDEM), named Y-Geo, and illustrated two of its geomechanical applications [13].Wei Wu put forward a discontinuous deformation analysis (DDA)-key block theory that can alleviate some shortcomings of the key block theory and illustrated an example [14].Leandro L. improved the classic lattice spring model by considering the Synthetic Rock Mass technique and the Barton-Bandis joint constitutive model to analyze the stability of rock tunnels under low in situ stress conditions [15].Fei Zheng developed a Newmark-based predictor corrector solution (NPC) approach in order to improve the performance of the DDA solution module in modeling discontinuous problems [16].Moreover, in comparison with the elastoplastic model, the stain-softening model can describe the post-peak mechanical behavior of rocks well.In order to describe the strain-softening characteristics more accurately, researchers have investigated them using theoretical analyses and laboratory tests, such as the strain-softening model considering the quantitative relationship between the plastic-softening parameters and the confining pressure of Cui, and the strain-softening model considering the variation in the dilation angle of Alejano [17,18].However, few models have considered the following three conditions simultaneously: (1) the reduction in the elastic constant caused by rock failure; (2) rock failure caused by damage accumulation; and (3) the heterogeneity of rocks and their application in practical engineering.
FLAC3D is a type of finite-difference software that cannot simulate micro-cracks, so the strain-softening constitutive is used to investigate this mechanical behavior of rock mass.However, strain-softening constituted according to the FLAC3D only takes the reduction of strength parameters, such as c and ϕ, into account but not the deformation parameters, such as the elastic modulus and Poisson s ratio.Therefore, the research in this paper aims at a secondary development of FLAC3D on the basis of strain-softening.In the improved constitutive model, statistical damage (SSD) is considered and the heterogeneity of engineering rock mass is also taken into account.Then, a numerical simulation calculation is carried out on a test roadway in the Jinchuan nickel mining area in China, and the obtained results are compared with the field-monitoring results.Finally, the rationality of the proposed constitutive model is discussed, and the case study provides a reference for the simulation calculation of other similar projects.

Strain-Softening Behavior of Rocks
The research shows that the strength parameters of surrounding rock will be obviously reduced in the nonlinear stage during the excavation process of deep-buried projects due to the connections among the macro-cracks in the so-called strain-softening rock [19,20].When the stress of the rock reaches peak strength, the deterioration will become more apparent as it continues to increase.The material parameters in the commonly-used Mohr-Coulomb (MC) model are preset to be constant, which cannot reflect the softening effect of the surrounding rock caused by the variation in the surrounding rock parameters.Thus, many scholars are aiming to establish a constitutive relationship that can reflect the nonlinear mechanical response of rock based on a numerical simulation and the strain-softening theory [21,22].
The strain-softening (SS) model considers the cohesion, friction angle, and dilatancy angle of the material, which change with the plastic strain, and we can define the relationship through a piecewise linear function.In addition, the yield criterion, potential function, plastic flow rule and stress correction in the SS model are exactly the same as in the MC model.It can be seen from Figure 1 that the strain of the element is only elastic strain (ε e ) in the elastic phase, that is, The strain of the element consists of elastic strain (ε e ) and plastic strain (ε p ) after the element yields, that is, The yield function of the element is, where f s = 0 is the shear failure; f t = 0 is the tensile failure; σ 1 , σ 3 are the maximum and minimum principal stresses; σ t is the tensile strength; C is the cohesion; ϕ is the internal friction angle; and surrounding rock parameters.Thus, many scholars are aiming to establish a constitutive relationship that can reflect the nonlinear mechanical response of rock based on a numerical simulation and the strain-softening theory [21,22].The strain-softening (SS) model considers the cohesion, friction angle, and dilatancy angle of the material, which change with the plastic strain, and we can define the relationship through a piecewise linear function.In addition, the yield criterion, potential function, plastic flow rule and stress correction in the SS model are exactly the same as in the MC model.It can be seen from Figure 1 that the strain of the element is only elastic strain (εe) in the elastic phase, that is, The strain of the element consists of elastic strain (εe) and plastic strain (εp) after the element yields, that is, The yield function of the element is, where 0 s f  is the shear failure; 0 t f  is the tensile failure; (1 sin ) / (1 sin ) Stress-strain relation of an element.
In the numerical calculation process, after the unit yields, the new shear strength parameters of the unit are calculated according to the plastic strain at each iteration step.The parameters are updated by the non-linear equation between the shear strength parameter and the plastic strain and then used in the next iteration step.Through this cycle, the strain-softening behavior of the rocks can be reflected.

Statistical Damage Constitutive Model
For deep-buried tunnels, due to the high ground stress, the stress redistribution of surrounding rock caused by excavation leads to damage in the surrounding rock of the tunnel, and when the damage accumulates to a certain degree, the surrounding rock fails.A variety of micro-defects are distributed inside the rock in its natural state, and following engineering disturbances, micro-defects will occur and develop under load.Characterizing the micro-defects affecting the physical and mechanical properties of rock and then applying them to the calculation or evaluation of rock In the numerical calculation process, after the unit yields, the new shear strength parameters of the unit are calculated according to the plastic strain at each iteration step.The parameters are updated by the non-linear equation between the shear strength parameter and the plastic strain and then used in the next iteration step.Through this cycle, the strain-softening behavior of the rocks can be reflected.

Statistical Damage Constitutive Model
For deep-buried tunnels, due to the high ground stress, the stress redistribution of surrounding rock caused by excavation leads to damage in the surrounding rock of the tunnel, and when the damage accumulates to a certain degree, the surrounding rock fails.A variety of micro-defects are distributed inside the rock in its natural state, and following engineering disturbances, micro-defects will occur and develop under load.Characterizing the micro-defects affecting the physical and mechanical properties of rock and then applying them to the calculation or evaluation of rock mechanical properties has always been a research hotspot of the rock mechanics field.The statistical damage model provides a feasible way to address this problem [23].
This method replaces the measurement of micro-defects by the assessment of the degree of damage.It assumes that these defects obey a distribution function and then establishes a statistical damage model of rock.The establishment of the rock statistical damage model is usually divided into two steps.The first step is to select a strength criterion of rock such as Mohr-Coulomb, Drucker-Prager, or the strain-softening measure used in this paper.The second step involves fitting an analytical probability distribution to characterize the stress level of rock elements.Several empirical distributions can be used, such as the function distribution and the normal distribution, with the Weibull distribution being the most popular method [24][25][26].
Regarding shear failure, this paper refers to Yang' assumption that the cohesion of rock micro-units obeys the three-parameter Weibull distribution, and the rock statistical damage model is expressed as [27]: where c is the cohesion of rock micro-units; c 0 is the initial cohesion; λ is the scale parameter whose value is in the same order of magnitude as c; m is the shape parameter; In the calculation of the ith incremental step, the maximum principal stress σ 1 and the minimum principal stress σ 3 are substituted into the damage factor D i−1 of the previous step, and c i is calculated by Equation (5) in the current step.The initial state of the rock is not damaged, so the D 0 is 0 in the first incremental step.If c i ≤ c 0 , there is no damage in the model, so D i is 0; If c i > c 0 , D i is calculated by Equation (4).In the calculation of step i plus the first incremental step, c i+1 can be calculated by substituting the damage factor D i into Equation (5), and then c i+1 is substituted into Equation (4) to calculate D i+1 , which is applied in the next incremental step.In every incremental step, the damage factor can be determined based on this loop computation.
Regarding the tensile failure, referring to the damage accumulation theory, it is considered that damage to the micro-units is gradually accumulated with the tensile failure.When the rock undergoes tensile yielding, the bearing capacity is not completely lost immediately.D 0 is defined as the initial damage, and then the calculation of the cumulative damage D is performed on the micro-units [28].
Under tension and shearing, there will be different degrees of damage inside the rock mass.According to the strain equivalence principle, the elastic parameters decrease at the stage of damage, following [29,30]: where E 0 is the initial elastic modulus of a micro-body.

Heterogeneity in Rocks
Rock is a kind of natural heterogeneous material.In general, when we study the mechanical properties and failure process of rocks, rocks are regarded as particles that should be analyzed within the framework of continuum mechanics.However, rocks are large enough to contain micro-cracks and their spatial distributions of physical and mechanical properties are discontinuous.
Heterogeneity has a great effect on the mechanical and deformational properties of rocks.Weibull described the heterogeneity of materials with probability and statistical methods.In this paper, the Weibull statistical distribution, which is widely used, is adopted to describe the heterogeneity of rock through the following equation [31]: where x is the attribute parameter of the element.This paper takes the elastic modulus as a representative to study the heterogeneity of rock.k is a scale parameter that reflects the average value of the elastic modulus.m is a shape parameter that reflects the homogeneity of the elastic modulus.
The heterogeneous assignment of the elastic modulus based on the Fish function can be accomplished by the following three steps: (1) Generate a (0, 1) uniformly distributed random number sequence F through the urand function in FLAC 3D ; (2) calculate the value of the elastic modulus by the formula x = λ(− ln(F)) 1/m ; (3) Assign values to individual units by command order.

Implementation of the SSD Model in FLAC 3D
In recent years, with the development of computer technology, more and more researchers have adopted numerical methods in their studies because of the convenience, repeatability, and economy.The numerical methods include the finite element method, the discrete element method, and the finite difference method.RFPA (Discontinuous Deformation Analysis), R-T 2D (Rock and Tool Interaction), and FLAC 3D are commonly used in geotechnical engineering, and FLAC 3D is used in this paper because it provides users with a flexible interface and has advantages in simulating rock failure processes and large deformations [32,33].
FLAC 3D uses the object-oriented language C++ for programming, and it provides users with a number of model examples, including header files (.h) and source files (.cpp).If there is no one appropriate model for the research problem, users can customize the constitutive model and compile it into a dynamic link library file (.dll) for software calls.Furthermore, the customized constitutive model and the software-owned model have the same execution efficiency.The main function of the constitutive model is to provide a strain increment and obtain a new stress value.It also has some auxiliary functions including the model name, material properties, and interactions between the model and code [34].
Based on the above analysis, a strain-softening constitutive model of heterogeneous jointed rock mass that considers statistical damage was developed, and the process of numerical implementation using software is shown in Figure 2.

Uniaxial Compression Test Based on SSD Model
In order to prove the reasonability of the model, a uniaxial compression test was carried out based on the MC, SS, and SSD models.The numerical model, as shown in Figure 3, was established according to the actual standard size of rock specimens which is Φ 5 × 10 cm, and the parameters are set based on the values for marble shown in Table 1.The mesh has an inevitable impact on the results of the calculations-the finer the mesh is, the closer the result is to the actual value-but the trend of the calculation results is consistent.The meshes are divided based on the calculation efficiency and the accuracy of the result in this paper.The stain-stress-acoustic emission ringing counts curves are shown in Figure 4.

Uniaxial Compression Test Based on SSD Model
In order to prove the reasonability of the model, a uniaxial compression test was carried out based on the MC, SS, and SSD models.The numerical model, as shown in Figure 3, was established according to the actual standard size of rock specimens which is Φ 5 × 10 cm, and the parameters are set based on the values for marble shown in Table 1.The mesh has an inevitable impact on the results of the calculations-the finer the mesh is, the closer the result is to the actual value-but the trend of the calculation results is consistent.The meshes are divided based on the calculation efficiency and the accuracy of the result in this paper.The stain-stress-acoustic emission ringing counts curves are shown in Figure 4.    Figure 4a shows an ideal elastoplastic curve under the MC model-the stress increases as the strain increases until the yield strength is reached.Under the SS model shown in Figure 4b, the stress can be divided into three stages: monotonic increasing before the critical strain, decline, and basic stability.The yield stresses under the MC model and the SS model are almost the same-about 45 MPa.It can be seen from Figure 4c that the SSD model proposed in this paper can reflect the full process of failure for rock, especially the strain-softening property.The stress-strain curve is basically the same as the previous two in the elastic stage, except that the yield stress is slightly smaller.However, after the rock yield, a cliff fall appears because the statistical damage is considered.The damage on the weak part of the rock accumulates and the mechanical parameters are continuously reduced until rock failure, which explains why although the intact rock is very strong, the rock mass is destroyed along the weak structural plane during practical engineering.

Engineering Application
To validate the rationality of the SSD model proposed above, a comparative numerical calculation was carried out to analyze the deformation of surrounding rock at a test roadway in the Jinchuan II mining area based on the MC model, SS model, and SSD model, respectively.

Geological Background
Jinchuan Mine, the largest nickel production base in China, is located in Jinchang City, Gansu Province.The ore deposit is about 6.5 km long, 200 m wide, and over 1000 m deep.It is a large and complex deposit with high ground stress, deep buried ore bodies, and broken rocks.The engineering    Figure 4a shows an ideal elastoplastic curve under the MC model-the stress increases as the strain increases until the yield strength is reached.Under the SS model shown in Figure 4b, the stress can be divided into three stages: monotonic increasing before the critical strain, decline, and basic stability.The yield stresses under the MC model and the SS model are almost the same-about 45 MPa.It can be seen from Figure 4c that the SSD model proposed in this paper can reflect the full process of failure for rock, especially the strain-softening property.The stress-strain curve is basically the same as the previous two in the elastic stage, except that the yield stress is slightly smaller.However, after the rock yield, a cliff fall appears because the statistical damage is considered.The damage on the weak part of the rock accumulates and the mechanical parameters are continuously reduced until rock failure, which explains why although the intact rock is very strong, the rock mass is destroyed along the weak structural plane during practical engineering.

Engineering Application
To validate the rationality of the SSD model proposed above, a comparative numerical calculation was carried out to analyze the deformation of surrounding rock at a test roadway in the Jinchuan II mining area based on the MC model, SS model, and SSD model, respectively.

Geological Background
Jinchuan Mine, the largest nickel production base in China, is located in Jinchang City, Gansu Province.The ore deposit is about 6.5 km long, 200 m wide, and over 1000 m deep.It is a large and complex deposit with high ground stress, deep buried ore bodies, and broken rocks.The engineering  Figure 4a shows an ideal elastoplastic curve under the MC model-the stress increases as the strain increases until the yield strength is reached.Under the SS model shown in Figure 4b, the stress can be divided into three stages: monotonic increasing before the critical strain, decline, and basic stability.The yield stresses under the MC model and the SS model are almost the same-about 45 MPa.It can be seen from Figure 4c that the SSD model proposed in this paper can reflect the full process of failure for rock, especially the strain-softening property.The stress-strain curve is basically the same as the previous two in the elastic stage, except that the yield stress is slightly smaller.However, after the rock yield, a cliff fall appears because the statistical damage is considered.The damage on the weak part of the rock accumulates and the mechanical parameters are continuously reduced until rock failure, which explains why although the intact rock is very strong, the rock mass is destroyed along the weak structural plane during practical engineering.

Engineering Application
To validate the rationality of the SSD model proposed above, a comparative numerical calculation was carried out to analyze the deformation of surrounding rock at a test roadway in the Jinchuan II mining area based on the MC model, SS model, and SSD model, respectively.

Geological Background
Jinchuan Mine, the largest nickel production base in China, is located in Jinchang City, Gansu Province.The ore deposit is about 6.5 km long, 200 m wide, and over 1000 m deep.It is a large and complex deposit with high ground stress, deep buried ore bodies, and broken rocks.The engineering geology of this area is shown in Figure 5. Due to its extremely complex regional geological background, special engineering geological environment, and distinctive rock mechanics, Jinchuan Mine has attracted the attention of scholars and experts at home and abroad.With the mining depth increasing, the safety of deep roadways has become a prominent problem in the Jinchuan mining area [35][36][37].geology of this area is shown in Figure 5. Due to its extremely complex regional geological background, special engineering geological environment, and distinctive rock mechanics, Jinchuan Mine has attracted the attention of scholars and experts at home and abroad.With the mining depth increasing, the safety of deep roadways has become a prominent problem in the Jinchuan mining area [35][36][37].The simulated test roadway in this paper is located in mining area II at a depth of about 750 m, as shown in Figure 6 The surrounding rock is mainly gray-white marble with a massive structure and a medium-fine grain crystalloblastic texture.The main mineral compositions include dolomite, calcite, serpentine, and a small amount of quartz stone.There is only Quaternary pore water and a small amount of structural fissure water in this mining area, so the hydrogeological conditions have little impact on the project.The intact rock and rock mass mechanical parameters in the tentative roadway are shown in Table 1.The value of the in situ stress near the test roadway is shown in Table 2, and there is a linear relationship between the three principal stresses and the depth as follows [38,39]:  The simulated test roadway in this paper is located in mining area II at a depth of about 750 m, as shown in Figure 6 The surrounding rock is mainly gray-white marble with a massive structure and a medium-fine grain crystalloblastic texture.The main mineral compositions include dolomite, calcite, serpentine, and a small amount of quartz stone.There is only Quaternary pore water and a small amount of structural fissure water in this mining area, so the hydrogeological conditions have little impact on the project.The intact rock and rock mass mechanical parameters in the tentative roadway are shown in Table 1.The value of the in situ stress near the test roadway is shown in Table 2, and there is a linear relationship between the three principal stresses and the depth as follows [38,39]: σ 1 = 8.620 + 0.018H σ 2 = 3.984 + 0.017H σ 3 = 0.077 + 0.00.024H(8) where σ 1 , σ 2 , σ 3 are the maximum principal stress, the intermediate principal stress, and the minimum principal stress (unit: MPa); H is the depth in m. geology of this area is shown in Figure 5. Due to its extremely complex regional geological background, special engineering geological environment, and distinctive rock mechanics, Jinchuan Mine has attracted the attention of scholars and experts at home and abroad.With the mining depth increasing, the safety of deep roadways has become a prominent problem in the Jinchuan mining area [35][36][37].The simulated test roadway in this paper is located in mining area II at a depth of about 750 m, as shown in Figure 6 The surrounding rock is mainly gray-white marble with a massive structure and a medium-fine grain crystalloblastic texture.The main mineral compositions include dolomite, calcite, serpentine, and a small amount of quartz stone.There is only Quaternary pore water and a small amount of structural fissure water in this mining area, so the hydrogeological conditions have little impact on the project.The intact rock and rock mass mechanical parameters in the tentative roadway are shown in Table 1.The value of the in situ stress near the test roadway is shown in Table 2, and there is a linear relationship between the three principal stresses and the depth as follows [38,39]:  The test roadway adopts a semi-circular arched shape because of the reasonable stress distribution and high utilization rate.The primary support is a composite supporting scheme including a double-layer bolt-mesh-anchor and U-shaped steel, as shown in Figure 7.The anchor rod with cement mortar is made of secondary rebar and has a length of 2250 mm, a diameter of 22 mm and spacing of 1000 × 1000 mm.The anchor plate is made of an ordinary steel plate with a thickness of 10 mm and a size of 200 × 200 mm.Concrete with a compressive strength is 25 MPa is used for shotcrete lining, and the thickness of a single layer is 100 mm.The U-shaped steel support uses U36 steel, and the specific parameters of each support structure are shown in Tables 3-5  The test roadway adopts a semi-circular arched shape because of the reasonable stress distribution and high utilization rate.The primary support is a composite supporting scheme including a double-layer bolt-mesh-anchor and U-shaped steel, as shown in Figure 7.The anchor rod with cement mortar is made of secondary rebar and has a length of 2250 mm, a diameter of 22 mm and spacing of 1000 × 1000 mm.The anchor plate is made of an ordinary steel plate with a thickness of 10 mm and a size of 200 × 200 mm.Concrete with a compressive strength is 25 MPa is used for shotcrete lining, and the thickness of a single layer is 100 mm.The U-shaped steel support uses U36 steel, and the specific parameters of each support structure are shown in Tables 3-5

Numerical Model Establishment
The numerical model, as shown in Figure 8, was established according to the actual standard size of the tentative roadway: the semicircle arch radius is 2 m, and the height of the straight wall is also 2 m.The dimensions of the entire model are 40 m in length, 20 m in width, and 40 m in height.The displacement boundary is set in six directions to control the normal displacement, and triaxial ground stresses are applied in the model.The primary support is simulated based on the structural units in the software: a beam unit for the U-shape steel, a cable unit for the anchor, and a Shell unit for the concrete spray layer.In this program, Equation ( 4) is an incremental distribution function between 0-1 and the value of the parameters only affects the shape of the curve, so the value is required to be within a reasonable range, namely, the value of λ is in the same order of magnitude as c, m and λ refers to the literature [27] and [42].Thus, λ and m were selected to be 1 × e 6 and 2 respectively, and m was taken as 5 in Equation (7).Six numerical simulations were carried out in this paper, using the MC, SS, and SSD models under no support and primary support, respectively.

Numerical Model Establishment
The numerical model, as shown in Figure 8, was established according to the actual standard size of the tentative roadway: the semicircle arch radius is 2 m, and the height of the straight wall is also 2 m.The dimensions of the entire model are 40 m in length, 20 m in width, and 40 m in height.The displacement boundary is set in six directions to control the normal displacement, and triaxial ground stresses are applied in the model.The primary support is simulated based on the structural units in the software: a beam unit for the U-shape steel, a cable unit for the anchor, and a Shell unit for the concrete spray layer.In this program, Equation ( 4) is an incremental distribution function between 0-1 and the value of the parameters only affects the shape of the curve, so the value is required to be within a reasonable range, namely, the value of λ is in the same order of magnitude as c, m and λ refers to the literature [27] and [42].Thus, λ and m were selected to be 1×e 6 and 2 respectively, and m was taken as 5 in Equation (7).Six numerical simulations were carried out in this paper, using the MC, SS, and SSD models under no support and primary support, respectively.

Unsupported Condition
The results of the horizontal and vertical displacement of the surrounding rock mass under the unsupported condition calculated using the MC model, SS model, and SSD model are shown in Figures 9-12.It can be seen from the Figures above that the excavation caused serious uneven deformation, including floor heave, roof subsidence and side shrinkage, in the surrounding rock of the roadway in the unsupported condition.The SSD model results of the maximum deformation of the sidewall, floor, and roof can be seen in Figure 11; the values are approximately 70, 100, and 100 cm, respectively.For horizontal displacement, the maximum value produced by the MC and SS models is uniformly distributed in a semi-elliptical way on the sidewall.However, the maximum value appeared near the shoulder of the arch, and the displacement distribution was layered in the SSD model, which coincides with the deformation and failure of the roadway near the study area, as shown in Figure 13(A) and (B).For the vertical displacement, the distribution characteristics are similar to the horizontal characteristics.Moreover, bigger deformations occurred on the roof, and the root of the wall in the SSD model was able to reflect the actual status of the roadway better, as shown in Figure 13(C) and (D).In order to compare the differences among the three models more directly, four monitoring points were set on the roof, floor, and two sides of the roadway.Figure 12 reveals that the displacement between two points increased the fastest in the SSD model, which illustrates that the strain-softening and statistical damage really work in this application.In summary, the SSD model can perform a better qualitative analysis of roadway deformation and failure.It can be seen from the Figures above that the excavation caused serious uneven deformation, including floor heave, roof subsidence and side shrinkage, in the surrounding rock of the roadway in the unsupported condition.The SSD model results of the maximum deformation of the sidewall, floor, and roof can be seen in Figure 11; the values are approximately 70, 100, and 100 cm, respectively.For horizontal displacement, the maximum value produced by the MC and SS models is uniformly distributed in a semi-elliptical way on the sidewall.However, the maximum value appeared near the shoulder of the arch, and the displacement distribution was layered in the SSD model, which coincides with the deformation and failure of the roadway near the study area, as shown in Figure 13(A) and (B).For the vertical displacement, the distribution characteristics are similar to the horizontal characteristics.Moreover, bigger deformations occurred on the roof, and the root of the wall in the SSD model was able to reflect the actual status of the roadway better, as shown in Figure 13(C) and (D).In order to compare the differences among the three models more directly, four monitoring points were set on the roof, floor, and two sides of the roadway.Figure 12 reveals that the displacement between two points increased the fastest in the SSD model, which illustrates that the strain-softening and statistical damage really work in this application.In summary, the SSD model can perform a better qualitative analysis of roadway deformation and failure.

Under Primary Support
The computation results of the displacement of surrounding rock mass under the primary support based on the MC, SS, and SSD models were obtained and are shown in Figures 14-16 in which the displacement is the result of the horizontal and vertical displacement.

Under Primary Support
The computation results of the displacement of surrounding rock mass under the primary support based on the MC, SS, and SSD models were obtained and are shown in Figures 14-16 in which the displacement is the result of the horizontal and vertical displacement.

Under Primary Support
The computation results of the displacement of surrounding rock mass under the primary support based on the MC, SS, and SSD models were obtained and are shown in Figures 14-16 in which the displacement is the result of the horizontal and vertical displacement.As can be seen from the figures above, the current support method effectively controls the convergence deformation of the roadway, and the overall deformation is reduced by about an order of magnitude.The deformation of the roadway in the first two models mainly involves horizontal convergence, while the deformation of the roadway in the SSD model is larger at the roof and floor.The values of displacement in the MC and SS models are bigger than in the other two models.
At the same time, field monitoring was carried out in the test roadway in order to measure the deformation of the roadway and compare it with the numerical simulation results.Five monitoring sections with spacings of 5 m were selected for the test roadway, and eight monitoring points were set on each section, as shown in Figure 17.The deformation of each section was periodically monitored using a laser range finder (Figure 18).In the past year, 17 measurements were made on the test roadway, which is sufficient to accurately quantify the development of the surrounding rock deformation.As can be seen from the figures above, the current support method effectively controls the convergence deformation of the roadway, and the overall deformation is reduced by about an order of magnitude.The deformation of the roadway in the first two models mainly involves horizontal convergence, while the deformation of the roadway in the SSD model is larger at the roof and floor.The values of displacement in the MC and SS models are bigger than in the other two models.
At the same time, field monitoring was carried out in the test roadway in order to measure the deformation of the roadway and compare it with the numerical simulation results.Five monitoring sections with spacings of 5 m were selected for the test roadway, and eight monitoring points were set on each section, as shown in Figure 17.The deformation of each section was periodically monitored using a laser range finder (Figure 18).In the past year, 17 measurements were made on the test roadway, which is sufficient to accurately quantify the development of the surrounding rock deformation.As can be seen from the figures above, the current support method effectively controls the convergence deformation of the roadway, and the overall deformation is reduced by about an order of magnitude.The deformation of the roadway in the first two models mainly involves horizontal convergence, while the deformation of the roadway in the SSD model is larger at the roof and floor.The values of displacement in the MC and SS models are bigger than in the other two models.
At the same time, field monitoring was carried out in the test roadway in order to measure the deformation of the roadway and compare it with the numerical simulation results.Five monitoring sections with spacings of 5 m were selected for the test roadway, and eight monitoring points were set on each section, as shown in Figure 17.The deformation of each section was periodically monitored using a laser range finder (Figure 18).In the past year, 17 measurements were made on the test roadway, which is sufficient to accurately quantify the development of the surrounding rock deformation.The average displacement curves of the monitoring points in the five sections are shown in Figure 19.I can be observed that the deformation process of the surrounding rock can be divided into three stages.Using 15 days and 200 days as the boundary points, the deformation experienced three stages: rapid growth, slow growth, and basic stability.At the last measurement, the maximum deformation occurred between points 4 and 8, which is the of the top and bottom plates, and that was about 15 cm.The average displacement curves of the monitoring points in the five sections are shown in Figure 19.I can be observed that the deformation process of the surrounding rock can be divided into three stages.Using 15 days and 200 days as the boundary points, the deformation experienced three stages: rapid growth, slow growth, and basic stability.At the last measurement, the maximum deformation occurred between points 4 and 8, which is the deformation of the top and bottom plates, and that was about 15 cm.The average displacement curves of the monitoring points in the five sections are shown in Figure 19.I can be observed that the deformation process of the surrounding rock can be divided into three stages.Using 15 days and 200 days as the boundary points, the deformation experienced three stages: rapid growth, slow growth, and basic stability.At the last measurement, the maximum deformation occurred between points 4 and 8, which is the deformation of the top and bottom plates, and that was about 15 cm.The average displacement curves of the monitoring points in the five sections are shown in Figure 19.I can be observed that the deformation process of the surrounding rock can be divided into three stages.Using 15 days and 200 days as the boundary points, the deformation experienced three stages: rapid growth, slow growth, and basic stability.At the last measurement, the maximum deformation occurred between points 4 and 8, which is the deformation of the top and bottom plates, and that was about 15 cm.Comparing the displacement curves of surrounding rock mass based on the three models with the monitoring results, we can see that the curves of the SSD model are most similar to the monitored ones, regardless of the distribution law or the value.In terms of the distribution law of displacement, both the SSD model and the field monitoring are characterized by vertical convergence rather than the horizontal convergence shown by the MC model and the SS model.Furthermore, the displacement values of the SSD model agree well with the measured values-the vertical convergence (4-8) values are about 18 and 15 cm, respectively, and the horizontal convergence (1-7 and 2-6) values are about 9 and 10 cm, Thus, the SSD model can perform a better quantitative analysis of roadway deformation and failure.

Conclusions
In order to study the stability of deep tunnels and to address the shortcomings of the current constitutive model, a strain-softening model of heterogeneous jointed rock mass considering statistical damage was developed and implemented through FLAC 3D simulation software.Then, based on a test roadway in the Jinchuan nickel mining area in China, numerical simulation calculations were carried out using the MC, SS, and SSD models; the computation results were compared with the field-monitoring results.
The results show that the SSD model simulates the behavior of the surrounding rocks well, and it can perform a better qualitative and quantitative analysis of roadway deformation and failure.Furthermore, the case study in the test roadway of the Jinchuan mine shows that the deformations of the roof and floor are larger under the primary support area, which may serve as a reference for the support pattern of deep roadways and provide a reference for simulation calculations for other similar projects.

1 3 ,
  are the maximum and minimum principal stresses; t  is the tensile strength;c C is the cohesion;  is the internal friction angle; and

Figure 1 .
Figure 1.Stress-strain relation of an element.

Figure 2 .
Figure 2. Program flowchart of the statistical damage (SSD) model.

Figure 3 .
Figure 3. Numerical model of the rock specimen.

Figure 4 .
Figure 4. Stain-stress-acoustic emission ringing count curves of different models.

Figure 4 .
Figure 4. Stain-stress-acoustic emission ringing count curves of different models.

Figure 4 .
Figure 4. Stain-stress-acoustic emission ringing count curves of different models.
are the maximum principal stress, the intermediate principal stress, and the minimum principal stress (unit: MPa); H is the depth in m.
are the maximum principal stress, the intermediate principal stress, and the minimum principal stress (unit: MPa); H is the depth in m.

Figure 6 .
Figure 6.Location of the test roadway.

Figure 7 .
Figure 7.Primary supporting method of the tentative roadway.

Figure 8 .
Figure 8. Numerical model of primary supporting.

Figure 8 .
Figure 8. Numerical model of primary supporting.

3. 3 .
Results Analysis 3.3.1.Unsupported Condition The results of the horizontal and vertical displacement of the surrounding rock mass under the unsupported condition calculated using the MC model, SS model, and SSD model are shown in Figures 9-12.

Figure 9 .
Figure 9. Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the MC model.

Figure 10 .
Figure 10.Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the SS model.

Figure 11 .
Figure 11.Horizontal and vertical displacement of surrounding rock mass under no support calculated by the SSD model.

Figure 9 . 19 Figure 9 .
Figure 9. Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the MC model.

Figure 10 .
Figure 10.Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the SS model.

Figure 11 .
Figure 11.Horizontal and vertical displacement of surrounding rock mass under no support calculated by the SSD model.

Figure 10 . 19 Figure 9 .
Figure 10.Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the SS model.

Figure 10 .
Figure 10.Horizontal and vertical displacement of the surrounding rock mass under no support calculated by the SS model.

Figure 11 .
Figure 11.Horizontal and vertical displacement of surrounding rock mass under no support calculated by the SSD model.

Figure 11 .
Figure 11.Horizontal and vertical displacement of surrounding rock mass under no support calculated by the SSD model.

Figure 12 .
Figure 12.Displacement curve of monitoring points in the roadway under no support calculated by the three models.

Figure 12 .
Figure 12.Displacement curve of monitoring points in the roadway under no support calculated by the three models.

Figure 13 .
Figure 13.Deformation and failure of the roadway near the study area.

Figure 14 .
Figure 14.Displacement of the surrounding rock mass under current support calculated by the MC model.

Figure 13 .
Figure 13.Deformation and failure of the roadway near the study area: (A,B), surface cracking and layered peeling; (C) roof subsidence; (D) floor heave.

Sustainability 2019 , 19 Figure 13 .
Figure 13.Deformation and failure of the roadway near the study area.

Figure 14 .
Figure 14.Displacement of the surrounding rock mass under current support calculated by the MC model.

Figure 14 .
Figure 14.Displacement of the surrounding rock mass under current support calculated by the MC model.

Figure 15 .
Figure 15.Displacement of the surrounding rock mass under current support calculated by the SS model.

Figure 16 .
Figure 16.Displacement of the surrounding rock mass under current support calculated by the SSD model.

Figure 15 . 19 Figure 15 .
Figure 15.Displacement of the surrounding rock mass under current support calculated by the SS model.

Figure 16 .
Figure 16.Displacement of the surrounding rock mass under current support calculated by the SSD model.

Figure 16 .
Figure 16.Displacement of the surrounding rock mass under current support calculated by the SSD model.

Figure 19 .
Figure 19.Displacement curves of the monitoring points at the roadway surface under current support.

Figure 19 .
Figure 19.Displacement curves of the monitoring points at the roadway surface under current support.

Figure 19 .
Figure 19.Displacement curves of the monitoring points at the roadway surface under current support.

Figure 19 .
Figure 19.Displacement curves of the monitoring points at the roadway surface under current support.

Figure 20 .Figure 21 .
Figure 20.Displacement curves of surrounding rock mass under current support calculated by the MC

Figure 22 .
Figure 22.Displacement curves of surrounding rock mass under current support calculated by the SSD model.

Figure 20 .
Figure 20.Displacement curves of surrounding rock mass under current support calculated by the MC model.

Figure 20 .
Figure 20.Displacement curves of surrounding rock mass under current support calculated by the MC model.

Figure 21 .
Figure 21.Displacement curves of surrounding rock mass under current support calculated by the SS model.

Figure 22 .
Figure 22.Displacement curves of surrounding rock mass under current support calculated by the SSD model.

Figure 21 .
Figure 21.Displacement curves of surrounding rock mass under current support calculated by the SS model.

Figure 20 .
Figure 20.Displacement curves of surrounding rock mass under current support calculated by the MC model.

Figure 21 .
Figure 21.Displacement curves of surrounding rock mass under current support calculated by the SS model.

Figure 22 .
Figure 22.Displacement curves of surrounding rock mass under current support calculated by the SSD model.

Figure 22 .
Figure 22.Displacement curves of surrounding rock mass under current support calculated by the SSD model.

Table 1 .
Intact rock properties and rock mass properties of the tentative roadway.

Table 1 .
Intact rock properties and rock mass properties of the tentative roadway.

Table 1 .
Intact rock properties and rock mass properties of the tentative roadway.

Table 2 .
Ground stress in the test section.

Table 3 .
Parameters of the concrete lining.

Table 2 .
Ground stress in the test section.

Table 3 .
Parameters of the concrete lining.

Table 4 .
Parameters of the anchor bolt.

Table 5 .
Parameters of the U-shaped steel supports.