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

Authors: Shanpo Jia; Zhenyun Zhao; Guojun Wu; Bisheng Wu; Caoxuan Wen
Journal: Geofluids
 (2020)

---

## Abstract

This paper presents a new constitutive model for describing the strain-hardening and strain-softening behaviors of clayey rock. As the conventional Mohr-Coulomb (CMC) criterion has its limitation in the tensile shear region, a modified Mohr-Coulomb (MMC) criterion is proposed for clayey rock by considering the maximal tensile stress criterion. Based on the results of triaxial tests, a coupled elastoplastic damage (EPD) model, in which the elastic and plastic damage laws are introduced to describe the nonlinear hardening and softening behaviors, respectively, is developed so as to fully describe the mechanical behavior of clayey rock. Starting from the implicit Euler integration algorithm, the stress-strain constitutive relationships and their numerical formulations are deduced for finite element implementation in the commercial package ABAQUS where a user-defined material subroutine (UMAT) is provided for clayey rock. Finally, the proposed model is used to simulate the triaxial tests and the results validate the proposed model and numerical implementation.


---

## Body

## 1. Introduction

Clayey rock is always selected as the geological barriers for radioactive waste disposal and considered the effective caprock for CO2 sequestration and oil and gas reservoir [1–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–8]. Therefore, a suitable constitutive model, which can characterize accurately the mechanical behavior of clayey rock, is very important for repository design [9–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 strain-softening behavior for clayey rock, the modelling approaches mainly include the following [7, 14–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–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–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 Section2, 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.


## 2. 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 ABAQUS, 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]:
(1)F=σmsinϕ+σ¯K′−ccosϕ=0,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 ϕ.
(2)K′=cosθ−13sinϕsinθ.In Equations (1) and (2), σm, σ¯, and θ are defined by
(3)σm=σ1+σ2+σ33,σ¯=J2,θ=13sin−1−332J3σ¯3,−30°≤θ≤30°,where σ1,σ2, and σ3 are the three principal stresses and J2 and J3 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 Figure1, σm∗ is the three-dimensional tensile strength by the CMC criterion, which is larger than the actual uniaxial tensile strength of rock ft. Assuming that the three-dimensional tensile strength of rock is equal to the uniaxial tensile strength, σt=ft, the tensile Mohr-Coulomb criterion is defined as [4]
(4)σ1≥ft,where ft is the uniaxial tensile strength of geomaterials and σ1=2/3σ¯sinθ+120°+σm.Figure 1
Modified Mohr-Coulomb criterion in the meridional plane.Thus, Equation (4) is rewritten in terms of stress invariants:
(5)F=23σ¯sinθ+120°+σm−ft=0.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 failure criteria: shear failure and tensile failure. The yield function is defined as follows [11] (as shown in Figure 1):
(6)F=σmsinϕ+σ¯K′2+m2c2cos2ϕ−ccosϕ=0,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:
(7)K′=A¯−B¯sin3θ,θ>θT,cosθ−13sinϕsinθ,θ≤θT,(8)K′=A¯−B¯sin3θ,θ>θT,cosθ−13sinϕsinθ,θ≤θT,where
(9)A¯=13cosθT3+tanθTtan3θT+13signθtan3θT−3tanθTsinϕ,B¯=13cos3θTsignθsinθT+13sinϕcosθT,signθ=+1,θ≥0°,−1,θ<0°.Figure 2
Smooth processing of modified Mohr-Coulomb criterion in the octahedral plane.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:(10)G=σmsinφ+σ¯K′2+m2c2cos2φ,where φ is the dilatancy angle. Similarly, K′ is a function of Lode angle θ and dilatancy angle φ. Associated flow occurs if φ=ϕ, while unassociated flow occurs if 0≤φ<ϕ.


## 3. Coupled Elastoplastic Damage Model

### 3.1. Model Description

Figure3 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 stress-strain 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.Figure 3
The stress-strain curve in different stages.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:
(11)σ~=σ1−Ω,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(12)F=σ~msinϕ+σ¯~K′2+m2c2cos2ϕ−ccosϕ.By substituting Equation (11) into Equation (6), the yield function in terms of nominal stress can be determined by
(13)F=σmsinϕ+σ¯K′2+1−Ω2m2c2cos2ϕ−1−Ωccosϕ,Similarly, the potential function can be defined in terms of nominal stress as(14)G=σmsinφ+σ¯K′2+1−Ω2m2c2cos2φ.


### 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 strain-hardening stage and the total strain dεij is decomposed into the following two parts:
(15)dεij=dεijR=dεije+dεijed,where dεijR is the reversible strain, dεije is the elastic strain, and dεijedis the elastic damage strain.The elastic damage straindεijed is deduced from the change of elastic parameters during the strain hardening [17]:
(16)dεijed=∂CijklΩ∂ΩdΩ⋅σkl,where Cijkl 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 straindεijIR is decomposed into plastic and damage parts:
(17)dεijIR=dεijp+dεijd,where dεijp is the plastic strain and dεijd is the plastic damage strain.According to the potential function in Equation (14), the irreversible strain dεijIR can be defined as follows [8, 10]:
(18)dεijIR=dλ∂Gσij,Ω∂σij,where dλ is the plastic multiplier.The general stress and strain relationship can be expressed as follows:(19)dσij=Cijkldεkl−dεklIR+εkl−εklIR∂Cijkl∂ΩdΩ.


### 3.3. Definition of Damage

According to stress-strain relationship of clayey rock as shown in Figure3, 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
(20)Ωe=β1e¯−e¯0e,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]
(21)e¯=εijCijkl0εkl,where Cijkl0 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:(22)Ωp=e¯−e¯0pα2+β2e¯−e¯0p,where e¯0p is the energy factor corresponding to the plastic damage starting point and α2β2 are the plastic damage parameters.


### 3.4. Evolution of Model Parameters

According to the continuum damage mechanics, the degradation of elastic modulus during loading is defined as follows [4, 29]:
(23)E=1−ΩE0,where E is the damaged elastic modulus, E0 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:
(24)E=E0−E0−ErΩΩlim,0≤Ω≤Ωlim,E=Er,Ωlim≤Ω≤1,where Er 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:
(25)c=c0−c0−cr⋅Ωpη,where c0 and cr denote the initial and residual cohesions, respectively; η∈0,1 is a model parameter controlling the slope of the strain-softening curve.


## 3.1. Model Description

Figure3 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 stress-strain 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.Figure 3
The stress-strain curve in different stages.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:
(11)σ~=σ1−Ω,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(12)F=σ~msinϕ+σ¯~K′2+m2c2cos2ϕ−ccosϕ.By substituting Equation (11) into Equation (6), the yield function in terms of nominal stress can be determined by
(13)F=σmsinϕ+σ¯K′2+1−Ω2m2c2cos2ϕ−1−Ωccosϕ,Similarly, the potential function can be defined in terms of nominal stress as(14)G=σmsinφ+σ¯K′2+1−Ω2m2c2cos2φ.


## 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 strain-hardening stage and the total strain dεij is decomposed into the following two parts:
(15)dεij=dεijR=dεije+dεijed,where dεijR is the reversible strain, dεije is the elastic strain, and dεijedis the elastic damage strain.The elastic damage straindεijed is deduced from the change of elastic parameters during the strain hardening [17]:
(16)dεijed=∂CijklΩ∂ΩdΩ⋅σkl,where Cijkl 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 straindεijIR is decomposed into plastic and damage parts:
(17)dεijIR=dεijp+dεijd,where dεijp is the plastic strain and dεijd is the plastic damage strain.According to the potential function in Equation (14), the irreversible strain dεijIR can be defined as follows [8, 10]:
(18)dεijIR=dλ∂Gσij,Ω∂σij,where dλ is the plastic multiplier.The general stress and strain relationship can be expressed as follows:(19)dσij=Cijkldεkl−dεklIR+εkl−εklIR∂Cijkl∂ΩdΩ.


## 3.3. Definition of Damage

According to stress-strain relationship of clayey rock as shown in Figure3, 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
(20)Ωe=β1e¯−e¯0e,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]
(21)e¯=εijCijkl0εkl,where Cijkl0 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:(22)Ωp=e¯−e¯0pα2+β2e¯−e¯0p,where e¯0p is the energy factor corresponding to the plastic damage starting point and α2β2 are the plastic damage parameters.


