# A Damaged Constitutive Model for Rock under Dynamic and High Stress State

Authors: Yan-Long Li; Zong-Yuan Ma
Journal: Shock and Vibration
 (2017)

---

## Abstract

The main research work of this paper focuses on the theoretical prediction of the constitutive relationship for rock, concrete, and other quasi-brittle materials under dynamic and complex stress state and the influence of dynamic mechanical behavior of rock on practical engineering problems was studied. A damaged elastoplastic model (DEPM) is established for the investigation and prediction of static or dynamic mechanical behavior of rock material. The mechanical behavior (brittleness or plasticity) and dynamic response (due to underground impact pressure and high-velocity impact of projectile) of rock under high in situ stress were investigated via the DEPM combined with the explicit finite element method. This paper suggests the influence of the brittle or plastic mechanical behavior of rock material on deep underground rock engineering.


---

## Body

## 1. Introduction

The strength of the rock or concrete is studied by many researchers for a long period. The strength of rock or concrete (quasi-brittle materials) is sensitive to the hydrostatic pressure(σm=(σ1+σ2+σ3)/3) under the high stress state [1, 2]. The high in situ stress and the quality of rock mass are also the key factors of deep underground rock engineering [3–5]. However, the regularity of brittle and plastic mechanical behavior for intact rock or concrete material under the high confining pressure has seldom been investigated. The high confining pressure σ3 has significant influence on the mechanical behavior (brittleness or plasticity) of intact rock or concrete material [6–8]. The high confining pressure σ3 also has significant influence on the stability and failure mechanism of deep underground rock engineering (Cai 2008). The constitutive relationship for quasi-brittle materials subjected to large strains, high strain rates, and high pressures was investigated by some researchers [9–13], but the brittle or plastic mechanical behavior for rock or concrete under high confining pressure has still not been clearly defined. At present, the theoretical methods and experimental techniques are still limited for the research of dynamic mechanical behavior of deep underground rock engineering [14–16]. The main research work of this study focuses on the theoretical prediction of the constitutive relationship for rock under high confining pressure condition, and the influence of dynamic mechanical behavior of rock on deep-buried engineering problems will be investigated.


## 2. Constitutive Modeling

The general objective of this study is to establish a damaged elastoplastic model (DEPM) for the investigation and prediction of static or dynamic mechanical behavior of intact rock materials. The plastic and damage behavior of rock under high in situ stress and dynamic response (due to earthquake, underground impact pressure, and high-velocity impact) can be investigated via this model combined with the explicit finite element method.


### 2.1. Damage Part for DEPM

The stress-strain relations of damaged mechanics are governed by the scalar damage which is expressed as follows [17, 18]:(1)σ=1-d·σ-=1-dD0e:ε-εp,where D0e is the initial (undamaged) elastic stiffness of the material; 1-dD0e is the degraded elastic stiffness; and d is the scalar stiffness degradation variable, which can take values in the range from zero (undamaged material) to one (fully damaged material). Damage associated with the failure mechanisms of the rock or concrete (cracking and crushing) therefore results in a reduction in the elastic stiffness. Within the context of the scalar damage, the stiffness degradation is isotropic and characterized by the degradation variable d. Following the usual notions of continuum damage mechanics, the Cauchy stress is related to the effective stress through the scalar degradation relation. The damage evolution law of damage variable d is shown as follows [19]:(2)d=0εmax<ε01-η-λη-1ε0εmax+1-λη-1ε0<εmax≤εr1-λε0εmaxεr<εmax≤εu1εmax>εu,where ε0 is the initial damage strain (micro crack initiation, i.e., damage threshold), ε0=σt/E, or ε0=σc/E (σt and σc are the uniaxial tension and compression peck strength, and E is Young’s modulus for rock), εmax is the maximum strain during the deformation of rock material, and it is evaluated from the incremental form of the maximum principal strain ε1 in tension and minimum principal strain ε3 in compression. The maximum strain εmax cannot be decreased. The relationship between damage variable d and maximum strain εmax for damage evolution is shown in Figure 1, where σ0 is the peck values of the stress and strain relationship yield by the uniaxial tension or compression test and σr is the residual values of the stress and strain relationship yield by the uniaxial tension or compression test. εr is the residual strain corresponding to the residual values of the stress and strain curve yield by the uniaxial tension or compression test. εu is the ultimate strain corresponding to the fully failure of rock material. The relationship of damage evolution is bilinear and divided by three stages: firstly, εmax<ε0 and the damage variable d is equal to zero; secondly, ε0<εmax<εr is the first damage stage, and the brittle failure behavior of rock material occurred; εr<εmax<εu is the second damage stage; a residual strength and deformation of rock occurred; if εmax>εu, the damage variable d is equal to one, and the rock material is fully failure. λ is the coefficient of residual strength (σr=λσ0), η is the coefficient of residual strain (εr=ηε0), and ξ is the coefficient of ultimate strain (εu=ξε0).Figure 1
Damage evolution relationship between damage variabled and maximum strain εmax.A stress weight factor is defined that is equal to one if all principal stressesσ-i (i=1,2,3) are positive and equal to zero if they are negative [20]. The Macaulay bracket · is defined by x=x+x/2.(3)rσ-=∑i=13σ-i∑i=13σ-i,0≤rσ-≤1.The scalar damage degradation variable d must be consistent with the damage variable in tension dt or in compression dc. The damage degradation variable d is assumed as follows: (4)d=1-1-rσ-dt1-1-rσ-dc.The values of the initial damage strain ε0 are increased with the values of confining pressure σ3 increased, and the relationship for initial ε0 and σ3 is expressed as follows:(5)ε0=ε0′hσ3σc+1,where ε0′ is initial damage strain obtained by uniaxial tension or compression test.


