# Determination of Constitutive Parameters of Crystalline Limestone Based on Improved RHT Model

Authors: Hao-Fan Tian; Tai Bao; Zhuo Li; Hao Peng; Shuai You; Shi-li Xiao
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

This study improves the RHT constitutive model based on some shortcomings, and the main improvements work are as follows: (1) hyperbolic function is used to make up for the deficiency of tensile strain rate enhancement factor of RHT constitutive model at high strain rate; (2) considering the constant change of hydrostatic pressure, the initial damage value of rock is introduced into the constitutive model, and the improved damage accumulation expression is derived based on this value. (3) In the equation of residual failure surface, the influence of the Lode angle factor is considered, the RHT constitutive parameters of crystalline limestone are determined by the SHPB experiment, and the difference in rock failure forms before and after the correction of constitutive parameters is simulated by ANSYS/LS-DYNA software. The results show that the improved constitutive parameters can effectively represent the yield strength surface and residual strength surface of the rock, which verifies the rationality of the improved RHT constitutive model, and also shows that it is an effective and simple method to obtain RHT constitutive parameters by SHPB experiment.


---

## Body

## 1. Introduction

Compared with the HJC [1] and K & C [2] model, the RHT model [3] fully considers the compression effect, strain rate effect, and damage accumulation on the failure intensity of the rock under explosive impact. The elastic limit surface equation, failure surface equation, and residual strength surface equation related to pressure are embedded in this model, which is used to describe the variation law of yield strength, failure strength, and residual strength of concrete and brittle materials under high strain rate. With the rapid development of numerical simulation technology, the model is widely used in the numerical simulation of high strain rates such as explosion impact. When Tu, Hansson, and Skoglund [4–8] et al. applied the RHT model, they found that the model was incomplete for tensile damage description, so the model was modified for tensile and high strain rate effects. Zhang ruoqi and Pavlovic A [9–14] et al. adjusted and improved the failure surface equation and residual surface equation in RHT constitutive parameters based on concrete experiments and verified them by AUTODYN element simulation, all of which achieved good results. A total of 34 constitutive parameters of the RHT model need to be determined. Most scholars can only refer to Riedel’s [3, 15–17] research results because of its complexity when applying the RHT model. Li [18] used marble as an example to explain in detail the method for determining the constitutive parameters of RHT. Xie [19, 20] and Wang [21–23] et al. studied the constitutive parameters of rock under the action of high ground stress and successfully applied them to the numerical simulation analysis of rock crack propagation, cut hole blasting, and cyclic blasting. For the convenience of parameter determination, Xu and Ye[24], Liu et al. [25], and Li et al. [26–28] used an orthogonal experimental method and optimized Latin hypercube design method to analyze the sensitivity of the parameters. The analysis results show that the strength parameters of failure surface and residual surface have a sensitive influence on the fracture morphology and damage accumulation of rock. Therefore, some parameters that have little influence on rock failure can be neglected when determining rock constitutive parameters, which makes the RHT model more convenient to use.The study of the above scholars for the RHT constitutive model mainly includes the following: the modification of the constitutive equation, the sensitivity analysis of the model parameters, and the method of determining the constitutive parameters. In the improvement of the constitutive equation, most scholars use the linear interpolation method for the residual strength of the rock in the damage softening stage, while the initial damage of rock has not been considered. In fact, for rocks in the karst area, due to the development of karst phenomena, the rock has certain initial damage, which will reduce the strength of the rock. In addition, when rock is subjected to dynamic loads such as explosion, vibration, and impact, the hydrostatic pressure on its failure surface often changes with the change of load, so it is necessary to update the expression of damage accumulation, considering that many complex experiments are usually needed in the process of determining RHT constitutive parameters. Based on this, the residual surface equation of the RHT constitutive model is improved, and the constitutive parameters of crystalline limestone are determined by the Hopkinson bar impact test.


## 2. RHT Constitutive Model

### 2.1. Failure Surface Equation

In the RHT model, the equivalent stress intensityσeq∗ of the failure surface is a function of normalized pressure p∗=p/fc, Lode angle θ, and strain rate ε˙.(1)σeq∗p,θ,ε˙=YTXC∗pR3θFrateε˙.In formula (1), the ∗ in the upper right corner of the formula indicates the normalized uniaxial compressive strength fc (unit: MPa) of the material, and YTXC∗p is the compression meridian strength.(2)YTXCp∗=Ap∗−pspall∗Frateε˙N.In formula (2), A is the failure surface constant; N is the failure surface index; and pspall∗ is the normalized layer crack strength.Frateε˙ is the strain rate dynamic enhancement factor:(3)Frateε˙=ε˙ε˙0cα,p≥fc3,p+ft/3fc/3+ft/3ε˙ε˙0cα+p−fc/3−ft/3−fc/3ε˙ε˙0cδ,−ft3<p<fc3ε˙ε˙0tδ.p≤−ft3.,In formula (3), ft is uniaxial tensile strength (unit: MPa)； ε˙0c=30×10−6s−1； ε˙0t=3×10−6s−1； ε˙ is the current strain rate; α=4/20+3fc is a compression strain rate index; and δ=2/20+fc is the tensile strain rate index.R3θ is the Lode angle factor: it is a function of the Lode angle θ and the tension-compression meridian ratio Q.(4)R3θ=21−Q2cosθ+2Q−141−Q2cos2θ+5Q2−4Q41−Q2cos2θ+2Q−12.In formula (4), θ=cos−133/2J3/J23/2/3θ∈0,π/3； J2 and J3 are deviatoric stress invariants; and Q is the meridian ratio of tension and compression, and in order to consider the effect of pressure on strength, it is expressed as follows:(5)Q=Qp∗=Q0+Bp∗0.51≤Q≤1.In formula (5), Q0 is the initial meridian ratio of tension and compression, and B is the material constant.


### 2.2. Elastic Ultimate Surface Equation

In the RHT model, the elastic limit surface equation is derived from the failure surface equation:(6)σel∗p,θ,ε˙=σeq∗FelasticFcapp∗,In formula (6), Felastic is the elastic scaling function and Fcapp∗ is the cap function, which can effectively solve the problem that materials are always in the elastic state without yielding under high hydrostatic pressure. The expressions for Felastic, Fcapp∗ are as follows:(7)Felastic=gc∗,p≥fc,el3,p+ft,el/3fc,el/3+ft,el/3gc∗+p−fc,el/3−ft,el/3−fc,el/3gt∗,−ft,el3<p<fc,el3,gt∗,p≤−ft,el3.In equation (7), ft,el is the uniaxial tensile elastic ultimate stress, fc,el is the uniaxial compressive elastic ultimate stress, gc∗=fc,el/fc is the compressive yield surface parameter, and gt∗=ft,el/ft is the tensile yield surface parameter.(8)Fcapp∗=1,p∗⩽13,1−p∗−1/3p0∗−1/32,13<p∗<p0∗,0,p⩾p0∗.In equation (8), p0 is the pressure at which the void material begins to crush.


