SPACE WAVE EQUATION CONSIDERING DAMAGE-INDUCED WEAKENING AND STRAIN RATE DEPENDENCY OF ROCK

Original scientific paper Split Hopkinson pressure bar (SHPB) apparatus has been used to characterize the dynamic properties of rock. The effect of the geometry-induced stress wave dispersion, especially under a high strain rate due to impact loading, deserves attention. In this study, the space wave equation which can simulate the stress wave propagating through rock was established, with consideration of the damage-induced weakening and strain rate dependency. The specific parameters of the wave equation were determined using a hybrid genetic algorithm based on the SHPB test results. Furthermore, stress wave propagation in the SHPB was simulated with a finite difference method involving a developed constitutive model in order to verify the space wave equation. The results obtained from numerical simulation show good agreement with experimental test results. In summary, the proposed space wave equation can be effectively used to investigate the dynamic properties of rock.


Introduction
The diameter of a heterogeneous brittle rock specimen (i.e., granite, marble, concrete) for the highspeed impact loading test should be enlarged to minimize the material heterogeneity and to characterize the performance of the sample under the expected strain rate.A number of numerical and theoretical research papers have been made using the Split Hopkinson Pressure Bar test (SHPB) to characterize rock material [1 ÷ 7].SHPB apparatus with a typical pressure bar diameter of 100 mm is also used to characterize concrete materials under high strain rate [8 ÷ 10].The stress wave propagation in the apparatus was analysed with numerical simulation.From the dynamic failure process and dynamic properties obtained from the large diameter SHPB test, it is found that the geometry-induced dispersion effect becomes very important as the diameter of the SHPB enlarges due to radial inertia.When the appropriate waveform is applied, the stress wave dispersion effect can be reduced [11].
Therefore, a large diameter SHPB system is essential in order to characterize the dynamic properties of brittle rock materials.However, the pressure bar diameter increases and the impact energy magnifies the geometryinduced dispersion significantly.In such case, the reliability of the test results can be seriously harmed.Furthermore, development of a three dimensional wave equation is crucial to analyse the stress wave propagating through the large diameter brittle rock materials in numerical terms considering the influence of friction between tip and pressure bar on the SHPB experimental test as well as the confining pressure acting on the specimens [12,13].

Theoretical and experimental study 2.1 Establishment of wave equation
Generally, stress wave propagation problems in cylindrical specimens have been simplified into two dimensional problems, since it is difficult to solve three dimensional problems directly.In order to simplify three dimensional problems, rock materials are assumed to be an elastic medium which can consider the effect of the damage-induced weakening and the strain rate hardening as follows: where is the elastic modulus considering the influence of the strain rate and is the elastic modulus which can consider the influence of damageinduced weakening and strain rate hardening.Where α and m are distribution ratio.E 0 is the initial elastic modulus, k is the material parameter reflecting strain rate effect of rock materials under impact load, and σ s is the crack damage threshold.If we assume the problem as an axisymmetric condition considering the SHPB and specimen for the analytical derivation of the wave equation along the cylindrical elastic bar, the following boundary conditions can be applied to the wave equation (e.g., 0 , 0 , 0 ).Then, the two dimensional expressions of tensor Δ can be as follows: where r u is radial displacement, and z u is axial displacement.
According to the three dimensional wave equation in a cylindrical coordinate system [14], a two dimenisonal differential equation can be derived in an axisymmetric condition as follows: where λ is Lame's constant and μ is the rigidity modulus.If we take the axis of the specimen as the z-axis, σ zz and ε zz are the axial impact stress and strain applied to the axis of symmetry, and Young's modulus can be defined as Considering the correlation between damage and strain rate, it can be assumed that E = E 1 and E = E 2 .
According to the generalized Hooke's law, elastic modulus E and Poisson's ratio ν can be expressed with λ and μ as follows: ( ) Therefore, the elastic wave velocity of material can be expressed as the wave velocity considering damageinduced weakening and strain rate hardening as follows: , , where c 0 is material wave velocity; ρ is material density.
where c 1 is longitudinal wave velocity, and c 2 is shear wave velocity.
The above equations can be simplified as,

( )( )
Substituting Eq. ( 7) and Eq. ( 8) into Eq.( 3), the two dimensional axisymmetric wave equation of elastic material considering damage-induced weakening and strain rate dependency can be obtained as follows: ( ) This can be expressed in matrix form as follows:

