A Coupled Elastoplastic Damage Model for Clayey Rock and Its Numerical Implementation and Validation

Institute of Unconventional Oil & Gas, Northeast Petroleum University, Daqing 163318, China Key Laboratory of Oil & Gas Reservoir and Underground Gas Storage Integrity Evaluation of Heilongjiang Province, Daqing 163318, China School of Urban Construction, Yangtze University, Jingzhou, 434023 Hubei, China State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan, 430071 Hubei, China State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China


Introduction
Clayey rock is always selected as the geological barriers for radioactive waste disposal and considered the effective caprock for CO 2 sequestration and oil and gas reservoir [1][2][3]. Clayey rock in its natural state exhibits good plastic deformation ability and very low permeability as geological barriers. One of the main concerns is that, due to underground excavation or fluid injection, the favorable properties of clayey formation change and the host rock loses part of its barrier function, thus negatively influencing the performance of a repository [4,5]. In addition, the damage associated with excavation or injection causes obvious degradation of mechanical and hydraulic properties of clayey rock and thus plays an important role in changing the stress and seepage fields in the practical engineering.
The damage and integrity behaviors of clayey rock are the key factors in studying the storage and environmental safety of geological disposal in deep repository. The underground excavation can redistribute the stress and generate irreversible damage around a gallery. The safety evaluation of a repository in clayey rock requires accurate prediction of the mechanical perturbations and damage zone, which is called as excavation damaged zone (EDZ) [6][7][8]. Therefore, a suitable constitutive model, which can characterize accurately the mechanical behavior of clayey rock, is very important for repository design [9][10][11].
Clayey rock is one type of sedimentary rock, and some characteristics of its mechanical behavior are similar to those of soft rock [4]. Laboratory experiments show that clayey rock has significant strain-softening behavior and good plastic deformation ability [2,12,13]. The strain-softening behavior of clayey rock has been studied from micro-and macromechanical viewpoints [2,4]. The micromechanical model can reproduce not only the crack opening, interaction between crack surfaces, and crack closure but also the process of strain localization or shear banding. However, the micromechanical model has limited use in practical engineering due to large number of excessive preoccupied cracks and complex modelling approaches. Therefore, the macromechanical model has been mainly used in underground engineering and can be easily applied to simulate engineering problems, for example, determining the distribution of EDZ and extent of failure. As for plastic deformation and strainsoftening behavior for clayey rock, the modelling approaches mainly include the following [7,[14][15][16][17]: (1) the strength parameters are considered to be weakened during strain softening based on laboratory tests and field observations and (2) the damage concept is introduced into the classical continuum theory to describe the degradation of stiffness and strength of clayey rock.
A large number of models, such as Cam-Clay, Drucker-Prager, Mohr-Coulomb, Cap model, and bounding surface model [18][19][20][21], have been developed to describe the failure behavior of clayey rock. Among these models, the Mohr-Coulomb failure criterion has been widely used for geotechnical applications due to its easy implementation. Many researchers have noted that modelling the failure behavior of clayey rock should take into account the damage that governs the complex mechanical behavior [21,22]. The laboratory tests show that plastic deformation and damage are coupled [2,4,11,21,23], indicating that the plastic yield criterion should be used simultaneously with the damage criteria to consider the changes of mechanical properties of clayey rock. The damage of rock has obvious effect on the stress, strain, elastic stiffness, strength, yield function, potential function, and softening parameters. In the coupled damage model developed recently for clayey rock [4,6,[8][9][10], in order to build the complex relationship between plastic flow and damage evolution, some assumptions are made on the basis of laboratory tests. To the best of the authors' knowledge, there is little available research on proposing a coupled damage constitutive model for clayey rock based on laboratory tests. This is also the motivation of the present work.
The paper is organized as follows. In Section 2, the MMC criterion, which combines the CMC criterion with the maximal tensile stress criterion, is proposed to describe the yield and plastic flow behavior of clayey rock. Section 3 presents a damage constitutive model, which considers the stiffness degradation, strain-hardening, and plastic softening effect. Numerical implementation of the proposed model in finite element software ABAQUS is presented in Section 4. Finally, in Section 5, numerical results are provided to show the effectiveness of the proposed model.

Modified Mohr-Coulomb Criterion
The Mohr-Coulomb criterion is still widely used for a large number of routine design in the geotechnical engineering, which is the lower limit of all the linear strength criterion [24]. According to many experimental studies [4,11,25,26], the Mohr-Coulomb criterion can well reflect the elastoplastic mechanical behavior of clayey rock. However, among the popular finite element software, such as ABA-QUS, ANSYS, and COMSOL, the built-in Mohr-Coulomb model cannot well describe the dilatancy characteristics and the truly associated flow properties of geomaterial because of the different expressions of yield and potential functions [4].
The CMC criterion is generally written in terms of stress invariants [11,27]: where σ m , σ, c, and ϕ denote the average stress, equivalent stress, cohesion, and friction angle, respectively.K ′ is a function of Lode angle θ and friction angle ϕ.
In Equations (1) and (2), σ m , σ, and θ are defined by where σ 1 , σ 2 , and σ 3 are the three principal stresses and J 2 and J 3 are the second and third invariant, respectively, of the deviatoric stress.
The CMC criterion does not consider the real tensile strength of rock, and it has the limitation for the actual tensile shear zone. In Figure 1, σ * m is the three-dimensional tensile strength by the CMC criterion, which is larger than the actual uniaxial tensile strength of rock f t . Assuming that the threedimensional tensile strength of rock is equal to the uniaxial tensile strength, σ t = f t , the tensile Mohr-Coulomb criterion is defined as [4] where f t is the uniaxial tensile strength of geomaterials and σ 1 = ð2/ ffiffi ffi 3 p Þ σ sin ðθ + 120°Þ + σ m . Thus, Equation (4) is rewritten in terms of stress invariants: Because the proposed MMC criterion considers both tensile strength and shear strength of geomaterials, the yield surface is a result of combination of two different 2 Geofluids failure criteria: shear failure and tensile failure. The yield function is defined as follows [11] (as shown in Figure 1): where 0 ≤ m ≤ 1 is a parameter reflecting the tensile strength of geomaterials. When m = 0 indicates a higher tensile strength, the MMC yield criterion in Equation (6) becomes the CMC criterion; m = 1 denotes that there is no tensile strength. In addition, the parameter m can smooth the vertex of the yield surface and avoid the numerical divergence and slow convergence. The Mohr-Coulomb failure surface is a hexagonal cone in the principal stress space with six sharp corners in the octahedral plane [11,28]. The six singularities occur at θ = ±30°in the octahedral plane (as shown in Figure 2), making the numerical convergence difficult. To approach the CMC yield surface, the parameter K ′ is expressed by using as a piecewise function: where A = 1 3 cos θ T 3 + tan θ T tan 3θ T ð The variable θ T denotes the transition angle in the vicinity of the six singularities, and its value is in the range of 0°≤ θ T ≤ 30°. In the present calculations θ T = 25°is adopted. The yield function rounded in both the meridional and octahedral planes is modified by adjusting the parameters m and K ′ in Equations (6) and (7), where the yield surface is continuous and differentiable for all stress states.
The plastic potential function can be defined from the yield function as follows: where φ is the dilatancy angle. Similarly, K ′ is a function of Lode angle θ and dilatancy angle φ. Associated flow occurs  Figure 3 schematically shows the stress-strain curve of clayey rock observed in laboratory tests. The mechanical properties of clayey rock are affected by fractures growth, cementation degree, and stress state [4]. When the microcracks grow, the clayey rock undergoes nonlinear deformation, leading to deterioration of its macroscopic elastic properties and strength. The interaction between the nonlinear deformation and the damage can be dealt with by a damage model. For the purpose of simplicity, the stressstrain curve is divided into four stages, i.e., elastic deformation, strain hardening, strain softening, and plastic flow. Due to the complexity of stress-strain behavior under triaxial stress state, the conventional elastoplastic constitutive model cannot accurately describe this process. Thus, the present work uses a piecewise method to fully characterize the mechanical behavior of clayey rock. Based on the damage mechanics theory, an EPD constitutive model, which contains elastic law, elastic damage law for strain hardening, and plastic damage law for strain softening, is developed. According to the position of peak stress point, the complete stress-strain curve is divided into two zones, i.e., elastic and plastic. For the plastic zone, the Mohr-Coulomb failure criterion is adopted to describe the progressive plastic deformation and strain-softening process.

Model Description.
Although anisotropic stress-induced damage has been observed in some laboratory tests on clayey rock, it is not obvious in this study [4,18,21]. Therefore, a scalar isotropic damage variable is used in the present coupled EPD model. Following the convention of continuum damage mechanics [29], the effective stressσ is defined as follows from the Lemaitre strain equivalent hypothesis: where σ is the nominal stress and Ω is the damage variable.
According to the MMC failure criterion, the yield function in terms of effective stress can be written as By substituting Equation (11) into Equation (6), the yield function in terms of nominal stress can be determined by Similarly, the potential function can be defined in terms of nominal stress as 3.2. General Stress and Strain Relationship. Generally, the assumption of small strain is suitable for clayey rock. The stress-strain relationship of clayey rock shows obvious nonlinear characteristics before peak strength [2,4,21]. Although the strain during the hardening after unloading cannot be restored completely, the plastic strain in this stage is so small that it can be neglected. For the purpose of simplicity, the elastic damage occurs in the strainhardening stage and the total strain dε ij is decomposed into the following two parts: where dε R ij is the reversible strain, dε e ij is the elastic strain, and dε ed ij is the elastic damage strain.

Geofluids
The elastic damage strain dε ed ij is deduced from the change of elastic parameters during the strain hardening [17]: where C ijkl is the elastic stiffness tensor and σ kl is the stress tensor.
When the stress in clayey rock reaches the peak strength, the rock will undergo plastic deformation which is coupled with the damage. For the stages of strain softening and plastic flow, the irreversible strain dε IR ij is decomposed into plastic and damage parts: where dε p ij is the plastic strain and dε d ij is the plastic damage strain.
According to the potential function in Equation (14), the irreversible strain dε IR ij can be defined as follows [8,10]: where dλ is the plastic multiplier. The general stress and strain relationship can be expressed as follows: 3.3. Definition of Damage. According to stress-strain relationship of clayey rock as shown in Figure 3, the curve is linear in the elastic stage and nonlinear in the strain-hardening stage. The transition point delineating these two stages is considered the elastic damage starting point. The elastic damage stage ends and the plastic damage commences when the peak strength is reached. The evolution equation of elastic damage is where β 1 is the elastic damage parameter and e 0e is the energy factor of the elastic damage starting point. The energy factor e is defined as [4] where C 0 ijkl is the undamaged elastic stiffness tensor. The peak stress of the stress-strain curve is defined as the threshold of plastic damage. Actually, plastic damage is dependent on the history of plastic deformation and its rate tends to stabilize gradually with the increase of accumulated deformation. The evolution equation of plastic damage is written as follows: where e 0p is the energy factor corresponding to the plastic damage starting point and α 2 β 2 are the plastic damage parameters.

Evolution of Model Parameters.
According to the continuum damage mechanics, the degradation of elastic modulus during loading is defined as follows [4,29]: where E is the damaged elastic modulus, E 0 is the initial undamaged elastic modulus, and Ω is the sum of elastic and plastic damage variables, i.e., Ω = Ω e + Ω p . Equation (23) shows that the elastic modulus decreases gradually during the damage deformation. When the rock approaches a complete damage, the damage variable Ω is close to 1 and thus, the elastic modulus tends to be 0. However, this is not in accordance with engineering application, where the elastic modulus of clayey rock retains its residual value even after a significant damage. In order to overcome this problem, here the elastic modulus is expressed by using a bilinear function of the damage variable: where E r is the residual elastic modulus of clayey rock and Ω lim is the critical damage value and Ω lim = 0:99 is adopted in this study. During the strain-softening stage, the strength of clayey rock decreases while the plastic damage variable and plastic deformation increase. As the friction angle has small change during the strain-softening process [4,30], it is assumed that the strength parameter evolution is only described by the cohesion, which is expressed as follows: where c 0 and c r denote the initial and residual cohesions, respectively; η ∈ ð0, 1 is a model parameter controlling the slope of the strain-softening curve. The fully BEII method is used in this study to return the stresses for predicting the yield surface when the stresses exceed the yield strength [31][32][33]. Given the responses at time t n , i.e., stress σ n , damage variable Ω n , and a total strain increment Δε n+1 , the objective is to determine these state variables σ n+1 , ε n+1 , and Ω n+1 at time t n+1 . The BEII method 5 Geofluids includes mainly three steps, i.e., elastic predictor, plastic corrector, and damage corrector.

Numerical Implementation
In the elastic predictor step, the scalar damage Ω n is assumed to be a constant. A trial damage is defined as and a trial stress σ trial n+1 is defined based on the elastic predictor where C is the elastic stiffness matrix, which is a function of the trial damage Ω trial n+1 . The onset of plastic flow and deformation is determined by the loading condition, which is expressed as follows: The trial stress at point B, σ trial B , is obtained by the elastic predictor (as shown in Figure 4) where Δσ e is the increment of elastic stress. In order to return the trial stress σ trial n+1 for predicting the yield surface, the stress increment is defined as where b = ð∂Gðσ trial n+1 , Ω trial n+1 Þ/∂σÞ n+1 . The derivatives of the potential function and yield function are given in the appendix.
Thus, the trial stress at point C is expressed as The first-order Taylor expansion of the yield function at point B gives In the plastic corrector step, the scalar damage Ω trial n+1 is regarded as constant. Then, Equation (32) reduces to be where a = ∂F/∂σ and F B = Fðσ trial B , Ω trial n+1 Þ.
Here the plastic multiplier Δλ can be determined explicitly as The backward Euler algorithm is based on the equation where σ B is the elastic trial stresses and the variables with subscript C are related to the current configuration. A starting estimate for σ C is defined as Generally, this starting value of stress σ C does not satisfy the yield function and further iterations will be required because the normal at the trial position B is not equal to the final normal (as shown in Figure 4). To control the iterative loop, a vector r is used to represent the difference between the current stresses and the backward Euler calculations, i.e., The iterations continue until the norm of vector r is small enough (almost zero) while the final stresses should satisfy the yield criterion.
With the trial stresses σ B being kept fixed, a truncated Taylor expansion can be applied to Equation (37). A new residual vector r N is expressed as follows: where r 0 is the residual vector at point B, _ σ is the change in σ, _ λ is the change of Δλ, and N is the number of iteration.

Geofluids
Setting r N to zero gives where I is the unit matrix; Q = I + ΔλCð∂b/∂σÞ. Similarly, a truncated Taylor series on the yield function is By substituting Equation (39) into Equation (40), the change of Δλ is determined by and the plastic multiplier after N times of iteration is obtained: The modified stress vector σ C after N times of iterations is The strain vector ε n+1 at time t n+1 is expressed as The third step in the BEII method is the damage corrector. The updated damage at time t n+1 is where ΔΩ = ΔΩ e + ΔΩ p . After the damage corrector is completed, the final stress at time t n+1 is determined:  (35), the standard back-Euler algorithm is expressed as It should be noted that the variables without subscript "C" in Equation (47) are related to the current configuration.
When the change of damage in plastic corrector is omitted, differentiation of Equation (47) gives which is simplified to be where R = Q −1 C.
To retain the current stress σ on the yield surface, _ F should be zero. The consistency conditions of yield function F is expressed as follows: By substituting Equation (49) into Equation (50), _ λ is determined by The consistent tangent stiffness matrix,D ct , can be obtained by substituting _ λ into Equation (49): 4.2. Procedure of the Constitutive Integration Algorithm. In the fully implicit backward Euler algorithm [31,34], the increment of irreversible strain and damage variable is calculated at the end of increment step n. The constitutive model integration algorithm can be expressed as Equation (53) is a system of nonlinear equations of ε ðn+1Þ , ε IR ðn+1Þ , and Ω ðn+1Þ . At any time t n , a set of values ðε ðnÞ , ε IR ðnÞ , Ω ðnÞ Þ and strain increment Δε = Δt _ ε are given. The updated variables come from the convergence value at the end of the previous time step, which achieves the effect of avoiding nonphysical meaning. The solution of the 7 Geofluids nonlinear equations is solved by Newton-Raphson iterative method. For the updated variables, the superscript k represents the number iteration and the subscript n represents the increment step.
The stress update algorithm flow is as follows: Step 1. Set the initial value. The initial value of irreversible strain and damage variable is the convergence value at the end of the last load step. The incremental value of plastic parameter is set to zero, and the elastic trial stress is calculated.
Step 2. Check the yield condition and convergence at the iteration number k.
For a given stress error tolerance TOL, if F ðkÞ ðn+1Þ < TOL, convergence; otherwise, go to Step 3.
Step 3. Calculate the increment of plastic parameters, then the increment of stress and damage variables can be obtained.

Secondary Development of UMAT Subroutine.
The numerical formulations discussed above are programmed as a UMAT subroutine in ABAQUS with FORTRAN language [28]. The material properties are defined in this UMAT subroutine, in which the stress and other related state variables are updated at the end of each time step. The UMAT subroutine can also work with a user subroutine to redefine field variables at a material point (USDFLD) to update the field variables [4,21]. The stress, elastic, and plastic strain are obtained by iteration in the UMAT subroutine, and then, the state variables such as plastic parameters and damage variable are updated. The change rate of _ σ with respect to _ ε is provided by the constitutive Jacobian (DDSDDE) in the UMAT subroutine as the Jacobian matrix of the material constitutive model. The flow chart of the UMAT subroutine for the proposed model is shown in Figure 5.

Application of the Proposed Model
5.1. Model Validation. The proposed constitutive model was developed and implemented in the UMAT subroutine of ABAQUS. In order to test the calculation ability, accuracy, and efficiency of this UMAT subroutine, the performance of the constitutive formulation was evaluated for uniaxial tension and confined compression scenarios. In the special case when no damage is mobilized, the proposed EPD model reduces to an ideal elastoplastic model. For the purpose of numerical illustration, the numerical results using the UMAT subroutine based on the MMC criterion were compared with those by using the built-in Mohr-Coulomb criterion in ABAQUS. Figure 6 shows the geometrical model for one cylinder rock sample with a diameter of 38 mm and a height of 76 mm. The bottom of the sample is fixed, and a vertical displacement load is applied on the top surface of the sample. The elastic modulus E = 700 MPa and Poisson's ratio μ = 0:25. The plastic properties of rock include the friction angle, cohesion, and dilatancy angle, whose values are φ = 18°, c = 0:3 MPa, and ϕ = 18°, respectively. For the nonassociated flow rule, the dilatancy angle is ϕ = 0°. The parameter m = 0:05 is used in the MMC criterion for avoiding the numerical divergence in this study. Figures 7 and 8 compare the stress-strain relationship for uniaxial and triaxial compression, respectively, between the present MMC model and the built-in Mohr-Coulomb model in ABAQUS. It can be found that the numerical results of deviatoric stress and volumetric strain predicted from both models are in good agreement. For a given axial strain under a specific confining pressure, the yield stress obtained from the MMC model (with the UMAT subroutine) is slightly smaller than that obtained from the built-in model. This is due to the fact that the yield surface of the MMC criterion is inside that of built-in Mohr-Coulomb criterion. In particular, the volumetric strain predicted from both models agree very well. This indicates that the proposed MMC model is feasible to describe the ideal plastic response under compression conditions.
The mechanical response of the sample under uniaxial tension is shown in Figure 9 for different values of m reflecting the tensile strength of rock. It can be seen that when m = 0:05, the numerical results calculated from the present model are in agreement with those obtained by

Triaxial Tests and Numeric Simulation on Clayey Rock.
In this subsection, the mechanical characteristics of clayey rock are studied by using the proposed EPD model.
The clayey rock in this study is a stiff clay. The undisturbed samples (Φ38 × 76 mm) used in the undrained triaxial tests are extracted from clayey rock formation at a depth of 247 m. Considering the low permeability of clayey rock (with an order of magnitude 10 −19 m 2 ), these tests can be used to estimate the undrained shear strength by analogy with quick excavation and tunneling for repository [18,35].
The tests were performed in a triaxial apparatus with a confining pressure between 0.89 MPa and 5.42 MPa according to following sequences [4,36]. First, samples were loaded in increment of approximately 0.5 MPa to the scheduled confined pressure, and then, a valve is opened to apply a back pressure (nominally 50% of cell pressure) allowing consolidation to commence. Once consolidation and pore pressure     dissipation were completed, the specimens were sheared by the application of small increments in axial stress. The finite element model is shown in Figure 6. The boundary conditions and simulation sequences are the same as those in laboratory tests. Simulation of undrained triaxial test involves solving a coupled hydromechanical problem, and the computational precision depends mainly on the mesh size and time step, especially the stage at the beginning of calculation. In order to ensure a uniform excess pore pressure in the sample, the initial time step Δt is estimated based on the following formula [28,37]: where E ′ is the drained elastic modulus, γ w is the density of pore water, k is the permeability of clayey rock, and Δh is the characteristic length of the mesh. The elastic modulus, Poison's ratio, peak cohesion, friction angle, and two energy factors of the damage starting point can be obtained directly from tested results, which are listed in Tables 1 and 2. The unknown parameters in the proposed EPD model can be determined by using a back analysis [4], which are listed in Table 3.
The stress-strain curves under varying confining pressures from 0.89 to 5.42 MPa are shown in Figure 10. The undrained triaxial tests show that clayey rock exhibits obvious strain hardening and strain softening during the shearing deformation, and the stress-strain curves can be classified into four stages. Moreover, large plastic deformation during strain softening is the main feature of the mechanical response of clayey rock. It can be observed that the numerical results from the proposed model are in good agreement with the test results. This indicates that the proposed EPD model is capable of effectively capturing the strain-hardening and strainsoftening characteristics of the clayey rock. When the deformation of clayey rock reaches the initial elastic damage point, the elastic damage increases slowly and then remains constant during plastic deformation stage. From Figure 10, it can also be seen that the total damage variable increases gradually with the increase of axial strain and the plastic damage behavior plays an important role in the irreversible deformation.

Conclusions
In the present work, a new coupled elastoplastic damage (EPD) model was developed by using a modified Mohr-Coulomb criterion which considers the maximal tensile strength. Elastic and damages were introduced to describe the strain-hardening and strain-softening processes, respectively.
Based on the implicit Euler stress integration algorithm, the constitutive integrating formulations and the consistent stiffness matrix for the coupled EPD model were deduced. For engineering application, the proposed model was coded as a subroutine UMAT in ABAQUS.
To validate the EPD model, the triaxial compression test and uniaxial tension test were first simulated without considering the damage effect. The comparison of the numerical results between the EPD model (with the UMAT subroutine) and built-in Mohr-Coulomb model in ABAQUS indicates that the proposed model is capable of effectively describing the tensile and compression responses of rocks. Then, the coupled EPD model is used to simulate the undrained triaxial tests of clayey rock. Comparisons between numerical simulations and experimental data show that the proposed model can characterize the strain hardening, strain softening, and plastic flow of clayey rock very well.
Further experimental studies are necessary to obtain more data for a better understanding of the plastic and damage mechanical behavior of clayey rock, especially when some parameters in the proposed model have strong dependence on the confining pressure. In the future, the proposed model will be extended to consider the effect of confining pressure.