### 2.3. Residual Strength Surface Equation

When the equivalent stress strength of the material is greater than the failure stress strength, the damage of the material begins to accumulate and enters the damage softening phase. Damage variableD is the ratio of cumulative equivalent plastic strain increment to final failure equivalent plastic strain. The plastic strain of the material at ultimate failure is as follows:(9)εpf=D1p∗−1−Dpt∗D2p∗⩾1−Dpt∗+εpmD11/D2εpmp∗<1−Dpt∗+εpmD11/D2,In formula (9), εpm is the minimum plastic strain at the time of material damage, and pt∗ is the destruction of the cutoff pressure. D1 and D2 are material parameters. In fact, when the pressure of the material exceeds the yield stress, the friction and confining pressure still exist in the broken part of the material, which makes it have the ability to resist shear. Therefore, the residual strength surface equation is introduced into the RHT model.(10)σr∗p∗=Afp∗nf,In formula (10), Af is the residual stress intensity parameter and nf is the residual stress intensity index. The equivalent stresses between the failure strength surface and the residual strength surface are obtained by linear interpolation.


## 2.1. Failure Surface Equation

In the RHT model, the equivalent stress intensityσeq∗ of the failure surface is a function of normalized pressure p∗=p/fc, Lode angle θ, and strain rate ε˙.(1)σeq∗p,θ,ε˙=YTXC∗pR3θFrateε˙.In formula (1), the ∗ in the upper right corner of the formula indicates the normalized uniaxial compressive strength fc (unit: MPa) of the material, and YTXC∗p is the compression meridian strength.(2)YTXCp∗=Ap∗−pspall∗Frateε˙N.In formula (2), A is the failure surface constant; N is the failure surface index; and pspall∗ is the normalized layer crack strength.Frateε˙ is the strain rate dynamic enhancement factor:(3)Frateε˙=ε˙ε˙0cα,p≥fc3,p+ft/3fc/3+ft/3ε˙ε˙0cα+p−fc/3−ft/3−fc/3ε˙ε˙0cδ,−ft3<p<fc3ε˙ε˙0tδ.p≤−ft3.,In formula (3), ft is uniaxial tensile strength (unit: MPa)； ε˙0c=30×10−6s−1； ε˙0t=3×10−6s−1； ε˙ is the current strain rate; α=4/20+3fc is a compression strain rate index; and δ=2/20+fc is the tensile strain rate index.R3θ is the Lode angle factor: it is a function of the Lode angle θ and the tension-compression meridian ratio Q.(4)R3θ=21−Q2cosθ+2Q−141−Q2cos2θ+5Q2−4Q41−Q2cos2θ+2Q−12.In formula (4), θ=cos−133/2J3/J23/2/3θ∈0,π/3； J2 and J3 are deviatoric stress invariants; and Q is the meridian ratio of tension and compression, and in order to consider the effect of pressure on strength, it is expressed as follows:(5)Q=Qp∗=Q0+Bp∗0.51≤Q≤1.In formula (5), Q0 is the initial meridian ratio of tension and compression, and B is the material constant.


## 2.2. Elastic Ultimate Surface Equation

In the RHT model, the elastic limit surface equation is derived from the failure surface equation:(6)σel∗p,θ,ε˙=σeq∗FelasticFcapp∗,In formula (6), Felastic is the elastic scaling function and Fcapp∗ is the cap function, which can effectively solve the problem that materials are always in the elastic state without yielding under high hydrostatic pressure. The expressions for Felastic, Fcapp∗ are as follows:(7)Felastic=gc∗,p≥fc,el3,p+ft,el/3fc,el/3+ft,el/3gc∗+p−fc,el/3−ft,el/3−fc,el/3gt∗,−ft,el3<p<fc,el3,gt∗,p≤−ft,el3.In equation (7), ft,el is the uniaxial tensile elastic ultimate stress, fc,el is the uniaxial compressive elastic ultimate stress, gc∗=fc,el/fc is the compressive yield surface parameter, and gt∗=ft,el/ft is the tensile yield surface parameter.(8)Fcapp∗=1,p∗⩽13,1−p∗−1/3p0∗−1/32,13<p∗<p0∗,0,p⩾p0∗.In equation (8), p0 is the pressure at which the void material begins to crush.


## 2.3. Residual Strength Surface Equation

When the equivalent stress strength of the material is greater than the failure stress strength, the damage of the material begins to accumulate and enters the damage softening phase. Damage variableD is the ratio of cumulative equivalent plastic strain increment to final failure equivalent plastic strain. The plastic strain of the material at ultimate failure is as follows:(9)εpf=D1p∗−1−Dpt∗D2p∗⩾1−Dpt∗+εpmD11/D2εpmp∗<1−Dpt∗+εpmD11/D2,In formula (9), εpm is the minimum plastic strain at the time of material damage, and pt∗ is the destruction of the cutoff pressure. D1 and D2 are material parameters. In fact, when the pressure of the material exceeds the yield stress, the friction and confining pressure still exist in the broken part of the material, which makes it have the ability to resist shear. Therefore, the residual strength surface equation is introduced into the RHT model.(10)σr∗p∗=Afp∗nf,In formula (10), Af is the residual stress intensity parameter and nf is the residual stress intensity index. The equivalent stresses between the failure strength surface and the residual strength surface are obtained by linear interpolation.


## 3. Improvement of the RHT Model

### 3.1. Correction of Strain Rate Enhancement Factor

In the RHT model, the dynamic strain rate enhancement factor (DIF) is sensitive to the high strain rate region of materials in tensile, which is different from most impact tests. Qi and Qian [29] research suggests the following: the strength of the material is not always infinitely enhanced with the increase in strain rate, when the strain rate is in the high strain rate region, and if the strain rate continues to increase, the strength of the material will slowly increase. Malvar and Ross [30] proposed a modified CEB model, which is consistent with the results of most experiments. However, the model still cannot show the characteristics of the explosion impact problem under a high strain rate. Based on this, this study uses the hyperbolic function (tanh) proposed by Gebbeken and Greulich [31] to describe the tensile strain rate enhancement factor of rock materials under a high strain rate. This function divides the strength change of rock into low strain rate region, medium strain rate region, and high strain rate region. As shown in Figure 1, the modified formula DIFt is as follows:(11)Ftε˙=ftdft=tanhlgε˙/ε˙0−WxSFmWy−1+1Wy.Figure 1
Law of rock material with tensile strain rate changes.In formula (11), ε˙0 is the reference strain rate. According to the data of Tedesco and Ross [32] experiments, the fitting formula (11) gives Fm=10, Wx=1.6, S=0.8, and Wy=5.5


