# A Coupled Elastoplastic Damage Dynamic Model for Rock

Authors: Xuelong Hu; Ming Zhang; Xiangyang Zhang; Min Tu; Zhiqiang Yin; Haifeng Ma; Minke Duan
Journal: Shock and Vibration
 (2021)

---

## Abstract

Rock dynamic constitutive model plays an important role in understanding dynamic response and addressing rock dynamic problems. Based on elastoplastic mechanics and damage mechanics, a dynamic constitutive model of rock coupled with elastoplastic damage is established. In this model, unified strength theory is taken as the yield criterion; to reflect the different damage evolution law of rocks under tension and pressure conditions, the effective plastic strain and volumetric plastic strain are used to represent the compressive damage variable and the equivalent plastic strain is used to represent the tensile damage variable; the plastic hardening behavior and strain rate effect of rocks are characterized by piecewise function and dynamic increase factor function, respectively; Fortran language and LS-DYNA User-Defined Interface (Umat) are used to numerically implement the constitutive model; the constitutive model is verified by three classical examples of rock uniaxial and triaxial compression tests, rock uniaxial tensile test, and rock ballistic test. The results show that the constitutive model can describe the dynamic and static mechanical behavior of rock comprehensively.


---

## Body

## 1. Introduction

Rock is a brittle material widely involved in the fields of mining engineering, civil engineering, and protection engineering. In these fields, rock is often subjected to dynamic load such as explosion and impact, which will bear high pressure and high strain rate [1–4]. How to accurately describe and comprehensively predict the dynamic mechanical behavior of rock materials is always one of the common goals of scholars at home and abroad. Studies have shown [5–7] that an appropriate rock materials model is essential for understanding and predicting the dynamic response and failure behavior of rock materials.At present, with the joint efforts of scholars at home and abroad, a series of dynamic constitutive models of rock materials have been proposed, among which HJC model [8], KCC model [9], CSC model [10], and RHT model [11] are widely used, and these models are known as classical models. All the classical models have certain defects. For example, the HJC model only considers the compression damage and ignores the tensile damage [12]. The KCC model cannot accurately capture the tensile behavior of rock materials [13]. The CSC model calculates plastic strain through the associated flow law, which may lead to the failure to accurately capture the volumetric expansion behavior of rock materials [14]. The RHT model ignores the Lode angle effect [15].In order to overcome the defects of the classical model, some modified models and new models have been established. Polanco-Loria et al. [12] modified the HJC model by considering the sensitivity of the third deviatorial stress affecting the strain rate and the tensile damage variable. Ling et al. [15] proposed a bilinear tension-softening model based on the RHT model, introduced Lode angle factor, and modified the calculation formulas of strain rate enhancement factor and tension-compression meridian-to-midday ratio. Yang et al. [16] established a new rock damage-plastic constitutive model, which took into account the effects of compressive pressure and strain rate on rock strength and postpeak characteristics under dynamic loading. Li and Shi [17] established a dynamic model of rock materials based on the extended Drucker Prager strength criterion and Johnson Cook material model, which considered the mechanical behavior of rock materials under high confining pressures and high strain rates. Huang et al. [1] proposed a rock constitutive model based on the dynamic model of concrete established by Kong et al. [18], which fully considered the effect of high confining pressure and strain rate of rock. Hu et al. [19] established a rock dynamic damage model based on the unified strength theory. Both the classical model and the newly established model have their own advantages and disadvantages [1]. Therefore, the improvement of the classical model and the establishment of the new model are still one of our joint efforts.The unified strength theory is used as the yield criterion in this paper; on the basis of fully considering the strain rate effect of rock materials, different hardening behaviors and damage laws under tension and compression were treated separately; an elastoplastic damage coupling dynamic model is established; and the constitutive model is verified by uniaxial and triaxial compression tests, uniaxial tensile tests, and rock ballistic tests.


## 2. Establishment of Model

### 2.1. Basic Relationship

According to elastic mechanics, the relationship between stress and strain can be expressed as(1)σ=Dεe,where σ is the stress vector; D is the elastic stiffness matrix; and εe is the elastic strain vector.According to elastoplastic mechanics, the relationship between total strain (increment), elastic strain (increment), and plastic strain (increment) is(2)ε=εe+εp,(3)dε=dεe+dεp,where ε is the total strain vector; εp is the plastic strain vector; dε is the total strain increment vector; dεe is an elastic strain increment vector; and dεp is the plastic strain increment vector.According to damage mechanics, the relation between the elastic stiffness matrix of damaged materials and that of undamaged materials is(4)D=1−ωD0,where ω is the damage variable and D0 is the elastic stiffness matrix of the undamaged materials.Combining equations (1)–(4), the relationship between stress increment and strain increment can be obtained as(5)dσ=1−ωD0dε−dεp−D0ε−εpdω.


### 2.2. Strength Criterion

