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

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. )e 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). )e results of numerical simulation using the modified Mohr–Coulomb yield criterion are closer to the measured data.


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. e subway stations and other underground structures were critically damaged [2,3] in the Kobe earthquake in Japan, which were considered aseismic structures [4][5][6][7] in highintensity earthquake area. e 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 highintensity earthquakes. erefore, studies of the seismic performance of underground structures should receive more attention.
ere are three main types of approaches for the analysis of seismic performance of underground structures: field investigations, model tests, and numerical simulations [8][9][10][11][12][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. e 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]. e 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][16][17][18][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. e 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.

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: where and σ m is the mean stress, σ is the equivalent stress, and σm J 2 and J 3 are the second and the third stress invariants, respectively. e Lode angel equation is as follows: e 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): When m � 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: e yield function of the maximum tensile stress theory is given in (6), where T c is the tensile strength: Zienkiewicz and Pande suggested replacing the classical Mohr-Coulomb yield surface function with a hyperbola [19]: To fit Zienkiewicz and Pande's hyperbola (7) using the maximum tensile stress yield function (6), the parameters in In the following, the yield function of the modified Mohr-Coulomb is set as (8).
e classical Mohr-Coulomb yield surface is not smooth, and it has a singular vertex point. e 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). e result is obtained in K(θ) is calculated using (10) and (12), as follows:

2
Shock and Vibration Based on Jia's research, θ T � 27°is recommended. Limestone was used as an example. e 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. e x-axis is the mean stress (σ m ), and the y-axis is the equivalent stress (σ).
As shown in Figure 1, 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]. e uniaxial tensile strength is substituted for the parameter m. e 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. e results are consistent with the experimental results under uniaxial tension. e plastic potential function is consistent with the yield function: e parameters in (13) are identical as those in (9).

Constitutive Model Used in ABAQUS Software
ABAQUS TM can employ FORTRAN TM programs to develop user-defined materials (UMATs). e calculation process of the general user subroutine [20] is shown in Figure 2.

Calculation of Plastic Parameters.
When the stress exceeds the yield surface, it needs to adjust the stress and return to the updated yield surface ( Figure 3). e method used is called the constitutive integration algorithm [20]. As shown in Figure 3, the stress state at point A (σ A ) is in the yield plane, while that at point B (σ B ) is outside of the yield plane: In (14), D e 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 � zG/zσ: e stress at point C is calculated in e function of the backward Euler algorithm is to gradually iterate point C onto the yield surface. e yield function F is expanded using a first-order Taylor expansion at point B: e backward Euler algorithm is used to pull the stress back to the yield surface, as follows: e initial estimate of σ C is based on the stress at point B: e vector r is defined as the difference between the current stress and the backward Euler stress: When r � 0, the stress is back on the yield plane. e Taylor expansion of r with fixed σ B is When r n � 0, I is the identity matrix. e Taylor expansion of the yield function F at point C is as follows: Substituting (22) into (23) gives e equivalent plastic strain rate ′ /ε pl is as follows: For the von Mises yield criteria, B(σ) � 1, whereas for the yield criterion of rock and soil, B(σ) ≠ 1.

Solution of the Stiffness Matrix of the Same Tangent D ct .
According to (19), if the subscript of the stress state C of the current configuration in the following formula is ignored Shock and Vibration and the subscript B represents the elastic probing stress, (19) can be expressed as e differential form of (26) is e simplified form of (27) is as follows: where 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: Substituting (28) into (29) gives Substituting (30) into (27) gives e stiffness matrix of the same tangent is D ct .

Newton-Raphson Iterative Algorithm for Plastic Parameters under the n th Loading
Step. In the fully implicit backward Euler method, the plastic strain and the equivalent plastic strain increment are calculated at the end of the nth incremental step. e integral algorithm can be expressed as follows: , and the strain increment Δε � Δt_ ε is given at t n . ese nonlinear equations are solved using the Newton-Raphson iterative method, and the stresses are updated as follows.
e initial values of the plastic strain and the equivalent plastic strain are the convergence values at the  end of the last loading step. e incremental value of the plastic parameter is set to 0 to calculate the elastic stress as follows: Step 2. Check the yield condition and convergence of the k-th iteration.
If F (n) (k) < tolerance, the calculation is convergent. Otherwise, proceed to Step 3.
e increment of the plastic parameter is calculated, and the increments of the stress and the equivalent plastic strain are obtained.
Step 4. e plastic strain, stress, and plastic multiplier are updated.
Set k � 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.

e First Derivative of Yield Function.
e vector a is the derivative of the yield function with stress: where e other parameters are calculated as follows:

e First Derivative of Potential Function.
e vector b is the derivative of the potential function with stress: e other parameters are calculated as follows: e calculation of dK/dθis shown in (39). e calculation of α is shown in Shock and Vibration

Verification of the Modified Model
rough 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. e calculation accuracy of the improved constitutive model was analysed. Table Test. e primary objective of the shaking table test was to test how the vibration-absorptive material protects a tunnel in a fault zone [21]. e overall layout of the shaking table is shown in Figures 4 and 5. e detailed accelerometer monitoring plan is presented in Table 1. e 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 cmthick foam board to absorb the reflected earthquake waves. e 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.

Introduction to the Shaking
e prototype of the model is the Longxi tunnel which was highly damaged in the Wenchuan earthquake. e 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]. e original seismic wave lasts for 180 s, and the peak acceleration value is 0.98 g. e input wave occurs at 20 s to 110 s. e time history curve of the input wave is shown in Figure 8. e 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.
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 Table 2. e prototype model for the rock material is classified as Grade V according to the Chinese railway tunnel design specifications (TB10003-2005) [23]. e prototype for the surrounding rock is strongly weathered rock. e prototype lining material is C30 concrete. e prototype for the fault zone is gravel. e material similarity relationships are detailed in Tables 3-5.

Introduction to the Numerical Simulation.
e numerical simulation was conducted using ABAQUS TM software. e overall model layout is shown in Figure 10. e 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. e physical parameters are shown in Table 6. e prototype surrounding rock is a strongly weathered rock mass, and the tensile strength of the model material is 3.0 kPa [24]. e stiffness damping coefficient of the foam board and the damping material is 0.2. e input earthquake wave shown in Figure 11 is the first 26 s acceleration of the shaking table test. e maximum horizontal acceleration of the input wave is 0.97 g, and the sampling is 200 Hz. e 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. e boundary filtering effect is negligible. Rayleigh damping was introduced in the dynamic analysis, and the damping matrix is C � ηM + β: 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. e damping ratio was obtained from the test and is ζ � 0.15. When running the modal analysis in ABAQUS TM software, the 1 st to 6 th orders of the natural vibration frequency are acquired (Table 7).
e Rayleigh damping parameters of the surrounding rock were calculated as follows: η � 1.951 and β � 0.009.

Acceleration Analysis of the Modified Mohr-Coulomb
Model.
e peak accelerations are introduced as the measurements [25][26][27][28][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 Figures 12-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.
Based on Figures 13, 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. e amplitudes of the modified Mohr-Coulomb criterion are closer to the test data than those of the classical ones. e 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 Table 8, 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 S9 S10 S11 S12 S13 S14 S15 S16 crown invert (b)    Figure 20 shows that the average peak acceleration error of the modified Mohr-Coulomb model is +2.71% and the maximum relative error is +8.91%. e average error of the classical Mohr-Coulomb model is −6.98%, and the maximum relative error is −8.87%. e 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.

Displacement Analysis of the Modified Mohr-Coulomb
Model.
e peak horizontal displacements are introduced as the measurements [30][31][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 Figures 21-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.
As reported in Table 9, the modified Mohr-Coulomb yield criterion has more similar maximum displacements      Figure 25 shows that the average displacement error of the modified Mohr-Coulomb model is +3.79% and the maximum relative error is -6.41%. e average displacement error of the classical Mohr-Coulomb model is −8.47%, and        Amplitude (g) Figure 13: : Frequency-domain analysis of point A4H.     the maximum relative error is −12.68%. e modified model has both lower average and relative errors than the classical model.

Axial Force Analysis of the Modified Mohr-Coulomb
Model. e 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. e axial forces of the shaking table test could be calculated from 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.42 MPa, b � 0.064 m, and h � 0.01 m.
As shown in Figures 26 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.
In Table 10, 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. Figures 28 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%. e average axial force error of the classical Mohr-Coulomb model is −23.93%, and the maximum relative error is −27.0%. e modified model has both lower average and relative errors than the classical   Amplitude (g) Figure 19: Frequency-domain analysis of point A7H.       model. In conclusion, the modified Mohr-Coulomb model exhibits a better performance in simulating the dynamic response of rocks during high-intensity horizontal earthquake.

Conclusions
e modified Mohr-Coulomb model was successfully developed and used to analyse the dynamic response of tunnels during high-intensity earthquakes: (    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) e 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) e 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. e 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. e modified Mohr-Coulomb model reduces the yield stresses of surrounding rocks and brings the simulation results closer to the shaking table test data.
e 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 (Figure 9), 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.

Data Availability
e data used to support the findings of this study are included within the article. e data are also available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.