# A Statistical Constitutive Model considering Deterioration for Brittle Rocks under a Coupled Thermal-Mechanical Condition

Authors: Meiben Gao; Tianbin Li; Tao Wei; Lubo Meng
Journal: Geofluids
 (2018)

---

## Abstract

Due to active actions of groundwater and geothermal, the stability of underground engineering is important during geological structure active area. The damage mechanical theory and statistical mesoscopic strength theory based on Weibull distribution are widely used to discuss constitutive behaviors of rocks. In these theories, a statistical method is used to capture mesoscopic properties of rocks in order to generate a realistic behavior at a macroscopic scale. Based on the above theories, this paper aims at establishing a constitutive relation of brittle rocks under thermal-mechanical coupling conditions. First, a statistical damage constitutive model was established by considering the thermal effects and crack initiation strength. Subsequently, the parameters of the model were determined and expressed according to the characteristics of stress-strain curve. Third, the model was verified by conventional triaxial experiments under thermal-mechanical actions, and the experimental data and theoretical results were compared and analyzed in the case study. Finally, the physical meaning of the parameters and their effects on the model performance were discussed.


---

## Body

## 1. Introduction

With this longer, larger, and deeper underground engineering’s construction, the rock deterioration caused by groundwater and geothermal has become increasingly prominent. Constitutive modeling of rocks, as a key in the theory development and numerical analysis of rock mechanics, has become a topic of everlasting interest to researchers and practitioners working on geosciences and geotechnics. In addition, many factors affect rock behaviors, such as water-rock interactions and temperature, which would lead to its deterioration. With the deterioration of rocks and soil, the state of the geoenvironment would be changed resulting in landslide [1, 2], rock avalanche [3, 4], and rock burst [5]. For brittle rocks, the deterioration process in laboratory samples is a consequence of microscale and macroscale fracturing that occurs in several stages. Recent research suggests that crack initiation stress can be used as an estimate for in situ spalling strength, which is commonly observed in brittle rocks around underground openings [6–10]. Additionally, under high-temperature and high-pressure conditions, the mechanical characteristics of deep rocks exhibit different behaviors compared to those at more mild temperatures at lower depths. A reasonable constitutive model considering thermal action and crack initiation is the key to accurate prediction and judgment of the reliability and stability of those works such as the exploitation of deep mining resources and geothermal resources. Therefore, the constitutive behaviors of rocks considering thermal effects and crack initiation strength are significant to be studied and explored in order to better solve rock engineering problems in thermal-mechanical coupling conditions.During the preceding four decades, various rock constitutive models have been established from theoretical, experimental approaches (e.g., [11–15]). Based on the traditional continuum mechanics and damage theory, thermoplastic and thermoelastic brittle models were proposed [16, 17]. However, these models failed to reflect the damage features of brittle rocks. For rock materials, numerous microscopic cracks, ranging from 0.01 to 1.0 mm in length, are statistically distributed, which have significant influences on the damage processes and failure characteristics of rocks. The nucleation and growth of microcracks would lead to a concentration of these microcracks into a narrow zone and produce a visible macroscopic fissure wider than 1.0 mm [18], so the process from damage to fracture could be studied on a mesoscopic scale. Since the statistical damage-based approach has been used successfully to address rock constitutive behaviors [19–21], as a quite attractive tool for investigating deformation processes and failure mechanisms in the mechanics of geomaterial community, it has been especially favored by many researchers. However, these previous statistical models did not reflect the residual strength of rocks. By introducing a coefficient into the damage variable, the statistical damage-based model will be able to describe the residual strength [22, 23]. For thermal effects, Zhang et al. [24] proposed a three-parameter Weibull distribution to express the rock constitutive behavior under a uniaxial compression test, which considers the thermal-mechanical coupling conditions.In the present study, a mesoscopic element is considered to be isotropically elastic and its properties are defined by Young’s modulus or Poisson’s ratio. The stress-strain relationship is linearly elastic until given damage threshold (crack initiation strength) is attained [17]. Thus, the macroscopic strength and the properties of rock depend on the statistical mechanical properties of individual mesoscopic elements, which could be described by a phenomenological model through a statistical method. A previous research work showed that continuum damage models can effectively simulate the elastic degradation caused by preexisting microcracks in rocks. Although rocks always exhibit anisotropy after macroscopic fissures occur, an isotropic damage model is still an effective method to estimate the gross damage of rocks subjected to external loading.This paper aims at exploring a statistical damage constitutive model considering crack initiation strength, deterioration, and thermal effects. Based on a statistical damage approach, the statistical constitutive model is established by introducing a three-parameter Weibull distribution, which describes a brittle constitutive behavior of rocks under thermal-mechanical coupling conditions.


## 2. Constitutive Model

### 2.1. Thermal-Mechanical-Damage Evolution Equations

#### 2.1.1. Thermal Damage

Under the thermal-mechanical coupling effect, numerous microcracks and considerable damage occur in the rock, which changes the mechanical properties. Macroscopically, the elastic modulus is usually chosen as the damage variable. According to Liu et al. [17], the thermal damage DT can be given by
(1)DT=1−ETE0,where ET and E0 are the elastic moduli at T°C and room temperature, respectively.


