Experimental Investigation and Micromechanical Modeling of Hard Rock in...

[Full title: Experimental Investigation and Micromechanical Modeling of Hard Rock in Protective Seam Considering Damage–Friction Coupling Effect]

Citation: Zhang, C.; Liu, W.; Shi, C.; Hu, S.; Zhang, J. Experimental Investigation and Micromechanical Modeling of Hard Rock in Protective Seam Considering Damage–Friction Coupling Effect Sustainability 2022 , 14 , 16296. https://doi.org/10.3390/ su 142316296 Academic Editors: Danqing Song, Zhuo Chen, Mengxin Liu and Yutian Ke Received: 5 October 2022 Accepted: 2 December 2022 Published: 6 December 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations Copyright: © 2022 by the authors Licensee MDPI, Basel, Switzerland This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/) sustainability Article Experimental Investigation and Micromechanical Modeling of Hard Rock in Protective Seam Considering Damage–Friction Coupling Effect Chuangye Zhang 1 , Wenyong Liu 2 , Chong Shi 2 , Shaobin Hu 2 and Jin Zhang 2, * 1 Pingdingshan Tianan Coal. Mining Co., Ltd., No. 12 Mine, Pingdingshan 467099, China 2 Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing 210024, China * Correspondence: zhangjin 90@hhu.edu.cn Abstract: The hard rock in the protective coal seam of the Pingdingshan Mine in China is a typical quasi-brittle material exhibiting complex mechanical characteristics. According to available research on the mechanical property, the inelastic deformation and development of damage are considered related with crack initiation and propagation, which are main causes of the material degradation. In the present study, an original experimental investigation on the rock sample of the Pingdingshan coal mine is firstly carried out to obtain the basic mechanical responses in a conventional triaxial compression test. Based on the homogenization method and thermodynamic theory, a damage– friction coupled model is proposed to simulate the non-linear mechanical behavior. In the framework of micromechanics, the hard rock in a protective coal seam is viewed as a heterogeneous material composed of a homogeneous solid matrix and a large number of randomly distributed microcracks, leading to a Representative Elementary Volume (REV), i.e., the matrix–cracks system. By the use of the Mori–Tanaka homogenization scheme, the effective elastic properties of cracked material are obtained within the framework of micromechanics. The expression of free energy on the characteristic unitary is derived by homogenization methods and the pairwise thermodynamic forces associated with the inelastic strain and damage variables. The local stress tensor is decomposed to hydrostatic and deviatoric parts, and the effective tangent stiffness tensor is derived by considering both the plastic yield law and a specific damage criterion. The associated generalized Coulomb friction criterion and damage criterion are introduced to describe the evolution of inelastic strain and damage, respectively. Prepeak and postpeak triaxial response analysis is carried out by coupled damage–friction analysis to obtain analytical expressions for rock strength and to clarify the basic characteristics of the damage resistance function. Finally, by the use of the returning mapping procedure, the proposed damage– friction constitutive model is applied to simulate the deformation of Pingdingshan hard rock in triaxial compression with respect to different confining pressures. It is observed that the numerical results are in good agreement with the experimental data, which can verify the accuracy and show the obvious advantages of the micromechanic-based model Keywords: micromechanical model, damage-friction coupling, rock mechanics, experimental investigation 1. Introduction An increasing number of large rock projects such as the Sichuan-tibet Railway tunnel project [ 1 ] and Beishan high-level nuclear waste geological disposal underground project [ 2 ] are currently being planned and developed worldwide. In order to guarantee the safety and stability of the rock mass engineering, the mechanical properties of the engineering rocks need to be mastered. The high ground stress, high gas and low permeability conditions that exist in protruding mines at a depth of 1000 meters make protrusion accidents a frequent problem [ 3 , 4 ]. The rock (as shown in Figure 1 ) here is drilled and blasted out of the Sustainability 2022 , 14 , 16296. https://doi.org/10.3390/su 142316296 https://www.mdpi.com/journal/sustainability

Sustainability 2022 , 14 , 16296 2 of 17 coal mine tunnels and made into standard rock samples. The corresponding mechanical parameters can be obtained by experiment and theoretical simulation. This method can provide reference for the calculation of deformation and stability of the coal mine [ 5 , 6 ]. In this paper, the hard rock in the protective coal seam of the Pingdingshan mine is studied experimentally and analytically in order to provided a unified constitutive model for predicting its mechanical behaviors based on the obtained experimental data Figure 1. Hard rock samples obtained from the protective coal seam of the Pingdingshan Mine The horizontal shafts of the Pingdingshan coal mine are located a thousand meters deep underground. The hard rock is collected from the protective coal seam, which is a typical discontinuous, anisotropic and quasi-brittle material. It has a complex microstructure, with a high number of fractures and microcracks in its internal structure. In recent years, to obtain the basic mechanical properties of various rocks, the conventional triaxial compression tests are accepted as a powerful method by many scholars on this subject [ 7 – 9 ]. Over the past three decades, significant progress has been made in the modeling of plastic degradation in quasi-brittle materials [ 10 – 12 ]. Several theoretical frameworks have also been proposed for plastic damage models, including [ 13 – 17 ], only to name a few. For concrete and other related materials, isotropic and anisotropic damage models have been developed with or without the plastic coupling, for example, [ 18 – 23 ]. For rock materials, certain models [ 24 – 27 ] have been developed. On the other hand, certain discrete plasticity–damage models [ 28 – 31 ] have been developed with the help of the micro-plane theory and the discrete thermodynamics formulation in order to better understand the consequences of anisotropic distribution of microcracks in brittle materials [ 32 – 34 ]. Furthermore, certain discrete plasticity–damage models have been developed with the help of the micro-plane theory and the discrete thermodynamics formulation [ 35 – 37 ] in order to better understand the consequences of anisotropic distribution of microcracks in brittle materials [ 38 – 40 ]. Recently, within the framework of micromechanics, elastoplastic damage modeling has been successfully applied on rock materials by considering penny-shaped cracks [ 41 , 42 ]. In this paper, based on the available micromechanic-based modeling method [ 43 – 45 ] , a new multiscale constitutive model simulating the mechanical properties of the hard rock in the protective coal seam of the Pingdingshan coal mine is constructed by combining homogenization theory and irreversible thermodynamics. The discussion focuses on strength prediction and parameter determination based on coupled damage–friction analysis, and new damage criteria are proposed based on the characterization of the damage evolution resistance function. The model is illustrated by a returning mapping procedure to simulate the mechanical behavior of the hard rock in conventional triaxial compression tests with respect to different confining pressures. It is emphasized that for the evolution of inelastic strains, the deformation can be still accurately predicted by the constructed model despite the application of the associated flow law. As one of the prominent advantages of the multiscale model, by using the homogenization method to obtain the full expression of

