Our research group previously proposed a simple two-dimensional (2D) constitutive model for rocks to simulate not only the axial stress–axial strain relationship, but also the axial stress–lateral strain relationship, with few complicated equations. However, the failure envelope that it predicted was linear, and it did not consider the effect of the intermediate principal stress (⁠σ2⁠). In the present study, the authors modify this simple 2D model to have a convex failure criterion. Then, the model is extended to a simple three-dimensional (3D) model that well approximates true triaxial stress–strain curves for real rocks under specific values of σ2 and σ3 and uses only four parameters. However, the predicted peak stress–σ2 relationship is linear. Finally, a modified 3D model was developed, which exhibited the true triaxial convex failure criterion. The equations in this model are simpler than the conventional true triaxial failure criteria. The proposed models can be implemented with a finite element method to improve the design of rock structures.

A 3D constitutive model is required for the analysis of rock structures. However, as shown below, complex equations are used just to represent the true triaxial strength. This paper aims to develop a 3D constitutive model which is simple and can fairly represent the deformation and failure of rocks under the effects of the intermediate principal stress. The model would contribute to better rock structure design if combined with 3D numerical stress software.

The simplest three-dimensional (3D) constitutive model for rocks is
(1)σ1=λ+2με1+λε2+λε3,(2)σ2=λε1+λ+2με2+λε3,(3)σ3=λε1+λε2+λ+2με3,
where ε1⁠, ε2⁠, and ε3 are the maximum (axial), intermediate, and minimum (lateral) principal strains, respectively; σ1⁠, σ2⁠, and σ3 are the maximum (axial), intermediate, and minimum (lateral) principal stresses, respectively, and λ and μ are Lame’s constants (Figures 1(a) and 1(b)).

Elasto-brittle models can be derived by introducing peak and residual strengths (Figures 1(c) and 1(d)). However, the resulting stress–strain curves, which should be curvilinear for real rocks, are unrealistic polylines. More realistic stress–strain curves, which are needed for precise analyses of rock structures, can be obtained using various two-dimensional (2D) constitutive models [1–14].

A variable-compliance-type constitutive model [7] is (Figure 2)
(4)σ1=ωb1/n−3ε1−m/n−m,
where ω=dψ/dt⁠, ψ=ε1/σ1⁠, and n⁠, m⁠, and b are constants. This equation simulates the axial stress–axial strain relationships for classes I and II. However, no reasonable methods to evaluate the lateral strain have yet been provided. Actual rocks under compression show lateral dilatancy around peak stress in most cases, due to the growth of axial microcracks [15–18]; this growth is essential to the deformation and failure of rocks [19]. Therefore, the axial stress–lateral strain relationship is more important than the axial stress–axial strain relationship [20, 21].

By contrast, a simple 2D constitutive model for rocks was proposed by our research group [22] based on the finding that the trace of the axial strain, lateral strain, and axial stress is located on a plane in a 3D coordinate system for rocks such as Paleogene Kamisunagawa sandstone (Figure 3, [22]), Cretaceous Pombetsu sandstone, Paleogene Bibai sandstone, and Paleogene Inada granite [23]. The equations in the model are very simple, and class I and II strain-softening behaviors can be obtained with appropriate lateral strain behavior by introducing the strain-dependent elastic modulus.

For the 3D case, the strength of rocks is affected by σ2 ([24–27], and [28–32], Figure 4; [33]), and various criteria have been proposed to describe these effects ([25–27], and [30, 34–43]). These criteria use complicated equations and many parameters. For example, a modified Lade criterion proposed by Ewy [44] (Figure 5(a)) is represented as
(5)I1′3I3′=27+η,(6)I1′=σ1+s1−p0+σ2+s1−p0+σ3+s1−p0,(7)I3′=σ1+s1−p0+σ2+s1−p0+σ3+s1−p0,
where p0 is pore pressure, η is related to the internal friction, and s1 represents the cohesion of the rock. The peak strength first increases and then decreases with σ2 according to the modified Lade criterion.
As another example, a 3D failure criterion for rocks proposed by Li et al. [38] based on the Hoek–Brown criterion [45, 46] is expressed as
(8)σ1=bσ2+σ3b+1+σcmbσcbσ2+σ3b+1+sαlowσ2,(9)bσ2+σ1b+1=σ3+σcmbσcσ3+sαhighσ2,
where σc and b are the unconfined compressive strength and a constant, respectively, and
(10)mb=expGSI−10028−14Dmi,(11)s=expGSI−1009−3D,(12)α=12+16e−GSI/15−e−20/3,
where mi⁠, GSI⁠, and D are the material constant for intact rock, geological strength index, and disturbance factor, respectively. The effects of σ2 on rock strength are represented well in this model (Figure 5(b)); however, the parameter selection is complex, and the resulting peak stress curves are unrealistic polylines.

