A Strain Hardening and Softening Constitutive Model for Hard Brittle Rocks

Abstract

To study the strain hardening and softening mechanism for hard brittle rocks, a strain hardening and softening constitutive model for hard brittle rocks is developed. First, the normalised hardening and softening factors are defined, which characterise the yield state of rock at the stages of pre-peak hardening and post-peak softening, respectively. Then, a unified strength parameter evolution model is established that can describe the nonlinear characteristics of cohesion and the internal friction angle under different confining pressures. Based on the Mohr–Coulomb criterion, a strain hardening and softening constitutive model is proposed. Finally, the proposed model was implemented in FLAC^3D, and triaxial compression numerical tests of granite and diabase were conducted. The results show that the constitutive model can characterise the nonlinear mechanical behaviour of the pre-peak hardening stage and post-peak softening stage of hard brittle rock. The model was also able to satisfactorily capture the transition from brittle failure to plastic failure for hard brittle rock under high confining pressures.

1. Introduction

With the advancement of technology, human society’s demand for deep underground rock engineering construction has increased [1,2]. Many deep underground structures are built in hard brittle rock masses such as granite and diabase [3,4]. The fracture or even failure of hard brittle rock is the most important challenge in engineering construction. Hard brittle rocks often exhibit obvious nonlinear mechanical properties such as strain hardening/softening and plastic dilatation under high in situ stress [5]. Therefore, the in-depth study of the nonlinear mechanical behaviours of hard brittle rocks and the establishment of a hardening and softening constitutive model is of great significance for engineering design and engineering safety assessment [6].

The phenomenon in which a significant increase or decrease in stress accompanies increasing strain after a rock has yielded under compression conditions is known as strain hardening or strain softening of rocks. In view of this nonlinear mechanical behaviour, many constitutive models have been proposed. Based on the traditional strength yield criterion, the damage constitutive models have been proposed [7]. These models can describe the nonlinear behaviour of the rock under loading but do not consider the true changes in the subsequent yield surface of the rock at the strain hardening and strain softening stages [8]. In recent years, some scholars have used plastic parameters as independent variables to characterise the strain hardening/softening behaviour of rocks and represent the evolution of subsequent yield surfaces and have proposed elastic-plastic strain hardening and softening constitutive models [9,10]. Because the Mohr–Coulomb yield criterion can satisfactorily describe the shear failure mechanism, it has been widely adopted as the yield surface equation. The following two strength parameters are needed to determine the evolution law of the yield surface: cohesion and the internal friction angle [11]. Yielding occurs once the shear stress of the rock reaches a stress threshold.

Based on the Mohr–Coulomb yield criterion, various strength parameter evolution models are proposed. According to the strength parameters at peak stress and residual stress, the strength parameters decrease linearly with the increase in the plastic parameter [12,13]. Based on the degradation index concept proposed by Fang and Harrison [14], the cohesion and internal friction angle strength degradation index was proposed to describe the process of linear reduction of the post-peak strength parameters [15], and Alejano et al. [16] proposed three strain softening models by adopting the drop modulus variation. Lee and Pietruszczak [17] and Fahimifar and Zareifard [18] obtained a finite strain solution for circular underground tunnels in a Mohr–Coulomb elastoplastic coupling strain softening rock mass, assuming that cohesion and the internal friction angle decrease linearly with increasing plastic parameters.

The application of the linear reduction model of post-peak strength parameters is relatively simple, but the actual softening process of the post-peak stage is ignored. To accurately describe the post-peak mechanical behaviour of rock, Ranjbarnia et al. [19] used multiple piecewise linear functions to approximate the nonlinear softening behaviour after the peak stress was reached. Li et al. [20] proposed that the decreasing law of the internal friction angle conforms to a quadratic function and that cohesion decreases linearly after the peak stress is reached. Based on the triaxial test data of sandstone, mudstone and argillaceous sandstone, O. Pourhosseini et al. [21] proposed a strength parameter evolution model of constant friction and cohesion weakening, assuming the internal friction angle remains constant during the post-peak softening process. The studies mentioned above consider the peak point as the starting point of the softening process, simplify the pre-peak process to one of linear elasticity, ignore the pre-peak rock hardening process and only study the post-peak strength parameter evolution law [22].

To overcome the abovementioned shortcomings, relevant scholars took the initial yield point of rock as the starting point and established strength parameter evolution models to describe the strain hardening and softening process. V. Hajiabdolmajid et al. [23] proposed a linear CWFS model (cohesion weakening-friction strengthening) based on experimental data derived from granite and implemented this model in FLAC2D numerical software to analyse the “V” damage caused by tunnel excavation. Walton [24] and Rafiei and Martin [25] proposed a nonlinear CWFS model based on the results of triaxial compression tests performed on granite under different confining stresses, and their model curves were in good agreement with the test curves. Li et al. [26] analysed the results of triaxial cyclic loading and unloading of Beishan granite and found that the cohesion gradually decreased as an exponential function with the plastic parameter increasing after damage occurred in the rock, while the internal friction angle increased as a log-normal distribution function until the peak was reached and then gradually decreased. On the basis of Li’s [26] research, Min et al. [27] proposed a strain softening constitutive model considering the elastic modulus degradation effect. Meng et al. [28] revealed the relationship between the softening parameter, cohesion and the internal friction angle through triaxial cyclic loading and unloading tests. Wang et al. [29] established a two-dimensional strength parameter model considering the confining pressure and plastic parameters and revealed the relationship between the parameters, cohesion and the internal friction angle.

