# A Statistical Damage Constitutive Model of Anisotropic Rock: Development and Validation

Authors: Tenglong Rong; Can Guan; Keliu Liu; Shuai Heng; Wenlong Shen; Ruiyong Mou
Journal: Geofluids
 (2021)

---

## Abstract

The damage constitutive model is of great significance to research the stress-strain relationship and damage evolution of rock under loading in engineering. In order to investigate the effect of anisotropic characteristic on the stress-strain relationship and damage evolution, a statistical damage constitutive model of anisotropic rock under true triaxial condition was developed. In this study, the plane which existed perpendicular to the coordinate axis was extracted from representative volume element (RVE) of rock. The extracted plane was assumed to be composed of abundant mesoscopic elements whose failure strength satisfied the Weibull distribution. According to the number of failure elements on the plane in each direction under loading, the anisotropic damage variable was established based on the proposed concept of areal damage. A statistical damage constitutive model of anisotropic rock was developed by using strain equivalent hypothesis and generalized Hooke constitutive model. Subsequently, the parameters in the anisotropic damage constitutive model were determined by the method of total differential. Thus, the damage evolution of anisotropic rock under various stress conditions can be conveniently evaluated by the anisotropic damage model. The model was validated based on the tests of rocks under the stress conditions of conventional triaxial and true triaxial, respectively. Moreover, for the purpose of studying the influence of parameters on the model, sensitivity analyses of mechanical parameters and model parameters were carried out. The results of statistical damage constitutive clearly demonstrate the stress-strain and damage evolution of anisotropic rock under various stress conditions.


---

## Body

## 1. Introduction

The stress-strain relationship of rock affected by damage under loading has always been the focus of rock mechanics. It is of great significance to accurately predict the mechanical state and deformation of rock for geotechnical [1], coal mining [2, 3], shale gas exploitation [4], etc. Due to the long-term geological effects in history, rock is a heterogeneous material which includes flaws and geological structure, i.e. pores, voids, defects, joints, fractures, and bedding. Therefore, anisotropy is an obvious characteristic of rock. It is very valuable to investigate the damage and deformation of anisotropic rock under loading. The stress-strain relationship and damage evolution of anisotropic rock have been studied by laboratory tests in previous researches [4–7].The damage evolution and stress-strain relationship can also be studied by means of theoretical analysis. The damage constitutive model considering anisotropic characteristics is the theoretical basis of analyzing the anisotropic rock deformation. Based on the previous researches, there are two main theoretical approaches to establish the constitutive relationship of rock under the influence of damage. One approach to establish constitutive model is continuum damage mechanics (CDM). In this approach, rock is regarded as continuous medium. The relationship between damage variable and loading condition is established according to CDM. The modeling approach based on CDM is developed very early. Some of the models based on CDM involve the anisotropic characteristic [8–12], but the anisotropic damage model based on this approach is complicated.The other approach to establish constitutive relationship of rock is the statistical damage model (SDM), which is according to the statistics and continuum damage mechanics. In the process of modeling, rock is regarded as clusters composed of abundant mesoscopic elements. The mechanical parameters of mesoscopic elements satisfy a certain mathematical function distribution. Then, the number of destructive mesoscopic elements in rock under loading is calculated by the statistical method. Consequently, damage variable is put forward according to the number of destructive mesoscopic elements to evaluate the damage evolution. Krajcinovic and Silva first proposed a simple statistical damage model to study the material damage evolution under uniaxial tension [13]. Tang and Xu assumed that rock was composed of a large number of mesoscopic elements, which were elastic before failure and useless after failure. The failure strength of mesoscopic elements was based on the normal statistical distribution. Then, the stress-strain of rock under uniaxial condition was analyzed by the continuum damage mechanics [14]. Subsequently, the theory of SDM was improved by Tang and coworkers [15–18]. Cao and coworkers [19–21] made many attempts to improve SDM from different aspects. The attempts included different statistical distribution functions (Weibull distribution and normal distribution) and different strength criterions (Mohr-Coulomb and Drucker-Prager). Meanwhile, the following issues were taken into consideration in these attempts: damage threshold, initial compression, strain softening, and residual strength. Considering scale effect in rock shearing, some statistical damage constitutive models of rock under shear stress were also proposed [22, 23]. In order to research the stress-strain of rock under damage effect more accurately, many mathematical and mechanical theories have been introduced into the SDM [20, 24–27].In addition to rock damage under loading (compression, tension and shear), the SDM has also been widely used in research of rock damage in the multiphysical environment [28–32]. The factors of porosity, temperature, or freeze-thaw cycles were introduced into the damage variables of SDM under the multiphysical conditions. Moreover, in order to investigate the effect of joints on damage, the factors of crack propagation length, joint friction effect, and joint orientation were drawn into the SDM [33, 34]. There are many statistical damage models that can be used to analyze the progressive damage process of rock under various actions, such as load (compression, tension and shear), high temperature, freeze-thaw, chemistry, and joint. It can be seen that most of the presented SDM are based on isotropic properties. Therefore, the SDM of rock which is based on anisotropic properties needs to be further studied. It is great valuable to develop a SDM which can represent the anisotropic characteristic of rock.In this study, based on the assumption that any plane of rock was composed of mesoscopic elements whose strength conformed to the Weibull distribution. The concept of areal damage was proposed to establish the expression of anisotropic damage variable. Then, a statistical damage constitutive model of anisotropic rock (SDAR model) was developed by using strain equivalent hypothesis. Subsequently, the determination of parameters in the anisotropic damage constitutive model was conducted by the method of total differential. Finally, the model was validated based on the tests of two kinds of rocks under different loading stress conditions. One was transversely isotropic shale under conventional triaxial stress condition. And the other was orthogonal anisotropic coal under true triaxial stress condition. In addition, the sensitivity study of parameters was conducted by the developed damage model.


## 2. Development of the Statistical Damage Constitutive Model of Anisotropic Rock

### 2.1. Anisotropic Damage Variable

Rock damage refers to the progressive deterioration of mechanical properties caused by the defect accumulation in rock under the action of stress and other factors. In the research on creep behavior of metal under uniaxial tension, Kachanov found that the propagation of microdefects was the main reason for creep damage in material [35]. Consequently, the definition of “continuity” was proposed to quantitatively show the damage behavior of material under loading in Eq. (1), i.e.,
(1)ψ=S~S,where ψ is the continuity of material, S is the initial total area of the material,and S~ is the actual bearing area, which is the result of deducting the failure area from the initial total area.Rabotnov [36] further proposed the concept of “damage factor ω,” which is presented as Eq. (2). Therefore, the damage factor is the ratio of the failure area which is invalid caused by the defects to the initial total area.
(2)ω=S−S~S.Damage variable is an important index to measure the damage degree of rock. The evolution of damage variable has an important relationship with the microdefects and fracture of rock. For the purpose of investigating the damage evolution of anisotropic rock in orthogonal directions, representative volume element (RVE) is taken for analysis. One plane which is perpendicular to the coordinatei-axis is extracted from the RVE (Figure 1(a)). The extracted plane is assumed to be composed of abundant initial mesosopic elements. The number of initial mesosopic elements is Ni0. The damage evolution of the plane can be evaluated under three-dimensional stress condition. The load in the direction of the i-axis is axial compression, and the load in the other two orthogonal directions is confining pressure. Meanwhile, the load in three directions can be changed independently. The plane can also be applied to true triaxial stress condition. Under the effect of three-dimensional stress, the elements are progressive failure, which results in the rock damage. The number of damaged elements in the plane under loading stress is Nid (Figure 1(b)).Figure 1
The procedures of proposing the areal damage. (a) Extracting the plane perpendicular to thei-axis from RVE. (b) The damaged elements in the plane under loading stress. (c) The failure probability density function of the mesoscopic elements.In this way, the area of rock plane under loading stress can be divided into two parts: damaged area and undamaged area. The undamaged area is regarded as the effective area which supports the loading. According to the concept of damage factor proposed by Rabotnov [36], the definition of areal damage in different directions which is based on the ratio of the failure element number (Nid) to the initial element number (Ni0) in the extracted plane is given. Because the notation of the i-axis can traverse the coordinate axis of x, y, and z, the extracted plane which is perpendicular to the coordinate i-axis can present the orthogonal anisotropic properties of rocks. Then, the damage evolution of rock in orthogonal directions can be analyzed by the areal damage in orthogonal directions, because the planes are extracted perpendicular to the orthogonal coordinate axes. The areal damage variable Di of the plane which is perpendicular to the coordinate i-axis can be given as
(3)Di=NidNi0.In recent years, on the basis of the assumption that the material was composed of elastic mesoscopic elements and the material properties satisfied the Weibull distribution, some statistical damage constitutive models were established [15, 20, 32]. It is also assumed that the strength of all elements in the plane of each direction are adopted as Weibull distribution in this study. The element would break and fail after the stress reaches to the peak strength. Therefore, the failure probability density function of the mesoscopic elements which are distributed in the plane perpendicular to the i-axis PFi in Figure 1(c) can be expressed as
(4)PFi=wiFi0FiFi0wi−1−FiFi0wi,where Fi is the strength of mesoscopic element, and wi and Fi0 are the parameters of the Weibull distribution function.When the element reaches to any stress level in loading procedure, the number of failure elements in the plane can be calculated as(5)Nid=Ni0∫0FiwiFi0FiFi0wi−1−FiFi0wi=Ni01−exp−FiFi0wi.Combining Eqs. (3) and (5), the areal damage variable of the plane which is perpendicular to the i-axis can be given by
(6)Di=1−exp−FiFi0wi.


