# Application of the Modified Mohr–Coulomb Yield Criterion in Seismic Numerical Simulation of Tunnels

Authors: Chuan-Yi Sui; Yu-Sheng Shen; Yu-Min Wen; Bo Gao
Journal: Shock and Vibration
 (2021)

---

## Abstract

To solve the classical problem that the Mohr–Coulomb yield criterion overestimates the tensile properties of geotechnical materials, a modified Mohr–Coulomb yield criterion that includes both maximum tensile stress theory and smooth processing was established herein. The modified Mohr–Coulomb constitutive model is developed using the user-defined material subroutine (UMAT) available in finite element software ABAQUS, and the modified Mohr–Coulomb yield criterion is applied to construct a numerical simulation of a shaking table model test. Compared with the measured data from the shaking table test, the accuracies of the classical Mohr–Coulomb yield criterion and the modified Mohr–Coulomb yield criterion are assessed. Compared to the shaking table test, the classical Mohr–Coulomb model has a relatively large average error (−6.98% in peak acceleration values, −8.47% in displacement values, −23.93% in axial forces), while the modified Mohr–Coulomb model has a smaller average error (+2.71% in peak accelerations value, +3.19% in displacements value, +7.56% in axial forces). The results of numerical simulation using the modified Mohr–Coulomb yield criterion are closer to the measured data.


---

## Body

## 1. Introduction

Since the Wenchuan earthquake in 2008 (Ms = 8.0), China has suffered from numerous serious earthquakes, like the Kaohsiung earthquake (Ms = 6.7, 2016) and the Yushu earthquake (Ms = 7.1, 2010) [1]. Most earthquakes pose potential threats to underground structures in mountainous areas. Taking the Longmen Mountain Fault Zone in the north-western margin of the Sichuan Basin as an example, an earthquake in this fault zone would impact a large number of traffic tunnels and underground structures. The subway stations and other underground structures were critically damaged [2, 3] in the Kobe earthquake in Japan, which were considered aseismic structures [4–7] in high-intensity earthquake area. The Lushan earthquake (Ms = 7.0, 2013), the Jiuzhaigou earthquake (Ms = 7.0, 2017), and the Changning earthquake (Ms = 6.0, 2019) show that the tunnels within the Sichuan Basin are threatened by high-intensity earthquakes. Therefore, studies of the seismic performance of underground structures should receive more attention.There are three main types of approaches for the analysis of seismic performance of underground structures: field investigations, model tests, and numerical simulations [8–13]. Numerical simulation methods have been widely used to analyse the dynamic response of tunnel structure in some bad conditions, but the computational ability of computers and the mechanical models needs to be studied further in the existing software. Most numerical analysis results should be verified using field tests or model tests, and improvements in the accuracy for numerical analysis would be beneficial.The classical Mohr–Coulomb yield criterion is the most widely used constitutive model in geotechnical mechanics. It focuses on the response of brittle materials to shear stress and normal stress. However, classical Mohr–Coulomb theory overestimates the tensile properties of geotechnical materials. Jia and Chen’s article shows that there are up to 18.7% differences in yield stress [14]. The dynamic response of the surrounding rocks in underground structures exhibits a clear tension-compression cycle during strong earthquake motion. When calculating the dynamic response of an underground structure, the classical Mohr–Coulomb yield criterion should be corrected.Zienkiewicz and Pande suggested using a hyperbola to approximate the two Mohr–Coulomb straight lines in the meridional plane. Deng et al. proposed the equal-area-cone principle, which uses a series of modified Drucker–Prager yield criteria to replace the Mohr–Coulomb yield criterion. But, the replaceable Drucker–Prager criteria are related to the stress states of rocks, which means the criteria of every element used may be changed during a dynamic simulation and need more computing resources [15–19]. Jia et al. combined the classical Mohr–Coulomb theory with the maximum tensile stress theory, which uses a parameter m to reduce the yield stress but leaves no definition of the parameter m [20]. Based on these previous studies, a constitutive model with the tensile strength replacing the parameter m in Jia’s model was developed. The modified model was applied to the dynamic numerical simulation of underground structures, and its advantages were compared with those of the classical Mohr–Coulomb model.


## 2. Modified Mohr–Coulomb Yield Criterion

For the constitutive model of geotechnical materials, the yield function is expressed in the form of stress invariants for convenience, as shown in the following:(1)σm=σ1+σ2+σ33,σ¯=12sx2+sy2+sz2+τxy2+τxz2+τyz2=J2,J3=sxsysz+2τxyτxzτyz−sxτyz2−syτxz2−szτxy2,where sx=σx−σm,sy=σy−σm,sz=σz−σm,andσm is the mean stress, σ¯ is the equivalent stress, and σmJ2andJ3 are the second and the third stress invariants, respectively.The Lode angel equation is as follows:(2)θ=13sin−1−332J3σ¯2,−30°≤θ≤30°.The classical Mohr–Coulomb model overestimates the tensile strength of materials. Jia et al. developed a modified Mohr–Coulomb yield criterion [20]; the yield function is given in (3). Kθ represents the yield curves with different Lode angles in the π plane.Kθ is calculated using (4):(3)F=σmsinϕ+σ¯2K2θ+mccosϕ2−ccosϕ=0,(4)Kθ=cosθ−13sinϕsinθ.Whenm=0, (3) is identical to the yield function of Mohr–Coulomb model (5). However, Jia and Chen did not provide a physical definition or recommend a value for the parameter m:(5)σmsinϕ+σ¯Kθ−ccosϕ=0.The yield function of the maximum tensile stress theory is given in (6), where Tc is the tensile strength:(6)23σ¯sinθ+120°+σm−Tc≥0.Zienkiewicz and Pande suggested replacing the classical Mohr–Coulomb yield surface function with a hyperbola [19]:(7)σm−d2a2−σ¯2b2=1.To fit Zienkiewicz and Pande’s hyperbola (7) using the maximum tensile stress yield function (6), the parameters in (7) can be set as a=ccotϕ−Tc, b=ccotϕ−Tcsinϕ/Kθ, and d=ccotϕ. Kθ is defined in (4). (7) can be simplified to(8)F=σmsinϕ+σ¯2K2θ+Tcsinϕ−ccosϕ2−ccosϕ=0.In the following, the yield function of the modified Mohr–Coulomb is set as (8).The classical Mohr–Coulomb yield surface is not smooth, and it has a singular vertex point. The corners and singular vertex point make the numerical calculations difficult. Shi and Yang suggest using a curve fit to approximate it [21]. Jia and Chen used a piecewise function to smooth the yield surface [20]. Since (8) is a special case of (3), Jia’s Kθ is introduced in (8). The result is obtained in(9)F=σmsinϕ+σ¯2K2θ+Tcsinϕ−ccosϕ2−ccosϕ=0.Kθ is calculated using (10) and (12), as follows:(10)Kθ=cosθ−13sinϕsinθ,whenθ≤θT,(11)Kθ=13cosθT3+tanθTtan3θT+13tan3θT−3tanθTsinϕ,(12)−13cosθTsinθT+13sinϕcosθTsin3θ,whenθ>θT.Based on Jia’s research,θT=27° is recommended. Limestone was used as an example. The cohesion is 60 kPa, the friction angle is 25°, and the tensile strength is 50 kPa. Figure 1 shows the curves of the constitutive model used in this paper and Jia’s model in the stress space. The x-axis is the mean stress (σm), and the y-axis is the equivalent stress (σ¯).Figure 1
Yield criterion function used in this article.As shown in Figure1, the model proposed by Jia is sensitive to the value of parameter m and no recommended value of m is given in Jia’s paper [20]. The uniaxial tensile strength is substituted for the parameter m. The constitutive model in this paper is similar to the classical Mohr–Coulomb constitutive model in the compression stage, and the control criterion of the tensile stress is approximated in the tensile section. The results are consistent with the experimental results under uniaxial tension.The plastic potential function is consistent with the yield function:(13)G=σmsinϕ+σ¯2K2θ+Tcsinϕ−ccosϕ2.The parameters in (13) are identical as those in (9).


