Experimental Investigation and Micromechanical Modeling of Elastoplastic Damage Behavior of Sandstone

The mechanical behavior of the sandstone at the dam site is important to the stability of the hydropower station to be built in Southwest China. A series of triaxial compression tests under different confining pressures were conducted in the laboratory. The critical stresses were determined and the relationship between the critical stress and confining pressure were analyzed. The Young’s modulus increases non-linearly with the confining pressure while the plastic strain increment Nϕ and the dilation angle ϕ showed a negative response. Scanning electron microscope (SEM) tests showed that the failure of the sandstone under compression is a coupled process of crack growth and frictional sliding. Based on the experimental results, a coupled elastoplastic damage model was proposed within the irreversible thermodynamic framework. The plastic deformation and damage evolution were described by using the micromechanical homogenization method. The plastic flow is inherently driven by the damage evolution. Furthermore, a numerical integration algorithm was developed to simulate the coupled elastoplastic damage behavior of sandstone. The main inelastic properties of the sandstone were well captured. The model will be implemented into the finite element method (FEM) to estimate the excavation damaged zones (EDZs) which can provide a reference for the design and construction of such a huge hydropower project.


Introduction
Thanks to their good geological and mechanical properties, sandstones serve as the privileged candidate materials for many rock engineering applications, such as hydropower engineering, petroleum engineering, road engineering, and other engineering applications [1]. Some examples of projects involving excavation in sandstones include the Xiang Jiaba hydropower station in Southwest China and railway tunnels in Northwest China. In this context, attaining a deeper understanding of the mechanical properties of sandstone is of crucial importance to the design and construction of rock engineering projects in such host rock.
an up-scaling method [23][24][25][26][27]. These models contain generally far fewer parameters compared to phenomenological ones. However, most of the micromechanical models were not formulated within the framework of irreversible thermodynamics. For practical rock mechanics problems such as the dam site sandstone, the predictive capability of these micromechanical based models is still under discussion [28].
In this study, the mechanical properties of sandstone collected from a dam foundation were systematically investigated. The basic inelastic mechanical behavior of the sandstone will be outlined. Based on the experimental results, a comprehensive analysis was performed on building a coupled elastoplastic damage model within the framework of irreversible thermodynamics. The damage evolution and plastic flow rules were developed according to the micromechanical based homogenization method. A new computational integration algorithm was proposed to deal with the coupled elastoplastic damage model. After the identification of the model's parameters, the proposed model was applied to simulate the experimental results of the sandstone under uniaxial/triaxial compression conditions. The proposed model would help with the design and construction of a huge hydropower project using sandstone.

Description of Sandstone
The sandstone used throughout this study was collected from the dam site of a hydropower project in Southwest China. The X-ray diffraction (XRD) and optical microscopy (OM) tests suggested that the mineralogical composition of this sandstone is about quartz (55%), feldspar (25%), sandy and clay detritus (20%). The quartz is of monocrystals of about 0.1-1.0 mm in size while the feldspar has granular structural and large phenocryst crystals. The mean density and initial porosity are 2.23 g/cm 3 and 8.43%, respectively. A thin section of sandstone was prepared for the SEM test and the result is given in Figure 1. From the 2000× enlarged image, we can see that the sandstone is a blocky structure. The grains of sandstone are tightly cemented together. Many pores and microcracks are uniformly distributed. Besides, most of the sizes of the initial defects are less than 5 µm. dam site sandstone, the predictive capability of these micromechanical based models is still under discussion [28].
In this study, the mechanical properties of sandstone collected from a dam foundation were systematically investigated. The basic inelastic mechanical behavior of the sandstone will be outlined. Based on the experimental results, a comprehensive analysis was performed on building a coupled elastoplastic damage model within the framework of irreversible thermodynamics. The damage evolution and plastic flow rules were developed according to the micromechanical based homogenization method. A new computational integration algorithm was proposed to deal with the coupled elastoplastic damage model. After the identification of the model's parameters, the proposed model was applied to simulate the experimental results of the sandstone under uniaxial/triaxial compression conditions. The proposed model would help with the design and construction of a huge hydropower project using sandstone.

Description of Sandstone
The sandstone used throughout this study was collected from the dam site of a hydropower project in Southwest China. The X-ray diffraction (XRD) and optical microscopy (OM) tests suggested that the mineralogical composition of this sandstone is about quartz (55%), feldspar (25%), sandy and clay detritus (20%). The quartz is of monocrystals of about 0.1-1.0 mm in size while the feldspar has granular structural and large phenocryst crystals. The mean density and initial porosity are 2.23 g/cm 3 and 8.43%, respectively. A thin section of sandstone was prepared for the SEM test and the result is given in Figure 1. From the 2000× enlarged image, we can see that the sandstone is a blocky structure. The grains of sandstone are tightly cemented together. Many pores and microcracks are uniformly distributed. Besides, most of the sizes of the initial defects are less than 5 μm.
Because of its significant impact on the safety of the dam foundation, the permeability evolution of this sandstone in the triaxial and hydrostatic compression has been studied sufficiently [29,30]. The change of permeability with crack growth under different pore pressures was studied. The inert gas test technique was developed to measure the permeability of this sandstone. This study is devoted to the experimental and numerical investigations on the mechanical properties of this sandstone. All the cylindrical specimens were drilled and polished from the same block of material to a diameter of 50 mm and approximately 100 mm in length. The tolerance of straightness and flatness of the samples meets the requirement of the International Society for Rock Mechanics and Rock Engineering (ISRM) suggested method [31].  Because of its significant impact on the safety of the dam foundation, the permeability evolution of this sandstone in the triaxial and hydrostatic compression has been studied sufficiently [29,30]. The change of permeability with crack growth under different pore pressures was studied. The inert gas test technique was developed to measure the permeability of this sandstone. This study is devoted to the experimental and numerical investigations on the mechanical properties of this sandstone. All the cylindrical specimens were drilled and polished from the same block of material to a diameter

