# Numerical Simulation Study on the Strengthening Mechanism of Rock Materials under Impact Loads

Authors: Xin Liu; Kai Wang; Chunan Tang; Xikun Qian
Journal: Shock and Vibration
 (2022)

---

## Abstract

Based on an analysis of real stress in rock specimens, the mechanisms of inertia and cracks on the dynamic properties of concrete material are discussed in this paper. A dynamic enhancement factor is introduced on the basis of static mesodamage constitutive theory, and a rate-dependent three-dimensional rock dynamic constitutive model is established. The orthogonal design and range method are used to analyze the sensitivity of the dynamic damage constitutive parameters, and a differential evolution algorithm is used to invert the relevant parameters of the constitutive model. Based on the drop hammer impact model established by the rock dynamic damage constitutive model, the AE distribution, crack propagation mode, and load-displacement curve of the rock sample are obtained. The numerical simulation results show that rock specimens under dynamic loading have different failure forms, mechanical response mechanisms, and energy evolution than those under static loading. Under dynamic loading, the rock specimen forms more cracks and releases more elastic energy. Furthermore, it shows greater strength characteristics. Inhomogeneities and inertial effects of rock materials are the fundamental reasons for the increase in the dynamic strength of rock materials. A reduction in inertial effects is the main reason for strength reduction in rock samples with hole defects.


---

## Body

## 1. Introduction

Rock materials are typically inhomogeneous materials with significant strain rate effects [1–3]. The strength, damage mode, fracture surface, and energy characteristics of rock materials exhibit various damage characteristics under different strain rates. When the rock loading strain is low, i.e., when it reaches a static or quasistatic loading range, the rock material is damaged according to the least energy dissipation principle, with the crack expansion path along the weakest face of the rock [4, 5]. As the strain rate increases, the energy inside the rock rapidly accumulates, causing the material to reach its storage limit in a short period of time and then rapidly release its energy. In this case, the crack expansion of the rock material follows the fastest energy release principle [6, 7], and the specimen reaches a stable state in the shortest time.Li et al. [8] investigated the energy characteristics of rock specimens during dynamic damage with an SHPB device. It was shown that, for higher rock loading strain rates, more energy was absorbed, stored, and dissipated from the rock specimen, resulting in greater fragmentation. The fracture surfaces of rock materials also exhibit distinct morphologies at different strain rates. Liang et al. [9] used an electron microscope laser scanner to examine fractures in granite at different strain rates. The results showed that as the strain rate increased, the fracture in the rock changed from an along-crystal fracture to a deconvolution fracture, to a through-crystal fracture with a zigzag fracture surface, and finally to a through-crystal fracture with a flat fracture surface. The results of Zhang et al. [10] also show that, at high strain rates, the rock has a large rupture damage surface that is significantly different from the rupture damage surface under static loading. Hull et al. [11] used scanning electron microscopy to observe the fracture pattern of rock specimens at different strain rates and to investigate the damage mechanism of rock specimens at different strain rates. Li et al. [12] used a similar approach to analyze the dynamic damage mechanism of soft rocks. Liang et al. [13] studied the relationship between rock energy storage and the release mechanisms of rock fractures by loading rocks at low strain rates. The study showed that the higher the rock loading strain rate is, the more energy the rock absorbs and stores, the more severe the rock fragmentation is, and the higher the measured rock strength is. Therefore, when loading at high strain rates, multiple cracks can penetrate areas of higher material strength, which differs significantly from the way weak surface expansion occurs during static loading [14]. When cracks penetrate high strength areas at higher strain rates, the ability of the rock specimen to resist deformation improves, thus increasing the measured rock strength. In other words, the strength of the rock is related to the strength of the grain through which the crack penetrates [15–18], and the greater the strength of the grain through which the crack penetrates, the greater the strength of the rock. When the rock specimen is subjected to a high strain rate, the inertial force [19–22] is the main reason the dynamic strength of the rock specimen increases. This mechanism of action is also consistent with the “S“ pattern of change in the dynamic strength of the rock. When the rock is subjected to a range of static and quasistatic strain rates, crack expansion in the rock specimen follows the least energy dissipation principle; cracks form along the weak surface of the specimen, where there is little change in strength. When the rock specimen is subjected to a higher strain rate, the dynamic strength of the rock is mainly related to the length of the crack and the strength of the area [23–26] through which the crack penetrates; the greater the strength of the area, the greater the strength of the rock specimen.However, the contribution of the inertial force and the crack penetration to the strength of the rock are not invariant for different models and are both closely related to the structural form of the model. Therefore, this paper constructs a dynamic damage intrinsic structure model based on the fine-scale damage intrinsic structure. The relevant parameters of the model were inverted using a differential evolutionary algorithm to analyze the effects of different inhomogeneities that cause inertia effects and cracks on the dynamic tensile strength with a finite element model of the falling weight impact test. First, on the basis of previous studies [9, 10], the necessary theoretical analysis was carried out to investigate the analytical expression for the dynamic strength with inhomogeneous effects. Second, the damage modes, tensile strength, and energy characteristics of different inhomogeneous numerical models were investigated under dynamic and static loading. Various models, including conventional and specially designed models, were used. Type I crack extensions and their dynamic tensile strength calculation values were first investigated for a case without particle inclusions. Subsequently, type I crack extensions and their dynamic tensile strength values were calculated for cases with rounded hard particle inclusions and rounded holes, and the static and dynamic damage modes were compared and analyzed for rocks with particles. Furthermore, on the basis of previous studies [27–33], the necessary theoretical analysis was carried out to investigate the dynamic macroscopic strength analytical expressions due to inhomogeneity. The hard particle parameters used in this paper were from data on sand conglomerates collected during the construction of a large reservoir at Sinopec's Shengli Oilfield [34].


## 2. Numerical Calculation Model

### 2.1. Dynamic Damage Criterion considering Rate Effects