## 3. Constitutive Model Used in ABAQUS Software

ABAQUSTM can employ FORTRANTM programs to develop user-defined materials (UMATs). The calculation process of the general user subroutine [20] is shown in Figure 2.Figure 2
Flow chart of the UMAT used in ABAQUS software.


### 3.1. Calculation of Plastic Parameters

When the stress exceeds the yield surface, it needs to adjust the stress and return to the updated yield surface (Figure3). The method used is called the constitutive integration algorithm [20].Figure 3
Diagram of the backward Euler integration algorithm.As shown in Figure3, the stress state at point A (σA) is in the yield plane, while that at point B (σB) is outside of the yield plane:(14)σB=σA+DeΔε=σA+Δσe.In (14), De is the elastic matrix, Δσe is the elastic stress increment, and σB is the probing stress. To bring the probing stress σB into the yield plane, the change in the stress Δσ is calculated in (15), where b=∂G/∂σ:(15)Δσ=DeΔε−ΔλDeb=Δσe−ΔλDeb.The stress at point C is calculated in(16)σC=σB−ΔλDeb.The function of the backward Euler algorithm is to gradually iterate point C onto the yield surface. The yield function F is expanded using a first-order Taylor expansion at point B:(17)F=FB+∂F∂σTΔσ+∂F∂εpl¯Δεpl¯=FB−ΔλaBTDebB−ΔλAB′=0,Δλ=FBaBTDebB+AB′.The backward Euler algorithm is used to pull the stress back to the yield surface, as follows:(18)σC=σB−ΔλDebC.The initial estimate ofσC is based on the stress at point B:(19)σC=σB−ΔλDebB.The vectorr is defined as the difference between the current stress and the backward Euler stress:(20)r=σ−σB−ΔλDebC=σ−σB+ΔλDebC.Whenr=0, the stress is back on the yield plane. The Taylor expansion of r with fixed σB is(21)rn=r0+σ˙+λ˙Deb+ΔλDe∂b∂σσ˙.Whenrn=0,(22)σ˙=−I+ΔλDe∂b∂σ−1r0+λ˙Deb=−Q−1r0−λ˙Q−1Deb.I is the identity matrix.The Taylor expansion of the yield function F at point C is as follows:(23)FCn=FC0+∂F∂σTσ˙+∂F∂εpl¯εpl¯˙=FC0+aCTσ˙+AC′λ˙=0.Substituting (22) into (23) gives(24)λ˙=FC−aCTQ−1r0aCTQ−1Deb+AC′.The equivalent plastic strain rate′/εpl is as follows:(25)εpl¯˙=λ˙Bσ.For the von Mises yield criteria,Bσ=1, whereas for the yield criterion of rock and soil, Bσ≠1.


### 3.2. Solution of the Stiffness Matrix of the Same TangentDct

According to (19), if the subscript of the stress state C of the current configuration in the following formula is ignored and the subscript B represents the elastic probing stress, (19) can be expressed as(26)σ=σB−ΔλDeb.The differential form of (26) is(27)σ˙=Deε˙−λ˙Deb−ΔDe∂b∂σσ˙.The simplified form of (27) is as follows:(28)σ˙=I+ΔλDe∂b∂σ−1Deε˙−λ˙b=Q−1Deε˙−λ˙b=Rε˙−λ˙b,where R=Q−1De.To bring the stressσ back to the yield plane, F˙=0 is needed. According to the yield function F, the consistency conditions are as follows:(29)F˙=aTσ˙−A′λ˙=0.Substituting (28) into (29) gives(30)λ˙=aTRε˙aTRb+A′.Substituting (30) into (27) gives(31)σ˙=R−RbaTRTaTRb+A′ε˙=Dctε˙.The stiffness matrix of the same tangent isDct.


### 3.3. Newton–Raphson Iterative Algorithm for Plastic Parameters under thenth Loading Step

In the fully implicit backward Euler method, the plastic strain and the equivalent plastic strain increment are calculated at the end of thenth incremental step. The integral algorithm can be expressed as follows:(32)εn+1=εn+Δε,εpn+1=εpn+Δλn+1bn+1,ε¯pln+1=ε¯pln+Δλn+1Bn+1,σn+1=De:εn+1−εpn+1,Fn+1=Fn+1σn+1,ε¯pln+1.(32) includes nonlinear equations about εn+1,εpn+1,ε¯pln+1 and εn,εpn,ε¯pln, and the strain increment Δε=Δtε˙ is given at tn. These nonlinear equations are solved using the Newton–Raphson iterative method, and the stresses are updated as follows.Step 1.
The initial values of the plastic strain and the equivalent plastic strain are the convergence values at the end of the last loading step. The incremental value of the plastic parameter is set to 0 to calculate the elastic stress as follows:(33)k=0:εp0n=εpn−1,ε¯pl0n=ε¯pln−1,Δλ0n=0,σ0n=De:εn−εp0n.Step 2.
Check the yield condition and convergence of the k-th iteration.(34)Fkn=Fσkn,ε¯plkn.
IfFkn<tolerance, the calculation is convergent. Otherwise, proceed to Step 3.Step 3.
The increment of the plastic parameter is calculated, and the increments of the stress and the equivalent plastic strain are obtained.(35)δλk=Fk−akTQ−1r0akTQ−1Deb+Ak',δσk=−I+ΔλDe∂b∂σ−1r0+λ˙Deb,δε¯plk=δλkB.Step 4.
The plastic strain, stress, and plastic multiplier are updated.(36)Δλk+1=Δλk+δλk,ε¯plk+1n=ε¯plkn+δε¯plk,εplk+1n=εplkn+De−1:δσk,σk+1n=σkn+De:εn−εplk+1n.
Setk=k+1, and return to Step 3.
When using the implicit backward Euler algorithm, the first derivative of the yield function, the first derivative of the potential function, and the second derivative of the potential function are needed.


### 3.4. Derivations of Yield Function and Potential Function

#### 3.4.1. The First Derivative of Yield Function

The vectora is the derivative of the yield function with stress:(37)a=∂F∂σ=C1′∂σm∂σ+C2′∂σ¯∂σ+C3′∂J3∂σ.where C1′=∂F/∂σm;C2′=∂F/∂σ¯−tan3θ/σ¯∂F/∂θ;andC3′=−3/2cos3θσ¯3∂F/∂θ. The other parameters are calculated as follows:(38)∂J3∂σ=sysz−τyz2sxsz−τxz2sxsy−τxy22τyzτxz−szτxy2τxyτyz−syτxz2τxzτxy−sxτyz+σ¯23111000,∂σm∂σ=131,1,1,0,0,0T,∂σ¯∂σ=12σ¯sx,sy,sz,2τxy,2τxz,2τyzT,C1′=sinϕ,C2′=α′K−tan3θdKdθ,C3′=α′−32cos3θσ¯2dKdθ,α′=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(39)dKdθ=−3Bcos3θ,θ>θT,sinθ−13sinϕcosθ,θ≤θT.


