A Micromechanical Anisotropic Damage Model for Brittle Rocks With Non-Associated Plastic Flow Rule Under True Triaxial Compressive Stresses

A micromechanical anisotropic damage model with a non-associated plastic flow rule is developed for describing the true triaxial behaviors of brittle rocks. We combine the Eshelby’s solution to the inclusion problem with the framework of irreversible thermodynamics. The main dissipative mechanisms of inelastic deformation due to the frictional sliding and damage by microcrack propagation are strongly coupled to each other. A Coulomb-type friction criterion is formulated in terms of the local stress applied onto the microcracks as the yielding function. The back-stress term contained in this local stress plays a critical role in describing the material’s hardening/softening behaviors. With a non-associated flow rule, a potential function is involved. Some analytical analysis of the non-associated micromechanical anisotropic damage model are conducted, which are useful for the model parameters calibration. The proposed model is used to simulate the laboratory tests on Westerly granite under true triaxial stresses. Comparing the numerical simulation results provided by the models with associated/non-associated plastic flow rule and experimental results, it is clear that the proposed non-associated model gives a better prediction than the previous associated model.

1 Introduction

Constitutive model and simulation of the mechanical behaviors of heterogeneous rocks under general stress conditions have a great significance on the investigation of the safety and stability of rock engineering. Since host rocks in underground engineering are most often in a true triaxial stress state, researchers have conducted many true triaxial compression tests [1–7]. In 1967, [8] firstly developed an experimental apparatus that combined torsion and triaxial compression. It provided general triaxial stress states but did not produce the conditions of homogeneous triaxial stress. Then [9, 10] implemented homogeneous triaxial stress using triaxial cells. [11] used Mogi’s cell design to fabricate a cell to test larger sample sizes. After that, [4,12] did a series of experiments on different rocks using the true triaxial loading system at the University of Wisconsin. Recently, Feng et al. [13,14] designed a novel Mogi type true triaxial testing apparatus and utilized it to obtain complete stress-strain curves of many hard Rocks.

In terms of constitutive model, a number of researchers [15–18] have developed numerical models or used commercial software to study the failure processes of rocks under polyaxial stress conditions. [19] established a damage softening statistical constitutive model with the assumption that the rock micro-unit failure obeys the Weibull random distribution under the true triaxial stress state. Through linear fitting method and elastoplastic mechanical analysis, [20] built a true triaxial constitutive model of coal rock under horizontal stress loading. At present, it is very common to use the isotropic assumption in most constitutive models. However, geomaterials, such as rocks and concrete, display stress-induced anisotropy in their mechanical properties. Considering it difficult in modeling the anisotropy, few researchers have established anisotropic theoretical models with physical meaning for simulating the true triaxial stress-strain curve of rocks.

This paper will extend the previous work [21] for modeling the true triaxial mechanical behavior of hard rocks. The micromechanical anisotropic damage model [21] has been established with an associated plastic flow rule. However, a large number of experimental and theoretical results show that a non-associated plastic flow rule must be adopted to more accurately describe the inelastic deformation of rocks. To this end, a micromechanical anisotropic damage model with a non-associated plastic flow rule will be developed, which will be used to simulate laboratory tests on Westerly granite under true triaxial stresses.

The following tensorial product notations are used throughout this paper: ${\left(\mathit{A}\otimes \mathit{B}\right)}_{ijkl}=A_{ij}{B}_{kl}$, ${\left(\mathit{A}\bar{\underline{\otimes }}\mathit{B}\right)}_{ijkl}=\left(A_{ik}{B}_{jl}+A_{il}{B}_{jk}\right)/2$, ${\left(\mathit{a}\otimes \mathit{b}\right)}_{ij}=a_i{b}_j$, ${\left(\mathit{a}\stackrel{s}{\otimes }\mathit{b}\right)}_{ij}=\left(a_i{b}_j+a_j{b}_i\right)/2$.

2 Formulation of Non-assocaited Micromechanical Anisotropic Damamge Model

In this section, a micromechanical anisotropic damage model with a non-associated plastic flow rule for brittle rocks is formulated. A representative elementary volume (REV), defined by a geometrical domain Ω and its boundary surface ∂Ω, is shown in Figure 1. The relevant REV as a matrix-inclusion system for microcracked rocks is composed of an isotropic linearly elastic matrix with a stiffness tensor ${\mathbb{C}}^s$ (or compliance tensor ${\mathbb{S}}^s$) and a large number of randomly distributed microcracks with the elasticity tensor ${\mathbb{C}}^c$. Based on the Eshelby’s solution, microcracks with a unit normal vector of $\bar{n}$ are treated as penny-shaped and microcracks with the same normal vector $\bar{n}$ are placed into the same family. The volume fraction of the rth family of microcracks is expressed as

${\phi }^{c,r}=\frac{4}{3}\pi {(a^r)}^2{c}^r\mathcal{N=}\frac{4}{3}\pi {\zeta }^r{d}^r(1)$