#### 2.1.2. Damage Evolution Equation considering the Thermal Effect

The randomly distributed preexisting microcracks in rocks are the main factors leading to the damage of the rock.Suppose that a microunit within the rock is sufficiently large to contain numerous cracks and defects and adequately small in dimension compared with the whole rock structure. Obviously, the strengths of these microunits vary randomly according to a Weibull distribution. The rock damage level can be expressed by the ratio of the number of damaged microunits to total units. On the microscopic level, the rock damage evolution process includes three stages, that is, crack initiation, propagation, and coalescence. Based on this, the coupled thermal-mechanical constitutive model of the complete failure process of the rock is proposed from a damage statistical point.Next, before the damage evolution equation presented, some assumptions are given in the following [22]:
(a)
Rock is an isotropic homogenous geomaterial on a macroscale.
(b)
Microunits conform to Hooke’s law, prior to failure.
(c)
Rock damage occurs continuously and is the gradual accumulation of failures on the mesoscopic level.
(d)
Heat diffusion in a rock is only in the form of heat transfer, without considering convection and radiation.
(e)
Microunit strength follows a Weibull distribution, where the three-parameter density and distribution functions are given by the following:(2)Densityfunction:fk=mFk−γFm−1exp−k−γFm,k∈γ,∞,0,k∈−∞,γ,Distributionfunction:Fk=1−exp−k−γFm,k∈γ,∞,where k is the microunit strength and m, F, and γ are the mean uniformity, peak strength, and damage evolution threshold, respectively, which represent the shape, scale, dimension, and position (as shown in Figure 1).Figure 1
Density function of Weibull distributions for each parameter.
(a)
Density function of Weibull distributions for differentm
(b)
Density function of Weibull distributions for differentF
(c)
Density function of Weibull distributions for differentγConsidering the thermal effects, the Weibull parameters can be given by(3)mT=m01−DT,FT=F01−DT,γT=γ01−DT,where m0, F0, and γ0 are the Weibull parameters at room temperature and m, F, and γ are as defined previously.Rock damage is caused by the continuous failure of the microunits within the rock. Assuming that there areNf failed microunits at a single load level, the statistic damage is defined as the ratio of the failure number to the total number (N) of microunits, namely,
(4)D=NfN.Then, at stressfσ, the failure number Nf is
(5)Nf=∫0fσNfσdσ.Substituting (2) and (4) into (3) yields
(6)D=NfN=∫0fσNfσdσN=∫0fσfσdσ=1−exp−fσ−γ01−DTF01−DTm01−DT,Using the Drucker-Prager criterion as the unit’s yield criterion,(7)fσ=k=α0I1+J2,where α0=sinφ/9+3sin2φ, I1 is the first stress invariable, and J2 is the second invariant of deviator stress. These are expressed as follows:
(8)I1=σx∗+σy∗+σz∗=σ1∗+σ2∗+σ3∗,J2=16σ1∗−σ2∗2+σ2∗−σ3∗2+σ3∗−σ1∗2.Therefore, (5) can be written as follows:
(9)D=1−exp−α0I1+J2−γ01−DTF01−DTm01−DT.Rock damage evolution does not occur until the damage threshold is reached, at which time the rock damage value is considered to be zero. When the stress surpasses the threshold, the rock damage value is given by (8), and for the whole stress state, the hard rock damage evolution equation considering the thermal effect can be expressed as follows:
(10)Dσ,T=0,σ<σci,1−exp−α0I1+J2−γ01−DTF01−DTm01−DT,σ≥σci,where Dσ,T is the damage value considering the thermal effect and σci is the crack initiation strength.


### 2.2. Rock Damage Deterioration Statistical Constitutive Model

According to the concept of the theories of Lemaitre [25] strain equivalence and effective stress principle, the damage variable D is the ratio of the number of the failure units to total units; then, the relationship between nominal stress σ and effective stress σ∗ in three-dimensional isotropy damage can be expressed by
(11)σ∗=σ1−D=Cε1−D,where σ is the nominal stress matrix, σ∗ is the effective stress matrix, ε is the strain matrix, C is the elastic flexibility matrix, and D is the damage matrix.During the compression process, rock retains the ability to pass a certain level of compressive stress and shear stress after the peak strength owing to a change in friction and confining pressure. Residual strength can be observed over the whole stress-strain curve. In most cases, the experimental curve of the residual strength is approximately horizontal after the peak strength. However, (10) does not take into account the rock residual strength; therefore, it has been revised by introducing the damage variable correction coefficient δ, and its value is between 0 and 1.By introducing the damage variable correction coefficient, the effective stress can be expressed by(12)σ∗=σ1−δD.Then,(13)σi∗=σi1−δD,i=1,2,3.According to the generalized Hooke’s law, we have(14)εi=1Eσi∗−μσj∗+σk∗,i,j,k=1,2,3.Substituting (12) into (15), we obtain
(15)εi=1E1−δDσi−μσj+σk,i,j,k=1,2,3.SubstitutingD=1−exp−α0I1+J2−γ01−DT/F01−DTm01−DT into (14) leads to
(16)σi=Eεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk.According to (1), the elastic modulus of hard rock under different temperatures is given as follows:
(17)DT=1−ETE0.Substituting (16) into (15), the hard rock damage constitutive equation considering thermal effects is obtained (after the damage threshold point):
(18)σi=E01−DTεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk.The hard rock constitutive equation considering the thermal-mechanical coupling is given by (17) after the damage threshold point. The constitutive equation can be fitted by a polynomial before the damage threshold point. Because the rock stress-strain curve covers the coordinate origin, its function is given as follows:
(19)σi=Aεi2+Bεi.Integrating (17) and (18), the hard rock constitutive equation with the thermal-mechanical coupling is given as follows:
(20)σi=Aεi2+Bεi,0≤ε≤εci,E01−DTεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk,ε>εci,where εci is the strain value at the damage threshold point.