This paper chooses the unified strength theory as the plastic yield criterion of rock materials. The unified strength theory is established on the basis of the twin shear orthogonal octahedral element model, which stipulates that the material will fail when a certain function of the two larger principal shear stresses and their corresponding normal stress reaches a certain threshold [20]. The unified strength theory takes into account all stress components and their different effects on material failure, which can be applied to various geotechnical materials, and has been widely used in geotechnical engineering and other fields [21–23]. The unified strength theory is expressed in terms of stress invariant and Lode angle as(6)F=131−aI1+13J2a−b−2ab−11+bsinθL+J2a+b+11+bcosθL−κθL∈−π6,arctan31−a31+a131−aI1+13J2a+2b+ab−11+bsinθL+J2a+ab+11+bcosθL−κθL∈arctan31−a31+a,π6,where I1 is the invariant of the first principal stress tensor; J2 is the invariant of the second deviator stress tensor; θL is the Lode angle; a is the ratio of tensile and compression strength of the rock materials; b is a parameter reflecting the influence of intermediate principal stress on material failure; in the convexity theory, the value range of b is 0 ≤ b ≤ 1; and κ is the strength function. In order to reflect the strength degradation of the rock materials during the deformation process and the hardening behavior of the rock materials, the damage variable ω and the hardening function are introduced into the strength function, and then κ can be expressed as(7)κ=h1−ω2ccosφ1+sinφ,where h is the hardening function; φ is the internal friction angle of the rock materials; and c is the cohesion of the rock materials.Rock materials generally show hardening mechanical behavior under compressive load but do not have this mechanical behavior under tensile load. In order to distinguish this behavior of rock materials, this paper introduces the following function to express the hardening behavior of rock materials:(8)h=β0+βm−β0γpb1+γp,p≤0,1,p>0,where β0 and βm are the initial critical value and maximum critical value of plastic hardening, respectively; b1 is the parameter that controls the hardening rate; γp is the effective plastic strain; and p is the hydrostatic pressure. In this paper, in order to be consistent with elastic mechanics, the tensile stress is taken as positive. The effective plastic strain increment can be calculated by the following formula:(9)dγp=23dεijp−13δijdεvpdεijp−13δijdεvp,where dγp is the effective plastic strain increment; dεijp is the plastic strain increment; dεvp is the volumetric plastic strain increment; and δij is the Kronecker symbol.The unified strength theory is composed of a series of piecewise linear strength criteria on the deviating plane [20], and its changes are shown in Figure 1. The specific form of the unified strength theory expression depends on the choice of parameter b. When b takes different values, the unified strength theory can be simplified to many current strength criteria. For example, when b = 0 (0 <a < 1), the unified strength theory becomes the Mohr–Coulomb criterion; when b = 1 (0 < a < 1), the unified strength theory becomes the generalized twin shear stress criterion.Figure 1
Locus of unified strength theory in deviatoric plane.


### 2.3. Strain Rate Effect

Rock is a strain rate-dependent geological material, and the strain rate has a significant influence on the mechanical behavior of the rock [24, 25]. The most direct effect of strain rate on rock mechanical behavior is to increase the strength of rock. In rock dynamics, dynamic increase factor is generally used to characterize the influence of strain rate on rock strength [26], and dynamic increase factor can be defined as(10)DIF=σdσs,where DIF is the dynamic increase factor; σs is the quasi-static strength of rock materials; and σd is the dynamic strength of the rock materials.According to references [26–28], it can be seen that the strain rate has little effect on the friction strength of the rock and can be ignored; the strain rate has a greater influence on the cohesive strength of the rock, so the dynamic cohesive strength of rock can be expressed as a function of strain rate and quasi-static cohesive strength, namely,(11)cd=DIFc,where cd is the dynamic cohesion of rock materials. In this paper, the dynamic increase factor DIF is expressed as(12)DIF=1+AeBlgε−C,where A, B, and C are strain rate parameters, which can be obtained by fitting test data.Substituting the quasi-static cohesionc in (7) with the dynamic cohesion cd in (11) can obtain the strength function considering the strain rate effect, namely,(13)κ=h1−ω2cdcosφ1+sinφ.


### 2.4. Law of Damage Evolution

In order to reflect the degradation of the stiffness and strength of rock materials in the process of deformation and failure, the damage variableω needs to be expressed by appropriate variables. Because the damage of rock materials has different evolution laws under the state of compressive and tensile stress, this paper uses different functions to express the damage variables.In the state of compressive and shear stress, the expression of damage variable in the HJC model [8] was introduced, namely,(14)ωc=∑dεvp+dγpD1P∗+T∗D2,where ωc is the compressive damage variable; P∗ is the standardized hydrostatic pressure, P∗ = P/σc, P is the actual hydrostatic pressure and σc is the uniaxial compressive strength of rock materials; T∗ is the standardized tensile strength, T∗ = σt/σc, σt is the uniaxial tensile strength of rock materials; and D1 and D2 are the damage constants of rock materials.The relationship betweenP and volumetric strain μ can be divided into three stages: linear elastic stage, transition stage, and compaction stage, as shown in Figure 2. In the linear elastic stage (0≤μ<μcrush), P can be expressed as P=Kμ (K=Pcrush/μcrush); in the transition stage, P can be expressed as P= (Plock−Pcrush)/(μplock−μcrush) (μ−μcrush)+Pcrush; in the compaction stage, P can be expressed as P = K1μ¯ + K2μ¯2 + K3μ¯3, where K is the volume modulus; Pcrush is the hydrostatic pressure when the rock reaches its elastic limit; μcrush is the volumetric strain corresponding to P=Pcrush; Plock is the hydrostatic pressure when the rock is compacted; μlock is the volumetric strain corresponding to P = Plock; μlock is the plastic deformation that occurs when the rock is compacted; K1, K2, and K3 are rock materials parameters; and μ¯ is the corrected volumetric strain.Figure 2
Relationship betweenP and μ.When rock materials are in tension, refer to [16], and the tensile damage variable of rock materials can be expressed as(15)ωt=c11−exp−c2ε¯p,where ωt is the tensile damage variable; c1 and c2 are the tensile damage parameters of the rock materials; and ε¯p is the equivalent plastic strain.


