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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 ${\sigma }_{\text{eq}}^{\ast }$ of the failure surface is a function of normalized pressure $p^{\ast }=p/f_c$, Lode angle $\theta$, and strain rate $\dot{\epsilon }$.$$\begin{array}{}\text{(1)}& {\sigma }_{\text{eq}}^{\ast }p,\theta ,\dot{\epsilon }=Y_{\text{TXC}}^{\ast }p{R}_3\theta F_{\text{rate}}\dot{\epsilon }.\end{array}$$

In formula (1), the ${}^{\ast }$ in the upper right corner of the formula indicates the normalized uniaxial compressive strength $f_c$ (unit: $\text{MPa}$) of the material, and $Y_{\text{TXC}}^{\ast }p$ is the compression meridian strength.$$\begin{array}{}\text{(2)}& Y_{\text{TXC}}{p}^{\ast }=A{p^{\ast }-p_{\text{spall}}^{\ast }{F}_{\text{rate}}\dot{\epsilon }}^N.\end{array}$$

In formula (2), $A$ is the failure surface constant; $N$ is the failure surface index; and $p_{\text{spall}}^{\ast }$ is the normalized layer crack strength.

  $F_{\text{rate}}\dot{\epsilon }$ is the strain rate dynamic enhancement factor:

$$\begin{array}{}\text{(3)}& F_{\text{rate\hspace{0.17em}}}\dot{\epsilon }=\begin{array}{cc}{\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0^c}}^{\alpha },\hfill & p\ge \frac{f_c}{3},\hfill \\ \frac{p+f_t/3}{f_c/3+f_t/3}{\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0^c}}^{\alpha }+\frac{p-f_c/3}{-f_t/3-f_c/3}{\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0^c}}^{\delta },\hfill & -\frac{f_t}{3}<p<\frac{f_c}{3}\hfill \\ {\frac{\dot{\epsilon }}{{\dot{\epsilon }}_0^t}}^{\delta }.\hfill & p\le -\frac{f_t}{3}.\hfill \end{array},\end{array}$$

In formula (3), $f_t$ is uniaxial tensile strength (unit: $\text{MPa}$)； ${\dot{\epsilon }}_0^c=30\times {10}^{-6}{s}^{-1}$； ${\dot{\epsilon }}_0^t=3\times {10}^{-6}{s}^{-1}$； $\dot{\epsilon }$ is the current strain rate; $\alpha =4/20+3f_c$ is a compression strain rate index; and $\delta =2/20+f_c$ is the tensile strain rate index.

$R_3\theta$ is the Lode angle factor: it is a function of the Lode angle $\theta$ and the tension-compression meridian ratio $Q$.$$\begin{array}{}\text{(4)}& R_3\theta =\frac{21-Q^2\cos \theta +2Q-1\sqrt{41-Q^2{\cos }^2\theta +5Q^2-4Q}}{41-Q^2{\cos }^2\theta +{2Q-1}^2}.\end{array}$$

In formula (4), $\theta ={\cos }^{-1}3\sqrt{3}/2J_3/J_2^{3/2}/3\hspace{1em}\theta \in 0,\pi /3$； $J_2$ and $J_3$ 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:$$\begin{array}{}\text{(5)}& Q=Q{p}^{\ast }=Q_0+B{p}^{\ast }0.51\le Q\le 1.\end{array}$$

In formula (5), $Q_0$ 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:$$\begin{array}{}\text{(6)}& {\sigma }_{\text{el}}^{\ast }p,\theta ,\dot{\epsilon }={\sigma }_{\text{eq}}^{\ast }{F}_{\text{elastic\hspace{0.17em}}}{F}_{\text{cap\hspace{0.17em}}}{p}^{\ast },\end{array}$$

