A Three-Dimensional Elastoplastic Constitutive Model for Geomaterials

Abstract

The Mohr-Coulomb (M-C) failure criterion has been a popular choice for geotechnical analysis because of its simplicity and ease of use. The fact that the M-C criterion disregards the intermediate principal stress’s impact is a significant drawback. As a result, the M-C criterion is only applied to materials under biaxial stress. This paper presents a three-dimensional version of the M-C criterion. The proposed criterion, called the Generalized Mohr-Coulomb (GMC) criterion, considers the intermediate principal stress’s effect, in addition to inheriting the original M-C criterion’s benefits. We obtained the conditions that the strength parameters must satisfy when the GMC criterion fulfills the π plane’s convexity. The GMC criterion can better describe geotechnical materials’ strengths under general stress conditions. Based on an implicit algorithm, the user material subroutine (UMAT) of the three-dimensional GMC model was developed in ABAQUS using the Fortran programming language. The established elastoplastic model’s validity and the program’s accuracy were examined using numerical simulation. Finally, a numerical simulation of a three-dimensional tunnel excavation under various working conditions was performed. The calculation results from the GMC model are precise and have some engineering-related practical significance.

1. Introduction

Including soils’ behaviors into a mathematical framework is required to calculate the bearing capacity. A failure criterion is typically used to define a material’s strength capabilities under particular circumstances [1]. Here, the simple Mohr-Coulomb material model is by far the most frequently used [2].

Geotechnical materials belong to the collection of scattered particles, which are widely distributed and are closely related to engineering construction. A crucial mechanical characteristic of construction materials is their shear strength [3]. Hence, the shear strength must be easy to obtain from the failure criteria. Additionally, the M-C criterion directly provides the shear strength. Since the limit surface of the M-C criterion is composed of piecewise linear equations, the corner points are formed at the intersection of the two equations [4]. If the associated flow rule is utilized in deformation analysis, the presence of corner points will result in the plastic flow direction not being unique [5]. To overcome this problem, many scholars have proposed a variety of smooth shape functions (ridge functions) to replace the M-C failure criterion [6,7]. Although ridge functions eliminate corner points, their mathematical expressions are cumbersome and lack clear physical meaning, making them not easy to apply. According to Zheng et al. [8,9], however, the corner problem has been solved.

The geotechnical materials’ failure criteria are generally divided into two categories, from the perspective of whether they have physical meaning. One category is hypothesis criteria, based on physical concepts, and the other is statistical abstract criteria, based on the findings of countless experiments. The former criteria include the M-C, the SMP [10], and the twin shear criteria [11]. The latter criteria include the Lade-Duncan [12], Hoek-Brown [13], and Yoshida criteria [14]. Due to its strong applicability, obvious physical meaning, and ease in calculating strength parameters using direct shear or traditional triaxial testing, the M-C failure criterion is frequently employed in geotechnical engineering. However, it ignores the intermediate principal stress’s influence.

Some researchers use elliptic, hyperbolic, and spatially moving plane deviatoric functions to extend the M-C criterion to three dimensions while taking the effect of ${\sigma }_2$ into account [15,16,17]. However, the existing 3D M-C failure criterion has no obvious advantage in either the mechanical mechanism or mathematical form. For example, the Mogi-Coulomb criterion [18] has significantly improved the true triaxial data’s fitting accuracy but does not meet the convexity condition in the deviatoric plane. Singh et al. [19] proposed the nonlinear M-C criterion formula but found that this criterion’s applicability was limited.

A popular calculation technique for solving a variety of engineering problems is the finite element method. It can be combined with various numerical approaches [20], which is just one of its numerous benefits. The key part of numerically calculating geotechnical material plasticity is the stress update. Since the stress update is performed multiple times in each loading step, the process must be fast and accurate to ensure an effective numerical solution. In general, schemes based on explicit and implicit procedures constitute the two categories of constitutive equation-integration algorithms [21]. The explicit integration algorithm is simple and easy to implement, but there may be cases where the predicted stress deviates from the current yield surface, thus resulting in inaccurate calculation results. Implicit algorithms have attracted increasing attention because of the progress. The return mapping algorithm was first put forth by Krieg and Krieg [22], and it has since been rapidly refined and used. Most implicit algorithms fall within the return mapping scheme category. However, the implicit return mapping algorithm is difficult to develop and is actually locally convergent. Based on Koiter, De Borst [23] presented a non-smooth yield surface plastic equation integration approach that is only useful when two yield surfaces cross. Adhikary [24] proposed a convex programming-based technique for multi-surface plastic calculation that, in order to calculate the stress at an integral point, needs the Hessian matrices for each effective yield function. Haji Aghajanpour et al. and Zheng et al. [8] proposed the Gauss-Seidel projection-contraction (GSPC) algorithm, based on the Gauss-Seidel idea and projection-contraction algorithm. The GSPC algorithm converges for any strain increment size, the computational efficiency is higher than the return mapping algorithm, and it does not need to calculate the potential function’s Hessian matrix.

ABAQUS has been extensively utilized in civil engineering computation because it is a powerful generic nonlinear finite element analysis program that can be used to analyze and compute several difficult linear and nonlinear engineering issues. Through the user subroutine UMAT, Liu et al. [25] integrated the extended Drucker Prager model into ABAQUS and applied it to tunnel construction. Based on the return mapping method, Zhang et al. [26] incorporated the Tresca yield criterion, plastic flow, and hardening rules into a UMAT in the finite element analysis program ABAQUS. Kossa and Szabó [27] proposed a novel precise stress integration strategy based on the combined linear hardening von Mises elastoplastic model, which is implemented in ABAQUS through its UMAT material interface. Some researchers [28,29,30,31,32] have also used the UMAT material interface in ABAQUS to numerically implement the existing model.

The aforementioned research shows that it is challenging for the current failure criteria to strike a balance between a straightforward mathematical form and a distinct physical meaning. In this paper, a generalized M-C (GMC) failure criterion that considers the intermediate principal stress’s impact is established. The suggested criterion simply regards its parameters of cohesion c and friction factor f as functions of the coefficient b of the intermediate principal stress, written as c(b) and f(b), respectively. In this way, the effect of σ₂ can be naturally reflected in the failure criterion. A true triaxial failure criterion can be easily established simply using a set of traditional triaxial compression tests and a set of triaxial extension tests. The GMC failure criterion’s predictions for several geomaterials were found to be congruent with experimental findings when they were compared to data from real triaxial testing. Based on the GSPC algorithm in the principal stress space, the UMAT subroutine of the three-dimensional GMC model is developed in ABAQUS using the Fortran programming language. The numerical simulation of a uniaxial compression test and the simulation of an underground chamber under hydrostatic pressure were performed to assess the established elastoplastic model’s viability and the program’s correctness. Finally, the GMC model was more accurate in determining the extent of the plastic zone and the surrounding rock stress by simulating tunnel excavation, which has a certain practical significance.