#### 3.4.2. The First Derivative of Potential Function

The vectorb is the derivative of the potential function with stress:(40)b=∂G∂σ=C1∂σm∂σ+C2∂σ¯∂σ+C3∂J3∂σ.where C1=∂G/∂σm;C2=∂G/∂σ¯−tan3θ/σ¯∂G/∂θ;andC3=−3/2cos3θσ¯3∂G/∂θ. The other parameters are calculated as follows:(41)C1=sinϕ,C2=αC2mc,C3=αC3mc.The calculation ofdK/dθis shown in (39). The calculation of α is shown in(42)α=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(43)C2mc=K−tan3θdKdθ,C3mc=−32cos3θσ¯2dKdθ.


#### 3.4.3. Quadratic Derivation of Plastic Potential Function

(44)∂b∂σ=∂2G∂σ2=∂C2∂σ∂σ¯∂σ+C2∂2σ¯∂σ2+∂C3∂σ∂J3∂σ+C2∂2J3∂σ2,∂C2∂σ=α∂C2mc∂σ+C2mc∂α∂σ,∂C3∂σ=α∂C3mc∂σ+C3mc∂α∂σ,where(45)∂C2mc∂σ=∂θ∂σ∂K∂θ+∂2K∂θ2tan3θ−3∂K∂θsec23θ∂C3mc∂σ=32cos3θσ¯2∂θ∂σ∂2K∂θ2−3∂K∂θtan3θ+2σ¯∂K∂θ∂σ¯∂σ∂α∂σ=1−α2σ¯2K2+Tcsinϕ−ccosϕ2∂σ¯∂σK+σ¯∂K∂θ∂θ∂σ∂θ∂σ=−32σ¯3cos3θ∂J3∂σ−3J3σ¯∂σ¯∂σ,∂2J3∂σ2=13sx−sy−sz2szsy−sx−szsym2sy2sxsz−sx−sy2τxy2τxy−4τxy−6sz2τxz−4τxz2τxz6τyz−6sy−4τyz2τyz2τyz6τxz6τxy−6sz,∂2σ¯∂σ2=1σ13−sx24σ¯2−16−sxsy4σ¯213−sy24σ¯2sym−16−sxsz4σ¯2−16−sysz4σ¯213−sz24σ¯2−τxysx2σ¯2−τxysy2σ¯2−τxysz2σ¯21−τxy2σ¯2−τxzsx2σ¯2−τxzsy2σ¯2−τxzsz2σ¯2−τxzτxyσ¯21−τxz2σ¯2−τyzsx2σ¯2−τxzsy2σ¯2−τyzsz2σ¯2−τyzsxyσ¯2−τyzτxzσ¯21−τyz2σ¯2.


## 3.1. Calculation of Plastic Parameters

When the stress exceeds the yield surface, it needs to adjust the stress and return to the updated yield surface (Figure3). The method used is called the constitutive integration algorithm [20].Figure 3
Diagram of the backward Euler integration algorithm.As shown in Figure3, the stress state at point A (σA) is in the yield plane, while that at point B (σB) is outside of the yield plane:(14)σB=σA+DeΔε=σA+Δσe.In (14), De is the elastic matrix, Δσe is the elastic stress increment, and σB is the probing stress. To bring the probing stress σB into the yield plane, the change in the stress Δσ is calculated in (15), where b=∂G/∂σ:(15)Δσ=DeΔε−ΔλDeb=Δσe−ΔλDeb.The stress at point C is calculated in(16)σC=σB−ΔλDeb.The function of the backward Euler algorithm is to gradually iterate point C onto the yield surface. The yield function F is expanded using a first-order Taylor expansion at point B:(17)F=FB+∂F∂σTΔσ+∂F∂εpl¯Δεpl¯=FB−ΔλaBTDebB−ΔλAB′=0,Δλ=FBaBTDebB+AB′.The backward Euler algorithm is used to pull the stress back to the yield surface, as follows:(18)σC=σB−ΔλDebC.The initial estimate ofσC is based on the stress at point B:(19)σC=σB−ΔλDebB.The vectorr is defined as the difference between the current stress and the backward Euler stress:(20)r=σ−σB−ΔλDebC=σ−σB+ΔλDebC.Whenr=0, the stress is back on the yield plane. The Taylor expansion of r with fixed σB is(21)rn=r0+σ˙+λ˙Deb+ΔλDe∂b∂σσ˙.Whenrn=0,(22)σ˙=−I+ΔλDe∂b∂σ−1r0+λ˙Deb=−Q−1r0−λ˙Q−1Deb.I is the identity matrix.The Taylor expansion of the yield function F at point C is as follows:(23)FCn=FC0+∂F∂σTσ˙+∂F∂εpl¯εpl¯˙=FC0+aCTσ˙+AC′λ˙=0.Substituting (22) into (23) gives(24)λ˙=FC−aCTQ−1r0aCTQ−1Deb+AC′.The equivalent plastic strain rate′/εpl is as follows:(25)εpl¯˙=λ˙Bσ.For the von Mises yield criteria,Bσ=1, whereas for the yield criterion of rock and soil, Bσ≠1.


## 3.2. Solution of the Stiffness Matrix of the Same TangentDct

According to (19), if the subscript of the stress state C of the current configuration in the following formula is ignored and the subscript B represents the elastic probing stress, (19) can be expressed as(26)σ=σB−ΔλDeb.The differential form of (26) is(27)σ˙=Deε˙−λ˙Deb−ΔDe∂b∂σσ˙.The simplified form of (27) is as follows:(28)σ˙=I+ΔλDe∂b∂σ−1Deε˙−λ˙b=Q−1Deε˙−λ˙b=Rε˙−λ˙b,where R=Q−1De.To bring the stressσ back to the yield plane, F˙=0 is needed. According to the yield function F, the consistency conditions are as follows:(29)F˙=aTσ˙−A′λ˙=0.Substituting (28) into (29) gives(30)λ˙=aTRε˙aTRb+A′.Substituting (30) into (27) gives(31)σ˙=R−RbaTRTaTRb+A′ε˙=Dctε˙.The stiffness matrix of the same tangent isDct.


## 3.3. Newton–Raphson Iterative Algorithm for Plastic Parameters under thenth Loading Step

