Modelling three-dimensional stress-dependent failure of hard rocks

Abstract

Hard rocks exhibit three-dimensional (3D) stress-dependent failure under true triaxial compression. The deformability and strength of hard rocks under true triaxial compression differ from those under traditional loading schemes of conventional triaxial compression tests. For the purpose of characterizing these distinctive features, including 3D stress-dependent brittleness and the failure process and 3D stress-induced anisotropy, a new model suited for hard rocks was proposed in this paper. In the new model, the 3D stress-dependent brittleness under true triaxial compression tests was reflected by a determination method of the internal variable considering the influences of both ${\sigma }_{2}$ and ${\sigma }_{3}$. For the failure process, different evolutions in cohesion and the friction angle were applied to realize different failure processes under different 3D stresses, and their 3D stress dependency was identified by the 3D brittleness index defined in this paper. The 3D stress-induced anisotropy involved in the deformation and failure of hard rocks under true triaxial compression was described via the deformation modulus evolution. The formulation for the stress-induced stiffness matrix and the model framework is fully thermodynamically consistent. The model was implemented in the 3D elasto-plastic cellular automaton system, and good agreement was achieved between the numerical simulation and experimental results, indicating that the new model can be applied to describe the failure behaviours of hard rocks under 3D stress.

Appendix: Thermodynamic aspects of the model framework and stress-induced anisotropy formulation

The constitutive equations have been constructed without details of thermodynamic aspects in this paper. Therefore, it is useful to check the satisfaction of the proposed model and stress-induced anisotropic method, as the internal variable determination and several evolution laws are simply postulated instead of being derived from the dissipation potential.

The thermodynamic formalism considered here is mainly based on the following hypotheses:

• Small transformation hypothesis

  The non-symmetrical part of the displacement gradient is not considered, indicating the rotation is ignored.

• Strain partition hypothesis

  The total strain contains the reversible part (elastic strain) and the irreversible part (plastic strain in this paper), which means:

  $${\varvec{\varepsilon}}={{\varvec{\varepsilon}}}^{e}+{{\varvec{\varepsilon}}}^{p}$$

  (30)

• Isothermal condition

  Dissipation is purely mechanical without any thermal dissipation.

The general Clausius–Duhem inequality is written as follows:

$$\mathcal{L}={\varvec{\sigma}}:\dot{{\varvec{\varepsilon}}}-\rho \left(\dot{\boldsymbol{\Psi }}+s\dot{T}\right)-\frac{1}{T}{\varvec{q}}\Delta T\ge 0$$

(31)

where ${\varvec{\sigma}}$ and ${\varvec{\varepsilon}}$ are stress and strain tensor. $\rho$ represents the density, and $\boldsymbol{\Psi }$ is the free energy. $s$, $T$ and ${\varvec{q}}$ are specific entropy, temperature and heat flux, respectively. Considering the isothermal condition, the dissipation inequality condition can be simplified as follows:

$$\mathcal{L}={\varvec{\sigma}}:\dot{{\varvec{\varepsilon}}}-\rho \dot{\boldsymbol{\Psi }}\ge 0$$

(32)

In our model, the Helmholtz free energy was applied as the potential. Here, the most natural choice for specific free energy is listed as follows:

$$\rho \boldsymbol{\Psi }=\frac{1}{2}\left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right):{\varvec{E}}\left({\varvec{D}}\right):\left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right)$$

(33)

Then, the following equations are acquired:

$${\varvec{\rho}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{\partial{\varvec{\varepsilon}}}=\boldsymbol{ }{\varvec{E}}\left({\varvec{D}}\right): \left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right)={\varvec{\sigma}}$$

(34)

$${\varvec{\rho}}\frac{\partial \boldsymbol{\Psi }}{\partial {{\varvec{\varepsilon}}}^{p}}=\boldsymbol{ }-{\varvec{E}}\left({\varvec{D}}\right): \left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right)=-{\varvec{\sigma}}$$

(35)

The inequality can be calculated as:

$${\varvec{\mathcal{L}}}=\left({\varvec{\sigma}}-{\varvec{\rho}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{{\varvec{\partial}}{\varvec{\varepsilon}}}\right) :{\dot{{\varvec{\varepsilon}}}}-{\varvec{\rho}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{\partial {{\varvec{\varepsilon}}}^{p}} :{\dot{{{\varvec{\varepsilon}}}^{p}}-{\varvec{\rho}}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{\partial {\varvec{D}}}:{\dot{{\varvec{D}}}}={\varvec{\sigma}} :{\dot{{{\varvec{\varepsilon}}}^{p}}-{\varvec{\rho}}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{{\varvec{\partial}} {\varvec{D}}} :{\dot{{\varvec{D}}}}={\varvec{\sigma}}:{\dot{{{\varvec{\varepsilon}}}^{p}}}+{\varvec{Y}} :{\dot{{\varvec{D}}}}\ge 0$$

(36)

Equation 36 means that the dissipative thermodynamic force conjugate to the plastic strain is the nominal stress ${\varvec{\sigma}}$. The damage conjugate force is defined as follows:

$${\varvec{Y}}=-{\varvec{\rho}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{{\varvec{\partial}} {\varvec{D}}}$$

(37)

1.1 Stress-induced anisotropic damage dissipation

The damage conjugate force can be calculated as follows:

$${\varvec{Y}}=-{\varvec{\rho}}\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{\partial {\varvec{E}}\left({\varvec{D}}\right)}:\frac{{\varvec{\partial}} {\varvec{E}}\left({\varvec{D}}\right)}{\partial {\varvec{D}}}$$

(38)

According to Eq. 33, $\frac{{\varvec{\partial}} \boldsymbol{\Psi }}{\partial {\varvec{E}}\left({\varvec{D}}\right)}$ will never be negative as the total strain ${\varvec{\varepsilon}}$ cannot be smaller than the plastic strain ${{\varvec{\varepsilon}}}^{p}$. It is equal to the partial derivative of the specific free energy ($\frac{1}{2}\left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right):{{\varvec{E}}}_{0}:\left({\varvec{\varepsilon}}-{{\varvec{\varepsilon}}}^{p}\right)$) for the undamaged material under the same elastic strain. Now, only $\frac{{\varvec{\partial}} {\varvec{E}}\left({\varvec{D}}\right)}{{\varvec{\partial}} {\varvec{D}}}$ is considered.

Considering the stress-induced anisotropic matrix in the damage principal axes space, it can be expressed as [14]:

$${E}_{ijkl}={\alpha }_{ijpq}{E}_{pqrs}^{0}{\alpha }_{klrs}$$

(39)

where ${E}_{pqrs}^{0}$ is the component of the initial stiffness matrix. For the definition of the integrity in this paper, the $\boldsymbol{\alpha }$ matrix can be expressed as [14]:

$${\boldsymbol{\alpha }}=\left[\begin{array}{cccccc}{\beta }_{(1)}& 0& 0& 0& 0& 0\\ 0& {\beta }_{(2)}& 0& 0& 0& 0\\ 0& 0& {\beta }_{(3)}& 0& 0& 0\\ 0& 0& 0& \sqrt{{\beta }_{(1)}{\beta }_{(2)}}& 0& 0\\ 0& 0& 0& 0& \sqrt{{\beta }_{(2)}{\beta }_{(3)}}& 0\\ 0& 0& 0& 0& 0& \sqrt{{\beta }_{(3)}{\beta }_{(1)}}\end{array}\right]$$

(40)

Consider the stiffness matrix as a function shown below:

$${\varvec{E}}\left({\varvec{D}}\right)=f({{\varvec{E}}}_{0},{\varvec{\beta}})$$

(41)

Hence, $\frac{\partial {\varvec{E}}\left({\varvec{D}}\right)}{\partial {\varvec{D}}}$ can be expressed as:

$$\frac{\partial {\varvec{E}}\left({\varvec{D}}\right)}{\partial {\varvec{D}}}=\frac{\partial {\varvec{E}}\left({\varvec{D}}\right)}{\partial {{\varvec{E}}}_{0}}\frac{\partial {{\varvec{E}}}_{0}}{\partial {\varvec{D}}}+\frac{\partial {\varvec{E}}\left({\varvec{D}}\right)}{\partial \boldsymbol{\alpha }}\frac{\partial \boldsymbol{\alpha }}{\partial{\varvec{\beta}}}\frac{\partial{\varvec{\beta}}}{\partial {\varvec{D}}}=2\boldsymbol{\alpha }\frac{\partial \boldsymbol{\alpha }}{\partial{\varvec{\beta}}}\frac{\partial{\varvec{\beta}}}{\partial {\varvec{D}}}$$

(42)

In the above equation, $\frac{\partial {{\varvec{E}}}_{0}}{\partial {\varvec{D}}}$ is zero as ${{\varvec{E}}}_{0}$ has no relation to ${\varvec{D}}$. As shown in Eq. 40, $\boldsymbol{\alpha }$ and $\frac{\partial \boldsymbol{\alpha }}{\partial{\varvec{\beta}}}$ are non-negative values since ${\beta }_{(j)}$ varies from 0 to 1. $\frac{\partial{\varvec{\beta}}}{\partial {\varvec{D}}}$ will be negative as the integrity and damage vary oppositely. In addition, $\dot{{\varvec{D}}}$ in Eq. 36 is the damage evolution ratio, and it is positive according to the definition in Eq. 13 in this paper. Therefore, the damage dissipation ${\varvec{Y}}:\dot{{\varvec{D}}}$ is non-negative, which means:

$${\mathcal{L}}_{{\varvec{D}}}={\varvec{Y}}:\dot{{\varvec{D}}}\ge 0$$

(43)

1.2 Plastic dissipation

For the plastic dissipation, the inequality can be calculated as follows:

$${\mathcal{L}}_{{\varvec{P}}}={\varvec{\sigma}}:\dot{{{\varvec{\varepsilon}}}^{p}}={\varvec{\sigma}}:\frac{\partial {\varvec{g}}}{\partial {{\varvec{\sigma}}}^{\mathrm{eff}}}$$

(44)

where $g$ is the plastic potential as defined in Sect. 3.1. According to the framework in this paper, the nominal stress can be expressed as follows:

$${\varvec{\sigma}}=\boldsymbol{\alpha }:{{\varvec{\sigma}}}^{\mathrm{eff}}$$

(45)

$\boldsymbol{\alpha }$ is the same as in Eqs. 10 and 11. Substituting Eq. 45 into Eq. 44, the following equation is acquired:

$${\mathcal{L}}_{{\varvec{P}}}={\varvec{\sigma}}:\dot{{{\varvec{\varepsilon}}}^{p}}=\boldsymbol{\alpha }:{{\varvec{\sigma}}}^{\mathrm{eff}}:\frac{\partial {\varvec{g}}}{\partial {{\varvec{\sigma}}}^{\mathrm{eff}}}$$

(46)

The plastic potential surface defined in this paper is convex; therefore, ${{\varvec{\sigma}}}^{\mathrm{eff}}:\frac{\partial {\varvec{g}}}{\partial {{\varvec{\sigma}}}^{\mathrm{eff}}}$ will be non-negative. Since the principal values of $\boldsymbol{\alpha }$ are between 0 and 1, it is easy to find that $\boldsymbol{\alpha }:{{\varvec{\sigma}}}^{\mathrm{eff}}:\frac{\partial {\varvec{g}}}{\partial {{\varvec{\sigma}}}^{\mathrm{eff}}}$ will be non-negative.

Then, the total inequality including stress-induced anisotropic damage and the plastic part can be calculated following Eq. 36:

$$\mathcal{L}={\mathcal{L}}_{{\varvec{D}}}+{\mathcal{L}}_{{\varvec{P}}}\ge 0$$

(47)

which indicates that the proposed model framework and stress-induced anisotropy formulation are fully thermodynamically consistent, and they are rigorous.