### 2.5. Constitutive Relationship

When the external load exceeds the yield strength of the rock materials, the plastic deformation of the rock materials will appear. In order to reflect the plastic behavior of the rock materials, it is necessary to calculate the plastic strain through the plastic potential function. According to the plastic potential theory, there is(16)dεp=dλ∂G∂σ,where dεp is the plastic strain increment; dλ is a plastic multiplier, dλ ≥ 0; and G is the plastic potential function. The plastic potential function of the unified strength theory can be expressed as follows [29]:(17)G=131−a∗I1+13J2a∗−b−2a∗b−11+bsinθL+J2a∗+b+11+bcosθLθL∈−π6,arctan31−a∗31+a∗131−a∗I1+13J2a∗+2b+a∗b−11+bsinθL+J2a∗+a∗b+11+bcosθLθL∈arctan31−a∗31+a∗,π6a∗=1−sinψ1+sinψ,where ψ is the expansion angle of rock materials. The plastic flow of rock materials is generally unassociated, and the flow rule can be controlled by the value of expansion angle.According to plastic mechanics, the loading condition of the yield surface can be expressed as(18)Fσ,ω,γp=0,dλ≥0,Fσ,ω,γpdλ=0.According to the plastic consistency condition,(19)dF=∂F∂σdσ+∂F∂γpdγp+∂F∂ωdω=0.According to equations (5), (9), (16), and (19), the plastic multiplier can be calculated:(20)dλ=∂F/∂σTDdε+∂F/∂ω−∂F/∂σTD0ε−εpdω∂F/∂σTD∂G/∂σ−2/3∂G/∂σij−1/3δijδmn∂G/∂σmn∂G/∂σij−1/3δijδmn∂G/∂σmn∂F/∂γp.According to equations (5), (16), and (20) the stress increment can be obtained:(21)dσ=D−D∂G/∂σ∂F/∂σTD∂F/∂σTD∂G/∂σ−Q∂F/∂γpdε−D∂G/∂σ∂F/∂ω−∂F/∂σTD0ε−εp∂F/∂σTD∂G/∂σ−Q∂F/∂γp+D0ε−εpdω,where Q=2/3∂G/∂σij−1/3δijδmn∂G/∂σmn∂G/∂σij−1/3δijδmn∂G/∂σmn.


## 2.1. Basic Relationship

According to elastic mechanics, the relationship between stress and strain can be expressed as(1)σ=Dεe,where σ is the stress vector; D is the elastic stiffness matrix; and εe is the elastic strain vector.According to elastoplastic mechanics, the relationship between total strain (increment), elastic strain (increment), and plastic strain (increment) is(2)ε=εe+εp,(3)dε=dεe+dεp,where ε is the total strain vector; εp is the plastic strain vector; dε is the total strain increment vector; dεe is an elastic strain increment vector; and dεp is the plastic strain increment vector.According to damage mechanics, the relation between the elastic stiffness matrix of damaged materials and that of undamaged materials is(4)D=1−ωD0,where ω is the damage variable and D0 is the elastic stiffness matrix of the undamaged materials.Combining equations (1)–(4), the relationship between stress increment and strain increment can be obtained as(5)dσ=1−ωD0dε−dεp−D0ε−εpdω.


## 2.2. Strength Criterion