In formula (6), $F_{\text{elastic\hspace{0.17em}}}$ is the elastic scaling function and $F_{\text{cap\hspace{0.17em}}}{p}^{\ast }$ 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 $F_{\text{elastic\hspace{0.17em}}}$, $F_{\text{cap\hspace{0.17em}}}{p}^{\ast }$ are as follows:$$\begin{array}{}\text{(7)}& F_{\text{elastic\hspace{0.17em}}}=\begin{array}{cc}{g}_c^{\ast },\hfill & p\ge \frac{f_{c,el}}{3},\hfill \\ \frac{p+f_{t,el}/3}{f_{c,el}/3+f_{t,el}/3}{g}_c^{\ast }+\frac{p-f_{c,el}/3}{-f_{t,el}/3-f_{c,el}/3}{g}_t^{\ast },\hfill & \frac{-f_{t,el}}{3}<p<\frac{f_{c,el}}{3},\hfill \\ g_t^{\ast },\hfill & p\le -\frac{f_{t,el}}{3}.\hfill \end{array}\end{array}$$

In equation (7), $f_{t,el}$ is the uniaxial tensile elastic ultimate stress, $f_{c,el}$ is the uniaxial compressive elastic ultimate stress, $g_c^{\ast }=f_{c,el}/f_c$ is the compressive yield surface parameter, and $g_t^{\ast }=f_{t,el}/f_t$ is the tensile yield surface parameter.$$\begin{array}{}\text{(8)}& F_{\text{cap\hspace{0.17em}}}{p}^{\ast }=\begin{array}{cc}1,\hfill & p^{\ast }\mathrm{⩽}\frac{1}{3},\hfill \\ \sqrt{1-{\frac{p^{\ast }-1/3}{p_0^{\ast }-1/3}}^2,}\hfill & \frac{1}{3}<p^{\ast }<p_0^{\ast },\hfill \\ 0,\hfill & p\mathrm{⩾}{p}_0^{\ast }.\hfill \end{array}\end{array}$$

In equation (8), $p_0$ 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 variable $D$ 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:$$\begin{array}{}\text{(9)}& {\epsilon }_p^{\text{f}}=\begin{array}{cc}{D}_1{p^{\ast }-1-D{p}_{\text{t}}^{\ast }}^{D_2}\hfill & p^{\ast }\mathrm{⩾}1-D{p}_{\text{t}}^{\ast }+{\frac{{\epsilon }_{\text{p}}^m}{D_1}}^{{1/D}_2}\hfill \\ {\epsilon }_{\text{p}}^m\hfill & p^{\ast }<1-D{p}_{\text{t}}^{\ast }+{\frac{{\epsilon }_{\text{p}}^m}{D_1}}^{{1/D}_2}\hfill \end{array},\end{array}$$

In formula (9), ${\epsilon }_{\text{p}}^m$ is the minimum plastic strain at the time of material damage, and $p_{\text{t}}^{\ast }$ is the destruction of the cutoff pressure. $D_1$ and $D_2$ 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.$$\begin{array}{}\text{(10)}& {\sigma }_r^{\ast }{p}^{\ast }=A_f{p^{\ast }}^{n_f},\end{array}$$

In formula (10), $A_f$ is the residual stress intensity parameter and $n_f$ 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 $DI{F}_t$ is as follows:$$\begin{array}{}\text{(11)}& F_t\dot{\epsilon }=\frac{f_{td}}{f_t}=\tanh \lg \dot{\epsilon }/{\dot{\epsilon }}_0-W_{x}S\frac{F_m}{W_y}-1+1W_y.\end{array}$$

[figure(s) omitted; refer to PDF]

In formula (11), ${\dot{\epsilon }}_0$ is the reference strain rate. According to the data of Tedesco and Ross [32] experiments, the fitting formula (11) gives $F_m=10$, $W_x=1.6$, $S=0.8$, and $W_y=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 value $D_0$ to modify the damage accumulation expression of the RHT model, and the improved damage value $D$ is as follows:$$\begin{array}{}\text{(12)}& 0\le D=D_0+{\displaystyle \sum \frac{\Delta {\epsilon }_p}{{\epsilon }_p^{\text{f\hspace{0.17em}}}}}={\displaystyle {\int }_0^{{\epsilon }_p}\frac{1}{{\epsilon }_p^{\text{f\hspace{0.17em}}}}}d{\epsilon }_p\le 1.\end{array}$$