## 3.4. Evolution of Model Parameters

According to the continuum damage mechanics, the degradation of elastic modulus during loading is defined as follows [4, 29]:
(23)E=1−ΩE0,where E is the damaged elastic modulus, E0 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:
(24)E=E0−E0−ErΩΩlim,0≤Ω≤Ωlim,E=Er,Ωlim≤Ω≤1,where Er 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:
(25)c=c0−c0−cr⋅Ωpη,where c0 and cr denote the initial and residual cohesions, respectively; η∈0,1 is a model parameter controlling the slope of the strain-softening curve.


## 4. Numerical Implementation

### 4.1. Implicit Constitutive Integration Algorithm

#### 4.1.1. Backward Euler Implicit Integration (BEII) Algorithm

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–33]. Given the responses at time tn, 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 tn+1. The BEII method includes mainly three steps, i.e., elastic predictor, plastic corrector, and damage corrector.In the elastic predictor step, the scalar damageΩn is assumed to be a constant. A trial damage is defined as
(26)Ωn+1trial=Ωn,and a trial stress σn+1trial is defined based on the elastic predictor
(27)σn+1trial=σn+CΔεn+1,where C is the elastic stiffness matrix, which is a function of the trial damage Ωn+1trial.The onset of plastic flow and deformation is determined by the loading condition, which is expressed as follows:(28)Fσn+1trial,Ωn+1trial=0,dλ≥0,Fσn+1trial,Ωn+1trial⋅dλ=0.The trial stress at point B,σBtrial, is obtained by the elastic predictor (as shown in Figure 4)
(29)σBtrial=σn+CΔεn+1=σn+Δσe,where Δσeis the increment of elastic stress.Figure 4
Sketch map of the backward Euler integration algorithm.In order to return the trial stressσn+1trial for predicting the yield surface, the stress increment is defined as
(30)Δσn+1=CΔεn+1−Δεn+1IR=CΔεn+1−ΔλCb=Δσe−ΔλCb,where b=∂Gσn+1trial,Ωn+1trial/∂σ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(31)σCtrial=σBtrial−ΔλCb.The first-order Taylor expansion of the yield function at point B gives(32)F=FB+∂F∂σTΔσ+∂F∂c∂c∂ΩpΔΩp+∂F∂ΩΔΩ=0.In the plastic corrector step, the scalar damageΩn+1trial is regarded as constant. Then, Equation (32) reduces to be
(33)F=FB+∂F∂σTΔσ=FB−ΔλaBTCbB=0,where a=∂F/∂σ and FB=FσBtrial,Ωn+1trial.Here the plastic multiplierΔλ can be determined explicitly as
(34)Δλ=FBaBTCbB.The backward Euler algorithm is based on the equation(35)σC=σB−ΔλCbC.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
(36)σC=σB−ΔλCbB.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.,
(37)r=σ−σB−ΔλCbC=σ−σB+ΔλCbC.The iterations continue until the norm of vectorr 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 rN is expressed as follows:
(38)rN=r0+σ̇+λ̇Cb+ΔλC∂b∂σσ̇,where r0 is the residual vector at point B, σ̇ is the change in σ, λ̇ is the change of Δλ, and N is the number of iteration.SettingrN to zero gives
(39)σ̇=−I+ΔλC∂b∂σ−1r0+λ̇Cb=−Q−1r0−λ̇Q−1Cb,where I is the unit matrix; Q=I+ΔλC∂b/∂σ.Similarly, a truncated Taylor series on the yield function is(40)FCN=FC0+∂F∂σTσ̇=0.By substituting Equation (39) into Equation (40), the change of Δλ is determined by
(41)λ̇=FC−aCTQ−1r0aCTQ−1Cb,and the plastic multiplier after N times of iteration is obtained:
(42)ΔλN=ΔλN−1+λ̇N−1.The modified stress vectorσC after N times of iterations is
(43)σCN=σB−ΔλN−1CbCN−1.The strain vectorεn+1 at time tn+1 is expressed as
(44)εn+1=εn+Δεn+1.The third step in the BEII method is the damage corrector. The updated damage at timetn+1 is
(45)Ωn+1=Ωn+ΔΩ,where ΔΩ=ΔΩe+ΔΩp.After the damage corrector is completed, the final stress at timetn+1 is determined:
(46)σn+1=Cn+1:εn+1e=Cn+1:εne+Δεn+1−Δεn+1IR=Cn+1Cnσn+1trial−Cn:Δεn+1IR=1−Ωn+11−ΩnσCN.