This paper chooses the unified strength theory as the plastic yield criterion of rock materials. The unified strength theory is established on the basis of the twin shear orthogonal octahedral element model, which stipulates that the material will fail when a certain function of the two larger principal shear stresses and their corresponding normal stress reaches a certain threshold [20]. The unified strength theory takes into account all stress components and their different effects on material failure, which can be applied to various geotechnical materials, and has been widely used in geotechnical engineering and other fields [21–23]. The unified strength theory is expressed in terms of stress invariant and Lode angle as(6)F=131−aI1+13J2a−b−2ab−11+bsinθL+J2a+b+11+bcosθL−κθL∈−π6,arctan31−a31+a131−aI1+13J2a+2b+ab−11+bsinθL+J2a+ab+11+bcosθL−κθL∈arctan31−a31+a,π6,where I1 is the invariant of the first principal stress tensor; J2 is the invariant of the second deviator stress tensor; θL is the Lode angle; a is the ratio of tensile and compression strength of the rock materials; b is a parameter reflecting the influence of intermediate principal stress on material failure; in the convexity theory, the value range of b is 0 ≤ b ≤ 1; and κ is the strength function. In order to reflect the strength degradation of the rock materials during the deformation process and the hardening behavior of the rock materials, the damage variable ω and the hardening function are introduced into the strength function, and then κ can be expressed as(7)κ=h1−ω2ccosφ1+sinφ,where h is the hardening function; φ is the internal friction angle of the rock materials; and c is the cohesion of the rock materials.Rock materials generally show hardening mechanical behavior under compressive load but do not have this mechanical behavior under tensile load. In order to distinguish this behavior of rock materials, this paper introduces the following function to express the hardening behavior of rock materials:(8)h=β0+βm−β0γpb1+γp,p≤0,1,p>0,where β0 and βm are the initial critical value and maximum critical value of plastic hardening, respectively; b1 is the parameter that controls the hardening rate; γp is the effective plastic strain; and p is the hydrostatic pressure. In this paper, in order to be consistent with elastic mechanics, the tensile stress is taken as positive. The effective plastic strain increment can be calculated by the following formula:(9)dγp=23dεijp−13δijdεvpdεijp−13δijdεvp,where dγp is the effective plastic strain increment; dεijp is the plastic strain increment; dεvp is the volumetric plastic strain increment; and δij is the Kronecker symbol.The unified strength theory is composed of a series of piecewise linear strength criteria on the deviating plane [20], and its changes are shown in Figure 1. The specific form of the unified strength theory expression depends on the choice of parameter b. When b takes different values, the unified strength theory can be simplified to many current strength criteria. For example, when b = 0 (0 <a < 1), the unified strength theory becomes the Mohr–Coulomb criterion; when b = 1 (0 < a < 1), the unified strength theory becomes the generalized twin shear stress criterion.Figure 1
Locus of unified strength theory in deviatoric plane.


## 2.3. Strain Rate Effect

Rock is a strain rate-dependent geological material, and the strain rate has a significant influence on the mechanical behavior of the rock [24, 25]. The most direct effect of strain rate on rock mechanical behavior is to increase the strength of rock. In rock dynamics, dynamic increase factor is generally used to characterize the influence of strain rate on rock strength [26], and dynamic increase factor can be defined as(10)DIF=σdσs,where DIF is the dynamic increase factor; σs is the quasi-static strength of rock materials; and σd is the dynamic strength of the rock materials.According to references [26–28], it can be seen that the strain rate has little effect on the friction strength of the rock and can be ignored; the strain rate has a greater influence on the cohesive strength of the rock, so the dynamic cohesive strength of rock can be expressed as a function of strain rate and quasi-static cohesive strength, namely,(11)cd=DIFc,where cd is the dynamic cohesion of rock materials. In this paper, the dynamic increase factor DIF is expressed as(12)DIF=1+AeBlgε−C,where A, B, and C are strain rate parameters, which can be obtained by fitting test data.Substituting the quasi-static cohesionc in (7) with the dynamic cohesion cd in (11) can obtain the strength function considering the strain rate effect, namely,(13)κ=h1−ω2cdcosφ1+sinφ.


## 2.4. Law of Damage Evolution

In order to reflect the degradation of the stiffness and strength of rock materials in the process of deformation and failure, the damage variableω needs to be expressed by appropriate variables. Because the damage of rock materials has different evolution laws under the state of compressive and tensile stress, this paper uses different functions to express the damage variables.In the state of compressive and shear stress, the expression of damage variable in the HJC model [8] was introduced, namely,(14)ωc=∑dεvp+dγpD1P∗+T∗D2,where ωc is the compressive damage variable; P∗ is the standardized hydrostatic pressure, P∗ = P/σc, P is the actual hydrostatic pressure and σc is the uniaxial compressive strength of rock materials; T∗ is the standardized tensile strength, T∗ = σt/σc, σt is the uniaxial tensile strength of rock materials; and D1 and D2 are the damage constants of rock materials.The relationship betweenP and volumetric strain μ can be divided into three stages: linear elastic stage, transition stage, and compaction stage, as shown in Figure 2. In the linear elastic stage (0≤μ<μcrush), P can be expressed as P=Kμ (K=Pcrush/μcrush); in the transition stage, P can be expressed as P= (Plock−Pcrush)/(μplock−μcrush) (μ−μcrush)+Pcrush; in the compaction stage, P can be expressed as P = K1μ¯ + K2μ¯2 + K3μ¯3, where K is the volume modulus; Pcrush is the hydrostatic pressure when the rock reaches its elastic limit; μcrush is the volumetric strain corresponding to P=Pcrush; Plock is the hydrostatic pressure when the rock is compacted; μlock is the volumetric strain corresponding to P = Plock; μlock is the plastic deformation that occurs when the rock is compacted; K1, K2, and K3 are rock materials parameters; and μ¯ is the corrected volumetric strain.Figure 2
Relationship betweenP and μ.When rock materials are in tension, refer to [16], and the tensile damage variable of rock materials can be expressed as(15)ωt=c11−exp−c2ε¯p,where ωt is the tensile damage variable; c1 and c2 are the tensile damage parameters of the rock materials; and ε¯p is the equivalent plastic strain.


