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.
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 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:where is the mean stress, is the equivalent stress, and σm are the second and the third stress invariants, respectively.
The Lode angel equation is as follows:
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). represents the yield curves with different Lode angles in the π plane. is calculated using (4):
When , (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 :
The yield function of the maximum tensile stress theory is given in (6), where 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 (7) can be set as , , and . is defined in (4). (7) can be simplified to
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 is introduced in (8). The result is obtained in
is calculated using (10) and (12), as follows:
Based on Jia’s research, 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 (), 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]. 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:
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.
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 (Figure 3). The method used is called the constitutive integration algorithm [20].
As shown in Figure 3, the stress state at point A () is in the yield plane, while that at point B () is outside of the yield plane:
In (14), is the elastic matrix, is the elastic stress increment, and is the probing stress. To bring the probing stress into the yield plane, the change in the stress is calculated in (15), where :
The stress at point C is calculated in
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:
The backward Euler algorithm is used to pull the stress back to the yield surface, as follows:
The initial estimate of is based on the stress at point B:
The vector is defined as the difference between the current stress and the backward Euler stress:
When, the stress is back on the yield plane. The Taylor expansion of with fixed is
When ,
is the identity matrix.
The Taylor expansion of the yield function F at point C is as follows:
Substituting (22) into (23) gives
The equivalent plastic strain rate is as follows:
For the von Mises yield criteria, , whereas for the yield criterion of rock and soil, .
3.2. Solution of the Stiffness Matrix of the Same Tangent
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
The differential form of (26) is
The simplified form of (27) is as follows:where .
To bring the stress back to the yield plane, is needed. According to the yield function F, the consistency conditions are as follows:
Substituting (28) into (29) gives
Substituting (30) into (27) gives
The stiffness matrix of the same tangent is .
3.3. Newton–Raphson Iterative Algorithm for Plastic Parameters under the nth 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. The integral algorithm can be expressed as follows:
(32) includes nonlinear equations about and and the strain increment is given at . 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:
Step 2. Check the yield condition and convergence of the k-th iteration.If , 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.
Step 4. The plastic strain, stress, and plastic multiplier are updated.Set , 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 vector is the derivative of the yield function with stress:where . The other parameters are calculated as follows:
3.4.2. The First Derivative of Potential Function
The vector b is the derivative of the potential function with stress:where . The other parameters are calculated as follows:
The calculation of is shown in (39). The calculation of is shown in
3.4.3. Quadratic Derivation of Plastic Potential Function
where
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.
(a)
(b)
(c)
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.
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.
The 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.
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.
The input earthquake wave shown in Figure 11 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.
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 is :where is the ith- or jth-order damping ratio of the natural vibration frequency. and and are the ith- and jth- order natural vibration frequencies, respectively. The damping ratio was obtained from the test and is . When running the modal analysis in ABAQUSTM software, the 1st to 6th orders of the natural vibration frequency are acquired (Table 7).
The Rayleigh damping parameters of the surrounding rock were calculated as follows: .
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 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. 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 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 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.
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%. 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.
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 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 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.
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%. 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.
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 fromwhere is the elastic modulus of the lining, is the width of the element, is the thickness of the lining, and and are the strains of inner and outer surfaces at the watch point.
In this article, , , and .
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%. 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.
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 (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 with the tensile strength.
Data Availability
The data used to support the findings of this study are included within the article. The data are also available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors thank LetPub (https://www.letpub.com) for its linguistic assistance and scientific consultation during the preparation of this manuscript. This study was funded by the National Natural Science Foundation of China (Nos. 51778540, 5167850, and 52108361). This research was also supported by the Miaozi Engineering Key Project from Science and Technology Department of Sichuan Province (No. 2020JDRC0078) and the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (China) Open Fund (No. SKLGP2020K018).