Abstract

The obvious anisotropy of layered rock mass prevents the built-in constitutive model of commercial simulation software from accurately describing the elastoplastic behavior of the rock mass. In this paper, a transversely isotropic elastoplastic constitutive model (called the AN-MC model) that characterizes the elastoplastic failure behavior of the layered rock mass is developed by introducing the Mohr-Coulomb yield criterion with the tension cutoff in the model based on the transversely isotropic constitutive model and deriving its finite difference scheme. The development and programming of the transversely isotropic elastoplastic constitutive model are achieved on the basis of the secondary development platform of FLAC3D and the VC++ environment. The accuracy and rationality of the proposed constitutive model are verified in terms of consistency of the calculation results of the developed transversely isotropic elastoplastic constitutive model and the built-in constitutive models. In the numerical analysis of the Queerxi tunnel project, the calculation results of the AN-MC model are in good agreement with the field monitoring data. Thus, the deformation characteristics of the tunnel surrounding rock are well characterized. Comparative analysis of the deformation and failure laws of the layered surrounding rock and isotropic surrounding rock indicates that the plastic failure range of the layered surrounding rock tunnel is larger, and the arch foot and waist are prone to plastic penetration failure with butterfly-shaped characteristics. The research results are important for the scientific prediction and precise support of layered rock mass deformation.

1. Introduction

Layered rock mass has a stratified and rhythmic structure with obvious transverse isotropy characteristics. In addition, the stress distribution, deformation, and failure characteristics of surrounding rock in layered rock mass are complex. Therefore, constructing a set of constitutive models that can scientifically characterize the elastoplastic mechanical behavior of layered rock mass has important theoretical significance and engineering value for analyzing the mechanical behavior and the deformation and failure mechanisms of layered rock mass.

In recent years, extensive research has been conducted on the deformation and failure characteristics of layered rock mass [17]. Xia et al. [1] studied the bearing deformation and failure process of layered rock mass through a model test and revealed the eccentric pressure failure phenomenon of surrounding rock. Liu et al. [3] studied the correlation characteristics between the bedding angle of sandstone and the development rate of cracks based on a laboratory test and acoustic emission technology. However, the size effect of the laboratory test prevents objective and true reflection of the structure of rock mass. Numerical simulations provide an alternative approach, and some innovative conclusions were obtained [819]. Ding et al. [8] and Tang and Sun [9] used an isotropic body with a bedding plane to characterize the anisotropy of layered rock mass based on the Mohr-Coulomb constitutive model. They then simulated and calculated the stability of a layered rock mass tunnel. Liang et al. [14] and Xia [15] studied the fracture mode and failure criterion of layered rock mass with different bedding plane inclination angles through numerical experiments and found the included angle between rock formation and the maximum principal stress direction to be the key factor that affects the failure process and range. Wu et al. [19] found that the failure surface was mainly developed along the weak surface of soft and hard stratification, and the interlayer failure was mainly caused by the failure of weak stratification.

A constitutive model is the key to simulate the mechanical deformation and failure characteristics of rock mass. Considerable research has been conducted on the constitutive model of layered rock mass [2026]. Han et al. [23] and Xie et al. [24] used the Mohr-Coulomb yield criterion to describe the failure characteristics of rock mass and established a transversely isotropic elastoplastic constitutive relationship. However, they did not consider the tensile failure of rock mass. Liu [25] used a ubiquitous-joint constitutive model to describe the anisotropy of layered rock mass. Zhou et al. [26] developed an enhanced equivalent continuum constitutive model based on the ubiquitous-joint constitutive model considering the transverse isotropy characteristics and tensile shearing yield criterion of layered rock mass. However, the transversely isotropic characteristic model was adopted even for the failed rock mass. The secondary development and programming of the constitutive model is an important approach to analyze special geotechnical problems [2731]. Chen and Liu [27] realized the development of the Duncan Chang constitutive model based on the secondary development platform of FLAC3D. Jiang et al. [29] and Zuo et al. [30] adopted the Mohr-Coulomb yield criterion and the Zienkiewicz-Pande yield criterion, respectively, to construct rock elastoplastic damage constitutive models.