where $\mathcal{N}$ is the number of microcracks per unit volume of all microcracks in the considered family and a^r and c^r denote the average radius and the half opening of the rth family of microcracks, respectively (see Figure 2). ζ = c/a ≪ 1 is defined as the aspect ratio of the microcracks; $d^r=\mathcal{N}{(a^r)}^3$ is the microcrack density parameter initially introduced by [22] and widely used as an internal damage variable.

FIGURE 1

FIGURE 1. Representative elementary volume (REV) of brittle rocks.

FIGURE 2

FIGURE 2. A penny-shaped microcrack and its orientation angles.

2.1 Strain Decomposition

The domain occupied by a family of microcracks (with the same unit normal vector $\bar{n}$) is denoted $\mathcal{S}$ and the displacements of the upper surface ${\mathcal{S}}^{+}$ and lower surface ${\mathcal{S}}^{-}$ of the microcracks are designated u⁺ and u⁻, respectively. Then, the related displacement jump between the two surfaces is expressed by $\left[\bar{u}\right]=u^{+}-u^{-}$. The unilateral contact conditions on $\mathcal{S}$ are taken into consideration, i.e.,

$\left[u_n\right]\ge 0,{\sigma }_n^c\le 0,\left[u_n\right]{\sigma }_n^c=0(2)$

where ${\sigma }_n^c$ and u_n are the normal component of the local stress tensor and the displacement jump $\left[\bar{u}\right]$, respectively.

The displacement discontinuity is characterized by two variables:

i) a scalar, β, representing the microcrack opening degree

$\beta =\mathcal{N}{\int }_{{\mathcal{S}}^{+}}\left[u_n\right]\text{d}\mathcal{S}(3)$

ii) a sliding vector, $\bar{\gamma }$, quantifying the relative sliding along the microcrack plane

$\bar{\gamma }=\mathcal{N}{\int }_{{\mathcal{S}}^{+}}\left[{\bar{u}}_t\right]\text{d}\mathcal{S}(4)$

with $\left[{\bar{u}}_t\right]=\left[\bar{u}\right]-\left[u_n\right]\bar{n}$.

The inelastic deformation due to the displacement discontinuities [23] generated by a family of penny-shaped microcracks of normal $\bar{n}$, denoted ϵ^p, takes the following general form

${\mathit{ϵ}}^p=\mathcal{N}{\int }_{S^{+}}\bar{n}\stackrel{s}{\otimes }\left[\bar{u}\right]\text{d}S=\beta \left(\bar{n}\otimes \bar{n}\right)+\bar{\gamma }\stackrel{s}{\otimes }\bar{n}(5)$

All the microcracks in the REV are discontinuities dispersed in the solid matrix. Accordingly, the total strain of REV is decomposed into two terms: an elastic part ɛ^e, which is the result of the deformation of the matrix phase and an inelastic part ɛ^p, which is due to the existence of microcracks.

$\mathit{\epsilon }={\mathit{\epsilon }}^e+{\mathit{\epsilon }}^p(6)$

where ɛ^p can be expressed by the simple sum of each family’s microcrack contribution

${\mathit{\epsilon }}^p=\sum _{r=1}^{n_c}{\varrho }^r{\mathit{ϵ}}^{p,r}(7)$

and ϱ^r is the weight of the rth family of microcracks.

2.2 Effective Elastic Properties and Free Energy

The Mori-Tanaka scheme [24] is taken into account within the standard linear homogenization framework, and the effective stiffness tensor ${\mathbb{C}}^{\mathrm{hom}}$ has the following general form [21]:

${\mathbb{C}}^{\mathrm{hom}}={\mathbb{C}}^s:{\left(\mathbb{I+}\sum _{r=1}^{n_c}{\phi }^c{\left(\mathbb{I}-{\mathbb{S}}_e^r\right)}^{-1}\right)}^{-1}(8)$

where ${\mathbb{C}}^s$ is the stiffness tensor of the matrix phase; ${\mathbb{S}}_e^r$ is the classical Eshelby tensor, such that ${\mathbb{S}}_e^r={\mathbb{P}}_e^r:{\mathbb{C}}^s$ with ${\mathbb{P}}_e^r$ being the so-called fourth-order Hill tensor. It has been proved that when microcracks are open and the aspect ratio ζ tends to zero, ${\mathbb{S}}_e^r$ has an analytical expression. In this case, the effective compliance tensor of the REV is attained by the inverse of Eq. 8 and is written in the following form:

${\mathbb{S}}^{\mathrm{hom}}={\mathbb{S}}^s+\sum _{r=1}^{n_c}{d}^r\left(\frac{1}{k_n}{\mathbb{E}}^{2,r}+\frac{1}{k_t}{\mathbb{E}}^{4,r}\right)(9)$