### 2.2. Damage Constitutive Model of the Anisotropic Rock

There are many original defects in rock and some new microcracks will be observed under the damage effect. During the process of loading, the load can be only supported by the effective area without damage. Consequently, the concept of effective stress is put forward. Effective stress refers to the increased actual stress in the material, which is increased due to the decrease of the actual bearing area caused by damage and failure. According to this phenomenon, the strain equivalent hypothesis was proposed by Lemaitre et al. [37]. The constitutive model of strain equivalent hypothesis can be expressed as
(7)σ∗=σ/I−Di=Cε/I−Di,whereσ∗ is the effective stress matrix, σ is the apparent stress matrix, C is the elastic matrix of the material, ε is the strain matrix, and I is the identity matrix.In the laboratory triaxial compression test, the three principal stresses can be measured by the testing machine. The generalized Hooke constitutive model can be transformed to(8)εi=1Eiσi∗−∑i≠jνijEjσj∗,where εi is the component of three principal strains, and Ei and νij are the elastic modulus and Poisson’s ratio in different directions, respectively. The subscript i represents the direction of the coordinate i-axis. Furthermore, i,j,k=1,2,3 and i≠j≠k. The subscript presents the traversal procedure, but they do not comply with the summation convention.Substituting Eq. (8) into (7), the three dimensional elastic constitutive equation which includes damage effect can be obtained as
(9)σi=Eiεi1−Di+∑i≠jEiνijEjσj.Substituting Eq. (6) into (9), a new statistical damage constitutive model of anisotropic rock (SDAR model) is developed as
(10)σi=Eiεiexp−FiFi0wi+∑i≠jEiνijEjσj.The Drucker-Prager failure criterion is widely used in the researches of rock mechanics [38–40]. In view of this, the Drucker Prager criterion was introduced as the failure criterion of mesoscopic elements in each plane perpendicular to the coordinate axis. Then, the failure strength of element (Fi) in the plane which is perpendicular to the coordinate i-axis can be given by
(11)Fi=αdpI1+J2,where αdp refers to the parameter of material, I1 is the first invariant of stress, and J2 is the second invariant of stress deviator.The three parameters of the Drucker-Prager failure criterion can be calculated by(12)αdp=sinφ/9+3sin2φ,I1=σ1∗+σ2∗+σ3∗,J2=16σ1∗−σ2∗2+σ2∗−σ3∗2+σ1∗−σ3∗2,where φ refers to the friction angle, and σ1∗, σ2∗, and σ3∗ are the effective stress in three directions.


### 2.3. Determination of Parameters in the Anisotropic Damage Constitutive Model

Eqs. (6) and (10) denote the evolution equation of the damage and SDAR model, respectively. There are two important parameters from the Weibull distribution function in the equations, i.e., wi and Fi0, but the values of wi and Fi0 are not determined. The method of determining the two parameters needs to be further investigated.The stressσi in the SDAR model is differentiated from Eq. (10), i.e.,
(13)dσi=Ai1dεi+Ai2dFi+Ai3dwi+Ai4dFi0+∑i≠jEiνijEjσj,where Ai1=∂σi/∂εi, Ai2=∂σi/∂Fi, Ai3=∂σi/∂wi, and Ai4=∂σi/∂Fi0.The failure strength of element (Fi) can be regarded as function of σ1, σ2, σ3, and εi. The subscript i of εi corresponds to that of Fi. Then, the following Eq. (14) can be obtained by differential of Fi.
(14)dFi=∑Fijdσj+Fi∗dεi.The parameters ofw and F0 from the Weibull distribution function in the development of statistical damage constitutive model were assumed to be a function of confining pressure in previous research [41]. And the assumption was validated by comparing with test results. In order to investigate the anisotropic damage of rock, the parameters of wi and Fi0 in each direction are assumed to be functions of loading stress in the other two directions (σj, σk) by analogy with that assumption. But σj and σk are not equal in true triaxial compression, the parameters are considered as functions of its average value. By differentiation, we can get
(15)dwi=wi∗dσ¯i,dFi0=Fi0∗dσ¯i,where σ¯i=1/2σj+σk.Substituting Eqs. (14) and (15) into Eq. (13) yields
(16)∑Uijdσj=−Xidεi,where
(17)Uij=A12F11−1A12F12+12A13w1∗+12A14F10∗+E1ν12E2A12F13+12A13w1∗+12A14F10∗+E1ν13E3A22F21+12A23w2∗+12A24F20∗+E2ν21E1A22F22−1A22F23+12A23w2∗+12A24F20∗+E2ν23E3A32F31+12A33w3∗+12A34F30∗+E3ν31E1A32F32+12A33w3∗+12A34F30∗+E3ν32E2A32F33−1,Xi=A11+A12F1∗A21+A22F2∗A31+A32F3∗.Solving Eq. (16) based on the Kramer’s law yields
(18)dσi=−1i∑XjΔijΔdεj,where Δ=Uij=U11U12U13U21U22U23U31U32U33, and Δij is the remainder of Δ, i.e.,
(19)Δ11=U22U23U32U33Δ12=U21U23U31U33Δ13=U21U22U31U32,Δ21=U12U13U32U33Δ22=U11U13U31U33Δ23=U11U12U31U32,Δ31=U12U13U22U23Δ32=U11U13U21U23Δ33=U11U12U21U22.Based on the generalized Hooke constitutive model, the principle stressσi can be expressed as a function of three principal strains ε1, ε2, and ε3. Thus, the total differential of σi can be obtained as
(20)dσi=∑∂σi∂εjdεj.In the true triaxial stress test, the stress-strain curve of rock in each direction to be solved should reach the following boundary conditions at the peak point, i.e.,(21)εi=εicp,σi=σicp,dσidεiεi=εicp=0,where εicp refers to the peak strain in the loading direction along the i-axis, and σicp refers to the peak stress in the loading direction along the i-axis.Comparing Eqs. (18) and (20), it can be known that the coefficients of items at the same position should be the same. Subsequently, Eqs. (18) and (20) are substituted into Eq. (21), and the Eq. (22) can be obtained.
(22)∂σi∂εi=−1iXiΔiiΔ=0.For the reason thatΔii≠0, the equation that Xi=0 can be obtained. By solving Xi=0, the calculated expressions of wi and Fi0 can be derived as
(23)wi=1lnEiεicp/σicp−∑i≠jEiνij/Ejσj,(24)Fi0=Ficwi1/wi.Eqs. (23) and (24) denote the determination of parameters in the SDAR model, i.e., wi and Fi0. When the material properties of rock are assumed to be isotropic, the elastic modulus and Poisson’s ratio in all directions are equal, respectively. In the conventional triaxial test, the two principal stresses in the horizontal direction are equal (σ2=σ3≠0). On this basis, for the conventional triaxial compression test of isotropic rock, the calculated expressions of Eqs. (23) and (24) can be degenerated to
(25)w=1lnEεcp/σcp−2νσ3,(26)F0=Fcw1/w,where σcp is the axial peak stress, and σ3 is the confining pressure.For the uniaxial triaxial compression condition of isotropic rock, the calculated expressions of parameters illustrated as Eqs. (23) and (24) can be further degenerated for lack of confining pressure (σ2=σ3=0). At this time, Eq. (25) for solving parameter w is degenerated to Eq. (27), but Eq. (26) for solving parameter F0 remains unchanged.
(27)w=1lnEεcp/σcp.Equations (25)–(27) are the calculated expressions of w and F0 obtained by the degradation of the calculation expressions of wi and Fi0 for isotropic rock under conventional triaxial and uniaxial compression, respectively. The calculated expressions of w and F0 under conventional triaxial compression condition of isotropic rock (Eqs. (25) and (26)) are consistent with those derived by Wang et al. [42]. Furthermore, the expressions under uniaxial compression condition of isotropic rock (Eqs. (26) and (27)) are the same as the research results obtained by Xie and Zhao [43].


## 2.1. Anisotropic Damage Variable

