# Development of Transversely Isotropic Elastoplastic Constitutive Model in FLAC3D and Its Application in Tunnel Engineering

Authors: Haoran Li; Yajun Zhang; Hongwu Yin; Xinlei Wang; Weigang Zhao; Wenjiang Li
Journal: Geofluids
 (2022)

---

## 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.


---

## Body

## 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 [1–7]. 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 [8–19]. 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 [20–26]. 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 [27–31]. 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 Table1.Table 1
Basic equation of the AN-MC constitutive model.
AN-MC modelFunctionsElastic relationshipσ′=D′ε′Yield criteriaShear yield criterion:fs=σ1−σ3Nϕ+2c√NϕTensile yield criterion: ft=σ3−σtPlastic potential functionsShear plastic potential function:gs=σ1−σ3NψTensile plasticity potential function: gt=−σ3


#### 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,gs and gt, which represent shear and tensile plastic flow, respectively.The functiongs corresponds to a nonassociated law and has the following form:
(1)gs=σ1−σ3Nψ,where ψ is the dilation angle.The functiongt corresponds to an associated flow rule and is written as
(2)gt=−σ3.Because of the anisotropic characteristics of layered rock mass, cohesionc and friction ϕ in Table 1 are variables that vary with the direction change. The minimum cohesion cmin and minimum friction ϕmin of layered rock mass are in the plane parallel to the bedding plane, and the maximum cohesion cmax and maximum friction ϕmax are in the plane perpendicular to the bedding plane and increases as θ (0°≤θ≤90°) increases. The values of c and ϕ are calculated by the expression proposed by Zhang and Liu [33]
(3)c=cmin+cmax−cminsinθϕ=ϕmin+ϕmax−ϕminsinθ,where θ is the included angle between the plane of values of c and ϕ and the bedding plane.As shown in Figure1, oxyz is the global coordinate system. The local coordinate system (i.e., the normal coordinate system of the bedding plane) ox′y′z′ is obtained by rotating an angle counterclockwise with the x-axis as the central axis; then, ox′y′ is the bedding plane, z′-axis is the normal direction of the bedding plane, y′-axis is the dip direction angle, and α is the dip angle. Further, the dip direction angle and dip angle are expressed by dd and dip, respectively.Figure 1
Relationship between local and global coordinate systems.In the local coordinate systemox′y′z′, the elastic matrix D′ in Table 1 is expressed as follows:
(4)D′=M1−nμ22μ1+nμ22μ21+μ1000μ1+nμ221−nμ22μ21+μ1000μ21+μ1μ21+μ11−μ22n000000G1M000000G2M000000G2M,where n=E2/E1, M=E1/1+μ11−μ1−2nμ22, and G1=E1/21+μ1. E1, E2, μ1, μ2, and G2 are the five independent elastic constants of the transversely isotropic constitutive model. E1 and μ1 are the elastic modulus and Poisson’s ratio of the plane parallel to the bedding plane (the plane of isotropy), respectively; E2, μ2, and G2 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:(5)Δσi=SiΔεne,where Si is a linear function of elastic strain increments Δεen.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(6)Δεi=Δεie+Δεip,where Δεi, Δεei, and Δεpi represents the strain and elastic strain and plastic strain increments, respectively.The incremental expression of the stress strain relationshipσ=Dε 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
(7)Δσi=αijΔεje,where αi represents the material parameter.According to Equation (5) and symmetry of the elastic matrix, Equation (7) is rewritten as
(8)S1Δε1e,Δε2e,Δε3e=α11Δε1e+α12Δε2e+α13Δε3eS2Δε1e,Δε2e,Δε3e=α12Δε1e+α22Δε2e+α23Δε3eS3Δε1e,Δε2e,Δε3e=α13Δε1e+α23Δε2e+α33Δε3e.The plastic strain part follows the relevant flow law [34–36]:
(9)Δεip=λ∂g∂σi,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:(10)Δεnps=λs∂gs∂σnΔεnpt=λt∂gt∂σn,where λs and λt are shear and tensile plastic factors, respectively, and gs and gt 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:
(11)Δε1ps=λs∂gs∂σ1=λsΔε2ps=λs∂gs∂σ2=0Δε3ps=λs∂gs∂σ3=−Nψλs.The plastic shear strain increment can also be expressed as follows:(12)dγps=λs∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=λsΔε1ps,Δε2ps,Δε3ps.After substituting the partial derivative of Equation (2) into Equation (10), the plastic tensile strain in each principal direction can be obtained:
(13)Δε1pt=λt∂gt∂σ1=0Δε2pt=λt∂gt∂σ2=0Δε3pt=λt∂gt∂σ3=−λt.The plastic tensile strain increment can also be expressed as follows:(14)dγpt=λt∂gt∂σ1,∂gt∂σ2,∂gt∂σ3=λtΔε1pt,Δε2pt,Δε3pt.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(15)σiN=σi+Δσi,σiI=σi+SiΔεn,where σNi is the new stress component and σIi is the elastic guess.The stress plastic correction is expressed as follows:(16)σiN=σiI−λSi∂g∂σn.From Equations (8), (12), and (16) and replacing Δεe1, Δεe2, and Δεe3 with Δεps1, Δεps2, and Δεps3, stress of each principal direction after shear plastic correction is obtained as follows:
(17)σ1N=σ1I−λsα11−α13Nψσ2N=σ2I−λsα12−α23Nψσ3N=σ3I−λsα13−α33Nψ,λs=fsσ1I,σ3Iα11−α13Nψ−α13−α33NψNϕ.From Equations (8), (14), and (16) and replacing Δεe1, Δεe2, and Δεe3 with Δεpt1, Δεpt2, and Δεpt3, stress of each principal direction after tensile plastic correction is obtained as follows:
(18)σ1N=σ1I−λtα13σ2N=σ2I−λtα23σ3N=σ3I−λtα33,λt=σt−σ3Iα33.