In the fully implicit backward Euler method, the plastic strain and the equivalent plastic strain increment are calculated at the end of thenth incremental step. The integral algorithm can be expressed as follows:(32)εn+1=εn+Δε,εpn+1=εpn+Δλn+1bn+1,ε¯pln+1=ε¯pln+Δλn+1Bn+1,σn+1=De:εn+1−εpn+1,Fn+1=Fn+1σn+1,ε¯pln+1.(32) includes nonlinear equations about εn+1,εpn+1,ε¯pln+1 and εn,εpn,ε¯pln, and the strain increment Δε=Δtε˙ is given at tn. These nonlinear equations are solved using the Newton–Raphson iterative method, and the stresses are updated as follows.Step 1.
The initial values of the plastic strain and the equivalent plastic strain are the convergence values at the end of the last loading step. The incremental value of the plastic parameter is set to 0 to calculate the elastic stress as follows:(33)k=0:εp0n=εpn−1,ε¯pl0n=ε¯pln−1,Δλ0n=0,σ0n=De:εn−εp0n.Step 2.
Check the yield condition and convergence of the k-th iteration.(34)Fkn=Fσkn,ε¯plkn.
IfFkn<tolerance, the calculation is convergent. Otherwise, proceed to Step 3.Step 3.
The increment of the plastic parameter is calculated, and the increments of the stress and the equivalent plastic strain are obtained.(35)δλk=Fk−akTQ−1r0akTQ−1Deb+Ak',δσk=−I+ΔλDe∂b∂σ−1r0+λ˙Deb,δε¯plk=δλkB.Step 4.
The plastic strain, stress, and plastic multiplier are updated.(36)Δλk+1=Δλk+δλk,ε¯plk+1n=ε¯plkn+δε¯plk,εplk+1n=εplkn+De−1:δσk,σk+1n=σkn+De:εn−εplk+1n.
Setk=k+1, and return to Step 3.
When using the implicit backward Euler algorithm, the first derivative of the yield function, the first derivative of the potential function, and the second derivative of the potential function are needed.


## 3.4. Derivations of Yield Function and Potential Function

### 3.4.1. The First Derivative of Yield Function

The vectora is the derivative of the yield function with stress:(37)a=∂F∂σ=C1′∂σm∂σ+C2′∂σ¯∂σ+C3′∂J3∂σ.where C1′=∂F/∂σm;C2′=∂F/∂σ¯−tan3θ/σ¯∂F/∂θ;andC3′=−3/2cos3θσ¯3∂F/∂θ. The other parameters are calculated as follows:(38)∂J3∂σ=sysz−τyz2sxsz−τxz2sxsy−τxy22τyzτxz−szτxy2τxyτyz−syτxz2τxzτxy−sxτyz+σ¯23111000,∂σm∂σ=131,1,1,0,0,0T,∂σ¯∂σ=12σ¯sx,sy,sz,2τxy,2τxz,2τyzT,C1′=sinϕ,C2′=α′K−tan3θdKdθ,C3′=α′−32cos3θσ¯2dKdθ,α′=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(39)dKdθ=−3Bcos3θ,θ>θT,sinθ−13sinϕcosθ,θ≤θT.


### 3.4.2. The First Derivative of Potential Function

The vectorb is the derivative of the potential function with stress:(40)b=∂G∂σ=C1∂σm∂σ+C2∂σ¯∂σ+C3∂J3∂σ.where C1=∂G/∂σm;C2=∂G/∂σ¯−tan3θ/σ¯∂G/∂θ;andC3=−3/2cos3θσ¯3∂G/∂θ. The other parameters are calculated as follows:(41)C1=sinϕ,C2=αC2mc,C3=αC3mc.The calculation ofdK/dθis shown in (39). The calculation of α is shown in(42)α=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(43)C2mc=K−tan3θdKdθ,C3mc=−32cos3θσ¯2dKdθ.


### 3.4.3. Quadratic Derivation of Plastic Potential Function

(44)∂b∂σ=∂2G∂σ2=∂C2∂σ∂σ¯∂σ+C2∂2σ¯∂σ2+∂C3∂σ∂J3∂σ+C2∂2J3∂σ2,∂C2∂σ=α∂C2mc∂σ+C2mc∂α∂σ,∂C3∂σ=α∂C3mc∂σ+C3mc∂α∂σ,where(45)∂C2mc∂σ=∂θ∂σ∂K∂θ+∂2K∂θ2tan3θ−3∂K∂θsec23θ∂C3mc∂σ=32cos3θσ¯2∂θ∂σ∂2K∂θ2−3∂K∂θtan3θ+2σ¯∂K∂θ∂σ¯∂σ∂α∂σ=1−α2σ¯2K2+Tcsinϕ−ccosϕ2∂σ¯∂σK+σ¯∂K∂θ∂θ∂σ∂θ∂σ=−32σ¯3cos3θ∂J3∂σ−3J3σ¯∂σ¯∂σ,∂2J3∂σ2=13sx−sy−sz2szsy−sx−szsym2sy2sxsz−sx−sy2τxy2τxy−4τxy−6sz2τxz−4τxz2τxz6τyz−6sy−4τyz2τyz2τyz6τxz6τxy−6sz,∂2σ¯∂σ2=1σ13−sx24σ¯2−16−sxsy4σ¯213−sy24σ¯2sym−16−sxsz4σ¯2−16−sysz4σ¯213−sz24σ¯2−τxysx2σ¯2−τxysy2σ¯2−τxysz2σ¯21−τxy2σ¯2−τxzsx2σ¯2−τxzsy2σ¯2−τxzsz2σ¯2−τxzτxyσ¯21−τxz2σ¯2−τyzsx2σ¯2−τxzsy2σ¯2−τyzsz2σ¯2−τyzsxyσ¯2−τyzτxzσ¯21−τyz2σ¯2.


## 3.4.1. The First Derivative of Yield Function

The vectora is the derivative of the yield function with stress:(37)a=∂F∂σ=C1′∂σm∂σ+C2′∂σ¯∂σ+C3′∂J3∂σ.where C1′=∂F/∂σm;C2′=∂F/∂σ¯−tan3θ/σ¯∂F/∂θ;andC3′=−3/2cos3θσ¯3∂F/∂θ. The other parameters are calculated as follows:(38)∂J3∂σ=sysz−τyz2sxsz−τxz2sxsy−τxy22τyzτxz−szτxy2τxyτyz−syτxz2τxzτxy−sxτyz+σ¯23111000,∂σm∂σ=131,1,1,0,0,0T,∂σ¯∂σ=12σ¯sx,sy,sz,2τxy,2τxz,2τyzT,C1′=sinϕ,C2′=α′K−tan3θdKdθ,C3′=α′−32cos3θσ¯2dKdθ,α′=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(39)dKdθ=−3Bcos3θ,θ>θT,sinθ−13sinϕcosθ,θ≤θT.


## 3.4.2. The First Derivative of Potential Function

The vectorb is the derivative of the potential function with stress:(40)b=∂G∂σ=C1∂σm∂σ+C2∂σ¯∂σ+C3∂J3∂σ.where C1=∂G/∂σm;C2=∂G/∂σ¯−tan3θ/σ¯∂G/∂θ;andC3=−3/2cos3θσ¯3∂G/∂θ. The other parameters are calculated as follows:(41)C1=sinϕ,C2=αC2mc,C3=αC3mc.The calculation ofdK/dθis shown in (39). The calculation of α is shown in(42)α=σ¯Kσ¯2K2+Tcsinϕ−ccosϕ2,(43)C2mc=K−tan3θdKdθ,C3mc=−32cos3θσ¯2dKdθ.


## 3.4.3. Quadratic Derivation of Plastic Potential Function