Rock damage refers to the progressive deterioration of mechanical properties caused by the defect accumulation in rock under the action of stress and other factors. In the research on creep behavior of metal under uniaxial tension, Kachanov found that the propagation of microdefects was the main reason for creep damage in material [35]. Consequently, the definition of “continuity” was proposed to quantitatively show the damage behavior of material under loading in Eq. (1), i.e.,
(1)ψ=S~S,where ψ is the continuity of material, S is the initial total area of the material,and S~ is the actual bearing area, which is the result of deducting the failure area from the initial total area.Rabotnov [36] further proposed the concept of “damage factor ω,” which is presented as Eq. (2). Therefore, the damage factor is the ratio of the failure area which is invalid caused by the defects to the initial total area.
(2)ω=S−S~S.Damage variable is an important index to measure the damage degree of rock. The evolution of damage variable has an important relationship with the microdefects and fracture of rock. For the purpose of investigating the damage evolution of anisotropic rock in orthogonal directions, representative volume element (RVE) is taken for analysis. One plane which is perpendicular to the coordinatei-axis is extracted from the RVE (Figure 1(a)). The extracted plane is assumed to be composed of abundant initial mesosopic elements. The number of initial mesosopic elements is Ni0. The damage evolution of the plane can be evaluated under three-dimensional stress condition. The load in the direction of the i-axis is axial compression, and the load in the other two orthogonal directions is confining pressure. Meanwhile, the load in three directions can be changed independently. The plane can also be applied to true triaxial stress condition. Under the effect of three-dimensional stress, the elements are progressive failure, which results in the rock damage. The number of damaged elements in the plane under loading stress is Nid (Figure 1(b)).Figure 1
The procedures of proposing the areal damage. (a) Extracting the plane perpendicular to thei-axis from RVE. (b) The damaged elements in the plane under loading stress. (c) The failure probability density function of the mesoscopic elements.In this way, the area of rock plane under loading stress can be divided into two parts: damaged area and undamaged area. The undamaged area is regarded as the effective area which supports the loading. According to the concept of damage factor proposed by Rabotnov [36], the definition of areal damage in different directions which is based on the ratio of the failure element number (Nid) to the initial element number (Ni0) in the extracted plane is given. Because the notation of the i-axis can traverse the coordinate axis of x, y, and z, the extracted plane which is perpendicular to the coordinate i-axis can present the orthogonal anisotropic properties of rocks. Then, the damage evolution of rock in orthogonal directions can be analyzed by the areal damage in orthogonal directions, because the planes are extracted perpendicular to the orthogonal coordinate axes. The areal damage variable Di of the plane which is perpendicular to the coordinate i-axis can be given as
(3)Di=NidNi0.In recent years, on the basis of the assumption that the material was composed of elastic mesoscopic elements and the material properties satisfied the Weibull distribution, some statistical damage constitutive models were established [15, 20, 32]. It is also assumed that the strength of all elements in the plane of each direction are adopted as Weibull distribution in this study. The element would break and fail after the stress reaches to the peak strength. Therefore, the failure probability density function of the mesoscopic elements which are distributed in the plane perpendicular to the i-axis PFi in Figure 1(c) can be expressed as
(4)PFi=wiFi0FiFi0wi−1−FiFi0wi,where Fi is the strength of mesoscopic element, and wi and Fi0 are the parameters of the Weibull distribution function.When the element reaches to any stress level in loading procedure, the number of failure elements in the plane can be calculated as(5)Nid=Ni0∫0FiwiFi0FiFi0wi−1−FiFi0wi=Ni01−exp−FiFi0wi.Combining Eqs. (3) and (5), the areal damage variable of the plane which is perpendicular to the i-axis can be given by
(6)Di=1−exp−FiFi0wi.


## 2.2. Damage Constitutive Model of the Anisotropic Rock

There are many original defects in rock and some new microcracks will be observed under the damage effect. During the process of loading, the load can be only supported by the effective area without damage. Consequently, the concept of effective stress is put forward. Effective stress refers to the increased actual stress in the material, which is increased due to the decrease of the actual bearing area caused by damage and failure. According to this phenomenon, the strain equivalent hypothesis was proposed by Lemaitre et al. [37]. The constitutive model of strain equivalent hypothesis can be expressed as
(7)σ∗=σ/I−Di=Cε/I−Di,whereσ∗ is the effective stress matrix, σ is the apparent stress matrix, C is the elastic matrix of the material, ε is the strain matrix, and I is the identity matrix.In the laboratory triaxial compression test, the three principal stresses can be measured by the testing machine. The generalized Hooke constitutive model can be transformed to(8)εi=1Eiσi∗−∑i≠jνijEjσj∗,where εi is the component of three principal strains, and Ei and νij are the elastic modulus and Poisson’s ratio in different directions, respectively. The subscript i represents the direction of the coordinate i-axis. Furthermore, i,j,k=1,2,3 and i≠j≠k. The subscript presents the traversal procedure, but they do not comply with the summation convention.Substituting Eq. (8) into (7), the three dimensional elastic constitutive equation which includes damage effect can be obtained as
(9)σi=Eiεi1−Di+∑i≠jEiνijEjσj.Substituting Eq. (6) into (9), a new statistical damage constitutive model of anisotropic rock (SDAR model) is developed as
(10)σi=Eiεiexp−FiFi0wi+∑i≠jEiνijEjσj.The Drucker-Prager failure criterion is widely used in the researches of rock mechanics [38–40]. In view of this, the Drucker Prager criterion was introduced as the failure criterion of mesoscopic elements in each plane perpendicular to the coordinate axis. Then, the failure strength of element (Fi) in the plane which is perpendicular to the coordinate i-axis can be given by
(11)Fi=αdpI1+J2,where αdp refers to the parameter of material, I1 is the first invariant of stress, and J2 is the second invariant of stress deviator.The three parameters of the Drucker-Prager failure criterion can be calculated by(12)αdp=sinφ/9+3sin2φ,I1=σ1∗+σ2∗+σ3∗,J2=16σ1∗−σ2∗2+σ2∗−σ3∗2+σ1∗−σ3∗2,where φ refers to the friction angle, and σ1∗, σ2∗, and σ3∗ are the effective stress in three directions.


## 2.3. Determination of Parameters in the Anisotropic Damage Constitutive Model

Eqs. (6) and (10) denote the evolution equation of the damage and SDAR model, respectively. There are two important parameters from the Weibull distribution function in the equations, i.e., wi and Fi0, but the values of wi and Fi0 are not determined. The method of determining the two parameters needs to be further investigated.The stressσi in the SDAR model is differentiated from Eq. (10), i.e.,
(13)dσi=Ai1dεi+Ai2dFi+Ai3dwi+Ai4dFi0+∑i≠jEiνijEjσj,where Ai1=∂σi/∂εi, Ai2=∂σi/∂Fi, Ai3=∂σi/∂wi, and Ai4=∂σi/∂Fi0.The failure strength of element (Fi) can be regarded as function of σ1, σ2, σ3, and εi. The subscript i of εi corresponds to that of Fi. Then, the following Eq. (14) can be obtained by differential of Fi.
(14)dFi=∑Fijdσj+Fi∗dεi.The parameters ofw and F0 from the Weibull distribution function in the development of statistical damage constitutive model were assumed to be a function of confining pressure in previous research [41]. And the assumption was validated by comparing with test results. In order to investigate the anisotropic damage of rock, the parameters of wi and Fi0 in each direction are assumed to be functions of loading stress in the other two directions (σj, σk) by analogy with that assumption. But σj and σk are not equal in true triaxial compression, the parameters are considered as functions of its average value. By differentiation, we can get
(15)dwi=wi∗dσ¯i,dFi0=Fi0∗dσ¯i,where σ¯i=1/2σj+σk.Substituting Eqs. (14) and (15) into Eq. (13) yields
(16)∑Uijdσj=−Xidεi,where
(17)Uij=A12F11−1A12F12+12A13w1∗+12A14F10∗+E1ν12E2A12F13+12A13w1∗+12A14F10∗+E1ν13E3A22F21+12A23w2∗+12A24F20∗+E2ν21E1A22F22−1A22F23+12A23w2∗+12A24F20∗+E2ν23E3A32F31+12A33w3∗+12A34F30∗+E3ν31E1A32F32+12A33w3∗+12A34F30∗+E3ν32E2A32F33−1,Xi=A11+A12F1∗A21+A22F2∗A31+A32F3∗.Solving Eq. (16) based on the Kramer’s law yields
(18)dσi=−1i∑XjΔijΔdεj,where Δ=Uij=U11U12U13U21U22U23U31U32U33, and Δij is the remainder of Δ, i.e.,
(19)Δ11=U22U23U32U33Δ12=U21U23U31U33Δ13=U21U22U31U32,Δ21=U12U13U32U33Δ22=U11U13U31U33Δ23=U11U12U31U32,Δ31=U12U13U22U23Δ32=U11U13U21U23Δ33=U11U12U21U22.Based on the generalized Hooke constitutive model, the principle stressσi can be expressed as a function of three principal strains ε1, ε2, and ε3. Thus, the total differential of σi can be obtained as
(20)dσi=∑∂σi∂εjdεj.In the true triaxial stress test, the stress-strain curve of rock in each direction to be solved should reach the following boundary conditions at the peak point, i.e.,(21)εi=εicp,σi=σicp,dσidεiεi=εicp=0,where εicp refers to the peak strain in the loading direction along the i-axis, and σicp refers to the peak stress in the loading direction along the i-axis.Comparing Eqs. (18) and (20), it can be known that the coefficients of items at the same position should be the same. Subsequently, Eqs. (18) and (20) are substituted into Eq. (21), and the Eq. (22) can be obtained.
(22)∂σi∂εi=−1iXiΔiiΔ=0.For the reason thatΔii≠0, the equation that Xi=0 can be obtained. By solving Xi=0, the calculated expressions of wi and Fi0 can be derived as
(23)wi=1lnEiεicp/σicp−∑i≠jEiνij/Ejσj,(24)Fi0=Ficwi1/wi.Eqs. (23) and (24) denote the determination of parameters in the SDAR model, i.e., wi and Fi0. When the material properties of rock are assumed to be isotropic, the elastic modulus and Poisson’s ratio in all directions are equal, respectively. In the conventional triaxial test, the two principal stresses in the horizontal direction are equal (σ2=σ3≠0). On this basis, for the conventional triaxial compression test of isotropic rock, the calculated expressions of Eqs. (23) and (24) can be degenerated to
(25)w=1lnEεcp/σcp−2νσ3,(26)F0=Fcw1/w,where σcp is the axial peak stress, and σ3 is the confining pressure.For the uniaxial triaxial compression condition of isotropic rock, the calculated expressions of parameters illustrated as Eqs. (23) and (24) can be further degenerated for lack of confining pressure (σ2=σ3=0). At this time, Eq. (25) for solving parameter w is degenerated to Eq. (27), but Eq. (26) for solving parameter F0 remains unchanged.
(27)w=1lnEεcp/σcp.Equations (25)–(27) are the calculated expressions of w and F0 obtained by the degradation of the calculation expressions of wi and Fi0 for isotropic rock under conventional triaxial and uniaxial compression, respectively. The calculated expressions of w and F0 under conventional triaxial compression condition of isotropic rock (Eqs. (25) and (26)) are consistent with those derived by Wang et al. [42]. Furthermore, the expressions under uniaxial compression condition of isotropic rock (Eqs. (26) and (27)) are the same as the research results obtained by Xie and Zhao [43].