### 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 t and the total strain increment for a timestep Δt, the stress state at time t+Δt 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. Figure2 shows the calculation cycle diagram. The constitutive model in the third step is obtained by the secondary development in this paper.Figure 2
FLAC3D calculation process.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”


## 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 Table1.Table 1
Basic equation of the AN-MC constitutive model.
AN-MC modelFunctionsElastic relationshipσ′=D′ε′Yield criteriaShear yield criterion:fs=σ1−σ3Nϕ+2c√NϕTensile yield criterion: ft=σ3−σtPlastic potential functionsShear plastic potential function:gs=σ1−σ3NψTensile plasticity potential function: gt=−σ3


### 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,gs and gt, which represent shear and tensile plastic flow, respectively.The functiongs corresponds to a nonassociated law and has the following form:
(1)gs=σ1−σ3Nψ,where ψ is the dilation angle.The functiongt corresponds to an associated flow rule and is written as
(2)gt=−σ3.Because of the anisotropic characteristics of layered rock mass, cohesionc and friction ϕ in Table 1 are variables that vary with the direction change. The minimum cohesion cmin and minimum friction ϕmin of layered rock mass are in the plane parallel to the bedding plane, and the maximum cohesion cmax and maximum friction ϕmax are in the plane perpendicular to the bedding plane and increases as θ (0°≤θ≤90°) increases. The values of c and ϕ are calculated by the expression proposed by Zhang and Liu [33]
(3)c=cmin+cmax−cminsinθϕ=ϕmin+ϕmax−ϕminsinθ,where θ is the included angle between the plane of values of c and ϕ and the bedding plane.As shown in Figure1, oxyz is the global coordinate system. The local coordinate system (i.e., the normal coordinate system of the bedding plane) ox′y′z′ is obtained by rotating an angle counterclockwise with the x-axis as the central axis; then, ox′y′ is the bedding plane, z′-axis is the normal direction of the bedding plane, y′-axis is the dip direction angle, and α is the dip angle. Further, the dip direction angle and dip angle are expressed by dd and dip, respectively.Figure 1
Relationship between local and global coordinate systems.In the local coordinate systemox′y′z′, the elastic matrix D′ in Table 1 is expressed as follows:
(4)D′=M1−nμ22μ1+nμ22μ21+μ1000μ1+nμ221−nμ22μ21+μ1000μ21+μ1μ21+μ11−μ22n000000G1M000000G2M000000G2M,where n=E2/E1, M=E1/1+μ11−μ1−2nμ22, and G1=E1/21+μ1. E1, E2, μ1, μ2, and G2 are the five independent elastic constants of the transversely isotropic constitutive model. E1 and μ1 are the elastic modulus and Poisson’s ratio of the plane parallel to the bedding plane (the plane of isotropy), respectively; E2, μ2, and G2 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.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,gs and gt, which represent shear and tensile plastic flow, respectively.The functiongs corresponds to a nonassociated law and has the following form:
(1)gs=σ1−σ3Nψ,where ψ is the dilation angle.The functiongt corresponds to an associated flow rule and is written as
(2)gt=−σ3.Because of the anisotropic characteristics of layered rock mass, cohesionc and friction ϕ in Table 1 are variables that vary with the direction change. The minimum cohesion cmin and minimum friction ϕmin of layered rock mass are in the plane parallel to the bedding plane, and the maximum cohesion cmax and maximum friction ϕmax are in the plane perpendicular to the bedding plane and increases as θ (0°≤θ≤90°) increases. The values of c and ϕ are calculated by the expression proposed by Zhang and Liu [33]
(3)c=cmin+cmax−cminsinθϕ=ϕmin+ϕmax−ϕminsinθ,where θ is the included angle between the plane of values of c and ϕ and the bedding plane.As shown in Figure1, oxyz is the global coordinate system. The local coordinate system (i.e., the normal coordinate system of the bedding plane) ox′y′z′ is obtained by rotating an angle counterclockwise with the x-axis as the central axis; then, ox′y′ is the bedding plane, z′-axis is the normal direction of the bedding plane, y′-axis is the dip direction angle, and α is the dip angle. Further, the dip direction angle and dip angle are expressed by dd and dip, respectively.Figure 1
Relationship between local and global coordinate systems.In the local coordinate systemox′y′z′, the elastic matrix D′ in Table 1 is expressed as follows:
(4)D′=M1−nμ22μ1+nμ22μ21+μ1000μ1+nμ221−nμ22μ21+μ1000μ21+μ1μ21+μ11−μ22n000000G1M000000G2M000000G2M,where n=E2/E1, M=E1/1+μ11−μ1−2nμ22, and G1=E1/21+μ1. E1, E2, μ1, μ2, and G2 are the five independent elastic constants of the transversely isotropic constitutive model. E1 and μ1 are the elastic modulus and Poisson’s ratio of the plane parallel to the bedding plane (the plane of isotropy), respectively; E2, μ2, and G2 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:(5)Δσi=SiΔεne,where Si is a linear function of elastic strain increments Δεen.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(6)Δεi=Δεie+Δεip,where Δεi, Δεei, and Δεpi represents the strain and elastic strain and plastic strain increments, respectively.The incremental expression of the stress strain relationshipσ=Dε 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
(7)Δσi=αijΔεje,where αi represents the material parameter.According to Equation (5) and symmetry of the elastic matrix, Equation (7) is rewritten as
(8)S1Δε1e,Δε2e,Δε3e=α11Δε1e+α12Δε2e+α13Δε3eS2Δε1e,Δε2e,Δε3e=α12Δε1e+α22Δε2e+α23Δε3eS3Δε1e,Δε2e,Δε3e=α13Δε1e+α23Δε2e+α33Δε3e.The plastic strain part follows the relevant flow law [34–36]:
(9)Δεip=λ∂g∂σi,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:(10)Δεnps=λs∂gs∂σnΔεnpt=λt∂gt∂σn,where λs and λt are shear and tensile plastic factors, respectively, and gs and gt 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:
(11)Δε1ps=λs∂gs∂σ1=λsΔε2ps=λs∂gs∂σ2=0Δε3ps=λs∂gs∂σ3=−Nψλs.The plastic shear strain increment can also be expressed as follows:(12)dγps=λs∂gs∂σ1,∂gs∂σ2,∂gs∂σ3=λsΔε1ps,Δε2ps,Δε3ps.After substituting the partial derivative of Equation (2) into Equation (10), the plastic tensile strain in each principal direction can be obtained:
(13)Δε1pt=λt∂gt∂σ1=0Δε2pt=λt∂gt∂σ2=0Δε3pt=λt∂gt∂σ3=−λt.The plastic tensile strain increment can also be expressed as follows:(14)dγpt=λt∂gt∂σ1,∂gt∂σ2,∂gt∂σ3=λtΔε1pt,Δε2pt,Δε3pt.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(15)σiN=σi+Δσi,σiI=σi+SiΔεn,where σNi is the new stress component and σIi is the elastic guess.The stress plastic correction is expressed as follows:(16)σiN=σiI−λSi∂g∂σn.From Equations (8), (12), and (16) and replacing Δεe1, Δεe2, and Δεe3 with Δεps1, Δεps2, and Δεps3, stress of each principal direction after shear plastic correction is obtained as follows:
(17)σ1N=σ1I−λsα11−α13Nψσ2N=σ2I−λsα12−α23Nψσ3N=σ3I−λsα13−α33Nψ,λs=fsσ1I,σ3Iα11−α13Nψ−α13−α33NψNϕ.From Equations (8), (14), and (16) and replacing Δεe1, Δεe2, and Δεe3 with Δεpt1, Δεpt2, and Δεpt3, stress of each principal direction after tensile plastic correction is obtained as follows:
(18)σ1N=σ1I−λtα13σ2N=σ2I−λtα23σ3N=σ3I−λtα33,λt=σt−σ3Iα33.