2. Geometric Representation of Stress Space

The three-dimensional space defined by principal stresses ${\sigma }_1,{\sigma }_2,{\sigma }_3$ is usually called Haigh-Westergaard stress space. Figure 1 provides a traditional failure surface for an isotope pressure-sensitive material. In geotechnical materials, the mechanical compression is usually considered to be positive. The stress point P in the Haigh-Westergaard stress space is shown in Figure 1. A line with the same angle as all the coordinate axes is called a hydrostatic line. The π plane is also called a deviatoric plane that is perpendicular to the hydrostatic line. Three parameters ($\xi ,\rho ,\theta$) represent the stress point in the plane. The following are these parameters’ expressions [33]:

$$\xi =\frac{I_1}{\sqrt{3}}$$

(1)

$$\rho =\sqrt{2J_2}$$

(2)

$$\theta =\frac{1}{3}{\cos }^{-1}\left(\frac{3\sqrt{3}}{2}\frac{J_3}{J_2^{3/2}}\right)$$

(3)

where $I_1$ = first invariant of the stress tensor, $I_1={{\displaystyle \sigma }}_{ii}$; $J_2$ = second invariant of the stress deviator tensor, ${{\displaystyle J}}_2=s_{ij}{s}_{ij}$/2; ${{\displaystyle J}}_3$ = third invariant of the stress deviator tensor, ${{\displaystyle J}}_3=s_{ij}{s}_{jk}{s}_{ki}/3$; $s_{ij}={\sigma }_{ij}-{\sigma }_{ij}{\delta }_{ij}/3$; and ${\delta }_{ij}$ = Kronecker delta. The stress point’s location on the hydrostatic line is represented by the symbol $\xi$, its distance from the line, by the symbol $\rho$, and its angular direction by the symbol $\theta$. The range of the Lode angle is [0, π/3].

The principal stress and $\xi ,\rho ,\theta$ have the following relationship [16]:

$$\left\{\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\end{array}\right\}=\left\{\begin{array}{l}\frac{\xi }{\sqrt{3}}\\ \frac{\xi }{\sqrt{3}}\\ \frac{\xi }{\sqrt{3}}\end{array}\right\}+\sqrt{\frac{2}{3}}\rho \left\{\begin{array}{l}\cos \theta \\ \cos (\theta -\frac{2}{3}\pi )\\ \cos (\theta +\frac{2}{3}\pi )\end{array}\right\}$$

(4)

3. Three-Dimensional Generalized Mohr-Coulomb Plastic Model

The effect of ${\sigma }_2$ on the shear strength is not considered by the M-C failure criterion. Alternatively, this condition might be seen as just being applicable in the case where ${\sigma }_2={\sigma }_3$. However, it is impossible to overlook the impact of ${\sigma }_2$ on the geotechnical materials’ shear strength [18]. For over-consolidated clays, plane strain angles, for example, are typically 10 percent higher than the comparable value obtained using traditional triaxial tests. Additionally, the peak value of $\varphi$ in plane strain for dense sands may be 4° or 5° higher [34].

${\sigma }_2$ always varies in the interval of $\left[{\sigma }_3,{\sigma }_1\right]$. From past experiences, the geotechnical materials’ shear strength gradually increases with the increase in ${\sigma }_2$. That is, when ${\sigma }_2$ = ${\sigma }_3$, the shear strength is the minimum; when ${\sigma }_2$ = ${\sigma }_1$, it is the maximum. Mogi [35] also observed similar situations.

To reflect the influence of ${\sigma }_2$ on the shear strength, we assume the M-C criterion still holds in form, but the cohesion $c$ and friction factor f are related to the coefficient $b$ of the intermediate principle stress ${\sigma }_2$ by

$$c_b=c\left(b\right)$$

and

$$f_b=f\left(b\right),$$

respectively; with $b$ defined as

$$b=\frac{{\sigma }_2-{\sigma }_3}{{\sigma }_1-{\sigma }_3}$$

(5)

Since ${\sigma }_1\ge {\sigma }_2\ge {\sigma }_3$ is specified, $0\le b\le 1$.

Therefore, the generalized Mohr-Coulomb criteria (GMC) adopts this form in the Mohr stress space:

$${\tau }_f=c\left(b\right)+f\left(b\right){\sigma }_n$$

(6)

The simplest way to obtain functions $c\left(b\right)$ and $f\left(b\right)$ is to apply the linear Lagrange interpolation to a series of data, $\left(c_i,f_i;b_i\right)$, $i=0,\cdots$, with $f_i=tg{\varphi }_i$, which are obtained from conventional or true triaxial tests.

Particularly, we can interpolate the shear strength parameter of $\left(c_0,f_0\right)$ from the conventional triaxial compression test with $b$ = 0, and that of $\left(c_1,f_1\right)$ from the conventional triaxial extension test with $b$ = 1, to obtain $c\left(b\right)$ and $f\left(b\right)$, reading

$$c\left(b\right)=c_0{l}_0\left(b\right)+c_1{l}_1\left(b\right)$$

(7)

and

$$f\left(b\right)=f_0{l}_0\left(b\right)+f_1{l}_1\left(b\right)$$

(8)

where $l_0\left(b\right)$ and $l_1\left(b\right)$ are the linear Lagrange interpolation functions, i.e.,

$$l_0\left(b\right)=1-b,l_1\left(b\right)=b$$

(9)

In this way, a three-dimensional failure criterion that considers the intermediate principal stress’s influence may be developed using just a piece of traditional triaxial test equipment.

We can obtain the yield surface of the GMC-L in the principal stress space:

$$\left({\sigma }_1-{\sigma }_3\right)-\left({\sigma }_1+{\sigma }_3\right)\sin {\varphi }_b-2c_b\cos {\varphi }_b=0$$

(10)

Each ${\varphi }_b$ corresponds to a unique b value.