No constitutive equations can fully represent true triaxial stress–strain curves considering the effects of σ2⁠. Therefore, this paper proposes a 3D constitutive model with fewer parameters and consideration of these true triaxial stress-strain curves and convex failure envelopes. Herein, we briefly explain the simple 2D constitutive model of Fujii and Ishijima [22]. Then, this simple model is modified to include the convex failure envelope by introducing stress dependency for an elastic modulus. We give an example of its application to a pressurized thick-walled cylinder under the plane strain condition. Then, the simple 2D model is extended to a simple 3D one. Finally, this 3D model is modified to represent true triaxial stress–strain curves and convex failure envelopes.

The simple 2D model developed by our research group should be briefly explained before the 3D models because they were developed based on the 2D model, and it has not been published as a full paper in English in a peer-reviewed journal.

A plane in the axial stress–axial strain–lateral strain space (Figure 3) can be represented as
(13)σ1=Aε1−ε0+kAε3,
where σ1⁠, ε1⁠, and ε3 are the maximum (axial) principal stress, the maximum (axial) principal strain, and the minimum (lateral) principal strain. A is an elastic modulus, and ε0 and k are constants; these are equal to λ+2 μ⁠, 0, and ν/1−ν⁠, respectively, where ν is Poisson’s ratio for a linear elastic medium. Using a function of normal strain, A’ ε⁠, to replace the constant A⁠, and ignoring the constant ε0 in Equation (13) for simplicity, we get
(14)σ1=A′ε1ε1+kAε3.
By analogy, the following equation can be written for σ3⁠:
(15)σ3=A′ε3ε3+kAε1,
where σ3 is the minimum (lateral) principal stress. A’ε should increase and converge to a specific value as ε tends to +∞ (compression). The increase represents the increase in axial elastic modulus (⁠AA⁠) due to microcrack closure (Figure 6), as illustrated in Figure 7(a). A’ε should also decrease with expansion and should converge to a specific value as ε tends to −∞ (expansion). The decrease represents a decrease in the lateral elastic modulus (⁠AL⁠) under extension due to the initiation and growth of axial microcracks (Figure 6), as illustrated in Figure 7(b). These assumptions have already been confirmed numerically by a boundary element method combining displacement discontinuity [47] with the body force [49] elements in Fujii and Ishijima [23] and explained by various mechanisms (e.g., Figure 8, [17, 48]).
The following function (Figure 9) was chosen as A’ε to satisfy the above requirements:
(16)A′ε=Dtan−1Cεεs+1+π2+F,(17)D=1−k2Aπ,(18)F=k2A,
where εs is a constant whose value is approximately half the absolute value of the critical extensile strain (lateral strain value at the peak load point; [20]), C is a positive constant, and D and F are determined so that the residual strength becomes constant as ε1 and ε3 tend to ∞ and −∞⁠, respectively.
From Equations (14) to (18), an ultimate residual strength, σur⁠, for infinite strains can be calculated as
(19)σur=σ1ε1=∞,ε3=−∞=σ3k.

However, the residual strength σr at −4.5% lateral strain (Figure 10(f)), which is higher than the ultimate residual strength, will be used hereafter for convenience.

The simple 2D model requires only four parameters (⁠A⁠, k⁠, C⁠, and εs⁠) to simulate classes I and II nonlinear stress–strain behaviors (Figures 10(a)–10(e)), and tensile behavior can also be simulated seamlessly (Figure 10(a)). Curves representing such behavior can be obtained as follows, for example:

(a1) Negative ε3 increment (extensile) is assigned

(a2) ε1 is calculated by Equation (15)

(a3) A’ε1 and A’ε3 are calculated by Equation (16)

(a4) σ1 can be calculated by Equation (14)

(a5) Iterate steps (a2–a4) until the solution converges

(a6) Iterate steps (a1–a5) for loading until ε3 reaches -0.045

It is rather surprising that only a decrease in the lateral elastic modulus can induce strain-softening behavior. The mechanism of the strain softening is as follows. The elastic modulus A’ε3 decreases with lateral extensile strain, and this strain increases by Equation (15). According to Equation (14), the lateral extensile strain increment causes a decrease in axial stress. The peak load point appears when the stress increase due to the axial strain increment becomes the same as the stress decrease due to the lateral extensile strain increment. In step (a1), the axial strain increment can be assigned instead of the lateral extensile strain one to simulate class I behavior.

Considering the model’s parameters, A mainly affects strength and tangent modulus (Figure 10(b)). k mainly affects strength, the shape of the σ1–ε1 curve around the peak load, Poisson’s ratio, and the critical compressive strain (⁠ε1 value at the peak load point; [20]; Figure 10(c)). C mainly affects the shape of the stress–strain curves around the peak load, strength, and the critical strains (Figure 10(d)). εS mainly affects the critical strains and strength (Figure 10(e)).

One of the limitations of this model is that it predicts a linear failure envelope (Figure 10(f)). Therefore, we modify it in the next chapter. However, other results are similar to real rocks:

  • (i)

    Compressive strength increases with σ3 (Figure 10(a))

  • (ii)

    Compressive strength is larger than tensile strength (Figure 10(a))

  • (iii)

    Critical compressive strain increases with σ3 (Figure 10(a), [20])

  • (iv)

    Critical extensile strain is unaffected by σ3 (Figure 10(a), [20])

3.1. Introduction of k Dependent on Confining Pressure

We introduce k dependent on confining pressure (Figure 11(a)) to give a convex failure envelope based on the increase in the secant Poisson’s ratio with confining pressure under triaxial compression (Figure 11(b)):
(20)k=Δkπtan−11ασ3−σ0+π2+k0,
where k0 is the minimum value, ∆k is the amplitude, and σ0 and α are constants. The ultimate residual strength is unchanged from Equation (19).

The stress–strain curves obtained by setting default values of A=20 GPa⁠, C=1⁠, εs=0.002⁠, k0=0.20⁠, σ0=35 MPa⁠, α=20 MPa⁠, and ∆k=0.30 show strain-softening behaviors (Figure 12(a)) with a convex failure envelope (Figure 12(b)). The stress drop decreases, and the strain-softening behavior becomes mild as confining pressure increases; these are common features of rock deformation. Figure 12(b) also shows the failure envelope predicted by the Hoek–Brown strength criterion with σc=71.6 MPa⁠, mb=119.0⁠, s=1⁠, and α=0.104⁠, which is similar. Therefore, the equations are expected to approximate the stress–strain curves and convex failure envelope of various rocks. Considering the model parameters, k0 mainly affects the maximum stress (Figure 13(a)). ∆k (Figure 13(b)) and σ0 (Figure 13(c)) mainly affect the internal friction angle. α mainly affects the convexity of the failure envelope (Figure 13(d)).

3.2. Approximating Convex Failure Envelope for Test Results

The procedure for approximating ordinary triaxial experimental data using the modified 2D model is as follows:

(b1) Adjust A to approximate the gradient of the stress–strain curve

(b2) Adjust C to approximate the overall shape of each stress–strain curve

(b3) Adjust εs to approximate the peak load point and slope of the failure envelope

(b4) Adjust k to approximate the gradient of the axial stress–lateral strain curve and peak strength

(b5) Iterate steps (b1–b4) until the stress–strain curves for the rocks are approximated well

(b6) Adjust k to approximate the peak strength for the triaxial compression test

(b7) Use a nonlinear least squares method for the k vs. σ3 curve to obtain values of k0⁠, ∆k⁠, α⁠, and σ0 (e.g., Figures 14 and 15(b))