## 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 t and the total strain increment for a timestep Δt, the stress state at time t+Δt 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. Figure2 shows the calculation cycle diagram. The constitutive model in the third step is obtained by the secondary development in this paper.Figure 2
FLAC3D calculation process.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 is30m×30m×30m, 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 X, Y, and Z 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.Table 2
Calculation parameters of constitutive model [40].
ParameterAN-MC modelTransversely isotropic modelLevel occurrenceDip direction angle dd (°)4545Dip angle dip (°)4545Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)1e91e9The Poisson ratio of the plane parallel to the bedding planeμ10.260.26The elastic modulus in the plane perpendicular to the bedding planeE2 (Pa)0.5e90.5e9The plane Poisson ratio in the plane perpendicular to the bedding planeμ20.360.36The shear modulus in the plane perpendicular to the bedding planeG2 (Pa)0.34e90.34e9Plastic parametersThe cohesionc (Pa)20e9—The frictionϕ (°)40—The tensile strengthσt (Pa)20e9—Figure 3
Comparison of stress, displacement contour plot, and plastic zone under different constitutive models.
(a)
Stress contour plot of the AN-MC model(b)
Displacement contour plot of the AN-MC model(c)
The plastic zone of the AN-MC model(d)
Stress contour plot of the transversely isotropic constitutive model(e)
Displacement contour of the plot if the transversely isotropic constitutive model(f)
The plastic zone of the transversely isotropic constitutive model