with ${\mathbb{S}}^s={\left({\mathbb{C}}^s\right)}^{-1}$ being the compliance tensor. k_n and k_t are two constants expressed as $k_n=\frac{3E^s}{16\left(1-{\left({\nu }^s\right)}^2\right)}$ and $k_t=k_n\left(2-{\nu }^s\right)$ [25]. E^s and ν^s are the Young’s modulus and the Poisson’s ratio of the matrix, respectively. ${\mathbb{E}}^2$ and ${\mathbb{E}}^4$ are the fourth-order orientation-dependent normal and tensor operators, which are both functions of the unit normal vector $\bar{n}$. The specific applications of ${\mathbb{E}}^2$ and ${\mathbb{E}}^4$ are described in detail in literature [26].

${\mathbb{E}}^2=\bar{n}\otimes \bar{n}\otimes \bar{n}\otimes \bar{n},{\mathbb{E}}^4=\left(\bar{n}\otimes \bar{n}\right)\bar{\underline{\otimes }}\left(\mathrm{\delta }-\bar{n}\otimes \bar{n}\right)+\left(\mathrm{\delta }-\bar{n}\otimes \bar{n}\right)\bar{\underline{\otimes }}\left(\bar{n}\otimes \bar{n}\right)(10)$

Since microcrack propagation gives rise to the degradation of the material stiffness, it is possible to define an alternative macroscopic damage variable ω^r according to the relative variation of the Young’s modulus [27,28].

${\omega }^r\left(\bar{n}\right)=1-\frac{E^r\left(\bar{n}\right)}{E^s}(11)$

where $E\left(\bar{n}\right)$ is the longtitudinal Young’s modulus of the matrix in the direction $\bar{n}$. In this way, we can link the macro damage ω to the micro damage d.

Theoretical studies [29,30] have shown that $E\left(\bar{n}\right)$ can be formulated by means of the compliance tensor ${\mathbb{S}}^{\mathrm{hom}}$ with the following expression:

$E^r\left(\bar{n}\right)={\left[\left({\bar{n}}^r\otimes {\bar{n}}^r\right):{\mathbb{S}}^{\mathrm{hom}}:\left({\bar{n}}^r\otimes {\bar{n}}^r\right)\right]}^{-1}(12)$

Along with the previous work [21], for any opening/closure combination of the microcrack families, the free energy takes the general form

$W=\frac{1}{2}\mathit{\sigma }:{\mathbb{S}}^s:\mathit{\sigma }+\mathit{\sigma }:\sum _{r=1}^{n_c}{\mathit{ϵ}}^{p,r}-\sum _{r=1}^{n_c}\frac{1}{2}{\mathit{ϵ}}^{p,r}:{\mathbb{C}}^{n,r}:{\mathit{ϵ}}^{p,r}(13)$

with ${\mathbb{C}}^{n,r}=\frac{k_n}{d^r}{\mathbb{E}}^{2,r}+\frac{k_t}{d^r}{\mathbb{E}}^{4,r}$.

2.3 Damage Criterion

The damage criterion is a function of the thermodynamic force F^d,r associated with the internal damage variable d^r, which can be derived by applying the standard differentiation of the macroscopic free energy W.

$F^{d,r}=\frac{\partial W}{\partial d^r}=\frac{1}{2{(d^r)}^2}{\mathit{ϵ}}^{p,r}:{\mathbb{C}}^{n,r}:{\mathit{ϵ}}^{p,r}(14)$

In regards to the damage evolution law, a strain energy release rate-based damage criterion is largely adopted for all microcrack families:

$g^r\left(F^{d,r},d^r\right)=F^{d,r}-\mathcal{R}\left(d^r\right)\le 0(15)$

where $\mathcal{R}\left(d^r\right)$ represents the current resistance to further damage propagation for the rth family of microcracks. It is usually assumed that the damage resistance function should be dependent on the damage level. In order to describe the strain hardening/softening of materials induced by the microcracks’ coalescence, the following power form for $\mathcal{R}\left(d^r\right)$ [31] is adopted:

$\mathcal{R}\left(d^r\right)=\mathcal{R}\left(d_c\right)G\left({\kappa }^r\right),G\left({\kappa }^r\right)=\frac{2{\kappa }^r}{{({\kappa }^r)}^2+1}(16)$

where κ^r = d^r/d_c is defined as a dimensionless parameter. Physically, d_c represents the critical damage values corresponding to the peak stress, and $\mathcal{R}\left(d_c\right)$ is the maximum resistance value to the damage propagation.

When the damage criterion Eq. 15 is satisfied, the damage evolution rate is determined by using the normality rule:

$\dot{d^r}={\lambda }^{d,r}\frac{\partial g^r}{\partial F^{d,r}}={\lambda }^{d,r}(17)$

with λ^d,r being a non-negative damage multiplier for the rth family of microcracks.

2.4 Friction Criterion With Non-associated Flow Rule