In formula (12), $\mathrm{\Delta }{\epsilon }_p$ is equivalent plastic strain increment, $D_0=1-{\overline{V_c}}^2/V_c^2$ is the initial injury of the rock mass, ${\overline{V_c}}^2$ is the wave velocity of the rock mass, and $V_c^2$ 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 $p_{pre}^{\ast }$ corresponding to the current stress state surface, and finally, the corrected residual surface strength can be obtained as follows:$$\begin{array}{}\text{(13)}& {\sigma }_{r-pre}^{\ast }{p}^{\ast }=A_f{p_{pre}^{\ast }}^{n_f}{R}_3\theta .\end{array}$$

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.$$\begin{array}{}\text{(14)}& \Delta {\epsilon }_P=\Delta \epsilon \cdot {{\displaystyle \sum \frac{\Delta \epsilon }{D_1{p_{pre}^{\ast }-1-D{p}_{\text{t}}^{\ast }}^{D_2}}}}^{-1}.\end{array}$$

From formula (14), it can be inferred that$$\begin{array}{}\text{(15)}& D_{\text{pre}}={\displaystyle {\int }_0^{{\epsilon }_{pre}-{\epsilon }_{fail}}\frac{\text{d}\epsilon }{D_1{p_{pre}^{\ast }-1-D{p}_{\text{t}}^{\ast }}^{D_2}}}\times {{\displaystyle {\int }_0^{{\epsilon }_{\max }-{\epsilon }_{fail}}\frac{\text{d}\epsilon }{D_1{p_{pre}^{\ast }-1-D{p}_{\text{t}}^{\ast }}^{D_2}}}}^{-1}.\end{array}$$

In equation (15), ${\epsilon }_{fail}$ is the strain corresponding to the residual strength and ${\epsilon }_{\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:$$\begin{array}{}\text{(16)}& {\sigma }_{\text{damage}}=1-D_{\text{pre}}{\sigma }_{\text{fail}}+D_{\text{pre}}{\sigma }_{r-pre}^{\ast }.\end{array}$$

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 $\mathrm{\Phi }$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(s) omitted; refer to PDF]

The $f_c$ and $f_t$ of crystalline limestone can be obtained by conventional uniaxial compressive test, Brazilian splitting test, and rock wave velocity test. The rock density ${\rho }_0$ and initial porosity ${\alpha }_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 rock	$f_c$ (MPa)	$f_t$ (MPa)	$E\text{GPa}$	$\mu$	${\rho }_{0}g\cdot {\text{cm}}^{-3}$	${\alpha }_0$	$v_p\text{m}\cdot {\text{s}}^{-1}$	$v_s\text{m}\cdot {\text{s}}^{-1}$	$G\text{GPa}$
Crystal limestone	44.2	5.34	47.5	0.36	2.68	1.025	4235	2439	17.46

4.1.2. $p-\alpha$ 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-\alpha$, and its compression process is shown in Figure 3:

[figure(s) omitted; refer to PDF]

From Figure 3, it can be seen that when the pressure value is $p=0$, the initial density of the brittle material is ${\rho }_0={\rho }_{initial}$, and when the pressure linearly increases to $p_{cursh}$, then the density of the material is compacted to ${\rho }_p$. When the pressure increases in a nonlinear way until the internal voids of the rock are completely compressed, the corresponding pressure is $p_{lock}$ and compaction density is ${\rho }_{lock}$. Generally, the dense relationship is expressed in Mie-Greisen [34] form as follows:$$\begin{array}{}\text{(17)}& A_1={\rho }_0{c}_0^2.\text{(18)}& A_2=A_{1}2s-1.\text{(19)}& A_3=A_{1}3s^2-4s+1.\text{(20)}& B_0=B_1=2s-1.\text{(21)}& T_1=A_1{T}_2=0.\end{array}$$