### 3.2. Model Modification Considering Initial Damage

In the RHT model, the initial damage value of brittle materials such as rock is 0 by default, but in practical engineering, the rock will be affected by different degrees of excavation disturbance, weathering, and groundwater dissolution, which leads to the initial damage value of rock not being zero. Therefore, it is necessary to introduce the initial damage valueD0 to modify the damage accumulation expression of the RHT model, and the improved damage value D is as follows:(12)0≤D=D0+∑Δεpεpf=∫0εp1εpfdεp≤1.In formula (12), Δεp is equivalent plastic strain increment, D0=1−Vc¯2/Vc2 is the initial injury of the rock mass, Vc¯2 is the wave velocity of the rock mass, and Vc2 is the wave velocity of the rock test piece.


### 3.3. Correction of Residual Strength Surfaces

The Lode angle factor is considered in the failure surface equation, which can transform the compression meridian into the tension meridian. Drawing on this idea, the rudder angle factor can be taken into account in the calculation formula of residual stress. In the RHT model, the equivalent stress strength between the failure stress surface and the residual stress surface is obtained by linear interpolation, which is inconsistent with the actual stress state in the rock mass, because the hydrostatic pressure in this interval is usually constantly changed. Therefore, the normalized hydrostatic pressure in equation (10) can be corrected to the normalized hydrostatic pressure ppre∗ corresponding to the current stress state surface, and finally, the corrected residual surface strength can be obtained as follows:(13)σr−pre∗p∗=Afppre∗nfR3θ.When the strain exceeds the strain corresponding to the yield stress, it can be assumed that there is a linear relationship between the strain increment and the plastic strain increment.(14)ΔεP=Δε⋅∑ΔεD1ppre∗−1−Dpt∗D2−1.From formula (14), it can be inferred that(15)Dpre=∫0εpre−εfaildεD1ppre∗−1−Dpt∗D2×∫0εmax−εfaildεD1ppre∗−1−Dpt∗D2−1.In equation (15), εfail is the strain corresponding to the residual strength and εmax is the strain corresponding to the failure strength when a brittle material such as rock loses its load-bearing capacity, and the new equivalent stress intensity between the failure stress surface and the residual stress surface can be obtained by linear interpolation:(16)σdamage=1−Dpreσfail+Dpreσr−pre∗.


## 3.1. Correction of Strain Rate Enhancement Factor

In the RHT model, the dynamic strain rate enhancement factor (DIF) is sensitive to the high strain rate region of materials in tensile, which is different from most impact tests. Qi and Qian [29] research suggests the following: the strength of the material is not always infinitely enhanced with the increase in strain rate, when the strain rate is in the high strain rate region, and if the strain rate continues to increase, the strength of the material will slowly increase. Malvar and Ross [30] proposed a modified CEB model, which is consistent with the results of most experiments. However, the model still cannot show the characteristics of the explosion impact problem under a high strain rate. Based on this, this study uses the hyperbolic function (tanh) proposed by Gebbeken and Greulich [31] to describe the tensile strain rate enhancement factor of rock materials under a high strain rate. This function divides the strength change of rock into low strain rate region, medium strain rate region, and high strain rate region. As shown in Figure 1, the modified formula DIFt is as follows:(11)Ftε˙=ftdft=tanhlgε˙/ε˙0−WxSFmWy−1+1Wy.Figure 1
Law of rock material with tensile strain rate changes.In formula (11), ε˙0 is the reference strain rate. According to the data of Tedesco and Ross [32] experiments, the fitting formula (11) gives Fm=10, Wx=1.6, S=0.8, and Wy=5.5


## 3.2. Model Modification Considering Initial Damage

In the RHT model, the initial damage value of brittle materials such as rock is 0 by default, but in practical engineering, the rock will be affected by different degrees of excavation disturbance, weathering, and groundwater dissolution, which leads to the initial damage value of rock not being zero. Therefore, it is necessary to introduce the initial damage valueD0 to modify the damage accumulation expression of the RHT model, and the improved damage value D is as follows:(12)0≤D=D0+∑Δεpεpf=∫0εp1εpfdεp≤1.In formula (12), Δεp is equivalent plastic strain increment, D0=1−Vc¯2/Vc2 is the initial injury of the rock mass, Vc¯2 is the wave velocity of the rock mass, and Vc2 is the wave velocity of the rock test piece.


## 3.3. Correction of Residual Strength Surfaces

The Lode angle factor is considered in the failure surface equation, which can transform the compression meridian into the tension meridian. Drawing on this idea, the rudder angle factor can be taken into account in the calculation formula of residual stress. In the RHT model, the equivalent stress strength between the failure stress surface and the residual stress surface is obtained by linear interpolation, which is inconsistent with the actual stress state in the rock mass, because the hydrostatic pressure in this interval is usually constantly changed. Therefore, the normalized hydrostatic pressure in equation (10) can be corrected to the normalized hydrostatic pressure ppre∗ corresponding to the current stress state surface, and finally, the corrected residual surface strength can be obtained as follows:(13)σr−pre∗p∗=Afppre∗nfR3θ.When the strain exceeds the strain corresponding to the yield stress, it can be assumed that there is a linear relationship between the strain increment and the plastic strain increment.(14)ΔεP=Δε⋅∑ΔεD1ppre∗−1−Dpt∗D2−1.From formula (14), it can be inferred that(15)Dpre=∫0εpre−εfaildεD1ppre∗−1−Dpt∗D2×∫0εmax−εfaildεD1ppre∗−1−Dpt∗D2−1.In equation (15), εfail is the strain corresponding to the residual strength and εmax is the strain corresponding to the failure strength when a brittle material such as rock loses its load-bearing capacity, and the new equivalent stress intensity between the failure stress surface and the residual stress surface can be obtained by linear interpolation:(16)σdamage=1−Dpreσfail+Dpreσr−pre∗.


## 4. Crystalline Limestone RHT Parameter Determination

### 4.1. Parameter Determination of the Standard RHT Model

#### 4.1.1. Determination of Static Load Mechanic Parameters

The rock sample used in the test is crystalline limestone in a mining area in Guizhou province. The size of the specimen isΦ50 × 50 mm cylinder, the aspect ratio is 1.0, and core taking, cutting, and polishing of rock samples are carried out by core-taking machine, cutting machine, double-face grinding machine, and sandpaper. The nonparallelism and nonperpendicularity of rock specimens are limited below 0.02 mm in order to make the subsequent tests go on smoothly. The rock sample is shown in Figure 2.Figure 2
Experimental samples of crystalline limestone. (a) Standard sample size. (b) Partial rock specimen.
(a)(b)Thefc and ft of crystalline limestone can be obtained by conventional uniaxial compressive test, Brazilian splitting test, and rock wave velocity test. The rock density ρ0 and initial porosity α0 were measured by the balance method and the mass method. The mechanical parameters of crystalline limestone are shown in Table 1:Table 1
Static load mechanic parameters.
Type of rockfc (MPa)ft (MPa)EGPaμρ0g⋅cm−3α0vpm⋅s−1vsm⋅s−1GGPaCrystal limestone44.25.3447.50.362.681.0254235243917.46