To study the dynamic constitutive behavior of rocks, a combination of approaches is usually used: the static damage factor and dynamic constitutive model, the dynamic damage factor and static constitutive model, and the dynamic damage factor and dynamic constitutive model. As an inhomogeneous material with damage, rocks are subjected to dynamic loading due to the interaction of rate effects and damage weakening effects.Numerous laboratory tests on rocks have shown that the strength and elastic moduli of rock materials have a linear relationship with the logarithm of the strain rate [35]. Zhao et al. [35] have shown through dynamic tests on various rock materials that the relationship between the strength and the strain rate conforms to a certain logarithmic relationship, as follows:(1)ηDIF=1+Alogε˙ε˙s,where ηDIF is the dynamic strain rate factor, ε˙ represents the dynamic strain rate, ε˙s represents the quasi-static strain rate, and A is a material parameter.In this paper, to establish a dynamic damage intrinsic model for rocks [34], the strengthening effect of the strain rate is considered to decouple the dynamic damage, and a dynamic damage constitutive model that considers the damage weakening effect and rate correlation of rock materials is proposed as follows:(2)σ=E01+Alnε˙ε˙sε<ε0σ=E0exp−ε−εsma1+Alnε˙ε˙sε<ε0,where E0 denotes the initial Young's modulus, ε0 represents the strain when the deformation reaches the elastic limit, and m represents the heterogeneity coefficient.


### 2.2. Differential Evolution Algorithm Based Parameter Inversion

In this paper, a numerical calculation model was established based on the falling weight impact test in the literature [28], as shown in Figure 1. The length lb, height b, and span distance l of the numerical model were 200 mm, 50 mm, and 180 mm, respectively. The thickness of the model was 5 mm. In the middle of the lower boundary of the rock specimen, a prefabricated crack with a length c of 5 mm was set perpendicular to the boundary.Figure 1
The numerical model.The rock specimen was fixed on both sides of the bottom surface, and an impact load consisting of a 4 kg hammer with a 2 m drop was applied to the central axis of the top surface. The peak stress wave was 43 kN, and the duration of the stress wave was 1.2 ms. The model boundaries were all free reflective boundaries. The number of meshes was reduced to ensure accuracy, and the rock specimens were meshed to a 0.5 mm square hexagonal mesh.Due to the large scale of the model, to speed up the calculation and make full use of the computing hardware, OpenMP technology was used to parallelize the finite element stiffness calculation and the unit stress calculations. In addition, a reasonable superrelaxation factor was used to speed up the iterative convergence of the system of equations through a large number of trial calculations. With these improvements, the software has very good computational efficiency and parallelism for some large-scale problems, significantly reducing the computational time and improving the computational efficiency.The selection of the numerical model parameters is crucial because they determine the accuracy of the calculation. In this paper, a differential evolution algorithm with fast convergence and no initial values was used to implement an optimized inverse analysis of the parameters of the falling weight impact model.The specific inversion steps are as follows:(1)
The numerical models were developed based on the literature, and the model parameters in the present constitutive model were identified as the optimization variables for the differential evolution algorithm: the heterogeneity coefficientm, the modulus of elasticity E, Poisson's ratio υ, the residual strength σtr, the maximum tensile strain coefficient μ, the dynamic strength enhancement factor η and the tensile to compressive strength ratio γ. The orthogonal experimental design method was used to construct the parameter combination scheme for the numerical simulation. Seven relevant parameters were used in this constitutive model, and the parameter selection of the orthogonal experimental scheme was arranged using the Q31(57) orthogonal table.(2)
A sensitivity analysis of the parameters was carried out using the extreme difference method. Thirty-one sets of orthogonal test protocols were used to calculate the ultimate loads and their corresponding displacement values under different parameter combinations using the falling hammer impact model; based on the difference between the numerical results and the test results, the sensitivity of the relevant parameters wasη>m>γ>E>μ>υ>σcr. The first five parameters were selected as the final objects to be inversely analyzed, and the values of the remaining parameters were given in accordance with the experimental parameters and empirical values.(3)
The final inverse analysis parameters were determined to be the optimization variables; they were used to prepare the difference evolution inverse analysis algorithm program in MATLAB, set the initial values and parameters, and determine the fitness function. In this paper, the ultimate load and its corresponding displacement value were used as the target values, and the mean squared difference of the calculated data of the numerical model was used as the fitness function(3)gx=1n∑i=1n−1lx−lili2+Fx−FiFi2.where F(x) and l(x) are the experimentally measured peak stress and displacement data, respectively;Fi and li are the calculated peak load and displacement, respectively; and n is the number of experimental data points.The parametric inverse analysis program in MATLAB was used, and the corresponding initial parameters were set. The number of optimization variables was 4, the population size QT was 100, the variation factor a was 0.7, the crossover factorc was 0.8, and the maximum number of iterations dmax was 1500. By attempting different parameters in the optimization search test, the ultimate load and peak displacement were obtained under different parameter combinations.(4)
The inversion parameter results were output, and the results were analyzed and validated.Table1 shows the optimization results of the model parameters. Figure 2 shows the load–displacement curves of the numerical model and the physical results for the falling weight impact test. The results show that the calculated inverse parameters are consistent with the results of the laboratory falling weight impact experiment and that the inverse results are reliable.Table 1
Result of the back analysis of parameters.
ParametersE(GPa)η (10-2)mμγValue33.523.5812Figure 2
Comparison between numerical simulation and laboratory experiment.To further verify the accuracy of the selected parameters of the numerical model, the vertical displacements under loading were extracted, and the inertial force and vertical displacement data calculated according to the equations in the paper were compared and analyzed with the experimental results in the literature. The comparison results are shown in Figure3, which indicates that the results of the numerical simulation and the laboratory experiments are in good agreement, and the results calculated by the parameter inversion are consistent with the physical test results in terms of both the vertical displacement variation pattern and the inertia load curve of the rock specimen, thus proving the reliability of the parameter simulation of the falling weight impact test obtained by inversion in this paper.Figure 3
Comparison between numerical simulation and laboratory experiment. (a) Inertia force; (b) vertical displacement.
(a)(b)


### 2.3. Numerical Models

To analyze the effects of macroscopic inhomogeneities caused by inclusions of granular materials on the dynamic tensile strength in different model structural situations, a numerical model was established, as shown in Figure4. A circular hole region in the specimen, with a radius R of 7.5 mm, was considered for type I crack extensions and their calculated dynamic tensile strength values for three cases: a structural material with no inhomogeneities, a material with round hard particle inclusions, and a material with a circular hole. The stress wave shown in Figure 5 was applied for 1.4 ms with a crest of 42.5 kN. After the calculation was completed, the tensile strengths of the different structural forms were calculated, and the static and dynamic damage modes were compared and analyzed for rocks with particles. The tensile strength of a three-point bent rock specimen beam can be determined by the following equation:(4)σT=3Pl2bd,where P is the peak loading value, and l, b and d are the span, height, and thickness of the three-point bending beam, respectively.Figure 4
The numerical model.Figure 5
Shock waveform diagram.