In the above formula, $A_1$, $A_2$, and $A_3$ are the Rankine-Hugoniot polynomial coefficients. $c_0$ 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$, $B_0$, $B_1$, $T_1$, and $T_2$ 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 $p_{crush}=f_c/3$. Generally, the pressure at the time of crushing of the material is $p_{lock}=6000\text{MPa}$ 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:$$\begin{array}{}\text{(22)}& {\sigma }_1={\sigma }_3+\sqrt{m_b{\sigma }_c{\sigma }_3+s{\sigma }_c^2}.\end{array}$$

${\sigma }_1$ and ${\sigma }_3$ in equation (22) are the maximum and minimum principal stresses at the failure of the rock mass, respectively, $m_b$ reflects the soft and hard degrees of rock mass, $m_b=24$, and $s$ reflects the degree of fragmentation of the rock mass and is taken as $s=1$. ${\sigma }_f=\sqrt{1/2{{\sigma }_1-{\sigma }_3}^2+{{\sigma }_2-{\sigma }_3}^2+{{\sigma }_3-{\sigma }_1}^2}$, and $p=-{\sigma }_1+\sigma {}_2+{\sigma }_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.

${\sigma }_2={\sigma }_3$/ $MPa$	${\sigma }_1$ (MPa)	$p$ (MPa)	${\sigma }_f$	$p^{\ast }$	${\sigma }_f^{\ast }$
0	5.34	−1.78	5.34	−0.04	0.12
0	44.2	14.73	14.73	0.33	0.33
5	90.20	33.40	85.20	0.76	1.93
20	172.22	70.74	152.22	1.60	3.44
50	284.51	128.17	234.51	2.90	5.31
100	428.68	209.56	328.68	4.74	7.44

The Mohr-Coulomb criterion was chosen to calculate the shear strength $f_s$ and the cohesion $c$ of the rock mass.$$\begin{array}{}\text{(23)}& f_s=\tau =\frac{{\sigma }_1-{\sigma }_3}{2}\sin 2\theta ,\end{array}$$

$2\theta =\pi /2+\phi$ in equation (23), and $\phi$ is the internal friction angle of rock mass, which is determined by Mohr-Coulomb criterion:$$\begin{array}{}\text{(24)}& {\sigma }_1=\frac{{\sigma }_{3}1+\sin \phi }{1-\sin \phi }+\frac{2\cos \phi }{1-\sin \phi }.\end{array}$$

The data brought into Table 2 are calculated as $\phi ={36.6}^{\circ }$ and $f_s=17.69\text{MPa}$.