#### 4.1.2.p−α State Equation Parameter Determination

Brittle materials such as rock and concrete contain a large number of voids in the interior. When the material is under a strong dynamic load, the material will be subjected to the combined action of shear stress and high hydrostatic pressure, which makes the mechanical response of the material more complex. Based on this, Herrman [33] puts forward an equation of state considering the internal voids of brittle materials in 1969, which is called equation p−α, and its compression process is shown in Figure 3:Figure 3
Pressure versus density diagram.From Figure3, it can be seen that when the pressure value is p=0, the initial density of the brittle material is ρ0=ρinitial, and when the pressure linearly increases to pcursh, then the density of the material is compacted to ρp. When the pressure increases in a nonlinear way until the internal voids of the rock are completely compressed, the corresponding pressure is plock and compaction density is ρlock. Generally, the dense relationship is expressed in Mie-Greisen [34] form as follows:(17)A1=ρ0c02.(18)A2=A12s−1.(19)A3=A13s2−4s+1.(20)B0=B1=2s−1.(21)T1=A1T2=0.In the above formula,A1, A2, and A3 are the Rankine-Hugoniot polynomial coefficients. c0 is the sound velocity in rock mass when the pressure is zero, and s is an empirical parameter. Calcium is the main component of crystalline limestone. Available through Meyers [35] literature s=0.95, B0, B1, T1, and T2 are the parameters of the equation of state. According to the literature [36], the pressure at the beginning of crushing of the void can be obtained as pcrush=fc/3. Generally, the pressure at the time of crushing of the material is plock=6000MPa and the compression index is n=3


#### 4.1.3. Determination of RHT Constitutive Parameters

Hoek and Brown [37] empirical formula is used to estimate rock strength under different confining pressures, and its expression is as follows:(22)σ1=σ3+mbσcσ3+sσc2.σ1 and σ3 in equation (22) are the maximum and minimum principal stresses at the failure of the rock mass, respectively, mb reflects the soft and hard degrees of rock mass, mb=24, and s reflects the degree of fragmentation of the rock mass and is taken as s=1. σf=1/2σ1−σ32+σ2−σ32+σ3−σ12, and p=−σ1+σ+2σ3/3. The mechanical parameters of the rock mass under different loading conditions can be obtained according to equation (22), as shown in Table 2.Table 2
mechanical parameters of crystalline limestone under different confining pressures.
σ2=σ3/ MPaσ1 (MPa)p (MPa)σfp∗σf∗05.34−1.785.34−0.040.12044.214.7314.730.330.33590.2033.4085.200.761.9320172.2270.74152.221.603.4450284.51128.17234.512.905.31100428.68209.56328.684.747.44The Mohr-Coulomb criterion was chosen to calculate the shear strengthfs and the cohesion c of the rock mass.(23)fs=τ=σ1−σ32sin2θ,2θ=π/2+φ in equation (23), and φ is the internal friction angle of rock mass, which is determined by Mohr-Coulomb criterion:(24)σ1=σ31+sinφ1−sinφ+2cosφ1−sinφ.The data brought into Table2 are calculated as φ=36.6∘ and fs=17.69MPa.When the rock is under quasi-static loading conditions,Fr=1, ε˙=3.0E−6, and the strength of the rock under 5MPa and 100MPa confining pressures is selected to determine the failure surface parameters. The equations for the compression meridian are as follows:(25)A0.76−13+A−1/NN=1.93,A4.74−13+A−1/NN=7.44.According to equation (25), the failure surface parameters A=2.47 and N=0.714 can be obtained. So far, most of the parameters in the constitutive model have been determined. The determination of parameters B and Q0 is complicated, and the sensitivity analysis of Li et al. [26, 38] to the parameters of the RHT model shows that it has little influence on the calculation results. Therefore, the value determined by the Riedel experiment can be used. The RHT constitutive parameters of crystalline limestone are shown in Table 3.Table 3
Values of RHT constitutive parameters.
Parameter symbolsParameter descriptionValueParameter symbolsParameter descriptionValueρ0Density (g·cm−3)2.68NFailure surface index0.714pelPressure in gap compression (GPa)0.0147GShear modulus17.46pcompPressure in gap compression (GPa)6Q0Pull pressure meridian ratio0.6805AFailure surface parameters2.47αCompression strain rate index0.0262A1Rankine-Hugoniot coefficient (GPa)48.06δTensile strain rate index0.0311A2Rankine-Hugoniot coefficient (GPa)43.26ε˙0cReference compressive strain rate (ms−1)3 × 10−8A3Rankine-Hugoniot coefficient (GPa)4.44ε˙0tReference tensile strain rate (ms−1)3 × 10−8BLode angle-related parameters0.0105ε˙cFailure compression strain rate3 × 1022B0Equation of state parameters0.9ε˙tFailure tensile strain rate3 × 1022B1Equation of state parameters0.9εpmMinimum residual strain of damage0.01T1Equation of state parameters (GPa)48.06D1Initial damage parameters0.04T2Equation of state parameters (GPa)0D2Damage parameters1nPorosity index3ξShear modulus reduction factor0.5α0Initial porosity1.1884gc∗Compression yield surface parameters0.88fcUniaxial compressive strength (GPa)0.0442gt∗Tensile yield surface parameters0.72ft∗Relative tensile strength0.12AfResidual stress strength parameters1.62fs∗Relative shear strength0.40nfResidual stress strength index0.6


### 4.2. Parameter Determination of the Modified RHT Model

#### 4.2.1. SHPB Experiment