### 2.3. Model Parameters

The determination of the model parameters is closely related to the rock macroparameters. Thus, the determined model parameters have certain physical meaning and a wider applicability. Therefore, the strain values at the damage threshold pointσci,εci and the peak strength σC,εC are introduced as the macroparameters to determine the model parameters.To simplify the calculation, (19) is rewritten in the following form:
(21)σi=Aεi2+Bεi,0≤ε≤εci,Eεi1−δ+δexp−α0I1+J2−γFm+μσj+σk,ε>εci.


#### 2.3.1. Determination ofm, F, and γ

In the rock triaxial tests, the nominal stressσ1,σ2,σ3σ2=σ3 and strain ε1 can be measured. Based on Hooke’s Law and the effective stress theory, we have
(22)ε1=1Eσ1∗−2μσ3∗,σ2∗=σ3∗=σ21−δD,σ1∗=σ11−δD.Solving (21), σ1∗,σ2∗,σ3∗ can be obtained. Then, I1,J2,J2 is solved:
(23)I1=σ1+2σ3Eε1σ1−2μσ3,J2=σ1−σ3Eε123σ1−2μσ32,J2=σ1−σ3Eε13σ1−2μσ3.In the conventional triaxial test (σ2=σ3) and when ε>εD, (20) can be abbreviated as follows:
(24)σ1=Eε11−δ+δexp−α0I1+J2−γFm+2μσ3.Similarly, parameterγ determines the starting point of the rock damage evolution, that is, the damage evolution threshold point σci,εci. According to fσ−γ≥0 and assuming the damage value is zero here, we obtain
(25)γ=limσ→σciα0I1+J2=α0σci+2σ3Eεciσci−2μσ3+σci−σ3Eεci3σci−2μσ3.By referring the rock stress-strain curve geometry (as shown in Figure2), we know that ①σ=σc,ε=εc;②σ=σc,dσc/dεc=0.Figure 2
Typical stress-strain curve in conventional triaxial compressive test (E is the elastic modulus of rock, σ is the axial stress, and ε is the axial strain).Combining boundary①σ=σc,ε=εc, we obtain
(26)σc=Eεc1−δ+δexp−α0I1+J2−γFm+2μσ3.Substituting boundary②σ=σc,dσc/dεc=0 into the partial derivative of (23), we have
(27)∂σ1∂ε1ε=εc=0.Integrating (26) and (27) yields
(28)α0I1+J2−γFm=InδEεcσc−2μσ3+δ−1Eεc,(29)mα0I1+J2−γm−1Fm=σc−2μσ3εcYσc−2μσ3+δ−1Eεc.Solving (28) and (29), the parameters m and F can be obtained:
(30)m=Xσc−2μσ3εcYσc−2μσ3+δ−1EεcInδEεc/σc−2μσ3+δ−1Eεc,F=XmInδEεc/σc−2μσ3+δ−1Eεc,where X=α0σc+2σ3Eεc/σc−2μσ3+σc−σ3Eεc/3σc−2μσ3−γ and Y=α0σc+2σ3E/σc−2μσ3/+σc−σ3E/3σc−2μσ3.


#### 2.3.2. Determination ofA and B

As the stress is continuous in the compression and damage evolution section of the rock stress-strain curve, we obtain(31)Eεci1−δ+δexp−α0I1+J2−γFm+2μσ3=Aεci2+Bεci.In addition, the first derivative is also continuous at these two sections; then,(32)∂σ1∂ε1ε=εci=2Aεci+B.Dividing (31) by εci on both sides and integrating,
(33)E1−δ+δexp−α0I1+J2−γFm+2μσ3εci=Aεci+B.By supposingR∣ε=εci=E1−δ+δexp−α0I1+J2−γ/Fm+2μσ3/εci, Z=∂σ1/∂ε1∣ε=εci.Equations (32) and (33) can be rewritten as follows:
(34)Z=2Aεci+B,R=Aεci+B.Solving (34), we have
(35)A=Z−Rεci,B=2R−Z.


#### 2.3.3. Determination of Damage Threshold Point (Crack Initiation Point)

In the compression, the original cracks in the rock sample close first under loading. This refers to the key issue in the crack evolution section of the crack initiation and propagation, that is, the damage threshold point. Obviously, it is not easy to obtain the crack initiation stress on a rock stress-strain curve. By researching Lac du Bonnet granite, Martin and Chandler [26] proposed a method to determine the granite initiation point (shown in Figure 3) assuming a crack initiation strength of approximately 40% of the peak strength (σci=40%σc).Figure 3
Crack volumetric strain method (Martin and Chandler [26]).