Sustainability 2022 , 14 , 16296 3 of 17 the free energy on the characteristics of a unit cell to establish the basic pattern of damage friction coupling, the strengthening of the local stress contains a similar back stress weakening function, and the model parameters are greatly reduced. The physical meaning is clear and provides a convenient process for the parameter calibration and the engineering application, which are also demonstrated in this paper 2. Experimental Investigation of the Hard Rock in Protective Coal Seam of Pingdingshan Mine 2.1. Sample Preparation and Testing Producure The tested samples in this study are obtained from the protective coal seam of the Pingdingshan Mine in China. The type of rock is diorite. The overall rock samples are light gray surrounded by black particles (as shown in Figure 1 ). Rock test specifications require standard specimens to be cylindrical, with a diameter of 50 mm and a range of 48 to 52 mm allowed. The height is 100 mm and the allowable range is 95 to 105 mm. For rocks with a heterogeneous coarse-grained structure, non-standard samples are allowed, but the ratio of height to diameter should be within 2 to 2.5. After mining by blasting method, the rock is processed into standard samples with a diameter of 50 mm and a height of 100 mm The non-parallelism error of the two end faces of the rock sample is less than 0.05 mm, the end faces are perpendicular to the axis of the rock sample, and the deviation is less than 0.25 ◦ . No obvious cracks are observed on the outer surface, and a good homogeneity is seen. The average density of the samples is 2.6 g/cm 3 The conventional triaxial compression tests were carried out at the Multiscale Multi- Field Coupled Rock Mechanics Laboratory of Hohai University, using a triaxial rheometer manufactured by TOP Industrie France (see Figure 2 ). The equipment mainly consists of a triaxial pressure chamber, an axial pressure servo pump, a perimeter pressure servo pump and a computer control system, which can realize conventional triaxial compression tests on rocks, with a wide range of application and high measurement accuracy. Pressure control adopts a high-precision electronic control servo high pressure pump, and the measurement accuracy can reach 0.01 MPa; two highly sensitive displacement sensors are used for axial displacement measurement, which can directly output the axial displacement value of the tested sample. The utilized loading technique is strain-prescribed load. The strainprescribed load requires the instrument’s indenter to press down on the rock at a rate of 0.02 mm/min. The lateral strain deformation measuring device (the right subfigure of Figure 2 ) here is directly placed on the rubber sleeve outside the rock sample when in use. During the test, the computer can directly read the corresponding deformation of the measuring device. The measurement range is 20 mm, and the measurement accuracy is 10 − 3 mm. The system can produce the pressure in many ways, among which the axial pressure can be carried out by axial displacement control, pressure control, flow control, lateral displacement control and other loading methods Figure 2. Autonomous and auto-compensated multi-field coupling testing system 2.2. Experimental Results Figure 3 illustrates the curves of deviatoric stress ( σ 1 − σ 3 ) versus the axial and lateral strains at various confining pressures (confining pressure P c = 0, 10, 20 and 30 MPa).