### 2.2. Plastic Part for DEPM

For plastic part, the shear strength envelope for rock material follows the Hoek–Brown criterion [1].(6)fs=σ1-σ3-σcimbσ3σci+sn=0,where fs yield functions of Hoek–Brown criterion, σci is the uniaxial compressive strength of intact rock material, and mb, n, and s are the material parameter. The order of the three principal stresses is as follows, σ1≥σ2≥σ3, and positive in compression. The tension strength σt of the rock material suggested by Hoek and Brown [2] is given by:(7)σt=σci2mb-mb2+4s.The volumetric strength envelope for rock material is as follows:(8)fV=I12-2I21/2-pc=0,where I1 and I2 are the first and second stress invariant, pc is the cap pressure, and the plastic volume changed is controlled by this volumetric yield function. The shape of the volumetric strength envelope is sphere. The common boundary between the shear strength envelope and volumetric strength envelope is pc=σcimbσ3/σci+sn. The plastic potential function g for shear loading surface can be written as follows:(9)gs=σ1-σ3-Hγspσcimbσ3σci+sn=0and the volumetric loading surface(10)gV=I12-2I21/2-Hεvppc=0,where H is the strain hardening function and γsp and εvp are the shear and volumetric hardening parameter. The strain hardening function for the shear and volumetric loading surface can be written as follows:(11)Hγsporεvp=1+σ3σc1-e-γsporεvpn.The plastic deformation and damage failure behavior of rock under high stress state can be predicated by this model. The DEPM constitutive model is written in FORTRAN as user subroutine to define material behavior (VUMAT) and loaded into the ABAQUS FEM code and computed by the explicit FEM scheme.The strain rate effect is reflected by the viscoplasticity in the constitutive model [21]. The viscoplastic plastic flow rule is as follows:(12)ε˙pl=λ˙∂gσ∂σ,where ε˙pl is the plastic strain rate, λ˙ is the nonnegative plastic multiplier, g is the plastic potential function and defined in the former contents, and associated plasticity is used in DEPM.


### 2.3. Verification

The plastic deformation and brittle failure behavior of Yamaguchi Marble under triaxial compression is simulated by the DEPM model via ABAQUS explicit scheme. The values of parameters in DEPM for triaxial test simulation are shown in Table1. Only one hexahedron element is used to simulate the triaxial test, and the maximum principal stress σ1 is loaded by displacement boundary condition, and confining pressure σ3 is stress boundary condition. The stress-strain relationships of Yamaguchi Marble predicted by DEPM model and suggestd by Mogi [7] are shown in Figure 2.Table 1
Values of parameters for numerical simulation.
Parameters for damage part Parameters for plastic part Young’s modulusE (GPa) Poisson’s ratiov Residual strength coefficientλ Residual strain coefficientη Ultimate strain coefficientξ σ c i (MPa) m b a s 40 0.25 0.1 10 20 80 5 0.5 1Figure 2
Stress-strain relationships of Yamaguchi Marble under triaxial compression by DEPM prediction (this study) and test [7].Another verification case is the drop hammer test and its numerical simulation. The drop hammer test was progressed in this study with a concrete beam in three-point bending condition. Three steel bars were placed in the position the same as the three-point bending test, and the hammer (1.0 kg) fell down from the height of 10 meters and hit the steel bar on the top of the concrete beam. The results obtained by FEM simulation using DEPM model and comparison of the drop hammer test photograph are shown in Figure3. The results of the drop hammer test suggest that the damage zone evolution and the crack propagation are similar between the test results and the numerical simulation.Figure 3
Concrete beam fracture due to impact load.
(a)
FEM simulation using DEPM (b)
Test results