When the rock is under quasi-static loading conditions, $F_r=1$, $\dot{\epsilon }=3.0E-6$, and the strength of the rock under $5\text{MPa}$ and $100\text{MPa}$ confining pressures is selected to determine the failure surface parameters. The equations for the compression meridian are as follows:$$\begin{array}{}\text{(25)}& \begin{array}{c}A{0.76-\frac{1}{3}+A^{-1/N}}^N=1.93,\hfill \\ A{4.74-\frac{1}{3}+A^{-1/N}}^N=\mathrm{7.44.}\hfill \end{array}\end{array}$$

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 $Q_0$ 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 symbols	Parameter description	Value	Parameter symbols	Parameter description	Value
${\rho }_0$	Density (g·cm−3)	2.68	$N$	Failure surface index	0.714
$p_{el}$	Pressure in gap compression (GPa)	0.0147	$G$	Shear modulus	17.46
$p{}_{comp}$	Pressure in gap compression (GPa)	6	$Q_0$	Pull pressure meridian ratio	0.6805
$A$	Failure surface parameters	2.47	$\alpha$	Compression strain rate index	0.0262
$A_1$	Rankine-Hugoniot coefficient (GPa)	48.06	$\delta$	Tensile strain rate index	0.0311
$A_2$	Rankine-Hugoniot coefficient (GPa)	43.26	${\dot{\epsilon }}_0^c$	Reference compressive strain rate (ms−1)	3 × 10−8
$A_3$	Rankine-Hugoniot coefficient (GPa)	4.44	${\dot{\epsilon }}_0^t$	Reference tensile strain rate (ms−1)	3 × 10−8
$B$	Lode angle-related parameters	0.0105	${\dot{\epsilon }}^c$	Failure compression strain rate	3 × 1022
$B_0$	Equation of state parameters	0.9	${\dot{\epsilon }}^t$	Failure tensile strain rate	3 × 1022
$B_1$	Equation of state parameters	0.9	${\epsilon }_p^m$	Minimum residual strain of damage	0.01
$T_1$	Equation of state parameters (GPa)	48.06	$D_1$	Initial damage parameters	0.04
$T_2$	Equation of state parameters (GPa)	0	$D_2$	Damage parameters	1
$n$	Porosity index	3	$\xi$	Shear modulus reduction factor	0.5
${\alpha }_0$	Initial porosity	1.1884	$g_c^{\ast }$	Compression yield surface parameters	0.88
$f_c$	Uniaxial compressive strength (GPa)	0.0442	$g_t^{\ast }$	Tensile yield surface parameters	0.72
$f_t^{\ast }$	Relative tensile strength	0.12	$A_f$	Residual stress strength parameters	1.62
$f_s^{\ast }$	Relative shear strength	0.40	$n_f$	Residual stress strength index	0.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⁻¹, 49.64 s⁻¹, 87.05 s⁻¹, 124.75 s⁻¹, and 138.18 s⁻¹ respectively, and the resulting stress-strain curve is shown in Figure 4.

[figure(s) omitted; refer to PDF]

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 ${\sigma }_1$, the corresponding strain is ${\epsilon }_{eff}$, and the strain rate is ${\dot{\epsilon }}_{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 pressure	Mean strain rate (s−1)	$f_{c,el}$	${\sigma }_f$	${\sigma }_r$	$p$	$p^{\ast }$	${\epsilon }_{eff}$ (s−1)
0.2	30.37	46.2	61.02	39.02	20.34	0.46	3.32 × 10−3
0.3	49.64	65.36	82.42	49.19	27.46	0.62	2.94 × 10−3
0.4	87.05	89.57	109.18	57.24	36.39	0.82	2.82 × 10−3
0.5	124.75	111.12	121.81	60.20	40.36	0.91	4.29 × 10−3
0.6	138.18	133.7	148.56	68.11	49.52	1.12	4.02 × 10−3

It can be seen from Figure 3 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 Table 4, it is possible to plot the relationship between ${\sigma }_r{}^{\ast }$ and $p^{\ast }$, the curve equation ${\sigma }_{r-pre}^{\ast }=1.94{p_{pre}^{\ast }}^{0.42}0.68+0.05p_{pre}^{\ast }$ is fitted by the modified residual surface equation of equation (13), and decision factor $R^2=0.954$. The fitting curve is shown in Figure 4:

From Figure 5 we can obtain $A_f=1.94$ and $n_f=0.42$, $A=2.95$ and $N=0.76$ are calculated by taking $p^{\ast }$ 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(s) omitted; refer to PDF]

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 $\ast$ 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\mu \text{s}$, the peak value of incident stress is 134.26 MPa, and the stress wave curve is shown in Figure 7.

[figure(s) omitted; refer to PDF]

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(s) omitted; refer to PDF]

As can be seen from Figure 8, 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⁻⁴ 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(s) omitted; refer to PDF]

From Figure 9, 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 variable $D_0$ 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 variable $D_{pre}$ 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.

Acknowledgments

The author’s research work was supported by the National Natural Science Foundation of China (no. 51664007).

References

[1] G. R. Johnson, W. H. Cook, "A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures," Engineering Fracture Mechanics, vol. 21, pp. 541-548, 1983.