Hopkinson pressure bar test was carried out, and the average impact pressures of the projectile in the experiment were 0.2, 0.3, 0.4, 0.5, and 0.6 MPa, respectively. Based on the theory of one-dimensional elastic waves and the assumption of stress uniformity, the three-wave method is used to calculate the stress-strain and the average strain rates of the specimen based on the experimentally data measured, the strain rates corresponding to different shock pressures are 30.37 s−1, 49.64 s−1, 87.05 s−1, 124.75 s−1, and 138.18 s−1 respectively, and the resulting stress-strain curve is shown in Figure 4.Figure 4
Stress-strain curve under different impact pressure.In order to determine the elastic ultimate strength, failure strength, and residual strength of rock, it is necessary to smooth the stress-strain curves obtained from experiments. Elastic ultimate strength is the point at which the slope of the stress-strain curve begins to change, failure strength is the point at which the maximum stress value is on the stress-strain curve, and residual strength is the first turning point after the peak value of the stress-strain curve. The stress in the direction of the compression bar isσ1, the corresponding strain is εeff, and the strain rate is ε˙eff. The strength indexes are shown in Table 4.Table 4
State quantities such as elastic limit, failure strength, and residual strength in the SHPB test curve (unit: MPa).
The impact pressureMean strain rate (s−1)fc,elσfσrpp∗εeff (s−1)0.230.3746.261.0239.0220.340.463.32 × 10−30.349.6465.3682.4249.1927.460.622.94 × 10−30.487.0589.57109.1857.2436.390.822.82 × 10−30.5124.75111.12121.8160.2040.360.914.29 × 10−30.6138.18133.7148.5668.1149.521.124.02 × 10−3It can be seen from Figure3 that when the impact pressure is 0.2 MPa, the corresponding yield strength is 61.02 MPa, and with the gradual increase in impact pressure, the failure strength of the rock also gradually increases. When the impact pressure is 0.6 MPa, the failure strength of rock reaches 148.56 MPa, which is 2.43 times higher than the yield strength when the impact pressure is 0.2 MPa and 3.36 times higher than the uniaxial compressive strength. It can be seen that with the increase in impact pressure, the yield strength of the rock is also gradually improved.Based on the experimental results in Table4, it is possible to plot the relationship between σr∗ and p∗, the curve equation σr−pre∗=1.94ppre∗0.420.68+0.05ppre∗ is fitted by the modified residual surface equation of equation (13), and decision factor R2=0.954. The fitting curve is shown in Figure 4:From Figure5 we can obtain Af=1.94 and nf=0.42, A=2.95 and N=0.76 are calculated by taking p∗ corresponding to 0.2 MPa and 0.6 MPa, respectively, and the rest of the parameters are the same as those taken in Table 3.Figure 5
Fitted curve of improved normalized pressure versus residual strength.


## 4.1. Parameter Determination of the Standard RHT Model

### 4.1.1. Determination of Static Load Mechanic Parameters

The rock sample used in the test is crystalline limestone in a mining area in Guizhou province. The size of the specimen isΦ50 × 50 mm cylinder, the aspect ratio is 1.0, and core taking, cutting, and polishing of rock samples are carried out by core-taking machine, cutting machine, double-face grinding machine, and sandpaper. The nonparallelism and nonperpendicularity of rock specimens are limited below 0.02 mm in order to make the subsequent tests go on smoothly. The rock sample is shown in Figure 2.Figure 2
Experimental samples of crystalline limestone. (a) Standard sample size. (b) Partial rock specimen.
(a)(b)Thefc and ft of crystalline limestone can be obtained by conventional uniaxial compressive test, Brazilian splitting test, and rock wave velocity test. The rock density ρ0 and initial porosity α0 were measured by the balance method and the mass method. The mechanical parameters of crystalline limestone are shown in Table 1:Table 1
Static load mechanic parameters.
Type of rockfc (MPa)ft (MPa)EGPaμρ0g⋅cm−3α0vpm⋅s−1vsm⋅s−1GGPaCrystal limestone44.25.3447.50.362.681.0254235243917.46


### 4.1.2.p−α State Equation Parameter Determination

Brittle materials such as rock and concrete contain a large number of voids in the interior. When the material is under a strong dynamic load, the material will be subjected to the combined action of shear stress and high hydrostatic pressure, which makes the mechanical response of the material more complex. Based on this, Herrman [33] puts forward an equation of state considering the internal voids of brittle materials in 1969, which is called equation p−α, and its compression process is shown in Figure 3:Figure 3
Pressure versus density diagram.From Figure3, it can be seen that when the pressure value is p=0, the initial density of the brittle material is ρ0=ρinitial, and when the pressure linearly increases to pcursh, then the density of the material is compacted to ρp. When the pressure increases in a nonlinear way until the internal voids of the rock are completely compressed, the corresponding pressure is plock and compaction density is ρlock. Generally, the dense relationship is expressed in Mie-Greisen [34] form as follows:(17)A1=ρ0c02.(18)A2=A12s−1.(19)A3=A13s2−4s+1.(20)B0=B1=2s−1.(21)T1=A1T2=0.In the above formula,A1, A2, and A3 are the Rankine-Hugoniot polynomial coefficients. c0 is the sound velocity in rock mass when the pressure is zero, and s is an empirical parameter. Calcium is the main component of crystalline limestone. Available through Meyers [35] literature s=0.95, B0, B1, T1, and T2 are the parameters of the equation of state. According to the literature [36], the pressure at the beginning of crushing of the void can be obtained as pcrush=fc/3. Generally, the pressure at the time of crushing of the material is plock=6000MPa and the compression index is n=3


### 4.1.3. Determination of RHT Constitutive Parameters

Hoek and Brown [37] empirical formula is used to estimate rock strength under different confining pressures, and its expression is as follows:(22)σ1=σ3+mbσcσ3+sσc2.σ1 and σ3 in equation (22) are the maximum and minimum principal stresses at the failure of the rock mass, respectively, mb reflects the soft and hard degrees of rock mass, mb=24, and s reflects the degree of fragmentation of the rock mass and is taken as s=1. σf=1/2σ1−σ32+σ2−σ32+σ3−σ12, and p=−σ1+σ+2σ3/3. The mechanical parameters of the rock mass under different loading conditions can be obtained according to equation (22), as shown in Table 2.Table 2
mechanical parameters of crystalline limestone under different confining pressures.
σ2=σ3/ MPaσ1 (MPa)p (MPa)σfp∗σf∗05.34−1.785.34−0.040.12044.214.7314.730.330.33590.2033.4085.200.761.9320172.2270.74152.221.603.4450284.51128.17234.512.905.31100428.68209.56328.684.747.44The Mohr-Coulomb criterion was chosen to calculate the shear strengthfs and the cohesion c of the rock mass.(23)fs=τ=σ1−σ32sin2θ,2θ=π/2+φ in equation (23), and φ is the internal friction angle of rock mass, which is determined by Mohr-Coulomb criterion:(24)σ1=σ31+sinφ1−sinφ+2cosφ1−sinφ.The data brought into Table2 are calculated as φ=36.6∘ and fs=17.69MPa.When the rock is under quasi-static loading conditions,Fr=1, ε˙=3.0E−6, and the strength of the rock under 5MPa and 100MPa confining pressures is selected to determine the failure surface parameters. The equations for the compression meridian are as follows:(25)A0.76−13+A−1/NN=1.93,A4.74−13+A−1/NN=7.44.According to equation (25), the failure surface parameters A=2.47 and N=0.714 can be obtained. So far, most of the parameters in the constitutive model have been determined. The determination of parameters B and Q0 is complicated, and the sensitivity analysis of Li et al. [26, 38] to the parameters of the RHT model shows that it has little influence on the calculation results. Therefore, the value determined by the Riedel experiment can be used. The RHT constitutive parameters of crystalline limestone are shown in Table 3.Table 3
Values of RHT constitutive parameters.
Parameter symbolsParameter descriptionValueParameter symbolsParameter descriptionValueρ0Density (g·cm−3)2.68NFailure surface index0.714pelPressure in gap compression (GPa)0.0147GShear modulus17.46pcompPressure in gap compression (GPa)6Q0Pull pressure meridian ratio0.6805AFailure surface parameters2.47αCompression strain rate index0.0262A1Rankine-Hugoniot coefficient (GPa)48.06δTensile strain rate index0.0311A2Rankine-Hugoniot coefficient (GPa)43.26ε˙0cReference compressive strain rate (ms−1)3 × 10−8A3Rankine-Hugoniot coefficient (GPa)4.44ε˙0tReference tensile strain rate (ms−1)3 × 10−8BLode angle-related parameters0.0105ε˙cFailure compression strain rate3 × 1022B0Equation of state parameters0.9ε˙tFailure tensile strain rate3 × 1022B1Equation of state parameters0.9εpmMinimum residual strain of damage0.01T1Equation of state parameters (GPa)48.06D1Initial damage parameters0.04T2Equation of state parameters (GPa)0D2Damage parameters1nPorosity index3ξShear modulus reduction factor0.5α0Initial porosity1.1884gc∗Compression yield surface parameters0.88fcUniaxial compressive strength (GPa)0.0442gt∗Tensile yield surface parameters0.72ft∗Relative tensile strength0.12AfResidual stress strength parameters1.62fs∗Relative shear strength0.40nfResidual stress strength index0.6