Figure 2a shows the yield image of the GMC-L criterion and the M-C failure criterion on the principal stress space. The GMC-L and M-C criteria are similar in shape and are hexagonal pyramids. Figure 2b shows the comparison of the GMC-L and M-C criteria on the π plane. The GMC-L surface contains the M-C surface inside but coincides with the M-C surface at vertex A (b = 0).

4. Convexity of GMC-L Failure Criterion

According to the Drucker the stabilized materials’ failure surface must be convex [36]. Therefore, to ensure the failure surface’s convexity, the GMC-L criterion needs to meet certain parameter conditions. Since the GMC-L criterion is linear on the meridian plane, it is only necessary to consider its convexity on the π plane.

Similarly, combined with Equations (4) and (10), we can obtain the GMC criterion’s expression using parameters $\xi ,\rho$, and $\theta$:

$$\begin{array}{l}\xi \sin {\varphi }_b+\frac{\sqrt{2}}{4}\rho [(\cos \theta -\sqrt{3}\sin \theta )\sin {\varphi }_b-(3\cos \theta +\sqrt{3}\sin \theta )]\\ +\sqrt{3}{c}_b\cos {\varphi }_b=0\end{array}$$

(11)

When $\theta =0°$, the geotechnical material is in the triaxial compression state; while b = 0, Equation (11) is as follows:

$$\xi \sin {\varphi }_0+\frac{\sqrt{2}}{4}\rho (\sin {\varphi }_0-3)+\sqrt{3}{c}_0\cos {\varphi }_0=0$$

(12)

When $\theta =60°$, the material is in the triaxial extension state; while b = 1, Equation (11) can be written as follows:

$$\xi \sin {\varphi }_1-\frac{\sqrt{2}}{4}\rho (\sin {\varphi }_1+3)+\sqrt{3}{c}_1\cos {\varphi }_1=0$$

(13)

Assuming the cohesion c_b = 0, according to Equation (12), the triaxial compression strength vector ${\rho }_c$ of GMC-L on the π plane can be written as follows:

$${\rho }_c=\frac{2\sqrt{2}\sin {\varphi }_0}{3-\sin {\varphi }_0}\xi$$

(14)

According to Equation (13), the triaxial extension strength vector ${\rho }_e$ of GMC-L can be provided in a similar manner:

$${\rho }_e=\frac{2\sqrt{2}\sin {\varphi }_1}{3+\sin {\varphi }_1}\xi$$

(15)

As shown in Figure 3, OP = ${\rho }_c$ and OC = ${\rho }_e$. To ensure the hexagon failure surface’s convexity on the π plane, according to the geometric relationship, OB ≤ OP ≤ OA is obtained. According to the relationship, the strength parameters ${\varphi }_0$ and ${\varphi }_1$ must meet the conditions:

$$\frac{\sin {\varphi }_1}{3+\sin {\varphi }_1}\le \frac{2\sin {\varphi }_0}{3-\sin {\varphi }_0}\le \frac{4\sin {\varphi }_1}{3+\sin {\varphi }_1}$$

(16)

Expanding and rearranging, Equation (16) becomes

$$\frac{\sin {\varphi }_1}{2+\sin {\varphi }_1}\le \sin {\varphi }_0\le \frac{2\sin {\varphi }_1}{1+\sin {\varphi }_1}$$

(17)

Some studies [37,38,39] also show that the material’s strength envelope surface should meet the convexity requirements. Therefore, when ${\varphi }_0$ and ${\varphi }_1$ satisfy Equation (17), the GMC-L criterion satisfies convexity.

5. Strength Parameters Determined on Meridian Plane

5.1. Failure Function on Meridian Plane

In Section 3, we obtained the GMC-L criterion’s expression in the principal stress space. Here, the criterion’s form in the triaxial compression and extension are considered.

When considering the compression state (b = 0, ${\sigma }_1>{\sigma }_2={\sigma }_3$), Formula (10) can be written as shown:

$$\left({\sigma }_1-{\sigma }_3\right)=\left({\sigma }_1+{\sigma }_3\right)\sin {\varphi }_0+2c_0\cos {\varphi }_0$$

(18)

which can be expressed using p and q:

$$q=\frac{6\sin {\varphi }_0}{3-\sin {\varphi }_0}p+\frac{6c_0\cos {\varphi }_0}{3-\sin {\varphi }_0}$$

(19)

where p = mean stress, $p=\frac{1}{3}({\sigma }_1+{\sigma }_2+{\sigma }_3)$ and q = deviator stress, $q=\sqrt{\frac{1}{2}\left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2\right]}$.

Equation (19) is the GMC-L’s failure function on the meridian plane of the triaxial compression. Similarly, we can obtain the criteria’s formula on the p-q plane under triaxial extension (b = 1, ${\sigma }_1={\sigma }_2>{\sigma }_3$):

$$q=\frac{6\sin {\varphi }_1}{3+\sin {\varphi }_1}p+\frac{6c_1\cos {\varphi }_1}{3+\sin {\varphi }_1}$$

(20)

According to the GMC-L criterion, Figure 4 depicts the stress state of the triaxial compression and extension. GMC-L is a straight line on the meridian plane, as is shown in Figure 4. The triaxial compression and extension lines intersect at point T, and the stress state of the true triaxial is situated in the grey area. The strength parameters of $\varphi$ (${\varphi }_0$ and ${\varphi }_1$) and $c$($c_0$ and $c_1$) can be simply obtained using the least square method with the triaxial compression and extension test data.

5.2. Determination and Verification of Strength Parameters

The GMC-L criterion can be established only using triaxial compression and extension tests. We selected the triaxial compression and extension test data in the literature [40] to verify the GMC-L criterion’s applicability.

Figure 5 shows the strength test data of Monterey sand on the p-q plane and the corresponding GMC-L failure surfaces. It can be seen from Figure 5 that the GMC-L criterion’s fitting effect on Monterey sand is satisfactory, and that the correlation coefficient R² > 0.96. According to Equations (19) and (20), the internal friction angle under triaxial compression is 37° and under triaxial extension is 46°.

We also calculated the internal friction angle of triaxial compression and extension of other geotechnical materials [41,42,43] using the GMC-L criterion, as shown in

Table 1. Strength parameters of soils under GMC-L criterion.

S.No	Soil Name	${\mathit{\varphi }}_0 (°)$	${\mathit{\varphi }}_1 (°)$	Source
1	Monterey sand	37	46	Reddy and Saxena (1993)
2	Santa Monica sand	40.92	46.96	Lade and Wang (2001)
3	Shanghai fine sand	45.26	52.25	Hu et al. (2011)
4	Normally-consolidated Fujinomori clay	35.19	37.31	Nakai et al. (1986)