### 2.4. Comparison with Similar Models

A damaged constitutive model for brittle materials subjected to large strains, high strain rates, and high pressures was suggested by Johnson and Holmquist [9, 10], and this model is implemented in LS-DYNA produced by Livermore Software Technology Corporation. Two yield surfaces (shear and volumetric) are used in Johnson and Holmquist’s model, and the hydrostatic pressure-volume response is taking into account during the impact or perforation process. A plastic-damage constitutive model for cyclic loading of concrete structure was suggested by Lee and Fenves [20], and this model is implemented in ABAQUS produced by SIMULIA Corporation. Only one yield surface (shear) is used in Lee and Fenves’s model, and two scalar degradation variables are used for axial tension and compression conditions, respectively. However, the dynamic mechanical behavior for brittle materials (rock or concrete) under high confining pressure has not been fully considered in such constitutive models. The dynamic mechanical behavior under high confining pressure condition can be taken into account and well predicted by the constitutive models DEPM suggested in this paper.


## 2.1. Damage Part for DEPM

The stress-strain relations of damaged mechanics are governed by the scalar damage which is expressed as follows [17, 18]:(1)σ=1-d·σ-=1-dD0e:ε-εp,where D0e is the initial (undamaged) elastic stiffness of the material; 1-dD0e is the degraded elastic stiffness; and d is the scalar stiffness degradation variable, which can take values in the range from zero (undamaged material) to one (fully damaged material). Damage associated with the failure mechanisms of the rock or concrete (cracking and crushing) therefore results in a reduction in the elastic stiffness. Within the context of the scalar damage, the stiffness degradation is isotropic and characterized by the degradation variable d. Following the usual notions of continuum damage mechanics, the Cauchy stress is related to the effective stress through the scalar degradation relation. The damage evolution law of damage variable d is shown as follows [19]:(2)d=0εmax<ε01-η-λη-1ε0εmax+1-λη-1ε0<εmax≤εr1-λε0εmaxεr<εmax≤εu1εmax>εu,where ε0 is the initial damage strain (micro crack initiation, i.e., damage threshold), ε0=σt/E, or ε0=σc/E (σt and σc are the uniaxial tension and compression peck strength, and E is Young’s modulus for rock), εmax is the maximum strain during the deformation of rock material, and it is evaluated from the incremental form of the maximum principal strain ε1 in tension and minimum principal strain ε3 in compression. The maximum strain εmax cannot be decreased. The relationship between damage variable d and maximum strain εmax for damage evolution is shown in Figure 1, where σ0 is the peck values of the stress and strain relationship yield by the uniaxial tension or compression test and σr is the residual values of the stress and strain relationship yield by the uniaxial tension or compression test. εr is the residual strain corresponding to the residual values of the stress and strain curve yield by the uniaxial tension or compression test. εu is the ultimate strain corresponding to the fully failure of rock material. The relationship of damage evolution is bilinear and divided by three stages: firstly, εmax<ε0 and the damage variable d is equal to zero; secondly, ε0<εmax<εr is the first damage stage, and the brittle failure behavior of rock material occurred; εr<εmax<εu is the second damage stage; a residual strength and deformation of rock occurred; if εmax>εu, the damage variable d is equal to one, and the rock material is fully failure. λ is the coefficient of residual strength (σr=λσ0), η is the coefficient of residual strain (εr=ηε0), and ξ is the coefficient of ultimate strain (εu=ξε0).Figure 1
Damage evolution relationship between damage variabled and maximum strain εmax.A stress weight factor is defined that is equal to one if all principal stressesσ-i (i=1,2,3) are positive and equal to zero if they are negative [20]. The Macaulay bracket · is defined by x=x+x/2.(3)rσ-=∑i=13σ-i∑i=13σ-i,0≤rσ-≤1.The scalar damage degradation variable d must be consistent with the damage variable in tension dt or in compression dc. The damage degradation variable d is assumed as follows: (4)d=1-1-rσ-dt1-1-rσ-dc.The values of the initial damage strain ε0 are increased with the values of confining pressure σ3 increased, and the relationship for initial ε0 and σ3 is expressed as follows:(5)ε0=ε0′hσ3σc+1,where ε0′ is initial damage strain obtained by uniaxial tension or compression test.