## 4.1.1. Determination of Static Load Mechanic Parameters

The rock sample used in the test is crystalline limestone in a mining area in Guizhou province. The size of the specimen isΦ50 × 50 mm cylinder, the aspect ratio is 1.0, and core taking, cutting, and polishing of rock samples are carried out by core-taking machine, cutting machine, double-face grinding machine, and sandpaper. The nonparallelism and nonperpendicularity of rock specimens are limited below 0.02 mm in order to make the subsequent tests go on smoothly. The rock sample is shown in Figure 2.Figure 2
Experimental samples of crystalline limestone. (a) Standard sample size. (b) Partial rock specimen.
(a)(b)Thefc and ft of crystalline limestone can be obtained by conventional uniaxial compressive test, Brazilian splitting test, and rock wave velocity test. The rock density ρ0 and initial porosity α0 were measured by the balance method and the mass method. The mechanical parameters of crystalline limestone are shown in Table 1:Table 1
Static load mechanic parameters.
Type of rockfc (MPa)ft (MPa)EGPaμρ0g⋅cm−3α0vpm⋅s−1vsm⋅s−1GGPaCrystal limestone44.25.3447.50.362.681.0254235243917.46


## 4.1.2.p−α State Equation Parameter Determination

Brittle materials such as rock and concrete contain a large number of voids in the interior. When the material is under a strong dynamic load, the material will be subjected to the combined action of shear stress and high hydrostatic pressure, which makes the mechanical response of the material more complex. Based on this, Herrman [33] puts forward an equation of state considering the internal voids of brittle materials in 1969, which is called equation p−α, and its compression process is shown in Figure 3:Figure 3
Pressure versus density diagram.From Figure3, it can be seen that when the pressure value is p=0, the initial density of the brittle material is ρ0=ρinitial, and when the pressure linearly increases to pcursh, then the density of the material is compacted to ρp. When the pressure increases in a nonlinear way until the internal voids of the rock are completely compressed, the corresponding pressure is plock and compaction density is ρlock. Generally, the dense relationship is expressed in Mie-Greisen [34] form as follows:(17)A1=ρ0c02.(18)A2=A12s−1.(19)A3=A13s2−4s+1.(20)B0=B1=2s−1.(21)T1=A1T2=0.In the above formula,A1, A2, and A3 are the Rankine-Hugoniot polynomial coefficients. c0 is the sound velocity in rock mass when the pressure is zero, and s is an empirical parameter. Calcium is the main component of crystalline limestone. Available through Meyers [35] literature s=0.95, B0, B1, T1, and T2 are the parameters of the equation of state. According to the literature [36], the pressure at the beginning of crushing of the void can be obtained as pcrush=fc/3. Generally, the pressure at the time of crushing of the material is plock=6000MPa and the compression index is n=3


## 4.1.3. Determination of RHT Constitutive Parameters

Hoek and Brown [37] empirical formula is used to estimate rock strength under different confining pressures, and its expression is as follows:(22)σ1=σ3+mbσcσ3+sσc2.σ1 and σ3 in equation (22) are the maximum and minimum principal stresses at the failure of the rock mass, respectively, mb reflects the soft and hard degrees of rock mass, mb=24, and s reflects the degree of fragmentation of the rock mass and is taken as s=1. σf=1/2σ1−σ32+σ2−σ32+σ3−σ12, and p=−σ1+σ+2σ3/3. The mechanical parameters of the rock mass under different loading conditions can be obtained according to equation (22), as shown in Table 2.Table 2
mechanical parameters of crystalline limestone under different confining pressures.
σ2=σ3/ MPaσ1 (MPa)p (MPa)σfp∗σf∗05.34−1.785.34−0.040.12044.214.7314.730.330.33590.2033.4085.200.761.9320172.2270.74152.221.603.4450284.51128.17234.512.905.31100428.68209.56328.684.747.44The Mohr-Coulomb criterion was chosen to calculate the shear strengthfs and the cohesion c of the rock mass.(23)fs=τ=σ1−σ32sin2θ,2θ=π/2+φ in equation (23), and φ is the internal friction angle of rock mass, which is determined by Mohr-Coulomb criterion:(24)σ1=σ31+sinφ1−sinφ+2cosφ1−sinφ.The data brought into Table2 are calculated as φ=36.6∘ and fs=17.69MPa.When the rock is under quasi-static loading conditions,Fr=1, ε˙=3.0E−6, and the strength of the rock under 5MPa and 100MPa confining pressures is selected to determine the failure surface parameters. The equations for the compression meridian are as follows:(25)A0.76−13+A−1/NN=1.93,A4.74−13+A−1/NN=7.44.According to equation (25), the failure surface parameters A=2.47 and N=0.714 can be obtained. So far, most of the parameters in the constitutive model have been determined. The determination of parameters B and Q0 is complicated, and the sensitivity analysis of Li et al. [26, 38] to the parameters of the RHT model shows that it has little influence on the calculation results. Therefore, the value determined by the Riedel experiment can be used. The RHT constitutive parameters of crystalline limestone are shown in Table 3.Table 3
Values of RHT constitutive parameters.
Parameter symbolsParameter descriptionValueParameter symbolsParameter descriptionValueρ0Density (g·cm−3)2.68NFailure surface index0.714pelPressure in gap compression (GPa)0.0147GShear modulus17.46pcompPressure in gap compression (GPa)6Q0Pull pressure meridian ratio0.6805AFailure surface parameters2.47αCompression strain rate index0.0262A1Rankine-Hugoniot coefficient (GPa)48.06δTensile strain rate index0.0311A2Rankine-Hugoniot coefficient (GPa)43.26ε˙0cReference compressive strain rate (ms−1)3 × 10−8A3Rankine-Hugoniot coefficient (GPa)4.44ε˙0tReference tensile strain rate (ms−1)3 × 10−8BLode angle-related parameters0.0105ε˙cFailure compression strain rate3 × 1022B0Equation of state parameters0.9ε˙tFailure tensile strain rate3 × 1022B1Equation of state parameters0.9εpmMinimum residual strain of damage0.01T1Equation of state parameters (GPa)48.06D1Initial damage parameters0.04T2Equation of state parameters (GPa)0D2Damage parameters1nPorosity index3ξShear modulus reduction factor0.5α0Initial porosity1.1884gc∗Compression yield surface parameters0.88fcUniaxial compressive strength (GPa)0.0442gt∗Tensile yield surface parameters0.72ft∗Relative tensile strength0.12AfResidual stress strength parameters1.62fs∗Relative shear strength0.40nfResidual stress strength index0.6