## 2.5. Constitutive Relationship

When the external load exceeds the yield strength of the rock materials, the plastic deformation of the rock materials will appear. In order to reflect the plastic behavior of the rock materials, it is necessary to calculate the plastic strain through the plastic potential function. According to the plastic potential theory, there is(16)dεp=dλ∂G∂σ,where dεp is the plastic strain increment; dλ is a plastic multiplier, dλ ≥ 0; and G is the plastic potential function. The plastic potential function of the unified strength theory can be expressed as follows [29]:(17)G=131−a∗I1+13J2a∗−b−2a∗b−11+bsinθL+J2a∗+b+11+bcosθLθL∈−π6,arctan31−a∗31+a∗131−a∗I1+13J2a∗+2b+a∗b−11+bsinθL+J2a∗+a∗b+11+bcosθLθL∈arctan31−a∗31+a∗,π6a∗=1−sinψ1+sinψ,where ψ is the expansion angle of rock materials. The plastic flow of rock materials is generally unassociated, and the flow rule can be controlled by the value of expansion angle.According to plastic mechanics, the loading condition of the yield surface can be expressed as(18)Fσ,ω,γp=0,dλ≥0,Fσ,ω,γpdλ=0.According to the plastic consistency condition,(19)dF=∂F∂σdσ+∂F∂γpdγp+∂F∂ωdω=0.According to equations (5), (9), (16), and (19), the plastic multiplier can be calculated:(20)dλ=∂F/∂σTDdε+∂F/∂ω−∂F/∂σTD0ε−εpdω∂F/∂σTD∂G/∂σ−2/3∂G/∂σij−1/3δijδmn∂G/∂σmn∂G/∂σij−1/3δijδmn∂G/∂σmn∂F/∂γp.According to equations (5), (16), and (20) the stress increment can be obtained:(21)dσ=D−D∂G/∂σ∂F/∂σTD∂F/∂σTD∂G/∂σ−Q∂F/∂γpdε−D∂G/∂σ∂F/∂ω−∂F/∂σTD0ε−εp∂F/∂σTD∂G/∂σ−Q∂F/∂γp+D0ε−εpdω,where Q=2/3∂G/∂σij−1/3δijδmn∂G/∂σmn∂G/∂σij−1/3δijδmn∂G/∂σmn.


## 3. Numerical Implementation of Model

In order to verify the correctness of the constitutive model and facilitate the application of the model, the constitutive model needs to be numerically implemented to achieve the purpose of solution. Since the constitutive model established in this paper is a model of plasticity and damage coupling, it is more complicated to solve the model. This paper uses a semi-implicit stress return algorithm to numerically implement the model [30], a schematic diagram of the semi-implicit stress return algorithm as shown in Figure 3. The algorithm mainly includes the following steps as follows:(1)
Calculate trial stress:σn+1tri=σn+DnΔεn+1;(2)
Determine whetherFσn+1tri,γnp,ωn≤0; if yes, then σn+1=σn+1tri; if no, continue to step (3);(3)
Stress return: whenFσn+1tri,γnp,ωn>0, the trial stress is located outside the yield surface and needs to be pulled back to the yield surface. According to plastic mechanics, on the yield surface, there is(22)Fn+1σn+1,ωn+1,γn+1p=0.For equation (22), perform a Taylor expansion at the trial stress and there is(23)Fn+1=Fn+1tri+∂F∂σTΔσcre|tri+∂F∂γpΔγp|tri+∂F∂ωΔω|tri=0.TheΔσcre in the above equation can be calculated by (5)(24)Δσcre=σn+1−σn+1tri=−1−ωnD0Δεp−σ1−ωnΔω|tri.According to equations (16), (20) (23) and (24),(25)dλ=Ftri∂F/∂σT1−ωnD0∂G/∂σ|tri+∂F/∂σTσ/1−ωn−∂F/∂ωQ1|tri−∂F/∂γpQ|tri,whereQ1=Q+δij∂G/∂σij/D1P∗+T∗D2,p≤0c1c2exp−c2dε¯p∂G/∂σij∂G/∂σij,p>0.(4)
Update stressσn+1, plastic strain εn+1p, effective plastic strain γn+1p, and damage variables ωn+1.(26)σn+1=σn+Δσn+1tri+Δσcre︸Δσn+1,εn+1p=εnp−Cn+1Δσcre,γn+1p=γnp+Δγn+1p,ωn+1=ωn+Δωn+1.According to the above calculation process, the Fortran language is used to program the solution of the constitutive model, and the programmed constitutive model is imported into the LS-DYNA software through the user material custom interface UMAT.Figure 3
Schematic diagram of the semi-implicit stress return algorithm.


## 4. Validation of Model

The correctness of a model needs to be verified by tests. In this part, the rock dynamic model established in this paper is verified by uniaxial and triaxial compression tests, uniaxial tensile tests, and ballistic tests.


### 4.1. Uniaxial Compression and Triaxial Compression Tests of Rock