.

Table 1 shows that the triaxial compression friction angle ${\varphi }_0$ of the geotechnical materials is smaller than that of the triaxial extension friction angle ${\varphi }_1$, which is consistent with the results observed in other studies [44,45]. With the increase in ${\sigma }_2$, the internal friction angle of the geotechnical materials also gradually increases. Therefore, the influence of ${\sigma }_2$ on the soils’ strengths can be referred to as ${\varphi }_0$ and ${\varphi }_1$. For different soil types, the influence of ${\sigma }_2$ is different. For example, for Santa Monica sand, the internal friction angle of the triaxial extension is about 14.8% higher than that of the compression, but, for Fujinomori clay, this value is about 6%. The GMC-L criterion has good applicability to geotechnical materials.

We also verified the GMC criterion’s applicability for rock materials by calculating the failure surfaces corresponding to the GMC-L criterion of different rocks on the p-q plane and comparing them with the experimental data in the literature [35,46,47].

The strength test data of different rocks on the p-q plane and the corresponding GMC-L failure surfaces are shown in Figure 6. The GMC-L criterion clearly has good applicability to rocks.

Table 2. Strength parameters of soils under GMC-L criterion.

S.No	Rock Name	${\mathit{\varphi }}_0 (°)$	${\mathit{\varphi }}_1 (°)$	c0 (MPa)	c1 (MPa)	Source
1	Darley sandstone	35.3	39.6	26.8	31.3	Murrell (1965)
2	Westerly granite	53.4	59.7	41.4	52.7	Mogi (1967)
3	Dulham dolomite	46.8	51.2	49.5	57.6	Mogi (1967)
4	Red Wildmoor sandstone	28.5	35.6	6.0	7.9	Papamichos et al. (2000)

displays the GMC-L strength parameters for rocks.

The internal friction angle and rock cohesion are smaller in triaxial compression than they are in extension, as can be determined from Table 2; therefore, as ${\sigma }_2$ rises, the parameters governing the strength will gradually increase. The experimental outcomes produced by [48] were consistent. We could also observe that the intermediate principal stress’s effects on the strength metrics of various rocks vary. For example, for Red Wildmoor sandstone, the internal friction angle of the triaxial extension is about 25% higher than that of compression, but, for Darley sandstone, this value is about 12.2%.

The research above also reveals that the strength parameters required to create the GMC criterion for geotechnical materials can be acquired by fitting the test results on the p-q plane. A true triaxial failure criterion can be easily established simply using a set of traditional triaxial compression tests and a set of triaxial extension tests, which greatly reduces the workload of the test and has certain engineering practicability.

6. Experimental Verification of GMC-L Criterion

This section describes the comparison of the failure criterion on the π plane with the previous test results of the geotechnical materials. True triaxial testing was carried out by Lade and Duncan [49] on cube specimens of Monterey sand. The principal stress component ${\sigma }_1$ was increased and ${\sigma }_3$ was kept constant (${\sigma }_3=58.8 \mathrm{kPa}$), while another principal stress ${\sigma }_2$ was adjusted so that the intermediate principal stress coefficient b was kept constant. The test results of the loose Monterey sand were plotted on the π plane in Figure 7. The actual test results only exist in the range from θ = 0° to 60°, and the plots in the range of 60° < θ < 360° are the mirror images assuming the material’s isotropy. In Figure 7, the failure envelopes generated using the L–D, SMP, M-C, and GMC-L are also indicated.

In Figure 7, for loose Monterey sand, the triaxial compression strengths predicted using the four failure criteria are equal, but the triaxial extension strength predicted using the GMC-L criterion is significantly higher than those predicted using the SMP and M-C criteria. The M-C criteria is at the innermost position because it disregards the intermediate principal stress’s impact. In a three-dimensional stress condition, it is discovered that the GMC-L may fairly characterize the geotechnical materials’ properties.

Figure 8 shows the test data and GMC-L prediction results of the Dense Monterey sand and Fujinomori clay [42]. The proposed criteria clearly agree with the experiment’s results.

7. Numerical Implementation

7.1. Elastoplastic Constitutive Relations’ General Manifestation

The current elastic zone’s limit is set by the yield surface in the stress space. A stress point within this yield surface is in an elastic state and only has elastic properties. The point’s stress state on the yield surface is in a plastic state with elastic or elastoplastic characteristics. A loading process that is determined by the loading criterion can be obtained if the stress state tends to move away from the yield surface.

The loading condition in the proposed constitutive model is established using the loading criterion put out by Qn [50] and Chen et al. [51]. The stress tensor is determined using the generalized Hooke’s law as follows:

$${\sigma }_{ij}=D_{ijkl}^e\left({\epsilon }_{kl}-{\epsilon }_{kl}^p\right)$$

(21)

where ${\sigma }_{ij}$ represents the stress tensor; ${\epsilon }_{kl}$ represents the strain tensor; ${\epsilon }_{kl}^p$ represents the plastic strain tensor; and $D_{ijkl}^e$ is the fourth-order elastic material tensor, which has the form [52]:

$$D_{ijkl}^{\mathrm{e}}=\frac{E}{2(1+v)}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)+\frac{vE}{(1+v)(1-2v)}{\delta }_{ij}{\delta }_{kl}$$

(22)

where $v$ is Poisson’s ratio, E is Young’s modulus, and ${\delta }_{ij}$ is the Kronecker symbol.

One may write the yield function f as follows:

$$f\left({\sigma }_{ij},H\right)=F\left({\epsilon }_{kl},{\epsilon }_{kl}^p,H\right)$$

(23)

where H is hardening parameter. Taking the partial differential of Equation (21) gives the following:

$$\frac{\partial F}{\partial {\epsilon }_{ij}}=\frac{\partial f}{\partial {\sigma }_{kl}}{D}_{ijkl}^{\mathrm{e}}$$

(24)

Equation (25), which provides the loading criteria, provides their physical meanings (which are shown in Figure 9):