When microcracks are open, there is no friction effect between the microcracks. Therefore, we only consider the friction criterion for closed microcracks. The generalized Coulomb criterion is used as the yielding function to describe the friction sliding occurring along the closed microcracks. Given the macroscopic free energy in Eq. 13, the thermodynamic force associated with the local inelastic strain ϵ^p, denoted by the local stress σ^c, is deduced by the standard derivation of W with respect to ϵ^p:

${\mathit{\sigma }}^{c,r}=\frac{\partial W}{\partial {\mathit{ϵ}}^{p,r}}=\mathit{\sigma }-{\mathbb{C}}^{n,r}:{\mathit{ϵ}}^{p,r}(18)$

where ${\mathbb{C}}^{n,r}:{\mathit{ϵ}}^{p,r}$ plays a critical role in describing the materials’ hardening/softening behaviors.

At the microscopic scale, the Coulomb criterion is determined by the normal and tangential components represented by ${\sigma }_n^c$ and τ^c of the local stress σ^c:

$f\left({\mathit{\sigma }}^{c,r}\right)=‖{\mathit{\tau }}^{c,r}‖+c_f{\sigma }_n^{c,r}\le 0(19)$

where c_f is the friction coefficient of the microcracks, which is related to the asperity of microcracks’ surfaces [32].

For most brittle rocks, a large number of true triaxial tests indicated that the associated plastic flow rule was generally not suitable to describe the volumetric deformation. Therefore, we propose the following plastic potential function. To be consistent, this function is similar to the friction criterion:

$G\left({\mathit{\sigma }}^{c,r}\right)=‖{\mathit{\tau }}^{c,r}‖+c_v{\sigma }_n^{c,r}=0(20)$

where c_v is termed as the current volumetric dilation coefficient. When c_v = c_f, the friction function is completely consistent with the potential function. At this time, the model adopts the associated plastic flow rule.

In the classical plasticity theory, the evolution rate of the local inelastic strain ε^p is given by the normality rule:

${\dot{\mathit{ϵ}}}^{p,r}={\lambda }^{p,r}\frac{\partial G}{\partial {\mathit{\sigma }}^{c,r}}={\lambda }^{p,r}{\mathit{D}}^{n,r}(21)$

where λ^p,r is a non-negative friction multiplier for the rth family of microcracks and D^n,r is served as the plastic flow direction by a second-order tensor,

${\mathit{D}}^{n,r}=\left(\frac{{\mathbb{E}}^{4,r}}{2‖{\mathit{\tau }}^{c,r}‖}+c_v\frac{{\mathbb{E}}^{2,r}}{{\sigma }_n^{c,r}}\right):{\mathit{\sigma }}^{c,r}(22)$

3 Model Parameters Calibration Method

The proposed model only contains six material constants or model parameters, each having a clear physical meaning. Using a series of conventional triaxial compression tests under different confining pressures, the model parameters can be determined. Before discussing the model parameters calibration, some analytical analyses of the non-associated micromechanical anisotropic damage model under conventional triaxial compression are first conducted.

3.1 Analytical Expression of Peak Stress and Crack Damage Stress

Under the loading path of conventional triaxial compression, the local friction criterion (19) for the critical plane with θ = θ_c which satisfies the condition $\tan {\theta }_c=c_f+\sqrt{c_f^2+1}$ has the following form [21]:

$f={\sigma }_1-{\tan }^2{\theta }_c{\sigma }_3-2\tan {\theta }_c\left(\frac{k_t\mathit{\gamma }.\mathit{t}}{2d}+\frac{c_f{k}_n\beta }{d}\right)=0(23)$

with $\mathit{t}=\frac{{\mathit{\tau }}^{c,r}}{‖{\mathit{\tau }}^{c,r}‖}$ being the shear flow direction within the microcrack plane.

It is possible to define the plastic multiplier in such a way that Λ^p = ∫λ^p, then $\bar{\gamma }$ and β can be calculated as follows:

$\bar{\gamma }={\mathrm{\Lambda }}^p\mathit{t},\beta ={\mathrm{\Lambda }}^p{c}_v(24)$

By substituting Eq. 24 into Eq. 23, we can derive the following expression by defining ${\xi }_1=\frac{k_t}{2}+k_n{c}_f{c}_v$:

$f={\sigma }_1-{\tan }^2{\theta }_c{\sigma }_3-2\tan {\theta }_c\frac{{\mathrm{\Lambda }}^p}{d}{\xi }_1=0(25)$

By combining Eqs. 14, 21, 22, the damage criterion (15) is reformulated as

$g=\frac{1}{2}{\left(\frac{{\mathrm{\Lambda }}^p}{d}\right)}^2\left(\frac{k_t}{2}+c_v^2{k}_n\right)-\mathcal{R}\left(d\right)=0(26)$

with ${\xi }_2=\frac{k_t}{2}+k_n{c}_v^2$. The following relation can be derived from the damage criterion (26):