## 3. Model Validation

### 3.1. Stress-Strain Relationship

#### 3.1.1. Transversely Isotropic Shale under Conventional Triaxial Stress Condition

As a kind of sedimentary rock, shale is observed to be transversely isotropic. Shale has the same mechanical properties in the horizontal direction, but the mechanical properties in the horizontal direction are different from those in the vertical direction. Heng [44] used shale which is from the Longmaxi Formation in Chongqing to drill tested samples along two orthognoal directions (Figure 2). One is parallel to bedding (drilling angle 0°). The other is perpendicular to bedding (drilling angle 90°). Shale samples drilled from different directions were tested under conventional triaxial compression with confining pressures of 10, 20, and 30 MPa, respectively. Then, the stress-strain relationships of shale which were drilled in orthogonal directions under different confining pressures that are obtained. The parameters used for validation of shale are shown in Table 1.Figure 2
Directional coring diagram of shale.Table 1
Mechanical parameters of shale in different directions [44].
Confining pressure (MPa)Elastic modulus parallel to bedding (GPa)Elastic modulus perpendicular to bedding (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)1022.9142.940.3840.32530.32022.7752.510.4160.26330.33024.0456.450.4220.25830.3Comparing the theoretical results of the SDAR model and test data of shale on stress-strain in Figure3, it can be seen that the theoretical results are very close to the test data in the same drilling direction under confining pressures of 10, 20, and 30 MPa, respectively. The comparison results show that all the stress-strain curves of shale under conventional triaxial compression can well be evaluated by the SDAR model. In the prepeak zone, the curve of theoretical calculation is in good agreement with that of test data. After the peak stress point, the curve of theoretical results present strain softening in the postpeak region. The presentation of strain softening in the SDAR model is also consistent with the mechanical behavior of shale.Figure 3
Stress-strain curves of shale under different confining pressures: (a) confining pressure 10 MPa, (b) confining pressure 20 MPa, and (c) confining pressure 30 MPa.
(a)(b)(c)Comparing the peak strength of shale in the same drilling direction under each confining pressure (Figure3), the peak strength of samples in the same drilling direction increases with the increase of confining pressure. Moreover, the peak strength under different confining pressures can be accurately calculated by the SDAR model. Then, comparing the results of samples drilled in different directions under the same confining pressure, it can be found that the peak strength of sample drilled in the direction of parallel to bedding is always bigger than that of perpendicular to bedding. Meanwhile, the peak strain of sample drilled in the direction of parallel to bedding is smaller than that of perpendicular to bedding. The stress-strain of transversely isotropic shale under conventional triaxial stress can be well calculated by the SDAR model.


#### 3.1.2. Orthogonal Anisotropic Coal under True Triaxial Stress Condition

Under the influence of sedimentation, there is an obvious approximate horizontal bedding plane structure in coal. In addition, there are a lot of cleats in coal, which can be divided into face cleat and butt cleat [45]. These two kinds of cleats are approximately orthogonal to each other and perpendicular to the bedding plane (Figure 4(a)) [46, 47]. Due to the existence of cleat and bedding (referred to as fractures), coal is generally divided into matrix and fracture. The coal matrix is also usually simplified to shape of cubes. Therefore, the mechanical properties of coal are considered to be nearly orthotropic [48].Figure 4
Stress loading path and loading direction variations in the true triaxial tests. (a) Stress loading path. (b) Variations of stress loading direction. [49].
(a)(b)According to the geological structure characteristics of coal seam, the raw cubic coal samples were prepared in the orthogonal fracture directions in Liu’s research [49]. The orthogonal fracture directions were bedding plane direction (BD), face cleat plane direction (FD), and butt cleat plane direction (UD), respectively. Then, the true triaxial loading test was carried out with the cubic coal samples in two true triaxial stress paths. The two true triaxial stress paths were conducted as follows (Figure 4(a)). Firstly, the coal sample was loaded to the hydrostatic pressure (σ1=σ2=σ3) of 10 MPa. Next, σ3 was kept constant while the σ1 and σ2 were increased to the design values of σ2 simultaneously. The design values of σ2 were 20 and 30 MPa, respectively. Finally, the σ2 and σ3 were kept constant, but the σ1 was increased until the coal sample failure.In this way, the complete stress-strain of coal samples under the true triaxial stress conditon was tested. The loading directions of three principal stresses (σ1, σ2, σ3) were parallel or perpendicular to the original fractures (bedding plane, face cleat, and butt cleat) in the tests. During the whole tests, the directions of the three principal stresses were kept unchanged, but the directions of original fractures in cubic coal sample were changed in turn (Figure 4(b)). Based on the true triaxial test results under different stress paths, the developed SDAR model is validated. The parameters used for validation of coal are illustrated in Table 2.Table 2
Mechanical parameters of coal in different directions [49, 50].
Stress loading pathElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)Path 1 (σ2=20MPa, σ3=10MPa)4.113.283.240.140.2228Path 2 (σ2=30MPa, σ3=10MPa,)4.553.672.870.140.2228The theoretical results of the SDAR model and test data of coal on stress-strain are shown in Figure5. By comparing the stress-strain in each stress path, the peak strength of coal samples prepared in three directions is in the same order. The peak strength of coal sample along the BD is the largest. And the peak strength of coal sample along UD is the smallest. Moreover, the peak strain also decreased in the order of BD, FD, and UD.Figure 5
Stress-strain of coal in orthogonal directions in different stress loading paths. (a) Path 1 (σ2=20MPa, σ3=10MPa). (b) Path 2 (σ2=30MPa, σ3=10MPa).
(a)(b)Compared with the test data, the results of the SDAR model can basically represent the stress-strain relationship of coal in different directions under true triaxial condition. In the phase of prepeak, the theoretical results are in good agreement with the test data. The peak stress of the SDAR model is very close to the experimental value of coal in each case. In the postpeak phase, all the theoretical curves obviously present variation of strain softening, which is the same as the test data. All the values of the SDAR model decrease rapidly after the peak stress point. The corresponding tested results also show a decrease trend, but the decrease rate is relatively slow. There is a certain deviation between the SDAR model and the test results in the postpeak phase. In the postpeak phase, there are many reasons for this deviation, including the properties of rock, test conditions, loading procedure, and failure mode of rock [20]. On the whole, stress-strain relationship of orthotropic coal under true triaxial loading can be well evaluated by the SDAR model.


### 3.2. Evolution of Damage Variable

As a heterogeneous geological material, rock has many original defects. Under the influence of various stresses, the defects continuously emerge and expand, which result in damage accumulation in rock. Damage accumulation is the main reason for progressive failure of rock. It is also the important factor determining the stress-strain relationship of rock. It is of great significance to investigate the rock damage evolution under loading for developing the constitutive model.Based on Liu’s research [49], the damage evolution of orthotropic coal samples drilled in different directions under true triaxial loading is analyzed by the SDAR model. Figure 6 demonstrates the evolution of coal samples prepared in each direction (BD, FD, and UD) under different stress paths. It can be seen that when the axial strain is less than a certain value, the value of damage variable in all directions is 0. In other words, there is no damage in the early phase of loading. In the process of test, the axial strain increases continuously with the increase of loading stress. After the peak stress point, the stress-strain curves present the strain softening with the decrease of stress. But during the whole loading process, both the damage variable and the axial strain are monotonically increase. The value of damage variable increases from 0 to1.Figure 6
Damage evolution of anisotropic coal under different stress loading paths. (a) Path 1 (σ2=20MPa, σ3=10MPa). (b) Path 2 (σ2=30MPa, σ3=10MPa).
(a)(b)In stress path 1 (σ2=20MPa, σ3=10MPa), the damage evolution of coal samples loading along the FD and UD is similar. The initial damage first emerges in the coal sample loaded in FD, but the damage variable of coal sample loaded in UD first reaches 1. The initial damage of coal sample loaded in BD emerges the latest. Meanwhile, the damage accumulation rate of coal loaded in BD is the slowest. Therefore, the increase rate of damage variable in BD is the slowest, and the value of axial strain corresponding to damage variable reaching 1 is the largest.In stress path 2 (σ2=30MPa, σ3=10MPa), the initial damage first emerges in the coal sample loaded in UD, but the damage variable of coal sample loaded in FD first reaches 1. The development of damage in the BD is still the slowest. The SDAR model established in this study can well evaluate the damage evolution of anisotropic rock. And the SDAR model also provides a theoretical analysis method for the damage evolution of anisotropic rock under true triaxial condition.