It is worth noting that by solving for cohesion and the internal friction angle based on the Mohr–Coulomb yield criterion [20,26,27,28,29], the rock was assumed to experience the same plastic deformation when reaching peak stress and residual stress under different confining stresses, and on this basis, the relationship between confining pressure and stress under the same plastic deformation was established, and then the cohesion and internal friction angle evolution was obtained. However, the experimental results show that the plastic deformation of the rock differs under different confining pressures when the stress value of the yield state (such as peak stress and residual stress) remains constant [30]. Therefore, the strength parameter evolution law that is based on the same plastic deformation will produce results inconsistent with the experimental results.

In summary, on the one hand, there are few methods to determine the stress value of the same yield state under different confining pressures; on the other hand, the current model cannot simultaneously calculate the cohesion and internal friction angle at the initial yield stress, peak stress and residual stress and reflect the evolution law of the strain hardening and softening process. Therefore, the main aim of this paper is to develop a strain hardening and softening constitutive model that can consider the nonlinear characteristics of rocks under different confining pressures. The constitutive model proposed in this paper differs from other models based on the following three characteristics:

(i) The normalised hardening factor and softening factor are defined to characterise the yield state of the rock at the initial yield stress point, peak stress point, residual stress point and any yield stress points.

(ii) A unified strength parameter evolution model is proposed, and it can reflect the variations in cohesion and the internal friction angle in the subsequent yield process.

(iii) The strain hardening and softening constitutive model can simultaneously capture the mechanical nonlinear behaviours of hard brittle rocks in the pre-peak hardening and post-peak softening stages under different confining pressures. The results show that the nonlinear characteristics of stress-strain curves can be reproduced well under different confining pressures.

In the following sections, the hardening and softening factors are defined based on the simplified stress-strain curve. Then, the relationship between the factors, the strength parameters and the axial plastic strain was analysed, and a unified strength parameter evolution model considering confining pressures was proposed. Then, the evolution model is implemented in FLAC^3D to realise the secondary development of the strain hardening and softening constitutive model. Finally, numerical triaxial compression tests are conducted on granite and diabase, and the numerical results are verified with the test results.

2. Simplification of the Complete Stress-Strain Curve

To analyse the strain hardening and softening characteristics of hard brittle rock before and after the peak, a reasonable simplification of the complete stress-strain curve for hard brittle rocks is needed. The characteristic stress at the reversal point of the volume–axial strain curve is taken as the initial yield point [31], as shown in Figure 1a, which is denoted by the initial yield stress σ_s. The initial yield stress σ_s, peak stress σ_p, and residual stress σ_r are regarded as segmentation points, and the typical complete stress-strain curve was idealised as a four-stage simplified complete stress-strain curve, which is a prerequisite for establishing a strain hardening and softening constitutive model. The simplified complete stress-strain curve covers the linear elastic stage, pre-peak hardening stage, post-peak softening stage and residual strength stage [32]. The simplified curve is shown in Figure 1b. Note that the OA-AB-BC segment of the curve in Figure 1a is simplified by a straight line (OC).

In addition, the following statement is made in this study:

The Mohr–Coulomb yield criterion, which has been widely adopted to describe the shear failure mechanism of rock [33,34,35,36], and based on this criterion, the strain hardening and softening constitutive model is studied. That is,

$${\sigma }_1=\frac{1+\sin \phi }{1-\sin \phi }{\sigma }_3+\frac{2c\cdot \cos \phi }{1-\sin \phi }$$

(1)

where σ₁ is the maximum principal stress; σ₃ is the minimum principal stress; c is cohesion; and φ is the internal friction angle.

3. Mechanical Behaviour of Brittle Hard Rock

3.1. The Test Data of Granite and Diabase

The test data of granite and diabase are obtained from published references [26] and [37], respectively. The granite samples [26] were collected from the borehole of BS19 in Beishan of Gansu province, China, with a sampling depth of 450–455 m. Beishan granite has good uniformity and compactness. The mineral composition of the rock is analysed by an optical microscope. Its main components are quartz, alkali feldspar, plagioclase, and biotite, and the rock is named monzogranite. The average dry density is 2690 kg/m³, and the longitudinal wave velocity is 5019 m/s. According to the standard requirements of the ISRM, the samples were prepared into standard specimens with a diameter of 50 mm and a high of 100 mm. The experiments were carried out on the dry rock sample at room temperature. The triaxial compression test was carried out on the MTS-815 rock mechanics test machine, and six series of tests at confining pressures of 0.2 MPa, 1 MPa, 3.5 MPa, 10 MPa, 20 MPa and 40 MPa were performed.