$\frac{{\mathrm{\Lambda }}^p}{d}=\sqrt{\frac{2\mathcal{R}\left(d\right)}{{\xi }_2}}(27)$

For the critical sliding plane, the friction criterion can be expressed in terms of ξ₁ and ξ₂:

$f={\sigma }_1-{\tan }^2{\theta }_c{\sigma }_3-2\tan {\theta }_c\sqrt{\frac{2\mathcal{R}\left(d_c\right)}{{\xi }_2}}{\xi }_1=0,\text{with}d=d_c(28)$

Finally, the analytical expression of peak stress can be derived as:

${\sigma }_1={\tan }^2{\theta }_c{\sigma }_3+2\tan {\theta }_c\sqrt{\frac{2\mathcal{R}\left(d_c\right)}{{\xi }_2}}{\xi }_1(29)$

On the other hand, in the volume strain-deviatoric stress curve, the volume strain will reverse with the increase of deviatoric stress, and this point is called the volume compressibility dilatancy (C/D) transition point. A series of researches [33,34] have shown that the crack damage stress σ_cd is defined as the volume compressibility dilatancy (C/D) transition stress. Inspired by the strength prediction, we assume that there also exists a critical damage value d_cd corresponding to σ_cd. On the basis of Eq. 29, the following crack damage stress can be expressed as:

${\sigma }_{cd}={\tan }^2{\theta }_c{\sigma }_3+2\tan {\theta }_c\sqrt{\frac{2\mathcal{R}\left(d_c\right)\varpi }{{\xi }_2}}{\xi }_1\text{with}\varpi =\frac{2{\kappa }_{cd}}{{({\kappa }_{cd})}^2+1},{\kappa }_{cd}=\frac{d_{cd}}{d_c}(30)$

According to the classical elastic-plastic theory, the volume strain ɛ_v is the sum of elastic volume strain ${\epsilon }_v^e$ and plastic volume strain ${\epsilon }_v^p$:

${\epsilon }_v={\epsilon }_v^e+{\epsilon }_v^p(31)$

By ignoring the volumetric strain generated by the confining pressure and the inelastic strain generated by the initial crack closure and using Eqs. 21, 31, the volumetric strain generated in the axial loading phase can be rewritten into the following form:

${\epsilon }_v=\mathit{\delta }:{\mathbb{S}}^s:\left(\mathit{\sigma }-{\sigma }_3\mathit{\delta }\right)-{\mathrm{\Lambda }}^p{c}_v(32)$

Inserting Eq. 27 into the above formula, the volumetric strain ϵ_v becomes:

${\epsilon }_v=\frac{1}{3k^s}\left(\left({\tan }^2{\theta }_c-1\right){\sigma }_3+2\tan {\theta }_c\sqrt{\frac{2\mathcal{R}\left(d\right)}{{\xi }_2}}{\xi }_1\right)-c_{v}d\sqrt{\frac{2\mathcal{R}\left(d\right)}{{\xi }_2}}(33)$

The increment of volume strain is equal to 0 at the volume C/D transition point, i.e. $\text{d}{ϵ}_v=\frac{\partial {ϵ}_v}{\partial \kappa }\frac{\partial \kappa }{\partial d}\text{d}d=0$. At C/D transition point, d = d_cd and κ = κ_cd, we can get the following characteristic equation:

$c_v{d}_c{\kappa }_{cd}^3+\vartheta {\xi }_1{\kappa }_{cd}^2+3c_v{d}_c{\kappa }_{cd}-\vartheta {\xi }_1=0(34)$

with $\vartheta =\frac{2\tan {\theta }_c}{3k^s}$.

3.2 Model Parameters Calibration

With the above analytical analysis at hand, we here discuss the model parameters calibration procedure as follows:

• The means of Young’s modulus E^s and Poisson’s ratio ν^s of can be determined using the linear part of the stress-strain curves.

• The friction coefficient c_f can be obtained by comparing strength criterion Eq. 29 with the peak stress envelope.

• Parameter d_c has no influence on peak strength of materials and is related to the deformation at peak strength state and the post-peak stress-strain curve. With the increase of d_c, the deformation at peak strength is increasing [21].

• The parameter c_v can be identified by solving Eq.30, 34. In this process, the intermediate variable ϖ must first be determined by comparing the peak stress envelope with the crack damage stress envelope.

• The critical damage resistance $\mathcal{R}\left(d_c\right)$ can be calibrated with the analytical expression of peak stress Eq. 29 and c_v.

If the associated flow rule is adopted to describe the inelastic strain of the material, there are only 5 parameters in the model. The parameters calibration method above is still available.

4 Model Application on Westerly Granite

In this section, the proposed model is applied to modeling the mechanical behavior of Westerly granite subjected to true triaxial compression loading. Firstly conventional triaxial compression test results on Westerly granite conducted by [4] are adopted to determine the model parameters. With the parameters calibration method in Section 4, we obtained the associated/non-associated micromechanical anisotropic damage models parameters and listed as follows:

• Associated micromechanical anisotropic damage model parameters: E^s = 68000MPa, ν^s = 0.21, d_c = 18, c_f = 1.27, $\mathcal{R}(d_c)=2.59\times 10^{-2}$.

• Non-associated micromechanical anisotropic damage model parameters: E^s = 68000MPa, ν^s = 0.21, d_c = 18, c_f = 1.27, c_v = 0.45, $\mathcal{R}(d_c)=3.32\times 10^{-2}$.

4.1 Numberical Simulations

With the above parameters and using the plastic-damage decoupled correction (PDDC) numerical algorithm [35], numerical simulations of assocaited/non-associated micromechanical anisotropic damage models on true triaxial comoression tests of Westerly granite are conducted.

Figure 3 shows the curves of the deviator stress $\left({\sigma }_1-{\sigma }_3\right)$ versus three principal strains from the series of tests in which the σ₃ magnitude is kept at 60 MPa and σ₂ is varied from 60 to 249 MPa. The numerical simulation results provided by the associated and non-associated micromechanical anisotropic damage model are compared. The dotted lines in Figure 3 are the simulated results using the associated flow rule, and the solid lines are the results with the non-associated flow rule. The comparison of the two simulated results shows the non-associated model gives a better prediction than the associated one. Besides, one can see generally a good agreement between the non-associated model’s predictions and experimental data. The proposed non-associated model presents the ability to capture the main characteristics of mechanical behavior of the rock, such as the peak stress, pressure sensitivity, transition from volumetric compaction to dilatancy.

FIGURE 3

FIGURE 3. Comparisons of the models’ predictions with the experimental data for the true triaxial compression tests on Westerly granite (The dotted lines are the associated model’s results and the solid lines are the non-associated model’s results). (A) intermediate principal stress σ₂ = 60 MPa, (B) intermediate principal stress σ₂ = 113 MPa, (C) intermediate principal stress σ₂ = 180 MPa, and (D) intermediate principal stress σ₂ = 249 MPa.

Figure 4 shows the damage density distribution with the different intermediate principal stress values. These 2D coupled plans ${\bar{e}}_1{\bar{e}}_2,{\bar{e}}_1{\bar{e}}_3$ and ${\bar{e}}_2{\bar{e}}_3$ are subsequently approximated by the distribution function Eq. 11, which is defined on a unit sphere and centered on a material point. If o is denoted as the original point in the considered space and p is denoted as a point on the surfaces of the distribution functions, the orientation of the vector $\overrightarrow{op}$ amounts to the family of microcracks with the unit normal $\overrightarrow{op}/‖\overrightarrow{op}‖$. Then, the damage magnitude is evaluated by $‖\overrightarrow{op}‖$. Theoretically, XY and XZ plans will coincide under conventional triaxial loading path (σ₂ = σ₃ = 60 MPa). This result is confirmed in Figure 4A. We can see that the growth of the damage in the plan XY is progressively blocked with the increase in σ₂ from 60 to 249 MPa.

FIGURE 4

FIGURE 4. Damage distribution at the peak stresses in different intermediate principal stress compression test (A) intermediate principal stress σ₂ = 60 MPa, (B) intermediate principal stress σ₂ = 113 MPa, (C) intermediate principal stress σ₂ = 180 MPa, and (D) intermediate principal stress σ₂ = 249 MPa.

5 Conclusion

A new micromechanical anisotropic damage model with a non-associated plastic flow rule has been developed for describing the true triaxial compression behaviors for brittle rocks. Unlike the previous models, the potential function plays a critical role in the integration of the model. Importantly, the non-associated flow rule works quite well to describe the inelastic deformations in both the axial and lateral directions. We have discussed the model parameters calibration procedure based on some analytical analysis of the proposed model under conventional triaxial compression loading. The model has been finally applied to simulate true triaxial stress-strain curves of Westerly granite. One can see that a good agreement between the non-associated model’s predictions and the experimental data for under true triaxial compression loading path. Moreover, the non-associated model could provide a better description of true triaxial mechanical behaviors than that of our previous associated micromechanical anisotropic damage model [21].

It is worth noting that our proposed model only contains six parameters, each having a clear physical meaning. Each of the parameters can be easily identified from a series of conventional triaxial tests. We will present important extensions related to the time-dependent behaviors of brittle rocks under true triaxial compression in future work.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

SY: Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Writing—original draft. QZ: Conceptualization, Methodology. WZ: Software, Visualization, Writing—review and editing, Resources. JZ: Supervision, Writing—review and editing. LZ: Supervision, Writing—review and editing.

Funding

This work has been jointly supported by the Natural Science Foundation of Nanjing Vocational University of Industry Technology (YK20-02-12).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Mogi K. Effect of the Intermediate Principal Stress on Rock Failure. J Geophys Res (1967) 72:5117–31. doi:10.1029/jz072i020p05117