Sustainability 2022 , 14 , 16296 4 of 17 The rock of the protective coal seam is noted to exhibit the characteristic mechanical properties of brittle solids. According to the linear portion of the curves, the mean Young’s modulus E = 29 GPa and the Poisson’s ratio ν = 0.08 are calculated -1.5 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 250 Deviatoric stress (MPa) Figure 3. The stress–strain curves of the hard rock in a protective coal seam with respect to different confining pressures 3. Construction of the Damage–Friction Coupling Constitutive Model 3.1. Theoretical Framework of Elastic Damage Model The initiation, propagation and connection of microcracks within the solid matrix are the main mechanisms leading to the regression of mechanical property and material failure. Quasi-brittle rocks with microcracks can be regarded as a kind of composite material and can be therefore studied as a Representative Elementary Volume (REV) containing a solid matrix and a large number of penny-shaped microcracks. By experimental observation, the fracture of the rock sample in conventional compression tests is mainly caused by the crack initiation and propagation. Thus, only the reduction of effective properties by cracks is considered and that by pores is ignored. Based on this, the REV is viewed as a heterogeneous system with scattered microcracks and a solid matrix that has been weakened by pores (Figure 4 ). The following derivations are studied within the framework of homogenization methods S.-L. Liu, H.-R. Chen and S.-S. Yuan et al. International Journal of Mechanical Sciences 179 (2020) 105654 Fig. 5. Representative elementary volume of microcracked solid of 60 to 120 MPa, the macroscopic dilatant cataclastic flow with no shear bands was observed, resulting from the interactions between microcracking and dislocation activities The stress concentration at crack tips is relaxed by dislocation activities within the neighboring grains and the ability of cracks is insufficient to propagate to localized into shear bands In this scenario, shear-oriented interfacial microcracks nucleate and form the grain-scale faults, leading to macroscopic distributed deformation In addition, with the increase of confining pressure, the dominant mechanism of microcracking changes from tensile cracks to compressive-shear cracks 4. Constitutive modelling and numerical simulation In view of the operative deformation mechanism related to microcracking, the damaged sandstone can be modeled as a heterogeneous system composed of an elastic matrix and a number of randomly oriented and distributed microcracks For upscaling analyses, let us consider the representative elementary volume in Fig 5 to represent a macroscopic material point Within the framework of linear homogenization, the elasticity tensors for the matrix and microcracks are denoted respectively by ℂ m and ℂ c all microcracks are geometrically assumed to be penny-shaped and each microcrack is identified by its unit normal vector n and aspect ratio {?} = {?} ∕ {?}, where a and c are the average crack radius and crack aperture, respectively The dimensionless crack density parameter, {?} =  {?} 3 is adopted as internal damage variable, where  denotes the number of cracks per unit volume The focus of this work is on the mechanical behaviors of the sandstone under triaxial compression with the increase of deviatoric stress, the rock will experience inelastic deformation due to frictional sliding along closed crack surfaces This deformation can be treated as plastic strain The total strain {?} can thus be decomposed into an elastic reversible strain in the matrix {?} e and a plastic strain due to frictional sliding {?} p {?} = {?} e + {?} p (1) 4.1. Free energy and state equations For the heterogeneous material described above, its free energy has been formulated by following the linear homogenization procedure [32,49] {?} = 1 2 {?} e ∶ ℂ m ∶ {?} e + 1 2 {?} p ∶ ℂ p ∶ {?} p (2) Under isotropic assumption, ℂ m = 3 {?} m {?} + 2 {?} m {?} with k m and {?} m being the bulk modulus and shear modulus of the solid matrix, respectively The full determination of the fourth order modulus tensor ℂ p depends on the choice of homogenization scheme In the literature, various homogenization schemes have been proposed to approximate the effective properties of composite solids When microcracking problems are concerned, the Mori-Tanaka (MT) homogenization method shows close links with the results established in the Linear Elastic Fracture Mechanics [23,32] According to the Mori-Tanaka homogenization method [5,29] , ℂ p involved in the stored energy takes the following explicit expression [49] ℂ p = 1 {?} 1 {?} 3 {?} m {?} + 1 {?} 2 {?} 2 {?} m {?} (3) where {?} 1 = 16 9 1−( {?} m ) 2 1−2 {?} m and {?} 2 = 32 45 (1− {?} m )(5− {?} m ) 2− {?} m are two constants only dependent on the Poisson’s ratio {?} m of the solid matrix Given the strain free energy (2) , one can derive the relation between the macroscopic stress {?} and the macroscopic strain {?} {?} = {?}{?} {?} {?} = ℂ m ∶ ( {?} − {?} p ) , (4) and the thermodynamic force associated with the plastic strain {?} p {?} c = − {?}{?} {?} {?} c = {?} − ℂ p ∶ {?} p (5) as well as the damage-conjugated thermodynamic force {?} {?} = − {?}{?} {?}{?} = − 1 2 {?} p ∶ {?} ℂ p {?}{?} ∶ {?} p (6) 4.2. Plasticity-damage coupling analyses A generalized Coulomb-type friction criterion in terms of the local stress {?} c is adopted to characterize the frictional sliding along microcracks, such that {?} = ‖‖ {?} c ‖‖ + {?}{?} c m = 0 (7) where {?} c = {?} ∶ {?} c and {?} c m = 1 3 {?} c ∶ {?} are the deviatoric part and the mean value of the local stress {?} c , respectively {?} denotes the generalized coefficient of friction In the framework of plasticity theory, the evolution rate of {?} p is determined by following the normality rule ̇ {?} p = {?} p {?} where {?} p is the non-negative plastic multiplier and {?} = {?}{?} {?} {?} c defines the plastic flow direction It has been importantly noted that for the particular case of conventional triaxial compression loading path, the flow direction tensor D is unchanged in the whole loading process and takes the form {?} = − 1 √ 6 diag ( 2 , −1 , −1 ) + 1 3 {?} diag ( 1 , 1 , 1 ) [49] It is then possible to define the accumulated plastic multiplier Λ p = ∫ {?} p and the total plastic strain is rewritten in the form {?} p = {?} p {?} (8) By benefiting from this salient feature, the friction criterion is cast into the form {?} = ‖ {?} ‖ + {?}{?} m − 2 Λ p {?} {?} = 0 (9) with {?} = {?} m {?} 2 2 {?} 1 + {?} m {?} 2 {?} = {?} ∶ {?} and {?} m = 1 3 {?} ∶ {?} are the deviatoric part and the mean part of the macroscopic stress {?} , respectively Figure 4. Representative Elementary Volume of hard rock considering microcracks Consider an isotropic rock matrix where the fourth-order elastic stiffness tensor is expressed as: C m = 3 k m J + 2 µ m K (1)

Sustainability 2022 , 14 , 16296 5 of 17 where k m and µ m are the compressive bulk modulus and shear modulus of the matrix, respectively. By introducing the second-order unit tensor δ , the fourth-order tensor operators J and K , they are denoted as J ijkl = δ ij δ kl /3, K ijkl = ( δ ik δ jl + δ il δ jk ) /2 − δ ij δ kl /3 The existence of microcracks leads to the discontinuity of the displacement field in rock material. Therefore, the macroscopic strain ε can be decomposed into two components, namely, the elastic strain ε m and the inelastic strain ε c , corresponding to the solid matrix component and the microcrack component, respectively: ε = ε m + ε c (2) Then, the relationship between macroscopic stress σ and strain ε is obtained as follows: σ = C m : ( ε − ε c ) (3) The degree of shear expansion of the microcrack and the relative slip of the microcrack face are expressed in terms of the scalar β and the second-order tensor T , respectively Thereby, one has the inelastic strain tensor ε c as the sum of the hydrostatic and deviatoric parts: ε c = T + 1 3 β δ , β = tr ε c (4) The construction of a constitutive model in terms of damage mechanics generally includes three steps: firstly, select a suitable damage variable ϕ to describe the damage state of the material. Budiansky and Connell [ 46 ] proposed that the damage variable is related to the fracture density, namely, ϕ = Nd 3 (where N is the number of microcracks per unit volume and d is the radius of the coin crack surface). Secondly, an effective elastic tensor or the expression for free energy of the REV is established, and the thermodynamic forces associated with the internal variables are derived; finally, a suitable damage criterion is proposed and the evolution equation for the damage variables is determined [ 2 ]. According to Zhu et al. [ 40 ], for damaged solids with microcracks, the free energy W of a single unit can be expressed in the following general form: W = 1 2 ( ε − ε c ) : C m : ( ε − ε c ) + 1 2 ε c : C b : ε c (5) where W represents the elastic free energy of a solid matrix, and the term of the right side of the equation is the free energy stored in a solid matrix caused by inelastic strain related to crack, and C b is the fourth-order back stress modulus. We are taking into account the internal relationship [ 37 ] between the Mori–Tanaka homogenization method (MT) and linear elastic fracture mechanics in dealing with crack problems. In the case of isotropic and open cracks, the following effective elastic tensors can be obtained by the application of the MT method: C hom = 1 1 + α 1 ϕ 3 k m J + 1 1 + α 2 ϕ 2 µ m K (6) where α 1 , α 2 are constants only related to the Poisson coefficient ν m of the rock matrix, α 1 = 16 9 1 − ( ν m ) 2 1 − 2 ν m , and α 2 = 32 45 ( 1 − ν m )( 5 − ν m )) 2 − ν m . Considering the expression of C hom , the system free energy is expressed as W = ε : C hom ε /2. By combining Equation ( 5 ), we can obtain: C b = 1 α 1 ϕ 3 k m J + 1 α 2 ϕ 2 µ m K (7) It should be pointed out that according to the research of Zhu et al. [ 40 ], Equation ( 5 ) is also suitable for crack closure. In this case, the energy dissipation mechanism of crack propagation and sliding friction coupling exist in the REV. The analytical relationship between inelastic strain and macroscopic strain caused by microcracks is no longer valid.