The diabase rock samples [37] were collected from the fresh rock mass of Danjiangkou Reservoir in Central China. The dry density is 2970 kg/m³, and the average porosity is 0.37%. Meanwhile, the water content is 0.059%, and the sample dimensions are Φ = 50 mm × h = 100 mm. All parameters were obtained by standard methods recommended in ISRM. Triaxial compression tests were carried out on an MTS testing machine, and four series of tests at confining pressures of 5 MPa, 10 MPa, 20 MPa and 40 MPa were performed, and all the experiments were carried out on the dry rock sample at room temperature.

3.2. Young’s Modulus and Poisson’s Ratio

Young’s modulus E and Poisson’s ratio μ of the rock under different confining pressures can be formulated by Equations (2) and (3) according to the linear elastic stage of the simplified stress-strain curve.

$$E=\frac{{\sigma }_s}{{\epsilon }_{1s}}$$

(2)

$$\mu =\left|\frac{{\epsilon }_{3s}}{{\epsilon }_{1s}}\right|$$

(3)

where σ_s is the initial yield stress of the rock and ε_1s and ε_3s are the axial and radial strains of the initial yield stress of the rock.

For most hard brittle rocks such as granite and diabase, the elastic stage is less affected by the confining pressure. Therefore, the average Young’s modulus E_m and the average Poisson’s ratio μ_m are usually used to characterise Young’s modulus and Poisson’s ratio of rock under various confining pressures [38], and the average values are shown in

Table 1. Average values of Young’s modulus and Poisson’s ratio for Beishan granite [26] and Danjiangkou diabase [37].

Rock Type	Em/GPa	μm
Granite	65.47	0.26
Diabase	50.28	0.24

.

3.3. Plastic Parameter

In FLAC^3D, the equivalent plastic strain increment Δε^ps is used as the plastic parameter to describe the evolution process of the strength parameters in the softening regime and is defined as follows:

$$\Delta {\epsilon }^{ps}=\frac{1}{\sqrt{2}}\sqrt{{\left(\Delta {\epsilon }_1^{ps}-\Delta {\epsilon }_m^{ps}\right)}^2+{\left(\Delta {\epsilon }_m^{ps}\right)}^2+{\left(\Delta {\epsilon }_3^{ps}-\Delta {\epsilon }_m^{ps}\right)}^2}$$

(4)

where Δ${\epsilon }_1^{ps}$ and Δ${\epsilon }_3^{ps}$ are the maximum and minimum principal plastic strain increments, and

$$\Delta {\epsilon }_m^{ps}=\frac{1}{3}\left(\Delta {\epsilon }_1^{ps}+\Delta {\epsilon }_3^{ps}\right)$$

(5)

The plastic potential function g^s corresponds to the noncorrelated flow rule, which can be expressed as follows:

$$g^s={\sigma }_1-{\sigma }_3{N}_{\psi }$$

(6)

where

$$N_{\psi }=\frac{1+\sin \psi }{1-\sin \psi }$$

(7)

and ψ is the dilation angle of the rock. The principal plastic strain increment Δ${\epsilon }_i^{ps}$ can be calculated as follows:

$$\Delta {\epsilon }_i^{ps}={\lambda }^s\frac{\partial g^s}{\partial {\sigma }_i},i=1,2,3$$

(8)

where λ^s is the plastic multiplier, which cannot be calculated in this process. Substituting Equation (6) into Equation (26), the principal plastic strain increment Δ${\epsilon }_i^{ps}$ can be written as follows:

$$\Delta {\epsilon }_1^{ps}={\lambda }^s$$

(9)

$$\Delta {\epsilon }_3^{ps}=-{\lambda }^s{N}_{\psi }$$

(10)

$$\Delta {\epsilon }_3^{ps}=-N_{\psi }\Delta {\epsilon }_1^{ps}$$

(11)

Using Equation (10) above, Equation (4) becomes the following:

$$\Delta {\epsilon }^{ps}=\frac{\sqrt{3}}{3}\sqrt{1+N_{\psi }+N_{\psi }^2}\Delta {\epsilon }_1^{ps}$$

(12)

If dilation angle ψ is constant, N_ψ is also constant, and relation between ε^ps and ${\epsilon }_1^{ps}$ can be established in the following form:

$${\epsilon }^{ps}=\frac{\sqrt{3}}{3}\sqrt{1+N_{\psi }+N_{\psi }^2}\sum \Delta {\epsilon }_1^{ps}=\frac{\sqrt{3}}{3}\sqrt{1+N_{\psi }+N_{\psi }^2}{\epsilon }_1^{ps}$$

(13)

According to Equation (13), the axial plastic strain ${\epsilon }_1^{ps}$ can be used as a plastic parameter in FlAC^3D to control the strength parameter evolution process in the strain hardening and softening stages [31].

3.4. Normalised Hardening and Softening Factors

To obtain the stress value with the same yield state under each confining pressure and then discuss the relationship between the strength parameters and the axial plastic strain, the normalised hardening factor and softening factor are proposed to characterise a certain yield state corresponding to a certain stress value in the strain hardening and softening stages under different confining pressures. The initial yield stress, peak stress and residual stress are regarded as the stress values corresponding to the initial yield state, peak yield state and residual yield state, respectively. Then, a certain stress value of one yield point under each confining pressure can also be regarded as the stress value corresponding to a certain yield state. According to reference [11,39], the normalised hardening factor and softening factor are defined to describe the yield state.