In this paper, based on the transversely isotropic constitutive model, a transversely isotropic elastoplastic constitutive model is developed by introducing the Mohr-Coulomb (MC) yield criterion with the tension cutoff in the model; this constitutive model is called the AN-MC model. Compared with the transversely isotropic constitutive model and the MC model, the AN-MC model can both express the anisotropic characteristics of rock mass and be used for elastoplastic analysis of rock mass. The secondary development and programming of the constitutive model are realized on the basis of the FLAC3D platform. The proposed AN-MC constitutive model is applied to the Queerxi tunnel project to further verify the models’ accuracy, and the deformation and failure characteristics of layered surrounding rock and stability of layered surrounding rock are deeply analyzed. The research results have universal applicability for analyzing the deformation and failure characteristics of layered rock mass.

2. Implementation of the AN-MC Model in FLAC3D

2.1. Introduction to the AN-MC Model

Based on the transversely isotropic constitutive model, the AN-MC model is proposed through introducing the MC yield criterion with the tension cutoff in the model. Compared with the transversely isotropic model, which considers only the elastic analysis function, and the MC model, which assumes an isotropic body with the elastoplastic analysis function, the AN-MC model has the following two remarkable advantages: The constitutive model can ① express the anisotropic characteristics of rock mass and ② be used for elastoplastic analysis of rock mass. The equations of the AN-MC model are described in Table 1.

2.1.1. Elastic Relationship and Yield Criterion

The stress-strain relationship in the elastic stage of the AN-MC model is consistent with the elastic relationship of the transversely isotropic constitutive model [15]. The proposed model adopts the MC yield criterion with the tension cutoff, which is described by the shear and tensile yield criteria [32].

2.1.2. Plastic Potential Functions

The plastic potential function is described using two functions, and , which represent shear and tensile plastic flow, respectively.

The function corresponds to a nonassociated law and has the following form: where is the dilation angle.

The function corresponds to an associated flow rule and is written as

Because of the anisotropic characteristics of layered rock mass, cohesion and friction in Table 1 are variables that vary with the direction change. The minimum cohesion and minimum friction of layered rock mass are in the plane parallel to the bedding plane, and the maximum cohesion and maximum friction are in the plane perpendicular to the bedding plane and increases as () increases. The values of and are calculated by the expression proposed by Zhang and Liu [33] where is the included angle between the plane of values of and and the bedding plane.

As shown in Figure 1, is the global coordinate system. The local coordinate system (i.e., the normal coordinate system of the bedding plane) is obtained by rotating an angle counterclockwise with the -axis as the central axis; then, is the bedding plane, -axis is the normal direction of the bedding plane, -axis is the dip direction angle, and is the dip angle. Further, the dip direction angle and dip angle are expressed by dd and , respectively.

In the local coordinate system , the elastic matrix in Table 1 is expressed as follows: where , , and . , , , , and are the five independent elastic constants of the transversely isotropic constitutive model. and are the elastic modulus and Poisson’s ratio of the plane parallel to the bedding plane (the plane of isotropy), respectively; , , and are the elastic modulus, Poisson ratio, and shear modulus in the plane perpendicular to the bedding plane (normal to the plane of isotropy), respectively.

2.2. Finite Difference Scheme of the AN-MC Model

According to the elastoplastic theory, the stress increment in elastic deformation can be solved by the strain increment: this process is called elastic guess. A linear relationship exists between the strain increment and the stress increment. The stress increment is expressed as follows: where is a linear function of elastic strain increments .

When the rock mass enters the plastic stage, the strain increment of the rock mass is composed of elastic and plastic strain increments. The strain increment is expressed in where , , and represents the strain and elastic strain and plastic strain increments, respectively.

The incremental expression of the stress strain relationship of the elastic part of the AN-MC model is obtained by the extension of the incremental expression of Hooke’s law. It is expressed in where represents the material parameter.

According to Equation (5) and symmetry of the elastic matrix, Equation (7) is rewritten as

The plastic strain part follows the relevant flow law [3436]: where is the plastic factor and g is the plastic potential function.

The proposed model includes the shear plastic flow law and tensile plastic flow law which are expressed as follows: where and are shear and tensile plastic factors, respectively, and and are the shear and tensile plastic potential functions, respectively.

After substituting the partial derivative of Equation (1) into Equation (10), the plastic shear strain in each principal direction can be obtained:

The plastic shear strain increment can also be expressed as follows:

After substituting the partial derivative of Equation (2) into Equation (10), the plastic tensile strain in each principal direction can be obtained:

The plastic tensile strain increment can also be expressed as follows:

When plastic deformations of rock mass are involved, only the elastic part of the stress increment will contribute to the stress increment. In this case, the plastic correction must be made to the elastic stress increment as computed from the total strain increment in order to obtain the real stress state.

The expression of the elastic guess is shown in where is the new stress component and is the elastic guess.

The stress plastic correction is expressed as follows:

From Equations (8), (12), and (16) and replacing , , and with , , and , stress of each principal direction after shear plastic correction is obtained as follows:

From Equations (8), (14), and (16) and replacing , , and with , , and , stress of each principal direction after tensile plastic correction is obtained as follows:

2.3. Implementation of the AN-MC Model

All constitutive models in FLAC3D share the same incremental numerical algorithm. Given the stress state at time and the total strain increment for a timestep , the stress state at time is determined [37, 38].

The following four cyclic steps are involved in this process: ① Based on Newton’s second law, the nodal velocities are solved according to the initial stress and boundary conditions. ② The element strain rate is solved from the nodal velocities based on Gauss’s law. ③ Based on the material constitutive model, the element stress and state variables are updated according to the element strain rate. ④ The nodal force is obtained by integrating the updated stress on the element. Figure 2 shows the calculation cycle diagram. The constitutive model in the third step is obtained by the secondary development in this paper.

The AN-MC constitutive model is written in C++ and is compiled into a dynamic link library file for it to be loaded and called by FLAC3D. The development of the constitutive model needs to modify the header file (.H file) and the program file (.CPP file). The declaration of the derived class, the definition of the function name, and the modification of the private member of the derived class in the new constitutive model are performed in the header file. The new constitutive model is implemented in the program file [27, 39]. Based on Figure 2, the AN-MC program files are explained as follows. (1)The functions getProperties(), getProperty(), and setProperty() are used to assign model parameters. The function getProperties() returns the string of the parameter name recognized by the calculation command of FLAC3D. The function getProperty() returns the value of the model parameter. The function setProperty() establishes a one-to-one correspondence between the string of the parameter name and the value of the model parameter. The parameters of the AN-MC model, including young-plane, young-normal, Poisson-plane, Poisson-normal, shear-normal, mincohesion, maxcohesion, minfriction, maxfriction, dilation, tension, dip, and dd, are developed(2)The function Initialize() is mainly used to initialize the model parameters and the state variables. This function is called once for each element when the CYCLE command is executed(3)The function run() is the main part of the constitutive model, and its main function is to calculate the new stress increment according to the strain increment and to update the stress and state variables. This function realizes the functions of elastic guess, yield criterion, flow law, and plastic correction(4)The function getStates() is the yield state function of the constitutive model. Based on the MC yield criterion, four yield states are set for the element in this paper, namely, experiencing shear instability (shear-n), experiencing tension instability (tension-n), shear instability in the past (shear-p), and tension instability in the past (tension-p). If the element does not yield, the yield state will be shown as “None”

3. Model Verification of the AN-MC Model

The verification of the constitutive model of engineering materials involves the following three methods: comparison with the built-in constitutive model of commercial numerical analysis software, comparison with the laboratory mechanical test, and engineering data and numerical calculation results of self-developed constitutive model. The AN-MC model is used to calculate the stress of surrounding rock and the elastoplastic deformation of surrounding rock during the excavation of underground engineering. The results are compared with the calculation results of the built-in transversely isotropic constitutive model and the MC constitutive model in FLAC3D.

3.1. Introduction of Numerical Calculation Model

The size of the numerical calculation model is , and it is divided into 29160 finite elements. The excavation section is round with a radius of 1.5 m and an excavation depth of 30 m. The excavation method is full-section excavation, and the excavation position is the geometric center of the calculation area. First, gravity is applied to the calculation model (with the acceleration of gravity of 10 m/s2), and the displacement and velocity are set to zero. The AN-MC model, the transversely isotropic constitutive model, and the MC constitutive model are applied to the calculation model. A uniform load of 10 MPa is applied in the , , and directions. Excavation is performed after the displacement, velocity, and plastic zone state are set to zero.