The stress–strain curves and convex failure envelopes for two types rocks, limestone and red sandstone [12, 13], are approximated by the modified 2D model. For limestone, the parameters A=27 GPa⁠, C=3.5⁠, and εs=0.0013 were given, and the parameters k0=0.143⁠, ∆k=0.0239⁠, α=2.04 MPa⁠, and σ0=20.9 MPa were obtained by a nonlinear least-squares method (Figure 14(b)). For red sandstone, the parameters A=9 GPa⁠, C=2.5⁠, and εs=0.002 were given, and the parameters k0=0.18⁠, ∆k=0.0928⁠, α=20.6 MPa⁠, and σ0=29.0 MPa were obtained by a nonlinear least-squares method (Figure 15(b)). The stress–strain curves (Figures 14(a) and 15(a)) and convex failure envelopes (Figures 14(c) and 15(c)) for rocks are simulated well.

3.3. Comparison with Other 2D Models

For comparison, we discuss the examples of two conventional 2D elasto-plastic constitutive models (Figure 16(a), [12]; Figure 16(b), [13]). These models are complicated and need 11 and 10 parameters, respectively, to approximate the experimental stress–strain curves.

The results according to Zhang [12], including elastic (⁠E0=45.51 GPa and ν0=0.320⁠), strength (⁠α1=2.07⁠, κ=34.5 MPa⁠, and α2=2.07⁠), plastic hardening (⁠h0=0.58⁠, h1=1.40⁠, and b=900⁠), damage (⁠ωc=0.70⁠, β1=0.50 MPa−1⁠, and r=0.2⁠), and damage softening (⁠a=3.45 exp−0.047σ3⁠) parameters, are represented by the blue lines in Figure 14(a).

The results according to Zhang, et al. [13], including elastic (bulk modulus, K0=8 GPa⁠; shear modulus initial value, G0=6 GPa⁠; shear modulus residual value, Gr=G0×0.8⁠) and plastic (friction angle plastic internal variable threshold, Pϕ=0.4⁠; friction angle initial value, ϕ0=18°⁠; friction angle residual value, ϕr=34.5°⁠; shear dilatancy angle initial value, ψ0=25°⁠; shear dilatancy angle residual value, ψr=10°⁠; cohesion initial value c0=15 MPa⁠; cohesion residual value cr=9.0 MPa⁠) parameters, and a plastic internal variable (⁠2.5σ3/σc+0.27⁠) are represented by the black lines in Figure 15(a).

The two models simulate the experimental stress–strain curves well. However, the models can neither simulate tensile failure nor consider the critical strains, and they both approximate failure envelopes as linear.

Among these three models, the modified 2D model is the best because it can approximate the stress–strain curves as precisely as the others, but can approximate convex failure envelope and uses fewer equations and parameters.

3.4. Example Application

As an example of its application, we calculated the stress distribution of a pressurized thick-walled cylinder with inner and outer radii r1=1 m and r2=100 m⁠, respectively, under internal and external pressures P1=1 MPa and P2=80 MPa⁠, respectively, using the simple 2D model and assuming the plane strain condition. The calculation was carried out by an axisymmetric finite element method (FEM), dividing the radial axis from 1 to 19 m into 470 elements, and from 19 to 100 m into 320 elements. The results were compared to Bray’s analytical elasto-plastic solution [50] for a circular hole in an infinite elasto-brittle material (Figure 1(c)). The modified 2D model was not used because Bray’s solution was only for a linear failure envelope and thus cannot be compared.

Parameters that were the same as for Figures 10(a) and 10(f) (⁠A=20 GPa⁠, k=0.25⁠, C=1⁠, and εs=0.002⁠) and were used for the simple 2D model. The uniaxial compressive strength qu=71.4 MPa⁠, internal friction angle ϕ=36.5°⁠, friction angle for the residual strength ϕ’=36.9°⁠, E=14.5 GPa⁠, ν=0.270⁠, and residual cohesion Sj=9.89 MPa were used for Bray’s solution. The parameters are obtained from Figures 10(a)–10(f).