The hardening factor D_h is defined as the ratio of the difference between certain axial stress and the initial yield stress to the difference between the peak stress and the initial yield stress in the strain hardening process, i.e.,

$$D_h=\frac{{\sigma }_h-{\sigma }_s}{{\sigma }_p-{\sigma }_s}$$

(14)

where σ_p is the peak stress and σ_h is a certain axial stress of the strain hardening stage, σ_s < σ_h < σ_p.

The softening factor D_ss is defined as the ratio of the difference between the peak stress and a certain axial stress and the difference between the peak stress and residual stress in the strain softening process, i.e.,

$$D_{ss}=\frac{{\sigma }_p-{\sigma }_{\mathrm{ss}}}{{\sigma }_p-{\sigma }_r}$$

(15)

where σ_r is the residual stress and σ_ss is a certain axial stress of the strain softening stage, σ_r < σ_ss < σ_p.

When D_h = 0, the rock reaches the initial yield stress σ_s, plastic yielding occurs, and subsequently, the rock enters the strain hardening stage. The hardening factor increases with increasing axial stress. When the rock reaches the peak stress, the strain hardening stage ends and then enters the strain softening stage. The peak stress point is both the endpoint of the strain hardening stage and the starting point of the strain softening stage, that is, D_h = 1 and D_ss = 0. As the post-peak stress gradually decreases, the softening factor gradually increases until it reaches 1, that is, D_ss = 1. This indicates the end of the strain softening stage, and then the rock enters the residual stress stage, as shown in Figure 2. If the initial yield, peak and residual stress under each confining pressure are determined, the stress value corresponding to any hardening and softening factor under each confining pressure can be calculated accurately by Equations (14) and (15).

3.5. Strength Parameter Evolution Model

After the rock has reached the initial yield stress, plastic yielding begins to occur. Considering the axial strain ε_1s at the initial yield stress as the zero point of the plastic parameter, the strength parameter and plastic parameter values under the same yield state are determined by the normalised hardening and softening factors based on the Mohr–Coulomb yield criterion. Then, a cohesion and internal friction angle evolution model under different confining pressures with respect to the axial plastic strain ${\epsilon }_1^{ps}$ is established. The modelling process is as follows:

(1) Establish the corresponding relationship between the hardening and softening factors and the axial plastic strain under different confining pressures.

Since the hardening and softening factors are defined by the ratio of the stress difference, it is appropriate to consider more points near the peak and residual stress where the stress changes less so that the points are more uniformly distributed on the stress curve. In this way, the evolution law of the strength parameters and axial plastic strain can be discussed better. Considering that less plastic deformation occurs in the strain hardening process than in the strain softening process, five hardening factors D_h = 0, 0.2, 0.5, 0.8, 1.0 and eight softening factors D_ss = 0, 0.1, 0.2, 0.4, 0.6, 0.8, 0.9, 1.0 are selected and substituted into Equations (10) and (11), respectively, to obtain the twelve axial stresses under each confining pressure. The initial yield stress is taken as the starting point of the axial stress, and the relationship curve between the axial stress and axial plastic strain is established. The linear interpolation method is used for 12 axial stress values, and the data regarding the relationship between each hardening and softening factor and axial plastic strain under different confining pressures are approximately obtained. Beishan granite [26] data are shown as an example in Figure 3.

(2) Solve the cohesion c and internal friction angle φ corresponding to each hardening and softening factor under different confining pressures.

Based on the Mohr–Coulomb yield criterion, the axial stress and confining pressure of rock under different hardening and softening factors are substituted into the following equation as the maximum principal stress σ₁ and the maximum principal stress σ₃, respectively.

$${\sigma }_1=a{\sigma }_3+b$$

(16)

where a is the slope and b is the intercept.

The coefficients a and b obtained by linear fitting can be substituted into Equations (17) and (18), as follows:

$$\phi =\arcsin \left(\frac{a-1}{a+1}\right)$$

(17)

$$c=\frac{b}{2\sqrt{a}}$$

(18)

The cohesion c and internal friction angle φ of the rock under each hardening and softening factor can be calculated, as given in

Table 2. The cohesion c and internal friction angle φ values corresponding to each hardening and softening factor (Beishan granite [26] and Danjiangkou diabase [37]).

Dh/Dss	Granite			Diabase
	c/MPa	φ/°	R2	c/MPa	φ/°	R2
Dh = 0	27.29	41.44	0.94	43.91	28.54	0.97
Dh = 0.2	27.86	42.61	0.94	45.13	29.67	0.97
Dh = 0.5	28.69	44.20	0.94	46.91	31.26	0.96
Dh = 0.8	29.50	45.65	0.94	48.61	32.73	0.95
Dh = 1.0 (Dss = 0)	30.03	46.54	0.94	49.72	33.64	0.94
Dss = 0.1	28.00	45.70	0.95	48.24	32.29	0.95
Dss = 0.2	25.90	44.82	0.95	46.71	30.84	0.95
Dss = 0.4	21.44	42.87	0.97	43.5	27.57	0.96
Dss = 0.5	16.56	40.65	0.98	40.03	23.68	0.98
Dss = 0.8	11.17	38.07	0.99	36.23	18.94	0.99
Dss = 0.9	8.24	36.63	1.00	34.18	16.16	1.00
Dss = 1.0	5.12	35.05	1.00	31.99	13.00	1.00