3.2. Elastic Verification of the AN-MC Model

The comparison of calculation between the AN-MC model and the transversely isotropic constitutive model is realized by the method of model degradation [40]. By artificially increasing the plastic parameters of the AN-MC model, the plastic deformation of rock mass is prevented. The values of rock mass parameters of the AN-MC model and transversely isotropic constitutive model are shown in Table 2. The surrounding rock stress contour plot, displacement contour plot, and the plastic zone after tunnel excavation are presented in Figure 3. The numerical calculation results indicate that the compressive stress values of the AN-MC model, and the transversely isotropic constitutive model are in the ranges of 1.0419 MPa-4.9914 MPa and 1.0418 MPa-4.9914 MPa, respectively, and the maximum displacements are 7.6427 and 7.6421 mm, respectively. Moreover, no plastic zone is observed. The relative error values of maximum and minimum stress are 0 and 100 Pa between the AN-MC model and the transversely isotropic constitutive model. The relative error values of maximum and minimum displacement are 0.0006 and -0.000092 mm between the AN-MC model and the transversely isotropic constitutive model. The calculation results of the two constitutive models are quite similar and the errors are acceptable; therefore, the AN-MC model can degenerate into the transversely isotropic constitutive model. Thus, the elastic verification conclusion is correct and reasonable.

3.3. Plastic Verification of the AN-MC Model

Constant values of the elastic modulus and Poisson’s ratio are considered in all directions in the AN-MC model to verify whether the MC yield criterion is successfully embedded in the AN-MC model. The values of rock mass parameters of the AN-MC model and MC constitutive model are shown in Table 3. The surrounding rock stress contour plot, displacement contour plot, and plastic zone after tunnel excavation are shown in Figure 4. The numerical results indicate that the compressive stress values of the AN-MC model and the MC constitutive model ranges within 1.5137 MPa-7.2357 MPa and 1.5155 MPa-7.2361 MPa, respectively; the maximum displacements are 9.8666 and 9.8708 mm; and the shape of the plastic zone is similar. The relative error values of maximum and minimum stress are -400 and -1800 Pa between the AN-MC model and the MC constitutive model. The relative error values of maximum displacement and minimum displacement are -0.0042 mm and -0.00004757 mm between the AN-MC model and the MC constitutive model. The calculation results of the two constitutive models are very similar, and the errors are acceptable; therefore, the AN-MC model can degenerate into the MC constitutive model, and thus, the plastic verification of the model is completed.

4. Engineering Verification of the AN-MC Model and Its Application

4.1. Engineering Situations and Establishment of a Numerical Calculation Model

The Queerxi Tunnel in Huxi County, Hunan Province, China, is 1417 m long and passes through the low mountains. The mountain terrain has large fluctuation and deep cutting. No large geological structures such as faults and folds are found in the tunnel area, and the surface water system is undeveloped. The surface coverage of the tunnel area is extremely thin, with most of the rocks being exposed. The underlying bedrock is mainly interbedded with calcareous sandstone and calcareous argillaceous sandstone, with poor interlayer bonding.

The excavation method of this tunnel is full-section excavation, and the uneven deformation of tunnel surrounding rock is obvious. The surrounding rock deformation during the tunnel excavation has been studied by field monitoring and numerical simulation [41]. In this study, the AN-MC model is applied to analyze the deformation and failure characteristics of layered surrounding rock during tunnel excavation, and the accuracy and practicability of the AN-MC model are tested.

The sizes of the numerical model in the , , and directions are 150 m, 100 m, and 150 m, respectively. The tunnel cross-section has a horseshoe shape; the maximum excavation span from the arch vault to the inverted arch is 12 m, and the span between the left sidewall and the right sidewall is 12 m. The eastern, western, southern, northern, and lower boundaries of the numerical calculation model are set as velocity constraints, and the constraint value of each boundary is 0. Free boundary is set at the upper part of the numerical calculation model. The mechanical parameters related to the numerical model are shown in Table 4 [41]. The sizes and boundary conditions of the numerical model in the - plane are shown in Figure 5, where (a)–(e) are monitoring points for displacement.

4.2. Result Analysis
4.2.1. Displacement Results and Analysis

