## Abstract

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 ($\sigma 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 $\sigma 2$ and $\sigma 3$ and uses only four parameters. However, the predicted peak stress–$\sigma 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.

## 1. Introduction

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.

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].

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.

No constitutive equations can fully represent true triaxial stress–strain curves considering the effects of $\sigma 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.

## 2. The Simple 2D Constitutive Model

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.

*ε*tends to $\u2212\u221e$ (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]).

However, the residual strength $\sigma 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 $\epsilon 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 $\epsilon 3$ increment (extensile) is assigned

(a2) $\epsilon 1$ is calculated by Equation (15)

(a3) $A\u2019\epsilon 1$ and $A\u2019\epsilon 3$ are calculated by Equation (16)

(a4) $\sigma 1$ can be calculated by Equation (14)

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

(a6) Iterate steps (a1–a5) for loading until $\epsilon 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\u2019\epsilon 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 $\sigma 1\u2013\epsilon 1$ curve around the peak load, Poisson’s ratio, and the critical compressive strain ($\epsilon 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)). $\epsilon 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:

## 3. Modified 2D Model

### 3.1. Introduction of $k$ Dependent on Confining Pressure

The stress–strain curves obtained by setting default values of $A=20\u2009GPa$, $C=1$, $\epsilon s=0.002$, $k0=0.20$, $\sigma 0=35\u2009MPa$, $\alpha =20\u2009MPa$, and $\u2206k=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 $\sigma c=71.6\u2009MPa$, $mb=119.0$, $s=1$, and $\alpha =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. $\alpha $ 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 $\epsilon 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. $\sigma 3$ curve to obtain values of $k0$, $\u2206k$, $\alpha $, and $\sigma 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\u2009GPa$, $C=3.5$, and $\epsilon s=0.0013$ were given, and the parameters $k0=0.143$, $\u2206k=0.0239$, $\alpha =2.04\u2009MPa$, and $\sigma 0=20.9\u2009MPa$ were obtained by a nonlinear least-squares method (Figure 14(b)). For red sandstone, the parameters $A=9\u2009GPa$, $C=2.5$, and $\epsilon s=0.002$ were given, and the parameters $k0=0.18$, $\u2206k=0.0928$, $\alpha =20.6\u2009MPa$, and $\sigma 0=29.0\u2009MPa$ 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\u2009GPa$ and $\nu 0=0.320$), strength ($\alpha 1=2.07$, $\kappa =34.5\u2009MPa$, and $\alpha 2=2.07$), plastic hardening ($h0=0.58$, $h1=1.40$, and $b=900$), damage ($\omega c=0.70$, $\beta 1=0.50\u2009MPa\u22121$, and $r=0.2$), and damage softening ($a=3.45\u2009exp\u22120.047\sigma 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\u2009GPa$; shear modulus initial value, $G0=6\u2009GPa$; shear modulus residual value, $Gr=G0\xd70.8$) and plastic (friction angle plastic internal variable threshold, $P\varphi =0.4$; friction angle initial value, $\varphi 0=18\xb0$; friction angle residual value, $\varphi r=34.5\xb0$; shear dilatancy angle initial value, $\psi 0=25\xb0$; shear dilatancy angle residual value, $\psi r=10\xb0$; cohesion initial value $c0=15\u2009MPa$; cohesion residual value $cr=9.0\u2009MPa$) parameters, and a plastic internal variable ($2.5\sigma 3/\sigma 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\u2009m$ and $r2=100\u2009m$, respectively, under internal and external pressures $P1=1\u2009MPa$ and $P2=80\u2009MPa$, 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\u2009GPa$, $k=0.25$, $C=1$, and $\epsilon s=0.002$) and were used for the simple 2D model. The uniaxial compressive strength $qu=71.4\u2009MPa$, internal friction angle $\varphi =36.5\xb0$, friction angle for the residual strength $\varphi \u2019=36.9\xb0$, $E=14.5\u2009GPa$, $\nu =0.270$, and residual cohesion $Sj=9.89\u2009MPa$ 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.

## 4. Simple 3D Model

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

(c1) Assign negative $\epsilon 3$ increment

(c2) Calculate $A\u2019\epsilon $ by Equation (24)

(c3) Calculate $\epsilon 1$ by Equation (22)

(c4) Calculate $\sigma 1$ by Equation (21)

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

(c6) Calculate $\epsilon 2$ by Equation (23)

(c7) Iterate steps (c1–c6) until $\epsilon 3$ reaches -0.045

We approximated the true triaxial compression data under $\sigma 2=300$ (MPa) and $\sigma 3=100$ (MPa) for Cretaceous–Tertiary boundary (KTB) amphibolite [28] and under $\sigma 2=125\u2009MPa$ and $\sigma 3=50\u2009MPa$ 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\u2009GPa$, $C=2$, $\epsilon s=0.0035$, and $k=0.22$ for Figure 18(a); $A=37.5\u2009GPa$, $C=1.5$, $\epsilon s=0.0025$, and $k=0.217$ for Figure 18(b)). However, the peak stress–$\sigma 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–$\sigma 2$ relationship.

## 5. Modified 3D Model

### 5.1. Introduction of Stress-Dependent $k$

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

(d1) Assign negative $\epsilon 3$ increment

(d2) Calculate $A\u2019\epsilon $ by Equation (24)

(d3) Calculate $\epsilon 1$ by Equation (22)

(d4) Calculate $\sigma 1$ by Equation (21)

(d5) Calculate $\epsilon 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 $\epsilon 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 $\sigma 2$ and $\sigma 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 $\beta $ at 0.3, 0.5, 0.7, and 0.9

(e7) Obtain $k$ vs. $\sigma 3\u2212\beta \sigma 2$ curve for $\sigma 2=\sigma 3$

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

(e9) Obtain values of $\u2206k$, $\alpha $, and $\sigma 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–$\sigma 2$ curve under $\sigma 2>\sigma 3$

An example plot of $k$ vs. $\sigma 3\u2212\beta \sigma 2$ for Westerly granite with $\beta =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 $\sigma 2=\sigma 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 $\sigma 2>\sigma 3$ at high $\sigma 3$ was simulated well, while the modified 3D model underestimates the peak stress at low $\sigma 3$ (Figures 20(b) and 21(b)). For sandstone, this model simulated the failure envelope well for $\sigma 2>\sigma 3$ at low $\sigma 3$, while it overestimates the peak stress at high $\sigma 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–$\sigma 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.

## 6. Concluding Remarks

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 $\sigma 2$ and $\sigma 3$ were approximated reasonably. However, the predicted peak stress–$\sigma 2$ relationship was linear.

Finally, the 3D model was developed by introducing stress-dependent $k$. Under this, the peak stress–$\sigma 2$ relationship for true triaxial compression first increased and then decreased with $\sigma 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.

## Data Availability

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

## Conflicts of Interest

We declare there is no conflict of interest.

## Acknowledgments

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