$$\frac{\partial F}{\partial {\epsilon }_{ij}}d{\epsilon }_{ij}=\left\{\begin{array}{ll}\frac{\partial f}{\partial {\sigma }_{ij}}{D}_{ijkl}^{\mathrm{e}}d{\epsilon }_{kl}>0\hfill & \mathrm{Loading} \hfill \\ \frac{\partial f}{\partial {\sigma }_{ij}}{D}_{ijkl}^{\mathrm{e}}d{\epsilon }_{kl}=0\hfill & \mathrm{Neutral} \mathrm{loading} \hfill \\ \frac{\partial f}{\partial {\sigma }_{ij}}{D}_{ijkl}^{\mathrm{e}}d{\epsilon }_{kl}<0\hfill & \mathrm{Unloading} \hfill \end{array}\right.$$

(25)

According to Hooke’s law, the relationship between the elastic strain increases and the stress increment during the pure elastic deformation process is as follows:

$$d{\sigma }_{ij}=D_{ijkl}^{e}d{\epsilon }_{kl}^e$$

(26)

The whole strain increment can be split into elastic and plastic strain increments, according to the elastic-plastic theory:

$$d{\epsilon }_{ij}=d{\epsilon }_{ij}^e+d{\epsilon }_{ij}^p$$

(27)

where $d{\epsilon }_{ij}^p$ stands for an increment in the plastic strain ${\epsilon }_{ij}^p$, and can be computed through the associated flow rule with the normality condition, as

$$d{\epsilon }_{ij}^p=\lambda \frac{\partial f}{\partial {\sigma }_{ij}}$$

(28)

with $\lambda$ as a scalar representing the plastic strain increment’s magnitude, which is called the plastic scalar factor.

Combining Equations (26) and (27), the formula for the stress–strain increment relationship is as follows:

$$d{\sigma }_{ij}=D_{ijkl}^e\left(d{\epsilon }_{ij}-d{\epsilon }_{ij}^p\right)$$

(29)

7.2. The Partial Derivative of Principal Stress to Stress Component

The yield function is simply and intuitively expressed in the three-dimensional principal stress space. When calculating the first-order partial derivative of the yield function f to the stress component, it can be expressed as follows:

$$\frac{\partial f}{\partial {\sigma }_{ij}}=\frac{\partial f}{\partial {\sigma }_1}\frac{\partial {\sigma }_1}{\partial {\sigma }_{ij}}+\frac{\partial f}{\partial {\sigma }_2}\frac{\partial {\sigma }_2}{\partial {\sigma }_{ij}}+\frac{\partial f}{\partial {\sigma }_3}\frac{\partial {\sigma }_3}{\partial {\sigma }_{ij}}$$

(30)

In Equation (30), for the partial derivative of the principal stress to the stress component $\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}$ (where k = 1,2,3), it can be obtained from the characteristic equation:

$$\mathit{S}{\mathit{l}}_k={\sigma }_k{\mathit{l}}_k$$

(31)

The subscript k does not represent a summation symbol; ${\mathit{l}}_k$ represents the principal direction corresponding to the principal stress, and ${‖{\mathit{l}}_k‖}_2=1$. $\mathit{S}$ is a symmetric stress matrix:

$$\mathit{S}=\left[\begin{array}{ccc}{\sigma }_x& {\tau }_{xy}& {\tau }_{xz}\\ {\tau }_{yx}& {\sigma }_y& {\tau }_{yz}\\ {\tau }_{zx}& {\tau }_{zy}& {\sigma }_z\end{array}\right]$$

(32)

By multiplying both sides of Equation (31) by ${\mathit{l}}_k^T$,we can obtain the following:

$${\mathit{l}}_k^T\mathit{S}{l}_k={\mathit{l}}_k^T{\sigma }_k{\mathit{l}}_k$$

(33)

By simultaneously deriving the stress components ${\sigma }_{ij}$ on both sides of Equation (33), the following can be obtained:

$$\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}=\frac{\partial {\mathit{l}}_k^T}{\partial {\sigma }_{ij}}\mathit{S}{\mathit{l}}_k+{\mathit{l}}_k^T\frac{\partial \mathit{S}}{\partial {\sigma }_{ij}}{\mathit{l}}_k+{\mathit{l}}_k^T\mathit{S}\frac{\partial {\mathit{l}}_k}{\partial {\sigma }_{ij}}$$

(34)

Since $\mathit{S}={\mathit{S}}^T$, Equation (34) can be written as shown:

$$\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}=2\frac{\partial {\mathit{l}}_k^T}{\partial {\sigma }_{ij}}\mathit{S}{\mathit{l}}_k+{\mathit{l}}_k^T\frac{\partial \mathit{S}}{\partial {\sigma }_{ij}}{\mathit{l}}_k$$

(35)

Substituting Equation (31) into Equation (35), we can obtain the following:

$$\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}=2{\sigma }_k\frac{\partial {\mathit{l}}_k^T}{\partial {\sigma }_{ij}}{\mathit{l}}_k+{\mathit{l}}_k^T\frac{\partial \mathit{S}}{\partial {\sigma }_{ij}}{\mathit{l}}_k$$

(36)

Using ${\mathit{l}}_k^T{\mathit{l}}_k=1$, we can obtain $\frac{\partial {\mathit{l}}_k^T}{\partial {\sigma }_{ij}}{\mathit{l}}_k=0$; thus, for Equation (36), there is the following condition:

$$\frac{\partial {\sigma }_k}{\partial {\sigma }_{ij}}={\mathit{l}}_k^T\frac{\partial \mathit{S}}{\partial {\sigma }_{ij}}{\mathit{l}}_k=(2-{\delta }_{ij})l_k^i{l}_k^j$$

(37)

where $l_k^i$ represents the i-th element in the principal stress direction ${\mathit{l}}_k$. It can be seen from Equation (37) that the derived partial derivative of the principal stress to the stress component is simple and easy to program.

7.3. The Mixed Complementarity Problem and its GSPC Algorithm

The elastoplastic constitutive relation in the rate form can be expressed as:

$$\dot{\mathbf{\sigma }}={\mathit{D}}^e\left(\dot{\mathbf{\epsilon }}-{\dot{\mathbf{\epsilon }}}^p\right)$$

(38)

where $\dot{\mathbf{\sigma }}$, $\dot{\mathbf{\epsilon }}$, and ${\dot{\mathbf{\epsilon }}}^p$ represent the stress, strain, and plastic strain vector, respectively; ${\mathit{D}}^e=D_{ijkl}^e$ represents the elastic matrix.

Considering the time discretization on the interval [0, T], the relationship between the strain rate $\dot{\mathbf{\epsilon }}$ and T is as follows:

$$\dot{\mathbf{\epsilon }}=\frac{\Delta \mathbf{\epsilon }}{T}$$

(39)

where $\Delta \mathbf{\epsilon }$ denotes the incremental strain vector.