Determination of characteristic parameters
In order to verify the applicability of the space wave equation, the characteristic parameters of the constitutive model should be determined.Based on the stress and strain data obtained from a large diameter SHPB test on granite, the corresponding characteristic parameters can be obtained through inversion analysis by compiling an adaptive hybrid genetic algorithm.In this paper, the experiment was carried out at the structure laboratory in the School of Architectural and Environmental Engineering, Nanyang Technological University (NTU) in Singapore.The testing installation is the SHPB, as shown in Fig. 1.
The dynamic strain measurement instrument is a dynamic strain indicator containing eight channels; the power system is a high-pressure compressed air device; and the speed of the striker bar is obtained by laser speedometer.These mating devices are shown in Fig. 2. The granite used in the test is mainly composed of pleocrystalline and hypidiomorphic granular granodiorite which is similar to a porphyritic texture.The density of granite is 2670 kg/m 3 and the maximum grain size is 3 ÷ 6 mm.To ensure its parallelism, roughness and smoothness, both ends of the cylinder specimen are smoothed by grinding and polishing.The diameter and thickness are 50 mm and 30 mm respectively.In order to reduce the boundary confinement induced by friction at the contact surfaces between the pressure bar and the specimen, a lubricant is generously applied to both ends of the specimen and the pressure bar.The sample is then located between the input bar and the transmission bar.The specimen, compressive bar and the axis of the bullets are in alignment along a straight line.The test conducted in this study uses a conventional SHPB apparatus for the typical heterogeneous brittle characteristics of rock materials.A circular brass strip (C11000) 18 mm in diameter and 1 mm thick was used, and using grease it was positioned in the center of the impact surface in the input bar, as shown in Fig. 3.The pressure bar and bullet material of the SHPB system are made of high strength stainless steel.The lengths of the input bar, transforming bar and absorption bar are all 2000 mm and the length of the bullet is 300 mm.The diameter of the aforementioned components is 50 mm.A schematic diagram of the SHPB apparatus and its monitoring system are presented in Fig. 4.
The waveforms of the incident wave, reflecting wave and transmission wave are presented in Fig. 5, and in which the pulse shaper was modified.From the waveforms, it can be found that the improved SHPB apparatus effectively modifies the square wave with high frequency oscillation into a half-sine wave and provides relatively smooth waveform and gentle rising edges.
The stress-strain curve of the granite specimen is presented in Fig. 6 depending on the three different strain rates.It shows that uniaxial dynamic compressive strength improves obviously with the increase in strain rate, so strain rate-dependency is noticeable.In additon, as the strain rate increases, the strain corresponding with the peek stress changes slightly, so it is weakly correlative with strain rate.On the other hand, from the figure, it can be found that granite elastic modulus is little changed, so it also can be considered a strain rate independent quantity.Technical Gazette 22, 4(2015), 1035-1042   The granite specimen's broken degree and broken form at different high strain rates are shown in Fig. 7. From the broken degree, it can be seen that, with the increase in strain rate, fragment size is considerably smaller and fragment numbers are larger, which shows a strong strain-dependency.From the broken form, it can be seen that, at different high strain rates, dynamic compression damage of the granite specimen exhibits an axial splitting failure mode, and the fracture surface is just like the fracture surface of a tensile fracture.a) Specimen (strain rate 50 s −1 ) b) Specimen (strain rate 75 s −1 ) c) Specimen (strain rate 97 s −1 ) Figure 7 Broken granite specimen at different high strain rate In this paper, a new hybrid algorithm is established by embedding the downhill simplex method of multivariate function into an adaptive genetic algorithm [14 ÷ 16].As an inverse analysis objective function, Eq. ( 1), which can take into account the damage-induced weakening factor and strain rate hardening factor, was adopted to characterize the dynamic properties of rock.The inversion analysis algorithm is compiled to establish the undetermined parameters in the model as shown in Tab. 1.