## 2.2. Plastic Part for DEPM

For plastic part, the shear strength envelope for rock material follows the Hoek–Brown criterion [1].(6)fs=σ1-σ3-σcimbσ3σci+sn=0,where fs yield functions of Hoek–Brown criterion, σci is the uniaxial compressive strength of intact rock material, and mb, n, and s are the material parameter. The order of the three principal stresses is as follows, σ1≥σ2≥σ3, and positive in compression. The tension strength σt of the rock material suggested by Hoek and Brown [2] is given by:(7)σt=σci2mb-mb2+4s.The volumetric strength envelope for rock material is as follows:(8)fV=I12-2I21/2-pc=0,where I1 and I2 are the first and second stress invariant, pc is the cap pressure, and the plastic volume changed is controlled by this volumetric yield function. The shape of the volumetric strength envelope is sphere. The common boundary between the shear strength envelope and volumetric strength envelope is pc=σcimbσ3/σci+sn. The plastic potential function g for shear loading surface can be written as follows:(9)gs=σ1-σ3-Hγspσcimbσ3σci+sn=0and the volumetric loading surface(10)gV=I12-2I21/2-Hεvppc=0,where H is the strain hardening function and γsp and εvp are the shear and volumetric hardening parameter. The strain hardening function for the shear and volumetric loading surface can be written as follows:(11)Hγsporεvp=1+σ3σc1-e-γsporεvpn.The plastic deformation and damage failure behavior of rock under high stress state can be predicated by this model. The DEPM constitutive model is written in FORTRAN as user subroutine to define material behavior (VUMAT) and loaded into the ABAQUS FEM code and computed by the explicit FEM scheme.The strain rate effect is reflected by the viscoplasticity in the constitutive model [21]. The viscoplastic plastic flow rule is as follows:(12)ε˙pl=λ˙∂gσ∂σ,where ε˙pl is the plastic strain rate, λ˙ is the nonnegative plastic multiplier, g is the plastic potential function and defined in the former contents, and associated plasticity is used in DEPM.


## 2.3. Verification

The plastic deformation and brittle failure behavior of Yamaguchi Marble under triaxial compression is simulated by the DEPM model via ABAQUS explicit scheme. The values of parameters in DEPM for triaxial test simulation are shown in Table1. Only one hexahedron element is used to simulate the triaxial test, and the maximum principal stress σ1 is loaded by displacement boundary condition, and confining pressure σ3 is stress boundary condition. The stress-strain relationships of Yamaguchi Marble predicted by DEPM model and suggestd by Mogi [7] are shown in Figure 2.Table 1
Values of parameters for numerical simulation.
Parameters for damage part Parameters for plastic part Young’s modulusE (GPa) Poisson’s ratiov Residual strength coefficientλ Residual strain coefficientη Ultimate strain coefficientξ σ c i (MPa) m b a s 40 0.25 0.1 10 20 80 5 0.5 1Figure 2
Stress-strain relationships of Yamaguchi Marble under triaxial compression by DEPM prediction (this study) and test [7].Another verification case is the drop hammer test and its numerical simulation. The drop hammer test was progressed in this study with a concrete beam in three-point bending condition. Three steel bars were placed in the position the same as the three-point bending test, and the hammer (1.0 kg) fell down from the height of 10 meters and hit the steel bar on the top of the concrete beam. The results obtained by FEM simulation using DEPM model and comparison of the drop hammer test photograph are shown in Figure3. The results of the drop hammer test suggest that the damage zone evolution and the crack propagation are similar between the test results and the numerical simulation.Figure 3
Concrete beam fracture due to impact load.
(a)
FEM simulation using DEPM (b)
Test results


## 2.4. Comparison with Similar Models