Experimental Method
A servo-controlled rock mechanics experimental system was used to complete all the experiments ( Figure 2). This apparatus comprises a triaxial cell, three high-pressure servo-controlled pumps, and a data monitoring system. The confining pressure and pore pressure up to 60 MPa are loaded through separated pumps. The maximum deviatoric stress is 375 MPa. The axial strain ε 1 is measured by two linear variable displacement transducers (LVDT) with a resolution of ±1 µm, while the radial deformation ε 3 is monitored continuously using a ring radial displacement transducer chain wrapped tightly around the middle height of the specimen. All the stress and strain data are monitored and recorded continuously by an integrated acquisition system.

Experimental Method
A servo-controlled rock mechanics experimental system was used to complete all the experiments ( Figure 2). This apparatus comprises a triaxial cell, three high-pressure servo-controlled pumps, and a data monitoring system. The confining pressure and pore pressure up to 60 MPa are loaded through separated pumps. The maximum deviatoric stress is 375 MPa. The axial strain 1 is measured by two linear variable displacement transducers (LVDT) with a resolution of ±1 μm, while the radial deformation 3 is monitored continuously using a ring radial displacement transducer chain wrapped tightly around the middle height of the specimen. All the stress and strain data are monitored and recorded continuously by an integrated acquisition system.
Following the ISRM suggested method, the conventional triaxial compression tests are performed in two steps. The confining pressure is loaded to the desired value. Then, the deviatoric stress ( 1 − 3 ) is increased in the displacement control method until specimen failure while the confining pressure is kept constant during the whole process.  Following the ISRM suggested method, the conventional triaxial compression tests are performed in two steps. The confining pressure is loaded to the desired value. Then, the deviatoric stress (σ 1 − σ 3 ) is increased in the displacement control method until specimen failure while the confining pressure is kept constant during the whole process.

Experimental Results and Discussions
The conventional uniaxial/triaxial compression test results are presented in Figure 3. Consistently, the stress-strain curves of sandstone under different confining pressures can be divided into four symptomatic stages. (1) The curves are concave upward, a consequence of the closure of pre-existing defects. With the increase of confining pressure, the initial hardening is usually not marked because the pre-existing defects have already closed before the additional stress is applied. (2) The curves are approximately linear, as the elastic deformation of the grains is in the dominant role. (3) The curves depart from perfectly elastic behavior as a large number of microcracks proliferate and propagate throughout the specimens. (4) The failure of the specimens marks the significant drop of the curves due to the macroscopic fractures developed by the lining-up of microcracks. Furthermore, the frictional sliding capacity of fractures finally sustains the relative flattening of the curves.