## 2.1. Dynamic Damage Criterion considering Rate Effects

To study the dynamic constitutive behavior of rocks, a combination of approaches is usually used: the static damage factor and dynamic constitutive model, the dynamic damage factor and static constitutive model, and the dynamic damage factor and dynamic constitutive model. As an inhomogeneous material with damage, rocks are subjected to dynamic loading due to the interaction of rate effects and damage weakening effects.Numerous laboratory tests on rocks have shown that the strength and elastic moduli of rock materials have a linear relationship with the logarithm of the strain rate [35]. Zhao et al. [35] have shown through dynamic tests on various rock materials that the relationship between the strength and the strain rate conforms to a certain logarithmic relationship, as follows:(1)ηDIF=1+Alogε˙ε˙s,where ηDIF is the dynamic strain rate factor, ε˙ represents the dynamic strain rate, ε˙s represents the quasi-static strain rate, and A is a material parameter.In this paper, to establish a dynamic damage intrinsic model for rocks [34], the strengthening effect of the strain rate is considered to decouple the dynamic damage, and a dynamic damage constitutive model that considers the damage weakening effect and rate correlation of rock materials is proposed as follows:(2)σ=E01+Alnε˙ε˙sε<ε0σ=E0exp−ε−εsma1+Alnε˙ε˙sε<ε0,where E0 denotes the initial Young's modulus, ε0 represents the strain when the deformation reaches the elastic limit, and m represents the heterogeneity coefficient.


## 2.2. Differential Evolution Algorithm Based Parameter Inversion

In this paper, a numerical calculation model was established based on the falling weight impact test in the literature [28], as shown in Figure 1. The length lb, height b, and span distance l of the numerical model were 200 mm, 50 mm, and 180 mm, respectively. The thickness of the model was 5 mm. In the middle of the lower boundary of the rock specimen, a prefabricated crack with a length c of 5 mm was set perpendicular to the boundary.Figure 1
The numerical model.The rock specimen was fixed on both sides of the bottom surface, and an impact load consisting of a 4 kg hammer with a 2 m drop was applied to the central axis of the top surface. The peak stress wave was 43 kN, and the duration of the stress wave was 1.2 ms. The model boundaries were all free reflective boundaries. The number of meshes was reduced to ensure accuracy, and the rock specimens were meshed to a 0.5 mm square hexagonal mesh.Due to the large scale of the model, to speed up the calculation and make full use of the computing hardware, OpenMP technology was used to parallelize the finite element stiffness calculation and the unit stress calculations. In addition, a reasonable superrelaxation factor was used to speed up the iterative convergence of the system of equations through a large number of trial calculations. With these improvements, the software has very good computational efficiency and parallelism for some large-scale problems, significantly reducing the computational time and improving the computational efficiency.The selection of the numerical model parameters is crucial because they determine the accuracy of the calculation. In this paper, a differential evolution algorithm with fast convergence and no initial values was used to implement an optimized inverse analysis of the parameters of the falling weight impact model.The specific inversion steps are as follows:(1)
The numerical models were developed based on the literature, and the model parameters in the present constitutive model were identified as the optimization variables for the differential evolution algorithm: the heterogeneity coefficientm, the modulus of elasticity E, Poisson's ratio υ, the residual strength σtr, the maximum tensile strain coefficient μ, the dynamic strength enhancement factor η and the tensile to compressive strength ratio γ. The orthogonal experimental design method was used to construct the parameter combination scheme for the numerical simulation. Seven relevant parameters were used in this constitutive model, and the parameter selection of the orthogonal experimental scheme was arranged using the Q31(57) orthogonal table.(2)
A sensitivity analysis of the parameters was carried out using the extreme difference method. Thirty-one sets of orthogonal test protocols were used to calculate the ultimate loads and their corresponding displacement values under different parameter combinations using the falling hammer impact model; based on the difference between the numerical results and the test results, the sensitivity of the relevant parameters wasη>m>γ>E>μ>υ>σcr. The first five parameters were selected as the final objects to be inversely analyzed, and the values of the remaining parameters were given in accordance with the experimental parameters and empirical values.(3)
The final inverse analysis parameters were determined to be the optimization variables; they were used to prepare the difference evolution inverse analysis algorithm program in MATLAB, set the initial values and parameters, and determine the fitness function. In this paper, the ultimate load and its corresponding displacement value were used as the target values, and the mean squared difference of the calculated data of the numerical model was used as the fitness function(3)gx=1n∑i=1n−1lx−lili2+Fx−FiFi2.where F(x) and l(x) are the experimentally measured peak stress and displacement data, respectively;Fi and li are the calculated peak load and displacement, respectively; and n is the number of experimental data points.The parametric inverse analysis program in MATLAB was used, and the corresponding initial parameters were set. The number of optimization variables was 4, the population size QT was 100, the variation factor a was 0.7, the crossover factorc was 0.8, and the maximum number of iterations dmax was 1500. By attempting different parameters in the optimization search test, the ultimate load and peak displacement were obtained under different parameter combinations.(4)
The inversion parameter results were output, and the results were analyzed and validated.Table1 shows the optimization results of the model parameters. Figure 2 shows the load–displacement curves of the numerical model and the physical results for the falling weight impact test. The results show that the calculated inverse parameters are consistent with the results of the laboratory falling weight impact experiment and that the inverse results are reliable.Table 1
Result of the back analysis of parameters.
ParametersE(GPa)η (10-2)mμγValue33.523.5812Figure 2
Comparison between numerical simulation and laboratory experiment.To further verify the accuracy of the selected parameters of the numerical model, the vertical displacements under loading were extracted, and the inertial force and vertical displacement data calculated according to the equations in the paper were compared and analyzed with the experimental results in the literature. The comparison results are shown in Figure3, which indicates that the results of the numerical simulation and the laboratory experiments are in good agreement, and the results calculated by the parameter inversion are consistent with the physical test results in terms of both the vertical displacement variation pattern and the inertia load curve of the rock specimen, thus proving the reliability of the parameter simulation of the falling weight impact test obtained by inversion in this paper.Figure 3
Comparison between numerical simulation and laboratory experiment. (a) Inertia force; (b) vertical displacement.
(a)(b)