## 3.1. Stress-Strain Relationship

### 3.1.1. Transversely Isotropic Shale under Conventional Triaxial Stress Condition

As a kind of sedimentary rock, shale is observed to be transversely isotropic. Shale has the same mechanical properties in the horizontal direction, but the mechanical properties in the horizontal direction are different from those in the vertical direction. Heng [44] used shale which is from the Longmaxi Formation in Chongqing to drill tested samples along two orthognoal directions (Figure 2). One is parallel to bedding (drilling angle 0°). The other is perpendicular to bedding (drilling angle 90°). Shale samples drilled from different directions were tested under conventional triaxial compression with confining pressures of 10, 20, and 30 MPa, respectively. Then, the stress-strain relationships of shale which were drilled in orthogonal directions under different confining pressures that are obtained. The parameters used for validation of shale are shown in Table 1.Figure 2
Directional coring diagram of shale.Table 1
Mechanical parameters of shale in different directions [44].
Confining pressure (MPa)Elastic modulus parallel to bedding (GPa)Elastic modulus perpendicular to bedding (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)1022.9142.940.3840.32530.32022.7752.510.4160.26330.33024.0456.450.4220.25830.3Comparing the theoretical results of the SDAR model and test data of shale on stress-strain in Figure3, it can be seen that the theoretical results are very close to the test data in the same drilling direction under confining pressures of 10, 20, and 30 MPa, respectively. The comparison results show that all the stress-strain curves of shale under conventional triaxial compression can well be evaluated by the SDAR model. In the prepeak zone, the curve of theoretical calculation is in good agreement with that of test data. After the peak stress point, the curve of theoretical results present strain softening in the postpeak region. The presentation of strain softening in the SDAR model is also consistent with the mechanical behavior of shale.Figure 3
Stress-strain curves of shale under different confining pressures: (a) confining pressure 10 MPa, (b) confining pressure 20 MPa, and (c) confining pressure 30 MPa.
(a)(b)(c)Comparing the peak strength of shale in the same drilling direction under each confining pressure (Figure3), the peak strength of samples in the same drilling direction increases with the increase of confining pressure. Moreover, the peak strength under different confining pressures can be accurately calculated by the SDAR model. Then, comparing the results of samples drilled in different directions under the same confining pressure, it can be found that the peak strength of sample drilled in the direction of parallel to bedding is always bigger than that of perpendicular to bedding. Meanwhile, the peak strain of sample drilled in the direction of parallel to bedding is smaller than that of perpendicular to bedding. The stress-strain of transversely isotropic shale under conventional triaxial stress can be well calculated by the SDAR model.


### 3.1.2. Orthogonal Anisotropic Coal under True Triaxial Stress Condition

Under the influence of sedimentation, there is an obvious approximate horizontal bedding plane structure in coal. In addition, there are a lot of cleats in coal, which can be divided into face cleat and butt cleat [45]. These two kinds of cleats are approximately orthogonal to each other and perpendicular to the bedding plane (Figure 4(a)) [46, 47]. Due to the existence of cleat and bedding (referred to as fractures), coal is generally divided into matrix and fracture. The coal matrix is also usually simplified to shape of cubes. Therefore, the mechanical properties of coal are considered to be nearly orthotropic [48].Figure 4
Stress loading path and loading direction variations in the true triaxial tests. (a) Stress loading path. (b) Variations of stress loading direction. [49].
(a)(b)According to the geological structure characteristics of coal seam, the raw cubic coal samples were prepared in the orthogonal fracture directions in Liu’s research [49]. The orthogonal fracture directions were bedding plane direction (BD), face cleat plane direction (FD), and butt cleat plane direction (UD), respectively. Then, the true triaxial loading test was carried out with the cubic coal samples in two true triaxial stress paths. The two true triaxial stress paths were conducted as follows (Figure 4(a)). Firstly, the coal sample was loaded to the hydrostatic pressure (σ1=σ2=σ3) of 10 MPa. Next, σ3 was kept constant while the σ1 and σ2 were increased to the design values of σ2 simultaneously. The design values of σ2 were 20 and 30 MPa, respectively. Finally, the σ2 and σ3 were kept constant, but the σ1 was increased until the coal sample failure.In this way, the complete stress-strain of coal samples under the true triaxial stress conditon was tested. The loading directions of three principal stresses (σ1, σ2, σ3) were parallel or perpendicular to the original fractures (bedding plane, face cleat, and butt cleat) in the tests. During the whole tests, the directions of the three principal stresses were kept unchanged, but the directions of original fractures in cubic coal sample were changed in turn (Figure 4(b)). Based on the true triaxial test results under different stress paths, the developed SDAR model is validated. The parameters used for validation of coal are illustrated in Table 2.Table 2
Mechanical parameters of coal in different directions [49, 50].
Stress loading pathElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)Path 1 (σ2=20MPa, σ3=10MPa)4.113.283.240.140.2228Path 2 (σ2=30MPa, σ3=10MPa,)4.553.672.870.140.2228The theoretical results of the SDAR model and test data of coal on stress-strain are shown in Figure5. By comparing the stress-strain in each stress path, the peak strength of coal samples prepared in three directions is in the same order. The peak strength of coal sample along the BD is the largest. And the peak strength of coal sample along UD is the smallest. Moreover, the peak strain also decreased in the order of BD, FD, and UD.Figure 5
Stress-strain of coal in orthogonal directions in different stress loading paths. (a) Path 1 (σ2=20MPa, σ3=10MPa). (b) Path 2 (σ2=30MPa, σ3=10MPa).
(a)(b)Compared with the test data, the results of the SDAR model can basically represent the stress-strain relationship of coal in different directions under true triaxial condition. In the phase of prepeak, the theoretical results are in good agreement with the test data. The peak stress of the SDAR model is very close to the experimental value of coal in each case. In the postpeak phase, all the theoretical curves obviously present variation of strain softening, which is the same as the test data. All the values of the SDAR model decrease rapidly after the peak stress point. The corresponding tested results also show a decrease trend, but the decrease rate is relatively slow. There is a certain deviation between the SDAR model and the test results in the postpeak phase. In the postpeak phase, there are many reasons for this deviation, including the properties of rock, test conditions, loading procedure, and failure mode of rock [20]. On the whole, stress-strain relationship of orthotropic coal under true triaxial loading can be well evaluated by the SDAR model.


## 3.1.1. Transversely Isotropic Shale under Conventional Triaxial Stress Condition

As a kind of sedimentary rock, shale is observed to be transversely isotropic. Shale has the same mechanical properties in the horizontal direction, but the mechanical properties in the horizontal direction are different from those in the vertical direction. Heng [44] used shale which is from the Longmaxi Formation in Chongqing to drill tested samples along two orthognoal directions (Figure 2). One is parallel to bedding (drilling angle 0°). The other is perpendicular to bedding (drilling angle 90°). Shale samples drilled from different directions were tested under conventional triaxial compression with confining pressures of 10, 20, and 30 MPa, respectively. Then, the stress-strain relationships of shale which were drilled in orthogonal directions under different confining pressures that are obtained. The parameters used for validation of shale are shown in Table 1.Figure 2
Directional coring diagram of shale.Table 1
Mechanical parameters of shale in different directions [44].
Confining pressure (MPa)Elastic modulus parallel to bedding (GPa)Elastic modulus perpendicular to bedding (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)1022.9142.940.3840.32530.32022.7752.510.4160.26330.33024.0456.450.4220.25830.3Comparing the theoretical results of the SDAR model and test data of shale on stress-strain in Figure3, it can be seen that the theoretical results are very close to the test data in the same drilling direction under confining pressures of 10, 20, and 30 MPa, respectively. The comparison results show that all the stress-strain curves of shale under conventional triaxial compression can well be evaluated by the SDAR model. In the prepeak zone, the curve of theoretical calculation is in good agreement with that of test data. After the peak stress point, the curve of theoretical results present strain softening in the postpeak region. The presentation of strain softening in the SDAR model is also consistent with the mechanical behavior of shale.Figure 3
Stress-strain curves of shale under different confining pressures: (a) confining pressure 10 MPa, (b) confining pressure 20 MPa, and (c) confining pressure 30 MPa.
(a)(b)(c)Comparing the peak strength of shale in the same drilling direction under each confining pressure (Figure3), the peak strength of samples in the same drilling direction increases with the increase of confining pressure. Moreover, the peak strength under different confining pressures can be accurately calculated by the SDAR model. Then, comparing the results of samples drilled in different directions under the same confining pressure, it can be found that the peak strength of sample drilled in the direction of parallel to bedding is always bigger than that of perpendicular to bedding. Meanwhile, the peak strain of sample drilled in the direction of parallel to bedding is smaller than that of perpendicular to bedding. The stress-strain of transversely isotropic shale under conventional triaxial stress can be well calculated by the SDAR model.