#### 4.1.2. Consistent Tangent Stiffness Matrix

The consistent tangent stiffness matrix is related to the convergence speed of global equilibrium iteration, and it does not influence the final results of stress updating. By dropping the subscript C in Equation (35), the standard back-Euler algorithm is expressed as
(47)σ=σB−ΔλCb.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
(48)σ̇=Cε̇−λ̇Cb−ΔλC∂b∂σσ̇,which is simplified to be
(49)σ̇=I+ΔλC∂b∂σ−1Cε̇−λ̇b=Q−1Cε̇−λ̇b=Rε̇−λ̇b,where R=Q−1C.To retain the current stressσ on the yield surface, Ḟ should be zero. The consistency conditions of yield function F is expressed as follows:
(50)Ḟ=∂F∂σTσ̇+∂F∂c∂c∂ΩpΩ̇p+∂F∂ΩΩ̇=aTσ̇=0.By substituting Equation (49) into Equation (50), λ̇ is determined by
(51)λ̇=aTRε̇aTRb.The consistent tangent stiffness matrix,Dct, can be obtained by substituting λ̇ into Equation (49):
(52)σ̇=R−RbaTRTaTRbε̇=Dctε̇.


### 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
(53)εn+1=εn+Δε,εn+1IR=εnIR+Δλn+1bn+1,Ωn+1=Ωn+ΔΩ,σn+1=C:εn+1−εn+1IR,Fn+1=Fn+1σn+1,Ωn+1.Equation (53) is a system of nonlinear equations of εn+1,εn+1IR, and Ωn+1. At any time tn, a set of values εn,εnIR,Ω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 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.(54)k=0:εn+1IR0=εnIR,Ωn+10=Ωn,Δλn+10=0,σn+10=C:εn+1−εn+1IR0.Step 2.
Check the yield condition and convergence at the iteration numberk.
(55)Fn+1k=Fσn+1k,Ωn+1k.
For a given stress error tolerance TOL, ifFn+1k<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.(56)δλn+1k=Fn+1k−an+1kTQ−1r0an+1kTQ−1Cb,δσn+1k=−I+ΔλC∂b∂σ−1r0+λ̇Cb,δΩn+1k=δΩen+1k+Ωpn+1k.Step 4.
Update plastic multiplier, strain, damage, and stress.(57)Δλn+1k+1=Δλn+1k+δλn+1k,Ωn+1k+1=Ωn+1k+δΩn+1k,εn+1IRk+1=εn+1IRk+C−1:δσn+1k,σn+1k+1=σn+1k+C:εn+1−εn+1IRk+1.
Letk=k+1; go to Step 3.


