A Coupled Elastoplastic Damage Dynamic Model for Rock

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, uniﬁed strength theory is taken as the yield criterion; to reﬂect the diﬀerent damage evolution law of rocks under tension and pressure conditions, the eﬀective 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 eﬀect of rocks are characterized by piecewise function and dynamic increase factor function, respectively; Fortran language and LS-DYNA User-Deﬁned Interface (Umat) are used to numerically implement the constitutive model; the constitutive model is veriﬁed 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.


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][2][3][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][6][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]. e KCC model cannot accurately capture the tensile behavior of rock materials [13]. e 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]. e 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]. erefore, the improvement of the classical model and the establishment of the new model are still one of our joint efforts. e 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.

Basic Relationship.
According to elastic mechanics, the relationship between stress and strain can be expressed as 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 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 where ω is the damage variable and D 0 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

Strength Criterion.
is paper chooses the unified strength theory as the plastic yield criterion of rock materials. e 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]. e 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][22][23]. e unified strength theory is expressed in terms of stress invariant and Lode angle as where I 1 is the invariant of the first principal stress tensor; J 2 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 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: 2 Shock and Vibration where β 0 and β m are the initial critical value and maximum critical value of plastic hardening, respectively; b 1 is the parameter that controls the hardening rate; c 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. e effective plastic strain increment can be calculated by the following formula: where dc p is the effective plastic strain increment; dε p ij is the plastic strain increment; dε p v is the volumetric plastic strain increment; and δ ij is the Kronecker symbol. e 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. e 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.

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]. e 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 where DIF is the dynamic increase factor; σ s is the quasistatic strength of rock materials; and σ d is the dynamic strength of the rock materials. According to references [26][27][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, where c d is the dynamic cohesion of rock materials. In this paper, the dynamic increase factor DIF is expressed as where A, B, and C are strain rate parameters, which can be obtained by fitting test data. Substituting the quasi-static cohesion c in (7) with the dynamic cohesion c d in (11) can obtain the strength function considering the strain rate effect, namely,

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, 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 D 1 and D 2 are the damage constants of rock materials. e relationship between P 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 � P crush /μ crush ); in the transition stage, P can be expressed as P � (P lock − P crush )/(μ plock − μ crush ) (μ − μ crush ) + P crush ; in the compaction stage, P can be expressed as P � K 1 μ + K 2 μ 2 + K 3 μ 3 , where K is the volume modulus; P crush is the hydrostatic pressure when the rock reaches its elastic limit; μ crush is the volumetric strain Figure 1: Locus of unified strength theory in deviatoric plane. Shock and Vibration 3 corresponding to P � P crush ; P lock is the hydrostatic pressure when the rock is compacted; μ lock is the volumetric strain corresponding to P � P lock ; μ lock is the plastic deformation that occurs when the rock is compacted; K 1 , K 2 , and K 3 are rock materials parameters; and μ is the corrected volumetric strain.
When rock materials are in tension, refer to [16], and the tensile damage variable of rock materials can be expressed as where ω t is the tensile damage variable; c 1 and c 2 are the tensile damage parameters of the rock materials; and ε p is the equivalent plastic strain.

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 where dε p is the plastic strain increment; dλ is a plastic multiplier, dλ ≥ 0; and G is the plastic potential function. e plastic potential function of the unified strength theory can be expressed as follows [29]: where ψ is the expansion angle of rock materials. e 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 According to the plastic consistency condition, According to equations (5)

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. is paper uses a semi-implicit stress return algorithm to numerically implement the model [30], a schematic diagram of the semiimplicit stress return algorithm as shown in Figure 3. e algorithm mainly includes the following steps as follows: (1) Calculate trial stress: σ tri n+1 � σ n + D n Δε n+1 ; (2) Determine whether F(σ tri n+1 , c p n , ω n ) ≤ 0; if yes, then σ n+1 � σ tri n+1 ; if no, continue to step (3); (3) Stress return: when F(σ tri n+1 , c p n , ω 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 For equation (22), perform a Taylor expansion at the trial stress and there is e Δσ cre in the above equation can be calculated by (5) According to equations (16), (20) (23) and (24), where 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.

Validation of Model
e 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.

Uniaxial Compression and Triaxial Compression Tests of Rock.
is experiment uses marble as the research object, and the basic physical and mechanical parameters of marble are as follows: density ρ � 2680 kg/m 3 , elastic modulus E � 71.91 GPa, Poisson's ratio υ � 0.20, cohesion c � 34.65 MPa, and internal friction angle φ � 34.25°. e marble finite element model in the numerical simulation Shock and Vibration 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]. e finite element model and boundary conditions are shown in Figure 4. Other model parameters used in numerical simulation are as follows: e detailed determination of model parameters can refer to [5]. Figure 5 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. e 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.

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/m 3 , 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: μ lock � 0.0245, P crush � 11 MPa, P lock � 0.54 GPa, K 1 � 9.76 GPa, K 2 � -21.80 GPa, K 3 � 401.85 GPa, c 1 � 0.98, and c 2 � 160. e finite element model is consistent with the test sample; that is, the shape is cylindrical and the size is Φ25 mm × 50 mm [32]. e tuff finite element model and boundary conditions are shown in Figure 6. Figure 7 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 stressstrain curve of tuff have a good consistency with the experimental data.

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. e steel bullet is 94.9 mm long, 20 mm in diameter, and 0.2 kg in mass. e detailed dimensions are shown in Figure 8. e 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. e 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. e element size of the bullet is 1 mm. e erosion surface contact algorithm is used to reflect the contact behavior between bullet and granite target. e basic physical and mechanical parameters of granite in this test are as follows: density ρ � 2626 kg/m 3 , 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: 33 MPa, P lock � 0.8 GPa, K1 � 85 GPa, K2 � -171 GPa，K3 � 208 GPa, c 1 � 0.98, and c 2 � 100. e model parameters A, B, and C can be determined by fitting experimental data [1], as shown in Figure 11. e bullet adopts a rigid body material model [17], and its material parameters are as follows: density ρ � 7850 kg/ m 3 , elastic modulus E � 210 GPa, Poisson's ratio υ � 0.2, and yield strength σ y � 1900 MPa. Figure 12 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. e 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%). is error may be caused by the microstructure of granite in the test. erefore, this model can well describe the mechanical behavior of granite under impact load.      Figure 11: Relationship between dynamic increase factor and strain rate of granite.

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) is 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.

Data Availability
e data used to support the findings of the study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.