(44)∂b∂σ=∂2G∂σ2=∂C2∂σ∂σ¯∂σ+C2∂2σ¯∂σ2+∂C3∂σ∂J3∂σ+C2∂2J3∂σ2,∂C2∂σ=α∂C2mc∂σ+C2mc∂α∂σ,∂C3∂σ=α∂C3mc∂σ+C3mc∂α∂σ,where(45)∂C2mc∂σ=∂θ∂σ∂K∂θ+∂2K∂θ2tan3θ−3∂K∂θsec23θ∂C3mc∂σ=32cos3θσ¯2∂θ∂σ∂2K∂θ2−3∂K∂θtan3θ+2σ¯∂K∂θ∂σ¯∂σ∂α∂σ=1−α2σ¯2K2+Tcsinϕ−ccosϕ2∂σ¯∂σK+σ¯∂K∂θ∂θ∂σ∂θ∂σ=−32σ¯3cos3θ∂J3∂σ−3J3σ¯∂σ¯∂σ,∂2J3∂σ2=13sx−sy−sz2szsy−sx−szsym2sy2sxsz−sx−sy2τxy2τxy−4τxy−6sz2τxz−4τxz2τxz6τyz−6sy−4τyz2τyz2τyz6τxz6τxy−6sz,∂2σ¯∂σ2=1σ13−sx24σ¯2−16−sxsy4σ¯213−sy24σ¯2sym−16−sxsz4σ¯2−16−sysz4σ¯213−sz24σ¯2−τxysx2σ¯2−τxysy2σ¯2−τxysz2σ¯21−τxy2σ¯2−τxzsx2σ¯2−τxzsy2σ¯2−τxzsz2σ¯2−τxzτxyσ¯21−τxz2σ¯2−τyzsx2σ¯2−τxzsy2σ¯2−τyzsz2σ¯2−τyzsxyσ¯2−τyzτxzσ¯21−τyz2σ¯2.


## 4. Verification of the Modified Model

Through comparison to the actual data from the shaking table test, the relative errors of both the improved Mohr–Coulomb model and the classical Mohr–Coulomb model were calculated. The calculation accuracy of the improved constitutive model was analysed.


### 4.1. Introduction to the Shaking Table Test

The primary objective of the shaking table test was to test how the vibration-absorptive material protects a tunnel in a fault zone [21]. The overall layout of the shaking table is shown in Figures 4 and 5. The detailed accelerometer monitoring plan is presented in Table 1. The rock surrounding the tunnel between the two fault zones is placed on a movable platform with 10 cm-thick foam pieces on both the left and right sides. Both the left and right sides of the model contain a 25 cm-thick foam board to absorb the reflected earthquake waves. The accelerometers used in the test were LC0113 M piezoelectric accelerometers [22] (Figures 6 and 7). Strain gauges used in the test were BX120-50AA resin-based strain gauges.Figure 4
Layout of the shaking table test (all dimensions are in cm). (a) Longitudinal plane of the test and layout of accelerometers. (b) Cross section of layout for strain gauges. (c) Layout of strain monitoring sections.
(a)(b)(c)Figure 5
Photograph of the shaking table test.Table 1
Accelerometers and test objects.
AccelerometersTest objectsA1HBoundary accelerationA2HRock acceleration near the damping layerA3HTunnel acceleration at the fault with a damping layerA4HTunnel acceleration between the faultsA5HRock acceleration at the fault without a damping layerA6HRock acceleration between the faultsA7HRock acceleration near the damping layerA8HSurface accelerationFigure 6
Accelerometers installed on the shaking table.Figure 7
Accelerometers buried in the simulation material.The prototype of the model is the Longxi tunnel which was highly damaged in the Wenchuan earthquake. The input wave of the shaking table is part of the east and west vector components of the Wolong seismic waves which is widely used in simulating tunnel dynamic response [21]. The original seismic wave lasts for 180 s, and the peak acceleration value is 0.98 g. The input wave occurs at 20 s to 110 s. The time history curve of the input wave is shown in Figure 8. The frequency-domain analysis shown in Figure 9 suggests that the predominant frequency is about 2.3 Hz and the majority of the frequency domain is below 10 Hz.Figure 8
Time history curve of the input Wuolong wave for the shaking table test.Figure 9
Frequency-domain analysis of acceleration of the Wuolong wave.In view of the limitations of the shaking table and considering the similarity criterion, the length similarity ratio was taken as 1 : 25, the density similarity ratio as 1 : 1.3, and the stress similarity ratio as 1 : 32.5. A detailed list of similarity relations is presented in Table2.Table 2
Similarity ratios of the shaking table test.
Physical quantityDimensionSimilarity relationSimilarity ratioLengthLlr1 : 25Elastic modulusFL−2Er1 : 32.5DensityFTL−4ρr1 : 1.3StressFL−2σr=Er1 : 32.5AccelerationLT−2ar=Erρr−1lr−11 : 1Poisson's ratio—μr=11 : 1Friction angle—θr=11 : 1The prototype model for the rock material is classified as Grade V according to the Chinese railway tunnel design specifications (TB10003-2005) [23]. The prototype for the surrounding rock is strongly weathered rock. The prototype lining material is C30 concrete. The prototype for the fault zone is gravel. The material similarity relationships are detailed in Tables 3–5.Table 3
Physical parameters of the surrounding rock.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material13001.070.253.0625Prototype material169034.80.2599.525Similarity ratio1 : 1.31 : 32.51 : 11 : 32.51 : 1Table 4
Physical parameters of the linings.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material1780571.420.260.7150.8Prototype material2314185710.2197350.8Similarity ratio1 : 1.31 : 32.51 : 11 : 32.51 : 1Table 5
Physical parameters of the faults.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material15000.710.25—35Prototype material1950230.25—35Similarity ratio1 : 1.31 : 32.51 : 1—1 : 1


### 4.2. Introduction to the Numerical Simulation

The numerical simulation was conducted using ABAQUSTM software. The overall model layout is shown in Figure 10. The model is based on the shaking table test, and the modeling is 1 : 1. A viscoelastic boundary was introduced to absorb the reflected earthquake waves. The physical parameters are shown in Table 6. The prototype surrounding rock is a strongly weathered rock mass, and the tensile strength of the model material is 3.0 kPa [24]. The stiffness damping coefficient of the foam board and the damping material is 0.2.Figure 10
Computational model of the shaking table test.Table 6
Physical parameters of the materials.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Surrounding rock13001.070.253.0625Linings1780571.420.260.7150.8Fault15000.710.25-35Q235 steel75002100000.25--Damping layer150700.05--The input earthquake wave shown in Figure11 is the first 26 s acceleration of the shaking table test. The maximum horizontal acceleration of the input wave is 0.97 g, and the sampling is 200 Hz.Figure 11
Time history of the input wave of the numerical simulation.The maximum frequency of the seismic wave is 100 Hz; the seismic wave velocity in the surrounding rock is 1.35 km/s; and the maximum length of the model units is 10 cm, which is less than one-tenth of the maximum frequency corresponding to the seismic wave wavelength. The boundary filtering effect is negligible. Rayleigh damping was introduced in the dynamic analysis, and the damping matrix isC=ηM+β:(46)η=2ωiωjωi+ωjζ,β=2ωi+ωjζ,where ζ is the ith- or jth-order damping ratio of the natural vibration frequency. ζi=ζj and ωi and ωj are the ith- and jth- order natural vibration frequencies, respectively. The damping ratio was obtained from the test and is ζ=0.15. When running the modal analysis in ABAQUSTM software, the 1st to 6th orders of the natural vibration frequency are acquired (Table 7).Table 7
First six orders of the natural vibration frequency of the rock material.
1st order2nd order3rd order4th order5th order6th orderNatural vibratory circular frequencies (rad/s)8.9812.5818.70519.30423.18523.58The Rayleigh damping parameters of the surrounding rock were calculated as follows:η=1.951andβ=0.009.