.

As Table 2 shows, the linear fitting results for axial stress and confining pressure under the same hardening and softening factors are good, and the fitting correlation coefficients are greater than 0.94, which conforms to the shear failure law of rock described by the Mohr–Coulomb yield criterion. It also shows that the hardening factor and softening factor can determine the stress values under the same yield state and different confining pressures.

(3) Strength parameter evolution model considering axial plastic strain and confining pressure.

After the above two steps, the one-to-one relationship between the strength parameters and the axial plastic strain for the same hardening and softening factors under different confining pressures can be obtained, so the relationship between the cohesion c and internal friction angle φ of two kinds of hard brittle rocks and the axial plastic strain considering the confining pressure can be captured. Figure 4a–d show the strength parameter relationship curves with both axial plastic strain and confining pressure for the granite and diabase samples.

As shown in Figure 4, the cohesion and internal friction angle evolution laws for granite and diabase are very similar under different confining pressures. On the one hand, with the increase in the axial plastic strain, the cohesion and internal friction angle values first increase and then decrease; on the other hand, with the increase in the confining pressure, the axial plastic strain at the peak and residual stress increase gradually. The reason is that the brittle failure of rock gradually transitions to plastic failure with increasing confining pressure. As far as one confining pressure is concerned, the cohesion and internal friction angle of both hard brittle rocks increase rapidly in the pre-peak hardening stage, decrease gradually in the post-peak strain softening stage and gradually become constant with increasing axial plastic strain. The reason is as follows: when the rock reaches the yield stress, unstable crack propagation occurs. Because the crack surface is rough and does not form a penetrating failure surface, the friction force and biting force between the cracks increase. At this time, the rock has not been damaged but continues to harden with the propagation of cracks [40]. As the load increases, more cracks form, grow, connect and slip in the rock sample; therefore, the cohesive strength and occlusal force between the cracks gradually decrease and finally tend to become stable [41].

According to the work of reference [42], a smooth and continuous Gaussian function can be adopted to describe the evolution law of the strength parameters of granite and diabase with respect to the axial plastic strain, as follows:

$$G\left({\epsilon }_1^{ps}\right)=\left(c_p-c_r\right)\exp \left[-{\left(\frac{{\epsilon }_1^{ps}-{\epsilon }_{1,p}^{ps}}{{\xi }_c}\right)}^2\right]+c_r$$

(19)

$$F\left({\epsilon }_1^{ps}\right)=\left({\phi }_p-{\phi }_r\right)\exp \left[-{\left(\frac{{\epsilon }_1^{ps}-{\epsilon }_{1,p}^{ps}}{{\xi }_{\phi }}\right)}^2\right]+{\phi }_r$$

(20)

where c_p and φ_p are the cohesion and internal friction angle at the peak stress of the rock and c_r and φ_r are the cohesion and internal friction angle at the residual stress of the rock, respectively. ${\epsilon }_{1,p}^{ps}$ is the axial plastic strain at the peak stress of the rock under each confining pressure. ξ_c and ξ_φ are the cohesion and internal friction angle fitting parameters, which control the width of the model curve.

Although the Gaussian function can adequately describe the evolution process of the strength parameters in the strain softening and residual stages, there is still a considerable difference in the description of the evolution process of the strength parameters in the hardening stage. As shown in Figure 5, the difference between the value calculated by the Gaussian function and the experimental value decreases continuously with increasing axial plastic strain at the pre-peak hardening stage, and the difference value is zero at the peak stress. Therefore, Boltzmann’s sigmoid mathematical function is introduced to properly modify the Gaussian function to eliminate the difference in the Gaussian function at the strain hardening stage.

After the Gaussian function was modified, a unified strength parameter evolution model of cohesion/internal friction angle–axial plastic strain for two hard brittle rocks under different confining pressures can be expressed as follows:

$$c\left({\epsilon }_1^{ps}\right)=\left(c_p-c_r\right)\exp \left[-{\left(\frac{{\epsilon }_1^{ps}-{\epsilon }_{1,p}^{ps}}{{\xi }_c}\right)}^2\right]-\frac{2\Delta X_c}{1+\exp {\left(\left|\Delta X_c\right|\cdot {\epsilon }_1^{ps}/{\epsilon }_{1,p}^{ps}\right)}^{0.9}}+c_r$$

(21)

$$\begin{array}{c}\phi \left({\epsilon }_1^{ps}\right)\end{array}=\begin{array}{c}\left({\phi }_p-{\phi }_r\right)\exp \left[-{\left(\frac{{\epsilon }_1^{ps}-{\epsilon }_{1,p}^{ps}}{{\xi }_{\phi }}\right)}^2\right]-\frac{2\Delta X_{\phi }}{1+\exp {\left(\left|\Delta X_{\phi }\right|\cdot {\epsilon }_1^{ps}/{\epsilon }_{1,p}^{ps}\right)}^{0.9}}+{\phi }_r\end{array}$$

(22)

where:

$$\Delta X_c=G\left(0\right)-c_s$$

(23)

$$\Delta X_{\phi }=F\left(0\right)-{\phi }_s$$