Experimental Results and Discussions
The conventional uniaxial/triaxial compression test results are presented in Figure 3. Consistently, the stress-strain curves of sandstone under different confining pressures can be divided into four symptomatic stages. (1) The curves are concave upward, a consequence of the closure of pre-existing defects. With the increase of confining pressure, the initial hardening is usually not marked because the pre-existing defects have already closed before the additional stress is applied.  Based on the recorded axial strain 1 and radial strain 3 , the volumetric strain can be calculated according to the relation: = 1 + 2 3 . The evolution of the volumetric strain during a typical experimental is presented in Figure 4. Additionally, a virtual elastic reference line is added based on the extrapolation of the linear part of the volumetric strain. According to the figure, the volumetric strain switches from the compaction-dominated phase to the dilatancy-dominated phase. Three typical points map the evolution process of the volumetric strain. The first is the onset of dilatancy ′ , which can be determined at the point where the volumetric strain departs from the approximately linear part [32]. This implies that at stress beyond the ′ the deviatoric stress induces the pore structure to dilate. The point CD is a reversal point where the volumetric strain changes from compaction to dilatancy. Point D represents the volume of the specimen returning to the initial value. Based on the recorded axial strain ε 1 and radial strain ε 3 , the volumetric strain can be calculated according to the relation: ε v = ε 1 + 2ε 3 . The evolution of the volumetric strain during a typical experimental is presented in Figure 4. Additionally, a virtual elastic reference line is added based on the extrapolation of the linear part of the volumetric strain. According to the figure, the volumetric strain switches from the compaction-dominated phase to the dilatancy-dominated phase. Three typical points map the evolution process of the volumetric strain. The first is the onset of dilatancy C , which can be determined at the point where the volumetric strain departs from the approximately linear part [32]. This implies that at stress beyond the C the deviatoric stress induces the pore structure to dilate. The point CD is a reversal point where the volumetric strain changes from compaction to dilatancy. Point D represents the volume of the specimen returning to the initial value.  The critical stresses at the three points and the corresponding peak stresses under different confining pressures are listed in Table 1. The results show a strong positive influence of confining pressure on the stresses. The influence of 3 on the corresponding stresses is given in Figure 5. The linear Mohr-Coulomb criterion and the nonlinear Hoek-Brown criterion are adapted to fit the experimental results: where c is the cohesion and is the frictional angle, UCS is the uniaxial compressive stress, m, s and a are all material constants. For the intact rock, the parameters s and a are fixed to 1.0 and 0.5, respectively. The fitting results are plotted in Figure 5 and the obtained strength parameters are presented in Table 2.  The critical stresses at the three points and the corresponding peak stresses under different confining pressures are listed in Table 1. The results show a strong positive influence of confining pressure on the stresses. The influence of σ 3 on the corresponding stresses is given in Figure 5. The linear Mohr-Coulomb criterion and the nonlinear Hoek-Brown criterion are adapted to fit the experimental results: where c is the cohesion and ϕ is the frictional angle, UCS is the uniaxial compressive stress, m, s and a are all material constants. For the intact rock, the parameters s and a are fixed to 1.0 and 0.5, respectively. The fitting results are plotted in Figure 5 and the obtained strength parameters are presented in Table 2. Table 1. Critical stresses during the dilation and the peak stress of sandstone. In general, the Mohr-Coulomb provides a better representation of the experimental data. The determined cohesion c and frictional angle all increase during the hardening process. Even though the fitting R 2 s are all high than 0.97, the determined USC at the peak point (25.295 MPa) is much less than the experimental data. The is identified using the local gradient of a third-order polynomial fitted to the axial stressstrain curve and is given by Similarly, the Poisson's ratio is calculated by = − ⁄ [33]. Figure 6 illustrates the influence of the confining pressure on the deformation parameters ( and ) of the sandstone. The Young's modulus increases non-linearly with the confining pressure. When the confining pressure is lower, the increases relatively fast. The evolution of the Poisson's ratio within the test confining pressure is not clear. The mean value of the Poisson's ratio is 0.198. As determined above, dilatancy, which is characterized by the dilation angle , is a significant property of the sandstone. To assess the dilation angle from uniaxial or triaxial tests, Vermeer and De Borst [34] proposed the equation In general, the Mohr-Coulomb provides a better representation of the experimental data. The determined cohesion c and frictional angle ϕ all increase during the hardening process. Even though the fitting R 2 s are all high than 0.97, the determined USC at the peak point (25.295 MPa) is much less than the experimental data.
The E s is identified using the local gradient of a third-order polynomial fitted to the axial stress-strain curve and is given by E s = ∂(σ 1 − σ 3 )/∂ε 1 Similarly, the Poisson's ratio is calculated by υ = −∂ε 3 /∂ε 1 [33]. Figure 6 illustrates the influence of the confining pressure on the deformation parameters (E s and υ) of the sandstone. The Young's modulus increases non-linearly with the confining pressure. When the confining pressure is lower, the E s increases relatively fast. The evolution of the Poisson's ratio within the test confining pressure is not clear. The mean value of the Poisson's ratio is 0.198. In general, the Mohr-Coulomb provides a better representation of the experimental data. The determined cohesion c and frictional angle all increase during the hardening process. Even though the fitting R 2 s are all high than 0.97, the determined USC at the peak point (25.295 MPa) is much less than the experimental data. The is identified using the local gradient of a third-order polynomial fitted to the axial stressstrain curve and is given by Similarly, the Poisson's ratio is calculated by = − ⁄ [33]. Figure 6 illustrates the influence of the confining pressure on the deformation parameters ( and ) of the sandstone. The Young's modulus increases non-linearly with the confining pressure. When the confining pressure is lower, the increases relatively fast. The evolution of the Poisson's ratio within the test confining pressure is not clear. The mean value of the Poisson's ratio is 0.198. As determined above, dilatancy, which is characterized by the dilation angle , is a significant property of the sandstone. To assess the dilation angle from uniaxial or triaxial tests, Vermeer and De Borst [34] proposed the equation As determined above, dilatancy, which is characterized by the dilation angle φ, is a significant property of the sandstone. To assess the dilation angle from uniaxial or triaxial tests, Vermeer and De Borst [34] proposed the equation Materials 2020, 13, 3414 8 of 22 In Equation (3), dε p 1 and dε p v are axial and volumetric plastic strain increments, respectively. For φ > 0, the irreversible radial strain increment is larger than that of the axial strain, which indicates that the plastic volumetric strain increases. While φ < 0, the plastic contraction occurs [35]. The dilation angle is approved to relate to the plastic strain and the plastic strain increment N φ is always introduced The relationship between dilation parameters (N φ , φ) and confining pressure is presented in Figure 7. Both the plastic strain increment N φ and the dilation angle φ show a negative response to the confining pressure. Also, the drop of N φ is more significant under the lower σ 3 . An exponential function is introduced to approach the correlation between N φ and σ 3 . The correlation coefficient R 2 is 0.999.
Materials 2020, 13, x FOR PEER REVIEW 8 of 24 In Equation (3), and are axial and volumetric plastic strain increments, respectively. For > 0, the irreversible radial strain increment is larger than that of the axial strain, which indicates that the plastic volumetric strain increases. While < 0, the plastic contraction occurs [35]. The dilation angle is approved to relate to the plastic strain and the plastic strain increment is always introduced ( ) The relationship between dilation parameters ( , ϕ) and confining pressure is presented in Figure 7. Both the plastic strain increment and the dilation angle ϕ show a negative response to the confining pressure. Also, the drop of is more significant under the lower 3 . An exponential function is introduced to approach the correlation between and 3 . The correlation coefficient 2 is 0.999. The stress-strain curves during the cyclic loading tests under the confining pressure of 10 MPa is shown in Figure 8 together with the data from the monotonic compression test. Three key observations arising from the similarity of the envelope curves of the tests during different loading histories. Firstly, the peak stress and residual stress of the sandstone are independent of how the specimen is loaded. Secondly, the axial and radial envelope curves during the cyclic loading coincide well with the monotonic loading before peak stress. Thirdly, the radial strain response in the postpeak region is quite larger than that of the monotonic loading due to the strain localization and the location and orientation of stress-induced fractures [36]. The stress-strain curves during the cyclic loading tests under the confining pressure of 10 MPa is shown in Figure 8 together with the data from the monotonic compression test. Three key observations arising from the similarity of the envelope curves of the tests during different loading histories. Firstly, the peak stress and residual stress of the sandstone are independent of how the specimen is loaded. Secondly, the axial and radial envelope curves during the cyclic loading coincide well with the monotonic loading before peak stress. Thirdly, the radial strain response in the post-peak region is quite larger than that of the monotonic loading due to the strain localization and the location and orientation of stress-induced fractures [36].
The SEM tests were investigated on the thin sections prepared from the fractured surfaces of the failed cylindrical specimen. With the test results shown in Figure 9, the micromechanics of damage is investigated. Microcracks with a width of approximately 1-2.5 µm are observed. Enlarging the microcracks to 5000×, we can see that the surfaces of the microcracks are relatively smooth while apparent dislocation can reach 5 µm. This irreversible deformation is originated from the growth of the cracks. The SEM tests were investigated on the thin sections prepared from the fractured surfaces of the failed cylindrical specimen. With the test results shown in Figure 9, the micromechanics of damage is investigated. Microcracks with a width of approximately 1-2.5 μm are observed. Enlarging the microcracks to 5000×, we can see that the surfaces of the microcracks are relatively smooth while apparent dislocation can reach 5 μm. This irreversible deformation is originated from the growth of the cracks.