This experiment uses marble as the research object, and the basic physical and mechanical parameters of marble are as follows: densityρ = 2680 kg/m3, elastic modulus E = 71.91 GPa, Poisson’s ratio υ = 0.20, cohesion c = 34.65 MPa, and internal friction angle φ = 34.25°. The marble finite element model in the numerical simulation has the same shape and size as the marble specimen in the experiment; that is, the shape is cylindrical, and the size is Φ50 mm × 100 mm [31]. The finite element model and boundary conditions are shown in Figure 4. Other model parameters used in numerical simulation are as follows: A = 0, B = 0, C = 0, β0 = 0.8, βm = 1.15, b1 = 1.8 × 10−4, D1 = 0.015, D2 = 1.0, μcrush = 0.8761 × 10−3, μlock = 0.011, Pcrush = 35 MPa, Plock = 0.8 GPa, K1 = 51.86 GPa, K2 = −54.12 GPa, K3 = 1359.66 GPa, c1 = 0.98, and c2 = 300. The detailed determination of model parameters can refer to [5].Figure 4
Finite element models and boundary conditions of marble (a) uniaxial compression and (b) triaxial compression.
(a)(b)Figure5 shows the stress-strain curves of marble under different confining pressures. It can be seen from the figure that the numerical simulation results of the stress-strain curve of the marble are in good agreement with the experimental data. The stress-strain curves of marble under different confining pressures have experienced three stages: linear elasticity, hardening, and softening; as the confining pressure increases, the marble gradually changes from brittleness to ductility.Figure 5
Stress-strain curves of marble under different confining pressures: (a)σ3 = 0 MPa; (b) σ3 = 3.45 MPa; (c) σ3 = 6.90 MPa; (d) σ3 = 13.80 MPa.
(a)(b)(c)(d)


### 4.2. Uniaxial Tensile Test of Rock

In this test, tuff is taken as the research object, and the basic physical and mechanical parameters of tuff are as follows: densityρ = 1920 kg/m3, elastic modulus E = 5.87 GPa, Poisson’s ratio υ = 0.23, cohesion c = 4.72 MPa, and internal friction angle φ = 58°. Other model parameters used in the numerical simulation of tuff uniaxial tensile tests are as follows: A = 0, B = 0，C = 0, β0 = 1.0, βm = 1.0，b1 = 1.8 × 10−4, D1 = 0.04, D2 = 1.0, μcrush = 3.04 × 10−3, μlock = 0.0245, Pcrush = 11 MPa, Plock = 0.54 GPa, K1 = 9.76 GPa, K2 = -21.80 GPa, K3 = 401.85 GPa, c1 = 0.98, and c2 = 160. The finite element model is consistent with the test sample; that is, the shape is cylindrical and the size is Φ25 mm × 50 mm [32]. The tuff finite element model and boundary conditions are shown in Figure 6.Figure 6
Uniaxial tension finite element model and boundary conditions of tuff.Figure7 shows the uniaxial tensile stress-strain curve of tuff. From the figure, the tuff has undergone two processes of linear elasticity and softening under uniaxial tension; the numerical simulation results of the uniaxial tensile stress-strain curve of tuff have a good consistency with the experimental data.Figure 7
Uniaxial tensile stress-strain curve of tuff.


### 4.3. Ballistic Test of Rock