## 2.1. Thermal-Mechanical-Damage Evolution Equations

### 2.1.1. Thermal Damage

Under the thermal-mechanical coupling effect, numerous microcracks and considerable damage occur in the rock, which changes the mechanical properties. Macroscopically, the elastic modulus is usually chosen as the damage variable. According to Liu et al. [17], the thermal damage DT can be given by
(1)DT=1−ETE0,where ET and E0 are the elastic moduli at T°C and room temperature, respectively.


### 2.1.2. Damage Evolution Equation considering the Thermal Effect

The randomly distributed preexisting microcracks in rocks are the main factors leading to the damage of the rock.Suppose that a microunit within the rock is sufficiently large to contain numerous cracks and defects and adequately small in dimension compared with the whole rock structure. Obviously, the strengths of these microunits vary randomly according to a Weibull distribution. The rock damage level can be expressed by the ratio of the number of damaged microunits to total units. On the microscopic level, the rock damage evolution process includes three stages, that is, crack initiation, propagation, and coalescence. Based on this, the coupled thermal-mechanical constitutive model of the complete failure process of the rock is proposed from a damage statistical point.Next, before the damage evolution equation presented, some assumptions are given in the following [22]:
(a)
Rock is an isotropic homogenous geomaterial on a macroscale.
(b)
Microunits conform to Hooke’s law, prior to failure.
(c)
Rock damage occurs continuously and is the gradual accumulation of failures on the mesoscopic level.
(d)
Heat diffusion in a rock is only in the form of heat transfer, without considering convection and radiation.
(e)
Microunit strength follows a Weibull distribution, where the three-parameter density and distribution functions are given by the following:(2)Densityfunction:fk=mFk−γFm−1exp−k−γFm,k∈γ,∞,0,k∈−∞,γ,Distributionfunction:Fk=1−exp−k−γFm,k∈γ,∞,where k is the microunit strength and m, F, and γ are the mean uniformity, peak strength, and damage evolution threshold, respectively, which represent the shape, scale, dimension, and position (as shown in Figure 1).Figure 1
Density function of Weibull distributions for each parameter.
(a)
Density function of Weibull distributions for differentm
(b)
Density function of Weibull distributions for differentF
(c)
Density function of Weibull distributions for differentγConsidering the thermal effects, the Weibull parameters can be given by(3)mT=m01−DT,FT=F01−DT,γT=γ01−DT,where m0, F0, and γ0 are the Weibull parameters at room temperature and m, F, and γ are as defined previously.Rock damage is caused by the continuous failure of the microunits within the rock. Assuming that there areNf failed microunits at a single load level, the statistic damage is defined as the ratio of the failure number to the total number (N) of microunits, namely,
(4)D=NfN.Then, at stressfσ, the failure number Nf is
(5)Nf=∫0fσNfσdσ.Substituting (2) and (4) into (3) yields
(6)D=NfN=∫0fσNfσdσN=∫0fσfσdσ=1−exp−fσ−γ01−DTF01−DTm01−DT,Using the Drucker-Prager criterion as the unit’s yield criterion,(7)fσ=k=α0I1+J2,where α0=sinφ/9+3sin2φ, I1 is the first stress invariable, and J2 is the second invariant of deviator stress. These are expressed as follows:
(8)I1=σx∗+σy∗+σz∗=σ1∗+σ2∗+σ3∗,J2=16σ1∗−σ2∗2+σ2∗−σ3∗2+σ3∗−σ1∗2.Therefore, (5) can be written as follows:
(9)D=1−exp−α0I1+J2−γ01−DTF01−DTm01−DT.Rock damage evolution does not occur until the damage threshold is reached, at which time the rock damage value is considered to be zero. When the stress surpasses the threshold, the rock damage value is given by (8), and for the whole stress state, the hard rock damage evolution equation considering the thermal effect can be expressed as follows:
(10)Dσ,T=0,σ<σci,1−exp−α0I1+J2−γ01−DTF01−DTm01−DT,σ≥σci,where Dσ,T is the damage value considering the thermal effect and σci is the crack initiation strength.


## 2.1.1. Thermal Damage

Under the thermal-mechanical coupling effect, numerous microcracks and considerable damage occur in the rock, which changes the mechanical properties. Macroscopically, the elastic modulus is usually chosen as the damage variable. According to Liu et al. [17], the thermal damage DT can be given by
(1)DT=1−ETE0,where ET and E0 are the elastic moduli at T°C and room temperature, respectively.


## 2.1.2. Damage Evolution Equation considering the Thermal Effect