Plastic zones can be seen for both models (Figure 17(a)). The decrease in strain-dependent elastic modulus (Figure 17(b)) around the opening induced a tangential stress decrease for the simple 2D model. The results from Bray’s solution appear somewhat angular when plotted, while the results of the simple 2D model are more rounded, reflecting the difference in the stress–strain relationship. The inward displacement by the simple 2D model (Figure 17(c)) rapidly increases near the inner wall. This rapid increase may represent the dilatancy of rock failure better than Bray’s solution, which is based on the rupture plane slip in the plastic zone.

We propose the following simple 3D constitutive model:
(21)σ1=Aε1ε1+kAε2+kAε3,(22)σ3=kAε1+kAε2+Aε3ε3,(23)ε2=σ1−σ2σ1−σ3ε3+σ2−σ3σ1−σ3ε1.
Equations (21) and (22) are based on Equations (14) and (15) for the 2D model. Equation (23) is just an elastic equation. The same strain-dependent elastic modulus A’ε is used as in Equation (16):
(24)A′ε=Dtan−1Cεεs+1+π2+F,
and the terms D and F are chosen so that the residual strength becomes constant for infinite axial compression and lateral expansion (similar to the 2D model) as follows.
Let us assume
(25)A′ε1=∞=A.
As ε3 tends to −∞⁠,
(26)A′ε3=−∞=F.
Substituting Equations (25) and (26) into Equations (21) and (22),
(27)σ3ε1=∞,ε3=−∞=Fε3+kAσ1−σ2σ1−σ3ε3+σ2−σ3σ1−σ3ε1+kAε1,(28)σ1ε1=∞,ε3=−∞=Aε1+kAσ1−σ2σ1−σ3ε3+σ2−σ3σ1−σ3ε1+kAε3.
Substituting Equation (27) into Equation (28),
(29)σ1ε1=∞,ε3=−∞=Hε1+kA1+σ1−σ2/σ1−σ3kAσ1−σ2/σ1−σ3+Fσ3,
where
(30)H=A+kAσ2−σ3σ1−σ3−k2A21+σ1−σ2/σ1−σ3kAσ1−σ2/σ1−σ3+F1+σ2−σ3σ1−σ3.
Assuming σ1ε1=∞,ε3=−∞=const.⁠, H should be 0,
(31)F=k1+σ1−σ2/σ1−σ31+σ2−σ3/σ1−σ3−σ1−σ2/σ1−σ31+kσ2−σ3/σ1−σ31+kσ1−σ2/σ1−σ3kA,(32)D=A−Fπ.

For example, the simulation procedure for stress–strain curves becomes

(c1) Assign negative ε3 increment

(c2) Calculate A’ε by Equation (24)

(c3) Calculate ε1 by Equation (22)

(c4) Calculate σ1 by Equation (21)

(c5) Iterate steps (c2–c4) until convergence

(c6) Calculate ε2 by Equation (23)

(c7) Iterate steps (c1–c6) until ε3 reaches -0.045

We approximated the true triaxial compression data under σ2=300 (MPa) and σ3=100 (MPa) for Cretaceous–Tertiary boundary (KTB) amphibolite [28] and under σ2=125 MPa and σ3=50 MPa for Coconino sandstone [29]. The approximation procedure is the same as steps (b1–b5). The results (Figure 18) is approximated well by the simple 3D model. This model requires only four parameters (⁠A=90 GPa⁠, C=2⁠, εs=0.0035⁠, and k=0.22 for Figure 18(a); A=37.5 GPa⁠, C=1.5⁠, εs=0.0025⁠, and k=0.217 for Figure 18(b)). However, the peak stress–σ2 relationship is not considered in this model; this is one of its limitations. We introduce stress-dependent k in the next chapter to approximate the convex peak stress–σ2 relationship.

5.1. Introduction of Stress-Dependent k

We now introduce stress-dependent k (Figure 19), with an additional term –βσ2 compared to Equation (20):
(33)k=Δkπtan−11ασ3−βσ2−σ0+π2+k0,
where β is a constant to adjust the effect of σ2⁠, taking a value between 0 (no effect) and 1 (the same effect as σ3⁠).

For example, the calculation procedure for stress–strain curves becomes

(d1) Assign negative ε3 increment

(d2) Calculate A’ε by Equation (24)