## 3.1.2. Orthogonal Anisotropic Coal under True Triaxial Stress Condition

Under the influence of sedimentation, there is an obvious approximate horizontal bedding plane structure in coal. In addition, there are a lot of cleats in coal, which can be divided into face cleat and butt cleat [45]. These two kinds of cleats are approximately orthogonal to each other and perpendicular to the bedding plane (Figure 4(a)) [46, 47]. Due to the existence of cleat and bedding (referred to as fractures), coal is generally divided into matrix and fracture. The coal matrix is also usually simplified to shape of cubes. Therefore, the mechanical properties of coal are considered to be nearly orthotropic [48].Figure 4
Stress loading path and loading direction variations in the true triaxial tests. (a) Stress loading path. (b) Variations of stress loading direction. [49].
(a)(b)According to the geological structure characteristics of coal seam, the raw cubic coal samples were prepared in the orthogonal fracture directions in Liu’s research [49]. The orthogonal fracture directions were bedding plane direction (BD), face cleat plane direction (FD), and butt cleat plane direction (UD), respectively. Then, the true triaxial loading test was carried out with the cubic coal samples in two true triaxial stress paths. The two true triaxial stress paths were conducted as follows (Figure 4(a)). Firstly, the coal sample was loaded to the hydrostatic pressure (σ1=σ2=σ3) of 10 MPa. Next, σ3 was kept constant while the σ1 and σ2 were increased to the design values of σ2 simultaneously. The design values of σ2 were 20 and 30 MPa, respectively. Finally, the σ2 and σ3 were kept constant, but the σ1 was increased until the coal sample failure.In this way, the complete stress-strain of coal samples under the true triaxial stress conditon was tested. The loading directions of three principal stresses (σ1, σ2, σ3) were parallel or perpendicular to the original fractures (bedding plane, face cleat, and butt cleat) in the tests. During the whole tests, the directions of the three principal stresses were kept unchanged, but the directions of original fractures in cubic coal sample were changed in turn (Figure 4(b)). Based on the true triaxial test results under different stress paths, the developed SDAR model is validated. The parameters used for validation of coal are illustrated in Table 2.Table 2
Mechanical parameters of coal in different directions [49, 50].
Stress loading pathElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Poisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingFriction angle (°)Path 1 (σ2=20MPa, σ3=10MPa)4.113.283.240.140.2228Path 2 (σ2=30MPa, σ3=10MPa,)4.553.672.870.140.2228The theoretical results of the SDAR model and test data of coal on stress-strain are shown in Figure5. By comparing the stress-strain in each stress path, the peak strength of coal samples prepared in three directions is in the same order. The peak strength of coal sample along the BD is the largest. And the peak strength of coal sample along UD is the smallest. Moreover, the peak strain also decreased in the order of BD, FD, and UD.Figure 5
Stress-strain of coal in orthogonal directions in different stress loading paths. (a) Path 1 (σ2=20MPa, σ3=10MPa). (b) Path 2 (σ2=30MPa, σ3=10MPa).
(a)(b)Compared with the test data, the results of the SDAR model can basically represent the stress-strain relationship of coal in different directions under true triaxial condition. In the phase of prepeak, the theoretical results are in good agreement with the test data. The peak stress of the SDAR model is very close to the experimental value of coal in each case. In the postpeak phase, all the theoretical curves obviously present variation of strain softening, which is the same as the test data. All the values of the SDAR model decrease rapidly after the peak stress point. The corresponding tested results also show a decrease trend, but the decrease rate is relatively slow. There is a certain deviation between the SDAR model and the test results in the postpeak phase. In the postpeak phase, there are many reasons for this deviation, including the properties of rock, test conditions, loading procedure, and failure mode of rock [20]. On the whole, stress-strain relationship of orthotropic coal under true triaxial loading can be well evaluated by the SDAR model.


## 3.2. Evolution of Damage Variable

As a heterogeneous geological material, rock has many original defects. Under the influence of various stresses, the defects continuously emerge and expand, which result in damage accumulation in rock. Damage accumulation is the main reason for progressive failure of rock. It is also the important factor determining the stress-strain relationship of rock. It is of great significance to investigate the rock damage evolution under loading for developing the constitutive model.Based on Liu’s research [49], the damage evolution of orthotropic coal samples drilled in different directions under true triaxial loading is analyzed by the SDAR model. Figure 6 demonstrates the evolution of coal samples prepared in each direction (BD, FD, and UD) under different stress paths. It can be seen that when the axial strain is less than a certain value, the value of damage variable in all directions is 0. In other words, there is no damage in the early phase of loading. In the process of test, the axial strain increases continuously with the increase of loading stress. After the peak stress point, the stress-strain curves present the strain softening with the decrease of stress. But during the whole loading process, both the damage variable and the axial strain are monotonically increase. The value of damage variable increases from 0 to1.Figure 6
Damage evolution of anisotropic coal under different stress loading paths. (a) Path 1 (σ2=20MPa, σ3=10MPa). (b) Path 2 (σ2=30MPa, σ3=10MPa).
(a)(b)In stress path 1 (σ2=20MPa, σ3=10MPa), the damage evolution of coal samples loading along the FD and UD is similar. The initial damage first emerges in the coal sample loaded in FD, but the damage variable of coal sample loaded in UD first reaches 1. The initial damage of coal sample loaded in BD emerges the latest. Meanwhile, the damage accumulation rate of coal loaded in BD is the slowest. Therefore, the increase rate of damage variable in BD is the slowest, and the value of axial strain corresponding to damage variable reaching 1 is the largest.In stress path 2 (σ2=30MPa, σ3=10MPa), the initial damage first emerges in the coal sample loaded in UD, but the damage variable of coal sample loaded in FD first reaches 1. The development of damage in the BD is still the slowest. The SDAR model established in this study can well evaluate the damage evolution of anisotropic rock. And the SDAR model also provides a theoretical analysis method for the damage evolution of anisotropic rock under true triaxial condition.


## 4. Sensitivity Study of Parameters

The sensitivity analysis was based on the test in BD under stress path 1 (σ2=20MPa, σ3=10MPa) of Liu’s research [49]. Sensitivity analysis is conducted from two aspects: mechanical parameters and model parameters. The mechanical parameters include anisotropic elastic modulus and anisotropic Poisson’s ratio. The model parameters refer to wi and F0i.


### 4.1. Sensitivity of Mechanical Parameters

#### 4.1.1. Anisotropic Elastic Modulus

In order to research the influence of elastic modulus in each direction on the stress-strain relationship and damage evolution, the sensitivity of elastic modulus in orthogonal directions (BD, FD, and UD) is analyzed on the basis of only changing the value of elastic modulus in each direction. The other parameters remain constant. The values of elastic modulus in all three directions are listed in Table3. Poisson’s ratio of coal in all directions is 0.3, and the friction angle of coal is 28°. The loading path of true triaxial stress in the analysis is based on the path 1 (σ2=20MPa, σ3=10MPa) in Figure 4(a). According to the above parameters and stress path, the axial strain and damage evolution of coal samples prepared along the butt cleat plane direction under true triaxial condition are analyzed.Table 3
Values of elastic modulus in orthogonal directions.
Serial numberElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Case 1333Case 23.533Case 3433Case 433.53Case 5343Case 6333.5Case 7334There are different stress-strain relationships and damage evolution under various combinations of elastic modulus. Figure7 shows the results of stress-strain and damage under different anisotropic elastic moduli.Figure 7
Stress-strain and damage evolution of different elastic moduli in orthogonal directions. (a) Stress-strain curves. (b) Damage variable.
(a)(b)Comparing the cases 1, 2, and 3 in Figure7(a), the results show that the larger the elastic modulus along the face cleat plane direction is, the higher the peak strength of coal is, and the faster the strain softening rate is after peak point. From the comparison among cases 1, 2, and 3 in Figure 7(b), it shows that the larger the elastic modulus along the face cleat plane direction is, the later the accelerated damage accumulation phase emerges. Meanwhile, the increasing rate of damage accumulation becomes faster, and the damage variable reaches 1 earlier.Then, comparing the cases 1, 4, and 5 in Figure7(a), it can be obtained that with the increase of elastic modulus along the butt cleat plane direction, the slope of the elastic deformation phase in the stress-strain curve increases while the peak strength of coal decreases. At the same time, the rate of stress dropping increases after peak point. The results of comparison among cases 1, 4, and 5 in Figure 7(b) demonstrate that the larger the elastic modulus along the butt cleat plane direction is, the earlier the accelerated damage accumulation phase emerges. Moreover, the increasing rate of damage accumulation becomes slower, and the value of damage variable reaches 1 later.By comparing the curves of cases 1, 6, and 7 in Figure7, we can know that the effects of elastic modulus along the bedding plane direction on the variation of peak stress, stress-strain relationship, and damage are same as the effect caused by the elastic modulus along the face cleat plane direction, but the variation degree of curves under the influence of elastic modulus along the bedding plane direction is smaller. In other words, the sensitivity of elastic modulus along the bedding plane direction is not as high as that of elastic modulus along the face cleat plane direction.Because the sensitivity analysis is focused on the case that the maximum principal stress (σ1) is loaded along the butt cleat plane direction, the variation of the elastic modulus along the butt cleat plane direction has a significant influence on the stress-strain and damage, but the change of elastic modulus in the other direction, as the face cleat plane or bedding plane, has little effect on the stress-strain and damage evolution.