## 2.3. Numerical Models

To analyze the effects of macroscopic inhomogeneities caused by inclusions of granular materials on the dynamic tensile strength in different model structural situations, a numerical model was established, as shown in Figure4. A circular hole region in the specimen, with a radius R of 7.5 mm, was considered for type I crack extensions and their calculated dynamic tensile strength values for three cases: a structural material with no inhomogeneities, a material with round hard particle inclusions, and a material with a circular hole. The stress wave shown in Figure 5 was applied for 1.4 ms with a crest of 42.5 kN. After the calculation was completed, the tensile strengths of the different structural forms were calculated, and the static and dynamic damage modes were compared and analyzed for rocks with particles. The tensile strength of a three-point bent rock specimen beam can be determined by the following equation:(4)σT=3Pl2bd,where P is the peak loading value, and l, b and d are the span, height, and thickness of the three-point bending beam, respectively.Figure 4
The numerical model.Figure 5
Shock waveform diagram.


## 3. Analysis of the Effect of Rock Inhomogeneity and Inertial Forces on Dynamic Strength

To determine the influence of inhomogeneities in rock materials on dynamic tensile strength, this paper integrates the effects of inhomogeneities and inertia on the dynamic strength based on the literature [9] and establishes a dynamic tensile strength calculation method based on the crack extension area and inertia effects. This method can be used to solve for the dynamic tensile strength when the crack passes through an inhomogeneous region. When a three-point bending beam is subjected to an external load, a type I fracture occurs in the middle section; the fracture area of the rupture surface is S = b×dc. When the three-point bending beam fractures completely, the fracture surface is the entire cross section, S0 = b×d. At this time, the specimen reaches a static fracture strength with an ultimate load P of σs = PL/bh2. When the crack passes through a region of strength σi, the crack area formed is Si. At this point, the magnitude of the strength can be expressed as(5)σ=σiSi+σsSsS0.According to the analysis of forces in concrete falling weight impact tests in the literature [10], when a rock is impacted by a falling weight, the forces can be divided into two parts: inertial forces Fit and deformation forces Fit(6)Ftt=Fdt+Fit,where Fit is the equivalent inertia force generated at the loading point.The damage process of a rock specimen subjected to an impact load consists mainly of elastic deformation and rock cracking stress relaxation. To simplify the analysis, the deflection curve of the rock specimen is assumed to be in the form of three deflection curves throughout the falling weight impact process, as shown in Figure6.Figure 6
The deflection curve of falling weight impact.When the falling weight impact is applied to the upper boundary of the specimen at the mid-axis point, a vertical inertial acceleration is generated. According to the principle of virtual work, the inertial force of the impact can be obtained in relation to the rock specimen as a whole as(7)Fidδ=2∫ρH¨x,ydHydV,where ρ is the density of the rock specimen, H(x, y) is the vertical displacement of the rock specimen, and V is the volume of the rock specimen.The expression for the inertial force of a falling weight impact can be obtained by integrating over the whole region of the monolithic rock formula while neglecting higher order infinitesimals:(8)Fit=64711340l3ρδ¨t.Taking into account the inhomogeneity of the rock and the inertial forces, an expression for the dynamic strength of the rock can be obtained:(9)σd=σ1S1+σ2S2+⋯+σnSnS0+Fil1bh.


## 4. Analysis of the Mechanism of Strength Enhancement by Rock Inhomogeneity and Inertial Forces

### 4.1. Static Tensile Strength for Inhomogeneity Numerical Models

The same displacement load of 0.01 mm per step was applied to the central axis of the upper boundary of the three models until a type I crack was generated in the middle section through the whole specimen; the damage process of the specimen is shown in Figure7. All three rock specimens formed fracture surfaces at 20 loading steps, and all fracture surfaces were located in the center of the bottom surface of the specimen. As the load was applied, the crack in Model I expanded along the central axis of the specimen, and the I-shaped tensile crack formed in Model I expanded along the middle section and truncated in the middle of the circular hole at step 50. At step 55, an I-shaped crack formed at the top of the hole and expanded to the central axis of the upper boundary of the specimen by step 80, as shown in Figure 7(a). The type I tensile crack formed in Model II extended along the central axis of the model to the top of the rock specimen and, by step 100, to the top of the model, causing integral instability damage to the specimen, as shown in Figure 7(b).Figure 7
Diagrams of damage process and minimum principal stress diagram under static load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)A short pause in crack extension occurred in Model III at step 50 when the extension reached the particle. Then, at step 55, the crack extended along the left-side edge of the hard particle; when the extension reached a radius length to the left of the center of the particle, the crack extended upward in a vertical direction to the upper boundary of the specimen within a radius length of the left of the center axis by step 100, as shown in Figure7(c).The numerical simulation results show that when the materials in the crack penetration area of the numerical model vary, there are some differences in the extended crack path and the damage process. When the intermediate region is a circular hole defect, the crack propagation is blocked, and the damage process of the specimen is the fastest. When the intermediate region is a hard particle with a higher strength than the matrix, the crack propagates first along the particle-matrix interface and then in the vertical direction.To study the acoustic emission pattern during the damage of the model under static loading [36–40], the distributions of the number of AEs and the energy with the loading step were extracted for each of the three models under static loading, as shown in Figure 8.Figure 8
Number of unit acoustic emissions and energy characteristics of the numerical models under static load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)The numerical simulation results show that the number of AEs and the energy of all three models were concentrated between loading steps 20 and 60. This is mainly due to the presence of hard particles in the crack expansion path in the model, which prevents the crack from penetrating the stronger particles and causes damage along the contact surface with the matrix, resulting in the destruction of the units and the release of more energy. Model I produced the lowest number of AEs and the smallest amount of energy, with values of 2071 and 1.13 J, respectively. This was mainly due to the presence of hole defects in the model, which caused the cracks to be truncated, resulting in fewer units breaking and less energy being released.Figure9 shows the load-displacement curves of Models I, II, and III under static loading. The simulation results show that the load peaks at 0.44 mm, 0.52 mm, and 0.66 mm were 74831 N, 94437 N, and 102852 N, respectively. The ultimate load and displacement were the largest for Model III, while the load and displacement were the smallest for Model I, which contained the circular hole defect. There was a rapid decrease in the stress, which was mainly due to the presence of hole defects in the model, and when the crack extended to the bottom of the round hole, a crack quickly formed at the top of the defect, reducing the resistance of the model. The ultimate load was substituted into (5) to calculate the tensile strength of the model and compared with the results of (4), as shown in Table 2.Figure 9
Load-displacement curve of model under static load.Table 2
Calculation results of tensile strength by different calculation methods.
Calculation methodModels IModels IIModels IIIEquation 48.07MPa10.21MPa11.93MPaEquation 57.93MPa10.01MPa10.90MPaThe calculations show that it is reasonable to use the calculated tensile strength of inhomogeneous materials of different strengths, with a maximum error of 0.95%.