Thermodynamic Framework
Based on the mechanical tests and microscopic observation, it has been well confirmed that the failure of brittle sandstone can be attributed to the coupling between the irreversible deformation and damage induced by microcracks. The basic physical process of damage is the initiation, propagation, and coalescence of microcracks. In the compression stress state, the frictional sliding along the rough surfaces of the microcracks produces the irreversible deformation. In this context, a coupled elastoplastic damage model is more adequate to reproduce the inherent coupling between the plastic and damage dissipation processes. For the brittle rock materials, the assumption of small strains is appropriate. In the isothermal conditions, the total strain tensor ε can be decomposed into an elastic strain and a plastic strain according to the plastic theory  The SEM tests were investigated on the thin sections prepared from the fractured surfaces of the failed cylindrical specimen. With the test results shown in Figure 9, the micromechanics of damage is investigated. Microcracks with a width of approximately 1-2.5 μm are observed. Enlarging the microcracks to 5000×, we can see that the surfaces of the microcracks are relatively smooth while apparent dislocation can reach 5 μm. This irreversible deformation is originated from the growth of the cracks.

Thermodynamic Framework
Based on the mechanical tests and microscopic observation, it has been well confirmed that the failure of brittle sandstone can be attributed to the coupling between the irreversible deformation and damage induced by microcracks. The basic physical process of damage is the initiation, propagation, and coalescence of microcracks. In the compression stress state, the frictional sliding along the rough surfaces of the microcracks produces the irreversible deformation. In this context, a coupled elastoplastic damage model is more adequate to reproduce the inherent coupling between the plastic and damage dissipation processes. For the brittle rock materials, the assumption of small strains is appropriate. In the isothermal conditions, the total strain tensor ε can be decomposed into an elastic strain and a plastic strain according to the plastic theory