Sustainability 2022 , 14 , 16296 6 of 17 The thermaldynamic force σ c associated with the inelastic strain ε c can be determined from the system free energy, that is, the local stress acting on the crack: σ c = − ∂ W ∂ ε c = σ − C b : ε c (8) Furthermore, σ c is decomposed into two portions: deviatoric stress s c and hydrostatic part: s c = K : σ c and p c = tr σ c /3 . In order to describe the inelastic strain due to the sliding friction of closed cracks and to capture the compressive shear damage pattern of quasi-brittle rock materials under compression, a Coulomb-type yielding criterion based on local stress is adopted in the present study: f s ( σ c ) = || s c || + α p c ≤ 0 (9) where α is the coefficient of friction of the cracked surface of the rock material. The local stress tensor can be decomposed into deviatoric and hydrostatic parts: s = K : σ = σ − C b : ε c : K = s − 1 α 2 ϕ 2 µ m T (10) p = tr σ /3 = 1 3 σ − C b : ε c : δ = p − 1 α 1 ϕ k m β (11) Combined with the above equations, Equation ( 9 ) can be rewritten as: f s = || s − 1 α 2 ϕ 2 µ m T || + α ( p − 1 α 1 ϕ k m β ) ≤ 0 (12) where s = K : σ and p = tr σ /3 Considering the expression of free energy and the theory of irreversible thermodynamics, the thermaldynamic force associated with the damage variable ϕ , namely the damage driving force F ϕ , can be derived by following equation: F ϕ = − ∂ W ∂ϕ = − 1 2 ε c : ∂ C b ∂ϕ : ε c (13) Substituting Equations ( 4 ), ( 5 ) and ( 7 ) into Equation ( 13 ), the explicit form of damage driving force can be obtained as follows: F ϕ = 1 2 α 1 ϕ 2 k m β 2 + 1 α 2 ϕ 2 µ m T : T (14) The damage evolution criterion based on the strain energy release rate can also be derived f ϕ ( F ϕ , ϕ ) = F ϕ − R ( ϕ ) ≤ 0 (15) where R ( ϕ ) is the current damage evolution resistance force 3.2. Coupled Friction–Damage Effect and the Strength Criterion As mentioned above, in the process of rock deterioration, there are two fundamental pathways for energy dissipation in compressive stress-dominated loading condition. One is damage evolution caused by the development of microcracks, and the other is friction caused by sliding along the fissure surface accompanied by dilatation (volume expansion) caused by the non-smooth crack surface. Damage and inelastic strain constantly rise upon loading in the damage–friction coupling process.

Sustainability 2022 , 14 , 16296 7 of 17 In this study, the evolution of damage variable ω and inelastic strain ε c are determined by the associated flow rule, and the directions of the evolution are determined by the orthogonalization criteria: ˙ ϕ = λ ϕ ∂ f ϕ ∂ F ϕ = λ ϕ (16) ˙ ε c = λ c ∂ f s ∂ σ c = λ c ( V + 1 3 α δ ) (17) where V = s c / || s c || (18) in which V represents the flow direction of the partial portion of inelastic strain ε c λ ϕ and λ c are damage and plasticity multipliers, respectively Comparing the expression of Equation ( 18 ) with that of Equation ( 4 ), one has ˙ T = λ c V , ˙ β = λ c α (19) Under standard triaxial loading circumstances, the evolution direction V of inelastic shear strain does not change, so the cumulative damage parameter is Λ ϕ = R ε λ ϕ and cumulative inelastic variable Λ c = R ε λ c . Thus, ε c = Λ c ( V + 1 3 α δ ) , ϕ = Λ ϕ (20) Rock material in engineering is mainly subjected to compression, and the strength needs to be determined within the damage–friction coupling framework. It is generally believed that inelastic strain causes the strengthening behavior of the material, while damage causes the strain softening after the peak stress. Therefore, two competing nonlinear mechanical mechanisms occur in the coupling process of inelastic strain and damage. It is difficult to determine the analytical form of material strength for plastic damage coupled models For the conventional triaxial compression loading path, in the principal stress space, the stress tensor is σ = σ 1 σ 2 σ 3 . At the same time, assuming that σ 1 < σ 2 = σ 3 , the deviatoric stress s is s = σ 1 − σ 3 3 2 − 1 − 1 (21) When the applied stress increases monotonously, the flow direction V can be written as V = s / || s || V = − 1 √ 6 2 − 1 − 1 (22) According to the relation s c = K : σ c and p c = tr σ c /3, we can obtain the following equations || s c || = − r 2 3 ( σ 1 − σ 3 ) − 2 µ m α 2 Λ c ϕ (23) p c = 1 3 ( σ 1 + 2 σ 3 ) − k m α α 1 Λ c ϕ (24) Substituting Equation ( 20 ) into Equations ( 13 ) and ( 15 ), the damage criteria in the case of crack closure can be obtained as follows: f ϕ = k m α 2 2 α 1 + µ m α 2 Λ c ϕ 2 − R ( ω ) ≤ 0 (25)