### 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 Table3. 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.Table 3
Calculation parameters of constitutive model [40].
ParameterAN-MC modelMohr-Coulomb modelLevel occurrenceDip direction angle dd (°)270—Dip angle dip (°)90—Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)1e91e9The Poisson ratio of the plane parallel to the bedding planeμ10.30The elastic modulus in the plane perpendicular to the bedding planeE2 (Pa)1e90.30The plane Poisson ratio in the plane perpendicular to the bedding planeμ20.30The shear modulus in the plane perpendicular to the bedding planeG2 (Pa)0.3846e9—Plastic parametersThe cohesionc (Pa)2e62e6The frictionϕ (°)4545The tensile strengthσt (Pa)10e610e6Figure 4
Comparison of stress, displacement contour plot, and plastic zone under different constitutive models.
(a)
Stress contour plot of the AN-MC model(b)
Displacement contour plot of the AN-MC model(c)
The plastic zone of the AN-MC model(d)
Stress contour plot of the MC constitutive model(e)
Displacement contour plot of the MC constitutive model(f)
The plastic zone of the MC constitutive model


## 3.1. Introduction of Numerical Calculation Model

The size of the numerical calculation model is30m×30m×30m, 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 X, Y, and Z 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.Table 2
Calculation parameters of constitutive model [40].
ParameterAN-MC modelTransversely isotropic modelLevel occurrenceDip direction angle dd (°)4545Dip angle dip (°)4545Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)1e91e9The Poisson ratio of the plane parallel to the bedding planeμ10.260.26The elastic modulus in the plane perpendicular to the bedding planeE2 (Pa)0.5e90.5e9The plane Poisson ratio in the plane perpendicular to the bedding planeμ20.360.36The shear modulus in the plane perpendicular to the bedding planeG2 (Pa)0.34e90.34e9Plastic parametersThe cohesionc (Pa)20e9—The frictionϕ (°)40—The tensile strengthσt (Pa)20e9—Figure 3
Comparison of stress, displacement contour plot, and plastic zone under different constitutive models.
(a)
Stress contour plot of the AN-MC model(b)
Displacement contour plot of the AN-MC model(c)
The plastic zone of the AN-MC model(d)
Stress contour plot of the transversely isotropic constitutive model(e)
Displacement contour of the plot if the transversely isotropic constitutive model(f)
The plastic zone of the transversely isotropic constitutive model


## 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 Table3. 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.Table 3
Calculation parameters of constitutive model [40].
ParameterAN-MC modelMohr-Coulomb modelLevel occurrenceDip direction angle dd (°)270—Dip angle dip (°)90—Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)1e91e9The Poisson ratio of the plane parallel to the bedding planeμ10.30The elastic modulus in the plane perpendicular to the bedding planeE2 (Pa)1e90.30The plane Poisson ratio in the plane perpendicular to the bedding planeμ20.30The shear modulus in the plane perpendicular to the bedding planeG2 (Pa)0.3846e9—Plastic parametersThe cohesionc (Pa)2e62e6The frictionϕ (°)4545The tensile strengthσt (Pa)10e610e6Figure 4
Comparison of stress, displacement contour plot, and plastic zone under different constitutive models.
(a)
Stress contour plot of the AN-MC model(b)
Displacement contour plot of the AN-MC model(c)
The plastic zone of the AN-MC model(d)
Stress contour plot of the MC constitutive model(e)
Displacement contour plot of the MC constitutive model(f)
The plastic zone of the MC constitutive model