(24)

where c_s and φ_s are the cohesion and friction angle at the initial yield stress, respectively. G(0) and F(0) represent the initial cohesion and internal friction angle values calculated by Gaussian Equations (19) and (20). ΔX_c and ΔX_φ are temporary parameters, which represent the differences between the initial values calculated by Equations (19) and (20) and the experimental values.

The modified strength parameter evolution model has the following two advantages:

(i) It has fewer fitting parameters and a unified function form. The fitting parameters ξ_c and ξ_φ are obtained by Equations (19) and (20), respectively, and the temporary parameters ΔX_c and ΔX_φ can be calculated by entering the fitting parameters ξ_c and ξ_φ into Equations (23) and (24), respectively. The peak parameter ${\epsilon }_{1,p}^{ps}$ can be obtained from the axial plastic strain at the peak stress. c_s and φ_s can be obtained by fitting the initial yield stress at each confining pressure. c_p and φ_p can be obtained by fitting the peak stress at each confining pressure. c_r and φ_r can be obtained by fitting the residual stress at each confining pressure. All fit coefficients for the cohesion/internal friction angle–axial plastic strain curves of granite and diabase are shown in

Table 3. Fit coefficients for cohesion/internal friction angle–axial plastic strain curves of two hard brittle rock types at different confining pressures.

Rock Type	σ3/MPa	${\mathit{\epsilon }}_{1,\mathit{p}}^{\mathit{p}\mathit{s}}$/%	c					φ
			ξc	cs/MPa	cp/MPa	cr/MPa	R2	ξφ	φs/°	φp/°	φr/°	R2
Granite	0.2	0.003	0.223	27.29	30.03	5.12	0.98	0.200	41.44	46.54	35.05	0.91
	1	0.057	0.282	27.29	30.03	5.12	0.99	0.282	41.44	46.54	35.05	0.98
	3.5	0.078	0.467	27.29	30.03	5.12	0.99	0.462	41.44	46.54	35.05	0.97
	10	0.120	0.379	27.29	30.03	5.12	0.99	0.364	41.44	46.54	35.05	0.98
	20	0.149	0.525	27.29	30.03	5.12	1.00	0.517	41.44	46.54	35.05	0.99
	40	0.239	0.672	27.29	30.03	5.12	0.99	0.663	41.44	46.54	35.05	0.98
Diabase	5	0.031	0.154	43.91	49.72	32.00	0.96	0.166	28.54	33.64	13.00	0.98
	10	0.060	0.212	43.91	49.72	32.00	0.99	0.214	28.54	33.64	13.00	0.98
	20	0.085	0.258	43.91	49.72	32.00	0.97	0.281	28.54	33.64	13.00	0.97
	40	0.108	0.283	43.91	49.72	32.00	0.99	0.305	28.54	33.64	13.00	0.97

.

(ii) The model can smoothly simulate the nonlinear cohesion and internal friction angle processes in the pre-peak hardening stage and the post-peak softening stage, as shown in Figure 4a–d. Compared with the Gaussian function before correction, the modified model can compensate for the disadvantage that the strength parameter evolution process in the pre-peak hardening stage cannot be accurately described, as shown in Figure 5. The fitting correlation coefficients of the strength parameter evolution model are mostly distributed in the range from 0.96 to 1.00 (Table 3), which illustrates that the established strength parameter evolution model is correct and reasonable.

As Table 3 shows, the peak parameter ${\epsilon }_{1,p}^{ps}$ and the fitting parameters ξ_c and ξ_φ gradually increase with increasing confining pressure. The reason is that the rock has more obvious plastic characteristics under high confining pressures, and more plastic deformation occurs after reaching the peak and residual stresses, so the peak parameter ${\epsilon }_{1,p}^{ps}$ becomes larger and the stress-strain curve appears wider, and the fitting parameters ξ_c and ξ_φ, which control the width of the model, also increase gradually with increasing confining pressure. The confining pressure dependence of ${\epsilon }_{1,p}^{ps}$, ξ_c and ξ_φ can be observed in Figure 6. The variations in these parameters with confining pressures can be described by a unified empirical function, as follows:

$$P\left({\sigma }_3\right)=P_0-A_1{e}^{-{\sigma }_3/t_1}-A_2{e}^{-{\sigma }_3/t_2}$$

(25)

where P represents one of the parameters ${\epsilon }_{1,p}^{ps}$, ξ_c or ξ_φ. P₀, A₁, A₂, t₁, and t₂ are fitting coefficients, and all values are positive. In addition, P₀ represents the range of the confining pressure function; A₁ and A₂ jointly control the increase in the amplitude of the confining pressure function; and t₁ and t₂ jointly control the growth rate of the confining pressure function. The fitting coefficients of each parameter are shown in

Table 4. Fitting coefficient of parameters ${\epsilon }_{1,p}^{ps}$, ξ_c and ξ_φ of two hard brittle rocks with the confining pressure function.

Parameter	Granite						Diabse
	P0	A1	A2	t1	t2	R2	P0	A1	A2	t1	t2	R2