### 4.2. Dynamic Tensile Strength of Inhomogeneity Numerical Models

A stress wave, as shown in Figure5, was applied to the central axis of the upper boundary of the three models to simulate the impact load of a 4 kg hammer dropped from 2 m, with a peak stress wave of 43 kN and a duration of 1.2 ms. Figure 10 shows the damage process of the numerical model under the action of the falling weight impact. The numerical simulation results show that the three different structural models exhibit various damage modes under dynamic loading.Figure 10
Diagrams of damage process and minimum principal stress diagram under dynamic load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)Type I fracture cracks appeared in all three models at the lower boundary of the mid-axis of the precast crack, with Models I and II showing crack initiation at the precast crack at the bottom of the specimen att = 0.5 ms and Model III only showing stress concentration at the tip of the precast crack at t = 0.5 ms, with no type I fracture crack initiation. This is because the stress waves in the circular hole in Model III were reflected through the hole, which prevented some of the stress waves from propagating, resulting in the horizontal tensile displacement in the central axis of Model III lagging behind the other two models. As the stress wave continued to propagate, the type I crack produced in Model I continued to expand along the mid-axis region of the model specimen toward the top of the specimen, reaching the middle of the model at the upper boundary at t = 2 ms. After a type I crack formed in Model II, the crack began to expand rapidly along the mid-axis region of the model, and at t = 1 ms, the crack expanded to the lower endpoint of the circular particle, where the crack bifurcated, and the two bifurcated cracks 1 and 2 expanded along the boundary of the circular hard particle. At t = 1.5 ms, the two bifurcated cracks stopped after expanding to the endpoints on both sides. A third, vertically oriented type I crack emerged at the lower endpoint of the rounded particle. At this point, cracks 1, 2, and 3 all expanded vertically. At t = 2 ms, cracks 1 and 2 stopped expanding and did not reach the upper end of the model, while crack 3 stopped expanding at the center of the round particle.Att = 1 ms, in Model III, the type I crack expansion stopped at the lower endpoint of the circular hole due to the presence of the circular hole, which hindered the propagation of the stress wave. At the same time, two type I cracks formed on the left and right sides of the prefabricated cracks at the lower boundary of the model, mainly due to tensile waves reflected from the lower boundary of the circular hole. At t = 1.5 ms, type I cracks formed in the upper region of the circular hole defect and, at t = 2 ms, extended to the upper boundary of the model, while the two type I cracks sprouted at the lower boundary stopped extending after a certain distance. The numerical models of the three different structures showed significant differences in the damage process under the falling weight impact, with distinct crack expansion paths and damage patterns due to the various materials in the middle circular region.Figure11 shows the spatial distribution of AEs for the three different structural numerical models. Model III has the most complex form of damage, with concentrated unit tensile damage at the left and right ends of the circular hole and shear damage units in the upper region of the hole, mainly due to the presence of the circular hole, which causes the stress waves to reflect back and forth in this region.Figure 11
The distribution of AE under dynamic load.Figure12 shows the distribution of the number of AEs and the amount of energy with time for Models I to III under dynamic loading. The numerical simulation results show that the distribution of AEs over time differs for the three models. For example, in Model I, the number of AEs gradually increases at t = 0.4 ms, reaches a peak at t = 1 ms, and then gradually decreases. This is mainly due to the AE from the vertical tensile crack at the bottom of the model. At t = 1.5 ms, a second peak in the AE occurs at the top of the circular hole in the model due to crack initiation. Model II has only one peak in its AE distribution due to its homogeneous structure, with only one type I crack produced.Figure 12
Number of unit acoustic emissions and energy characteristics of the numerical models under dynamic load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)In Model III, due to the presence of hard particles, the crack expands rapidly along the particle-matrix interface, causing the AE first peak to occur att = 0.8 ms; at t = 1.5 ms, as the peak stress wave gradually increases, the crack passes through the hard particle region, causing a second AE peak. The AE time distribution curve in Figure 12 shows that crack initiation causes an increase in the AE, while crack expansion and penetration cause a concentrated burst in the number of AEs. Furthermore, by observing the cumulative release energy pattern of the AEs, it is found that the released energy of the AEs increases with increasing loading step. The highest number of AEs and the highest cumulative release energy were 8068 and 6.58 J, respectively, for Model I. Model I had the most damaged units due to the presence of circular defects and multiple transverse reflections of stress waves inside the specimen, which resulted in more tensile cracks in the model.Compared with Model I, the number of AEs and the energy produced in Model III are 7898 and 4.62 J, respectively. Although the number of AEs is similar, there is a difference of 1.92 J in the energy released, which is mainly due to the presence of more shear damage in Model I, which causes shear damage at the edge of the defect, releasing more energy. It is worth noting that Model II shows an increase in the number of AEs and the amount of energy released at the point, where the crack reaches the lower endpoint of the round particle. In particular, the rapid increase in AE energy, when combined with the AE distribution diagram at this point, is found to be the main reason for the rapid increase in microcrack energy. The sheer damage is mainly due to the apparent stress concentration at the tip of the particle and crack, which causes shear damage to the unit in this region.To further investigate the force and deformation mechanisms of the different structural inhomogeneous models under dynamic loading, the load-displacement curves, inertial force curves, and vertical displacement curves with time for the three numerical models were calculated, as shown in Figure13.Figure 13
The results of the numerical model under dynamic load. (a) Load-displacement curve. (b) Model inertia force curve. (c) Vertical displacement curve.
(a)(b)(c)According to the load-displacement curve of the model in Figure13(a), the maximum ultimate load of the model with hard particles is 49.07 kN, the minimum ultimate load of the model with circular hole defects is 38.7 kN, and the maximum displacement of the ultimate load is 0.26 mm. It is worth mentioning that the presence of circular defects causes a significant load drop in the load-displacement curve, as shown in Figure 13(b). The three models also have different inertial effects under stress waves due to their structural inhomogeneities, which confirms the conclusions reached in Section 5 of this paper.Figure13(c) shows the vertical displacement over time for the three model monitoring points. The vertical displacement for the model with round particles has the fastest rate of change but the smallest value. In Model III, the rate of change of the vertical displacement is the smallest because it contains circular defects, but the ultimate displacement value is the largest. This is mainly due to the presence of circular defects, which obstruct the propagation of stress waves, causing significant hysteresis in crack initiation and expansion. The equivalent inertial forces calculated according to (8) are shown in Figure 13(b). The numerical simulation results show that the equivalent inertial force of Model III is the largest, while that of Model I is the smallest. Overall, there are no significant differences between Model II and Model III in any calculation aspect, but the two models do have considerable differences with Model I. The dynamic strength magnitudes derived from the three numerical models based on their rupture areas and (4) and (9) were calculated, as shown in Table 3.Table 3
Comparison of calculation methods for tensile strength of numerical models.
ModelModel IModel IIModel IIIInhomogeneity strength value13.12MPa22.35MPa15.20MPaInertia force Strength values29.09MPa30.24MPa25.53MPaTotal model Strength values42.21MPa52.59MPa40.73MPaTraditional calculation method44.26MPa54.65MPa41.22MPaError2.66％3.92％1.21％According to the calculation results in the table, the error between the method in this work and the traditional calculation method is within 4%, indicating that the method has good accuracy. It is worth noting that the calculated values of the method in this work are all less than the calculated values of the traditional method, which is mainly because the traditional calculation method only considers tensile damage, but there may be a small amount of shear damage during the dynamic loading process, resulting in an increase in the strength.