Sustainability 2022 , 14 , 16296 8 of 17 For simplicity, let us define χ = k m α 2 2 α 1 + µ m α 2 . Substituting Equations ( 23 ) and ( 24 ) into Equation ( 9 ), one obtains f s = − r 2 3 ( σ 1 − σ 3 ) + 1 3 α ( σ 1 + 2 σ 3 ) − 2 χ Λ c ϕ ≤ 0 (26) From f ϕ = 0, it follows that Λ c ϕ = s R ( ϕ ) χ (27) Substituting Equation ( 27 ) into Equation ( 26 ) and using the sign convention in geotechnics, it is derived: f s = σ 1 − 2 α + √ 6 √ 6 − α σ 3 − 6 p R ( ϕ ) χ √ 6 − α ≤ 0 (28) Noting that p = ( σ 1 + 2 σ 3 ) /3 and q = ( σ 1 − σ 3 ) , Equation ( 28 ) can be reformed into a loading function on the p - q surface as follows: f s = q − r 3 2 α p − q 6 R ( ϕ ) χ ≤ 0 (29) When the damage variable reaches its critical value ϕ = ϕ c , the damage resistance function R ( ϕ ) reaches its maximum value R ( ϕ c ) , at which point the axial stress reaches its maximum value and the material reaches its peak strength. Based on the above analysis, the rock strength envelope predicted by the coupled friction-damage model is expressed in the following form: f s = q − r 3 2 α p − q 6 R ( ϕ c ) χ = 0 (30) As mentioned above, geotechnical materials have hardening and softening properties So, for a given loading path, R ( ϕ ) reaches its maximum value R ( ϕ c ) when the damage factor reaches its critical value ϕ c , namely R ( ϕ ) < R ( ϕ c ) , and the rock is in the strainhardening phase. Meanwhile, for R ( ϕ ) > R ( ϕ c ) , the rock is in the strain-softening phase The mathematical description of the above properties is R 0 ( ϕ ) > 0 when 0 < ϕ < ϕ c ; R 0 ( ϕ ) = 0 when ϕ = ϕ c ; R 0 ( ϕ ) < 0 when ϕ < ϕ c . For this purpose, the dimensionless parameter ξ = ϕ / ϕ c is defined, and the following expression for the damage resistance R ( ϕ ) is used: R ( ϕ ) = R ( ϕ c ) 2 ξ 1 + ξ 2 (31) 4. Numerical Simulation of Damage–Friction Coupling Model and Validation by Experimental Results 4.1. Description of the Returning Mapping Procedure After the constitutive model and strength criterion are established, the inelastic strain ε c and damage variable ω are calculated iteratively according to the loading criterion The values of the k − 1 loading step variables ε k − 1 , T k − 1 , β k − 1 , ϕ k − 1 and σ k − 1 are known. The flow for calculating the kth loading step ε k , T k , β k , ϕ k , σ k using strain loading is shown in Figure 5 . 1 The strain increment d ε k was superimposed onto ε k − 1 to estimate the strain ε k at step k, and the macroscopic stress σ k was preliminarily calculated according to Equation ( 3 ). 2 The loading and unloading conditions are judged by the yield function f s ( σ c ) . If f s ( σ c ) > 0, then by the consistency condition ˙ f ( ε , β , T ) = 0 and Equation ( 19 ), we find λ c and hence T k , β k ; f s ( σ c ) ≤ 0; then, we update the stress σ k according to ε k only 3 Calculate the damage driving force F ϕ ( k ) according to the updated T k , β k , based on Equation ( 14 ).

Sustainability 2022 , 14 , 16296 9 of 17 4 Examine the damage yield function Equation ( 15 ). If f ϕ ( F ϕ , ϕ ) > 0, apply the Newton–Raphson algorithm to calculate ϕ k , and if f ϕ ( F ϕ , ϕ ) ≤ 0, then there is no damage increment 5 Update the stress σ k from Equation ( 3 ). 6 Substitute the updated ε k , T k , β k , ϕ k , σ k into the next loading loop Yes Yes No Figure 5. Flow chart of the numerical returning mapping procedure 4.2. Determination of Model Parameters The coupled damage–friction model proposed in this paper contains only five parameters, E m , ν m , α , R ( ϕ c ) and ϕ c , which can all be determined by a set of conventional triaxial mechanical tests. The Young’s modulus and Poisson’s coefficient are taken as E m = 30 GPa and ν m = 0.1 for the hard rock samples from the Pingdingshan coal mine. The remaining parameters were determined as follows: The parameter ϕ c included in the damage criterion is the damage threshold value, the value of which corresponds to the damage value at the peak stress. According to Lockner [ 47 ], the value of ϕ c is approximately linearly related to the perimeter pressure For simplicity, a constant value of ϕ c = 1.5 is taken here. The effect of the value of ϕ c on the numerical simulation results will be discussed through parametric sensitivity analysis As mentioned above, when ω = ω c that is, R ( ϕ c ) = R ( ϕ c ) , the material reaches its maximum axial stress (the intensity envelope on the p - q surface is shown in Figure 6 ). Based on the protected hard rock test data, the parameter values can be determined by applying Equation ( 30 ): α 0 = 1.3, R ( ϕ c ) = 5.34 × 10 − 2 MPa.

Sustainability 2022 , 14 , 16296 10 of 17 0 50 100 150 0 50 100 150 200 250 y = 1.8059 x + 40.8307 Figure 6. Strength envelope of the Pingdingshan hard rock under triaxial compressions 4.3. Numerical Simulations for Diorite in Pingdingshan Coal Mine Utilizing the values of the model parameters established in the preceding subsection, the stress–strain data from conventional triaxial compression tests are numerically simulated for different envelope pressure conditions (envelope pressure P c = 0, 10, 20 and 30 MPa), as shown in Figures 7 – 10 . The model more accurately simulates the main macromechanical behavior of the rock, and it better describes the strength and stress-softening characteristics of the rock under different circumferential pressure conditions, especially in the elastic phase. Both in the axial and lateral directions, the stress–strain relationship is more accurately described; after the material enters non-linear deformation, before the peak strength, the model also fits the lateral deformation of the rock relatively well -1.5 -1 -0.5 0 0.5 1 Strain (%) 20 40 60 80 100 120 140 Deviatoric stress (MPa) Confining pressure p c = 0 MPa Figure 7. Simulation of conventional triaxial compression test on hard rock protection layer rock in the Pingdingshan coal mine with the confining pressure of 0 MPa.