CrossRef Full Text | Google Scholar

2. Wawersik WR, Carlson LW, Holcomb DJ, Williams RJ. New Method for True-Triaxial Rock Testing. Int J Rock Mech Mining Sci (1997) 34:330. doi:10.1016/s1365-1609(97)00049-x

CrossRef Full Text | Google Scholar

3. Chang C, Haimson B. True Triaxial Strength and Deformability of the German continental Deep Drilling Program (Ktb) Deep Hole Amphibolite. J Geophys Res (2000) 105:18999–9013. doi:10.1029/2000jb900184

CrossRef Full Text | Google Scholar

4. Haimson B, Chang C. A New True Triaxial Cell for Testing Mechanical Properties of Rock, and its Use to Determine Rock Strength and Deformability of westerly Granite. Int J Rock Mech Mining Sci (2000) 37:285–96. doi:10.1016/s1365-1609(99)00106-9

CrossRef Full Text | Google Scholar

5. Liu HD, Cao J. Test and Study on Influence of Intermediate Main Stress to Characteristic of Rock Mechanics (2008) 30:56–60+80.

Google Scholar

6. Xiang TB, Feng XT, Chen BR, Jiang Q, Zhang CQ. Rock Failure Mechanism and True Triaxial Experimental Study of Specimens with Single Structure Plane under Three-Dimensional Stress. Rock Soil Mech (2009) 30:2908–16. doi:10.16285/j.rsm.2009.10.053

CrossRef Full Text | Google Scholar

7. Feng X-T, Zhang X, Kong R, Wang G. A Novel Mogi Type True Triaxial Testing Apparatus and its Use to Obtain Complete Stress-Strain Curves of Hard Rocks. Rock Mech Rock Eng (2015) 49:1649–62. doi:10.1007/s00603-015-0875-y

CrossRef Full Text | Google Scholar

8. Handin J, Heard HC, Magouirk JN. Effects of the Intermediate Principal Stress on the Failure of limestone, Dolomite, and Glass at Different Temperatures and Strain Rates. J Geophys Res Atmospheres (1967) 72:320–2. doi:10.1029/jz072i002p00611

CrossRef Full Text | Google Scholar

9. Hojem J, Cook N. The Design and Constructon of a Triaxial and Polyaxial Cell for Testing Rock Specimens. South Africa Mech.Engr. (1968) 18:57–61.

Google Scholar

10. Mogi K. Fracture and Flow of Rocks under High Triaxial Compression. J Geophys Res (1971) 76:1255–69. doi:10.1029/jb076i005p01255

CrossRef Full Text | Google Scholar

11. Takahashi M, Koide H, Takahashi M. Effect of the Intermediate Principal Stress on Strength and Deformation Behavior of Sedimentary Rocks at the Depth Shallower Than 200m. In: International Symposium on Rock at Great Depth; 1989 August 28–39; PAU (1989).

Google Scholar

12. Haimson BC, Chang C. True Triaxial Strength of the Ktb Amphibolite under Borehole wall Conditions and its Use to Estimate the Maximum Horizontal In Situ Stress. J Geophys Res Atmospheres (2002) 107:ETG 15–1–ETG 15–4. doi:10.1029/2001jb000647

CrossRef Full Text | Google Scholar

13. Feng X-T, Zhang X, Kong R, Wang G. A Novel Mogi Type True Triaxial Testing Apparatus and its Use to Obtain Complete Stress-Strain Curves of Hard Rocks. Rock Mech Rock Eng (2016) 49:1649–62. doi:10.1007/s00603-015-0875-y

CrossRef Full Text | Google Scholar

14. Zhao J, Feng X-T, Zhang X-W, Zhang Y, Zhou Y-Y, Yang C-X. Brittle-ductile Transition and Failure Mechanism of Jinping marble under True Triaxial Compression. Eng Geology (2018) 232:160–70. doi:10.1016/j.enggeo.2017.11.008

CrossRef Full Text | Google Scholar

15. Cai M. Influence of Intermediate Principal Stress on Rock Fracturing and Strength Near Excavation Boundaries-Insight from Numerical Modeling. Int J Rock Mech Mining Sci (2008) 45:763–72. doi:10.1016/j.ijrmms.2007.07.026

CrossRef Full Text | Google Scholar

16. Lu S, Xiao-Chun LI. Analysis of End Friction Effect in True Triaxial Test. Rock Soil Mech (2009) 30:1159–64. doi:10.16285/j.rsm.2009.04.057

CrossRef Full Text | Google Scholar

17. Hudson JA, Tang C. Rock Failure Mechanisms: Explained and Illustrated. London: Taylor & Francis (2010).

Google Scholar