## 4.2. Parameter Determination of the Modified RHT Model

### 4.2.1. SHPB Experiment

Hopkinson pressure bar test was carried out, and the average impact pressures of the projectile in the experiment were 0.2, 0.3, 0.4, 0.5, and 0.6 MPa, respectively. Based on the theory of one-dimensional elastic waves and the assumption of stress uniformity, the three-wave method is used to calculate the stress-strain and the average strain rates of the specimen based on the experimentally data measured, the strain rates corresponding to different shock pressures are 30.37 s−1, 49.64 s−1, 87.05 s−1, 124.75 s−1, and 138.18 s−1 respectively, and the resulting stress-strain curve is shown in Figure 4.Figure 4
Stress-strain curve under different impact pressure.In order to determine the elastic ultimate strength, failure strength, and residual strength of rock, it is necessary to smooth the stress-strain curves obtained from experiments. Elastic ultimate strength is the point at which the slope of the stress-strain curve begins to change, failure strength is the point at which the maximum stress value is on the stress-strain curve, and residual strength is the first turning point after the peak value of the stress-strain curve. The stress in the direction of the compression bar isσ1, the corresponding strain is εeff, and the strain rate is ε˙eff. The strength indexes are shown in Table 4.Table 4
State quantities such as elastic limit, failure strength, and residual strength in the SHPB test curve (unit: MPa).
The impact pressureMean strain rate (s−1)fc,elσfσrpp∗εeff (s−1)0.230.3746.261.0239.0220.340.463.32 × 10−30.349.6465.3682.4249.1927.460.622.94 × 10−30.487.0589.57109.1857.2436.390.822.82 × 10−30.5124.75111.12121.8160.2040.360.914.29 × 10−30.6138.18133.7148.5668.1149.521.124.02 × 10−3It can be seen from Figure3 that when the impact pressure is 0.2 MPa, the corresponding yield strength is 61.02 MPa, and with the gradual increase in impact pressure, the failure strength of the rock also gradually increases. When the impact pressure is 0.6 MPa, the failure strength of rock reaches 148.56 MPa, which is 2.43 times higher than the yield strength when the impact pressure is 0.2 MPa and 3.36 times higher than the uniaxial compressive strength. It can be seen that with the increase in impact pressure, the yield strength of the rock is also gradually improved.Based on the experimental results in Table4, it is possible to plot the relationship between σr∗ and p∗, the curve equation σr−pre∗=1.94ppre∗0.420.68+0.05ppre∗ is fitted by the modified residual surface equation of equation (13), and decision factor R2=0.954. The fitting curve is shown in Figure 4:From Figure5 we can obtain Af=1.94 and nf=0.42, A=2.95 and N=0.76 are calculated by taking p∗ corresponding to 0.2 MPa and 0.6 MPa, respectively, and the rest of the parameters are the same as those taken in Table 3.Figure 5
Fitted curve of improved normalized pressure versus residual strength.


## 4.2.1. SHPB Experiment

Hopkinson pressure bar test was carried out, and the average impact pressures of the projectile in the experiment were 0.2, 0.3, 0.4, 0.5, and 0.6 MPa, respectively. Based on the theory of one-dimensional elastic waves and the assumption of stress uniformity, the three-wave method is used to calculate the stress-strain and the average strain rates of the specimen based on the experimentally data measured, the strain rates corresponding to different shock pressures are 30.37 s−1, 49.64 s−1, 87.05 s−1, 124.75 s−1, and 138.18 s−1 respectively, and the resulting stress-strain curve is shown in Figure 4.Figure 4
Stress-strain curve under different impact pressure.In order to determine the elastic ultimate strength, failure strength, and residual strength of rock, it is necessary to smooth the stress-strain curves obtained from experiments. Elastic ultimate strength is the point at which the slope of the stress-strain curve begins to change, failure strength is the point at which the maximum stress value is on the stress-strain curve, and residual strength is the first turning point after the peak value of the stress-strain curve. The stress in the direction of the compression bar isσ1, the corresponding strain is εeff, and the strain rate is ε˙eff. The strength indexes are shown in Table 4.Table 4
State quantities such as elastic limit, failure strength, and residual strength in the SHPB test curve (unit: MPa).
The impact pressureMean strain rate (s−1)fc,elσfσrpp∗εeff (s−1)0.230.3746.261.0239.0220.340.463.32 × 10−30.349.6465.3682.4249.1927.460.622.94 × 10−30.487.0589.57109.1857.2436.390.822.82 × 10−30.5124.75111.12121.8160.2040.360.914.29 × 10−30.6138.18133.7148.5668.1149.521.124.02 × 10−3It can be seen from Figure3 that when the impact pressure is 0.2 MPa, the corresponding yield strength is 61.02 MPa, and with the gradual increase in impact pressure, the failure strength of the rock also gradually increases. When the impact pressure is 0.6 MPa, the failure strength of rock reaches 148.56 MPa, which is 2.43 times higher than the yield strength when the impact pressure is 0.2 MPa and 3.36 times higher than the uniaxial compressive strength. It can be seen that with the increase in impact pressure, the yield strength of the rock is also gradually improved.Based on the experimental results in Table4, it is possible to plot the relationship between σr∗ and p∗, the curve equation σr−pre∗=1.94ppre∗0.420.68+0.05ppre∗ is fitted by the modified residual surface equation of equation (13), and decision factor R2=0.954. The fitting curve is shown in Figure 4:From Figure5 we can obtain Af=1.94 and nf=0.42, A=2.95 and N=0.76 are calculated by taking p∗ corresponding to 0.2 MPa and 0.6 MPa, respectively, and the rest of the parameters are the same as those taken in Table 3.Figure 5
Fitted curve of improved normalized pressure versus residual strength.


