- Email: [email protected]
Numerical simulation of Brazilian disk rock failure under static and dynamic loading
ARTICLE IN PRESS
International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252 www.elsevier.com/locate/ijrmms
Numerical simulation of Brazilian disk rock failure under static and dynamic loading W.C. Zhu, C.A. Tang Center for Rock Instability and Seismicity Research, Box 138, Northeastern University, Shenyang 110004, PR China Accepted 27 June 2005 Available online 22 August 2005
Abstract A numerical simulator based on RFPA (Rock Failure Process Analysis) is used to study the deformation and failure process of a Brazilian disk of heterogeneous rock when subjected to static and dynamic loading conditions. In this simulator, the heterogeneity of rock is considered by assuming that the material properties of elements conform to a Weibull distribution; an elastic damagebased law that considers the strain-rate dependency is used to describe the constitutive law at mesoscopic scale; and a finite element program is employed as a basic stress analysis tool. The simulator is firstly validated by simulating the dynamic spalling of a homogeneous rock bar and by comparing with the theoretical and experimental results. Then, the failure process of a Brazilian disk of rock subjected to static and dynamic loading is numerically simulated, and the numerical results are compared with the available experimental results. Particular attention is given to the typical failure patterns of the rock disk when the incident compressive stress waves with different amplitudes are applied. The numerical simulation also identifies the failure mechanisms of rock during dynamic failure processes that are closely related to the propagation of the stress wave. r 2005 Elsevier Ltd. All rights reserved. Keywords: Rock failure; Brazilian disk; Dynamic loading; Numerical simulation
1. Introduction Better understanding of the dynamic fracture processes of rocks promises benefit in many areas from rock mechanics to mining engineering and earthquake prediction. It is essential to understand how fractures initiate and propagate under different loading conditions in order to provide understanding of rock fracture process that occur in the engineering fields. During the past several decades, researchers have carried out many dynamic tests of rock using the apparatus such as the Split Hopkinson pressure bar (SHPB). From these experimental investigations, it has been understood that the mechanical properties (or responses) and failure characteristics of rock are sensitive to the strain rate. It is also a well-established fact that rocks exhibit increases Corresponding author. Tel.: 86 24 83687705; fax: 86 24 83671626.
E-mail address: [email protected] (W.C. Zhu). 1365-1609/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2005.06.008
in strength as the loading rate or strain rate increases [1–11]. Based on these experimental results, failure criteria considering the strain rate effect of rocks are proposed [12,13]. However, not so well understood is the dynamic failure process and failure mechanism, which are closely related to the microstructural behavior in terms of activation, growth, and coalescence of cracks when the rock is subjected to different loading conditions. With regard to the SHPB experiments, they are widely used to determine strain-rate dependency of the dynamic strength of rock, but not to study the failure process and failure mechanisms that are closely related to the stress wave propagation in heterogeneous rock. Continuum damage mechanics is considered to be appropriate to describe the failure process resulting from dynamic loading. Based on this approach, many theoretical models were developed to study the dynamic damage evolution of brittle materials such as rock and
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Nomenclature
A b C D E, E0 f c0 , f t0 f cr , f tr f cs0 k K m M p
parameter to reflect the stress rate effect on strength width of load-bearing strip damping matrix diameter of rock disk damaged and undamaged (initial) elastic modulus of element dynamic compressive and tensile strengths, respectively, of element residual dynamic compressive and tensile strengths, respectively, of element static uniaxial compressive strength ratio of tensile and compressive strength stiffness matrix homogeneity index, shape parameter of Weibull distribution mass matrix applied incident stress
concrete with microstructures [12,14–20]. Most of these models are developed based on combining the theory of damage mechanics or fracture mechanics with some statistic treatment to account for the random distribution of microcracks. Because many micromechanical parameters such as crack density must be used to describe the crack distribution in these damage models, it is often difficult to implement such models in a numerical code [21]. Therefore, there is a real need for further developments in the material modeling of rock to reflect more closely the physical mechanisms responsible for failure. In the numerical simulation with the discrete element method, the basic fracture patterns of rock and concrete subject to dynamic tension are captured [22]. But in the discrete element model, there is no accepted general method of averaging the contact displacement between particles in order to obtain the value of resulting macroscopic strain or stress [22]. The explicit non-linear finite element code DYNA 2D/3D is used to analyze a variety of dynamic problems, which is generally concentrated on collision analysis, such as an elastic sphere on an elastic–plastic metal material and sheet metal forming. Although DYNA 2D had been also used to calculate the response of jointed rock subjected to blasting [14,23], it was not focused on capturing the dynamic failure mechanism of rock. In fact, the stress distribution and damage initiation in the rock is closely related to stress wave propagation; moreover, the local damage in rock can in turn affect the
237
pmax R u
amplitude of stress wave vector of externally applied loads parameter of elements that conforms to Weibull distribution, such as elastic modulus and strength _ Udisplacement, € U, U, velocity, and acceleration vectors of elements u0 scale parameter of Weibull distribution e strain e0 strain at peak stress e1, e2, e3 principal strain ec0 strain at the peak compressive stress et0 strain at the peak tensile stress etu ultimate tensile strain ¯ equivalent strain f internal friction angle Z ultimate tensile strain coefficient ðZ ¼ tu =t0 Þ, l residual strength coefficient ðl ¼ f tr =f t0 Þ m Poisson’s ratio s stress s1 , s2 , s3 principal stress o damage variable
propagation of the stress wave. Numerical simulation of rock subjected to dynamic loading conditions requires the implementation of a constitutive relation capable of capturing the major features of the failure process during the propagation of stress waves. The Rock Failure Process Analysis (RFPA) Code, which had been developed and introduced by Tang [24,25] in many documents, has been successfully used to simulate the fracture process and failure mechanism of rock under static loading. In RFPA, in order to numerically denote the heterogeneity of rock, the rock sample is discretized into many mesoscopic elements and their material properties are conformed to the Weibull distribution [26,27]. The constitutive law of these mesoscopic elements is formed based on elastic damage mechanics. In cooperation with Professor Chau [28] at the Department of Civil and Structural Engineering, The Hong Kong Polytechnic University, the author proposed a stress-rate dependent constitutive law and implemented it into the original RFPA code [24,25], and therefore the RFPA is extended for analyzing the failure process of rock under dynamic loading [28,29]. From January 2003 on, based on that, Zhu et al. [29] have made many further improvements in utilizing the faster finite element solver and introducing more boundary conditions, and all these have been also implemented into the original RFPA code [24,25]. Therefore, this upgraded simulator used in this study is also called RFPA. The contribution of this work is to validate the simulator by reproducing the dynamic spalling of
ARTICLE IN PRESS 238
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
homogeneous material bar, and then to investigate the failure processes and related failure mechanism of a Brazilian disk of rock under static and dynamic loading. Based on the numerical simulations, the different failure processes are captured, and the related failure mechanisms of rock under static and dynamic loading are identified.
2. Numerical simulator description Briefly, the simulator used in this work is a twodimensional (2D) code that can simulate the fracture and failure process of rock under static or dynamic loading conditions. To model the failure of rock material (or rock mass), the rock medium is assumed to be composed of many mesoscopic elements whose material properties are different from one to another and are specified according to a Weibull distribution [30]. The finite element method is employed to obtain the stress fields in the mesoscopic elements. Elastic damage mechanics is used to describe the constitutive law of the meso-scale elements when the maximum tensile strain criterion and the Mohr–Coulomb criterion that incorporate the effect of stress rate are utilized as damage thresholds [27].
(a)
m = 1.5 (The values range from 1.0 to 363.)
(b)
m = 5.0 (The values range from 25.1 to 147)
(c)
m = 10.0 (The values range from 50.1 to 121)
2.1. Assignment of material properties In RFPA, the solid or structure is assumed to be composed of many mesoscopic elements with the same size, and the mechanical properties of these elements are assumed to conform to a given Weibull distribution as defined by the following probability density function (pdf): m m u m1 u f ðuÞ ¼ exp , (1) u0 u0 u0 where u is the mechanical parameter of the element (such as strength or elastic modulus); the scale parameter u0 is related to the average of the element parameters and the parameter m defines the shape of the distribution function. From the properties of the Weibull distribution, a larger value of m implies a more heterogeneous material and vice versa. Therefore, the parameter m is called the homogeneity index in our numerical simulations. In the definition of the Weibull distribution, the value of the parameter m must be larger than 1.0. Using the pdf, in a computer simulation of a medium composed of many mesoscopic elements, one can produce numerically a heterogeneous material. The computationally produced heterogeneous medium is analogous to a real specimen tested in the laboratory, so in this investigation it is referred to as a numerical specimen. Fig. 1 presents some 2D numerical specimens, which are all composed of 50 50 elements, produced randomly by computer according to the Weibull distributions with
Fig. 1. Numerical specimens with different homogeneity index m. The specimen is composed of 2500 square elements and the Weibull parameter u0 is 100 MPa. The intensity of the gray scale in the figure denotes the magnitude of strength of elements.
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
different homogeneity indices. In Fig. 1, the different degree of gray color corresponds to a different magnitude of strength of the elements. In general, it has been assumed that the Young’s modulus and the strength of the mesoscopic elements that were used to simulate a rock specimen conform to two individual distributions with the same heterogeneity indices. The mesoscopic elements themselves are isotropic and homogeneous. The mesoscopic elements in the specimen must be sufficiently small to reflect the heterogeneous mechanical properties of materials at the meso-scale and still provide conditions under which the computer can perform the analysis efficiently. 2.2. The constitutive model for the mesoscopic element Initially an element is considered elastic, with elastic properties defined by Young’s modulus and Poisson’s ratio. The stress–strain relation for an element is considered linear elastic until the given damage threshold is attained, and then is modified by softening. Under a dynamic stress state, the elements undergo damage when one of the following damage criteria is satisfied at the element level: 1 ¼ kf c0 =E 0 ; s1
1 þ sin f s3 Xf c0 , 1 sin f
(2)
where f c0 is the dynamic uniaxial compressive strength of the element, which is closely related to the strain rate (or stress rate) of dynamic loading condition. E0 is initial elastic modulus of the element that is assumed to be not affected by the stress rate, k is the ratio of compressive and tensile strength and f is the internal frictional angle of the element. s1 and s3 are the major and minor principal stresses of the element. In this study, the following relation between dynamic uniaxial compressive strength and loading rate, which has been proposed by Zhao [13], is used to reflect the effect of stress rate on the dynamic strength: _ f_cs0 Þ þ f cs0 ; f c0 ¼ A logðs=
_ f_cs0 , s4
239
be damaged in either tension (corresponding to the maximum tensile strain criterion) or shear (corresponding to the Mohr–Coulomb criterion). Once Eq. (2) is satisfied at the element level, the elastic modulus of the element is reduced according to the following expression: E ¼ ð1 oÞE 0 ,
(4)
where o represents a damage variable, E and E0 are the elastic moduli of the damaged and the undamaged element, respectively. In the current method, the element as well as its damage is assumed isotropic, so the E, E0 and o are all scalar quantities. The sign convention used throughout this paper is that compressive stress and strain are positive. When the mesoscopic element is in a uniaxial stress state (both uniaxial compression and uniaxial tension), the constitutive relation of elements is as illustrated in Fig. 2. Initially, the stress–strain curve is linear elastic and no damage exists, i.e. o ¼ 0. When the maximum tensile strain criterion is satisfied, damage occurs in the element in a brittle mode. The constitutive relation of the mesoscopic element under uniaxial tension, as shown in the third quartile of Fig. 2, can be expressed as 8 4t0 ; > < 0; lt0 o ¼ 1 ; tu opt0 ; (5) > : 1; ptu ; where l is the residual strength coefficient, which is given as f tr ¼ lf t0 and f t0 is the tensile strength of element. The parameter t0 is the strain at the elastic limit, which is called the threshold strain and σ
fc0
(3)
where f c0 is also the dynamic uniaxial compressive strength (MPa), s_ is the stress rate (MPa/s), f cs0 is the uniaxial compressive strength at the quasi-static stress rate f_cs0 that is approximately 5 102 MPa/s, and A is a parameter depending on the material. In addition, the experimental results of Zhao [13] also indicated that the ratio of tensile and compressive strength ðk ¼ f c0 =f t0 Þ and internal frictional angle f are not influenced by the stress rate. In this respect, when the f c0 obtained from Eq. (3) is substituted into Eq. (2), the effect of strain rate on the strength of elements will be incorporated. The first part of Eq. (2) is the maximum tensile strain criterion, whilst the second part is the classical Mohr–Coulomb failure criterion for tensile and shear damage thresholds, respectively. Thus, an element may
fcr ε tu
ε t0
- ftr
ε c0
ε
- ft0 Fig. 2. Elastic damage constitutive law of element under uniaxial stress state (Here ft0 and ftr are dynamic uniaxial tensile strength and residual uniaxial tensile strength of element, respectively. And fc0 and fcr are dynamic uniaxial compressive strength and residual corresponding strength of element, respectively.)
ARTICLE IN PRESS 240
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
calculated as t0 ¼ f t0 =E 0 .
(6)
From Eq. (6), it can be found that t0 also depends on the stress rate of element because f c0 is related to the stress rate according to Eq. (3). tu is the ultimate tensile strain of the element, describing the state at which the element would be completely damaged. The ultimate tensile strain is defined as tu ¼ Zt0 , where Z is called the ultimate strain coefficient. Both the residual strength coefficient (l) and ultimate strain coefficient (Z) are assumed to be independent of the stress rate. Additionally, it is assumed that the damage of a mesoscopic element under multiaxial states of stress is also isotropic and elastic. The constitutive equation described above can be extended for application to 3D stress states when the tensile strain threshold is attained. In multiaxial stress states, the element is still subject to damage in the tensile mode when the equivalent major tensile strain ¯ attains the above threshold strain t0 . The equivalent principal strain ¯ is defined by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (7) ¯ ¼ h1 i2 þ h2 i2 þ h3 i2 , where 1 , 2 and 3 are three principal strains, and o4is a function defined as follows: ( x; xX0; ox4 ¼ (8) 0; xo0: For an element subject to a multiaxial state of stress, the constitutive law can be easily obtained by substituting the strain e in Eq. (5) with the equivalent strain ¯ defined by Eqs. (7) and (8). The damage variable is expressed as 8 ¯4t0 ; > < 1; lt0 (9) o ¼ 1 ¯ ; tu o¯pt0 ; > : 1; ¯ptu : In the same way as for uniaxial tension, when the element is under uniaxial compression but damaged according to the Mohr–Coulomb criterion, the expression for damage variable o can be described as ( 0; oc0 ; o¼ (10) 1 lc0 ; Xc0 ; where l is also the residual strength coefficient, and f cr =f c0 ¼ f tr =f t0 ¼ l is assumed to apply when the element is under either uniaxial compression or tension. Also, c0 is calculated as c0 ¼ f c0 =E 0 .
(11)
If an element is under a multiaxial stress state and its strength satisfies the Mohr–Coulomb criterion, damage occurs, and it is necessary to consider the effect of the other principal stress in the model during damage
evolution. When the Mohr–Coulomb criterion is satisfied, the maximum principal (compressive) strain c0 can be calculated at the peak value of the maximum principal (compressive) stress: 1 1 þ sin f s3 mðs1 þ s2 Þ . f c0 þ (12) c0 ¼ E0 1 sin f It is assumed that shear damage evolution is related only to the maximum compressive principal strain 1 . So the maximum compressive principal strain 1 of the damaged element is substituted for the uniaxial compressive strain e in Eq. (10), extending the equation to represent shear damage under triaxial stress states: ( 0; 1 oc0 ; o¼ (13) lc0 1 1 ; 1 Xc0 : Based on the above equations, it can be found that the damage variable o is related to the stress rate because it is calculated based on the values of c0 or t0 that are obtained directly according to the stress-rate dependent strength of the element ðf c0 or f t0 Þ according to Eq. (6) or (11). Of course, the above-mentioned constitutive law can be used to simulate the failure process of rock under static loading when the stress rate effect on the strength as expressed by Eq. (3) is not incorporated. Therefore, the stress-rate dependent constitutive law described above is based on the original one for static analysis and is also compatible with it. 2.3. Finite element implementation As described above, the rock specimen is composed of many square elements with the same size. These elements are also acting as the four-noded iso-parametric elements for finite element analysis. The equilibrium equations governing the linear dynamic response of a system of finite elements can be expressed in the following form: € þ CU _ þ KU ¼ R, MU
(14)
where M, C and K are the mass, damping, and stiffness matrices; R is the vector of externally applied loads; and _ and U € are the displacement, velocity, and U, U acceleration vectors of the finite elements. A direct step-by-step integration procedure is found suitable for solving the problem in which a body is subjected to a short duration impulse loading [31]. At initial condition, the elements are elastic, and their stresses can be calculated via step-by-step integration. At a time step, the principal stresses from the previous time step are subtracted from the current principal stresses and divided by the time step in order to calculate the stress rates of the elements. Similarly, when the current minor principal stress and corresponding stress
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
rate are both negative, the increased strength of an element due to the increase of absolute value of this stress rate can also be obtained. When the effect of stress rate on the dynamic compressive strength is considered according to Eq. (3) and the stresses of (or strains) of elements meet the maximum tensile strain criterion or Mohr–Coulomb criterion as given in Eq. (2), the elements are damaged according to the constitutive law given before, then their stress will be re-analyzed iteratively for the current boundary conditions in order to reflect the stress redistribution at this time step. The program does not consider the analysis of the next time step until no new damaged elements are found for the iterative step of this time step. The above model described in the two preceding sections is implemented into the RFPA [24,25], and therefore RFPA is extended for analysis of rock under dynamic loading.
241
which means that no strain rate effect is considered in this numerical simulation. The bar is 200 mm in length, 5 mm in width, and it is simulated with 200 5 elements in 2D plane stress problem. The incident compressive stress wave is applied at the top end of the bar and the bottom end of it is kept free. Two compressive incident stress waves as shown in Fig. 3(b) are applied on the top surface of bar, in order to illustrate the effect of stress amplitudes on spalling. In Fig. 4, the stress distributions during the fracture process of the bar subjected to two different stress waves are given, where the brightness of the gray scale indicates the relative magnitude of the maximum shear stress in the mesoscopic elements. From this Figure, the propagation of the stress wave can be seen clearly. For incident compressive stress wave I (pmax ¼ 20 MPa), as shown in Fig. 4(a), spalling occurs, induced when the incident compressive stress wave travels down and is reflected from the bottom free end of the bar, and
3. Simulator validation As discussed by many researchers [32–35], the dynamic spalling of quasi-brittle material such as rock occurs when an incident compressive wave is reflected into a tensile one on reaching a free end. If the incident compressive stress is lower than the compressive strength but its absolute value exceeds the material dynamic tensile strength, the specimen will spall during the reflection process. In this section, the dynamic spalling of a homogeneous quasi-brittle material bar is studied in order to validate RFPA in simulating the dynamic failure process. The geometries and loading conditions for the bar specimen are shown in Fig. 3(a). The mechanical properties of this bar are homogeneous and are specified as: the Young’s modulus is 62.0 GPa, uniaxial compressive strength is 200 MPa, direct tensile strength is 19.0 MPa, and Poisson’s ratio is 0.20. The parameter A for this homogeneous bar is specified to be zero,
(10µs)
(20µs)
(a)
Stress wave I (Peak value of incident stress is 50 MPa)
(30µs)
(40µs)
(50µs)
(60µs)
(70µs)
(80µs)
(90µs)
(100µs)
p
100
10 mm
p (MPa)
200 mm
75
II
50 I 25 (10µs)
(a)
(b)
5
10
(b)
(20µs)
(30µs)
(40µs) (42.5µs)
(43µs)
(45µs)
(50µs)
(60µs)
(70µs)
Stress wave II (Peak value of incident stress is 100 MPa)
t(µs)
Fig. 3. (a) Geometries and loading conditions for homogeneous bar subjected to dynamic loading. (b) The incident compressive stress waves with different peak amplitudes applied on the top surface of bar.
Fig. 4. Dynamic spalling of homogeneous bar is numerically simulated. In (a) the amplitude of incident compressive stress wave is 20 MPa; and in (b) the amplitude of incident compressive stress wave is 40 MPa.
ARTICLE IN PRESS 242
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
reaches the tensile strength of the material. The fracture face is at 22 mm above the bottom end of the bar. Finally, the stress waves in two parts of the bar only propagate individually with a certain velocity and no further spalling is induced. Two spallings occur for incident compressive stress wave II (pmax ¼ 40 MPa) when the level of the stress wave becomes more than twice the tensile strength of the material. As shown in Fig. 4(b), the first spall fracture, upon its abrupt creation at 12 mm above the bottom free end of bar, places a free surface in front of the trailing portion of the original wave. When the material fails, the fraction of the tensile wave that has passed the
failure point, continues traveling along the specimen, and it causes another spalling since the magnitude of tensile stress is also larger than the tensile strength of the bar. The second fracture occurs 25 mm away from the original bottom end of the bar. No more spallings occur because the intensity of the rear portion of the wave is no longer greater than the tensile strength of the material. This numerical result compares well with the theoretical results obtained by Rinehart [32] based on 1D wave propagation theory and the experimental observation given by Brara et al. [22]. From the above numerical simulation, it can be concluded that RFPA is effective in capturing the dynamic failure process of quasi-brittle material subjected to dynamic loading in these circumstances.
y
p
4. Numerical simulations of static and dynamic failure 4.1. Numerical rock disk rock Steel plate
D =150 mm
Rock disk x
25
2
b
Steel plate
Fig. 5. Geometries and loading conditions for disk specimen subjected to static and dynamic loading (D ¼ 150 mm, b=D ¼ 0:12).
The geometries and loading conditions for disk specimens are shown in Fig. 5. In our numerical simulation, the rock has a Young’s modulus of 62.2 GPa, a uniaxial compressive strength of 150.0 MPa, a direct tensile strength of 15.0 MPa, a mass density of 2500 kg/m3, and a Poisson’s ratio of 0.28. The parameter A of the material is specified to be 13.0, which is given according to the uniaxial compressive strength of this rock. When the Weibull distribution parameters as listed in Table 1, are applied to a standard numerical specimen that is 200 100 mm in dimensions and composed of 200 100 elements, the standard numerical specimen will have the above-mentioned mechanical properties of this rock when subjected to static loading. The rock disk is 150 mm in diameter, and in all cases, it is simulated with 150 150 elements as a 2D plane stress problem. Since we assume that the cushion materials (load-bearing strip) and the steel plate will not fail in numerical simulations, they are specified to be homogeneous and have a very high strength. Here the ratio between the width of load-bearing strip and diameter of disk (b/D) is selected to be of 0.12 (i.e. b=D ¼ 0:12) in order to allow the specimen to statically fail in tension along the loading line and not locally by compression [36].
Table 1 Weibull distribution parameters used for RFPA to simulate rock Materials
Homogeneity index m
Mean of Young’s modulus (GPa)
Mean of uniaxial compressive strength (MPa)
Poisson’s ratio
Ratio of compressive and tensile strength
Rock Steel plate Cardboard
5 Homogeneous Homogeneous
70 210 10
510 800 500
0.28 0.3 0.3
12 1 1
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
4 3.5
113
3 Load (kN)
For the static loading condition, displacement is imposed on the top loading paten step by step with the bottom platen fixed in the vertical direction, which simulates the compressive loading with displacement control. Similarly, for the dynamic loading conditions, the incident compressive stress is imposed at the top loading platen while the bottom platen is also fixed in the vertical direction.
2.5
77
2 1.5 1
180 23
0.5
4.2. Under static loading
243
114
0
Initially, the rock specimen is elastic and the stress distribution on the vertical diameter is given in Fig. 6, and it is compared to theoretical results for a specific width of load-bearing strips [37]. In the figure, the stress is normalized with 2P=ðpDtÞ (here P is the applied load and t the unit thickness for plane stress analysis). It can be seen that the tensile stress distribution at the center area of the rock specimen shows good agreement with the theoretical results. Because the Young’s moduli of elements in the numerical specimen are assumed to be heterogeneous, the stress distribution across the loaded diameter is also heterogeneous and fluctuates around the theoretical curves. Fig. 7 shows the load response of the rock specimen during the static failure process. The load-step curve is basically linear until peak loading is attained. Then the specimen loses its loading capacity substantially because this specimen is relatively homogeneous. The failure patterns at several typical loading steps (as indicated on the curve of Fig. 7) are shown in Fig. 8. In this figure, the distribution of damaged elements, distribution of maximum shear stress and cracking pattern in the rock specimen subjected to static loading are presented. In Fig. 8(a), the damaged elements are denoted with different color, i.e. white and dark gray for elements damaged in shear and tensile at current step, respectively, and black for damaged elements at all the
1 Theoretical (Hondros) σx 0. 5
-10
0
Numerical σx 0
10
20
σy σy
0
20
40
60
80
100 120 140 160 180 200 Step
Fig. 7. Load-step curves of specimens subjected to static loading (numerical results, the applied displacement at the top loading platen is 0.002 mm/step).
previous steps. In Fig. 8(b) and (c), the brightness indicates the relative magnitude of maximum shear stress and elastic modulus of elements, respectively. Since the damage of elements leads to degradation of its elastic modulus, the cracking patterns can be clearly given from the elastic modulus figure of elements. The first element is damaged in tensile mode near the vertical diameter of the specimen (as shown at Step 23 of Fig. 8). With the increase of the external loading, a damage zone near the vertical diameter is developed due to the tensile damage of many elements (as Steps 77 and 113 of Fig. 8). Just after the peak load the damaged elements in the damage zone coalescence with each other to form the primary crack along the vertical diameter of specimen. At the same loading step, the secondary diagonal cracks near two loading platens also develop. It should be mentioned that, for this specimen, fracture initiation does not mean collapse of the specimen. From Fig. 8, it can be seen that the initiation of fractures (at Step 23 as shown in Fig. 7) usually occurs when the load is 20% of the peak load (at Step 113). Also, a damage zone that is distributed in a wide range around the vertical diameter of disk is developed, which agrees well with the acoustic emission observation of Yanagidani et al. [38]. Besides, due to the heterogeneity, the path of the primary crack in the rock disk, which is not a straight line, is just like the common phenomena that are observed in many laboratory experiments [36–39].
30
4.3. Under dynamic loading -0.5
-1 Fig. 6. Theoretical and numerical results of stresses distribution across loaded diameter for disk specimen under indirect test.
In an attempt to more fully comprehend the response of the rock disk to the different amplitudes of incident compressive stress waves, the three incident compressive stress waves as shown in Fig. 9 are considered in the numerical simulations. Three incident compressive stress waves with different amplitudes (pmax) are applied on
ARTICLE IN PRESS 244
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 8. Failure process of rock disk subjected to static loading (The applied displacement at the top loading platen is 0.002 mm/step. Damage initiates at step 23, and load attains its peak value at step 113.)
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
III
p (MPa)
100 80
II
60
I
40 20 10
20
30
40 t (µs)
50
60
Fig. 9. Incident compressive stress waves imposed on the top loading platen above rock specimen: (I) pmax ¼ 50 MPa; (II) pmax ¼ 70 MPa; (III) pmax ¼ 100 MPa.
the top loading plates (as shown in Fig. 5). The amplitudes of the three incident compressive stress waves are 50, 70 and 100 MPa, respectively. 4.3.1. Elastic stress wave Firstly, the elastic stress field for the homogeneous specimen under the stress wave I is given, in order to illustrate the stress wave propagation in the homogeneous rock disk. The material parameters are specified for the material properties of this rock, i.e. the Young’s modulus is 62.2 GPa, Poisson’s ratio is 0.28 and mass density is 2500 kg/m3. The photoelastic fringes of maximum shear stress at different times are obtained from numerical simulation and are shown in Fig. 10. It takes t ¼ 3 ms for the stress wave to arrive the top end of the rock disk and it is partially reflected back into the top steel plate and partially transmitted into the rock specimen. Then, the stress wave propagates downwards and arrive at the bottom end of the rock disk around t ¼ 35 ms, at which time the stress wave in the rock specimens will also be partially reflected back into the rock specimen and partially transmitted into the bottom steel plate. About 3 ms later, the stress wave that is transmitted into the bottom steel plate reaches the bottom fixed end of the steel plate and is mostly reflected because the bottom surface of the loading plate is fixed in the vertical direction. Under the superposition of incident and reflected stress waves in the rock disk, the stress distribution in the homogeneous disk appear identical to a statically loaded specimen at t ¼ 85 ms, which means that the loads at both contacts points of the specimen become equal and the specimen is in equilibrium at this time. However, the stress deviates from this transient equilibrium condition as time passes because the stress in the lower half part of the disk become relatively high due to the superposition of the incident stress wave and the reflected one after the reflection at the bottom fixed end. However, another equilibrium condition is also achieved at about
245
t ¼ 179 ms. Therefore, it can be predicted that a periodical stress distribution in the rock disk is achieved under this dynamic loading condition. The ‘photoelastic fringe’ of stress distribution obtained from the numerical simulation indicates that the stress distribution in a rock disk is not so regularly distributed as that under the static loading, but changes periodically with time. Before the equilibrium condition at t ¼ 85 ms, the top end of the rock disk has the larger compressive stress, while the bottom end of the rock disk is unstressed. In contrast, after t ¼ 85 ms, the relative high stress is induced at the bottom end of the rock disk. This means that the locations with the maximum stresses are shifted to the bottom end of the rock disk. Therefore, the stress distribution in the rock disk changes as the time passes, which makes the deformation and failure process of rock under dynamic loading much complicated. It should also be mentioned that the rock disk under this dynamic loading conditions is quite different from that in the SHPB experiments because the rock specimen in the SHPB experiments remains in equilibrium until the time of failure. Hence, it can be seen that the numerical simulations of this work are most useful for studying the complex failure process and failure mechanisms of rock under a variety of dynamic loading conditions. 4.3.2. Failure process and failure mechanisms Since no damage of elements occurs when stress wave I is applied, the failure process analysis is only concentrated on the other two incident compressive stress waves. The maximum shear stress distribution, as well as cracking sequence from beginning to final failure, for the rock disk under stress wave II is shown in Fig. 11. Profiles for the normal stresses along the vertical diameter at selected times are also presented in Fig. 12. Around the time t ¼ 80 ms, the stress distribution appears similar to that in a statically loaded specimen (as shown in Fig. 12(d)). It can be seen that the damage occurs only after both the top and bottom ends of the specimen are loaded (at t ¼ 95 ms). Therefore, the failure pattern at this time (t ¼ 95 ms) is similar to that of the specimen under static loading. After t ¼ 95 ms, more cracks initiate at a location near the bottom side of the disk at time t ¼ 98 ms because the compressive stress wave is enhanced when the reflected wave also arrives this location. Then, the cracks propagate upwards around the vertical diameter of the disk. The failure pattern at time t ¼ 98 ms closely resembles to the experimental results of concrete obtained by Hughes et al. [40]. Besides, a large number of radially angled cracks that are located at the bottom half of the disk also initiate as time passes. As shown in Fig. 12, at t ¼ 1102130 ms, the stress concentration on the vertical diameter is released because the main crack has been
ARTICLE IN PRESS 246
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 10. The photoelastic fringes of maximum shear stress of homogeneous rock specimen subjected to incident compressive stress wave I as indicated in Fig. 9 (elastic stress field).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
247
Fig. 11. Failure process of rock specimen subjected to stress wave II (the brightness in the figures indicates the relative magnitude of maximum shear stress).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
248
formed in this zone (Fig. 12(i)-(l)). From the cracking process and the stress distribution as given in Fig. 11 and Fig. 12, it can be concluded that this failure pattern is mainly caused after the stress wave is reflected from the bottom fixed end of the loading plate. Under the stress wave III, as shown in Fig. 13, the cracks initiate much earlier than under stress wave II. Crack initiation occurs near the top loading platen at approximately t ¼ 50 ms because the high compressive stress concentration occurs directly when the applied incident compressive stress arrives at the top end of the rock disk, and the formed cracks grow downwards with a bifurcation occurring along the vertical diameter (t ¼ 60 ms). At this stage, many other cracks also initiate
at the top half of the specimen and propagate downwards. The failure of the specimen under this stress wave is initially caused by the incident compressive stress wave applied at the top loading plate. In addition, the stress wave reflected from the bottom end of loading plate also induces a large number of cracks at the bottom half of the specimen. Therefore, under this stress wave, the rock disk becomes much more fractured when it is compared to that under stress wave II. In Fig. 14, the variations of load and displacement at the top loading plate, as well as the distribution of damaged elements at different times, are shown. Under stress waves II and III, the peak loads applied on the top loading plates are 10.5 and 15.0 kN, respectively, which
65
65
45
σx
45
σx
25
σy
25
σy
5 -100
-15
5
0
100
200
300
400
500
-100
-35
-75
200
300
400
-75 Stress(MPa)
65
45
σx
45
σx
25
σy
25
σy
5
5
-15 0
100
200
300
400
500
-100 -15 0
-35
100
200
300
400
500
Step 950 (t = 95 µs)
-35
-55
-55
Step 400 (t = 40 µs)
-75
-75
Stress(MPa)
(b)
Stress(MPa)
(e)
65
65 45
σx
45
σx
25
σy
25
σy
5
5 -15
500
Step 800 (t = 80 µs)
(d)
Stress(MPa)
65
0
100
200
300
400
-100
500
-15 -35
-35 -55
0
100
200
300
400
500
Step 980 (t = 98 µs)
-55
Step 600 (t = 60 µs)
-75
-75 (c)
100
-55
Step 200 (t = 20 µs)
(a)
-100
0
-35
-55
-100
-15
Stress(MPa)
(f)
Stress(MPa)
Fig. 12. Profiles for the normal stresses along the vertical diameter at selected times when the rock specimen is subjected to incident stress wave II.
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
65
65
45
σx
45
σx
25
σy
25
σy
5 -100 -15 0
5 100
200
-35
300 400 500 Step 1000 (t = 100 µs)
-100 -15 0
100
200
-55
-75
400
Step 1200 (t = 120 µs)
(j)
65
Stress(MPa)
65
45
σx
45
σx
25
σy
25
σy
5 -100 -15 0
5 100
200
300
400
500
-100 -15 0
-35
100
200
300
400
Step 1100 (t = 110 µs)
-75
-55
Step 1250 (t = 125 µs)
-75
Stress(MPa)
(h)
(k)
Stress(MPa) 65
65 45
σx
45
σx
25
σy
25
σy
5
5 100
200
300
400
-100 -15 0
500
100
200
300
400
500
-35
-35 -55
500
-35
-55
-100 -15 0
500
-75
Stress(MPa)
(g)
300
-35
-55
-55
Step 1150 (t = 115 µs)
Step 1300 (t = 130 µs)
-75
-75 (i)
249
Stress(MPa)
(l)
Stress(MPa)
Fig. 12. (Continued)
are both 3 times larger than the peak load under the static loading condition. There has been a peak load of 7.5 kN applied at the top loading plate under the stress wave I, even though no damage is induced in the rock disk, which shows the high load-bearing capacity of rock under the dynamic loading condition. With regard to the displacement at the top loading plate, under the stress wave II the displacement has attained a peak value and begun to decline, while under stress III, it has not attained its peak value. Depending on the stress distribution that varies with time due to the propagation of the stress wave, the failure mechanisms of rock under the above two incident
compressive stress waves are different. When stress wave II is applied, as discussed above, the majority of the damaged elements are found after t ¼ 100 ms when the incident stress wave has arrived at the bottom loading plate and has been superposed with the reflected stress wave. Then another concentration of the damage occurs around t ¼ 100 ms at the bottom half of the rock disk. In contrast, under stress wave III, the large number of the damaged elements has already been found at about t ¼ 50 ms, and it remains stable until the t ¼ 130 ms. Compared to the numerical results of rock under static loading, the amount of cracking is significantly
ARTICLE IN PRESS 250
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 13. Failure process of rock specimen subjected to loading case III (the brightness in the figures indicates the relative magnitude of maximum shear stress).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252 1.2 Load
Displacement
40
1 35 0.8
30 25
0.6
20
0.4
15 10
0.2 5 0
Number of damaged elements
Displacement (mm) / load (*10 kN)
The following conclusions are drawn from our numerical simulations:
45 Damaged elements
0
0
10
20
30
40
50
60
70 80
90 100 110 120 130
Time (µs)
Stress wave II
(a)
60 Damaged elements
1.6
Load
Displacement
50 1.4 1.2
40
1 30 0.8 0.6
20
0.4 10 0.2 0 0
10
20
30
40
50
60
70 80
Number of damaged elements
Displacement (mm) / load (*10 kN)
1.8
0 90 100 110 120 130
Time (µs)
(b)
251
Stress wave III
Fig. 14. The variations of responded load, displacement at the top loading plate and the number of damaged elements during the failure process of the rock disk.
(1) Under static loading, numerical simulation reproduces well the primary crack at the disk center propagating along the diameter of specimen and the secondary cracks near the loading platens. (2) Under dynamic loading, different failure patterns and failure mechanisms are numerically captured when incident compressive stress waves with different amplitudes are applied to the rock specimen, which indicates that the failure pattern and failure mechanism of rock under dynamic loading is closely related to the propagation of stress waves. (3) Compared to the static loading condition, more cracks are developed in rock under dynamic loading, and the numerical results of this study are advantageous in capturing the failure pattern and failure mechanisms resulting from the propagation of stress wave in the rock specimen. In the future, care should be exercised in such simulations to accurately model the failure specimen subject to a variety of input stress pulses and to compare the results with experimental results. Although the work reported here has calculated on rock disks, such RFPA simulations are possible for a wide variety of specimen geometries.
Acknowledgements increased when rock is subjected to the dynamic loading. This is because much more damage is activated when the stress wave propagates in the rock. The cracks initiate, not only near the vertical diameter, but also in other parts of the rock disk. Under dynamic loading conditions, the rock specimen is much more fractured with the increase in magnitude of the stress waves. In contrast to the experimental results of the SHPB experiments, which were observed similar to a statically loaded specimen [31,40,41], the numerical simulation of this study is therefore advantageous in capturing the failure process and full failure mechanism during the propagation of a stress wave in the heterogeneous rock.
5. Summary and conclusions A numerical simulator, which considers the heterogeneity of material properties and the effect of strain rate on strength, was used to study the failure processes of a rock disk subjected to static and dynamic loading.
The authors would like to thank the financial support from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. PolyU 5044/99E) from 2001 to December 2002. They also thank the Natural Science Foundation of China (Grant Nos. 50504005 and 50204003) for their continuing support.
References [1] Lindholm US, Yeakley LM, Nagy A. The dynamic strength and fracture properties of Dreser basalt. Int J Rock Mech Min Sci Geomech Abstr 1974;11:181–91. [2] Shockey DA, Curran DR, Seaman L, Rosenberg JT, Petersen CF. Fragmentation of rock under dynamic loads. Int J Rock Mech Min Sci Geomech Abstr 1974;11:303–17. [3] Hanson ME. Some effects of macroscopic flaws on dynamic fracture patterns near a pressure-driven fracture. Int J Rock Mech Min Sci Geomech Abstr 1975;12:311–23. [4] Peng SS. A note on the fracture propagation and time-dependent behavior of rocks in uniaxial tension. Int J Rock Mech Min Sci 1975;12:125–7. [5] Janach W. The role of bulking in brittle failure of rock under rapid compression. Int J Rock Mech Min Sci 1976;13:177–86.
ARTICLE IN PRESS 252
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
[6] Chong KP, Boresi AP. Strain rate dependent mechanical properties of new Albany reference shale. Int J Rock Mech Min Sci Geomech Abstr 1990;27:199–205. [7] Ahrens TJ, Rubin AM. Impact-induced tensional failure in rock. J Geophys Res 1993;98(E1):1185–203. [8] He HL, Ahrens TJ. Mechanical properties of shock-damaged rocks. Int J Rock Mech Min Sci Geomech Abstr 1994;31(5): 525–33. [9] Zhao J, Li HB, Wu MB, Li TJ. Dynamic uniaxial compression test on granite. Int J Rock Mech Min Sci 1999;36(2):273–7. [10] Zhang ZX, Kou SQ, Yu J, Yu Y, Jiang LG, Lindqvist PA. Effects of loading rate on rock fracture. Int J Rock Mech Min Sci 1999;36(5):597–611. [11] Frew DJ, Forrestal MJ, Chen W. A split Hopkinson pressure bar technique to determine compressive stress-strain data for rock materials. Exp Mech 2001;41(1):40–6. [12] Grady DE, Kipp ME. The micromechanics of impact fracture of rock. Int J Rock Mech Min Sci 1979;16:293–302. [13] Zhao J. Application of Mohr–Coulomb and Hoek–Brown strength criteria to the dynamic strength of brittle rock. Int J Rock Mech Min Sci 2000;37:105–12. [14] Taylor LM, Chen EP, Kuszmaul JS. Microcrack-induced damage accumulation in brittle rock under dynamic loading. Comput Method Appl Mech Eng 1986;55:301–20. [15] Fahrenthold EP. A continuum damage model for fracture of brittle solids under dynamic loading. J Appl Mech ASME 1991;58:904–9. [16] Yang R, Bawden WF, Katsabanis PD. A new constitutive model for blast damage. Int J Rock Mech Min Sci Geomech Abstr 1996;33(3):245–54. [17] Yazdchi M, Valliappan S, Zhang WH. Continuum model for dynamic damage evaluation of anisotropic brittle materials. Int J Numer Methods Eng 1996;39(9):1555–83. [18] Liu LQ, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci 1997;34(2): 217–31. [19] Li HB, Zhao J, Li TJ, Yuan JX. Analytical simulation of the dynamic compressive strength of granite using the sliding crack model. Int J Numer Anal Methods Geomech 2001;25:853–69. [20] Huang CY, Subhash G, Vitton SJ. A dynamic damage growth model for uniaxial compressive response of rock aggregates. Mech Mater 2002;34:267–77. [21] Homand-Etienne F, Hoxha D, Shao JF. A continuum damage constitutive law of brittle rocks. Comput Geotechn 1998;22(2): 135–51. [22] Brara A, Camborde F, Klepaczko JR, Mariotti C. Experimental and numerical study of concrete at high strain rates in tension. Mech Mater 2001;33:33–45. [23] Sanchidrian JA, Pesquero JM, Garbayo E. Damage in rock under explosive loading: implementation in DYNA2D of a TCK model. Int J Surf Min Reclamat 1992;6(3):109–14.
[24] Tang CA. Numerical simulation of progressive rock failure and associated seismicity. Int J Rock Mech Min Sci 1997;34:249–61. [25] Tang CA, Liu H, Lee PKK, Tsui Y, Tham LG. Numerical tests on micro-macro relationship of rock failure under uniaxial compression, part I: effect of heterogeneity. Int J Rock Mech Min Sci 2000;37:555–69. [26] Zhu WC, Tang CA. Numerical simulation on shear fracture process of concrete using mesoscopic mechanical model. Construct Build Mater 2002;16(8):453–63. [27] Zhu WC, Tang CA. Micromechanical model for simulating the fracture process of rock. Rock Mech Rock Eng 2004;37(1):25–56. [28] Chau KT, Zhu WC, Tang CA. Numerical simulation of failure of brittle solids under dynamic impact using a new computer program—DIFAR. Key Eng Mater 2004;261–263:239–44. [29] Zhu WC, Tang CA, Huang ZP, Liu J. A numerical study of the effect of loading conditions on the dynamic failure of rock. Int J Rock Mech Min Sci 2004;41(3):424. [30] Weibull W. A statistical distribution function of wide applicability. J Appl Mech 1951;18:293–7. [31] Tedesco JW, Ross CA, Mcgill PB, O’Neil BP. Numerical analysis of high strain rate concrete direct tension tests. Comput Struct 1991;40(2):313–27. [32] Rinehart JS, Pearson J. Behavior of metals under impulsive loads. Scranton, Pennsylvania: The Haddon Craftsmen, Inc.; 1954. [33] Williams RE, Potter RM, Miska S. Experiments in thermal spallation of various rocks. J Energy Resour Technol ASME 1996;118(1):2–8. [34] Cho SH, Ogata Y, Kaneko K. Strain-rate dependency of the dynamic tensile strength of rock. Int J Rock Mech Min Sci 2003;40:763–77. [35] Dı´ az-Rubio FG, Pe´rez JR, Ga´lvez VS. The spalling of long bars as a reliable method of measuring the dynamic tensile strength of ceramics. Int J Impact Eng 2002;27:161–77. [36] Rocco C, Guinea GV, Planas J, Elices M. Mechanisms of rupture in splitting tests. ACI Mater J 1999;96(1):52–60. [37] Hondros G. The evaluation of Poisson’s ratio and the modulus of materials of a low tensile resistance by the Brazilian (indirect tensile) test with particular reference to concrete. Aust J Appl Sci 1959;10:243–68. [38] Yanagidani T, Sano O, Terada M, Ito I. The observation of cracks propagating in diametrically-compressed rock disks. Int J Rock Mech Min Sci Geomech Abstr 1978;15:225–35. [39] Hudson JA, Brown ET, Rummel F. The controlled failure of rock disks and rings loaded in diametral compression. Int J Rock Mech Min Sci 1972;9(2):241–8. [40] Hughes ML, Tedesco JW, Ross CA. Numerical analysis of high strain rate splitting-tensile tests. Comput Struct 1993;47(4/5): 653–71. [41] Gomez JT, Shukla A, Sharma A. Static and dynamic behavior of concrete and granite in tension with damage. Theoret Appl Fract Mech 2001;36:37–49.
International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252 www.elsevier.com/locate/ijrmms
Numerical simulation of Brazilian disk rock failure under static and dynamic loading W.C. Zhu, C.A. Tang Center for Rock Instability and Seismicity Research, Box 138, Northeastern University, Shenyang 110004, PR China Accepted 27 June 2005 Available online 22 August 2005
Abstract A numerical simulator based on RFPA (Rock Failure Process Analysis) is used to study the deformation and failure process of a Brazilian disk of heterogeneous rock when subjected to static and dynamic loading conditions. In this simulator, the heterogeneity of rock is considered by assuming that the material properties of elements conform to a Weibull distribution; an elastic damagebased law that considers the strain-rate dependency is used to describe the constitutive law at mesoscopic scale; and a finite element program is employed as a basic stress analysis tool. The simulator is firstly validated by simulating the dynamic spalling of a homogeneous rock bar and by comparing with the theoretical and experimental results. Then, the failure process of a Brazilian disk of rock subjected to static and dynamic loading is numerically simulated, and the numerical results are compared with the available experimental results. Particular attention is given to the typical failure patterns of the rock disk when the incident compressive stress waves with different amplitudes are applied. The numerical simulation also identifies the failure mechanisms of rock during dynamic failure processes that are closely related to the propagation of the stress wave. r 2005 Elsevier Ltd. All rights reserved. Keywords: Rock failure; Brazilian disk; Dynamic loading; Numerical simulation
1. Introduction Better understanding of the dynamic fracture processes of rocks promises benefit in many areas from rock mechanics to mining engineering and earthquake prediction. It is essential to understand how fractures initiate and propagate under different loading conditions in order to provide understanding of rock fracture process that occur in the engineering fields. During the past several decades, researchers have carried out many dynamic tests of rock using the apparatus such as the Split Hopkinson pressure bar (SHPB). From these experimental investigations, it has been understood that the mechanical properties (or responses) and failure characteristics of rock are sensitive to the strain rate. It is also a well-established fact that rocks exhibit increases Corresponding author. Tel.: 86 24 83687705; fax: 86 24 83671626.
E-mail address: [email protected] (W.C. Zhu). 1365-1609/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijrmms.2005.06.008
in strength as the loading rate or strain rate increases [1–11]. Based on these experimental results, failure criteria considering the strain rate effect of rocks are proposed [12,13]. However, not so well understood is the dynamic failure process and failure mechanism, which are closely related to the microstructural behavior in terms of activation, growth, and coalescence of cracks when the rock is subjected to different loading conditions. With regard to the SHPB experiments, they are widely used to determine strain-rate dependency of the dynamic strength of rock, but not to study the failure process and failure mechanisms that are closely related to the stress wave propagation in heterogeneous rock. Continuum damage mechanics is considered to be appropriate to describe the failure process resulting from dynamic loading. Based on this approach, many theoretical models were developed to study the dynamic damage evolution of brittle materials such as rock and
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Nomenclature
A b C D E, E0 f c0 , f t0 f cr , f tr f cs0 k K m M p
parameter to reflect the stress rate effect on strength width of load-bearing strip damping matrix diameter of rock disk damaged and undamaged (initial) elastic modulus of element dynamic compressive and tensile strengths, respectively, of element residual dynamic compressive and tensile strengths, respectively, of element static uniaxial compressive strength ratio of tensile and compressive strength stiffness matrix homogeneity index, shape parameter of Weibull distribution mass matrix applied incident stress
concrete with microstructures [12,14–20]. Most of these models are developed based on combining the theory of damage mechanics or fracture mechanics with some statistic treatment to account for the random distribution of microcracks. Because many micromechanical parameters such as crack density must be used to describe the crack distribution in these damage models, it is often difficult to implement such models in a numerical code [21]. Therefore, there is a real need for further developments in the material modeling of rock to reflect more closely the physical mechanisms responsible for failure. In the numerical simulation with the discrete element method, the basic fracture patterns of rock and concrete subject to dynamic tension are captured [22]. But in the discrete element model, there is no accepted general method of averaging the contact displacement between particles in order to obtain the value of resulting macroscopic strain or stress [22]. The explicit non-linear finite element code DYNA 2D/3D is used to analyze a variety of dynamic problems, which is generally concentrated on collision analysis, such as an elastic sphere on an elastic–plastic metal material and sheet metal forming. Although DYNA 2D had been also used to calculate the response of jointed rock subjected to blasting [14,23], it was not focused on capturing the dynamic failure mechanism of rock. In fact, the stress distribution and damage initiation in the rock is closely related to stress wave propagation; moreover, the local damage in rock can in turn affect the
237
pmax R u
amplitude of stress wave vector of externally applied loads parameter of elements that conforms to Weibull distribution, such as elastic modulus and strength _ Udisplacement, € U, U, velocity, and acceleration vectors of elements u0 scale parameter of Weibull distribution e strain e0 strain at peak stress e1, e2, e3 principal strain ec0 strain at the peak compressive stress et0 strain at the peak tensile stress etu ultimate tensile strain ¯ equivalent strain f internal friction angle Z ultimate tensile strain coefficient ðZ ¼ tu =t0 Þ, l residual strength coefficient ðl ¼ f tr =f t0 Þ m Poisson’s ratio s stress s1 , s2 , s3 principal stress o damage variable
propagation of the stress wave. Numerical simulation of rock subjected to dynamic loading conditions requires the implementation of a constitutive relation capable of capturing the major features of the failure process during the propagation of stress waves. The Rock Failure Process Analysis (RFPA) Code, which had been developed and introduced by Tang [24,25] in many documents, has been successfully used to simulate the fracture process and failure mechanism of rock under static loading. In RFPA, in order to numerically denote the heterogeneity of rock, the rock sample is discretized into many mesoscopic elements and their material properties are conformed to the Weibull distribution [26,27]. The constitutive law of these mesoscopic elements is formed based on elastic damage mechanics. In cooperation with Professor Chau [28] at the Department of Civil and Structural Engineering, The Hong Kong Polytechnic University, the author proposed a stress-rate dependent constitutive law and implemented it into the original RFPA code [24,25], and therefore the RFPA is extended for analyzing the failure process of rock under dynamic loading [28,29]. From January 2003 on, based on that, Zhu et al. [29] have made many further improvements in utilizing the faster finite element solver and introducing more boundary conditions, and all these have been also implemented into the original RFPA code [24,25]. Therefore, this upgraded simulator used in this study is also called RFPA. The contribution of this work is to validate the simulator by reproducing the dynamic spalling of
ARTICLE IN PRESS 238
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
homogeneous material bar, and then to investigate the failure processes and related failure mechanism of a Brazilian disk of rock under static and dynamic loading. Based on the numerical simulations, the different failure processes are captured, and the related failure mechanisms of rock under static and dynamic loading are identified.
2. Numerical simulator description Briefly, the simulator used in this work is a twodimensional (2D) code that can simulate the fracture and failure process of rock under static or dynamic loading conditions. To model the failure of rock material (or rock mass), the rock medium is assumed to be composed of many mesoscopic elements whose material properties are different from one to another and are specified according to a Weibull distribution [30]. The finite element method is employed to obtain the stress fields in the mesoscopic elements. Elastic damage mechanics is used to describe the constitutive law of the meso-scale elements when the maximum tensile strain criterion and the Mohr–Coulomb criterion that incorporate the effect of stress rate are utilized as damage thresholds [27].
(a)
m = 1.5 (The values range from 1.0 to 363.)
(b)
m = 5.0 (The values range from 25.1 to 147)
(c)
m = 10.0 (The values range from 50.1 to 121)
2.1. Assignment of material properties In RFPA, the solid or structure is assumed to be composed of many mesoscopic elements with the same size, and the mechanical properties of these elements are assumed to conform to a given Weibull distribution as defined by the following probability density function (pdf): m m u m1 u f ðuÞ ¼ exp , (1) u0 u0 u0 where u is the mechanical parameter of the element (such as strength or elastic modulus); the scale parameter u0 is related to the average of the element parameters and the parameter m defines the shape of the distribution function. From the properties of the Weibull distribution, a larger value of m implies a more heterogeneous material and vice versa. Therefore, the parameter m is called the homogeneity index in our numerical simulations. In the definition of the Weibull distribution, the value of the parameter m must be larger than 1.0. Using the pdf, in a computer simulation of a medium composed of many mesoscopic elements, one can produce numerically a heterogeneous material. The computationally produced heterogeneous medium is analogous to a real specimen tested in the laboratory, so in this investigation it is referred to as a numerical specimen. Fig. 1 presents some 2D numerical specimens, which are all composed of 50 50 elements, produced randomly by computer according to the Weibull distributions with
Fig. 1. Numerical specimens with different homogeneity index m. The specimen is composed of 2500 square elements and the Weibull parameter u0 is 100 MPa. The intensity of the gray scale in the figure denotes the magnitude of strength of elements.
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
different homogeneity indices. In Fig. 1, the different degree of gray color corresponds to a different magnitude of strength of the elements. In general, it has been assumed that the Young’s modulus and the strength of the mesoscopic elements that were used to simulate a rock specimen conform to two individual distributions with the same heterogeneity indices. The mesoscopic elements themselves are isotropic and homogeneous. The mesoscopic elements in the specimen must be sufficiently small to reflect the heterogeneous mechanical properties of materials at the meso-scale and still provide conditions under which the computer can perform the analysis efficiently. 2.2. The constitutive model for the mesoscopic element Initially an element is considered elastic, with elastic properties defined by Young’s modulus and Poisson’s ratio. The stress–strain relation for an element is considered linear elastic until the given damage threshold is attained, and then is modified by softening. Under a dynamic stress state, the elements undergo damage when one of the following damage criteria is satisfied at the element level: 1 ¼ kf c0 =E 0 ; s1
1 þ sin f s3 Xf c0 , 1 sin f
(2)
where f c0 is the dynamic uniaxial compressive strength of the element, which is closely related to the strain rate (or stress rate) of dynamic loading condition. E0 is initial elastic modulus of the element that is assumed to be not affected by the stress rate, k is the ratio of compressive and tensile strength and f is the internal frictional angle of the element. s1 and s3 are the major and minor principal stresses of the element. In this study, the following relation between dynamic uniaxial compressive strength and loading rate, which has been proposed by Zhao [13], is used to reflect the effect of stress rate on the dynamic strength: _ f_cs0 Þ þ f cs0 ; f c0 ¼ A logðs=
_ f_cs0 , s4
239
be damaged in either tension (corresponding to the maximum tensile strain criterion) or shear (corresponding to the Mohr–Coulomb criterion). Once Eq. (2) is satisfied at the element level, the elastic modulus of the element is reduced according to the following expression: E ¼ ð1 oÞE 0 ,
(4)
where o represents a damage variable, E and E0 are the elastic moduli of the damaged and the undamaged element, respectively. In the current method, the element as well as its damage is assumed isotropic, so the E, E0 and o are all scalar quantities. The sign convention used throughout this paper is that compressive stress and strain are positive. When the mesoscopic element is in a uniaxial stress state (both uniaxial compression and uniaxial tension), the constitutive relation of elements is as illustrated in Fig. 2. Initially, the stress–strain curve is linear elastic and no damage exists, i.e. o ¼ 0. When the maximum tensile strain criterion is satisfied, damage occurs in the element in a brittle mode. The constitutive relation of the mesoscopic element under uniaxial tension, as shown in the third quartile of Fig. 2, can be expressed as 8 4t0 ; > < 0; lt0 o ¼ 1 ; tu opt0 ; (5) > : 1; ptu ; where l is the residual strength coefficient, which is given as f tr ¼ lf t0 and f t0 is the tensile strength of element. The parameter t0 is the strain at the elastic limit, which is called the threshold strain and σ
fc0
(3)
where f c0 is also the dynamic uniaxial compressive strength (MPa), s_ is the stress rate (MPa/s), f cs0 is the uniaxial compressive strength at the quasi-static stress rate f_cs0 that is approximately 5 102 MPa/s, and A is a parameter depending on the material. In addition, the experimental results of Zhao [13] also indicated that the ratio of tensile and compressive strength ðk ¼ f c0 =f t0 Þ and internal frictional angle f are not influenced by the stress rate. In this respect, when the f c0 obtained from Eq. (3) is substituted into Eq. (2), the effect of strain rate on the strength of elements will be incorporated. The first part of Eq. (2) is the maximum tensile strain criterion, whilst the second part is the classical Mohr–Coulomb failure criterion for tensile and shear damage thresholds, respectively. Thus, an element may
fcr ε tu
ε t0
- ftr
ε c0
ε
- ft0 Fig. 2. Elastic damage constitutive law of element under uniaxial stress state (Here ft0 and ftr are dynamic uniaxial tensile strength and residual uniaxial tensile strength of element, respectively. And fc0 and fcr are dynamic uniaxial compressive strength and residual corresponding strength of element, respectively.)
ARTICLE IN PRESS 240
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
calculated as t0 ¼ f t0 =E 0 .
(6)
From Eq. (6), it can be found that t0 also depends on the stress rate of element because f c0 is related to the stress rate according to Eq. (3). tu is the ultimate tensile strain of the element, describing the state at which the element would be completely damaged. The ultimate tensile strain is defined as tu ¼ Zt0 , where Z is called the ultimate strain coefficient. Both the residual strength coefficient (l) and ultimate strain coefficient (Z) are assumed to be independent of the stress rate. Additionally, it is assumed that the damage of a mesoscopic element under multiaxial states of stress is also isotropic and elastic. The constitutive equation described above can be extended for application to 3D stress states when the tensile strain threshold is attained. In multiaxial stress states, the element is still subject to damage in the tensile mode when the equivalent major tensile strain ¯ attains the above threshold strain t0 . The equivalent principal strain ¯ is defined by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (7) ¯ ¼ h1 i2 þ h2 i2 þ h3 i2 , where 1 , 2 and 3 are three principal strains, and o4is a function defined as follows: ( x; xX0; ox4 ¼ (8) 0; xo0: For an element subject to a multiaxial state of stress, the constitutive law can be easily obtained by substituting the strain e in Eq. (5) with the equivalent strain ¯ defined by Eqs. (7) and (8). The damage variable is expressed as 8 ¯4t0 ; > < 1; lt0 (9) o ¼ 1 ¯ ; tu o¯pt0 ; > : 1; ¯ptu : In the same way as for uniaxial tension, when the element is under uniaxial compression but damaged according to the Mohr–Coulomb criterion, the expression for damage variable o can be described as ( 0; oc0 ; o¼ (10) 1 lc0 ; Xc0 ; where l is also the residual strength coefficient, and f cr =f c0 ¼ f tr =f t0 ¼ l is assumed to apply when the element is under either uniaxial compression or tension. Also, c0 is calculated as c0 ¼ f c0 =E 0 .
(11)
If an element is under a multiaxial stress state and its strength satisfies the Mohr–Coulomb criterion, damage occurs, and it is necessary to consider the effect of the other principal stress in the model during damage
evolution. When the Mohr–Coulomb criterion is satisfied, the maximum principal (compressive) strain c0 can be calculated at the peak value of the maximum principal (compressive) stress: 1 1 þ sin f s3 mðs1 þ s2 Þ . f c0 þ (12) c0 ¼ E0 1 sin f It is assumed that shear damage evolution is related only to the maximum compressive principal strain 1 . So the maximum compressive principal strain 1 of the damaged element is substituted for the uniaxial compressive strain e in Eq. (10), extending the equation to represent shear damage under triaxial stress states: ( 0; 1 oc0 ; o¼ (13) lc0 1 1 ; 1 Xc0 : Based on the above equations, it can be found that the damage variable o is related to the stress rate because it is calculated based on the values of c0 or t0 that are obtained directly according to the stress-rate dependent strength of the element ðf c0 or f t0 Þ according to Eq. (6) or (11). Of course, the above-mentioned constitutive law can be used to simulate the failure process of rock under static loading when the stress rate effect on the strength as expressed by Eq. (3) is not incorporated. Therefore, the stress-rate dependent constitutive law described above is based on the original one for static analysis and is also compatible with it. 2.3. Finite element implementation As described above, the rock specimen is composed of many square elements with the same size. These elements are also acting as the four-noded iso-parametric elements for finite element analysis. The equilibrium equations governing the linear dynamic response of a system of finite elements can be expressed in the following form: € þ CU _ þ KU ¼ R, MU
(14)
where M, C and K are the mass, damping, and stiffness matrices; R is the vector of externally applied loads; and _ and U € are the displacement, velocity, and U, U acceleration vectors of the finite elements. A direct step-by-step integration procedure is found suitable for solving the problem in which a body is subjected to a short duration impulse loading [31]. At initial condition, the elements are elastic, and their stresses can be calculated via step-by-step integration. At a time step, the principal stresses from the previous time step are subtracted from the current principal stresses and divided by the time step in order to calculate the stress rates of the elements. Similarly, when the current minor principal stress and corresponding stress
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
rate are both negative, the increased strength of an element due to the increase of absolute value of this stress rate can also be obtained. When the effect of stress rate on the dynamic compressive strength is considered according to Eq. (3) and the stresses of (or strains) of elements meet the maximum tensile strain criterion or Mohr–Coulomb criterion as given in Eq. (2), the elements are damaged according to the constitutive law given before, then their stress will be re-analyzed iteratively for the current boundary conditions in order to reflect the stress redistribution at this time step. The program does not consider the analysis of the next time step until no new damaged elements are found for the iterative step of this time step. The above model described in the two preceding sections is implemented into the RFPA [24,25], and therefore RFPA is extended for analysis of rock under dynamic loading.
241
which means that no strain rate effect is considered in this numerical simulation. The bar is 200 mm in length, 5 mm in width, and it is simulated with 200 5 elements in 2D plane stress problem. The incident compressive stress wave is applied at the top end of the bar and the bottom end of it is kept free. Two compressive incident stress waves as shown in Fig. 3(b) are applied on the top surface of bar, in order to illustrate the effect of stress amplitudes on spalling. In Fig. 4, the stress distributions during the fracture process of the bar subjected to two different stress waves are given, where the brightness of the gray scale indicates the relative magnitude of the maximum shear stress in the mesoscopic elements. From this Figure, the propagation of the stress wave can be seen clearly. For incident compressive stress wave I (pmax ¼ 20 MPa), as shown in Fig. 4(a), spalling occurs, induced when the incident compressive stress wave travels down and is reflected from the bottom free end of the bar, and
3. Simulator validation As discussed by many researchers [32–35], the dynamic spalling of quasi-brittle material such as rock occurs when an incident compressive wave is reflected into a tensile one on reaching a free end. If the incident compressive stress is lower than the compressive strength but its absolute value exceeds the material dynamic tensile strength, the specimen will spall during the reflection process. In this section, the dynamic spalling of a homogeneous quasi-brittle material bar is studied in order to validate RFPA in simulating the dynamic failure process. The geometries and loading conditions for the bar specimen are shown in Fig. 3(a). The mechanical properties of this bar are homogeneous and are specified as: the Young’s modulus is 62.0 GPa, uniaxial compressive strength is 200 MPa, direct tensile strength is 19.0 MPa, and Poisson’s ratio is 0.20. The parameter A for this homogeneous bar is specified to be zero,
(10µs)
(20µs)
(a)
Stress wave I (Peak value of incident stress is 50 MPa)
(30µs)
(40µs)
(50µs)
(60µs)
(70µs)
(80µs)
(90µs)
(100µs)
p
100
10 mm
p (MPa)
200 mm
75
II
50 I 25 (10µs)
(a)
(b)
5
10
(b)
(20µs)
(30µs)
(40µs) (42.5µs)
(43µs)
(45µs)
(50µs)
(60µs)
(70µs)
Stress wave II (Peak value of incident stress is 100 MPa)
t(µs)
Fig. 3. (a) Geometries and loading conditions for homogeneous bar subjected to dynamic loading. (b) The incident compressive stress waves with different peak amplitudes applied on the top surface of bar.
Fig. 4. Dynamic spalling of homogeneous bar is numerically simulated. In (a) the amplitude of incident compressive stress wave is 20 MPa; and in (b) the amplitude of incident compressive stress wave is 40 MPa.
ARTICLE IN PRESS 242
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
reaches the tensile strength of the material. The fracture face is at 22 mm above the bottom end of the bar. Finally, the stress waves in two parts of the bar only propagate individually with a certain velocity and no further spalling is induced. Two spallings occur for incident compressive stress wave II (pmax ¼ 40 MPa) when the level of the stress wave becomes more than twice the tensile strength of the material. As shown in Fig. 4(b), the first spall fracture, upon its abrupt creation at 12 mm above the bottom free end of bar, places a free surface in front of the trailing portion of the original wave. When the material fails, the fraction of the tensile wave that has passed the
failure point, continues traveling along the specimen, and it causes another spalling since the magnitude of tensile stress is also larger than the tensile strength of the bar. The second fracture occurs 25 mm away from the original bottom end of the bar. No more spallings occur because the intensity of the rear portion of the wave is no longer greater than the tensile strength of the material. This numerical result compares well with the theoretical results obtained by Rinehart [32] based on 1D wave propagation theory and the experimental observation given by Brara et al. [22]. From the above numerical simulation, it can be concluded that RFPA is effective in capturing the dynamic failure process of quasi-brittle material subjected to dynamic loading in these circumstances.
y
p
4. Numerical simulations of static and dynamic failure 4.1. Numerical rock disk rock Steel plate
D =150 mm
Rock disk x
25
2
b
Steel plate
Fig. 5. Geometries and loading conditions for disk specimen subjected to static and dynamic loading (D ¼ 150 mm, b=D ¼ 0:12).
The geometries and loading conditions for disk specimens are shown in Fig. 5. In our numerical simulation, the rock has a Young’s modulus of 62.2 GPa, a uniaxial compressive strength of 150.0 MPa, a direct tensile strength of 15.0 MPa, a mass density of 2500 kg/m3, and a Poisson’s ratio of 0.28. The parameter A of the material is specified to be 13.0, which is given according to the uniaxial compressive strength of this rock. When the Weibull distribution parameters as listed in Table 1, are applied to a standard numerical specimen that is 200 100 mm in dimensions and composed of 200 100 elements, the standard numerical specimen will have the above-mentioned mechanical properties of this rock when subjected to static loading. The rock disk is 150 mm in diameter, and in all cases, it is simulated with 150 150 elements as a 2D plane stress problem. Since we assume that the cushion materials (load-bearing strip) and the steel plate will not fail in numerical simulations, they are specified to be homogeneous and have a very high strength. Here the ratio between the width of load-bearing strip and diameter of disk (b/D) is selected to be of 0.12 (i.e. b=D ¼ 0:12) in order to allow the specimen to statically fail in tension along the loading line and not locally by compression [36].
Table 1 Weibull distribution parameters used for RFPA to simulate rock Materials
Homogeneity index m
Mean of Young’s modulus (GPa)
Mean of uniaxial compressive strength (MPa)
Poisson’s ratio
Ratio of compressive and tensile strength
Rock Steel plate Cardboard
5 Homogeneous Homogeneous
70 210 10
510 800 500
0.28 0.3 0.3
12 1 1
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
4 3.5
113
3 Load (kN)
For the static loading condition, displacement is imposed on the top loading paten step by step with the bottom platen fixed in the vertical direction, which simulates the compressive loading with displacement control. Similarly, for the dynamic loading conditions, the incident compressive stress is imposed at the top loading platen while the bottom platen is also fixed in the vertical direction.
2.5
77
2 1.5 1
180 23
0.5
4.2. Under static loading
243
114
0
Initially, the rock specimen is elastic and the stress distribution on the vertical diameter is given in Fig. 6, and it is compared to theoretical results for a specific width of load-bearing strips [37]. In the figure, the stress is normalized with 2P=ðpDtÞ (here P is the applied load and t the unit thickness for plane stress analysis). It can be seen that the tensile stress distribution at the center area of the rock specimen shows good agreement with the theoretical results. Because the Young’s moduli of elements in the numerical specimen are assumed to be heterogeneous, the stress distribution across the loaded diameter is also heterogeneous and fluctuates around the theoretical curves. Fig. 7 shows the load response of the rock specimen during the static failure process. The load-step curve is basically linear until peak loading is attained. Then the specimen loses its loading capacity substantially because this specimen is relatively homogeneous. The failure patterns at several typical loading steps (as indicated on the curve of Fig. 7) are shown in Fig. 8. In this figure, the distribution of damaged elements, distribution of maximum shear stress and cracking pattern in the rock specimen subjected to static loading are presented. In Fig. 8(a), the damaged elements are denoted with different color, i.e. white and dark gray for elements damaged in shear and tensile at current step, respectively, and black for damaged elements at all the
1 Theoretical (Hondros) σx 0. 5
-10
0
Numerical σx 0
10
20
σy σy
0
20
40
60
80
100 120 140 160 180 200 Step
Fig. 7. Load-step curves of specimens subjected to static loading (numerical results, the applied displacement at the top loading platen is 0.002 mm/step).
previous steps. In Fig. 8(b) and (c), the brightness indicates the relative magnitude of maximum shear stress and elastic modulus of elements, respectively. Since the damage of elements leads to degradation of its elastic modulus, the cracking patterns can be clearly given from the elastic modulus figure of elements. The first element is damaged in tensile mode near the vertical diameter of the specimen (as shown at Step 23 of Fig. 8). With the increase of the external loading, a damage zone near the vertical diameter is developed due to the tensile damage of many elements (as Steps 77 and 113 of Fig. 8). Just after the peak load the damaged elements in the damage zone coalescence with each other to form the primary crack along the vertical diameter of specimen. At the same loading step, the secondary diagonal cracks near two loading platens also develop. It should be mentioned that, for this specimen, fracture initiation does not mean collapse of the specimen. From Fig. 8, it can be seen that the initiation of fractures (at Step 23 as shown in Fig. 7) usually occurs when the load is 20% of the peak load (at Step 113). Also, a damage zone that is distributed in a wide range around the vertical diameter of disk is developed, which agrees well with the acoustic emission observation of Yanagidani et al. [38]. Besides, due to the heterogeneity, the path of the primary crack in the rock disk, which is not a straight line, is just like the common phenomena that are observed in many laboratory experiments [36–39].
30
4.3. Under dynamic loading -0.5
-1 Fig. 6. Theoretical and numerical results of stresses distribution across loaded diameter for disk specimen under indirect test.
In an attempt to more fully comprehend the response of the rock disk to the different amplitudes of incident compressive stress waves, the three incident compressive stress waves as shown in Fig. 9 are considered in the numerical simulations. Three incident compressive stress waves with different amplitudes (pmax) are applied on
ARTICLE IN PRESS 244
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 8. Failure process of rock disk subjected to static loading (The applied displacement at the top loading platen is 0.002 mm/step. Damage initiates at step 23, and load attains its peak value at step 113.)
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
III
p (MPa)
100 80
II
60
I
40 20 10
20
30
40 t (µs)
50
60
Fig. 9. Incident compressive stress waves imposed on the top loading platen above rock specimen: (I) pmax ¼ 50 MPa; (II) pmax ¼ 70 MPa; (III) pmax ¼ 100 MPa.
the top loading plates (as shown in Fig. 5). The amplitudes of the three incident compressive stress waves are 50, 70 and 100 MPa, respectively. 4.3.1. Elastic stress wave Firstly, the elastic stress field for the homogeneous specimen under the stress wave I is given, in order to illustrate the stress wave propagation in the homogeneous rock disk. The material parameters are specified for the material properties of this rock, i.e. the Young’s modulus is 62.2 GPa, Poisson’s ratio is 0.28 and mass density is 2500 kg/m3. The photoelastic fringes of maximum shear stress at different times are obtained from numerical simulation and are shown in Fig. 10. It takes t ¼ 3 ms for the stress wave to arrive the top end of the rock disk and it is partially reflected back into the top steel plate and partially transmitted into the rock specimen. Then, the stress wave propagates downwards and arrive at the bottom end of the rock disk around t ¼ 35 ms, at which time the stress wave in the rock specimens will also be partially reflected back into the rock specimen and partially transmitted into the bottom steel plate. About 3 ms later, the stress wave that is transmitted into the bottom steel plate reaches the bottom fixed end of the steel plate and is mostly reflected because the bottom surface of the loading plate is fixed in the vertical direction. Under the superposition of incident and reflected stress waves in the rock disk, the stress distribution in the homogeneous disk appear identical to a statically loaded specimen at t ¼ 85 ms, which means that the loads at both contacts points of the specimen become equal and the specimen is in equilibrium at this time. However, the stress deviates from this transient equilibrium condition as time passes because the stress in the lower half part of the disk become relatively high due to the superposition of the incident stress wave and the reflected one after the reflection at the bottom fixed end. However, another equilibrium condition is also achieved at about
245
t ¼ 179 ms. Therefore, it can be predicted that a periodical stress distribution in the rock disk is achieved under this dynamic loading condition. The ‘photoelastic fringe’ of stress distribution obtained from the numerical simulation indicates that the stress distribution in a rock disk is not so regularly distributed as that under the static loading, but changes periodically with time. Before the equilibrium condition at t ¼ 85 ms, the top end of the rock disk has the larger compressive stress, while the bottom end of the rock disk is unstressed. In contrast, after t ¼ 85 ms, the relative high stress is induced at the bottom end of the rock disk. This means that the locations with the maximum stresses are shifted to the bottom end of the rock disk. Therefore, the stress distribution in the rock disk changes as the time passes, which makes the deformation and failure process of rock under dynamic loading much complicated. It should also be mentioned that the rock disk under this dynamic loading conditions is quite different from that in the SHPB experiments because the rock specimen in the SHPB experiments remains in equilibrium until the time of failure. Hence, it can be seen that the numerical simulations of this work are most useful for studying the complex failure process and failure mechanisms of rock under a variety of dynamic loading conditions. 4.3.2. Failure process and failure mechanisms Since no damage of elements occurs when stress wave I is applied, the failure process analysis is only concentrated on the other two incident compressive stress waves. The maximum shear stress distribution, as well as cracking sequence from beginning to final failure, for the rock disk under stress wave II is shown in Fig. 11. Profiles for the normal stresses along the vertical diameter at selected times are also presented in Fig. 12. Around the time t ¼ 80 ms, the stress distribution appears similar to that in a statically loaded specimen (as shown in Fig. 12(d)). It can be seen that the damage occurs only after both the top and bottom ends of the specimen are loaded (at t ¼ 95 ms). Therefore, the failure pattern at this time (t ¼ 95 ms) is similar to that of the specimen under static loading. After t ¼ 95 ms, more cracks initiate at a location near the bottom side of the disk at time t ¼ 98 ms because the compressive stress wave is enhanced when the reflected wave also arrives this location. Then, the cracks propagate upwards around the vertical diameter of the disk. The failure pattern at time t ¼ 98 ms closely resembles to the experimental results of concrete obtained by Hughes et al. [40]. Besides, a large number of radially angled cracks that are located at the bottom half of the disk also initiate as time passes. As shown in Fig. 12, at t ¼ 1102130 ms, the stress concentration on the vertical diameter is released because the main crack has been
ARTICLE IN PRESS 246
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 10. The photoelastic fringes of maximum shear stress of homogeneous rock specimen subjected to incident compressive stress wave I as indicated in Fig. 9 (elastic stress field).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
247
Fig. 11. Failure process of rock specimen subjected to stress wave II (the brightness in the figures indicates the relative magnitude of maximum shear stress).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
248
formed in this zone (Fig. 12(i)-(l)). From the cracking process and the stress distribution as given in Fig. 11 and Fig. 12, it can be concluded that this failure pattern is mainly caused after the stress wave is reflected from the bottom fixed end of the loading plate. Under the stress wave III, as shown in Fig. 13, the cracks initiate much earlier than under stress wave II. Crack initiation occurs near the top loading platen at approximately t ¼ 50 ms because the high compressive stress concentration occurs directly when the applied incident compressive stress arrives at the top end of the rock disk, and the formed cracks grow downwards with a bifurcation occurring along the vertical diameter (t ¼ 60 ms). At this stage, many other cracks also initiate
at the top half of the specimen and propagate downwards. The failure of the specimen under this stress wave is initially caused by the incident compressive stress wave applied at the top loading plate. In addition, the stress wave reflected from the bottom end of loading plate also induces a large number of cracks at the bottom half of the specimen. Therefore, under this stress wave, the rock disk becomes much more fractured when it is compared to that under stress wave II. In Fig. 14, the variations of load and displacement at the top loading plate, as well as the distribution of damaged elements at different times, are shown. Under stress waves II and III, the peak loads applied on the top loading plates are 10.5 and 15.0 kN, respectively, which
65
65
45
σx
45
σx
25
σy
25
σy
5 -100
-15
5
0
100
200
300
400
500
-100
-35
-75
200
300
400
-75 Stress(MPa)
65
45
σx
45
σx
25
σy
25
σy
5
5
-15 0
100
200
300
400
500
-100 -15 0
-35
100
200
300
400
500
Step 950 (t = 95 µs)
-35
-55
-55
Step 400 (t = 40 µs)
-75
-75
Stress(MPa)
(b)
Stress(MPa)
(e)
65
65 45
σx
45
σx
25
σy
25
σy
5
5 -15
500
Step 800 (t = 80 µs)
(d)
Stress(MPa)
65
0
100
200
300
400
-100
500
-15 -35
-35 -55
0
100
200
300
400
500
Step 980 (t = 98 µs)
-55
Step 600 (t = 60 µs)
-75
-75 (c)
100
-55
Step 200 (t = 20 µs)
(a)
-100
0
-35
-55
-100
-15
Stress(MPa)
(f)
Stress(MPa)
Fig. 12. Profiles for the normal stresses along the vertical diameter at selected times when the rock specimen is subjected to incident stress wave II.
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
65
65
45
σx
45
σx
25
σy
25
σy
5 -100 -15 0
5 100
200
-35
300 400 500 Step 1000 (t = 100 µs)
-100 -15 0
100
200
-55
-75
400
Step 1200 (t = 120 µs)
(j)
65
Stress(MPa)
65
45
σx
45
σx
25
σy
25
σy
5 -100 -15 0
5 100
200
300
400
500
-100 -15 0
-35
100
200
300
400
Step 1100 (t = 110 µs)
-75
-55
Step 1250 (t = 125 µs)
-75
Stress(MPa)
(h)
(k)
Stress(MPa) 65
65 45
σx
45
σx
25
σy
25
σy
5
5 100
200
300
400
-100 -15 0
500
100
200
300
400
500
-35
-35 -55
500
-35
-55
-100 -15 0
500
-75
Stress(MPa)
(g)
300
-35
-55
-55
Step 1150 (t = 115 µs)
Step 1300 (t = 130 µs)
-75
-75 (i)
249
Stress(MPa)
(l)
Stress(MPa)
Fig. 12. (Continued)
are both 3 times larger than the peak load under the static loading condition. There has been a peak load of 7.5 kN applied at the top loading plate under the stress wave I, even though no damage is induced in the rock disk, which shows the high load-bearing capacity of rock under the dynamic loading condition. With regard to the displacement at the top loading plate, under the stress wave II the displacement has attained a peak value and begun to decline, while under stress III, it has not attained its peak value. Depending on the stress distribution that varies with time due to the propagation of the stress wave, the failure mechanisms of rock under the above two incident
compressive stress waves are different. When stress wave II is applied, as discussed above, the majority of the damaged elements are found after t ¼ 100 ms when the incident stress wave has arrived at the bottom loading plate and has been superposed with the reflected stress wave. Then another concentration of the damage occurs around t ¼ 100 ms at the bottom half of the rock disk. In contrast, under stress wave III, the large number of the damaged elements has already been found at about t ¼ 50 ms, and it remains stable until the t ¼ 130 ms. Compared to the numerical results of rock under static loading, the amount of cracking is significantly
ARTICLE IN PRESS 250
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
Fig. 13. Failure process of rock specimen subjected to loading case III (the brightness in the figures indicates the relative magnitude of maximum shear stress).
ARTICLE IN PRESS W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252 1.2 Load
Displacement
40
1 35 0.8
30 25
0.6
20
0.4
15 10
0.2 5 0
Number of damaged elements
Displacement (mm) / load (*10 kN)
The following conclusions are drawn from our numerical simulations:
45 Damaged elements
0
0
10
20
30
40
50
60
70 80
90 100 110 120 130
Time (µs)
Stress wave II
(a)
60 Damaged elements
1.6
Load
Displacement
50 1.4 1.2
40
1 30 0.8 0.6
20
0.4 10 0.2 0 0
10
20
30
40
50
60
70 80
Number of damaged elements
Displacement (mm) / load (*10 kN)
1.8
0 90 100 110 120 130
Time (µs)
(b)
251
Stress wave III
Fig. 14. The variations of responded load, displacement at the top loading plate and the number of damaged elements during the failure process of the rock disk.
(1) Under static loading, numerical simulation reproduces well the primary crack at the disk center propagating along the diameter of specimen and the secondary cracks near the loading platens. (2) Under dynamic loading, different failure patterns and failure mechanisms are numerically captured when incident compressive stress waves with different amplitudes are applied to the rock specimen, which indicates that the failure pattern and failure mechanism of rock under dynamic loading is closely related to the propagation of stress waves. (3) Compared to the static loading condition, more cracks are developed in rock under dynamic loading, and the numerical results of this study are advantageous in capturing the failure pattern and failure mechanisms resulting from the propagation of stress wave in the rock specimen. In the future, care should be exercised in such simulations to accurately model the failure specimen subject to a variety of input stress pulses and to compare the results with experimental results. Although the work reported here has calculated on rock disks, such RFPA simulations are possible for a wide variety of specimen geometries.
Acknowledgements increased when rock is subjected to the dynamic loading. This is because much more damage is activated when the stress wave propagates in the rock. The cracks initiate, not only near the vertical diameter, but also in other parts of the rock disk. Under dynamic loading conditions, the rock specimen is much more fractured with the increase in magnitude of the stress waves. In contrast to the experimental results of the SHPB experiments, which were observed similar to a statically loaded specimen [31,40,41], the numerical simulation of this study is therefore advantageous in capturing the failure process and full failure mechanism during the propagation of a stress wave in the heterogeneous rock.
5. Summary and conclusions A numerical simulator, which considers the heterogeneity of material properties and the effect of strain rate on strength, was used to study the failure processes of a rock disk subjected to static and dynamic loading.
The authors would like to thank the financial support from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. PolyU 5044/99E) from 2001 to December 2002. They also thank the Natural Science Foundation of China (Grant Nos. 50504005 and 50204003) for their continuing support.
References [1] Lindholm US, Yeakley LM, Nagy A. The dynamic strength and fracture properties of Dreser basalt. Int J Rock Mech Min Sci Geomech Abstr 1974;11:181–91. [2] Shockey DA, Curran DR, Seaman L, Rosenberg JT, Petersen CF. Fragmentation of rock under dynamic loads. Int J Rock Mech Min Sci Geomech Abstr 1974;11:303–17. [3] Hanson ME. Some effects of macroscopic flaws on dynamic fracture patterns near a pressure-driven fracture. Int J Rock Mech Min Sci Geomech Abstr 1975;12:311–23. [4] Peng SS. A note on the fracture propagation and time-dependent behavior of rocks in uniaxial tension. Int J Rock Mech Min Sci 1975;12:125–7. [5] Janach W. The role of bulking in brittle failure of rock under rapid compression. Int J Rock Mech Min Sci 1976;13:177–86.
ARTICLE IN PRESS 252
W.C. Zhu, C.A. Tang / International Journal of Rock Mechanics & Mining Sciences 43 (2006) 236–252
[6] Chong KP, Boresi AP. Strain rate dependent mechanical properties of new Albany reference shale. Int J Rock Mech Min Sci Geomech Abstr 1990;27:199–205. [7] Ahrens TJ, Rubin AM. Impact-induced tensional failure in rock. J Geophys Res 1993;98(E1):1185–203. [8] He HL, Ahrens TJ. Mechanical properties of shock-damaged rocks. Int J Rock Mech Min Sci Geomech Abstr 1994;31(5): 525–33. [9] Zhao J, Li HB, Wu MB, Li TJ. Dynamic uniaxial compression test on granite. Int J Rock Mech Min Sci 1999;36(2):273–7. [10] Zhang ZX, Kou SQ, Yu J, Yu Y, Jiang LG, Lindqvist PA. Effects of loading rate on rock fracture. Int J Rock Mech Min Sci 1999;36(5):597–611. [11] Frew DJ, Forrestal MJ, Chen W. A split Hopkinson pressure bar technique to determine compressive stress-strain data for rock materials. Exp Mech 2001;41(1):40–6. [12] Grady DE, Kipp ME. The micromechanics of impact fracture of rock. Int J Rock Mech Min Sci 1979;16:293–302. [13] Zhao J. Application of Mohr–Coulomb and Hoek–Brown strength criteria to the dynamic strength of brittle rock. Int J Rock Mech Min Sci 2000;37:105–12. [14] Taylor LM, Chen EP, Kuszmaul JS. Microcrack-induced damage accumulation in brittle rock under dynamic loading. Comput Method Appl Mech Eng 1986;55:301–20. [15] Fahrenthold EP. A continuum damage model for fracture of brittle solids under dynamic loading. J Appl Mech ASME 1991;58:904–9. [16] Yang R, Bawden WF, Katsabanis PD. A new constitutive model for blast damage. Int J Rock Mech Min Sci Geomech Abstr 1996;33(3):245–54. [17] Yazdchi M, Valliappan S, Zhang WH. Continuum model for dynamic damage evaluation of anisotropic brittle materials. Int J Numer Methods Eng 1996;39(9):1555–83. [18] Liu LQ, Katsabanis PD. Development of a continuum damage model for blasting analysis. Int J Rock Mech Min Sci 1997;34(2): 217–31. [19] Li HB, Zhao J, Li TJ, Yuan JX. Analytical simulation of the dynamic compressive strength of granite using the sliding crack model. Int J Numer Anal Methods Geomech 2001;25:853–69. [20] Huang CY, Subhash G, Vitton SJ. A dynamic damage growth model for uniaxial compressive response of rock aggregates. Mech Mater 2002;34:267–77. [21] Homand-Etienne F, Hoxha D, Shao JF. A continuum damage constitutive law of brittle rocks. Comput Geotechn 1998;22(2): 135–51. [22] Brara A, Camborde F, Klepaczko JR, Mariotti C. Experimental and numerical study of concrete at high strain rates in tension. Mech Mater 2001;33:33–45. [23] Sanchidrian JA, Pesquero JM, Garbayo E. Damage in rock under explosive loading: implementation in DYNA2D of a TCK model. Int J Surf Min Reclamat 1992;6(3):109–14.
[24] Tang CA. Numerical simulation of progressive rock failure and associated seismicity. Int J Rock Mech Min Sci 1997;34:249–61. [25] Tang CA, Liu H, Lee PKK, Tsui Y, Tham LG. Numerical tests on micro-macro relationship of rock failure under uniaxial compression, part I: effect of heterogeneity. Int J Rock Mech Min Sci 2000;37:555–69. [26] Zhu WC, Tang CA. Numerical simulation on shear fracture process of concrete using mesoscopic mechanical model. Construct Build Mater 2002;16(8):453–63. [27] Zhu WC, Tang CA. Micromechanical model for simulating the fracture process of rock. Rock Mech Rock Eng 2004;37(1):25–56. [28] Chau KT, Zhu WC, Tang CA. Numerical simulation of failure of brittle solids under dynamic impact using a new computer program—DIFAR. Key Eng Mater 2004;261–263:239–44. [29] Zhu WC, Tang CA, Huang ZP, Liu J. A numerical study of the effect of loading conditions on the dynamic failure of rock. Int J Rock Mech Min Sci 2004;41(3):424. [30] Weibull W. A statistical distribution function of wide applicability. J Appl Mech 1951;18:293–7. [31] Tedesco JW, Ross CA, Mcgill PB, O’Neil BP. Numerical analysis of high strain rate concrete direct tension tests. Comput Struct 1991;40(2):313–27. [32] Rinehart JS, Pearson J. Behavior of metals under impulsive loads. Scranton, Pennsylvania: The Haddon Craftsmen, Inc.; 1954. [33] Williams RE, Potter RM, Miska S. Experiments in thermal spallation of various rocks. J Energy Resour Technol ASME 1996;118(1):2–8. [34] Cho SH, Ogata Y, Kaneko K. Strain-rate dependency of the dynamic tensile strength of rock. Int J Rock Mech Min Sci 2003;40:763–77. [35] Dı´ az-Rubio FG, Pe´rez JR, Ga´lvez VS. The spalling of long bars as a reliable method of measuring the dynamic tensile strength of ceramics. Int J Impact Eng 2002;27:161–77. [36] Rocco C, Guinea GV, Planas J, Elices M. Mechanisms of rupture in splitting tests. ACI Mater J 1999;96(1):52–60. [37] Hondros G. The evaluation of Poisson’s ratio and the modulus of materials of a low tensile resistance by the Brazilian (indirect tensile) test with particular reference to concrete. Aust J Appl Sci 1959;10:243–68. [38] Yanagidani T, Sano O, Terada M, Ito I. The observation of cracks propagating in diametrically-compressed rock disks. Int J Rock Mech Min Sci Geomech Abstr 1978;15:225–35. [39] Hudson JA, Brown ET, Rummel F. The controlled failure of rock disks and rings loaded in diametral compression. Int J Rock Mech Min Sci 1972;9(2):241–8. [40] Hughes ML, Tedesco JW, Ross CA. Numerical analysis of high strain rate splitting-tensile tests. Comput Struct 1993;47(4/5): 653–71. [41] Gomez JT, Shukla A, Sharma A. Static and dynamic behavior of concrete and granite in tension with damage. Theoret Appl Fract Mech 2001;36:37–49.