18. Pan P-Z, Feng X-T, Hudson JA. The Influence of the Intermediate Principal Stress on Rock Failure Behaviour: A Numerical Study. Eng Geol (2012) 124:109–18. doi:10.1016/j.enggeo.2011.10.008

CrossRef Full Text | Google Scholar

19. Zhang H, Xie X, Zhang M, Yang G. Damage Constitutive Model of Rock under the True Triaxial Confinement State. Mech Eng (2015) 37:75–8.

Google Scholar

20. Cui C, Zhang L, Li H, Wang D. Three axis Failure Criterion and Constitutive Model of Coal Rock under Horizontal Stress Loading. In: International Conference on Oil and Gas Field Exploration and Development; 2018 September 18–20; Xi'an, China (2018).

Google Scholar

21. Zhu QZ, Shao JF. A Refined Micromechanical Damage-Friction Model with Strength Prediction for Rock-like Materials under Compression. Int J Sol Structures (2015) 60-61:75–83. doi:10.1016/j.ijsolstr.2015.02.005

CrossRef Full Text | Google Scholar

22. Budiansky B, O'connell RJ. Elastic Moduli of a Cracked Solid. Int J Sol Struct (1976) 12:81–97. doi:10.1016/0020-7683(76)90044-5

CrossRef Full Text | Google Scholar

23. Horii H, Nemat-Nasser S. Overall Moduli of Solids with Microcracks: Load-Induced Anisotropy. J Mech Phys Sol (1983) 31:155–71. doi:10.1016/0022-5096(83)90048-0

CrossRef Full Text | Google Scholar

24. Mori T, Tanaka K. Average Stress in Matrix and Average Elastic Energy of Materials with Misfitting Inclusions. Acta Metallurgica (1973) 21:571–4. doi:10.1016/0001-6160(73)90064-3

CrossRef Full Text | Google Scholar

25. Kachanov M. Effective Elastic Properties of Cracked Solids: Critical Review of Some Basic Concepts. Appl Mech Rev (1992) 45:304–35. doi:10.1115/1.3119761

CrossRef Full Text | Google Scholar

26. Zhu Q-Z, Yuan S-S, Shao J-f. Bridging Meso- and Microscopic Anisotropic Unilateral Damage Formulations for Microcracked Solids. Comptes Rendus Mécanique (2017) 345:281–92. doi:10.1016/j.crme.2017.02.003

CrossRef Full Text | Google Scholar

27. Lemaitre J, Dufailly J. Damage Measurements. Eng Fracture Mech (1987) 28:643–61. doi:10.1016/0013-7944(87)90059-2

CrossRef Full Text | Google Scholar

28. Murakami S. Continuum Damage Mechanics: A Continuum Mechanics Approach to the Analysis of Damage and Fracture, Vol. 41. London, New York: Springer Ebooks (2012). 4731–4755. doi:10.1007/978-94-007-2666-6

CrossRef Full Text | Google Scholar

29. Ladeveze P, Ladeveze P. Sur une theorie de l’endommagement anisotrope (1983). Laboratoire de Mechanique et Technologie, Cachan. Technical Report 34.

Google Scholar

30. He Q-C, Curnier A. A More Fundamental Approach to Damaged Elastic Stress-Strain Relations. Int J Sol Struct (1995) 32:1433–57. doi:10.1016/0020-7683(94)00183-w

CrossRef Full Text | Google Scholar

31. Zhu Qz. Strength Prediction of Dry and Saturated Brittle Rocks by Unilateral Damage-Friction Coupling Analyses. Comput Geotechn (2016) 73:16–23. doi:10.1016/j.compgeo.2015.11.015

CrossRef Full Text | Google Scholar

32. Zhao L-Y, Lai Y-M, Shao J-F, Zhang W-L, Zhu Q-Z, Niu F-J. Friction-damage Coupled Models and Macroscopic Strength Criteria for Ice-Saturated Frozen silt with Crack Asperity Variation by a Micromechanical Approach. Eng Geol (2021) 294:106405. doi:10.1016/j.enggeo.2021.106405

CrossRef Full Text | Google Scholar

33. Martin CD, Chandler NA. The progressive fracture of lac du bonnet granite. Int J Rock Mech Mining Sci Geomech Abstr (1994) 31:643–59. doi:10.1016/0148-9062(94)90005-1

CrossRef Full Text | Google Scholar

34. Xue L, Qin S, Sun Q, Wang Y, Lee LM, Li W. A Study on Crack Damage Stress Thresholds of Different Rock Types Based on Uniaxial Compression Tests. Rock Mech Rock Eng (2014) 47:1183–95. doi:10.1007/s00603-013-0479-3

CrossRef Full Text | Google Scholar

35. Zhu QZ, Zhao LY, Shao JF. Analytical and Numerical Analysis of Frictional Damage in Quasi Brittle Materials. J Mech Phys Sol (2016) 92:137–63. doi:10.1016/j.jmps.2016.04.002

CrossRef Full Text | Google Scholar