### 4.3. Mechanism of the Influence of the Radius of a Circular Hole on Strength

To investigate the effect of different defect radius sizes on rock specimen damage, three-point bent rock specimens with radiiR of 5 mm, 7.5 mm, 10 mm, 12.5 mm, and 15 mm were established. A shock wave, as shown in Figure 5, was applied to the central region of the upper boundary of the specimen. The final damage pattern of the rock specimen is shown in Figure 14.Figure 14
Failure pattern of different hole radius.As the radiusR of the circular hole increases, the pattern of the crack extension type of the rock specimen changes. When R = 5 mm, the crack extends from the top of the precast crack at the bottom of the specimen to the top of the specimen through the top and bottom of the circular hole. When R = 7.5 mm, two short tensile cracks form on both sides of the precast crack due to the stress waves reflecting at the hole boundary. As the hole diameter increases, the stress wave reflection effect increases, and the tensile cracks that form on either side of the precast crack expand toward the left and right sides of the hole. When R = 15 mm, the rock specimen shows a similar damage pattern to R = 12.5 mm, but the cracks start to sprout and expand as prefabricated cracks above the hole rather than at the hole. This indicates that when the diameter of the hole exceeds a certain value, the rock specimen no longer expands as a three-point bending type I crack.Table4 shows the contributions of inertia effects and crack extension to the dynamic tensile strength values for rock specimens with different hole radii, as calculated by the above equations. The calculation results show that the dynamic tensile strength of the rock gradually decreases as the hole radius increases, and the larger the radius, the more significant the decrease. At this point, although the strength of the specimen gradually decreases, the contribution of the crack extension area to the strength first decreases, then increases, and finally decreases. The influence of inertia effects during loading is always greater than the influence of inhomogeneities. In other words, the hole affects the integral tensile strength of the specimen mainly by reducing the inertia effects of the specimen.Table 4
Model strength of different hole radius.
Hole radius (mm)Inhomogeneous strength value (MPa)Inertial forces (MPa)Total strength (MPa)013.1229.0942.21512.9528.9241.877.515.2025.5340.731017.2521.2838.5312.517.7318.1235.851513.1215.3528.47


## 4.1. Static Tensile Strength for Inhomogeneity Numerical Models

The same displacement load of 0.01 mm per step was applied to the central axis of the upper boundary of the three models until a type I crack was generated in the middle section through the whole specimen; the damage process of the specimen is shown in Figure7. All three rock specimens formed fracture surfaces at 20 loading steps, and all fracture surfaces were located in the center of the bottom surface of the specimen. As the load was applied, the crack in Model I expanded along the central axis of the specimen, and the I-shaped tensile crack formed in Model I expanded along the middle section and truncated in the middle of the circular hole at step 50. At step 55, an I-shaped crack formed at the top of the hole and expanded to the central axis of the upper boundary of the specimen by step 80, as shown in Figure 7(a). The type I tensile crack formed in Model II extended along the central axis of the model to the top of the rock specimen and, by step 100, to the top of the model, causing integral instability damage to the specimen, as shown in Figure 7(b).Figure 7
Diagrams of damage process and minimum principal stress diagram under static load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)A short pause in crack extension occurred in Model III at step 50 when the extension reached the particle. Then, at step 55, the crack extended along the left-side edge of the hard particle; when the extension reached a radius length to the left of the center of the particle, the crack extended upward in a vertical direction to the upper boundary of the specimen within a radius length of the left of the center axis by step 100, as shown in Figure7(c).The numerical simulation results show that when the materials in the crack penetration area of the numerical model vary, there are some differences in the extended crack path and the damage process. When the intermediate region is a circular hole defect, the crack propagation is blocked, and the damage process of the specimen is the fastest. When the intermediate region is a hard particle with a higher strength than the matrix, the crack propagates first along the particle-matrix interface and then in the vertical direction.To study the acoustic emission pattern during the damage of the model under static loading [36–40], the distributions of the number of AEs and the energy with the loading step were extracted for each of the three models under static loading, as shown in Figure 8.Figure 8
Number of unit acoustic emissions and energy characteristics of the numerical models under static load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)The numerical simulation results show that the number of AEs and the energy of all three models were concentrated between loading steps 20 and 60. This is mainly due to the presence of hard particles in the crack expansion path in the model, which prevents the crack from penetrating the stronger particles and causes damage along the contact surface with the matrix, resulting in the destruction of the units and the release of more energy. Model I produced the lowest number of AEs and the smallest amount of energy, with values of 2071 and 1.13 J, respectively. This was mainly due to the presence of hole defects in the model, which caused the cracks to be truncated, resulting in fewer units breaking and less energy being released.Figure9 shows the load-displacement curves of Models I, II, and III under static loading. The simulation results show that the load peaks at 0.44 mm, 0.52 mm, and 0.66 mm were 74831 N, 94437 N, and 102852 N, respectively. The ultimate load and displacement were the largest for Model III, while the load and displacement were the smallest for Model I, which contained the circular hole defect. There was a rapid decrease in the stress, which was mainly due to the presence of hole defects in the model, and when the crack extended to the bottom of the round hole, a crack quickly formed at the top of the defect, reducing the resistance of the model. The ultimate load was substituted into (5) to calculate the tensile strength of the model and compared with the results of (4), as shown in Table 2.Figure 9
Load-displacement curve of model under static load.Table 2
Calculation results of tensile strength by different calculation methods.
Calculation methodModels IModels IIModels IIIEquation 48.07MPa10.21MPa11.93MPaEquation 57.93MPa10.01MPa10.90MPaThe calculations show that it is reasonable to use the calculated tensile strength of inhomogeneous materials of different strengths, with a maximum error of 0.95%.