[2] L. J. Malvar, J. E. Crawford, J. W. Wesevich, D. Simons, "A plasticity concrete material model for DYNA3D," International Journal of Impact Engineering, vol. 19 no. 9-10, pp. 847-873, DOI: 10.1016/s0734-743x(97)00023-7, 1997.

[3] W. Riedel, K. Thoma, S. Hiermaier, "Penetration of reinforced concrete by BETA-B-500 numerical analysis using a new macroscopic concrete model for hydrocodes," pp. 315-322, .

[4] P. Skoglund, "Simulation of concrete penetration in 2D and 3D with the RHT material model," Swedish Defense Research Agency, 2002.

[5] Z. Tu, Y. Lu, "Evaluation of typical concrete material models used in hydrocodes for high dynamic response simulations," International Journal of Impact Engineering, vol. 36 no. 1, pp. 132-146, DOI: 10.1016/j.ijimpeng.2007.12.010, 2009.

[6] Z. Tu, Y. Lu, "Modifications of RHT material model for improved numerical simulation of dynamic response of concrete," International Journal of Impact Engineering, vol. 37 no. 10, pp. 1072-1082, DOI: 10.1016/j.ijimpeng.2010.04.004, 2010.

[7] T. L. Ling, Y. T. Wang, D. S. Liu, "Application of modified RHT model in numerical simulation of rock blast response," Journal of Coal, vol. 43 no. S2, pp. 434-442, 2018.

[8] C. P. Kuang, X. K. Yuan, "Analysis and simulation of strength parameters of RHT concrete principal structure model," Quarterly Journal of Mechanics, vol. 33 no. 1, pp. 158-163, 2012.

[9] M. Abdel-Kader, "Modified settings of concrete parameters in RHT model for predicting the response of concrete panels to impact," International Journal of Impact Engineering, vol. 132, pp. 103312.1-103312.18, DOI: 10.1016/j.ijimpeng.2019.06.001, 2019.

[10] A. Pavlovic, C. Fragassa, A. Disic, "Comparative numerical and experimental study of projectile impact on reinforced concrete," Composites Part B: Engineering, vol. 108, pp. 122-130, DOI: 10.1016/j.compositesb.2016.09.059, 2017.

[11] D. S. Shi, G. Y. Shi, J. Y. Chen, "Improvement and numerical fitting of RHT model," Journal of Ningbo University (Science and Technology Edition), vol. 27 no. 03, pp. 93-96, 2014.

[12] Y. Wang, Research on Rock Blasting Fragmentation Model Based on RHT Ontology, 2015.

[13] Z. Wang, Y. Ni, J. Cao, "Progress in the study of dynamic mechanical properties of concrete under impact loading," Explosion and Impact, vol. 2005 no. 6, pp. 519-527, 2005.

[14] R. Zhang, Y. Ding, W. Tang, "Failure strength parameters of concrete HJC and RHT intrinsic structure models," Journal of High Pressure Physics, vol. 25 no. 1, pp. 15-22, 2011.

[15] A. By, A. Zl, A. Yc, "Mechanical and microstructural properties of recycling granite residual soil reinforced with glass fiber and liquid-modified polyvinyl alcohol polymer," Chemosphere, 2021.

[16] L. Wang, F. Guo, H. Yang, Y Wang, S Tang, "Comparison of fly ash, PVA fiber, MGO and shrinkage-reducing admixture on the frost resistance of face slab concrete via pore structural and fractal analysis," Fractals, vol. 29, 2020.

[17] L. Wang, R. Luo, W. Zhang, M. Jin, S. Tang, "Effects of fineness and content of phosphorus slag on cement hydration, permeability, pore structure and fractal dimension of concrete," Fractals, vol. 29,DOI: 10.1142/s0218348x21400041, 2020.

[18] H. Li, Research on the Theory of Rock RHT Model and the Method of Determining the Main Parameters, 2016.

[19] L. X. Xie, W. B. Lu, Q. B. Zhang, Q. H. Jiang, M. Chen, J. Zhao, "Analysis of damage mechanisms and optimization of cut blasting design under high in-situ stresses," Tunnelling and Underground Space Technology, vol. 66, pp. 19-33, DOI: 10.1016/j.tust.2017.03.009, 2017.