The randomly distributed preexisting microcracks in rocks are the main factors leading to the damage of the rock.Suppose that a microunit within the rock is sufficiently large to contain numerous cracks and defects and adequately small in dimension compared with the whole rock structure. Obviously, the strengths of these microunits vary randomly according to a Weibull distribution. The rock damage level can be expressed by the ratio of the number of damaged microunits to total units. On the microscopic level, the rock damage evolution process includes three stages, that is, crack initiation, propagation, and coalescence. Based on this, the coupled thermal-mechanical constitutive model of the complete failure process of the rock is proposed from a damage statistical point.Next, before the damage evolution equation presented, some assumptions are given in the following [22]:
(a)
Rock is an isotropic homogenous geomaterial on a macroscale.
(b)
Microunits conform to Hooke’s law, prior to failure.
(c)
Rock damage occurs continuously and is the gradual accumulation of failures on the mesoscopic level.
(d)
Heat diffusion in a rock is only in the form of heat transfer, without considering convection and radiation.
(e)
Microunit strength follows a Weibull distribution, where the three-parameter density and distribution functions are given by the following:(2)Densityfunction:fk=mFk−γFm−1exp−k−γFm,k∈γ,∞,0,k∈−∞,γ,Distributionfunction:Fk=1−exp−k−γFm,k∈γ,∞,where k is the microunit strength and m, F, and γ are the mean uniformity, peak strength, and damage evolution threshold, respectively, which represent the shape, scale, dimension, and position (as shown in Figure 1).Figure 1
Density function of Weibull distributions for each parameter.
(a)
Density function of Weibull distributions for differentm
(b)
Density function of Weibull distributions for differentF
(c)
Density function of Weibull distributions for differentγConsidering the thermal effects, the Weibull parameters can be given by(3)mT=m01−DT,FT=F01−DT,γT=γ01−DT,where m0, F0, and γ0 are the Weibull parameters at room temperature and m, F, and γ are as defined previously.Rock damage is caused by the continuous failure of the microunits within the rock. Assuming that there areNf failed microunits at a single load level, the statistic damage is defined as the ratio of the failure number to the total number (N) of microunits, namely,
(4)D=NfN.Then, at stressfσ, the failure number Nf is
(5)Nf=∫0fσNfσdσ.Substituting (2) and (4) into (3) yields
(6)D=NfN=∫0fσNfσdσN=∫0fσfσdσ=1−exp−fσ−γ01−DTF01−DTm01−DT,Using the Drucker-Prager criterion as the unit’s yield criterion,(7)fσ=k=α0I1+J2,where α0=sinφ/9+3sin2φ, I1 is the first stress invariable, and J2 is the second invariant of deviator stress. These are expressed as follows:
(8)I1=σx∗+σy∗+σz∗=σ1∗+σ2∗+σ3∗,J2=16σ1∗−σ2∗2+σ2∗−σ3∗2+σ3∗−σ1∗2.Therefore, (5) can be written as follows:
(9)D=1−exp−α0I1+J2−γ01−DTF01−DTm01−DT.Rock damage evolution does not occur until the damage threshold is reached, at which time the rock damage value is considered to be zero. When the stress surpasses the threshold, the rock damage value is given by (8), and for the whole stress state, the hard rock damage evolution equation considering the thermal effect can be expressed as follows:
(10)Dσ,T=0,σ<σci,1−exp−α0I1+J2−γ01−DTF01−DTm01−DT,σ≥σci,where Dσ,T is the damage value considering the thermal effect and σci is the crack initiation strength.


## 2.2. Rock Damage Deterioration Statistical Constitutive Model

According to the concept of the theories of Lemaitre [25] strain equivalence and effective stress principle, the damage variable D is the ratio of the number of the failure units to total units; then, the relationship between nominal stress σ and effective stress σ∗ in three-dimensional isotropy damage can be expressed by
(11)σ∗=σ1−D=Cε1−D,where σ is the nominal stress matrix, σ∗ is the effective stress matrix, ε is the strain matrix, C is the elastic flexibility matrix, and D is the damage matrix.During the compression process, rock retains the ability to pass a certain level of compressive stress and shear stress after the peak strength owing to a change in friction and confining pressure. Residual strength can be observed over the whole stress-strain curve. In most cases, the experimental curve of the residual strength is approximately horizontal after the peak strength. However, (10) does not take into account the rock residual strength; therefore, it has been revised by introducing the damage variable correction coefficient δ, and its value is between 0 and 1.By introducing the damage variable correction coefficient, the effective stress can be expressed by(12)σ∗=σ1−δD.Then,(13)σi∗=σi1−δD,i=1,2,3.According to the generalized Hooke’s law, we have(14)εi=1Eσi∗−μσj∗+σk∗,i,j,k=1,2,3.Substituting (12) into (15), we obtain
(15)εi=1E1−δDσi−μσj+σk,i,j,k=1,2,3.SubstitutingD=1−exp−α0I1+J2−γ01−DT/F01−DTm01−DT into (14) leads to
(16)σi=Eεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk.According to (1), the elastic modulus of hard rock under different temperatures is given as follows:
(17)DT=1−ETE0.Substituting (16) into (15), the hard rock damage constitutive equation considering thermal effects is obtained (after the damage threshold point):
(18)σi=E01−DTεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk.The hard rock constitutive equation considering the thermal-mechanical coupling is given by (17) after the damage threshold point. The constitutive equation can be fitted by a polynomial before the damage threshold point. Because the rock stress-strain curve covers the coordinate origin, its function is given as follows:
(19)σi=Aεi2+Bεi.Integrating (17) and (18), the hard rock constitutive equation with the thermal-mechanical coupling is given as follows:
(20)σi=Aεi2+Bεi,0≤ε≤εci,E01−DTεi1−δ+δexp−α0I1+J2−γ01−DTF01−DTm01−DT+μσj+σk,ε>εci,where εci is the strain value at the damage threshold point.