Sustainability 2022 , 14 , 16296 11 of 17 -1.5 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 Deviatoric stress (MPa) Confining pressure p c = 10 MPa Figure 8. Simulation of conventional triaxial compression test on hard rock protection layer rock in the Pingdingshan coal mine with the confining pressure of 10 MPa -1.5 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 Confining pressure p c = 20 MPa Figure 9. Simulation of conventional triaxial compression test on hard rock protection layer rock in the Pingdingshan coal mine with the confining pressure of 20 MPa According to the above simulated results, it can be seen that the constructed damage–friction coupling model in this paper can capture the main properties of the stress–strain curves of Pingdingshan hard rock under different confining pressures (envelope pressure P c = 0, 10, 20 and 30 MPa). Especially for confining pressures of 0 MPa, 10 MPa, and 20 MPa, this model can better simulate the peak pressure and related damage. However, when the confining pressure is 30 MPa, the simulation result of the height strength is smaller than the experimental data. This is because the rock sample in this experiment is obtained by the blasting method, so some microcracks exist inside the body, resulting in a certain dispersion of its peak strength. It must be pointed out that small but noticeable differences on the lateral strain are observed in the above figures. This is due to the fact that the tested tock is collected close to the coal seam, so there is more coal inclusions distributed in the rock solids, resulting in a large discrete property of the rock samples. Since the rock samples in the Pingdingshan Coal mine are obtained by an explosive method, unnatural cracks appear in the rock samples, which may lead to irregular lateral deformation. For most geomatrials, non-associated plastic flow rules are usually adopted to study the lateral and volumetric strains. However, in the present study, for the

Sustainability 2022 , 14 , 16296 12 of 17 sake of simplicity, an associated flow law is adopted, which may result in a large difference between the experimental and analytical results of lateral strain -1.5 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 250 Confining pressure p c = 30 MPa Figure 10. Simulation of conventional triaxial compression test on hard rock protection layer rock in the Pingdingshan coal mine with the confining pressure of 30 MPa 4.4. Sensitivity Analysis of Parameter ϕ c In order to determine the effect of the parameter ϕ c on the numerical simulation results, a conventional test with an envelope pressure of 10 MPa is used as the subject of the study. Figure 11 compares the results of the simulations for ϕ c values of 1.5, 2.5 and 3.5, respectively. It can be seen that the higher the ϕ c is, the greater the axial and lateral strains at the peak stress point are and the smoother the hardening and softening curve, namely, the slower the rate of hardening and softening, due to the larger the critical value of the damage. In addition, the number of microcracks within the characteristic cell or their radius is greater when the critical value is reached, which is the same for the material deterioration and the axial and lateral strains. The strengthening speed becomes smaller when the inelastic strain between the elastic stage and the peak stress point of REV become larger -1 -0.5 0 0.5 1 1.5 50 100 150 Figure 11. Sensitivity analysis for parameter ϕ c As mentioned before, the magnitude of the parameter ϕ c is approximately linearly related to the surrounding pressure. For the sake of simplicity, a constant value was taken in the numerical simulations, which need further improvement of the simulation results for different envelope pressures. At the same time, due to the associated flow rule,

Sustainability 2022 , 14 , 16296 13 of 17 the magnitude of lateral strain is measured when the stress peak point is reached at the stress-softening rate. The smaller ϕ c is, the faster the lateral strain reaches the peak, and the faster the stress-softening rate 4.5. Numerical Simulations for Dagangshan Diabase For further validation of the proposed model, comparisons between experimental results of the Dagangshan diabase and the numerical simulations are carried out in this subsection. The experimental data are provided by an earlier work [ 48 ]. Now, we will provide numerical simulations of standard triaxial compression tests conducted with confining pressures of p c = 10, 15, 20, 30, 50 MPa. The specific model parameters used are listed in Table 1 below. The results are illustrated in Figures 12 – 16 with respect to confining pressures. It can be seen that the predicted numerical solutions for both low and high confining pressure levels match quite well with the experimental data Table 1. Identified parameters in the numerical simulations for Dagangshan diabase E m (MPa) ν m α 0 ϕ c R ( ϕ c ) 55,000 0.2 1.65 5.5 1.3 × 10 − 2 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 Deviatoric stress (MPa) Confining pressure p c = 10 MPa Figure 12. Simulation of conventional triaxial compression test on Dagangshan diabase with the confining pressure of 10 MPa -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 250 Deviatoric stress (MPa) Confining pressure p c = 15 MPa Figure 13. Simulation of conventional triaxial compression test on Dagangshan diabase with the confining pressure of 15 MPa.

Sustainability 2022 , 14 , 16296 14 of 17 -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 250 300 Deviatoric stress (MPa) Confining pressure p c = 20 MPa Figure 14. Simulation of conventional triaxial compression test on Dagangshan diabase with the confining pressure of 20 MPa -1 -0.5 0 0.5 1 Strain (%) 50 100 150 200 250 300 350 Deviatoric stress (MPa) Confining pressure p c = 30 MPa Figure 15. Simulation of conventional triaxial compression test on Dagangshan diabase with the confining pressure of 30 MPa -1.5 -1 -0.5 0 0.5 1 Strain (%) 100 200 300 400 500 Deviatoric stress (MPa) Confining pressure p c = 50 MPa Figure 16. Simulation of conventional triaxial compression test on Dagangshan diabase with the confining pressure 50 MPa.