#### 4.1.2. Anisotropic Poisson’s Ratio

For the purpose of investigating the influence of Poisson’s ratio on the stress-strain and damage evolution, the values of Poisson’s ratio in different directions are changed while the other parameters are kept constant in the sensitivity analysis. The elastic modulus of coal in all directions is 3 GPa, and the friction angle of coal is 28°. It can be known that the strain in thei-axis direction is only influenced by Poisson’s ratio in the other two directions (vij and vik) from Eqs. (10), (23), and (24). Then, the calculation is carried out according to the case that the maximum principal stress (σ1) is loaded along the direction of butt cleat plane in Figure 4. The values of Poisson’s ratio in all directions are shown in Table 4.Table 4
Values of Poisson’s ratio in different directions.
Serial numberPoisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingCase 10.200.20Case 20.250.20Case 30.350.20Case 40.200.25Case 50.200.35The stress-strain and damage of coal under different Poisson’s ratio are shown in Figure8. The values of the slope concerning elastic deformation phase in the stress-strain curves under different Poisson’s ratios are equal (Figure 8(a)). With the independent increase of Poisson’s ratio in a certain direction (parallel to bedding or perpendicular to bedding), the peak strength, peak strain, and the rate of stress dropping decrease. Meanwhile, with the increase of Poisson’s ratio, both the axial strain corresponding to the initial damage and the increase rate of damage variable decrease (Figure 8(b)).Figure 8
Stress-strain and damage evolution of anisotropic Poisson’s ratio. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


### 4.2. Sensitivity of Model Parameters

#### 4.2.1. Parameter ofwi

In order to study the parameter ofwi, the value of wi is increased at intervals of 3. And the values of wi are shown in Table 5. On the premise that other parameters remain unchanged, with the increase of wi, the phase of linear elastic deformation becomes longer, and the peak strength and peak strain increase obviously (Figure 9(a)). Some changes in the damage variable corresponding to the stress-strain are also observed. With the increase of wi, the axial strain corresponding to initial damage and the rate of damage accumulation become larger (Figure 9(b)). The value of damage variable would reach up to 1 as soon as wi increases.Table 5
Increasing values ofwi.
Serial numberwiCase 16.32Case 29.32Case 312.32Case 415.32Case 518.32Figure 9
Stress-strain and damage evolution of differentwi. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


#### 4.2.2. Parameter ofFi0

In order to analyze the sensitivity ofFi0, the values of Fi0 are changed according to Table 6. The stress-strain and damage variable can be obtained in different values of Fi0. As it can be seen in Figure 10(a), with the increase of Fi0, the phase of elastic deformation becomes longer, and the peak strength and peak strain increase obviously. Moreover, the stress after peak point decreases based on the same slope at different values of Fi0. Meanwhile, the axial strain corresponding to the emergence of initial damage becomes larger with the increase of Fi0 (Figure 10(b)). It indicates that higher loading stress is required for emergence of coal damage. In addition, the larger the Fi0 is, the larger the axial strain corresponding to the damage variable reaching 1 becomes.Table 6
Increasing values ofFi0.
Serial numberFi0 (MPa)Case 153.57Case 263.57Case 373.57Case 483.57Case 593.57Figure 10
Stress strain and damage evolution of differentFi0. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 4.1. Sensitivity of Mechanical Parameters

### 4.1.1. Anisotropic Elastic Modulus

In order to research the influence of elastic modulus in each direction on the stress-strain relationship and damage evolution, the sensitivity of elastic modulus in orthogonal directions (BD, FD, and UD) is analyzed on the basis of only changing the value of elastic modulus in each direction. The other parameters remain constant. The values of elastic modulus in all three directions are listed in Table3. Poisson’s ratio of coal in all directions is 0.3, and the friction angle of coal is 28°. The loading path of true triaxial stress in the analysis is based on the path 1 (σ2=20MPa, σ3=10MPa) in Figure 4(a). According to the above parameters and stress path, the axial strain and damage evolution of coal samples prepared along the butt cleat plane direction under true triaxial condition are analyzed.Table 3
Values of elastic modulus in orthogonal directions.
Serial numberElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Case 1333Case 23.533Case 3433Case 433.53Case 5343Case 6333.5Case 7334There are different stress-strain relationships and damage evolution under various combinations of elastic modulus. Figure7 shows the results of stress-strain and damage under different anisotropic elastic moduli.Figure 7
Stress-strain and damage evolution of different elastic moduli in orthogonal directions. (a) Stress-strain curves. (b) Damage variable.
(a)(b)Comparing the cases 1, 2, and 3 in Figure7(a), the results show that the larger the elastic modulus along the face cleat plane direction is, the higher the peak strength of coal is, and the faster the strain softening rate is after peak point. From the comparison among cases 1, 2, and 3 in Figure 7(b), it shows that the larger the elastic modulus along the face cleat plane direction is, the later the accelerated damage accumulation phase emerges. Meanwhile, the increasing rate of damage accumulation becomes faster, and the damage variable reaches 1 earlier.Then, comparing the cases 1, 4, and 5 in Figure7(a), it can be obtained that with the increase of elastic modulus along the butt cleat plane direction, the slope of the elastic deformation phase in the stress-strain curve increases while the peak strength of coal decreases. At the same time, the rate of stress dropping increases after peak point. The results of comparison among cases 1, 4, and 5 in Figure 7(b) demonstrate that the larger the elastic modulus along the butt cleat plane direction is, the earlier the accelerated damage accumulation phase emerges. Moreover, the increasing rate of damage accumulation becomes slower, and the value of damage variable reaches 1 later.By comparing the curves of cases 1, 6, and 7 in Figure7, we can know that the effects of elastic modulus along the bedding plane direction on the variation of peak stress, stress-strain relationship, and damage are same as the effect caused by the elastic modulus along the face cleat plane direction, but the variation degree of curves under the influence of elastic modulus along the bedding plane direction is smaller. In other words, the sensitivity of elastic modulus along the bedding plane direction is not as high as that of elastic modulus along the face cleat plane direction.Because the sensitivity analysis is focused on the case that the maximum principal stress (σ1) is loaded along the butt cleat plane direction, the variation of the elastic modulus along the butt cleat plane direction has a significant influence on the stress-strain and damage, but the change of elastic modulus in the other direction, as the face cleat plane or bedding plane, has little effect on the stress-strain and damage evolution.


### 4.1.2. Anisotropic Poisson’s Ratio

For the purpose of investigating the influence of Poisson’s ratio on the stress-strain and damage evolution, the values of Poisson’s ratio in different directions are changed while the other parameters are kept constant in the sensitivity analysis. The elastic modulus of coal in all directions is 3 GPa, and the friction angle of coal is 28°. It can be known that the strain in thei-axis direction is only influenced by Poisson’s ratio in the other two directions (vij and vik) from Eqs. (10), (23), and (24). Then, the calculation is carried out according to the case that the maximum principal stress (σ1) is loaded along the direction of butt cleat plane in Figure 4. The values of Poisson’s ratio in all directions are shown in Table 4.Table 4
Values of Poisson’s ratio in different directions.
Serial numberPoisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingCase 10.200.20Case 20.250.20Case 30.350.20Case 40.200.25Case 50.200.35The stress-strain and damage of coal under different Poisson’s ratio are shown in Figure8. The values of the slope concerning elastic deformation phase in the stress-strain curves under different Poisson’s ratios are equal (Figure 8(a)). With the independent increase of Poisson’s ratio in a certain direction (parallel to bedding or perpendicular to bedding), the peak strength, peak strain, and the rate of stress dropping decrease. Meanwhile, with the increase of Poisson’s ratio, both the axial strain corresponding to the initial damage and the increase rate of damage variable decrease (Figure 8(b)).Figure 8
Stress-strain and damage evolution of anisotropic Poisson’s ratio. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 4.1.1. Anisotropic Elastic Modulus