In this part, the numerical simulation of ballistic penetration and perforation tests performed by Seah et al. [33] was carried out. In the test, a steel bullet was fired from a compressed gas gun, and its initial velocity was measured by a photovoltaic system before it finally hit the target. The steel bullet is 94.9 mm long, 20 mm in diameter, and 0.2 kg in mass. The detailed dimensions are shown in Figure 8. The size of the granite target is 600 mm × 600 mm × 100 mm, and the weight is about 100 kg. It is anchored by rigid steel plates. The granite ballistic test device is shown in Figure 9. In order to reduce calculation cost and improve calculation efficiency, only a quarter of the granite target and steel bullets was modeled, as shown in Figure 10. Steel bullets and granite targets use three-dimensional solid elements SOLID 164. In order to save calculation time and ensure calculation accuracy, local mesh refinement of granite target was carried out, namely, the element size of the granite target is 1 mm near the impact site and expands to 3 mm after exceeding 3 times the diameter of the bullet. The element size of the bullet is 1 mm. The erosion surface contact algorithm is used to reflect the contact behavior between bullet and granite target. The basic physical and mechanical parameters of granite in this test are as follows: density ρ = 2626 kg/m3, elastic modulus E = 54 GPa, Poisson’s ratio υ = 0.27, uniaxial compressive strength σc = 163 MPa, and uniaxial tensile strength σt = 7.1 MPa. Other model parameters used in the simulation are as follows: A = 0.84, B = 0.80, C = 1.50, β0 = 1.0, βm = 1.0, b1 = 1.8×10–4, D1 = 0.04，D2 = 1.0, μcrush = 1.39×10–3, μlock = 0.1, Pcrush = 54.33 MPa, Plock = 0.8 GPa, K1 = 85 GPa, K2 = –171 GPa，K3 = 208 GPa, c1 = 0.98, and c2 = 100. The model parameters A, B, and C can be determined by fitting experimental data [1], as shown in Figure 11. The bullet adopts a rigid body material model [17], and its material parameters are as follows: density ρ = 7850 kg/m3, elastic modulus E = 210 GPa, Poisson’s ratio υ = 0.2, and yield strength σy = 1900 MPa.Figure 8
Schematic diagram of steel bullet and its dimensions.Figure 9
Device of granite ballistic experiment [33].Figure 10
Finite element model (a quarter of the original test model) and boundary conditions of the granite ballistic test.Figure 11
Relationship between dynamic increase factor and strain rate of granite.Figure12 shows the test results (see Figure 12(a)) and numerical results (see Figure 12(b)) of the failure mode of the granite target. It can be seen from Figure 12(a) that the granite target plate first formed a crater under the impact of the bullet and finally formed a cone plug with the penetration of the bullet. The numerical simulation results of the granite failure form are highly consistent with the experimental results. When the bullet impacts the granite target plate at an initial velocity of 288 m/s, the height of the conical plug and the diameter of the crater on the back of the target plate are shown in Table 1. It can be seen from Table 1 that the numerical simulation results of conical plug height are highly consistent with the test results, with a relative error of only 1.4%, which is far less than the relative error between the analytical value and the test results (18.5%). Although the relative error between the numerical simulation results and the test results of the crater diameter on the back of the target plate is 21.1%, this result is far less than the relative error between the analytical value and the test results (45.8%). This error may be caused by the microstructure of granite in the test. Therefore, this model can well describe the mechanical behavior of granite under impact load.Figure 12
Comparison of failure form of granite target plate between (a) experiment [33] (typical failure modes of the granite ballistic test) and (b) numerical simulation (the initial velocity of a steel bullet v = 288 m/s).Table 1
The conical plug height of granite target plate impacted by bullets and the diameter of crater on the back of target plate.
Conical plug height (mm) (relative error)Diameter of the crater on the back of the target plate (mm) (relative error)Experiment [30]65312Analytical value [30]77 (18.5%)455 (45.8%)Numerical simulation65.9 (1.4%)246.1 (21.1%)


## 4.1. Uniaxial Compression and Triaxial Compression Tests of Rock

This experiment uses marble as the research object, and the basic physical and mechanical parameters of marble are as follows: densityρ = 2680 kg/m3, elastic modulus E = 71.91 GPa, Poisson’s ratio υ = 0.20, cohesion c = 34.65 MPa, and internal friction angle φ = 34.25°. The marble finite element model in the numerical simulation has the same shape and size as the marble specimen in the experiment; that is, the shape is cylindrical, and the size is Φ50 mm × 100 mm [31]. The finite element model and boundary conditions are shown in Figure 4. Other model parameters used in numerical simulation are as follows: A = 0, B = 0, C = 0, β0 = 0.8, βm = 1.15, b1 = 1.8 × 10−4, D1 = 0.015, D2 = 1.0, μcrush = 0.8761 × 10−3, μlock = 0.011, Pcrush = 35 MPa, Plock = 0.8 GPa, K1 = 51.86 GPa, K2 = −54.12 GPa, K3 = 1359.66 GPa, c1 = 0.98, and c2 = 300. The detailed determination of model parameters can refer to [5].Figure 4
Finite element models and boundary conditions of marble (a) uniaxial compression and (b) triaxial compression.
(a)(b)Figure5 shows the stress-strain curves of marble under different confining pressures. It can be seen from the figure that the numerical simulation results of the stress-strain curve of the marble are in good agreement with the experimental data. The stress-strain curves of marble under different confining pressures have experienced three stages: linear elasticity, hardening, and softening; as the confining pressure increases, the marble gradually changes from brittleness to ductility.Figure 5
Stress-strain curves of marble under different confining pressures: (a)σ3 = 0 MPa; (b) σ3 = 3.45 MPa; (c) σ3 = 6.90 MPa; (d) σ3 = 13.80 MPa.
(a)(b)(c)(d)


## 4.2. Uniaxial Tensile Test of Rock

In this test, tuff is taken as the research object, and the basic physical and mechanical parameters of tuff are as follows: densityρ = 1920 kg/m3, elastic modulus E = 5.87 GPa, Poisson’s ratio υ = 0.23, cohesion c = 4.72 MPa, and internal friction angle φ = 58°. Other model parameters used in the numerical simulation of tuff uniaxial tensile tests are as follows: A = 0, B = 0，C = 0, β0 = 1.0, βm = 1.0，b1 = 1.8 × 10−4, D1 = 0.04, D2 = 1.0, μcrush = 3.04 × 10−3, μlock = 0.0245, Pcrush = 11 MPa, Plock = 0.54 GPa, K1 = 9.76 GPa, K2 = -21.80 GPa, K3 = 401.85 GPa, c1 = 0.98, and c2 = 160. The finite element model is consistent with the test sample; that is, the shape is cylindrical and the size is Φ25 mm × 50 mm [32]. The tuff finite element model and boundary conditions are shown in Figure 6.Figure 6
Uniaxial tension finite element model and boundary conditions of tuff.Figure7 shows the uniaxial tensile stress-strain curve of tuff. From the figure, the tuff has undergone two processes of linear elasticity and softening under uniaxial tension; the numerical simulation results of the uniaxial tensile stress-strain curve of tuff have a good consistency with the experimental data.Figure 7
Uniaxial tensile stress-strain curve of tuff.


