Study on Damage Constitutive Model of High-Concentration Cemented Backfill in Coal Mine

In order to construct the damage constitutive model (DCM) of high-concentration cemented backfill (HCCB) in coal mine, the generalized Hoek-Brown strength criterion was used as the failure criterion. For the difference of theoretical derivation of constitutive relation, a new DCM based on residual strength was proposed. Combined with the conventional triaxial compression test, the correctness and rationality of the DCM were verified. The damage evolution characteristics of HCCB were analyzed, and the physical meaning of model parameters was clarified. The results show that (a) the theoretical curves of stress-strain relation are in good agreement with its experimental curves, which means DCM can simulate the deformation and failure process of HCCB. (b) The damage evolution curve of HCCB is S-shaped. To some extent, the confining pressure can inhibit the development of damage. (c) The parameter F0 reflects the position of the peak point of the DCM, and parameter n is the slope of the straight line segment in the postpeak strain softening stage, which are, respectively, used to characterize the strength level and brittleness of HCCB. The establishment of DCM of HCCB is helpful to reveal its deformation and failure mechanism and provides theoretical basis for its strength design.


Introduction
While coal mining has made great contribution to national economy, it has also brought serious environmental problems. Backfilling mining is one of the effective ways to solve the above problems. With the increasing popularity of environmental protection concept, backfilling mining will be applied more widely. High-concentration cemented backfill (HCCB) in coal mine is a goaf backfilling material composed of coal gangue, fly ash, cement, admixture, and water in a certain proportion. It has good flow characteristics before coagulating and high strength after coagulating. HCCB is mainly used to support the overlying strata in the coal mine goaf to effectively control the surface settlement and reduce all kinds of mining damage [1]. Constitutive relation is one of the key and hot issues in the study of mechanical properties of rock and other materials. The damage constitutive model (DCM) is established by combining the damage mechanics and statistical analysis theory, a stress-strain (σ 1 − ε 1 ) relation which can characterize the deformation and failure pro-cess of rock and other materials, used to analyze and solve the deformation and failure problems of materials [2][3][4][5]. The damage constitutive relation of HCCB is directly related to its stability in goaf and control effect of coal mine surface collapse.
At present, the research on damage constitutive relation is mainly focused on rock materials and has achieved fruitful research results. The reasonable measurement of rock microelements strength is one of the key problems in the establishment of DCM. Many experts and scholars regard the strength failure criterion as the random distribution variable of rock microelement strength. According to the randomness of the distribution of defects in rock materials, the damage in rock deformation and failure process was studied based on the maximum strain criterion and a DCM with simple form and parameter easy to obtain was established [6][7][8][9]. Using the Mohr-Coulomb (M-C) criterion as the expression method of microelement strength and assuming the microelement strength follows Weibull random distribution, a three-dimensional damage statistical constitutive model was established to reflect the postpeak softening characteristics of rock [10][11][12][13][14]. Based on the theory of probability statistics and continuous damage mechanics, the Drucker-Prager (D-P) criterion was introduced as the failure criterion of rock microelements, and the damage evolution equation of rock was deduced strictly, which greatly improved the degree of agreement between the theoretical curve of σ 1 − ε 1 relation and the test data [15][16][17][18][19]. Yan et al. [20], Pan et al. [21], and Zhang et al. [22] used the triple shear energy yield criterion, the unified strength criterion, and the double shear unified strength criterion to characterize the rock microelement strength, respectively, and constructed the DCM based on the above strength criterion. The failure criteria of rock materials involved in the above study include maximum strain criterion, M-C criterion, D-P criterion, triple shear energy yield criterion, unified strength criterion, and double shear unified strength criterion, among which M-C criterion and D-P criterion are in the majority.
As an artificial composite material, HCCB has a short diagenetic time and a large number of soft structural planes such as pores and cracks. The deformation and failure process after loading shows obvious nonlinear characteristics. The generalized Hoek-Brown strength criterion can reflect the nonlinear failure characteristics of rock mass better by comprehensively considering a variety of influencing factors, while D-P criterion is relatively conservative, and M-C criterion is more suitable to represent the linear relation [23][24][25][26][27]. Therefore, the generalized Hoek-Brown strength criterion was selected as the failure criterion of HCCB. In view of the difference of the theoretical derivation of constitutive relation, the theoretical derivation model was modified and a new DCM was proposed. The correctness and rationality of the model were verified by conventional triaxial compression tests. The damage evolution characteristics of HCCB were analyzed. The related studies are helpful to reveal the deformation and failure mechanism of HCCB under triaxial compression and provide a theoretical basis for its strength design.

