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

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.


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 waterrock 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][7][8][9][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][12][13][14][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][20][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 damagebased 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.

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 D T can be given by where E T and E 0 are the elastic moduli at T °C and room temperature, respectively.

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: 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).
Considering the thermal effects, the Weibull parameters can be given by where m 0 , F 0 , and γ 0 are the Weibull parameters at room temperature and m, F, and γ are as defined previously.

Geofluids
Rock damage is caused by the continuous failure of the microunits within the rock.Assuming that there are N f 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, Then, at stress f σ , the failure number N f is Substituting ( 2) and ( 4) into (3) yields Using the Drucker-Prager criterion as the unit's yield criterion, where α 0 = sin φ/ 9 + 3 sin 2 φ, I 1 is the first stress invariable, and J 2 is the second invariant of deviator stress.These are expressed as follows: Therefore, ( 5) can be written as follows: 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: where D σ, T is the damage value considering the thermal effect and σ ci is the crack initiation strength.

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 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 Then, According to the generalized Hooke's law, we have Substituting ( 12) into (15), we obtain 14) leads to According to (1), the elastic modulus of hard rock under different temperatures is given as follows: Substituting ( 16) into (15), the hard rock damage constitutive equation considering thermal effects is obtained (after the damage threshold point): 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: Integrating ( 17) and ( 18), the hard rock constitutive equation with the thermal-mechanical coupling is given as follows: where ε ci is the strain value at the damage threshold point.

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: 2.3.1.Determination of m, 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 3 can be obtained.Then, I 1 , J 2 , J 2 is solved: In the conventional triaxial test (σ 2 = σ 3 ) and when ε > ε D , (20) can be abbreviated as follows: 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 By referring the rock stress-strain curve geometry (as shown in Figure 2), we know that ①σ = σ c , ε = ε c ; ②σ = σ c , dσ c /dε c = 0.

Determination of A and B.
As the stress is continuous in the compression and damage evolution section of the rock stress-strain curve, we obtain In addition, the first derivative is also continuous at these two sections; then, Dividing (31) by ε ci on both sides and integrating, Equations ( 32) and ( 33) can be rewritten as follows: Solving (34), we have 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 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 ).

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/cm 3 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).
The 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 hightemperature box resistance furnace and resistance temperature controller, using a silicon carbon rod as the heating    6 Geofluids 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
the constitutive model parameter expressions (Tables 2  and 3).

Geofluids
The 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 (Figure 5).
As shown in Figure 5, 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.

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.
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.

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.

Figure 2 :
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).

Figure 4 :
Figure 4: Examples of samples and a generalized photomicrograph of the granite.

Figure 6 :
Figure 6: Stress-strain curves under different parameter values (25 MPa): (a) features of stress-strain curves with change in parameter F; (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).

Table 1 :
Mineralogical composition of the tested material (weight %).

Table 2 :
Values to calculate constitutive model parameters at different temperatures (25 MPa).