### 4.3. 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.Figure 5
Flow chart of the UMAT subroutine.


## 4.1. Implicit Constitutive Integration Algorithm

### 4.1.1. Backward Euler Implicit Integration (BEII) Algorithm

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–33]. Given the responses at time tn, 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 tn+1. The BEII method includes mainly three steps, i.e., elastic predictor, plastic corrector, and damage corrector.In the elastic predictor step, the scalar damageΩn is assumed to be a constant. A trial damage is defined as
(26)Ωn+1trial=Ωn,and a trial stress σn+1trial is defined based on the elastic predictor
(27)σn+1trial=σn+CΔεn+1,where C is the elastic stiffness matrix, which is a function of the trial damage Ωn+1trial.The onset of plastic flow and deformation is determined by the loading condition, which is expressed as follows:(28)Fσn+1trial,Ωn+1trial=0,dλ≥0,Fσn+1trial,Ωn+1trial⋅dλ=0.The trial stress at point B,σBtrial, is obtained by the elastic predictor (as shown in Figure 4)
(29)σBtrial=σn+CΔεn+1=σn+Δσe,where Δσeis the increment of elastic stress.Figure 4
Sketch map of the backward Euler integration algorithm.In order to return the trial stressσn+1trial for predicting the yield surface, the stress increment is defined as
(30)Δσn+1=CΔεn+1−Δεn+1IR=CΔεn+1−ΔλCb=Δσe−ΔλCb,where b=∂Gσn+1trial,Ωn+1trial/∂σ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(31)σCtrial=σBtrial−ΔλCb.The first-order Taylor expansion of the yield function at point B gives(32)F=FB+∂F∂σTΔσ+∂F∂c∂c∂ΩpΔΩp+∂F∂ΩΔΩ=0.In the plastic corrector step, the scalar damageΩn+1trial is regarded as constant. Then, Equation (32) reduces to be
(33)F=FB+∂F∂σTΔσ=FB−ΔλaBTCbB=0,where a=∂F/∂σ and FB=FσBtrial,Ωn+1trial.Here the plastic multiplierΔλ can be determined explicitly as
(34)Δλ=FBaBTCbB.The backward Euler algorithm is based on the equation(35)σC=σB−ΔλCbC.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
(36)σC=σB−ΔλCbB.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.,
(37)r=σ−σB−ΔλCbC=σ−σB+ΔλCbC.The iterations continue until the norm of vectorr 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 rN is expressed as follows:
(38)rN=r0+σ̇+λ̇Cb+ΔλC∂b∂σσ̇,where r0 is the residual vector at point B, σ̇ is the change in σ, λ̇ is the change of Δλ, and N is the number of iteration.SettingrN to zero gives
(39)σ̇=−I+ΔλC∂b∂σ−1r0+λ̇Cb=−Q−1r0−λ̇Q−1Cb,where I is the unit matrix; Q=I+ΔλC∂b/∂σ.Similarly, a truncated Taylor series on the yield function is(40)FCN=FC0+∂F∂σTσ̇=0.By substituting Equation (39) into Equation (40), the change of Δλ is determined by
(41)λ̇=FC−aCTQ−1r0aCTQ−1Cb,and the plastic multiplier after N times of iteration is obtained:
(42)ΔλN=ΔλN−1+λ̇N−1.The modified stress vectorσC after N times of iterations is
(43)σCN=σB−ΔλN−1CbCN−1.The strain vectorεn+1 at time tn+1 is expressed as
(44)εn+1=εn+Δεn+1.The third step in the BEII method is the damage corrector. The updated damage at timetn+1 is
(45)Ωn+1=Ωn+ΔΩ,where ΔΩ=ΔΩe+ΔΩp.After the damage corrector is completed, the final stress at timetn+1 is determined:
(46)σn+1=Cn+1:εn+1e=Cn+1:εne+Δεn+1−Δεn+1IR=Cn+1Cnσn+1trial−Cn:Δεn+1IR=1−Ωn+11−ΩnσCN.