Test Equipment.
The conventional triaxial compression test was carried out on the RTR-2000 triaxial dynamic test system for high-pressure rock (see Figure 1). The system is mainly used for testing physical and mechanical parameters in normal-and high-temperature and high-pressure environment. The system is equipped with a high stiffness loading frame, the load stiffness is up to 10 mN/mm, the maximum axial pressure is 2000 kN, the maximum confining pressure is 140 MPa, the maximum pore pressure is 140 MPa, and the maximum temperature is 200°C. Samples of different sizes can be tested according to the requirements of users. The sample is a cylinder with a diameter of 25~100 mm and a height of 50~200 mm. The change of test conditions and data acquisition can be completely controlled by a computer. The operator can control the operation of the test by means of load, displacement, strain, and stress.
Axial displacement control was adopted in the test. The axial loading rate was strictly in accordance with standard for test method of concrete physical and mechanical properties (GB/T50081-2019) [28]. During the test, the specimen was first installed in the pressure chamber. Specimens with insufficient end flatness should be polished in advance. Secondly, the specimen was preloaded to make the pressure head contact with the specimen. Then, the confining pressure was applied to the predetermined value according to the hydrostatic pressure condition. The loading rate of confining pressure was 0.25 MPa/min, and the predetermined values were 0 MPa, 1 MPa, 2 MPa, and 3 MPa, respectively. After the confining pressure loading, the specimen was under hydrostatic pressure. Finally, the load was applied along the axis at a constant displacement rate until the specimen was in the residual strength stage. In order to reduce the adverse effects of specimen discreteness on the test and ensure the success of the test, three groups of tests were carried out under each confining pressure, and one group of tests was selected for analysis.

Test Materials.
The coal gangue used in the test was taken from Xinyang Coal Mine in Xiaoyi City, the fly ash was selected from the coal-fired power plant around Xinyang Mine, the cement was commercially available ordinary Portland cement, the admixture was cellulose hydroxypropyl methyl ether, and the preparation water was laboratory tap water. Combined with the practice of backfilling mining in Xinyang Coal Mine, the mass fractions of coal gangue, fly ash, cement, cellulose hydroxypropyl methyl ether, and water were set as 49.94%, 18%, 12%, 0.06%, and 20%, respectively. The weighed raw materials were mixed in turn to prepare the backfilling slurry with mass concentration of 80%. The specimens (ϕ50 × 100 mm) were poured in time and maintained to the prescribed age (28 days). All specimens used the same batch of raw materials. Slurry preparation strictly referred to standard for test method of performance of ordinary fresh concrete (GB/T50080-2016) [29]. The pouring and maintaining of specimens strictly referred to standard for test method of concrete physical and mechanical properties (GB/T50081-2019) [28]. The specimens after maintaining are shown in Figure 2. 2.3. Test Results. The σ 1 − ε 1 curves of HCCB are similar under different confining pressures. Taking the specimen with confining pressure of 2 MPa as an example, the deformation and failure process of HCCB is analyzed. Similar to the deformation and failure process of materials such as 2 Geofluids rocks, the σ 1 − ε 1 curve of the specimens can be divided into five stages, as shown in Figure 3.
(1) Initial compaction stage (before point A): at the initial stage of compression, the original defects inside the HCCB are pressed and closed. The curve of σ 1 − ε 1 is approximately linear. The reason why there is no typical concave section is that the HCCB is so soft that most of the original defects have been compressed in the process of confining pressure loading. The deformation at this stage is close to elastic deformation and can be mainly restored after unloading (2) Elastic deformation stage (AB section): the HCCB bears the load stably. The curve σ 1 − ε 1 is close to a straight line. This stage is dominated by elastic deformation which can be recovered, but also contains a small amount of plastic deformation which cannot be recovered (3) Yield deformation stage (BC section): when the external load reaches the yield stress of the HCCB, new cracks initiate, expand, and gradually connect. The curve of σ 1 − ε 1 is convex, whose slope of tangent decreases until it drops to 0. With the increase of axial strain, the unrecoverable plastic deformation increases gradually (4) Strain softening stage (CD section): when the peak strength is reached, the HCCB is damaged and then enters the strain softening stage. The slope of the tangent of the σ 1 − ε 1 curve changes from zero to negative. The bearing capacity is gradually decreasing. The plastic deformation increases significantly. The specimens show obvious dilatancy characteristics (5) Residual strength stage (after point D): with the increase of axial strain, the axial stress does not change. The curve of σ 1 − ε 1 is a horizontal line. Although HCCB in this stage has been completely destroyed, it still maintains a certain bearing capacity