## 4.3. Ballistic Test of Rock

In this part, the numerical simulation of ballistic penetration and perforation tests performed by Seah et al. [33] was carried out. In the test, a steel bullet was fired from a compressed gas gun, and its initial velocity was measured by a photovoltaic system before it finally hit the target. The steel bullet is 94.9 mm long, 20 mm in diameter, and 0.2 kg in mass. The detailed dimensions are shown in Figure 8. The size of the granite target is 600 mm × 600 mm × 100 mm, and the weight is about 100 kg. It is anchored by rigid steel plates. The granite ballistic test device is shown in Figure 9. In order to reduce calculation cost and improve calculation efficiency, only a quarter of the granite target and steel bullets was modeled, as shown in Figure 10. Steel bullets and granite targets use three-dimensional solid elements SOLID 164. In order to save calculation time and ensure calculation accuracy, local mesh refinement of granite target was carried out, namely, the element size of the granite target is 1 mm near the impact site and expands to 3 mm after exceeding 3 times the diameter of the bullet. The element size of the bullet is 1 mm. The erosion surface contact algorithm is used to reflect the contact behavior between bullet and granite target. The basic physical and mechanical parameters of granite in this test are as follows: density ρ = 2626 kg/m3, elastic modulus E = 54 GPa, Poisson’s ratio υ = 0.27, uniaxial compressive strength σc = 163 MPa, and uniaxial tensile strength σt = 7.1 MPa. Other model parameters used in the simulation are as follows: A = 0.84, B = 0.80, C = 1.50, β0 = 1.0, βm = 1.0, b1 = 1.8×10–4, D1 = 0.04，D2 = 1.0, μcrush = 1.39×10–3, μlock = 0.1, Pcrush = 54.33 MPa, Plock = 0.8 GPa, K1 = 85 GPa, K2 = –171 GPa，K3 = 208 GPa, c1 = 0.98, and c2 = 100. The model parameters A, B, and C can be determined by fitting experimental data [1], as shown in Figure 11. The bullet adopts a rigid body material model [17], and its material parameters are as follows: density ρ = 7850 kg/m3, elastic modulus E = 210 GPa, Poisson’s ratio υ = 0.2, and yield strength σy = 1900 MPa.Figure 8
Schematic diagram of steel bullet and its dimensions.Figure 9
Device of granite ballistic experiment [33].Figure 10
Finite element model (a quarter of the original test model) and boundary conditions of the granite ballistic test.Figure 11
Relationship between dynamic increase factor and strain rate of granite.Figure12 shows the test results (see Figure 12(a)) and numerical results (see Figure 12(b)) of the failure mode of the granite target. It can be seen from Figure 12(a) that the granite target plate first formed a crater under the impact of the bullet and finally formed a cone plug with the penetration of the bullet. The numerical simulation results of the granite failure form are highly consistent with the experimental results. When the bullet impacts the granite target plate at an initial velocity of 288 m/s, the height of the conical plug and the diameter of the crater on the back of the target plate are shown in Table 1. It can be seen from Table 1 that the numerical simulation results of conical plug height are highly consistent with the test results, with a relative error of only 1.4%, which is far less than the relative error between the analytical value and the test results (18.5%). Although the relative error between the numerical simulation results and the test results of the crater diameter on the back of the target plate is 21.1%, this result is far less than the relative error between the analytical value and the test results (45.8%). This error may be caused by the microstructure of granite in the test. Therefore, this model can well describe the mechanical behavior of granite under impact load.Figure 12
Comparison of failure form of granite target plate between (a) experiment [33] (typical failure modes of the granite ballistic test) and (b) numerical simulation (the initial velocity of a steel bullet v = 288 m/s).Table 1
The conical plug height of granite target plate impacted by bullets and the diameter of crater on the back of target plate.
Conical plug height (mm) (relative error)Diameter of the crater on the back of the target plate (mm) (relative error)Experiment [30]65312Analytical value [30]77 (18.5%)455 (45.8%)Numerical simulation65.9 (1.4%)246.1 (21.1%)


## 5. Conclusion

Based on elastoplastic mechanics and damage mechanics, with unified strength theory as the yield criterion, a dynamic model of rock coupled with elastoplastic damage is established. Some conclusions can be made as follows:(1)
This model fully considers the strain rate effect of rock materials, the Lode angle effect, and the nonlinearity of the hydrostatic pressure of rock materials and reflects the different damage evolution laws of rock materials and different hardening behaviors under compression and tension conditions.(2)
Using Fortran language and using the semi-implicit stress return algorithm, the dynamic constitutive model of the rock established is numerically realized, thereby facilitating the application of the model.(3)
Through three classic examples of rock uniaxial and triaxial compression test, uniaxial tensile test, and rock ballistic test, the established constitutive model was verified. The results show that the model can well describe the mechanical behavior of rock materials under static and dynamic loads.



---

*Source: 5567019-2021-10-04.xml*