### 4.1.2. Consistent Tangent Stiffness Matrix

The consistent tangent stiffness matrix is related to the convergence speed of global equilibrium iteration, and it does not influence the final results of stress updating. By dropping the subscript C in Equation (35), the standard back-Euler algorithm is expressed as
(47)σ=σB−ΔλCb.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
(48)σ̇=Cε̇−λ̇Cb−ΔλC∂b∂σσ̇,which is simplified to be
(49)σ̇=I+ΔλC∂b∂σ−1Cε̇−λ̇b=Q−1Cε̇−λ̇b=Rε̇−λ̇b,where R=Q−1C.To retain the current stressσ on the yield surface, Ḟ should be zero. The consistency conditions of yield function F is expressed as follows:
(50)Ḟ=∂F∂σTσ̇+∂F∂c∂c∂ΩpΩ̇p+∂F∂ΩΩ̇=aTσ̇=0.By substituting Equation (49) into Equation (50), λ̇ is determined by
(51)λ̇=aTRε̇aTRb.The consistent tangent stiffness matrix,Dct, can be obtained by substituting λ̇ into Equation (49):
(52)σ̇=R−RbaTRTaTRbε̇=Dctε̇.


## 4.1.1. Backward Euler Implicit Integration (BEII) Algorithm

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–33]. Given the responses at time tn, 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 tn+1. The BEII method includes mainly three steps, i.e., elastic predictor, plastic corrector, and damage corrector.In the elastic predictor step, the scalar damageΩn is assumed to be a constant. A trial damage is defined as
(26)Ωn+1trial=Ωn,and a trial stress σn+1trial is defined based on the elastic predictor
(27)σn+1trial=σn+CΔεn+1,where C is the elastic stiffness matrix, which is a function of the trial damage Ωn+1trial.The onset of plastic flow and deformation is determined by the loading condition, which is expressed as follows:(28)Fσn+1trial,Ωn+1trial=0,dλ≥0,Fσn+1trial,Ωn+1trial⋅dλ=0.The trial stress at point B,σBtrial, is obtained by the elastic predictor (as shown in Figure 4)
(29)σBtrial=σn+CΔεn+1=σn+Δσe,where Δσeis the increment of elastic stress.Figure 4
Sketch map of the backward Euler integration algorithm.In order to return the trial stressσn+1trial for predicting the yield surface, the stress increment is defined as
(30)Δσn+1=CΔεn+1−Δεn+1IR=CΔεn+1−ΔλCb=Δσe−ΔλCb,where b=∂Gσn+1trial,Ωn+1trial/∂σ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(31)σCtrial=σBtrial−ΔλCb.The first-order Taylor expansion of the yield function at point B gives(32)F=FB+∂F∂σTΔσ+∂F∂c∂c∂ΩpΔΩp+∂F∂ΩΔΩ=0.In the plastic corrector step, the scalar damageΩn+1trial is regarded as constant. Then, Equation (32) reduces to be
(33)F=FB+∂F∂σTΔσ=FB−ΔλaBTCbB=0,where a=∂F/∂σ and FB=FσBtrial,Ωn+1trial.Here the plastic multiplierΔλ can be determined explicitly as
(34)Δλ=FBaBTCbB.The backward Euler algorithm is based on the equation(35)σC=σB−ΔλCbC.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
(36)σC=σB−ΔλCbB.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.,
(37)r=σ−σB−ΔλCbC=σ−σB+ΔλCbC.The iterations continue until the norm of vectorr 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 rN is expressed as follows:
(38)rN=r0+σ̇+λ̇Cb+ΔλC∂b∂σσ̇,where r0 is the residual vector at point B, σ̇ is the change in σ, λ̇ is the change of Δλ, and N is the number of iteration.SettingrN to zero gives
(39)σ̇=−I+ΔλC∂b∂σ−1r0+λ̇Cb=−Q−1r0−λ̇Q−1Cb,where I is the unit matrix; Q=I+ΔλC∂b/∂σ.Similarly, a truncated Taylor series on the yield function is(40)FCN=FC0+∂F∂σTσ̇=0.By substituting Equation (39) into Equation (40), the change of Δλ is determined by
(41)λ̇=FC−aCTQ−1r0aCTQ−1Cb,and the plastic multiplier after N times of iteration is obtained:
(42)ΔλN=ΔλN−1+λ̇N−1.The modified stress vectorσC after N times of iterations is
(43)σCN=σB−ΔλN−1CbCN−1.The strain vectorεn+1 at time tn+1 is expressed as
(44)εn+1=εn+Δεn+1.The third step in the BEII method is the damage corrector. The updated damage at timetn+1 is
(45)Ωn+1=Ωn+ΔΩ,where ΔΩ=ΔΩe+ΔΩp.After the damage corrector is completed, the final stress at timetn+1 is determined:
(46)σn+1=Cn+1:εn+1e=Cn+1:εne+Δεn+1−Δεn+1IR=Cn+1Cnσn+1trial−Cn:Δεn+1IR=1−Ωn+11−ΩnσCN.