## 2.3. Model Parameters

The determination of the model parameters is closely related to the rock macroparameters. Thus, the determined model parameters have certain physical meaning and a wider applicability. Therefore, the strain values at the damage threshold pointσci,εci and the peak strength σC,εC are introduced as the macroparameters to determine the model parameters.To simplify the calculation, (19) is rewritten in the following form:
(21)σi=Aεi2+Bεi,0≤ε≤εci,Eεi1−δ+δexp−α0I1+J2−γFm+μσj+σk,ε>εci.


### 2.3.1. Determination ofm, F, and γ

In the rock triaxial tests, the nominal stressσ1,σ2,σ3σ2=σ3 and strain ε1 can be measured. Based on Hooke’s Law and the effective stress theory, we have
(22)ε1=1Eσ1∗−2μσ3∗,σ2∗=σ3∗=σ21−δD,σ1∗=σ11−δD.Solving (21), σ1∗,σ2∗,σ3∗ can be obtained. Then, I1,J2,J2 is solved:
(23)I1=σ1+2σ3Eε1σ1−2μσ3,J2=σ1−σ3Eε123σ1−2μσ32,J2=σ1−σ3Eε13σ1−2μσ3.In the conventional triaxial test (σ2=σ3) and when ε>εD, (20) can be abbreviated as follows:
(24)σ1=Eε11−δ+δexp−α0I1+J2−γFm+2μσ3.Similarly, parameterγ determines the starting point of the rock damage evolution, that is, the damage evolution threshold point σci,εci. According to fσ−γ≥0 and assuming the damage value is zero here, we obtain
(25)γ=limσ→σciα0I1+J2=α0σci+2σ3Eεciσci−2μσ3+σci−σ3Eεci3σci−2μσ3.By referring the rock stress-strain curve geometry (as shown in Figure2), we know that ①σ=σc,ε=εc;②σ=σc,dσc/dεc=0.Figure 2
Typical stress-strain curve in conventional triaxial compressive test (E is the elastic modulus of rock, σ is the axial stress, and ε is the axial strain).Combining boundary①σ=σc,ε=εc, we obtain
(26)σc=Eεc1−δ+δexp−α0I1+J2−γFm+2μσ3.Substituting boundary②σ=σc,dσc/dεc=0 into the partial derivative of (23), we have
(27)∂σ1∂ε1ε=εc=0.Integrating (26) and (27) yields
(28)α0I1+J2−γFm=InδEεcσc−2μσ3+δ−1Eεc,(29)mα0I1+J2−γm−1Fm=σc−2μσ3εcYσc−2μσ3+δ−1Eεc.Solving (28) and (29), the parameters m and F can be obtained:
(30)m=Xσc−2μσ3εcYσc−2μσ3+δ−1EεcInδEεc/σc−2μσ3+δ−1Eεc,F=XmInδEεc/σc−2μσ3+δ−1Eεc,where X=α0σc+2σ3Eεc/σc−2μσ3+σc−σ3Eεc/3σc−2μσ3−γ and Y=α0σc+2σ3E/σc−2μσ3/+σc−σ3E/3σc−2μσ3.


### 2.3.2. Determination ofA and B

As the stress is continuous in the compression and damage evolution section of the rock stress-strain curve, we obtain(31)Eεci1−δ+δexp−α0I1+J2−γFm+2μσ3=Aεci2+Bεci.In addition, the first derivative is also continuous at these two sections; then,(32)∂σ1∂ε1ε=εci=2Aεci+B.Dividing (31) by εci on both sides and integrating,
(33)E1−δ+δexp−α0I1+J2−γFm+2μσ3εci=Aεci+B.By supposingR∣ε=εci=E1−δ+δexp−α0I1+J2−γ/Fm+2μσ3/εci, Z=∂σ1/∂ε1∣ε=εci.Equations (32) and (33) can be rewritten as follows:
(34)Z=2Aεci+B,R=Aεci+B.Solving (34), we have
(35)A=Z−Rεci,B=2R−Z.


### 2.3.3. Determination of Damage Threshold Point (Crack Initiation Point)

In the compression, the original cracks in the rock sample close first under loading. This refers to the key issue in the crack evolution section of the crack initiation and propagation, that is, the damage threshold point. Obviously, it is not easy to obtain the crack initiation stress on a rock stress-strain curve. By researching Lac du Bonnet granite, Martin and Chandler [26] proposed a method to determine the granite initiation point (shown in Figure 3) assuming a crack initiation strength of approximately 40% of the peak strength (σci=40%σc).Figure 3
Crack volumetric strain method (Martin and Chandler [26]).


## 2.3.1. Determination ofm, F, and γ