${\epsilon }_{1,p}^{ps}$	1.510	0.075	1.445	0.565	321.094	0.99	0.114	0.119	0	13.422	-	1.00
ξc	6.790	0.198	6.414	1.027	874.929	0.92	0.285	0.227	0	8.947	-	0.99
ξφ	7.075	0.203	6.709	1.006	914.635	0.91	0.313	0.237	0	10.757	-	0.99

, and their fitting curves are shown in Figure 6. The fitting correlation coefficients of each parameter are all greater than 0.9, which indicates that the confining pressure function Equation (25) can satisfactorily describe the relationship between parameters ${\epsilon }_{1,p}^{ps}$, ξ_c and ξ_φ with various confining pressures.

According to Equations (21)−(25), an overall fit for granite and diabase is created to illustrate the strength parameter evolution model that considers both axial plastic strain and confining pressure simultaneously. As shown in Figure 7a–d, the results of the overall fit using the strength parameter evolution model are in good agreement with the individual fit results using the experimental data. With increasing confining pressure, the strength parameters first increase and then decrease, accurately describing the initial, peak and residual strength values. The strength parameters at the peak points gradually shift to the right, and the model curve tends to increase in width, which means that the failure mode of hard brittle rock transitions from brittle failure to plastic failure under high confining pressures.

4. Establishment and Verification of the Constitutive Model

4.1. Establishment of the Constitutive Model

In FLAC^3D, the strain hardening and softening constitutive model considering confining pressures is coded by the C++ language. The new stress vector σ^N of the current step can be calculated from the stress vector σ^p of the previous step and the strain increment vector Δε of the current step.

First, the elastic guess stress ${\sigma }_i^I$ is defined as follows:

$${\sigma }_i^I={\sigma }_i^p+S_i\Delta {\epsilon }_i,i=1,2,3$$

(26)

$$S_i=\left[\begin{array}{ccc}K+\frac{4}{3}G& K-\frac{2}{3}G& K-\frac{2}{3}G\\ K-\frac{2}{3}G& K+\frac{4}{3}G& K-\frac{2}{3}G\\ K-\frac{2}{3}G& K-\frac{2}{3}G& K+\frac{4}{3}G\end{array}\right]$$

(27)

where S_i is the elastic stiffness matrix, K is the bulk modulus, and G is the shear modulus.

Δε_i is the total principal strain increment, which can be calculated as the sum of the principal elastic strain increment Δ${\epsilon }_i^e$ and the principal plastic strain increment Δ${\epsilon }_i^{ps}$, as follows:

$$\Delta {\epsilon }_i=\Delta {\epsilon }_i^e+\Delta {\epsilon }_i^{ps}$$

(28)

For the elastic part, the principal elastic strain increment Δ${\epsilon }_i^e$ and the principal stress increment Δσ_i follow Hooke’s law, as follows:

$$\Delta {\sigma }_i=S_i\left(\Delta {\epsilon }_i^e\right),i=1,2,3$$

(29)

Second, judge whether the simulation element is yielded. The maximum principal guess stress ${\sigma }_1^I$ and the minimum principal guess stress ${\sigma }_3^N$ are substituted into the yield criterion to judge whether the element is yielded. In FLAC^3D, the Mohr–Coulomb yield criterion is given as follows:

$$f^s={\sigma }_1-{\sigma }_3\frac{1+\sin \phi }{1-\sin \phi }+2c\sqrt{\frac{1+\sin \phi }{1-\sin \phi }}$$

(30)

Third, if the element is yielded, then its stress should be updated. Substituting Equations (28) and (29) into Equation (8), the principal stress increment Δσ_i would be recalculated as follows:

$$\Delta {\sigma }_i=S_i\left(\Delta {\epsilon }_i^e\right)=S_i\left(\Delta {\epsilon }_i\right)-{\lambda }^s{S}_i\left(\frac{\partial g^s}{\partial {\sigma }_i}\right),i=1,2,3$$

(31)

The principal stress ${\sigma }_i^N$ of the current step can be expressed as follows:

$${\sigma }_i^N={\sigma }_i^p+\Delta {\sigma }_i,i=1,2,3$$

(32)

Substituting Equations (26) and (31) into Equation (32), the principal stress ${\sigma }_i^N$ of the current step can be written as follows:

$${\sigma }_i^N={\sigma }_i^I-{\lambda }^s{S}_i\left(\frac{\partial g^s}{\partial {\sigma }_i}\right),i=1,2,3$$

(33)

In FLAC^3D, when the element yields, the position of the principal stress point of the current step is on the yield surface, which can be expressed as follows:

$$f^s\left({\sigma }_1^N,{\sigma }_3^N\right)=0$$

(34)

where ${\sigma }_i^N$ and ${\sigma }_3^N$ are the maximum and minimum principal stresses of the current step, respectively. Substituting Equations (8), (26)–(30) and (33) into Equation (34), a formula for plastic multiplier λ^s can be written as follows:

$${\lambda }^s=\frac{f^s\left({\sigma }_1^I,{\sigma }_3^I\right)}{\left({\alpha }_1-{\alpha }_2{N}_{\psi }\right)-\left(-{\alpha }_1{N}_{\psi }+{\alpha }_2\right)N_{\psi }}$$

(35)

$${\alpha }_1=K+\frac{4}{3}G$$