## 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 theX, Y, and Z 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 X-Z plane are shown in Figure 5, where (a)–(e) are monitoring points for displacement.Table 4
Calculation parameters of the AN-MC constitutive model.
ParameterValueLevel occurrenceDip direction angle dd (°)90Dip angle dip (°)0Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)6.24e9The Poisson ratio of the plane parallel to the bedding planeμ10.3The elastic modulus of in the plane perpendicular to the bedding planeE2 (Pa)4.76e9The Poisson ratio of in the plane perpendicular to the bedding planeμ20.18The shear modulus of in the plane perpendicular to the bedding planeG2 (Pa)2.08e9Plastic parametersThe cohesionc (Pa)0.80e6The friction angleϕ (°)27The tensile strengthσt (Pa)0.1e6Figure 5
Geometric sizes and boundary conditions of the model.


### 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.Figure 6
Displacement comparison of monitoring points.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 (Figures7(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.Figure 7
Development process of the plastic zone in layered surrounding rock.


## 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 theX, Y, and Z 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 X-Z plane are shown in Figure 5, where (a)–(e) are monitoring points for displacement.Table 4
Calculation parameters of the AN-MC constitutive model.
ParameterValueLevel occurrenceDip direction angle dd (°)90Dip angle dip (°)0Elastic parametersThe elastic modulus of the plane parallel to the bedding planeE1 (Pa)6.24e9The Poisson ratio of the plane parallel to the bedding planeμ10.3The elastic modulus of in the plane perpendicular to the bedding planeE2 (Pa)4.76e9The Poisson ratio of in the plane perpendicular to the bedding planeμ20.18The shear modulus of in the plane perpendicular to the bedding planeG2 (Pa)2.08e9Plastic parametersThe cohesionc (Pa)0.80e6The friction angleϕ (°)27The tensile strengthσt (Pa)0.1e6Figure 5
Geometric sizes and boundary conditions of the model.


## 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.Figure 6
Displacement comparison of monitoring points.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 (Figures7(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.Figure 7
Development process of the plastic zone in layered surrounding rock.


## 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.Figure 6
Displacement comparison of monitoring points.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 (Figures7(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.Figure 7
Development process of the plastic zone in layered surrounding rock.


## 5. Discussion

(1)
In the secondary development process of the AN-MC model, the transformation between the transversely isotropic material parametersαi (α1, α2, and α3) of rock mass in Equation (7) and material parameters α1 and α2 of the MC yield criterion is complex. In this paper, the material parameters α1 and α2 in the incremental expression of Hooke’s law are transformed into material parameters e1 and e2 of the AN-MC model, and the incremental expression of transversely isotropic of layered rock mass is obtained. The key steps are as follows: ① α1=K+4G/3 and α2=K−2G/3 [30] in the incremental expressions of Hooke’s law are converted into e1=K′+4G′/3 and e2=K′−2G′/3, where K and G are the bulk and shear moduli of the isotropic material, and K′ and G′ are the bulk and shear moduli of the transversely isotropic material, respectively; that is, α1 and α2 are replaced with e1 and e2 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 E1 and that in the direction perpendicular to the bedding plane E2, Poisson’s ratio in the direction parallel to the bedding plane μ1 and that in the direction perpendicular to the bedding plane μ2, and the shear modulus in the direction perpendicular to the bedding plane G2—are transformed into two elastic parameters K′ and G′. The transformed program expression is shown in references [42], and their logical relationships are shown in Figures 8 and 9. E1 and μ1 are obtained by the uniaxial compression test of layered rock mass with a dip angle of 0°, E2 and μ2 are obtained by the uniaxial compression test of layered rock mass with dip angle of 90°, and G2 is calculated by the empirical formula (1): G2=1/E1+1/E2+2μ2/E2(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 Table4. The elastic moduli E1 and E2 are 6.24×109, and Poisson’s ratios μ1 and μ2 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 10Figure 8
Logic relationship of the bulk modulus conversion.Figure 9
Logic relationship of the shear modulus conversion.Figure 10
Development process of the plastic zone in isotropic surrounding rock.After the tunnel is excavated, the sidewalls and inverted arch of the rock mass first enter the plastic state (Figures10(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 Figures7 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



---

*Source: 3264675-2022-08-09.xml*