### 4.3. Acceleration Analysis of the Modified Mohr–Coulomb Model

The peak accelerations are introduced as the measurements [25–29]. By comparing the acceleration watch points of the two Mohr–Coulomb models with the results of the shaking table test, the peak acceleration similarity advantage of the modified Mohr–Coulomb yield criterion was assessed.As shown in Figures12–19, both the classical and modified Mohr–Coulomb yield criteria have acceleration time histories similar to those of the shaking table test, which indicates that both the classical and modified Mohr–Coulomb yield criteria can relatively truly represent the transmitted earthquake wave observed in the shaking table test.Figure 12
Acceleration time history at point A4H.Figure 13
: Frequency-domain analysis of point A4H.Figure 14
Acceleration time history at point A5H.Figure 15
Frequency-domain analysis of point A5H.Figure 16
Acceleration time history at point A6H.Figure 17
Frequency-domain analysis of point A6H.Figure 18
Acceleration time history at point A7H.Figure 19
Frequency-domain analysis of point A7H.Based on Figures13, 15, 17, and 19, the modified Mohr–Coulomb yield criterion has more similar waveforms in the low frequency band (<10 Hz) than the classical Mohr–Coulomb yield criterion. The amplitudes of the modified Mohr–Coulomb criterion are closer to the test data than those of the classical ones. The frequencies of each peak in the modified model are closer than those of the classical model. In the middle frequency band (10–15 Hz), both the classical and modified Mohr–Coulomb yield criteria have lower accelerations than the shaking table test. However, further research is still required. In the high frequency band (>15 Hz), the two Mohr–Coulomb yield criteria and the shaking table test have little energy (<0.005).As reported in Table8, the modified Mohr–Coulomb yield criterion has more similar peak accelerations than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the acceleration of the modified Mohr–Coulomb yield criterion at A4H is 2.8% smaller, while that of the classical Mohr–Coulomb yield criterion is 8.8% smaller. At A5H, the acceleration of the modified Mohr–Coulomb yield criterion is 1.6% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 5.5% smaller.Table 8
Maximum accelerations at the monitoring points for the different yield criterion models (unit: g).
PointShaking table testClassical MC modelModified MC modelA4H1.0370.9451.017A5H0.9800.9441.016A6H1.0140.9451.016A7H1.2341.1281.344Figure20 shows that the average peak acceleration error of the modified Mohr–Coulomb model is +2.71% and the maximum relative error is +8.91%. The average error of the classical Mohr–Coulomb model is −6.98%, and the maximum relative error is −8.87%. The modified model has a lower average error than the classical model, and both models have similar maximum relative errors. In aseismic tunnel analysis, a larger acceleration peak value leads to a more conservative aseismic design.Figure 20
Histogram of the maximum accelerations at all monitoring points for the different yield criterion models.


### 4.4. Displacement Analysis of the Modified Mohr–Coulomb Model

The peak horizontal displacements are introduced as the measurements [30–32]. Displacement can be obtained by acceleration double integral for time. By comparing the displacement watch points of the two Mohr–Coulomb models with the results of the shaking table test, the maximum horizontal displacement similarity advantage of the modified Mohr–Coulomb yield criterion was assessed.As shown in Figures21–24, the displacements of the shaking table test and classical and modified Mohr–Coulomb yield criteria are in minus, which indicate that the models were moving towards the negative direction of the sensors. To simplify the following analysis, all displacements are taken as absolute values.Figure 21
Displacement time history at point A4H.Figure 22
Displacement time history at point A5H.Figure 23
Displacement time history at point A6H.Figure 24
Displacement time history at point A7H.As reported in Table9, the modified Mohr–Coulomb yield criterion has more similar maximum displacements than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the displacement of the modified Mohr–Coulomb yield criterion at A4H is 6.41% smaller, while that of the classical Mohr–Coulomb yield criterion is 10.68% smaller. At A5H, the displacement of the modified Mohr–Coulomb yield criterion is 5.37% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 12.68% smaller.Table 9
Maximum displacements at the monitoring points for the different yield criterion models (unit: m).
PointShaking table testClassical MC modelModified MC modelA4H0.2340.2090.219A5H0.2050.1790.216A6H0.3230.3070.324A7H0.2870.2710.300Figure25 shows that the average displacement error of the modified Mohr–Coulomb model is +3.79% and the maximum relative error is -6.41%. The average displacement error of the classical Mohr–Coulomb model is −8.47%, and the maximum relative error is −12.68%. The modified model has both lower average and relative errors than the classical model.Figure 25
Histogram of the maximum displacements at all monitoring points for the different yield criterion models.


### 4.5. Axial Force Analysis of the Modified Mohr–Coulomb Model

The peak axial forces are introduced as the measurements [33]. Axial forces can be obtained by calculating the stress or strain at monitoring points. In numerical simulations, axial forces can be obtained from the elemental solutions. The axial forces of the shaking table test could be calculated from(47)N=112Ebhε1+ε2,where E is the elastic modulus of the lining, b is the width of the element, h is the thickness of the lining, and ε1 and ε2 are the strains of inner and outer surfaces at the watch point.In this article,E=571.42MPa, b=0.064m, and h=0.01m.As shown in Figures26 and 27, there are both positive and negative numbers of the axial forces of the shaking table test and classical and modified Mohr–Coulomb yield criteria, which indicate that the models were tensioned (positive numbers) and compressed (negative numbers). To simplify the following analysis, all axial forces are taken as absolute values.Figure 26
Axial force time history of section1 spandrel.Figure 27
: Axial force time history of section2 arch springing.In Table10, the modified Mohr–Coulomb yield criterion has more similar maximum axial forces than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the maximum axial force of the modified Mohr–Coulomb yield criterion at the crown of section 1 is 6.03% higher, while that of the classical Mohr–Coulomb yield criterion is 16.96% smaller. At the invert of section 2, the maximum axial force of the modified Mohr–Coulomb yield criterion is 3.10% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 18.83% smaller.Table 10
Maximum axial forces at the monitoring points for the different yield criterion models (unit: N).
SectionPointShaking table testClassical MC modelModified MC model1Crown7.966.618.441Spandrel15.6111.3517.821Sidewall8.316.138.621Arch springing28.1319.6930.661Invert9.727.4810.502Crown5.764.325.882Spandrel6.604.826.732Sidewall4.103.164.512Arch springing7.356.398.642Invert1.291.061.33Figures28 and 29 show that the average axial forces error of the modified Mohr–Coulomb model is +7.56% and the maximum relative error is +17.55%. The average axial force error of the classical Mohr–Coulomb model is −23.93%, and the maximum relative error is −27.0%. The modified model has both lower average and relative errors than the classical model. In conclusion, the modified Mohr–Coulomb model exhibits a better performance in simulating the dynamic response of rocks during high-intensity horizontal earthquake.Figure 28
Histogram of the maximum axial forces at section1 for the different yield criterion models.Figure 29
Histogram of the maximum axial forces at section2 for the different yield criterion models.


## 4.1. Introduction to the Shaking Table Test