The relationship between the numerical simulation results and the field monitoring displacement [41] is shown in Figure 6. The settlement values of the tunnel vault calculated in reference [41] and this paper are 10.0 and 10.5 mm, respectively, and the settlement range obtained by field monitoring is 8.0 mm-10.0 mm; these values are in good agreement. However, considerable differences exist in the horizontal convergence displacement of the tunnel sidewall in three ways. The horizontal convergence displacement values of the tunnel sidewall calculated in reference [41] and this paper are 8.0 and 5.2 mm, respectively, and those obtained from field monitoring range within 4.5 mm-5.5 mm. The displacement calculated by the AN-MC model agrees well with the field monitoring data.

The calculation results of the two constitutive models are different, which is due to the different failure criteria. The combined discontinuity criterion based on the MC criterion, the weak plane strength theory, and the maximum axial strain theory is used to determine the failure of rock mass [41]. The MC criterion is only applicable to the elastoplastic analysis of isotropic rock mass. The anisotropic body is assumed to be an isotropic body with a set of directional weak planes in the weak plane strength theory. Further, the MC criterion is used to describe the shear failure of the bedrock and slip failure of the weak plane. The maximum axial strain theory is based on the weak plane strength theory and reflects the anisotropic degree of rock by anisotropic parameters. The premise of the modified transversely isotropic failure criterion is that the layered rock mass is isotropic at the nonweak plane, which makes most of the elements of the rock mass isotropic. Thus, the numerical simulation results are not precise, because isotropic elements cannot accurately express the difference between the mechanical properties in the parallel and perpendicular directions to the bedding plane. The AN-MC model assigns the mechanical parameters of layered rock mass to each element. Each element is transversely isotropic, which can express the different mechanical properties between the parallel and perpendicular directions to the bedding plane. It explains why there is a great difference between the settlement values of the tunnel vault and the horizontal convergence displacement calculated by the AN-MC model. The displacements of numerical results are consistent with the field monitoring displacements; this further verifies the accuracy and rationality of the AN-MC model and explains the engineering application value of the AN-MC model.

4.2.2. Failure Characteristics of Layered Surrounding Rock

This part discusses the development and failure state of the plastic zone of surrounding rock during tunnel excavation.

After the excavation of layered rock mass, the tunnel sidewalls and the inverted arch first enter the plastic state (Figures 7(a) and 7(b)). Subsequently, the plastic zone develops from the sidewalls to the arch vault, and the middle of the sidewalls and the arch feet on both sides extend to the internal surrounding rock in the form of the shear plastic zone (Figures 7(c)–7(f)). Then, the shear plastic zone extends symmetrically up and down along the middle of the sidewalls and reaches a heading through with the plastic zone at the arch waist and foot (Figures 7(g)–7(i)); this shows the butterfly-shaped failure characteristics. Figure 7 shows that the layered surrounding rock is dominated by the shear failure zone, and only tensile failure occurs in the arch vault to the inverted arch. Both sides of the arch waist present shear tensile composite failure. The initial shear failure of the arch foot on both sides is gradually transformed into shear tensile composite failure.

5. Discussion

(1)In the secondary development process of the AN-MC model, the transformation between the transversely isotropic material parameters (, , and ) of rock mass in Equation (7) and material parameters and of the MC yield criterion is complex. In this paper, the material parameters and in the incremental expression of Hooke’s law are transformed into material parameters and of the AN-MC model, and the incremental expression of transversely isotropic of layered rock mass is obtained. The key steps are as follows: ① and [30] in the incremental expressions of Hooke’s law are converted into and , where and are the bulk and shear moduli of the isotropic material, and and are the bulk and shear moduli of the transversely isotropic material, respectively; that is, and are replaced with and in the position of the MC criterion in the AN-MC model. ② The five independent elastic parameters of the transversely isotropic constitutive model—namely, the elastic modulus in the direction parallel to the bedding plane and that in the direction perpendicular to the bedding plane , Poisson’s ratio in the direction parallel to the bedding plane and that in the direction perpendicular to the bedding plane , and the shear modulus in the direction perpendicular to the bedding plane —are transformed into two elastic parameters and . The transformed program expression is shown in references [42], and their logical relationships are shown in Figures 8 and 9. and are obtained by the uniaxial compression test of layered rock mass with a dip angle of 0°, and are obtained by the uniaxial compression test of layered rock mass with dip angle of 90°, and is calculated by the empirical formula (1): (2)This part discusses the differences in deformation and failure characteristics of the tunnel surrounding rock between isotropic rock mass and layered rock mass. The transformation from anisotropic to isotropic is achieved by taking the same values of the elastic modulus and Poisson’s ratio of the AN-MC model. The parameter values of the isotropic rock mass are shown in Table 4. The elastic moduli and are , and Poisson’s ratios and are 0.3. Failure characteristics of surrounding rock during tunnel excavation are calculated on the basis of the isotropic constitutive model, as shown in Figure 10