Thermodynamic Framework
Based on the mechanical tests and microscopic observation, it has been well confirmed that the failure of brittle sandstone can be attributed to the coupling between the irreversible deformation and damage induced by microcracks. The basic physical process of damage is the initiation, propagation, and coalescence of microcracks. In the compression stress state, the frictional sliding along the rough surfaces of the microcracks produces the irreversible deformation. In this context, a coupled elastoplastic damage model is more adequate to reproduce the inherent coupling between the plastic and damage dissipation processes. For the brittle rock materials, the assumption of small strains is appropriate. In the isothermal conditions, the total strain tensor ε can be decomposed into an elastic strain ε e and a plastic strain ε p according to the plastic theory ε = ε e + ε p , and dε = dε e + dε p The propagation of microcracks in rock materials is generally with faces parallel to the applied stress which results in stress-induced anisotropic damage of materials [37]. In this paper, we emphasize the formulation of a coupled elastoplastic damage model. For the sake of simplicity, an isotropic damage variable is adopted in this work. Therefore, an internal scalar variable ω is introduced to describe the growth of microcracks. Assuming an isothermal process and small strains, the free energy Ψ, taken as the thermodynamic potential, can be expressed in terms of a set of state variables: elastic strain ε e , internal plastic variable κ, and damage variable ω ψ = ψ(ε e , κ, ω) For any dissipative process, the Clausius-Duhem inequality must be satisfied [38] σ : dε − dψ ≥ 0 where σ is the stress tensor. Substituting the differential of ε and Ψ into inequality (7), one gets This relation should be satisfied for any values of state variables (ε e , κ, and ω), and hence we have the state equation By defining the thermodynamic forces associated with the plasticity (K) and damage (Y), we have To describe the plastic flow, a plastic potential g p = g p (σ, K, ω) should be employed. Besides, a damage potential g ω = g ω (Y, ω) is introduced to describe the damage evolution. Finally, the rate of change of the internal variables can be characterized as where dλ p and dλ ω are the plastic and damage multipliers, respectively. In the general case, a plastic criterion f p = f p (σ, K, ω) is necessary to account for the pressure dependence of the brittle rock. Also, a damage criterion f ω = f ω (Y, ω) is introduced to describe the damage initiation. Therefore, the loading-unloading conditions can be represented by the Kuhn-Tucker conditions with the formulations The first inequality in both Equations (15) and (16) suggests that the thermodynamic forces are within or on the yield surface and the second one indicates that the multipliers are nonnegative. The third equation ensures that the stress states lie on the yield surface during the complete loading or unloading process [13]. The two consistency conditions can be expressed as In the general case ( f p > 0 and f ω > 0), plastic flow and damage evolution take place in a coupled process. If f p > 0 and f ω ≤ 0, only plastic flow occurs. If f p ≤ 0 and f ω > 0, the material is in the elastic damage loading condition.

Microscopic Prediction
Under the general loading conditions (plastic flow is coupled with damage evolution), if the yield criteria ( f p , f ω ) and potential functions (g p , g ω ) are given, the plastic and damage multipliers can be determined by solving the two consistency conditions Equations (15) and (16). Under a thermodynamic framework, the description of the damage conjugate force Y and plastic hard function K are related to the formation of thermodynamic potential (6). However, the assumption of thermodynamic potential when considering the coupled relationship between plasticity and damage has never been theoretically and experimentally justified in the phenomenological constitutive model. Moreover, many parameters are introduced, but have no physical significance. The micromechanical approaches provide an alternative way to deal with the coupled elastoplastic damage problems.

Effective Elastic Properties of Cracked Materials
The essential of micromechanical damage is to determine the macroscopic properties of materials from its microstructure (geometry, number, size, and spatial distribution of defects as well as their evolution laws) via homogenization schemes [39]. In this study, we consider a representative element volume (REV) Ω (with boundary ∂Ω) of a brittle material composed of an isotropic linear solid matrix and sets of defects in different orientations. According to the shapes, the defects in the solid matrix are assumed to be classified into two categories: penny-shaped microcracks and quasi-spherical pores. Then the ensemble of the solid matrix and pore therein is considered as a homogenized porous matrix with stiffness tensor D m . Microcracks are classified into N families with the stiffness tensor D c,r , r = 1, 2 . . . N according to the direction. According to the micromechanical method, both pore-weakened matrix and microcracks are considered as local elastic medium. The effective (homogenized) elasticity tensor of the microcrack-matrix rock system D hom is obtained by taking the average of the local strain over Ω [28,40] where ϕ r and A c,r are the volume fraction and average strain concentration tensor, respectively. The local strain is linear with the macroscopic one ε on ∂Ω, i.e., = A c : ε. is The stiffness tensor D m of rock matrix can be expressed as where k m is the bulk modulus and µ m denotes shear modulus of rock materials. J and K are fourth order symmetric tensors where δ is a second order unit tensor, and I ijkl = 1 2 δ ik δ jl + δ il δ jk is a fourth order unit symmetric tensor. The sets of microcracks are considered as flat ellipsoids with radius a r and aspect ratio ϑ r = c r /a r (Figure 10) where the subscript "r" stands for the rth family. The parameter c is the half-length of the small axis. The volume fraction ϕ r of the rth family can be expressed mathematically in the form where d r = N Ω a 3 r is the well-known crack density parameter of the rth family and can be treated as an internal damage variable [41]. If the macroscopic damage ω is defined as the degradation of the elastic modulus, the microscopic damage variable d is certainly related to ω as ω = ω(d). In this context, the goal to determinate the homogenized elasticity tensor D hom is to evaluate the fourth-order concentration tensor A c,r for each phase. Several homogenization schemes can be found in the literature to deal with the concentration tensor. Considering that the Mori-Tanaka is especially suitable for a microcrack-matrix rock system [42], the homogenization procedure is based on this scheme and the determined homogenized elasticity tensor can be expressed in the form [43] where η 1 and η 2 are parameters only related to the Poisson's ratio ν m of a matrix, namely η 1 = 16 . Finally, the free energy of the matrix-cracks system can be expressed as Materials 2020, 13, x FOR PEER REVIEW 13 of 24 Figure 10. Schematic representation of a penny-shaped crack.