The primary objective of the shaking table test was to test how the vibration-absorptive material protects a tunnel in a fault zone [21]. The overall layout of the shaking table is shown in Figures 4 and 5. The detailed accelerometer monitoring plan is presented in Table 1. The rock surrounding the tunnel between the two fault zones is placed on a movable platform with 10 cm-thick foam pieces on both the left and right sides. Both the left and right sides of the model contain a 25 cm-thick foam board to absorb the reflected earthquake waves. The accelerometers used in the test were LC0113 M piezoelectric accelerometers [22] (Figures 6 and 7). Strain gauges used in the test were BX120-50AA resin-based strain gauges.Figure 4
Layout of the shaking table test (all dimensions are in cm). (a) Longitudinal plane of the test and layout of accelerometers. (b) Cross section of layout for strain gauges. (c) Layout of strain monitoring sections.
(a)(b)(c)Figure 5
Photograph of the shaking table test.Table 1
Accelerometers and test objects.
AccelerometersTest objectsA1HBoundary accelerationA2HRock acceleration near the damping layerA3HTunnel acceleration at the fault with a damping layerA4HTunnel acceleration between the faultsA5HRock acceleration at the fault without a damping layerA6HRock acceleration between the faultsA7HRock acceleration near the damping layerA8HSurface accelerationFigure 6
Accelerometers installed on the shaking table.Figure 7
Accelerometers buried in the simulation material.The prototype of the model is the Longxi tunnel which was highly damaged in the Wenchuan earthquake. The input wave of the shaking table is part of the east and west vector components of the Wolong seismic waves which is widely used in simulating tunnel dynamic response [21]. The original seismic wave lasts for 180 s, and the peak acceleration value is 0.98 g. The input wave occurs at 20 s to 110 s. The time history curve of the input wave is shown in Figure 8. The frequency-domain analysis shown in Figure 9 suggests that the predominant frequency is about 2.3 Hz and the majority of the frequency domain is below 10 Hz.Figure 8
Time history curve of the input Wuolong wave for the shaking table test.Figure 9
Frequency-domain analysis of acceleration of the Wuolong wave.In view of the limitations of the shaking table and considering the similarity criterion, the length similarity ratio was taken as 1 : 25, the density similarity ratio as 1 : 1.3, and the stress similarity ratio as 1 : 32.5. A detailed list of similarity relations is presented in Table2.Table 2
Similarity ratios of the shaking table test.
Physical quantityDimensionSimilarity relationSimilarity ratioLengthLlr1 : 25Elastic modulusFL−2Er1 : 32.5DensityFTL−4ρr1 : 1.3StressFL−2σr=Er1 : 32.5AccelerationLT−2ar=Erρr−1lr−11 : 1Poisson's ratio—μr=11 : 1Friction angle—θr=11 : 1The prototype model for the rock material is classified as Grade V according to the Chinese railway tunnel design specifications (TB10003-2005) [23]. The prototype for the surrounding rock is strongly weathered rock. The prototype lining material is C30 concrete. The prototype for the fault zone is gravel. The material similarity relationships are detailed in Tables 3–5.Table 3
Physical parameters of the surrounding rock.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material13001.070.253.0625Prototype material169034.80.2599.525Similarity ratio1 : 1.31 : 32.51 : 11 : 32.51 : 1Table 4
Physical parameters of the linings.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material1780571.420.260.7150.8Prototype material2314185710.2197350.8Similarity ratio1 : 1.31 : 32.51 : 11 : 32.51 : 1Table 5
Physical parameters of the faults.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Model material15000.710.25—35Prototype material1950230.25—35Similarity ratio1 : 1.31 : 32.51 : 1—1 : 1


## 4.2. Introduction to the Numerical Simulation

The numerical simulation was conducted using ABAQUSTM software. The overall model layout is shown in Figure 10. The model is based on the shaking table test, and the modeling is 1 : 1. A viscoelastic boundary was introduced to absorb the reflected earthquake waves. The physical parameters are shown in Table 6. The prototype surrounding rock is a strongly weathered rock mass, and the tensile strength of the model material is 3.0 kPa [24]. The stiffness damping coefficient of the foam board and the damping material is 0.2.Figure 10
Computational model of the shaking table test.Table 6
Physical parameters of the materials.
Density (kg/m3)Elastic modulus (MPa)Poisson's ratioCohesion (kPa)Friction angle (°)Surrounding rock13001.070.253.0625Linings1780571.420.260.7150.8Fault15000.710.25-35Q235 steel75002100000.25--Damping layer150700.05--The input earthquake wave shown in Figure11 is the first 26 s acceleration of the shaking table test. The maximum horizontal acceleration of the input wave is 0.97 g, and the sampling is 200 Hz.Figure 11
Time history of the input wave of the numerical simulation.The maximum frequency of the seismic wave is 100 Hz; the seismic wave velocity in the surrounding rock is 1.35 km/s; and the maximum length of the model units is 10 cm, which is less than one-tenth of the maximum frequency corresponding to the seismic wave wavelength. The boundary filtering effect is negligible. Rayleigh damping was introduced in the dynamic analysis, and the damping matrix isC=ηM+β:(46)η=2ωiωjωi+ωjζ,β=2ωi+ωjζ,where ζ is the ith- or jth-order damping ratio of the natural vibration frequency. ζi=ζj and ωi and ωj are the ith- and jth- order natural vibration frequencies, respectively. The damping ratio was obtained from the test and is ζ=0.15. When running the modal analysis in ABAQUSTM software, the 1st to 6th orders of the natural vibration frequency are acquired (Table 7).Table 7
First six orders of the natural vibration frequency of the rock material.
1st order2nd order3rd order4th order5th order6th orderNatural vibratory circular frequencies (rad/s)8.9812.5818.70519.30423.18523.58The Rayleigh damping parameters of the surrounding rock were calculated as follows:η=1.951andβ=0.009.


## 4.3. Acceleration Analysis of the Modified Mohr–Coulomb Model

The peak accelerations are introduced as the measurements [25–29]. By comparing the acceleration watch points of the two Mohr–Coulomb models with the results of the shaking table test, the peak acceleration similarity advantage of the modified Mohr–Coulomb yield criterion was assessed.As shown in Figures12–19, both the classical and modified Mohr–Coulomb yield criteria have acceleration time histories similar to those of the shaking table test, which indicates that both the classical and modified Mohr–Coulomb yield criteria can relatively truly represent the transmitted earthquake wave observed in the shaking table test.Figure 12
Acceleration time history at point A4H.Figure 13
: Frequency-domain analysis of point A4H.Figure 14
Acceleration time history at point A5H.Figure 15
Frequency-domain analysis of point A5H.Figure 16
Acceleration time history at point A6H.Figure 17
Frequency-domain analysis of point A6H.Figure 18
Acceleration time history at point A7H.Figure 19
Frequency-domain analysis of point A7H.Based on Figures13, 15, 17, and 19, the modified Mohr–Coulomb yield criterion has more similar waveforms in the low frequency band (<10 Hz) than the classical Mohr–Coulomb yield criterion. The amplitudes of the modified Mohr–Coulomb criterion are closer to the test data than those of the classical ones. The frequencies of each peak in the modified model are closer than those of the classical model. In the middle frequency band (10–15 Hz), both the classical and modified Mohr–Coulomb yield criteria have lower accelerations than the shaking table test. However, further research is still required. In the high frequency band (>15 Hz), the two Mohr–Coulomb yield criteria and the shaking table test have little energy (<0.005).As reported in Table8, the modified Mohr–Coulomb yield criterion has more similar peak accelerations than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the acceleration of the modified Mohr–Coulomb yield criterion at A4H is 2.8% smaller, while that of the classical Mohr–Coulomb yield criterion is 8.8% smaller. At A5H, the acceleration of the modified Mohr–Coulomb yield criterion is 1.6% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 5.5% smaller.Table 8
Maximum accelerations at the monitoring points for the different yield criterion models (unit: g).
PointShaking table testClassical MC modelModified MC modelA4H1.0370.9451.017A5H0.9800.9441.016A6H1.0140.9451.016A7H1.2341.1281.344Figure20 shows that the average peak acceleration error of the modified Mohr–Coulomb model is +2.71% and the maximum relative error is +8.91%. The average error of the classical Mohr–Coulomb model is −6.98%, and the maximum relative error is −8.87%. The modified model has a lower average error than the classical model, and both models have similar maximum relative errors. In aseismic tunnel analysis, a larger acceleration peak value leads to a more conservative aseismic design.Figure 20
Histogram of the maximum accelerations at all monitoring points for the different yield criterion models.