In order to research the influence of elastic modulus in each direction on the stress-strain relationship and damage evolution, the sensitivity of elastic modulus in orthogonal directions (BD, FD, and UD) is analyzed on the basis of only changing the value of elastic modulus in each direction. The other parameters remain constant. The values of elastic modulus in all three directions are listed in Table3. Poisson’s ratio of coal in all directions is 0.3, and the friction angle of coal is 28°. The loading path of true triaxial stress in the analysis is based on the path 1 (σ2=20MPa, σ3=10MPa) in Figure 4(a). According to the above parameters and stress path, the axial strain and damage evolution of coal samples prepared along the butt cleat plane direction under true triaxial condition are analyzed.Table 3
Values of elastic modulus in orthogonal directions.
Serial numberElastic modulus along the face cleat plane direction (GPa)Elastic modulus along the butt cleat plane direction (GPa)Elastic modulus along the bedding plane direction (GPa)Case 1333Case 23.533Case 3433Case 433.53Case 5343Case 6333.5Case 7334There are different stress-strain relationships and damage evolution under various combinations of elastic modulus. Figure7 shows the results of stress-strain and damage under different anisotropic elastic moduli.Figure 7
Stress-strain and damage evolution of different elastic moduli in orthogonal directions. (a) Stress-strain curves. (b) Damage variable.
(a)(b)Comparing the cases 1, 2, and 3 in Figure7(a), the results show that the larger the elastic modulus along the face cleat plane direction is, the higher the peak strength of coal is, and the faster the strain softening rate is after peak point. From the comparison among cases 1, 2, and 3 in Figure 7(b), it shows that the larger the elastic modulus along the face cleat plane direction is, the later the accelerated damage accumulation phase emerges. Meanwhile, the increasing rate of damage accumulation becomes faster, and the damage variable reaches 1 earlier.Then, comparing the cases 1, 4, and 5 in Figure7(a), it can be obtained that with the increase of elastic modulus along the butt cleat plane direction, the slope of the elastic deformation phase in the stress-strain curve increases while the peak strength of coal decreases. At the same time, the rate of stress dropping increases after peak point. The results of comparison among cases 1, 4, and 5 in Figure 7(b) demonstrate that the larger the elastic modulus along the butt cleat plane direction is, the earlier the accelerated damage accumulation phase emerges. Moreover, the increasing rate of damage accumulation becomes slower, and the value of damage variable reaches 1 later.By comparing the curves of cases 1, 6, and 7 in Figure7, we can know that the effects of elastic modulus along the bedding plane direction on the variation of peak stress, stress-strain relationship, and damage are same as the effect caused by the elastic modulus along the face cleat plane direction, but the variation degree of curves under the influence of elastic modulus along the bedding plane direction is smaller. In other words, the sensitivity of elastic modulus along the bedding plane direction is not as high as that of elastic modulus along the face cleat plane direction.Because the sensitivity analysis is focused on the case that the maximum principal stress (σ1) is loaded along the butt cleat plane direction, the variation of the elastic modulus along the butt cleat plane direction has a significant influence on the stress-strain and damage, but the change of elastic modulus in the other direction, as the face cleat plane or bedding plane, has little effect on the stress-strain and damage evolution.


## 4.1.2. Anisotropic Poisson’s Ratio

For the purpose of investigating the influence of Poisson’s ratio on the stress-strain and damage evolution, the values of Poisson’s ratio in different directions are changed while the other parameters are kept constant in the sensitivity analysis. The elastic modulus of coal in all directions is 3 GPa, and the friction angle of coal is 28°. It can be known that the strain in thei-axis direction is only influenced by Poisson’s ratio in the other two directions (vij and vik) from Eqs. (10), (23), and (24). Then, the calculation is carried out according to the case that the maximum principal stress (σ1) is loaded along the direction of butt cleat plane in Figure 4. The values of Poisson’s ratio in all directions are shown in Table 4.Table 4
Values of Poisson’s ratio in different directions.
Serial numberPoisson’s ratio parallel to beddingPoisson’s ratio perpendicular to beddingCase 10.200.20Case 20.250.20Case 30.350.20Case 40.200.25Case 50.200.35The stress-strain and damage of coal under different Poisson’s ratio are shown in Figure8. The values of the slope concerning elastic deformation phase in the stress-strain curves under different Poisson’s ratios are equal (Figure 8(a)). With the independent increase of Poisson’s ratio in a certain direction (parallel to bedding or perpendicular to bedding), the peak strength, peak strain, and the rate of stress dropping decrease. Meanwhile, with the increase of Poisson’s ratio, both the axial strain corresponding to the initial damage and the increase rate of damage variable decrease (Figure 8(b)).Figure 8
Stress-strain and damage evolution of anisotropic Poisson’s ratio. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 4.2. Sensitivity of Model Parameters

### 4.2.1. Parameter ofwi

In order to study the parameter ofwi, the value of wi is increased at intervals of 3. And the values of wi are shown in Table 5. On the premise that other parameters remain unchanged, with the increase of wi, the phase of linear elastic deformation becomes longer, and the peak strength and peak strain increase obviously (Figure 9(a)). Some changes in the damage variable corresponding to the stress-strain are also observed. With the increase of wi, the axial strain corresponding to initial damage and the rate of damage accumulation become larger (Figure 9(b)). The value of damage variable would reach up to 1 as soon as wi increases.Table 5
Increasing values ofwi.
Serial numberwiCase 16.32Case 29.32Case 312.32Case 415.32Case 518.32Figure 9
Stress-strain and damage evolution of differentwi. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


### 4.2.2. Parameter ofFi0

In order to analyze the sensitivity ofFi0, the values of Fi0 are changed according to Table 6. The stress-strain and damage variable can be obtained in different values of Fi0. As it can be seen in Figure 10(a), with the increase of Fi0, the phase of elastic deformation becomes longer, and the peak strength and peak strain increase obviously. Moreover, the stress after peak point decreases based on the same slope at different values of Fi0. Meanwhile, the axial strain corresponding to the emergence of initial damage becomes larger with the increase of Fi0 (Figure 10(b)). It indicates that higher loading stress is required for emergence of coal damage. In addition, the larger the Fi0 is, the larger the axial strain corresponding to the damage variable reaching 1 becomes.Table 6
Increasing values ofFi0.
Serial numberFi0 (MPa)Case 153.57Case 263.57Case 373.57Case 483.57Case 593.57Figure 10
Stress strain and damage evolution of differentFi0. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 4.2.1. Parameter ofwi

In order to study the parameter ofwi, the value of wi is increased at intervals of 3. And the values of wi are shown in Table 5. On the premise that other parameters remain unchanged, with the increase of wi, the phase of linear elastic deformation becomes longer, and the peak strength and peak strain increase obviously (Figure 9(a)). Some changes in the damage variable corresponding to the stress-strain are also observed. With the increase of wi, the axial strain corresponding to initial damage and the rate of damage accumulation become larger (Figure 9(b)). The value of damage variable would reach up to 1 as soon as wi increases.Table 5
Increasing values ofwi.
Serial numberwiCase 16.32Case 29.32Case 312.32Case 415.32Case 518.32Figure 9
Stress-strain and damage evolution of differentwi. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 4.2.2. Parameter ofFi0

In order to analyze the sensitivity ofFi0, the values of Fi0 are changed according to Table 6. The stress-strain and damage variable can be obtained in different values of Fi0. As it can be seen in Figure 10(a), with the increase of Fi0, the phase of elastic deformation becomes longer, and the peak strength and peak strain increase obviously. Moreover, the stress after peak point decreases based on the same slope at different values of Fi0. Meanwhile, the axial strain corresponding to the emergence of initial damage becomes larger with the increase of Fi0 (Figure 10(b)). It indicates that higher loading stress is required for emergence of coal damage. In addition, the larger the Fi0 is, the larger the axial strain corresponding to the damage variable reaching 1 becomes.Table 6
Increasing values ofFi0.
Serial numberFi0 (MPa)Case 153.57Case 263.57Case 373.57Case 483.57Case 593.57Figure 10
Stress strain and damage evolution of differentFi0. (a) Stress-strain curves. (b) Damage variable.
(a)(b)


## 5. Conclusions

(1)
Based on the proposed concept of areal damage and statistical damage mechanics, a statistical damage constitutive model of anisotropic rock (SDAR model) is established. According to the test validation, it is proved that the model can well evaluate the stress-strain relationship and damage evolution of anisotropic rock in different directions(2)
The three principal stresses in orthogonal directions are independent of each other in the process of developing the SDAR model. Therefore, the SDAR model can be applied to predict the stress-strain and damage evolution of anisotropic rock under true triaxial loading condition(3)
The effect of elastic modulus in one direction on the strain of the same direction is more significant than that in other directions. With the increase of elastic modulus in one direction, the slope of elastic deformation in the stress-strain curve of the same direction becomes larger. Meanwhile, the accelerated damage accumulation phase begins. With the increase of Poisson’s ratio, the peak strength, peak strain, and the rate of stress dropping in the postpeak zone decrease(4)
The elastic phase in the stress-strain curve becomes longer with the increased value ofwi. Meanwhile, the peak strength, peak strain, and rate of stress dropping after peak point increase with the increased value of wi. Moreover, the rate of damage accumulation increases as well. On the basis of increasing Fi0, the elastic deformation phase in the stress-strain curve becomes longer. Both peak strength and peak strain tend to be larger. The rate of stress dropping after peak point is constant. And the axial strain corresponding to the emergence of initial damage becomes larger



---

*Source: 6307895-2021-08-13.xml*