Plastic Characterization
Based on the thermodynamic framework, the state Equation (9) can be expressed as (27) In the microcrack-matrix rock system, the plastic strain induced by the friction sliding along the crack surfaces is directly selected as the internal plastic variable . The thermodynamic force associated with the plastic strain is divided according to Equation (10) : (28) It shows that the thermodynamic force is part of the macroscopic stress σ. The part : in the right hand of Equation (28) is the self-equilibrated stress in the solid matrix according to the homogenization scheme. Therefore, the thermodynamic force is the local stress act on the closed microcracks.
Before formulating the plastic criterion, the plastic strain induced by the microcracks is decomposed into spherical strain tensor and deviatoric strain tensor. If the scalar variable and vector denote the shear dilation and friction sliding, the plastic strain is expressed as The local stress is also convenient to be decomposed into a deviatoric part and a spherical part : : Besides, the free energy (24) can also be expressed in the general form [27] where the tensor D b can be expressed in the form Compared with the phenomenological model [44], a specific form of the locked plastic energy ψ p = 1 2 ε p : D b : ε p can be determined with the micromechanical homogenization procedure. Besides, each parameter has clear physical meaning and can be determined from laboratory experiments.

Plastic Characterization
Based on the thermodynamic framework, the state Equation (9) can be expressed as In the microcrack-matrix rock system, the plastic strain ε p induced by the friction sliding along the crack surfaces is directly selected as the internal plastic variable κ. The thermodynamic force associated with the plastic strain ε p is divided according to Equation (10) It shows that the thermodynamic force σ p is part of the macroscopic stress σ. The part D b : ε p in the right hand of Equation (28) is the self-equilibrated stress in the solid matrix according to the homogenization scheme. Therefore, the thermodynamic force σ p is the local stress act on the closed microcracks. Before formulating the plastic criterion, the plastic strain induced by the microcracks is decomposed into spherical strain tensor and deviatoric strain tensor. If the scalar variable β and vector Γ denote the shear dilation and friction sliding, the plastic strain is expressed as The local stress σ p is also convenient to be decomposed into a deviatoric part and a spherical part where s p = σ p : K and σ p m = trσ p /3 are deviatoric part and spherical part of the local stress. According to the Equation (28), the local stress σ p can also be expressed with the form of macroscopic stress σ where s = σ − 1 3 trσδ and σ m = 1 3 trσ are the deviatoric part and spherical part of the macroscopic stress, respectively.
To describe the local frictional sliding along the surfaces of microcracks, a generalized Coulomb criterion [45] is adopted in terms of local stress σ p where η ϕ is the generalized friction coefficient due to the roughness of the crack surfaces. Equation (32) can also be rewritten in the form of macroscopic global stress σ in Equation (31) Under triaxial compression condition, the damage d at the peak point reaches d c . The relationship between maximum principal stress σ 1 and minimum principal stress can be determined from (33) (under the sign convention in geomechanics) where χ is a constant which can be written with the parameters related to the Yong's modulus and Poisson's modulus [43] In geomaterials, the micromechanical approach always introduces the non-associated plastic potential to describe the volumetric dilatancy during the failure process. However, the micromechanical analysis which uses the back stress term D b : ε p to realize the hardening/softening behavior gives an alternative method [46]. Therefore, an associated flow rule g p = f p is utilized to describe the evolution of plastic strain. According to the plastic theory (12), the rate of the plastic strain can be expressed as where V = s p ||s p || is the flow direction. Compare Equation (36) with the rating form of plastic strain in Equation (29), one gets dΓ = dλ p V, dβ = dλ p η ϕ The current value of plastic strain can be expressed in terms of the cumulated values of the plastic multiplier

Damage Characterization
The damage criterion is a function of the damage conjugate force Y within the framework of irreversible thermodynamics. According to the Equation (11), the damage conjugate force Y can be determined with the function of free energy ψ.
Inspired by the work of Pensee et al. [47], a liner damage criterion is considered where R(d) is the damage energy release threshold at a given value of the damage. Based on the previous work by Zhu and Shao [45], the rock strength is attained when R(d) takes its maximum value at d = d c . d c is the microscopic damage variable at the peak stress. Also, the model should reflect the strain hardening behavior when d ≤ d c while the model should predict the damage softening behavior when d > d c . The following power form is adopted for R(d) where R(d c ) is the maximum value of R(d) at d = d c . The ratio ξ is defined by ξ = d/d c . In this study, the associated damage potential is adopted which means that the damage potential is equal to the damage criterion, namely g d = f d .