Sustainability 2022 , 14 , 16296 15 of 17 5. Conclusions Based on homogenization and irreversible thermodynamics theories, a coupled damage–friction model for Pingdingshan hard rock in a protective coal seam and Dagangshan diabase are constructed by the use of an associated flow law. The model successfully simulates the main mechanical behaviors of the tested rock in a conventional triaxial compression test. In order to assess the proposed coupled plastic–damage model, numerical simulations of triaxial compression tests on Pingdingshan rock and Dagangshan diabase have been carried out. By comparing with the experimental data, the proposed model can well describe the peak strengths and the mechanical responses from low to high levels of confining pressure for two types of rocks. However, the lateral deformation around the peak stress is not accurately predicted by the model due to the associated plastic flow rule. Determination of the critical damage parameter is also crucial, which determines the peak strength Determination of the free energy expressions for the RVE of cracked rock based on homogenization theory is the key point of the research results, which are combined with thermodynamic theory to the closed damage-friction coupling and strain strengthening and weakening processes. The proposed multiscale intrinsic model has the advantage of having only five parameters, which have specific physical meanings, and can be easily determined by the experiments. Inelastic deformation and the propagation of cracks are two of the fundamental mechanisms of material damage and destruction. Therefore, the coupled damage–friction model is more suitable for describing the mechanical behavior of such quasi-brittle material. In geotechnics, non-associated flow law is generally considered to be necessary. This model is suitable for quasi-brittle rock and closed fracture. Secondly, the selection of damage parameters in this model is related to the fracture density. After rock is compressed, cracks in its interior will be derived and developed, which is consistent with the material failure mechanism described above. However, in this study, the associated flow law of inelastic strain is adopted to predict the deformation of the tested rock Author Contributions: Software, C.S.; Validation, S.H.; Investigation, S.H.; Data curation, W.L.; Writing—original draft, C.Z. and W.L.; Writing—review & editing, J.Z.; Visualization, C.Z. and C.S.; Supervision, J.Z. All authors have read and agreed to the published version of the manuscript Funding: This research was funded by the project “Research of high efficiency precracking technology for working face in rock protective coal seam of kilometer-deep mine” Institutional Review Board Statement: Not applicable Informed Consent Statement: Not applicable Data Availability Statement: The experimental data used to support the findings of this study are available from the corresponding author upon request Acknowledgments: The authors would like to thank the project "Research of high efficiency precracking technology for working face in rock protective coal seam of kilometer-deep mine" for the funding support Conflicts of Interest: The writers affirm that they do not have any conflict of interest, either financial or otherwise References 1 Zhang, D.; Sun, Z.; Fang, Q. Scientific problems and research proposals for Sichuan–Tibet railway tunnel construction Undergr. Space 2022 , 7 , 419–439. [ CrossRef ] 2 Zhu, Q.; Liu, H.; Wang, W.; Shao, J. A micromechanical constitutive damage model for beishan granite Chin. J. Rock Mech. Eng 2015 , 34 , 433–439 3 Kim, M.M.; Ko, H.Y. Multistage triaxial testing of rocks Geotech. Test. J 1979 , 2 , 98–105 4 Xie, H.; He, C. Study of the unloading characteristics of a rock mass using the triaxial test and damage mechanics Int. J. Rock Mech. Min. Sci 2004 , 41 , 74–80. [ CrossRef ] 5 Liu, A.; Liu, S. A fully-coupled water-vapor flow and rock deformation/damage model for shale and coal: Its application for mine stability evaluation Int. J. Rock Mech. Min. Sci 2021 , 146 , 104880. [ CrossRef ]

