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

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.


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][5][6][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][9][10][11][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][16][17][18]. Cao and coworkers [19][20][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][25][26][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][29][30][31][32]. The factors of porosity, temperature, or freezethaw 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.  [35]. Consequently, the definition of "continuity" was proposed to quantitatively show the damage behavior of material under loading in Eq. (1), i.e., where ψ is the continuity of material, S is the initial total area of the material,andS 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.
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 coordinate i-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 N i0 . The damage evolution of the plane can be evaluated under threedimensional 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 2 Geofluids are progressive failure, which results in the rock damage. The number of damaged elements in the plane under loading stress is N id (Figure 1(b)). 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 (N id ) to the initial element number (N i0 ) 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 D i of the plane which is perpendicular to the coordinate i -axis can be given as 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 PðF i Þ in Figure 1(c) can be expressed as Initial mesoscopic element Damaged mesoscopic element

Geofluids
where F i is the strength of mesoscopic element, and w i and F i0 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 Combining Eqs. (3) and (5), the areal damage variable of the plane which is perpendicular to the i-axis can be given by

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 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 where ε i is the component of three principal strains, and E i 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 Substituting Eq. (6) into (9), a new statistical damage constitutive model of anisotropic rock (SDAR model) is developed as The Drucker-Prager failure criterion is widely used in the researches of rock mechanics [38][39][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 (F i ) in the plane which is perpendicular to the coordinate i -axis can be given by where α dp refers to the parameter of material, I 1 is the first invariant of stress, and J 2 is the second invariant of stress deviator.
The three parameters of the Drucker-Prager failure criterion can be calculated by α dp = sin φ/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9 + 3 sin 2 φ q , where φ refers to the friction angle, and σ * 1 , σ * 2 , and σ * 3 are the effective stress in three directions.

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., w i and F i0 , but the values of w i and F i0 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., Drilling angle 90º Drilling angle 0º Geofluids The failure strength of element (F i ) can be regarded as function of σ 1 , σ 2 , σ 3 , and ε i . The subscript i of ε i corresponds to that of F i . Then, the following Eq. (14) can be obtained by differential of F i .
The parameters of w and F 0 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 w i and F i0 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 where σ i = 1/2ðσ j + σ k Þ. Substituting Eqs. (14) and (15) into Eq. (13) yields where Solving Eq. (16) based on the Kramer's law yields where Δ = jU ij j = , and Δ ij is the remainder of Δ, i.e., 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 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., 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), , For the reason that Δ ii ≠ 0, the equation that X i = 0 can be obtained. By solving X i = 0, the calculated expressions of w i and F i0 can be derived as Eqs. (23) and (24) denote the determination of parameters in the SDAR model, i.e., w i and F i0 . 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 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 F 0 remains unchanged.
Equations (25)- (27) are the calculated expressions of w and F 0 obtained by the degradation of the calculation expressions of w i and F i0 for isotropic rock under conventional triaxial and uniaxial compression, respectively. The calculated expressions of w and F 0 under conventional triaxial compression condition of isotropic rock (Eqs. (25) and (26)) are consistent with those derived by Wang et al. [42]. Further-more, 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]. 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.

Model Validation
Comparing the theoretical results of the SDAR model and test data of shale on stress-strain in Figure 3, 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.
Comparing the peak strength of shale in the same drilling direction under each confining pressure (Figure 3), 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         9 Geofluids 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.

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   (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]. 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.
The theoretical results of the SDAR model and test data of coal on stress-strain are shown in Figure 5. 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.
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.  11 Geofluids 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.
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.

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 w i and F 0i .  (σ 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.
There are different stress-strain relationships and damage evolution under various combinations of elastic modulus. Figure 7 shows the results of stress-strain and damage under different anisotropic elastic moduli.
Comparing the cases 1, 2, and 3 in Figure 7(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 Figure 7(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 Figure 7, we can know that the effects of elastic modulus along the bedding plane direction on the variation of peak stress, stressstrain 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 elas-tic 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.

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 the i-axis direction is only influenced by Poisson's ratio in the other two directions (v ij and v ik ) 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.
The stress-strain and damage of coal under different Poisson's ratio are shown in Figure 8. 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)).  Table 5. On the premise that other parameters remain unchanged, with the increase of w i , 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 w i , 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 w i increases.

4.2.2.
Parameter of F i0 . In order to analyze the sensitivity of F i0 , the values of F i0 are changed according to Table 6. The stressstrain and damage variable can be obtained in different values of F i0 . As it can be seen in Figure 10(a), with the increase of F i0 , 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 F i0 . Meanwhile, the axial strain corresponding to the emergence of initial damage becomes larger with the increase of F i0 (Figure 10(b)). It indicates that higher loading stress is required for emergence of coal damage. In addition, the larger the F i0 is, the larger the axial strain corresponding to the damage variable reaching 1 becomes.

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 stressstrain 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 of w i . Meanwhile, the peak strength, peak strain, and rate of stress dropping after peak point increase with the increased value of w i . Moreover, the rate of damage accumulation increases as well. On the basis of increasing F i0 , 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