In the rock triaxial tests, the nominal stressσ1,σ2,σ3σ2=σ3 and strain ε1 can be measured. Based on Hooke’s Law and the effective stress theory, we have
(22)ε1=1Eσ1∗−2μσ3∗,σ2∗=σ3∗=σ21−δD,σ1∗=σ11−δD.Solving (21), σ1∗,σ2∗,σ3∗ can be obtained. Then, I1,J2,J2 is solved:
(23)I1=σ1+2σ3Eε1σ1−2μσ3,J2=σ1−σ3Eε123σ1−2μσ32,J2=σ1−σ3Eε13σ1−2μσ3.In the conventional triaxial test (σ2=σ3) and when ε>εD, (20) can be abbreviated as follows:
(24)σ1=Eε11−δ+δexp−α0I1+J2−γFm+2μσ3.Similarly, parameterγ determines the starting point of the rock damage evolution, that is, the damage evolution threshold point σci,εci. According to fσ−γ≥0 and assuming the damage value is zero here, we obtain
(25)γ=limσ→σciα0I1+J2=α0σci+2σ3Eεciσci−2μσ3+σci−σ3Eεci3σci−2μσ3.By referring the rock stress-strain curve geometry (as shown in Figure2), we know that ①σ=σc,ε=εc;②σ=σc,dσc/dεc=0.Figure 2
Typical stress-strain curve in conventional triaxial compressive test (E is the elastic modulus of rock, σ is the axial stress, and ε is the axial strain).Combining boundary①σ=σc,ε=εc, we obtain
(26)σc=Eεc1−δ+δexp−α0I1+J2−γFm+2μσ3.Substituting boundary②σ=σc,dσc/dεc=0 into the partial derivative of (23), we have
(27)∂σ1∂ε1ε=εc=0.Integrating (26) and (27) yields
(28)α0I1+J2−γFm=InδEεcσc−2μσ3+δ−1Eεc,(29)mα0I1+J2−γm−1Fm=σc−2μσ3εcYσc−2μσ3+δ−1Eεc.Solving (28) and (29), the parameters m and F can be obtained:
(30)m=Xσc−2μσ3εcYσc−2μσ3+δ−1EεcInδEεc/σc−2μσ3+δ−1Eεc,F=XmInδEεc/σc−2μσ3+δ−1Eεc,where X=α0σc+2σ3Eεc/σc−2μσ3+σc−σ3Eεc/3σc−2μσ3−γ and Y=α0σc+2σ3E/σc−2μσ3/+σc−σ3E/3σc−2μσ3.


## 2.3.2. Determination ofA and B

As the stress is continuous in the compression and damage evolution section of the rock stress-strain curve, we obtain(31)Eεci1−δ+δexp−α0I1+J2−γFm+2μσ3=Aεci2+Bεci.In addition, the first derivative is also continuous at these two sections; then,(32)∂σ1∂ε1ε=εci=2Aεci+B.Dividing (31) by εci on both sides and integrating,
(33)E1−δ+δexp−α0I1+J2−γFm+2μσ3εci=Aεci+B.By supposingR∣ε=εci=E1−δ+δexp−α0I1+J2−γ/Fm+2μσ3/εci, Z=∂σ1/∂ε1∣ε=εci.Equations (32) and (33) can be rewritten as follows:
(34)Z=2Aεci+B,R=Aεci+B.Solving (34), we have
(35)A=Z−Rεci,B=2R−Z.


## 2.3.3. Determination of Damage Threshold Point (Crack Initiation Point)

In the compression, the original cracks in the rock sample close first under loading. This refers to the key issue in the crack evolution section of the crack initiation and propagation, that is, the damage threshold point. Obviously, it is not easy to obtain the crack initiation stress on a rock stress-strain curve. By researching Lac du Bonnet granite, Martin and Chandler [26] proposed a method to determine the granite initiation point (shown in Figure 3) assuming a crack initiation strength of approximately 40% of the peak strength (σci=40%σc).Figure 3
Crack volumetric strain method (Martin and Chandler [26]).


## 3. Verification of the Constitutive Model