(36)

$${\alpha }_2=K+\frac{2}{3}G$$

(37)

According to Equation (12), the maximum principal plastic strain increment $\Delta {\epsilon }_1^{ps}$ can be calculated. Meanwhile, the maximum principal plastic strain/axial plastic strain ${\epsilon }_1^{ps}$ can also be obtained as follows:

$${\epsilon }_1^{ps}={\epsilon }_1^{p,ps}+\Delta {\epsilon }_1^{ps}$$

(38)

$$\Delta {\epsilon }_1^{ps}=\frac{\sqrt{3}\Delta {\epsilon }^{ps}}{\sqrt{1+N_{\psi }+N_{\psi }^2}}$$

(39)

where ${\epsilon }_1^{p,ps}$ is the accumulated maximum principal plastic strain of the previous step. Substituting the accumulated maximum principal plastic strain ${\epsilon }_1^{ps}$ into Equations (21) and (22), we can update the strength parameter value and proceed to the next loop, as shown in Figure 8. To study the effect of the strength parameter evolution law of the strain hardening and softening constitutive model and to simplify the numerical calculation, a common consideration of the dilation angle results can be calculated by equation ψ = φ − 20° [43]. For the employed granite and diabase, the internal friction angle at peak stress can be used for calculation, and the dilation angle ψ = 26° and ψ = 13° is considered, respectively [27].

Finally, the prepared program file can be transformed into a dynamic link library file (.dll file) and copied to the FLAC^3D installation directory, so the model file can be invoked from the command line.

4.2. Verification

To verify the rationality of the strain hardening and softening constitutive model for hard brittle rocks under different confining pressures, triaxial compression tests on granite and diabase are simulated using FLAC^3D software. Because the numerical model carries a uniform load, the number of elements in the numerical model has little influence on the simulation results [44], and only a single-grid model (1 m³) is established, as shown in Figure 9. The boundary conditions are as follows:

(i) The displacement at the bottom boundary is restricted when the model runs;

(ii) A stress σ_x = σ_y is applied to the surface of the model to simulate the confining pressure;

(iii) A constant displacement velocity u_z (10⁻⁸ m/step) is applied to the upper surface of the model to simulate axial displacement loading.

Numerical triaxial compression tests are performed on Beishan granite [26] under six confining pressures (σ_x = 0.2 MPa, 1.0 MPa, 3.5 MPa, 10 MPa, 20 MPa and 40 MPa) and on diabase [37] under four confining pressures (σ_x = 5 MPa, 10 MPa, 20 MPa and 40 MPa). The complete stress-strain curves of the two hard brittle rocks under different confining pressures are shown in Figure 10.

Figure 10 illustrates that the simulated results agree well with the test data. The initial yield stress, peak stress, peak strain and residual stress values at different confining pressures, especially at 0.2 MPa, 1.0 MPa, 20 MPa and 40 MPa for granite and 5 MPa, 10 MPa and 40 MPa for diabase are satisfactorily captured by the strain hardening and softening model, which demonstrates that the proposed constitutive model can describe the nonlinear behaviour of rock. The constitutive model can adequately capture the transition from brittle failure to plastic failure under high confining pressures and well describe the confining pressure-dependent characteristics of hard brittle rocks, which verifies the rationality and correctness of the strain hardening and softening constitutive model for hard brittle rocks under different confining pressures. It should be noted that the simulated peak stress for granite at 10 MPa in Figure 10a differs significantly from the experimental peak stress. This occurs due to the deviation in fitting the peak stress of each confining pressure with the Mohr–Coulomb linear yield criterion. In addition, the discrete experimental data may result in deviations when fitting the width parameters ξ_c and ξ_φ with the discrete data of the confining pressures, so the post-peak softening behaviours at a few confining pressure values (e.g., granite at 3.5 MPa, diabase at 20 MPa) are inconsistent with the true experimental curve.

5. Conclusions

According to the triaxial compression test data of hard brittle rocks such as granite and diabase, a strain hardening and softening constitutive model for hard brittle rocks is proposed. The normalised hardening factor and softening factor are defined, which can satisfactorily characterise the same yield state under different confining pressures. The relationship between the hardening and softening factors, cohesion, the internal friction angle and axial plastic strain can be easily obtained by linear interpolation. Based on the Mohr–Coulomb yield criterion, a strength parameter evolution model modified by the sigmoid function is proposed to consider the influence of various confining pressures, which can well describe the real evolution of the strength parameter during the process of rock from yield to failure.

The code of the strain hardening and softening constitutive model is written, and the secondary development of the constitutive model is implemented. Triaxial compression simulation tests are conducted on granite and diabase. The simulation results show that the constitutive model can satisfactorily capture the nonlinear behaviour of rocks, including the strain hardening behaviour before the peak stress and strain softening behaviour after peak stress, the tendency to change from brittle failure to plastic failure, and the effect of different confining pressures on stress-strain curves. The results presented in this study can provide preliminary guidance for determining the strength and strain of hard brittle rocks in the hardening and softening stages. In addition, to simplify the numerical calculation, the dilation angle ψ is assumed to be constant [36]. In future research, a mobilised dilation angle model will be considered to predict the deformation of the tunnel's working face.