(d3) Calculate ε1 by Equation (22)

(d4) Calculate σ1 by Equation (21)

(d5) Calculate ε2 by Equation (23)

(d6) Calculate k by Equation (33)

(d7) Iterate steps (d2–d6) until k converges

(d8) Iterate steps (d1–d7) for loading until ε3 reaches -0.045

5.2. Prediction of σ2 Effects

In chapter 4, the true stress–strain curves of KTB amphibolite [28] and Coconino sandstone [29] were simulated well by the 3D constitutive model for specific values of σ2 and σ3⁠. We then approximated the stress–strain curves and peak stress–σ2 relationship for Westerly granite [51] and the two rocks above using the modified 3D model. The following steps are added to the approximation steps above (b1–b5):

(e6) Assume β at 0.3, 0.5, 0.7, and 0.9

(e7) Obtain k vs. σ3−βσ2 curve for σ2=σ3

(e8) Assume k0 (e.g., 0.075, 0.080, 0.085, 0.090, and 0.095 for Westerly granite)

(e9) Obtain values of ∆k⁠, α⁠, and σs by a nonlinear least-squares method (e.g., Figure 19)

(e10) Calculate stress–strain curves and failure envelopes

(e11) Iterate steps (e8–e10)

(e12) Iterate steps (e6–e11)

(e13) Select a set of parameters that can best simulate the peak stress–σ2 curve under σ2>σ3

An example plot of k vs. σ3−βσ2 for Westerly granite with β=0.7 is shown in Figure 19. The parameter values were obtained for each k0 by a nonlinear least-squares method. The failure envelopes of all rock types for σ2=σ3 (Figures 20(a), 21(a), and 22(a)) were approximated reasonably well. The true stress–strain curves of rock (Figures 20(c), 21(c), and 22(c)) under the modified 3D model were somewhat different from those under the simple 3D model (Figure 18), because the former prioritized approximating the failure envelope. For Westerly granite and KTB amphibolite, the failure envelope for σ2>σ3 at high σ3 was simulated well, while the modified 3D model underestimates the peak stress at low σ3 (Figures 20(b) and 21(b)). For sandstone, this model simulated the failure envelope well for σ2>σ3 at low σ3⁠, while it overestimates the peak stress at high σ3 (Figure 22(b)).

5.3. Comparison with Other 3D Criteria

In Figure 23, we show a comparison of 3D failure criteria ([36, 40, 41, 43, 52], and [38]) from Li et al. [38]. In this figure, we have also plotted the peak stress–σ2 relationship according to the modified 3D model in chapter 5.2. Some of these criteria give better results than the proposed model; however, they do not give a reasonable nonlinear stress–strain relationship. Therefore, it is better to use one of these criteria when the stress–strain relationship is almost linear, and it is important that the failure envelope be followed strictly. However, it is better to use our proposed model to represent nonlinear stress–strain relationships well while following the true triaxial failure envelope.

To propose a 3D constitutive model for rocks, with fewer parameters and consideration of true triaxial stress–strain curves and convex failure envelopes, we first modified the simple 2D model previously proposed by our research group to have a convex failure criterion. The effects of the parameters on the failure envelope were observed, and a procedure to approximate the stress–strain curves and failure envelopes was proposed. This model approximated these well (much better than conventional 2D models). The new 2D model was applied to a pressurized thick-walled cylinder under the plane strain condition.

The simple 2D constitutive model was then extended to a simple 3D model. True triaxial stress–strain curves under specific values of σ2 and σ3 were approximated reasonably. However, the predicted peak stress–σ2 relationship was linear.

Finally, the 3D model was developed by introducing stress-dependent k⁠. Under this, the peak stress–σ2 relationship for true triaxial compression first increased and then decreased with σ2⁠, as occurs in actual rocks. True triaxial stress–strain curves and failure envelopes were simulated reasonably.

This modified 3D model is the only one that can fairly represent the true triaxial nonlinear stress–strain curves and convex failure envelopes. Combining it with a 3D finite element method could improve the design of rock structures.

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

We declare there is no conflict of interest.

This research was supported by the Ministry of Education, Culture, Sports, Science, and Technology of Japan and the China Scholarship Council.