[20] C. Pan, L-X. Xie, X. Li, K. Liu, P.-F. Gao, L.-G. Tian, "Numerical investigation of the effect of eccentric decoupled charge structure on blasting-induced rock damage," Journal of Central South University, vol. 29, 2022.

[21] Z. Wang, H. Wang, J. Wang, N. Tian, "Finite element analyses of constitutive models performance in the simulation of blast-induced rock cracks," Computers and Geotechnics, vol. 135,DOI: 10.1016/j.compgeo.2021.104172, 2021.

[22] B. Yuan, Z. Li, Z. Zhao, H. Ni, Z. Su, Z. Li, "Experimental study of displacement field of layered soils surrounding laterally loaded pile based on transparent soil," Journal of Soils and Sediments, vol. 21 no. 4,DOI: 10.1007/s11368-021-03004-y, 2021.

[23] B. Yuan, M. Sun, Y. Wang, L. Zhai, Q. Luo, X. Zhang, "Full 3D displacement measuring system for 3D displacement field of soil around a laterally loaded pile in transparent soil," International Journal of Geomechanics, vol. 2019 no. 5, 2019.

[24] C. Xu, G. Ye, "Application of orthogonal experimental design for sensitivity analysis of numerical model parameters," Hydrogeology & Engineering Geology, vol. 4 no. 01, pp. 95-97, 2004.

[25] D. Liu, N. Liu, T. Gao, "Parameter sensitivity study of RHT model by applying orthogonal test method," Journal of Beijing University of Technology, vol. 39 no. 06, pp. 558-564, 2019.

[26] H. Li, Y. Chen, D. Liu, "Study on the sensitivity of main parameters and determination method of rock RHT model," Journal of Beijing University of Technology, vol. 38 no. 08, pp. 779-785, 2018.

[27] Z. Nie, Y. Peng, R. Chen, "Parameter sensitivity analysis of RHT model for rock-like materials under intrusion conditions," Vibration and Shock, vol. 40 no. 14, pp. 108-116, 2021.

[28] X. Jian, "Sensitivity analysis of RHT model parameters under the effect of explosion," Naval Electronics Engineering, vol. 39 no. 04, pp. 111-113+122, 2019.

[29] C. Qi, Q. Qian, "Physical mechanism of strain rate dependence of dynamic strength of brittle materials such as rocks," Journal of Rock Mechanics and Engineering, vol. 3 no. 02, pp. 177-181, 2003.

[30] L. J. Malvar, C. A. Ross, "Review of strain rate effects for concrete in tension," Materials Journal, vol. 95 no. 6, 1998.

[31] N. Gebbeken, S. Greulich, "A new material model for SFRC under high dynamic loadings," Proceedings of the 1 Ith International Symposium on Interaction of the Effects of Munitions with Structures, .

[32] J. W. Tedesco, J. C. Powell, C. A. Ross, "A strain-rate-dependent concrete material model for ADINA," Computers & Structures, vol. 64 no. 5,DOI: 10.1016/s0045-7949(97)00018-7, 1997.

[33] W. Herrmann, "Constitutive equation for the dynamic compaction of ductile porous materials," Journal of Applied Physics, vol. 40 no. 6, pp. 2490-2499, DOI: 10.1063/1.1658021, 1969.

[34] X. Zhang, The Point-of-Matter Method, 2013.

[35] M. A. Meyers, Dynamic Behavior of Materials, 1994.

[36] H. G. R. A. Johnson, "Computational constitutive model for glass subjected to large strains, high strain rates and high pressures," Journal of Applied Mechanics, vol. 78 no. 5, 2011.

[37] E. Hoek, E. T. Brown, Underground Excavations in Rock, 1980.

[38] M. Yu, L. He, L. Song, "Twin shear stress theory and its generalization," Ence in China Ser A, vol. 10, 1985.