To establish the complementary equation of the elastic-plastic problem, the backward Euler algorithm is used to rewrite the constitutive Equation (38) as follows:

$$\mathbf{\sigma }={\mathbf{\sigma }}^e-{\mathbf{\sigma }}^p$$

(40)

where ${\mathbf{\sigma }}^e$ represents the test stress, which can be written as:

$${\mathbf{\sigma }}^e={\mathbf{\sigma }}_0+{\mathit{D}}^e\Delta \mathbf{\epsilon }$$

(41)

where ${\mathbf{\sigma }}_0$ denotes the total stress at the beginning of the load phase; the ${\mathbf{\sigma }}^p$ expression is as follows:

$${\mathbf{\sigma }}^p={\mathit{D}}^e\lambda \frac{\partial f}{\partial \mathbf{\sigma }}$$

(42)

The mixed complementarity issue that the plastic constitutive integral defines may be expressed as follows:

$$\left\{\begin{array}{l}0\le \lambda \perp f_I\left(\mathbf{\sigma },\lambda \right)\ge 0\\ f_E\left(\mathbf{\sigma },\lambda \right)=0\end{array}\right.$$

(43)

This problem is abbreviated as MiCP($f_I$,$f_E$); $f_I$ and $f_E$ are defined as follows:

$$\left\{\begin{array}{l}{f}_I\left(\mathbf{\sigma },\lambda \right)=-f(\mathbf{\sigma })\\ f_E\left(\mathbf{\sigma },\lambda \right)=\mathbf{\sigma }+{\mathbf{\sigma }}_p-{\mathbf{\sigma }}_e\end{array}\right.$$

(44)

For the mixed complementarity problem, Zheng et al. [8] designed a GSPC algorithm based on the Gauss-Seidel idea and projection-contraction algorithm. The GSPC algorithm is easy to program, can converge for any size of strain increment, and does not need to calculate the potential function’s Hessian matrix. The computational efficiency is higher than that of the return mapping algorithm.

The GSPC algorithm’s code is as follows:

$$\begin{array}{l}\mathrm{Step}0:{\beta }_I={\beta }_E=1;k=0;\\ \mathrm{Step}1:\overline{\lambda }=\max \left[\lambda -{\beta }_I{f}_I(\lambda ,\mathit{\sigma }),0\right];\\ \overline{\mathit{\sigma }}=\mathit{\sigma }-{\beta }_E{f}_E(\overline{\lambda },\mathit{\sigma })\\ \text{If }{‖\overline{\lambda }-\lambda ‖}_{\infty }\le {\epsilon }_{\lambda }{‖\overline{\lambda }‖}_{\infty }\text{ and }{‖\overline{\mathit{\sigma }}-\mathit{\sigma }‖}_{\infty }\le {\epsilon }_{\mathit{\sigma }}{‖\overline{\mathit{\sigma }}‖}_{\infty },\text{then }\lambda =\overline{\lambda };\mathit{\sigma }=\overline{\mathit{\sigma }};\text{break};\\ r_{\lambda }=\frac{{\beta }_I{‖f_I(\lambda ,\mathit{\sigma })-f_I(\overline{\lambda },\mathit{\sigma })‖}_2}{{‖\lambda -\overline{\lambda }‖}_2}\\ \text{while }{r}_{\lambda }>v\\ {\beta }_I=\frac{2}{3}{\beta }_I\min \left(1,\frac{1}{r_{\lambda }}\right);\overline{\lambda }=\max \left[\lambda -{\beta }_I{f}_I(\lambda ,\mathit{\sigma }),0\right]\text{;}\\ r_{\lambda }=\frac{{\beta }_I{‖f_I(\lambda ,\mathit{\sigma })-f_I(\overline{\lambda },\mathit{\sigma })‖}_2}{{‖\lambda -\overline{\lambda }‖}_2};\\ \text{end(while)}\\ r_{\mathit{\sigma }}=\frac{{\beta }_E{‖f_E(\overline{\lambda },\mathit{\sigma })-f_E(\overline{\lambda },\overline{\mathit{\sigma }})‖}_2}{{‖\mathit{\sigma }-\overline{\mathit{\sigma }}‖}_2};\\ \text{while }{r}_{\mathit{\sigma }}>v\\ {\beta }_E=\frac{2}{3}{\beta }_E\min \left(1,\frac{1}{r_{\mathit{\sigma }}}\right);\overline{\mathit{\sigma }}=\mathit{\sigma }-{\beta }_E{f}_E(\overline{\lambda },\mathit{\sigma })\text{; }\\ r_{\mathit{\sigma }}=\frac{{\beta }_E{‖f_E(\overline{\lambda },\mathit{\sigma })-f_E(\overline{\lambda },\overline{\mathit{\sigma }})‖}_2}{{‖\mathit{\sigma }-\overline{\mathit{\sigma }}‖}_2};\\ \text{end(while)};\\ d_{\lambda }(\lambda ,\overline{\lambda })=(\lambda -\overline{\lambda })-{\beta }_I\left[f_I(\lambda ,\mathit{\sigma })-f_I(\overline{\lambda },\mathit{\sigma })\right];\\ \alpha =\frac{{(\lambda -\overline{\lambda })}^T{d}_{\lambda }(\lambda ,\overline{\lambda })}{{‖d_{\lambda }(\lambda ,\overline{\lambda })‖}_2^2};\lambda =\lambda -\gamma \alpha d_{\lambda }(\lambda ,\overline{\lambda }); \\ \text{If }{r}_{\lambda }\le \mu \text{ then }{\beta }_I=1.5{\beta }_I;\\ d_{\mathit{\sigma }}(\mathit{\sigma },\overline{\mathit{\sigma }})=(\mathit{\sigma }-\overline{\mathit{\sigma }})-{\beta }_E\left[f_E(\overline{\lambda },\mathit{\sigma })-f_E(\overline{\lambda },\overline{\mathit{\sigma }})\right]\text{; }\\ \alpha =\frac{{(\mathit{\sigma }-\overline{\mathit{\sigma }})}^T{d}_{\mathit{\sigma }}(\mathit{\sigma },\overline{\mathit{\sigma }})}{{‖d_{\mathit{\sigma }}(\mathit{\sigma },\overline{\mathit{\sigma }})‖}_2^2};\mathit{\sigma }=\mathit{\sigma }-\gamma \alpha d_{\mathit{\sigma }}(\mathit{\sigma },\overline{\mathit{\sigma }});\\ \text{If }{r}_{\mathit{\sigma }}\le \mu \text{ then }{\beta }_E=1.5{\beta }_E;\\ \mathrm{Step}2:k=k+1\text{; go to Step 1}\end{array}$$

The GSPC algorithm has the following characteristics: (1) The GSPC converges for any size of the strain increment $\Delta \epsilon$, and the implicit return mapping algorithm calculates the stress $\mathbf{\sigma }$ and the plastic multiplier $\mathit{\lambda }$ using the Newton technique. When the strain increment is small enough, the convergence can be ensured because the Newton technique is locally convergent. (2) The GSPC algorithm’s computational efficiency is higher than the return mapping algorithm because it is simple to program and the second derivative of the potential function does not need computing. The return mapping algorithm needs to first determine the stress path and yield surface corner points, as well as whether those yield surfaces are active, before forming the plastic potential function’s Jacobian matrix. The Hessian matrix of the plastic potential function must be determined to derive the Jacobian matrix. Usually, the return mapping algorithm requires a lot of trial and error for this method to provide accurate results.

Except for ${\varphi }_b$ being replaced with the expansion angle ${\psi }_b$, the expression of f and the plastic potential function g are consistent. When using the GSPC integral algorithm to update the stress and calculate the plastic multiplier, the first-order derivative of the yield function $F_{GMC}$ to stress is required (refer to Appendix A).

ABAQUS provides users with many user subroutine interfaces, thus allowing users to write the required user subroutines according to the corresponding problems encountered in the actual project. The UMAT subroutine can specify the material’s characteristics. The Jacobian matrix, the major component of UMAT, updates the stress increments and associated state variables in accordance with the strain increments provided by the ABAQUS core program [53,54].

The ABAQUS main program runs the UMAT subroutine at the element’s integration point at the start of each incremental loading step and inputs the strain increment, time step, and load increment, as well as the stress, strain, and other variables related to the solution process in the current known state. The UMAT subroutine [55] solves the stress increment according to the constitutive equation, updates the stress and other related variables, and provides the Jacobian matrix to the ABAQUS main program to generate the global stiffness matrix. The main software solves the displacement increment by combining the current load increment and then it determines whether it is balanced. If the specified error is not satisfied, ABAQUS will iterate until convergence and then solve the next incremental step. According to the above ideas, the main process of the UMAT subroutine in this paper is shown in Figure 10.

8. Numerical Examples

8.1. Uniaxial Compression

The GMC yield criterion-based elastoplastic constitutive model’s UMAT subroutine is studied and constructed, and the standard soil mechanics experiments conducted under uniaxial compression are analyzed using finite elements. When the strength parameters ϕ₀ = ϕ₁ and c₀ = c₁, the GMC model degenerates into the Mohr-Coulomb model. The established UMAT subroutine’s computational capability, correctness, and effectiveness are confirmed by comparing the calculations of the subroutine and the Mohr-Coulomb model.

The specimen’s element type is a completely integrated C3D8 solid element, and Figure 11 depicts the calculation model. The specimen is 100 mm × 100 mm × 100 mm in size. The established calculation model is a cube, and the applied uniform load is 6.0 MPa.

Table 3. Material parameters.

Model	Elastic Modulus (Mpa)	Poisson’s Ratio	${\mathit{\varphi }}_0 (°)$	${\mathit{\varphi }}_1 (°)$	c0 (Mpa)	c1 (Mpa)
M-C	1000	0.3	45	-	1.2	-
GMC	1000	0.3	45	45	1.2	1.2

displays the material properties using the associated flow rule.

The stress ${\sigma }_{33}$ and plastic strain ${\epsilon }_{33}^p$ nephogram of the GMC model and M-C model incorporated in ABAQUS under uniaxial compression are shown in Figure 12 and Figure 13, respectively.

It is evident from Figure 12 and Figure 13 that the stress nephograms produced using the M-C model integrated into ABAQUS and the UMAT subroutine are essentially identical; the plastic strain nephogram obtained using the two models also have good consistency, and the plastic strain calculated using the UMAT subroutine is slightly larger. The difference between the calculated results is mainly because the plastic potential function of the M-C model embedded in ABAQUS is a continuous and smooth elliptic function, which is not consistent with the yield function. Nevertheless, the accuracies of both the stress and strain meet the calculation requirements.

8.2. Uniaxial Compression

To check the UMAT subroutine’s accuracy in the multi-element calculation and practical application, an ideal underground cavern under hydrostatic pressure is used as an example. The subroutine’s numerical solution is compared with the classical Fenner formula’s rigorous solution for analyzing such problems.

The tunnel has a 10.0 m radius and is 90.0 m below ground. The model adopts a C3D8 solid element. Due to the structure’s symmetry, 1/4 of the model is calculated and a force of 30.0 MPa is applied in the x and y directions.

Table 4. Parameters used for numerical simulation.

Model	Elastic Modulus (MPa)	Poisson’s Ratio	${\mathit{\varphi }}_0 (°)$	${\mathit{\varphi }}_1 (°)$	c0 (MPa)	c1 (MPa)
GMC	2000	0.3	30	30	6.0	6.0

displays the surrounding rock-related parameters, and Figure 14 displays the computational model.

Under ideal elastoplasticity, the internal stress of the surrounding rock made of M-C material can be calculated using the conventional Fenner formula [56]. The expression is as follows.

The plastic zone’s stress expression:

$${\sigma }_r^p=c\cdot \cot \varphi \cdot \left[{\left(\frac{r}{R_0}\right)}^{\frac{2\sin \varphi }{1-\sin \varphi }}-1\right]$$

(45)

$${\sigma }_{\theta }^p=c\cdot \cot \varphi \cdot \left[\frac{1+\sin \varphi }{1-\sin \varphi }{\left(\frac{r}{R_0}\right)}^{\frac{2\sin \varphi }{1-\sin \varphi }}-1\right]$$

(46)

$$R_1=R_0{\left[\left(\frac{\left(p_0+c\cdot \cot \varphi \right)\left(1-\sin \varphi \right)}{c\cdot \cot \varphi }\right)\right]}^{\frac{1-\sin \varphi }{2\sin \varphi }}$$

(47)

The elastic zone’s stress expression:

$${\sigma }_r^e=p_0\cdot \left(1-\frac{R_1^2}{r^2}\right)+c\cdot \cot \varphi \left[{\left(\frac{R_1}{R_0}\right)}^{\frac{2\sin \varphi }{1-\sin \varphi }}-1\right]\frac{R_1^2}{r^2}$$

(48)

$${\sigma }_{\theta }^e=p_0\left(1+\frac{R_1^2}{r^2}\right)-c\cdot \cot \varphi \left[{\left(\frac{R_1}{R_0}\right)}^{\frac{2\sin \varphi }{1-\sin \varphi }}-1\right]\frac{R_1^2}{r^2}$$

(49)

where c represents cohesion; $\varphi$ represents the internal friction angle; r represents the radius from a point in the surrounding rock to the tunnel’s center; $R_1$ represents the plastic zone’s radius; $R_0$ represents the circular tunnel’s radius; ${\sigma }_r^p$ and ${\sigma }_{\theta }^p$ are the plastic zone’s radial and circumferential stress, respectively; ${\sigma }_r^e$ and ${\sigma }_{\theta }^e$ are the elastic zone’s radial and circumferential stress, respectively; and $p_0$ is ground stress.

The comparison between the numerical and analytical solutions is shown in Figure 15. According to Figure 15, the numerical solution calculated using the UMAT subroutine is consistent with the exact solution calculated using the theoretical formula. Both the radial and circumferential stress have average relative errors of 1.63% and 1.67%, respectively. Even with flaws, the UMAT subroutine’s computation results are accurate enough to fulfill engineering standards. This also proves the GMC constitutive model program’s correctness from the perspective of multiple elements, and that it can solve practical engineering problems.

8.3. Modeling of Tunnel Excavation

Some scholars [57,58,59] have carried out corresponding research on tunnel excavation. However, to simulate tunnel excavation, most of them are based on the M-C model. The GMC model, in contrast to the M-C model, considers the impact of ${\sigma }_2$ and is able to fully utilize the material’s strength. Studying the tunnel excavation problem using the GMC model is therefore practical.

The model has dimensions of 20.0 m in length, width, and height, and the radius of the initial tunnel is 1.0m, which is drilled in a formation with in situ stresses ${\sigma }_V$ = 15.0 MPa, ${\sigma }_H$ = 30.0 MPa, and ${\sigma }_h$ (varying between ${\sigma }_V$ and ${\sigma }_H$). The M-C model is also used to mimic tunnel excavation as a comparison. The elastic mechanical parameters used are E = 30,000 MPa and $\nu$ = 0.25, and the plastic parameters are shown in

Table 5. Parameters used for numerical simulation.

Model	Elastic Modulus (MPa)	Poisson’s Ratio	${\mathit{\varphi }}_0 (°)$	${\mathit{\varphi }}_1 (°)$	c0 (MPa)	c1 (MPa)	$\mathit{\psi } (°)$
GMC	30,000	0.25	28.5	35.6	6.0	6.3	0
M-C	30,000	0.25	28.5	-	6.0	-	0

. A three-dimensional model was created in ABAQUS, as shown in Figure 16, to simulate the analysis as closely as possible. The model adopts the C3D8 solid element and divides into 13,640 units, with a total of 15,666 nodes.

The intermediate principal stress’s effect on the stability of the rock surrounding the tunnel is investigated by simulating tunnel excavation under various intermediate ground stresses (${\sigma }_h$). The ${\sigma }_h$ changes between ${\sigma }_V$ and ${\sigma }_H$, and the four working conditions of ${\sigma }_h$ = 15.0, 20.0, 25.0, and 30.0 MPa are simulated. Each condition is calculated using the GMC and M-C models, and the equivalent plastic strain nephograms from these two models are displayed in Figure 17 and Figure 18.

Since the M-C model does not reflect the influence of ${\sigma }_h$, as is shown in Figure 17, the range of the plastic zones computed under various ${\sigma }_h$ is essentially the same. According to Figure 18, the range of the surrounding rock’s plastic zone increasingly shrinks for the GMC model as ${\sigma }_h$ increases. This shows that the GMC model can accurately represent the intermediate principal stress’s effect on the plastic zone, since the surrounding rock’s strength continuously increases as the intermediate principal stress rises. Figure 19 shows the maximum equivalent plastic strain corresponding to different ${\sigma }_h$. The M-C model’s maximum equivalent plastic strain is almost unaltered for varied ${\sigma }_h$ values, proving that it does not account for the intermediate principal stress’s influence and that the calculation results are conservative. With an increase in ${\sigma }_h$, the maximum equivalent plastic strain of the GMC model decreases, thus demonstrating that the GMC model accounts for the intermediate principal stress’s effect. Most of the geotechnical materials used in practical engineering are in the three-dimensional stress state, which has a specific enhancing effect on the material’s strength. The GMC model may thus fully utilize the material’s strength potential while saving on engineering costs, which has some significance for practical engineering issues.

Figure 20 shows surrounding rock’s maximum stresses under different ${\sigma }_h$. The maximum radial stress and maximum circumferential stress of the surrounding rock using the M-C model essentially remain unchanged as ${\sigma }_h$ increases, whereas the maximum radial stress obtained using the GMC model decreases as ${\sigma }_h$ increases and the maximum circumferential stress barely changes. Under different ${\sigma }_h$, the maximum radial stress obtained using the GMC model is lower than the maximum radial stress using the M-C model, while the maximum circumferential stresses obtained using the two models are opposite. As a result, the GMC model can consider how the intermediate principal stress affects the surrounding rock stress. The GMC model’s outcomes appear more credible when contrasted with those of the M-C model.

9. Conclusions

This paper presents a three-dimensional failure criterion that considers intermediate principal stress under complex stress states. The proposed criterion has a simple linear expression, which is convenient for calculation, analysis, and application. The strength parameters needed to establish the GMC criterion can be obtained by fitting the triaxial compression and triaxial extension test data on the p-q plane. The failure surface fitted by the GMC criterion was revealed to be in good agreement with the experimental data by comparison with the plane’s classical failure criterion.

Based on the three-dimensional GMC elastoplastic model, a UMAT subroutine was developed using the ABAQUS platform, which uses a GSPC integral algorithm to program the constitutive relation. Through the uniaxial compression numerical experiment and the numerical experiment of the underground cavern under hydrostatic pressure, the subroutine’s correctness and effectiveness were verified. Finally, the excavation of a three-dimensional tunnel under various working conditions was numerically simulated. The results demonstrate that the surrounding rock’s stability steadily increases with increasing intermediate principal stress. The GMC model accounts for the intermediate principal stress’s impact, meaning it can depict the material’s strength properties under three-dimensional stress conditions and can fully utilize the material’s potential, which has some engineering significance in practice.