## 4.1.2. Consistent Tangent Stiffness Matrix

The consistent tangent stiffness matrix is related to the convergence speed of global equilibrium iteration, and it does not influence the final results of stress updating. By dropping the subscript C in Equation (35), the standard back-Euler algorithm is expressed as
(47)σ=σB−ΔλCb.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
(48)σ̇=Cε̇−λ̇Cb−ΔλC∂b∂σσ̇,which is simplified to be
(49)σ̇=I+ΔλC∂b∂σ−1Cε̇−λ̇b=Q−1Cε̇−λ̇b=Rε̇−λ̇b,where R=Q−1C.To retain the current stressσ on the yield surface, Ḟ should be zero. The consistency conditions of yield function F is expressed as follows:
(50)Ḟ=∂F∂σTσ̇+∂F∂c∂c∂ΩpΩ̇p+∂F∂ΩΩ̇=aTσ̇=0.By substituting Equation (49) into Equation (50), λ̇ is determined by
(51)λ̇=aTRε̇aTRb.The consistent tangent stiffness matrix,Dct, can be obtained by substituting λ̇ into Equation (49):
(52)σ̇=R−RbaTRTaTRbε̇=Dctε̇.


## 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
(53)εn+1=εn+Δε,εn+1IR=εnIR+Δλn+1bn+1,Ωn+1=Ωn+ΔΩ,σn+1=C:εn+1−εn+1IR,Fn+1=Fn+1σn+1,Ωn+1.Equation (53) is a system of nonlinear equations of εn+1,εn+1IR, and Ωn+1. At any time tn, a set of values εn,εnIR,Ω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 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.(54)k=0:εn+1IR0=εnIR,Ωn+10=Ωn,Δλn+10=0,σn+10=C:εn+1−εn+1IR0.Step 2.
Check the yield condition and convergence at the iteration numberk.
(55)Fn+1k=Fσn+1k,Ωn+1k.
For a given stress error tolerance TOL, ifFn+1k<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.(56)δλn+1k=Fn+1k−an+1kTQ−1r0an+1kTQ−1Cb,δσn+1k=−I+ΔλC∂b∂σ−1r0+λ̇Cb,δΩn+1k=δΩen+1k+Ωpn+1k.Step 4.
Update plastic multiplier, strain, damage, and stress.(57)Δλn+1k+1=Δλn+1k+δλn+1k,Ωn+1k+1=Ωn+1k+δΩn+1k,εn+1IRk+1=εn+1IRk+C−1:δσn+1k,σn+1k+1=σn+1k+C:εn+1−εn+1IRk+1.
Letk=k+1; go to Step 3.