## 4.4. Displacement Analysis of the Modified Mohr–Coulomb Model

The peak horizontal displacements are introduced as the measurements [30–32]. Displacement can be obtained by acceleration double integral for time. By comparing the displacement watch points of the two Mohr–Coulomb models with the results of the shaking table test, the maximum horizontal displacement similarity advantage of the modified Mohr–Coulomb yield criterion was assessed.As shown in Figures21–24, the displacements of the shaking table test and classical and modified Mohr–Coulomb yield criteria are in minus, which indicate that the models were moving towards the negative direction of the sensors. To simplify the following analysis, all displacements are taken as absolute values.Figure 21
Displacement time history at point A4H.Figure 22
Displacement time history at point A5H.Figure 23
Displacement time history at point A6H.Figure 24
Displacement time history at point A7H.As reported in Table9, the modified Mohr–Coulomb yield criterion has more similar maximum displacements than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the displacement of the modified Mohr–Coulomb yield criterion at A4H is 6.41% smaller, while that of the classical Mohr–Coulomb yield criterion is 10.68% smaller. At A5H, the displacement of the modified Mohr–Coulomb yield criterion is 5.37% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 12.68% smaller.Table 9
Maximum displacements at the monitoring points for the different yield criterion models (unit: m).
PointShaking table testClassical MC modelModified MC modelA4H0.2340.2090.219A5H0.2050.1790.216A6H0.3230.3070.324A7H0.2870.2710.300Figure25 shows that the average displacement error of the modified Mohr–Coulomb model is +3.79% and the maximum relative error is -6.41%. The average displacement error of the classical Mohr–Coulomb model is −8.47%, and the maximum relative error is −12.68%. The modified model has both lower average and relative errors than the classical model.Figure 25
Histogram of the maximum displacements at all monitoring points for the different yield criterion models.


## 4.5. Axial Force Analysis of the Modified Mohr–Coulomb Model

The peak axial forces are introduced as the measurements [33]. Axial forces can be obtained by calculating the stress or strain at monitoring points. In numerical simulations, axial forces can be obtained from the elemental solutions. The axial forces of the shaking table test could be calculated from(47)N=112Ebhε1+ε2,where E is the elastic modulus of the lining, b is the width of the element, h is the thickness of the lining, and ε1 and ε2 are the strains of inner and outer surfaces at the watch point.In this article,E=571.42MPa, b=0.064m, and h=0.01m.As shown in Figures26 and 27, there are both positive and negative numbers of the axial forces of the shaking table test and classical and modified Mohr–Coulomb yield criteria, which indicate that the models were tensioned (positive numbers) and compressed (negative numbers). To simplify the following analysis, all axial forces are taken as absolute values.Figure 26
Axial force time history of section1 spandrel.Figure 27
: Axial force time history of section2 arch springing.In Table10, the modified Mohr–Coulomb yield criterion has more similar maximum axial forces than the classical Mohr–Coulomb yield criterion. Compared to the shaking table test, the maximum axial force of the modified Mohr–Coulomb yield criterion at the crown of section 1 is 6.03% higher, while that of the classical Mohr–Coulomb yield criterion is 16.96% smaller. At the invert of section 2, the maximum axial force of the modified Mohr–Coulomb yield criterion is 3.10% larger than that of the shaking table test, while that of the classical Mohr–Coulomb yield criterion is 18.83% smaller.Table 10
Maximum axial forces at the monitoring points for the different yield criterion models (unit: N).
SectionPointShaking table testClassical MC modelModified MC model1Crown7.966.618.441Spandrel15.6111.3517.821Sidewall8.316.138.621Arch springing28.1319.6930.661Invert9.727.4810.502Crown5.764.325.882Spandrel6.604.826.732Sidewall4.103.164.512Arch springing7.356.398.642Invert1.291.061.33Figures28 and 29 show that the average axial forces error of the modified Mohr–Coulomb model is +7.56% and the maximum relative error is +17.55%. The average axial force error of the classical Mohr–Coulomb model is −23.93%, and the maximum relative error is −27.0%. The modified model has both lower average and relative errors than the classical model. In conclusion, the modified Mohr–Coulomb model exhibits a better performance in simulating the dynamic response of rocks during high-intensity horizontal earthquake.Figure 28
Histogram of the maximum axial forces at section1 for the different yield criterion models.Figure 29
Histogram of the maximum axial forces at section2 for the different yield criterion models.


## 5. Conclusions

The modified Mohr–Coulomb model was successfully developed and used to analyse the dynamic response of tunnels during high-intensity earthquakes:(1)
Compared to the shaking table test, the classical Mohr–Coulomb model has a relatively large average error (−6.98% in peak accelerations, -8.47% in maximum displacements, and -23.93% in axial forces), while the modified Mohr–Coulomb model has a smaller average error (+2.71% in peak accelerations, +3.79% in maximum displacements, and +7.56% in axial forces).(2)
The absolute values of the maximum acceleration relative errors of the two models are similar. Compared to the shaking table test, the classical Mohr–Coulomb model has a +8.91% maximum relative error and the modified Mohr–Coulomb model has a −8.87% maximum relative error. The classical Mohr–Coulomb model has larger relative errors in displacements (-12.68%) and axial forces (-27.0%) than those of the modified Mohr–Coulomb model (-6.41% in displacements and +17.55% in axial forces).(3)
The new modified model has better waveforms in the low frequency band (0–10 Hz) than the classical Mohr–Coulomb model. In the higher frequency band (>15 Hz), the difference between the two models is negligible.(4)
The classical Mohr–Coulomb model overshoots the tensile properties of the geotechnical materials, which leads to higher yield stresses in simulating the surrounding rocks of the shaking table test. The greater the dynamic stress the rock bears, the less the dynamic stress the tunnel lining bears during the earthquake. Less dynamic stress of the lining leads to lower peak accelerations, lower maximum displacements, and fewer axial forces, which match the numerical simulation results. The modified Mohr–Coulomb model reduces the yield stresses of surrounding rocks and brings the simulation results closer to the shaking table test data.The new modified Mohr–Coulomb yield criterion has overall closer results in simulating the accelerations, displacements, and axial forces of the tunnel linings in the horizontal shaking table test than the classical Mohr–Coulomb yield criterion. Since most seismic wave energy is in the low frequency band (Figure9), the modified Mohr–Coulomb yield criterion works better than the classical one. Moreover, the modified criterion presented in this article improves Jia’s equation by replacing the undefined parameter m with the tensile strength.



---

*Source: 9968935-2021-11-25.xml*