After the tunnel is excavated, the sidewalls and inverted arch of the rock mass first enter the plastic state (Figures 10(a) and 10(b)). Subsequently, the plastic zone develops from the sidewalls to the arch vault, and the middle of the sidewalls extend to the internal surrounding rock in the form of the shear plastic zone (Figures 10(c)–10(f)). Then, the shear plastic zone continues to develop further to the surrounding rock behind the sidewalls, and the arch feet on both sides and the inverted arch extend slightly into the surrounding rock in the forms of shear and tensile plastic zones, respectively, (Figures 10(g)–10(i)); they show the bat-shaped failure characteristics. Figure 10 indicates that the surrounding rock of the isotropic tunnel is dominated by the shear failure zone, and only tensile failure occurs in the arch vault to the inverted arch. Both sides of the arch waist present shear tension composite failure. The initial tension failure of the arch foot on both sides is gradually transformed into shear and shear tensile composite failure.

Comparison of Figures 7 and 10 indicates that the development process and failure state changes of the plastic zone are consistent between isotropic surrounding rock and layered surrounding rock. The plastic zone of the sidewalls is large, whereas that of the arch vault is small and is the tensile plastic zone. Compared with the isotropic surrounding rock, the plastic zone behind the sidewalls of the tunnel has a larger expansion range; the plastic zone behind the sidewalls reaches a heading through with the plastic zone at the arch waist and foot. This phenomenon can be attributed to the fact that the mechanical parameters of the layered surrounding rock in the direction perpendicular to the bedding plane are weaker than those parallel to the bedding plane. This is the risk point of layered rock mass engineering as well as it is also the fundamental starting point of the research on the transversely isotropic constitutive model in this paper.

6. Conclusions

(1)Based on the transversely isotropic constitutive model, the transversely isotropic elastoplastic constitutive model named AN-MC model is developed through introducing the MC yield criterion with the tension cutoff in the model. The innovation of the AN-MC model is that it can be used for elastoplastic analysis for anisotropic rock, and the correctness of AN-MC model is verified by the method of model degradation(2)The results calculated by the AN-MC model were compared, respectively, with existing simulation calculation results and field monitoring data based on the Queerxi tunnel project, and the results show that the calculation results of the AN-MC model are closer to the actual engineering situation, which verifies the correctness and rationality of the proposed constitutive model(3)The isotropic surrounding rock and the layered surrounding rock have a consistent plastic characteristic in the excavation unloading process. However, the plastic zone behind the sidewalls of layered surrounding rock has a larger expansion range than isotropic rock. The plastic zone among tunnel sidewalls, tunnel arch foot, and tunnel arch waist of the surrounding rock is connected and preformed as the butterfly-shaped, and the isotropic surrounding rock tunnel preformed as the bat-shaped(4)The secondary development of rock constitutive model is an important research method to analyze the stress and deformation characteristics of geotechnical engineering. Optimization and development on the basis of AN-MC model need be used to analyze the failure development process of layered rock mass with the water environment

Data Availability

All data are available on request.

Conflicts of Interest

The authors declare no conflict of interests related to this article.

Acknowledgments

This research was supported by the State Key Laboratory of Structural Mechanics Behavior and System Safety of Traffic Engineering Open Subjects (No. ZZ2021-06), the National Natural Science Foundation of China (No. U2034207), the Hebei Province Key Research and Development Program of China (No. 21374101D and 19210804D), the Science and Technology Research Project of China National Railway Group Corporation (No. K2020G034), the Top Young Talents of Hebei Education Department (No. BJ2019011), and the Natural Science Foundation of Hebei Province (D2021210006). The authors are thankful for all of the support for this basic research.