Sustainability 2022 , 14 , 16296 16 of 17 6 Wei, X.; Sun, Q. Damage mechanism and constitutive model of surface paste disposal of a coal mine collapsed pit in a complicated environment Constr. Build. Mater 2021 , 304 , 124637. [ CrossRef ] 7 Chang, C.; Haimson, B. A failure criterion for rocks based on true triaxial testing. In The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007–2014 ; Springer: Berlin/Heidelberg, Germany, 2012; pp. 259–262 8 Höfer, K.; Thoma, K. Triaxial tests on salt rocks. In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts; Elsevier: Amsterdam, The Netherlands, 1968; Volume 5, pp. 195–196 9 Lane, K.; Heck, W. Triaxial testing for strength of rock joints. In Proceedings of the 6 th US Symposium on Rock Mechanics (USRMS). OnePetro, Rolla, MI, USA, 28–30 October 1964 10 Dragon, A.; Mroz, Z. A continuum model for plastic-brittle behaviour of rock and concrete Int. J. Eng. Sci 1979 , 17 , 121–137 [ CrossRef ] 11 Frantziskonis, G.; Desai, C. Elastoplastic model with damage for strain softening geomaterials Acta Mech 1987 , 68 , 151–170 [ CrossRef ] 12 Maugin, G.A The Thermomechanics of Plasticity and Fracture ; Cambridge University Press: Cambridge, UK, 1992; Volume 7 13 Hansen, N.; Schreyer, H. A thermodynamically consistent framework for theories of elastoplasticity coupled with damage Int. J Solids Struct 1994 , 31 , 359–389. [ CrossRef ] 14 Lubarda, V.; Krajcinovic, D. Some fundamental issues in rate theory of damage-elastoplasticity Int. J. Plast 1995 , 11 , 763–797 [ CrossRef ] 15 Lawn, B.R.; Marshall, D.B. Nonlinear stress-strain curves for solids containing closed cracks with friction J. Mech. Phys. Solids 1998 , 46 , 85–113. [ CrossRef ] 16 de Sciarra, F.M. Hardening plasticity with nonlocal strain damage Int. J. Plast 2012 , 34 , 114–138. [ CrossRef ] 17 Paliwal, B.; Ramesh, K. An interacting micro-crack damage model for failure of brittle materials under compression J. Mech Phys. Solids 2008 , 56 , 896–923. [ CrossRef ] 18 Yazdani, S.; Schreyer, H. Combined plasticity and damage mechanics model for plain concrete J. Eng. Mech 1990 , 116 , 1435–1450 [ CrossRef ] 19 Abu-Lebdeh, T.M.; Voyiadjis, G.Z. Plasticity-damage model for concrete under cyclic multiaxial loading J. Eng. Mech 1993 , 119 , 1465–1484. [ CrossRef ] 20 Luccioni, B.M.; Rougier, V.C. A plastic damage approach for confined concrete Comput. Struct 2005 , 83 , 2238–2256. [ CrossRef ] 21 Jason, L.; Huerta, A.; Pijaudier-Cabot, G.; Ghavamian, S. An elastic plastic damage formulation for concrete: Application to elementary tests and comparison with an isotropic damage model Comput. Methods Appl. Mech. Eng 2006 , 195 , 7077–7092 [ CrossRef ] 22 Grassl, P.; Jirásek, M. Plastic model with non-local damage applied to concrete Int. J. Numer. Anal. Methods Geomech 2006 , 30 , 71–90. [ CrossRef ] 23 Wu, J.Y.; Li, J.; Faria, R. An energy release rate-based plastic-damage model for concrete Int. J. Solids Struct 2006 , 43 , 583–612 [ CrossRef ] 24 Khan, A.S.; Xiang, Y.; Huang, S. Behavior of Berea sandstone under confining pressure part I: Yield and failure surfaces, and nonlinear elastic response Int. J. Plast 1991 , 7 , 607–624. [ CrossRef ] 25 Chiarelli, A.S.; Shao, J.F.; Hoteit, N. Modeling of elastoplastic damage behavior of a claystone Int. J. Plast 2003 , 19 , 23–45 [ CrossRef ] 26 Salari, M.; Saeb, S.a.; Willam, K.; Patchet, S.; Carrasco, R. A coupled elastoplastic damage model for geomaterials Comput. Methods Appl. Mech. Eng 2004 , 193 , 2625–2643. [ CrossRef ] 27 Conil, N.; Djeran-Maigre, I.; Cabrillac, R.; Su, K. Thermodynamics modelling of plasticity and damage of argillite Comptes Rendus Mec 2004 , 332 , 841–848. [ CrossRef ] 28 Shao, J.F.; Jia, Y.; Kondo, D.; Chiarelli, A.S. A coupled elastoplastic damage model for semi-brittle materials and extension to unsaturated conditions Mech. Mater 2006 , 38 , 218–232. [ CrossRef ] 29 Chen, L.; Shao, J.; Zhu, Q.Z.; Duveau, G. Induced anisotropic damage and plasticity in initially anisotropic sedimentary rocks Int. J. Rock Mech. Min. Sci 2012 , 51 , 13–23. [ CrossRef ] 30 Lai, Y.; Jin, L.; Chang, X. Yield criterion and elasto-plastic damage constitutive model for frozen sandy soil Int. J. Plast 2009 , 25 , 1177–1205. [ CrossRef ] 31 Parisio, F.; Samat, S.; Laloui, L. Constitutive analysis of shale: A coupled damage plasticity approach Int. J. Solids Struct 2015 , 75 , 88–98. [ CrossRef ] 32 Klisi ´nski, M.; Mroz, Z. Description of inelastic deformation and degradation of concrete Int. J. Solids Struct 1988 , 24 , 391–416 [ CrossRef ] 33 Lubliner, J.; Oliver, J.; Oller, S.; Oñate, E. A plastic-damage model for concrete Int. J. Solids Struct 1989 , 25 , 299–326. [ CrossRef ] 34 Zhu, Q.; Shao, J.F.; Mainguy, M. A micromechanics-based elastoplastic damage model for granular materials at low confining pressure Int. J. Plast 2010 , 26 , 586–602. [ CrossRef ] 35 Gambarotta, L.; Lagomarsino, S. A microcrack damage model for brittle materials Int. J. Solids Struct 1993 , 30 , 177–198 [ CrossRef ] 36 Prat, P.C.; Bažant, Z.P. Tangential stiffness of elastic materials with systems of growing or closing cracks J. Mech. Phys. Solids 1997 , 45 , 611–636. [ CrossRef ]

Sustainability 2022 , 14 , 16296 17 of 17 37 Pensée, V.; Kondo, D.; Dormieux, L. Micromechanical analysis of anisotropic damage in brittle materials J. Eng. Mech 2002 , 128 , 889–897. [ CrossRef ] 38 François, B.; Dascalu, C. A two-scale time-dependent damage model based on non-planar growth of micro-cracks J. Mech Phys. Solids 2010 , 58 , 1928–1946. [ CrossRef ] 39 Zhu, Q.; Shao, J.F. A refined micromechanical damage–friction model with strength prediction for rock-like materials under compression Int. J. Solids Struct 2015 , 60 , 75–83. [ CrossRef ] 40 Zhu, Q.Z.; Kondo, D.; Shao, J. Micromechanical analysis of coupling between anisotropic damage and friction in quasi brittle materials: Role of the homogenization scheme Int. J. Solids Struct 2008 , 45 , 1385–1405. [ CrossRef ] 41 Yuan, S.S.; Zhu, Q.Z.; Zhao, L.Y.; Chen, L.; Shao, J.F.; Zhang, J. Micromechanical modelling of short-and long-term behavior of saturated quasi-brittle rocks Mech. Mater 2020 , 142 , 103298. [ CrossRef ] 42 Hu, K.; Zhu, Q.; Chen, L.; Shao, J.; Liu, J. A micromechanics-based elastoplastic damage model for rocks with a brittle–ductile transition in mechanical response Rock Mech. Rock Eng 2018 , 51 , 1729–1737. [ CrossRef ] 43 Cicekli, U.; Voyiadjis, G.Z.; Al-Rub, R.K.A. A plasticity and anisotropic damage model for plain concrete Int. J. Plast 2007 , 23 , 1874–1900. [ CrossRef ] 44 Voyiadjis, G.Z.; Taqieddin, Z.N.; Kattan, P.I. Anisotropic damage–plasticity model for concrete Int. J. Plast 2008 , 24 , 1946–1965 [ CrossRef ] 45 Zhu, Q.; Kondo, D.; Shao, J.; Pensee, V. Micromechanical modelling of anisotropic damage in brittle rocks and application Int. J Rock Mech. Min. Sci 2008 , 45 , 467–477. [ CrossRef ] 46 Budiansky, B.; O’connell, R.J. Elastic moduli of a cracked solid Int. J. Solids Struct 1976 , 12 , 81–97. [ CrossRef ] 47 Lockner, D.A. A generalized law for brittle deformation of Westerly granite J. Geophys. Res. Solid Earth 1998 , 103 , 5107–5123 [ CrossRef ] 48 Chen, F. Rheological Strength Aging Model of Hard Brittle Fracture Rock in High Dam Zone and Its Engineering Application. Ph.D. Thesis, Shandong University, Jinan, China, 2009.