## 4.3. 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.Figure 5
Flow chart of the UMAT subroutine.


## 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.Figure6 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.Figure 6
Geometrical model for triaxial tests.Figures7 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.Figure 7
Stress-strain relationship under uniaxial compression: (a) nonassociated flow and (b) associated flow.
(a)
(b)Figure 8
Stress-strain relationship under triaxial compression with a confining pressure 3 MPa: (a) nonassociated cases and (b) associated flow.
(a)
(b)The mechanical response of the sample under uniaxial tension is shown in Figure9 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 the built-in model in ABAQUS. In this case, the calculated tension strength from the MMC model is 0.424 MPa, which is very close to the calculated value of 0.436 MPa from the built-in model. When m=0.2, the tension strength obtained by the MMC model is 0.353 MPa, and it is significantly lower than that predicted by the built-in model. Therefore, the actual tension strength of rock is dependent on the parameter m.Figure 9
Stress-strain relation in uniaxial tensile cases.


### 5.2. 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 m2), 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 Figure6. 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]:
(58)Δt≥γw6E′kΔh2,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 Tables1 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.Table 1
Initial damage points under different confining pressures.
Confining pressure (MPa)
e¯0e (MPa1/2)
e¯0p(MPa1/2)
0.89
0.273
0.354
2.50
0.227
0.413
2.85
0.197
0.322
5.42
0.191
0.553Table 2
The known parameters for clayey rock.
Parameters
E0 (MPa)
Poisson’s ratioμ
c0 (MPa)
ϕ (°)
k (m2)
Value
300
0.13
0.30
18
3×10−19Table 3
The results of unknown parameters by back analysis.
Confining pressure
β1
α2
β2
cr (KPa)
η
Er (MPa)
φ (°)
0.89
1.381
0.569
0.990
58.132
0.608
77.418
7.290
2.50
0.860
0.583
1.243
97.946
0.404
95.593
2.501
2.85
0.939
1.190
1.886
44.399
0.693
66.270
0.804
5.42
0.590
2.008
0.951
49.008
0.699
117.196
0.921The stress-strain curves under varying confining pressures from 0.89 to 5.42 MPa are shown in Figure10. 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 strain-softening 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.Figure 10
Comparison of the stress-strain curves under different confining pressures between the numerical analysis and experiments: (a) confining pressure 0.89 MPa, (b) confining pressure2.5 MPa, (b) confining pressure 2.85 MPa, and (d) confining pressure 5.42 MPa.
(a)
(b)
(c)
(d)