## 4.2. Dynamic Tensile Strength of Inhomogeneity Numerical Models

A stress wave, as shown in Figure5, was applied to the central axis of the upper boundary of the three models to simulate the impact load of a 4 kg hammer dropped from 2 m, with a peak stress wave of 43 kN and a duration of 1.2 ms. Figure 10 shows the damage process of the numerical model under the action of the falling weight impact. The numerical simulation results show that the three different structural models exhibit various damage modes under dynamic loading.Figure 10
Diagrams of damage process and minimum principal stress diagram under dynamic load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)Type I fracture cracks appeared in all three models at the lower boundary of the mid-axis of the precast crack, with Models I and II showing crack initiation at the precast crack at the bottom of the specimen att = 0.5 ms and Model III only showing stress concentration at the tip of the precast crack at t = 0.5 ms, with no type I fracture crack initiation. This is because the stress waves in the circular hole in Model III were reflected through the hole, which prevented some of the stress waves from propagating, resulting in the horizontal tensile displacement in the central axis of Model III lagging behind the other two models. As the stress wave continued to propagate, the type I crack produced in Model I continued to expand along the mid-axis region of the model specimen toward the top of the specimen, reaching the middle of the model at the upper boundary at t = 2 ms. After a type I crack formed in Model II, the crack began to expand rapidly along the mid-axis region of the model, and at t = 1 ms, the crack expanded to the lower endpoint of the circular particle, where the crack bifurcated, and the two bifurcated cracks 1 and 2 expanded along the boundary of the circular hard particle. At t = 1.5 ms, the two bifurcated cracks stopped after expanding to the endpoints on both sides. A third, vertically oriented type I crack emerged at the lower endpoint of the rounded particle. At this point, cracks 1, 2, and 3 all expanded vertically. At t = 2 ms, cracks 1 and 2 stopped expanding and did not reach the upper end of the model, while crack 3 stopped expanding at the center of the round particle.Att = 1 ms, in Model III, the type I crack expansion stopped at the lower endpoint of the circular hole due to the presence of the circular hole, which hindered the propagation of the stress wave. At the same time, two type I cracks formed on the left and right sides of the prefabricated cracks at the lower boundary of the model, mainly due to tensile waves reflected from the lower boundary of the circular hole. At t = 1.5 ms, type I cracks formed in the upper region of the circular hole defect and, at t = 2 ms, extended to the upper boundary of the model, while the two type I cracks sprouted at the lower boundary stopped extending after a certain distance. The numerical models of the three different structures showed significant differences in the damage process under the falling weight impact, with distinct crack expansion paths and damage patterns due to the various materials in the middle circular region.Figure11 shows the spatial distribution of AEs for the three different structural numerical models. Model III has the most complex form of damage, with concentrated unit tensile damage at the left and right ends of the circular hole and shear damage units in the upper region of the hole, mainly due to the presence of the circular hole, which causes the stress waves to reflect back and forth in this region.Figure 11
The distribution of AE under dynamic load.Figure12 shows the distribution of the number of AEs and the amount of energy with time for Models I to III under dynamic loading. The numerical simulation results show that the distribution of AEs over time differs for the three models. For example, in Model I, the number of AEs gradually increases at t = 0.4 ms, reaches a peak at t = 1 ms, and then gradually decreases. This is mainly due to the AE from the vertical tensile crack at the bottom of the model. At t = 1.5 ms, a second peak in the AE occurs at the top of the circular hole in the model due to crack initiation. Model II has only one peak in its AE distribution due to its homogeneous structure, with only one type I crack produced.Figure 12
Number of unit acoustic emissions and energy characteristics of the numerical models under dynamic load. (a) Model I; (b) model II; and (c) model III.
(a)(b)(c)In Model III, due to the presence of hard particles, the crack expands rapidly along the particle-matrix interface, causing the AE first peak to occur att = 0.8 ms; at t = 1.5 ms, as the peak stress wave gradually increases, the crack passes through the hard particle region, causing a second AE peak. The AE time distribution curve in Figure 12 shows that crack initiation causes an increase in the AE, while crack expansion and penetration cause a concentrated burst in the number of AEs. Furthermore, by observing the cumulative release energy pattern of the AEs, it is found that the released energy of the AEs increases with increasing loading step. The highest number of AEs and the highest cumulative release energy were 8068 and 6.58 J, respectively, for Model I. Model I had the most damaged units due to the presence of circular defects and multiple transverse reflections of stress waves inside the specimen, which resulted in more tensile cracks in the model.Compared with Model I, the number of AEs and the energy produced in Model III are 7898 and 4.62 J, respectively. Although the number of AEs is similar, there is a difference of 1.92 J in the energy released, which is mainly due to the presence of more shear damage in Model I, which causes shear damage at the edge of the defect, releasing more energy. It is worth noting that Model II shows an increase in the number of AEs and the amount of energy released at the point, where the crack reaches the lower endpoint of the round particle. In particular, the rapid increase in AE energy, when combined with the AE distribution diagram at this point, is found to be the main reason for the rapid increase in microcrack energy. The sheer damage is mainly due to the apparent stress concentration at the tip of the particle and crack, which causes shear damage to the unit in this region.To further investigate the force and deformation mechanisms of the different structural inhomogeneous models under dynamic loading, the load-displacement curves, inertial force curves, and vertical displacement curves with time for the three numerical models were calculated, as shown in Figure13.Figure 13
The results of the numerical model under dynamic load. (a) Load-displacement curve. (b) Model inertia force curve. (c) Vertical displacement curve.
(a)(b)(c)According to the load-displacement curve of the model in Figure13(a), the maximum ultimate load of the model with hard particles is 49.07 kN, the minimum ultimate load of the model with circular hole defects is 38.7 kN, and the maximum displacement of the ultimate load is 0.26 mm. It is worth mentioning that the presence of circular defects causes a significant load drop in the load-displacement curve, as shown in Figure 13(b). The three models also have different inertial effects under stress waves due to their structural inhomogeneities, which confirms the conclusions reached in Section 5 of this paper.Figure13(c) shows the vertical displacement over time for the three model monitoring points. The vertical displacement for the model with round particles has the fastest rate of change but the smallest value. In Model III, the rate of change of the vertical displacement is the smallest because it contains circular defects, but the ultimate displacement value is the largest. This is mainly due to the presence of circular defects, which obstruct the propagation of stress waves, causing significant hysteresis in crack initiation and expansion. The equivalent inertial forces calculated according to (8) are shown in Figure 13(b). The numerical simulation results show that the equivalent inertial force of Model III is the largest, while that of Model I is the smallest. Overall, there are no significant differences between Model II and Model III in any calculation aspect, but the two models do have considerable differences with Model I. The dynamic strength magnitudes derived from the three numerical models based on their rupture areas and (4) and (9) were calculated, as shown in Table 3.Table 3
Comparison of calculation methods for tensile strength of numerical models.
ModelModel IModel IIModel IIIInhomogeneity strength value13.12MPa22.35MPa15.20MPaInertia force Strength values29.09MPa30.24MPa25.53MPaTotal model Strength values42.21MPa52.59MPa40.73MPaTraditional calculation method44.26MPa54.65MPa41.22MPaError2.66％3.92％1.21％According to the calculation results in the table, the error between the method in this work and the traditional calculation method is within 4%, indicating that the method has good accuracy. It is worth noting that the calculated values of the method in this work are all less than the calculated values of the traditional method, which is mainly because the traditional calculation method only considers tensile damage, but there may be a small amount of shear damage during the dynamic loading process, resulting in an increase in the strength.