Derivation of boundary condition
Taking into account the radial and axial deformation of materials at the same time, the boundary conditions at the interface between the pressure bar and the specimen are quite complicated.Moreover, a large amount of computation time is necessary for the calculation of the compressive bar since the length of the pressure bar is nearly 80 times greater than the thickness of the specimen.Meanwhile, lateral and axial deformation of the pressure bar can be ignored since the pressure bar is high strength stainless steel that shows rigid behavior with very small deformation.With this in mind, in order to simplify the complicated boundary condition and improve computational efficiency, stress wave propagation in the pressure bar is not considered in the present study.Instead, waveforms measured at the incident bar end and transmission bar head are applied at both ends of the specimen directly in this numerical simulation.Hence, without consideration of the slender pressure bar, the space wave equation developed in this study can be implemented to simulate the stress wave propagation in the rock specimen.
If we equate radial step-size to axial step-size, z r ∆ = ∆ . For the boundary condition at both ends and the radial free surface of the specimen, we can set the following rules, respectively: (1) When r = 0, considering the characteristics of axisymmetry for the particle on the axis of the specimen, its radial displacement is 0 and it can be simplified as a one dimensional propagation.The wave equation can be simply expressed as: The formula can be changed into finite difference form as follows, is the differential form of axial displacement; i is the radial component; j is the axial component; k is the time component; f is the material surface ratio, and (2) When r = a, all stress components at the free surface of the bar are 0, where a is the radius of the bar.
( ) where rr σ is the radial stress; rr ε is the strain component due to radial loading; θ r σ is the angular stress component; and rz σ is the axial stress component.From above, Eq. (12a) can be expressed as: And the finite difference form of Eq. (12c) can be expressed as: (3) When z = 0, according to the continuity assumption at the interface between incident bar and specimen, (i.e., there is no relative displacement between the specimen and two compressive bars in the whole process of loading), the axial velocity of particles at the starting surface of the specimen equals the velocity of particles at the incident surface.
), ( Δ where Δt is the time component.Without consideration of the friction at the specimen ends, the radial contact surface is free from any constraints, namely σ r = 0.And in this situation, the finite difference form can be expressed as: (4) When z = l, according to the continuity assumption similarly at the interface between transmission bar and specimen, the particle velocity of the specimen equals the velocity of particles at the starting surface of the transmission bar.This is expressed as follows: ). ( Δ Similarly, if friction at the contact surface is not taken into account, radial stress can be expressed as σ r = 0, and the finite difference form is as follows: (16) (5) When t = k, if there is no relative displacement between the specimen and two compressive bars in the whole process of loading, the particle displacement at the initial side and the other side of the specimen should be equal to the corresponding location of the bar, and they can be expressed as follows: .
Boundary conditions can be determined based on the pre-defined rules for the numerical calculation program under the impact loading.
(6) The axial dynamic stress in the specimen can be expressed as follows: And the finite difference form can be expressed as:

Numerical analysis and verification
For numerical analysis, the finite difference (FD) method is adopted in this study [17,18].Differential equations of the wave equation and the boundary conditions derived in previous sections are fomulated using a central finite differentiation scheme.The developed FD code implementing the space wave equation is compiled with VC++.In the paper, the radial grid and axial grid are all set as 1mm.The specimen is round with the diameter 50 mm, and the thickness is 30 mm.Since the calculation of the stress wave propagating in the compressive bar is not considered in this program, enough computation time is ensured for the stress to reach equilibrium state in the specimen.Here, 180 µs is taken as operation time range.The selected parameters for numerical analysis are tabulated in Tab. 2.
The axial strain at the middle cross-section of the granite obtained from numerical simulation is compared to the axial strain obtained from the experimental test as presented in Fig. 8.The axial strain curve obtained at the surface of the specimen creates a lower boundary and the axial strain curve obtained at the center of the specimen forms an upper boundary.The axial strain curve obtained from the experimental test is located just between the two curves obtained from numerical study.Therefore, the calculation results are reasonable.Meanwhile, it can be found that the particle's axial strain is not exactly the same.This is because of the geometric dispersion effect.Radial deformation at the specimen surface transforms the flat cross-section of the specimen into a curved surface.Using FD-based numerical simulation, radial tensile strain can be calculated at the cross section of the specimen before failure as shown in Fig. 9. Through the figure, it can be found that the maximum tensile strain before failure is about 0,01 %.This result is essentially identical to the critical tensile strain of granite materials obtained from a dynamic test in Brazil under the same conditions where relative test conditions are the conditions of same material, same size, same experiment equipment and same load, and the detailed information has been described in [19].It shows that significant tensile failure of the specimen under a shock pressure load can be caused due to the rapid radial deformation that accompanies tensile stress exceeding the tensile strength of the specimen.This failure mechanism is consistent with the findings of Yue Zhai et al. [19].The FD code developed in this study, which can calculate analytical space wave equation, can be effectively used to estimate the dynamic properties of rock materials by simulating the SHPB test under various multidirectional states of stress.

Conclusions
In the analytical study, a two dimensional axisymmetric wave equation, which can simulate a stress wave propagating through the rock, was developed which takes into account damage-induced weakening and strain rate dependency.
In the experimental study, high strain rate tests were conducted on granite specimens with SHPB apparatus.The characteristic parameters of the wave equation were determined using a hybrid genetic algorithm based on the SHPB test results.
In the numerical study, stress wave propagation in the SHPB was simulated using the finite difference method with a central finite differentiation scheme.The established stress wave equation is implemented in the FD code.The axial strain curve obtained from the numerical simulation shows good agreement with the experimental test.The FD code developed in this study can be effectively used to estimate the dynamic properties of rock materials by simulating the SHPB test under various multidirectional states of stress in the future.

Figure 3
Figure 3 The position of pulse shaper

Figure 5 Figure 6
Figure 5 Pulse curve with pulse shaper

Figure 8 Figure 9
Figure 8Comparison between the middle cross-section axial strain of granite specimen and the testing strain (50 s −1 )

Table 1
Results of inverse analysis of elastic dynamic constitutive modelResults of inverse analysis of elastic dynamic constitutive model

Table 2
Selected calculation parameters for numerical simulation