Integration by Return Mapping Algorithm
Numerical implementations of the inelastic constitutive model require the stress state to be corrected and returned onto the yield surface. Generally, there are two types of numerical integration techniques: explicit algorithms and implicit algorithms. The explicit algorithms parameter updates at the beginning of a given time step. The disadvantage is that when the time step is decreased, a non-convergent and infinite loop may happen during iterations [48]. The return mapping algorithm, first proposed by Simo and Ortiz [49], is a typical implicit algorithm and is widely used for the numerical implementation of elastoplastic models in the programming of FEM. This method can reach an asymptotic global quadratic convergence rate when using the full Newton-Raphson iteration method [50].
In this study, this algorithm is extended to solve the coupled elastoplastic damage problem. If the implementations are based on the strain-controlled strategy, a trial calculation is tested with a new strain increment. If the yield condition is reached, the return mapping algorithm is utilized to drive the Finally, the stress, elastic strain, plastic strain, and damage can be exactly updated by the returning mapping algorithm.

Implementation Procedure
The previous numerical integration scheme applied to phenomenological model [44] is also implemented to the micromechanical based model in this work. The flowchart of the numerical algorithm is summarized in Table 3. Table 3. Flowchart of return mapping algorithm of the micromechanical based model.
Plastic correction Determine the plastic multiplier δλ p with Equation (47) while δλ d = 0. Update the variables: Update the macroscopic stress and strain: 3 η ϕ δ δλ p Elastoplastic damage coupling correction Calculate the multipliers {δλ p , δλ p } using the Equation (48) Update the variables: Update the macroscopic stress and strain:

Numerical Simulations
In this section, the comparison between the experimental results of the sandstone and the numerical simulations data are presented. The proposed micromechanical based elastoplastic damage model has been implemented as a user-defined in a home-made c language code. At present, only a Gauss integration point inside a finite element is studied which indicates that the simulation results are independent of the element type. In the future work, the subroutine can be easily embedded into finite element software or finite difference software to analyze the safety and stability of hydropower station in excavation and operation conditions.

Identification of Model Parameters
As mentioned above, the phenomenological models always introduce so many parameters to account for the coupled relationship between plasticity and damage. However, the proposed coupled elastoplastic damage model in this work only contains five parameters and all the parameters can be determined from experiments in the laboratory. The feature of the developed model makes it more suitable for the engineering application. The calibration of the parameters is discussed below.
Equation (38) shows that the major principal stress σ 1 increases linearly with the increase of minor principal stress σ 3 which indicates that the friction criterion is similar to the Mohr-Coulomb criterion in the formulation. Therefore, the generalized friction coefficient η ϕ and parameter R(d c ) can be expressed with the cohesion c and internal friction angle ϕ in Equation (1) With the parameters (c and ϕ) listed in Table 1, the determined η ϕ and R(d c ) are 1.56 and 0.009, respectively.
The relationship between bulk modulus k m and the shear modulus and µ m with the elastic constants E s and υ can be expressed as: The Young's modulus E s and Poisson's ratio υ of the sandstone is shown in Figure 6. For the sack of simplicity, the average values of the Young's modulus E s = 15.86 GPa and Poisson's ratio υ = 0.198 are selected.
According to Lockner [51], the maximum damage variable d c (refers to the crack density at peak point) can be determined from the acoustic emission tests. Unfortunately, we do not have such a device in our laboratory. To calibrate the parameter d c , parametric studies on the influence of d c upon material responses are carried out and the results are given in Figure 11. The parameter d c influences the shape of the stress-strain curve in the post-peak region. With the increases in d c , the ductile characteristic is more significant. Even the confining pressure increases to 40 MPa, the sandstone is with brittle failure. According to the sensitivity analysis results of d c and the shapes of the stress-strain curve of the test sandstone, we take approximately d c = 1.0 for the test sandstone.
The relationship between bulk modulus and the shear modulus and with the elastic constants Es and can be expressed as: The Young's modulus Es and Poisson's ratio of the sandstone is shown in Figure 6. For the sack of simplicity, the average values of the Young's modulus Es = 15.86 GPa and Poisson's ratio = 0.198 are selected.
According to Lockner [51], the maximum damage variable dc (refers to the crack density at peak point) can be determined from the acoustic emission tests. Unfortunately, we do not have such a device in our laboratory. To calibrate the parameter dc, parametric studies on the influence of dc upon material responses are carried out and the results are given in Figure 11. The parameter dc influences the shape of the stress-strain curve in the post-peak region. With the increases in dc, the ductile characteristic is more significant. Even the confining pressure increases to 40 MPa, the sandstone is with brittle failure. According to the sensitivity analysis results of dc and the shapes of the stressstrain curve of the test sandstone, we take approximately dc = 1.0 for the test sandstone.
Finally, the parameters of the model are listed in Table 4.  Figure 11. Sensitivity analysis on the parameter dc.