A damaged constitutive model for brittle materials subjected to large strains, high strain rates, and high pressures was suggested by Johnson and Holmquist [9, 10], and this model is implemented in LS-DYNA produced by Livermore Software Technology Corporation. Two yield surfaces (shear and volumetric) are used in Johnson and Holmquist’s model, and the hydrostatic pressure-volume response is taking into account during the impact or perforation process. A plastic-damage constitutive model for cyclic loading of concrete structure was suggested by Lee and Fenves [20], and this model is implemented in ABAQUS produced by SIMULIA Corporation. Only one yield surface (shear) is used in Lee and Fenves’s model, and two scalar degradation variables are used for axial tension and compression conditions, respectively. However, the dynamic mechanical behavior for brittle materials (rock or concrete) under high confining pressure has not been fully considered in such constitutive models. The dynamic mechanical behavior under high confining pressure condition can be taken into account and well predicted by the constitutive models DEPM suggested in this paper.


## 3. Application

The stability and failure mechanism of deep underground excavation influenced by the underground impact pressure was analyzed via the DEPM model and explicit finite element method. A brittle failure behavior was obtained by the uniaxial tension or compression stress state (σ2=σ3=0) for rock material. However, a large plastic deformation was obtained by the high confining pressure (σ1≠σ2=σ3) condition for rock material. Moreover, the deep underground excavation is often influenced by the underground impact pressure due to the high in situ stress condition. Two application cases were analyzed using DEPM model under the dynamic and high magnitude stress state. The values of the parameters for the computation of application cases are the same as the verification cases (shown in Table 1).The FEM model, damage zone, and dynamic failure mechanisms of surrounding rock for rectangular tunnel under different initial in situ stress are shown in Figure4. The initial in situ stress and underground impact pressure are uniform applied around the circle boundary of the surrounding rock. The calculation results shown that, after the impact loading, the surrounding rock in the corners and lateral walls of the tunnel was cracked (brittle failure mechanism) under low or no initial in situ stress condition. However, the large plastic deformation was yielded in the surrounding rock closed to the tunnel after the impact loading and high initial in situ stress.Figure 4
FEM model and failure mechanisms of surrounding rock for rectangular tunnel under different initial in situ stress.Due to the geologic conditions limit and low permeability, ultra-heavy oil is very difficult to exploitation. In general, the ultra-heavy oil reservoir must be perforated by the metallic shaped charge jet in the deeply buried rock stratum. The depth and effect of penetration of rock stratum will be influenced by the magnitude of the in situ stress significantly. The high-velocity impact of a rigid penetrator into the rock stratum under deep underground and high initial in situ condition was investigated via this model combined with the explicit finite element method. The FEM model and damage zone of the rock stratum by a rigid projectile perforation under different confining pressureσ3 are shown in Figure 5. The initial velocity of rigid projectile is about 1000 m/s. Large-scale plastic zone of the rock stratum was yielded around the head of the projectile under high confining stress condition, and the penetration depth of the projectile under high confining stress condition is smaller than the case of low or no confining stress condition significantly.Figure 5
High-velocity impact and perforation by a rigid projectile into a rock stratum under different confining pressure (initial velocity = 1000 m/s).


## 4. Conclusions

In comparison with similar models, the dynamic mechanical behavior for brittle materials (rock or concrete) under high confining pressure can be taken into account and well predicted by the constitutive models DEPM suggested in this paper. The strain rate effect is reflected by viscoplasticity in DEPM. Based on the results of this study, the following conclusions and recommendations for future work can be stated.(1)
In this paper, the results suggest that a great influence will be yielded on the deep underground rock engineering due to the different mechanical behavior (brittle or ductile) of rock material under the high pressure stress state. The research works of this study also aim to provide the regularity about the influences of mechanical behavior for rock material on the practical problems of rock engineering under the high pressure stress state and dynamic situation.(2)
The investigation and prediction of the static or dynamic behavior for rock material under the high pressure stress state will be very helpful to the deep underground excavation engineering, heavy oil exploitation, and so forth. The new insight or idea for the dynamic response and failure mechanism of deep underground rock engineering influenced by the high in situ stress, rock blast, or underground impact pressure may be obtained.(3)
In future studies, the stability and failure mechanisms of the surrounding rock for deep underground rock engineering due to the dynamic load will be further investigated using the constitutive model suggested in this paper.



---

*Source: 8329545-2017-02-23.xml*