Establishment of DCM
According to the general pattern of damage mechanics, in order to obtain the DCM of materials, the damage variables are defined in a certain way and damage evolution equation is derived accordingly. Then, the damage evolution equation is substituted into the constitutive relation to obtain the DCM.

Damage Evolution Equation.
Assuming that the number of damaged microelements of the HCCB under a certain level of load is N d , the damage variable D is defined as the ratio of the number of damaged microelements N d to the number of total microelements N [2, 15]; then, When the stress level S reaches the strength F of the microelement, the microelement is destroyed. Assuming that F obeys a certain probability distribution, the number of microelements that fail in any stress level interval ½S, S + dS is as follows: where p is the density function of the probability distribution satisfied by F. For F, the commonly used distribution types include Weibull distribution and normal distribution. When the stress level reaches a certain S, the number of damaged microelements is as follows: where P is the distribution function of the probability distribution satisfied by F. It can be obtained from Equations (1) and (3) that Equation (4) is the deduced damage evolution equation. For any probability distribution, the value of the distribution function changes from 0 to 1 as S increases, which is consistent with the change law of D.

Geofluids
A microelement is taken from any section of the specimen of HCCB, and it is assumed that (1) the microelement conforms to the generalized Hoek theorem, (2) the yield of microelements follows Hoek-Brown criterion, and (3) microelement strength strictly follows Weibull random distribution [8]. Then, the probability density function of the microelement strength is as follows: where F 0 and n are the Weibull distribution parameters. The distribution function of the microelement strength is as follows: It can be obtained from Equations (4) and (6) that According to the assumption of point (2) satisfied by the microelement, S can be determined by the Hoek-Brown criterion. The expression of the generalized Hoek-Brown strength criterion [24] is as follows where σ 1c is the maximum principal stress under rock failure, MPa; σ 3c is the minimum principal stress under rock failure, MPa; σ c is the uniaxial compressive strength of rock, MPa; M is the material parameter, which can be obtained according to the rock mass classification index GSI; S and α are constants related to material characteristics. For materials with good quality, α = 0:5.
The generalized Hoek-Brown strength criterion expressed by the effective stress invariants [24] is as follows: where f ðσ * Þ is the strength of microelement, MPa; I * 1 is the first invariant of effective stress, MPa; J * 2 is the second invariant of effective stress deviator, MPa; θ σ is Lord's angle. For conventional triaxial compression test, θ σ = 30°; σ * i is the effective stress in the direction i, i = 1, 2, 3 MPa.
According to the generalized Hoek law and the concept of effective stress, the corresponding effective stress can be obtained as follows: where σ i is the nominal stress in the direction i, MPa; E is the elastic modulus, GPa; ε 1 is the axial strain, %; μ is Poisson's ratio. Substituting Equation (10) into Equation (9), the following can be obtained.
According to Hoek's law, the principal strain can be as follows: It can be obtained from Equations (12) and (13) that According to Equation (14), when the HCCB is completely damaged, that is, D = 1, σ 1 = μðσ 2 + σ 3 Þ. Then, the residual strength of HCCB is μðσ 2 + σ 3 Þ. However, the actual residual strength of HCCB is σ r . There is an absolute difference in value between them. In order to eliminate this difference, make the established constitutive model more in line with the reality and better characterize the σ 1 − ε 1 relation of HCCB. Equation (14) is modified based on the σ r of HCCB. A new constitutive model of HCCB is proposed, that is, where σ m is the difference between σ r and μðσ 2 + σ 3 Þ, that is,

Geofluids
In the conventional triaxial compression test, σ 2 = σ 3 , the DCM of HCCB can be expressed as follows: It can be seen from Equations (16) and (7) that the determination of distribution parameters F 0 and n is one of the key issues in the establishment of DCM. The unknown parameters of the model can be determined according to the extreme value characteristics of the σ 1 − ε 1 curve. The specific determination process is as follows.
Assuming that the stress and strain at the peak point corresponding to σ 1 − ε 1 curve under different confining pressures are σ 1c and ε 1c , respectively, then, The derivative of Equation (16) can be obtained.
According to Equations (17) and (18), the following can be obtained.
where D 1c is D when σ 1 = σ 1c and ε 1 = ε 1c , and it can be determined from Equation (16), that is, Partial differentiation of Equation (7) can be obtained.
Partial differentiation of Equation (11) can be obtained.
The transformation of Equation (7) can be obtained.
It can also be expressed as follows:

By substituting Equations (22)-(24) into Equation (21) and then substituting Equation (21) into Equation (19), it can be obtained that
where S 1c is S when σ 1 = σ 1c and ε 1 = ε 1c , and it can be determined from Equation (11), that is, By substituting D = D 1c and S = S 1c into Equation (23), it can be obtained that

Model Verification
According to conventional triaxial compression test with σ 3 = 0 MPa, the uniaxial compressive strength of HCCB was 8.5 MPa, that is, σ c = 8:5 MPa. The material parameter m of the generalized Hoek-Brown strength criterion was fitted by the least square method and m = 5:58. Combined with conventional triaxial compression test results, DCM parameters under different confining pressures were obtained, as shown in Table 1.
In order to verify the correctness and rationality of the established DCM, the obtained model parameters were replaced into the DCM to draw the theoretical curve of σ 1 − ε 1 relation. The theoretical curves were compared with the test curves under different confining pressures, as shown in Figure 4.
As can be seen from Figure 4, under low confining pressure, the postpeak curve of σ 1 − ε 1 relation has an obvious strain softening stage. The postpeak axial stress decreases significantly with the increase of axial strain, as shown in Figures 4(a)-4(c). When the confining pressure reaches a certain value, the postpeak curve of σ 1 − ε 1 relation no longer has obvious strain softening stage. After the specimen reaches the peak strength, with the increase of axial strain, the axial stress decreases but not significantly, as shown in Figure 4(d). In Figure 4(d), the peak strength is 16.6 MPa and the residual strength is 16.3 MPa, which are very similar, but the postpeak strain softening stage still exists. For the σ 1 − ε 1 relation whose postpeak curve has obvious strain softening stage, its theoretical curve agrees well with the experimental curve, especially in the prepeak stage and the residual strength stage. In the strain softening stage, the theoretical curve and the test curve are in relatively low agreement, showing a certain difference. However, in Figures 4(a)-4(c), the maximum difference between theoretical and experimental values in the strain softening stage is 5 Geofluids 1.35 MPa. Compared to local differences of tens of megapascals in the rocks, the degree of match has been greatly improved. For the σ 1 − ε 1 relation whose postpeak curve has no obvious strain softening stage, there is a great difference between the theoretical curve and the experimental curve, especially in the prepeak stage. The maximum differ-ence between theoretical value and experimental value is 7 MPa. From the initial compaction to the peak strength, the difference gradually decreases between them. In the postpeak stage, the theoretical curve can simulate the experimental curve well. In view of the relatively large difference between the theoretical value and the experimental value in   6 Geofluids the prepeak stage, it is considered that the damage constitutive equation established can only characterize the σ 1 − ε 1 relation whose postpeak curve has no obvious strain softening stage to a certain extent. When a DCM with a higher degree of consistency is needed, its damage constitutive relation can be reconstructed. When the confining pressure is 0 MPa, the postpeak axial stress decreases rapidly with the increase of axial strain. The absolute value of the slope of the straight segment of the postpeak curve is large. The difference between peak strength and residual strength is large. HCCB shows obvious brittleness. With the increase of confining pressure, the axial stress after the peak decreases slowly with the increase of axial strain. The absolute value of the slope of the straight section of the postpeak curve is small. The difference between the peak strength and the residual strength is reduced. The brittleness of the HCCB weakens, showing a certain plasticity. When the confining pressure is 3 MPa, the HCCB shows the lowest brittleness. The DCM can reflect the trend of brittleness decreasing and plasticity increasing.