Simulation of Experimental Results of Test Sandstone
The results of each modeled test are compared to the experimental results of the sandstone ( Figure 12). There is a good concordance between the modeled results and experimental data under different confining pressures. However, some systematic differences for the tests are also observed. The difference between the numerical curve and the experimental data under uniaxial compression stems from the compelling concave upward stage of sandstone under low confinements. The experimental data show a marked decrease of the stress in the post-peak region while a smooth fall of stress is received from the numerical model. This could be attributed to the shear localization and stress-induced fractures resulting from the coalescence near the peak stress of clusters of oriented microcracks. This behavior of sandstone is behind the topic of this study and will be developed in our ongoing works. Finally, the parameters of the model are listed in Table 4.

Simulation of Experimental Results of Test Sandstone
The results of each modeled test are compared to the experimental results of the sandstone ( Figure 12). There is a good concordance between the modeled results and experimental data under different confining pressures. However, some systematic differences for the tests are also observed. The difference between the numerical curve and the experimental data under uniaxial compression stems from the compelling concave upward stage of sandstone under low confinements. The experimental data show a marked decrease of the stress in the post-peak region while a smooth fall of stress is received from the numerical model. This could be attributed to the shear localization and stress-induced fractures resulting from the coalescence near the peak stress of clusters of oriented microcracks. This behavior of sandstone is behind the topic of this study and will be developed in our ongoing works. Model results (black) compared to experimental responses (grey). The cyclic experimental results are shown in red, the left is stress-radial strain curves and the right is stress-axial strain curves.

Simulation of Experimental Results of Shandong Red Sandstone
To illustrate the utility and consistency of the proposed coupled elastoplastic damage model, the mechanical properties of Shandong sandstone are here investigated. Figure 13 shows the experimental results and numerical simulations of triaxial compression tests performed on Shandong sandstone. Following the procedure presented in Section 5.1, the model's parameters of Shandong sandstone are determined and listed in Table 5.
The comparisons show a good agreement between numerical predictions and experimental results for five levels of confining pressures, namely 5, 20, 35, 50, 65 MPa. The strain hardening in the Model results (black) compared to experimental responses (grey). The cyclic experimental results are shown in red, the left is stress-radial strain curves and the right is stress-axial strain curves.

Simulation of Experimental Results of Shandong Red Sandstone
To illustrate the utility and consistency of the proposed coupled elastoplastic damage model, the mechanical properties of Shandong sandstone are here investigated. Figure 13 shows the experimental results and numerical simulations of triaxial compression tests performed on Shandong sandstone. Following the procedure presented in Section 5.1, the model's parameters of Shandong sandstone are determined and listed in Table 5.
The comparisons show a good agreement between numerical predictions and experimental results for five levels of confining pressures, namely 5, 20, 35, 50, 65 MPa. The strain hardening in the pre-peak region and strain softening in the post-peak region are clearly produced. The dependence of mechanical properties of sandstone on confining pressure is well simulated. In conclusion, the developed model can simulate the coupled elastoplastic damage properties of sandstone at the scale of the rock sample.
pre-peak region and strain softening in the post-peak region are clearly produced. The dependence of mechanical properties of sandstone on confining pressure is well simulated. In conclusion, the developed model can simulate the coupled elastoplastic damage properties of sandstone at the scale of the rock sample.  [52]).

Conclusions
The mechanical properties of studied sandstone in this work is significantly important to the stability and safety of the hydropower station in Southwest China. A series of uniaxial and triaxial compression tests were carried out in our laboratory. The strength and deformation characteristics were investigated. The microstructure of the failure specimen was examined through SEM. The

Conclusions
The mechanical properties of studied sandstone in this work is significantly important to the stability and safety of the hydropower station in Southwest China. A series of uniaxial and triaxial compression tests were carried out in our laboratory. The strength and deformation characteristics were investigated. The microstructure of the failure specimen was examined through SEM. The results showed that the complete stress-strain curves of the sandstone can be divided into four stages.
The volume of the sandstone changes from compaction to dilatancy during the failure process. The critical stresses all depend on the confining pressure. The Young's modulus increases nonlinearly with the confining stress while the relationship between the Poisson's ratio and the confining stress is not clear. The plastic strain increment N φ and the dilation angle φ show a negative response to the confining pressure. The SEM images present the growth and frictional sliding of the cracks induced by the stress.
According to the experimental results, a coupled elastoplastic damage model was developed based on the irreversible thermodynamic framework. Specific plastic and damage criteria based on the micromechanics are proposed to describe the physical process of propagation and frictional sliding of microcracks. Compared with the phenomenological model, the model developed in this paper can reflect the physical mechanism of the failure of the sandstone that observed from the tests. Besides, the model only has five parameters and each one can be determined from laboratory tests. The general constitutive integrated formulations of the coupled elastoplastic damage model were deduced based on the return mapping algorithm. The model was validated through a comparison of the numerical simulation results to the experimental data. A good concordance between the modeled results and experimental data suggests that this model can provide a good representation of the nonlinear behavior of the sandstone.
In this study, we mainly present the experimental results of the sandstone, and a micromechanical based elastoplastic damage model was developed. In the near future, the proposed constitutive model will be implemented into the finite element method (FEM). Combining the laboratory tests with the in-situ tests, the mechanical parameters can be determined. Using the developed model, the excavation damaged zones (EDZs) can be estimated, and it could provide a reference for the design and construction of a huge hydropower project.