## 5. Validation of the Modified RHT Model

### 5.1. Numerical Simulation

In order to verify the correctness of the parameters of the improved RHT model, the numerical model is established by ANSYS/LS-DYNA software according to the experimental process of 1 : 1, in which the incident rod is 2.0 m long, the transmission rod is 1.5 m long, and the diameter of the rod is 0.05 m. The keyword∗ LOAD_SEGMENT_SET is used to load the sinusoidal stress wave on the end face of the incident rod, which not only ensures the consistency of the waveform but also plays a role in simplifying the model. The model is meshed into 255,400 cells and 274,373 nodes using SOLID 164 cells. The contact between the bar and the specimen is set as surface-to-surface erosion contact, the contact stiffness is the default value, and the dynamic and static friction coefficients are zero. The finite element analysis model is shown in Figure 6. The incident stress wave curve is selected when the impact pressure is 0.2 MPa for loading. When t=60.373μs, the peak value of incident stress is 134.26 MPa, and the stress wave curve is shown in Figure 7.Figure 6
Hopkinson compression bar model (a) and specimen meshing model (b).
(a)(b)Figure 7
Incident stress wave curve.


### 5.2. Analysis of Results

In order to verify the rationality of the improved RHT model, it is necessary to compare the difference between the rock failure patterns before and after the modification of constitutive parameters in numerical simulation and the experimental results, and the keyword∗MAT_ADD-EROSION in LS-DYNA is used to control the failure of the element in order to simulate the crack propagation in rock. The results of numerical simulation before correction and the experimental impact effect are shown in Figure 8.Figure 8
Numerical simulation results before modification of constitutive parameters of RHT.As can be seen from Figure8, the stress-strain curves obtained by numerical simulation show that the elastic ultimate strength is 6.98 MPa higher than the experimental value, and the yield strength is 1.45 MPA higher than the experimental value. At the same time, it is not difficult to find that the strain value corresponding to the elastic ultimate strength in the numerical simulation and test results is 5 × 10−4 more than the test value. When the strength exceeds the yield value, the residual strength surface described by the numerical simulation results is different from the experimental results. Although the surface crack propagation trend of numerical simulation is similar to the experimental results, only a few elements fail on the side of rock samples. This is because the unmodified constitutive parameters are conservative in describing the yield strength surface and residual strength surface, which leads to the failure mode of rock not fully displayed. In order to make the simulation results more consistent with the failure mode of rock, the improved constitutive parameters are used for calculation, and the results are shown in Figure 9.Figure 9
A comparison chart of numerical simulation effect after modification of constitutive parameters of RHT.From Figure9, it can be found that the description of elastic ultimate strength surface, yield strength surface, and residual failure surface by the improved RHT model is basically consistent with the test results. This is because the modified RHT model considers the constant change of hydrostatic pressure in the failure surface equation and the influence of the Lode angle factor in the residual strength surface equation, which makes the constitutive parameters of the RHT model more reasonable to describe the failure strength surface and the residual strength surface, which also verifies the correctness of the modified RHT model.


## 5.1. Numerical Simulation

In order to verify the correctness of the parameters of the improved RHT model, the numerical model is established by ANSYS/LS-DYNA software according to the experimental process of 1 : 1, in which the incident rod is 2.0 m long, the transmission rod is 1.5 m long, and the diameter of the rod is 0.05 m. The keyword∗ LOAD_SEGMENT_SET is used to load the sinusoidal stress wave on the end face of the incident rod, which not only ensures the consistency of the waveform but also plays a role in simplifying the model. The model is meshed into 255,400 cells and 274,373 nodes using SOLID 164 cells. The contact between the bar and the specimen is set as surface-to-surface erosion contact, the contact stiffness is the default value, and the dynamic and static friction coefficients are zero. The finite element analysis model is shown in Figure 6. The incident stress wave curve is selected when the impact pressure is 0.2 MPa for loading. When t=60.373μs, the peak value of incident stress is 134.26 MPa, and the stress wave curve is shown in Figure 7.Figure 6
Hopkinson compression bar model (a) and specimen meshing model (b).
(a)(b)Figure 7
Incident stress wave curve.


## 5.2. Analysis of Results

In order to verify the rationality of the improved RHT model, it is necessary to compare the difference between the rock failure patterns before and after the modification of constitutive parameters in numerical simulation and the experimental results, and the keyword∗MAT_ADD-EROSION in LS-DYNA is used to control the failure of the element in order to simulate the crack propagation in rock. The results of numerical simulation before correction and the experimental impact effect are shown in Figure 8.Figure 8
Numerical simulation results before modification of constitutive parameters of RHT.As can be seen from Figure8, the stress-strain curves obtained by numerical simulation show that the elastic ultimate strength is 6.98 MPa higher than the experimental value, and the yield strength is 1.45 MPA higher than the experimental value. At the same time, it is not difficult to find that the strain value corresponding to the elastic ultimate strength in the numerical simulation and test results is 5 × 10−4 more than the test value. When the strength exceeds the yield value, the residual strength surface described by the numerical simulation results is different from the experimental results. Although the surface crack propagation trend of numerical simulation is similar to the experimental results, only a few elements fail on the side of rock samples. This is because the unmodified constitutive parameters are conservative in describing the yield strength surface and residual strength surface, which leads to the failure mode of rock not fully displayed. In order to make the simulation results more consistent with the failure mode of rock, the improved constitutive parameters are used for calculation, and the results are shown in Figure 9.Figure 9
A comparison chart of numerical simulation effect after modification of constitutive parameters of RHT.From Figure9, it can be found that the description of elastic ultimate strength surface, yield strength surface, and residual failure surface by the improved RHT model is basically consistent with the test results. This is because the modified RHT model considers the constant change of hydrostatic pressure in the failure surface equation and the influence of the Lode angle factor in the residual strength surface equation, which makes the constitutive parameters of the RHT model more reasonable to describe the failure strength surface and the residual strength surface, which also verifies the correctness of the modified RHT model.


## 6. Conclusions

In this study, the shortcomings of the RHT constitutive model are studied, and the conclusions are as follows:(1)
It is practical to use the hyperbolic function to describe the tensile strain rate enhancement factor of rock at high strain rates.(2)
The introduction of the initial damage variableD0 makes the expression of damage variable D in the RHT constitutive model more complete.(3)
Considering the influence of the Lode angle factor and hydrostatic pressure, the updated damage variableDpre can better reflect the actual failure mode of the rock mass.(4)
It is a simple and effective method to determine the constitutive parameters of the RHT model by SHPB experiment.



---

*Source: 3794898-2022-04-13.xml*