## 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.Figure6 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.Figure 6
Geometrical model for triaxial tests.Figures7 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.Figure 7
Stress-strain relationship under uniaxial compression: (a) nonassociated flow and (b) associated flow.
(a)
(b)Figure 8
Stress-strain relationship under triaxial compression with a confining pressure 3 MPa: (a) nonassociated cases and (b) associated flow.
(a)
(b)The mechanical response of the sample under uniaxial tension is shown in Figure9 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 the built-in model in ABAQUS. In this case, the calculated tension strength from the MMC model is 0.424 MPa, which is very close to the calculated value of 0.436 MPa from the built-in model. When m=0.2, the tension strength obtained by the MMC model is 0.353 MPa, and it is significantly lower than that predicted by the built-in model. Therefore, the actual tension strength of rock is dependent on the parameter m.Figure 9
Stress-strain relation in uniaxial tensile cases.


## 5.2. 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 m2), 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 Figure6. 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]:
(58)Δt≥γw6E′kΔh2,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 Tables1 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.Table 1
Initial damage points under different confining pressures.
Confining pressure (MPa)
e¯0e (MPa1/2)
e¯0p(MPa1/2)
0.89
0.273
0.354
2.50
0.227
0.413
2.85
0.197
0.322
5.42
0.191
0.553Table 2
The known parameters for clayey rock.
Parameters
E0 (MPa)
Poisson’s ratioμ
c0 (MPa)
ϕ (°)
k (m2)
Value
300
0.13
0.30
18
3×10−19Table 3
The results of unknown parameters by back analysis.
Confining pressure
β1
α2
β2
cr (KPa)
η
Er (MPa)
φ (°)
0.89
1.381
0.569
0.990
58.132
0.608
77.418
7.290
2.50
0.860
0.583
1.243
97.946
0.404
95.593
2.501
2.85
0.939
1.190
1.886
44.399
0.693
66.270
0.804
5.42
0.590
2.008
0.951
49.008
0.699
117.196
0.921The stress-strain curves under varying confining pressures from 0.89 to 5.42 MPa are shown in Figure10. 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 strain-softening 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.Figure 10
Comparison of the stress-strain curves under different confining pressures between the numerical analysis and experiments: (a) confining pressure 0.89 MPa, (b) confining pressure2.5 MPa, (b) confining pressure 2.85 MPa, and (d) confining pressure 5.42 MPa.
(a)
(b)
(c)
(d)


## 6. 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.



---

*Source: 9853782-2020-03-20.xml*