Analysis of Damage Evolution Characteristics
and Parameter Effect 5.1. Damage Evolution Characteristics. The theoretical curves of σ 1 − ε 1 relation are in good agreement with the experimental curves, which indicates that the established DCM is correct and the damage evolution equations are reasonable. Based on damage evolution equations, the damage evolution characteristics in the process of deformation and failure of HCCB will be discussed. The evolution equation of damage variable is shown in Equation (28).
According to Equation (28), combined with the parameters of DCM, damage evolution curves of HCCB under different confining pressures can be drawn, as shown in Figure 5.
Most of the damage evolution curves of rock materials show S-shaped [2,24]. The damage evolution curves increase monotonically under different confining pressures. With the increase of axial strain, the damage variable increases from 0 and gradually approaches 1, completing the continuous accumulation of damage until it reaches saturation. The whole damage evolution process can be divided into no damage stage, damage beginning stage, damage accelerating stage, damage slowing stage, and damage stabilizing stage. The corresponding deformation and failure stages are initial compaction and elastic deformation stage, yield deformation stage, early strain softening stage, late strain softening stage, and residual strength stage. As can be seen from Figure 5, the damage evolution curves of HCCB are half S-shaped and are the upper part of S. The damage evolution curves mainly include the damage slowing stage and damage stabilizing stage. The reasons are as follows: due to the confining pressure, the starting point of the damage evolution curve and the stress-strain curve is not the origin, resulting in no damage stage, damage beginning stage, and damage accelerating stage that are not shown in Figure 5. The dashed line in Figure 5 is the dividing line between the damage slowing stage and the damage stabilizing stage. Among them, the damage stabilizing stage corresponds to the residual strength stage of HCCB, and the damage slowing stage almost corresponds to all deformation and failure stages except the residual strength stage. The reason for the above differences is that the strength of cemented matrix forming HCCB is low. The cemented matrix is rapidly damaged from the beginning of axial loading. In the damage slowing stage, the damage of HCCB increases rapidly at first and then slowly. The increase rate gradually slows down. In the early part of this stage, most of the damage of the HCCB is completed within a small axial strain range. In the later part of this stage, the damage increment is limited, but this part lasts for a long time, which is several times of the early part. In the damage stabilizing stage, the damage of HCCB does not increase any more. D is always 1, which means that the HCCB has been completely damaged. Figure 5 shows that the damage slowing stage of HCCB will last for a long time no matter under low or high confining pressure. The axial strain of HCCB is maintained at a high value when it is completely damaged. When D = 1, the HCCB is completely damaged. The axial strain under complete damage is defined as ε r . The curve of ε r changing with σ 3 is shown in Figure 6.
It can be seen from Figure 6 that ε r increases with the increase of σ 3 . There is a roughly linear relation between them; that is, the greater the confining pressure is, the greater the axial strain is when the HCCB is completely damaged. The main reason is that confining pressure can inhibit the development of damage to some extent. With the increase of confining pressure, when HCCB achieves the same damage amount, the axial deformation required is larger.

Parametric Effect.
Most of the parameters in the DCM have clear physical meanings. For example, E is the elastic modulus, reflecting the elastic level of the material. At the same time, there are also some parameters whose physical meanings are unclear, such as Weibull distribution parameters F 0 and n. The following is to determine the physical significance of F 0 and n by analyzing their influence on DCM. Weibull distribution parameters F 0 and n are the common reflection of material characteristics and confining pressure of HCCB. Under the same confining pressure, F 0 and n are bound to be different for HCCB with different proportions. Similarly, for HCCB with the same proportion, F 0 and n are different under different confining pressures and show certain regularity. It can be seen from Table 1 that with the increase of confining pressure of HCCB, both F 0 and n of the DCM decrease.      9 Geofluids stage. Therefore, F 0 reflects the location of the peak point of the DCM, which can be regarded as a physical quantity representing the strength level of HCCB.

Effect of n on DCM.
Based on the specimen with confining pressure of 2 MPa, other parameters remain unchanged and only parameter n is changed. When n varies approximately in amplitude of 0.10, the DCM was obtained under different values of n, as shown in Figure 8. Figure 8 shows that the linear elastic deformation stage of the DCM is almost not affected by n. The linear elastic deformation stages under different n are approximately coincident. With the increase of n, the peak strength, peak strain, and residual strength of the DCM all decrease. The peak point of the DCM moves inward and is closer to the origin in the horizontal and vertical directions. On the other hand, the shape of DCM changes obviously. N changes the softening and steepening degree of the straight line segment in the strain softening stage. As n increases, the line segment becomes steeper and the brittleness of HCCB increases. Therefore, n reflects the slope of the straight section in the strain softening stage, which can be regarded as a physical quantity representing the brittleness degree of HCCB.

Conclusions
The following conclusions can be drawn from this study: (a) For the σ 1 − ε 1 relation whose postpeak curve has obvious strain softening stage, the theoretical curve agrees well with the experimental curve. The DCM can characterize the deformation and failure process of HCCB. For the σ 1 − ε 1 relation whose postpeak curve has no obvious strain softening stage, the damage constitutive relation with higher degree of coincidence can be established according to the demand (b) The damage evolution curve of HCCB is S-shaped.
The damage evolution process only shows damage slowing stage and damage stabilizing stage because of the confining pressure. When the HCCB is completely damaged, the axial strain remains high. Confining pressure can inhibit the development of damage to some extent  10 Geofluids (c) The linear elastic deformation stage of the DCM is almost not affected by the distribution parameters F 0 and n. F 0 reflects the position of the peak point of DCM and n reflects the slope of the straight line segment in the strain softening stage, respectively, used to characterize the strength level and brittleness degree of HCCB

Data Availability
The data used to support the findings of this study are available from the corresponding author and first author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.