The validity of the constitutive model was verified by the biotite granite triaxial experiment results under a confining pressure of 25 MPa and at varying temperatures (i.e., 40°C, 60°C, and 130°C).The samples of granite were obtained from the Dali to Ruili railway at the depth 677–682 m (lithology is mainly biotite granite). The average density of the samples was 2.65 g/cm3 at room temperature. Cylindrical specimens were used (diameter: 50 mm, height: 100 mm, Figure 4). The constituent mineralogical components and their compositions were examined by an X-ray diffraction test, and the results are listed in Table 1. The samples did not have distinct, macroscopic heterogeneities (e.g., veins) at the specimen scale. The particle size distribution of various mineral particles is between 0 and 3 mm and mainly in 1 mm (Figure 4).Figure 4
Examples of samples and a generalized photomicrograph of the granite.Table 1
Mineralogical composition of the tested material (weight %).
Plagioclase
Quartz
Biotite
K-feldspar
Chlorite
Tremolite
Pyrite
29
25
20
15
4
4
3The granite triaxial compression tests were carried out on a MTS815 Teststar programmable servo stiffness test machine in the Material Mechanics Laboratory in the Institute of Water Resources and Hydropower, Sichuan University. The heating equipment is a KSW-5D-12 high-temperature box resistance furnace and resistance temperature controller, using a silicon carbon rod as the heating sensor and high-performance fiber as the insulation material. The heating and loading methods are described as follows:(a)
Heating method: heat the samples to the predetermined temperature at a speed of 2°C/min and maintain a constant temperature for 5 h.
(b)
Loading method: apply lateral pressure and axial pressure at the predetermined values simultaneously at a speed of 0.005 mm/s; then, use a 5 mm displacement sensor to measure the axial and circular displacements, and apply an axial load at a speed of 0.005 mm/s until the sample fails completely.The constitutive model parameter values are obtained by combining the granite triaxial compression tests results and the constitutive model parameter expressions (Tables2 and 3).Table 2
Values to calculate constitutive model parameters at different temperatures (25 MPa).
Temperature (T (°C))
Peak strength (σc (MPa))
Strain at peak strength (εc (10−6))
Crack initiation strength (σci (MPa))
Strain at crack initiation (εci (10−6)
Elastic modulus (E (GPa))
Poisson’s ratio (μ)
Constant (α0)
40°C
253.3
4.9
101.3
1.9
5.9
0.27
0.23
60°C
229.5
5.1
91.8
1.6
6.1
0.34
0.21
130°C
263.7
6.2
105.5
1.9
5.6
0.34
0.22Table 3
Parameter values at different temperatures (25 MPa).
Temperature (T (°C))
Quadratic function coefficient
Average strength (F (MPa))
Uniformity coefficient (m)
Damage evolution threshold point (γ)
Damage variable correction coefficient (δ)
A (106)
B
40°C
−3.87
73468.5
233.6
3.18
99.39
0.800
60°C
−6.87
82808.2
229.7
1.87
88.43
0.781
130°C
−4.48
74223.2
274.0
2.09
100.40
0.890The graphs are generated according to the constitutive model expression. Comparisons are made between the theoretical and experimental curves of the granite under different temperatures (40°C, 60°C, and 130°C) under a confining pressure of 25 MPa (Figure5).Figure 5
Experimental and theoretical stress-strain curves at different temperatures (25 MPa).
(a)
(b)
(c)As shown in Figure5, the theory curves fit well to the test curves. This effectively demonstrates the features of the whole process of rock stress-strain, with a good fitting to the curves of different temperatures. This reflects the characteristics of the deterioration and the residual strength of the rock. Adopting the effective stress, which better fits the rock damage rule, the peak strength and its corresponding strain in the theoretical curves are basically identical to those of the experimental curves.


## 4. Discussion

In this section, we would like to explore the relationship between the stress-strain curve and its parameters, as well as understanding the physical meaning of each parameter; we use the constitutive model of brittle granite at 25 MPa and 60°C as an example.A and B are the coefficients of the quadratic function. The larger the A, the greater the degree of curvature of the concave curve, which means that the compaction stage of the stress-strain curve is more obvious. In these curves, B is the initial modulus of the rock sample; γ is the rock damage evolution threshold point; F reflects the rock average strength; m determines the shape of the stress-strain curve and reflects the homogeneity of rock mineral particles, also known as the uniformity coefficient; and δ reflects the residual strength. The complete stress-strain curves under different values of F, m, and δ are shown in Figure 6. The peak stress increases with increasing F. When m=1, no peak point is observed on the curve, whereas both maximum and minimum values are observed at m=3. With increasing δ, the residual declines.Figure 6
Stress-strain curves under different parameter values (25 MPa): (a) features of stress-strain curves with change in parameterF; (b) features of stress-strain curves with change in parameter m; (c) features of stress-strain curves with change in parameter δ (m is the mean uniformity of rock, F is the peak strength of rock, and δ is the damage variable correction coefficient).
(a)
(b)
(c)The excavation of the deeply buried underground openings, the heat generated by decay of radioactive wastes after their emplacements, or the hydraulic fracturing stimulation in geothermal reservoirs may induce damage and progressive failure of the surrounding brittle rocks. The damage accumulation caused by microcracking is therefore considered the main mechanism that leads to the inelastic behavior of stress-strain relation, deterioration of strength and stiffness, and evolution of transport properties (e.g., permeability, diffusivity, and thermal conductivity) of the host rocks. The proper model of coupled thermo-mechanical processes (TM) considering deterioration is hence of paramount importance for analyzing the rock mass response to excavation, thermal loading, and fluid flow. The model could be widely used in underground engineering, such as nuclear waste isolation, production of fossil fuels, underground gas storage, geothermal energy, and deep-buried tunnel.


## 5. Conclusions

The following conclusions can be drawn from this study:(a)
Based on previous works, a statistical damage constitutive model was proposed by a three-parameter Weibull distribution, which considers the crack initiation strength, deterioration, and thermal effects.
(b)
By introducing the boundary conditions and macroscopic parameters, the parameters of the model were established and expressed. The physical significance of the parameters and its effects on the model performance were discussed.
(c)
The model is verified by conventional triaxial experiments under thermal-mechanical actions. A relatively good coincidence between experimental data and theoretical results was found in case studies, which confirms the validity of the model in considering brittle constitutive relations under thermal-mechanical action.It is necessary to point out that our work does not intend to replace any of the existing models because the selection and application of models usually depend on a specific area or purpose. We want to check the possibility of applying this statistical damage model to describe the constitutive behavior of brittle rocks under thermal-mechanical coupling action. The results provide an alternative way of describing similar thermal-mechanical coupling properties of brittle rocks.



---

*Source: 3269423-2018-07-30.xml*