## 4.3. Mechanism of the Influence of the Radius of a Circular Hole on Strength

To investigate the effect of different defect radius sizes on rock specimen damage, three-point bent rock specimens with radiiR of 5 mm, 7.5 mm, 10 mm, 12.5 mm, and 15 mm were established. A shock wave, as shown in Figure 5, was applied to the central region of the upper boundary of the specimen. The final damage pattern of the rock specimen is shown in Figure 14.Figure 14
Failure pattern of different hole radius.As the radiusR of the circular hole increases, the pattern of the crack extension type of the rock specimen changes. When R = 5 mm, the crack extends from the top of the precast crack at the bottom of the specimen to the top of the specimen through the top and bottom of the circular hole. When R = 7.5 mm, two short tensile cracks form on both sides of the precast crack due to the stress waves reflecting at the hole boundary. As the hole diameter increases, the stress wave reflection effect increases, and the tensile cracks that form on either side of the precast crack expand toward the left and right sides of the hole. When R = 15 mm, the rock specimen shows a similar damage pattern to R = 12.5 mm, but the cracks start to sprout and expand as prefabricated cracks above the hole rather than at the hole. This indicates that when the diameter of the hole exceeds a certain value, the rock specimen no longer expands as a three-point bending type I crack.Table4 shows the contributions of inertia effects and crack extension to the dynamic tensile strength values for rock specimens with different hole radii, as calculated by the above equations. The calculation results show that the dynamic tensile strength of the rock gradually decreases as the hole radius increases, and the larger the radius, the more significant the decrease. At this point, although the strength of the specimen gradually decreases, the contribution of the crack extension area to the strength first decreases, then increases, and finally decreases. The influence of inertia effects during loading is always greater than the influence of inhomogeneities. In other words, the hole affects the integral tensile strength of the specimen mainly by reducing the inertia effects of the specimen.Table 4
Model strength of different hole radius.
Hole radius (mm)Inhomogeneous strength value (MPa)Inertial forces (MPa)Total strength (MPa)013.1229.0942.21512.9528.9241.877.515.2025.5340.731017.2521.2838.5312.517.7318.1235.851513.1215.3528.47


## 5. Conclusions

This paper investigated the tensile damage mechanism of rock materials under dynamic and static loads, used theoretical and numerical methods to analyze the influence of inertia and rock material inhomogeneity on the strength of rock materials under dynamic loads, and obtained the following conclusions:(1)
Under static loading, crack expansion follows the least energy dissipation principle, with cracks expanding along the weak side of the material or at the strength-weakness interface, and the strength of the material is proportional to the area of crack expansion.(2)
Under dynamic loading, crack expansion exhibits a multipath approach to expansion, following the shortest energy path and penetrating different regions of strength and weakness in a short period of time. Furthermore, inertial forces are generated as the rock transitions from a static deformation state to a dynamic deformation state in a short period of time. These two aspects increase the macroscopic strength of the rock material under dynamic loading.(3)
Based on the damage pattern of rock materials under dynamic and static loads, a strength calculation method based on material inhomogeneity and inertial forces is established, which, when compared with traditional calculation methods, has good calculation accuracy.(4)
Inertia effects dominate under dynamic loading. The larger the radii of hole defects in a rock specimen, the smaller the inertial force, and the smaller the dynamic tensile strength of the rock specimen.



---

*Source: 4559169-2022-03-30.xml*