Microcrack-induced damage accumulation in brittle rock under dynamic loading

Microcrack-induced damage accumulation in brittle rock under dynamic loading

COMPUTER METHODS NORTH-HOLLAND IN APPLIED MECHANICS AND ENGINEERING 55 (1986) 301-320 MICROCRACK-INDUCED DAMAGE ACCUMULATION IN BRITTLE ROCK UNDE...

1MB Sizes 0 Downloads 108 Views

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

55 (1986) 301-320

MICROCRACK-INDUCED DAMAGE ACCUMULATION IN BRITTLE ROCK UNDER DYNAMIC LOADING*

Lee M. TAYLOR

and Er-Ping CHEN

Applied Mechanics III, Divkion 1523, Sandia National Laboratories, Albuquerque, NM 87185, U.S. A

Joel S. KUSZMAUL In Situ Technologies, Division 6258, Sand& National Laboratories, Albuquerque, NM 87185, U.S.A.

Received 6 March 1985 Revised manuscript received 12 August 1985

This paper is concerned with the development of a computational constitutive model which simulates the dynamic fracture behavior of brittle rock. The essential feature of this model is the treatment of the dynamic fracture process in rock as a continuous accrual of damage, where the damage mechanism has been attributed to the microcracking in the rock medium. Detailed descriptions of the model and its numerical implementation are given. A method is presented for determining the three additional material constants governing the damage accumulation. All of the numerical algorithms are amenable to vectorization for a vector processing computer and the model. is computationally very efficient. A numerical example of an oil-shale blasting experiment is presented and compared to field observations. Results demonstrate that the model is applicable to the prediction of dynamic rock-fracture behavior.

1. Introduction

The phenomenon of fracture in rock under rapidly applied loads is a topic of great current interest and importance because of its potential utilization in energy exploration and defenserelated applications. For example, in situ oil-shale retorting requires blasting of the rock in a of the rock play a confined region. The explosive loading and resulting fragmentation predominate role in creating a permeable broken rubble bed for subsequent retorting. Under dynamic loads, rock exhibits an inelastic brittle material response. This inelastic response is due principally to the stress-induced microcracks. The growth of these microcracks renders portions of the rock volume unable to carry load. This is then reflected in the decrease of the stiffness of the material. Consequently, the dynamic fracture process in rock may be modeled as a continuous accrual of damage, where the damage is defined to be the volume fraction of material that has been stress-relieved by multiple crack growth. This is the essential feature for the dynamic rock-fracture model developed in the present investigation. *This work performed at Sandia National Laboratories Contract Number DE-AC04-76D00789. 0045-7825/86/$3.50

@ 1986, Elsevier Science Publishers

is supported by the U.S. Department

B.V. (North-Holland)

of Energy under

302

L.M. Taylor et al., Microcrack-induced damage accumulation

The approach to damage modeling presented here falls under the theory of continuum damage mechanics. A comprehensive review of the subject was given by Krajcinovic [l] and will not be repeated here. Continuum damage models introduce internal state variables which represent, on a macroscopic scale, the distribution of microcracks within the material. The damage model described here is written in terms of a single scalar damage parameter, combining the theory of fracture mechanics governing the behavior of individual cracks with a statistical treatment to account for the random distribution of the microcracks. More elaborate damage models can be written using vector or tensor state variables which are capable of providing more refined data concerning the orientation, interaction, and geometry of the microcrack systems, but requiring a heavy computational penalty. For rock materials, where microcracking is identified as the damage mechanism under dynamic loading conditions, two models [2,3] exist in the literature. Both models combine the theory of fracture mechanics which governs the behavior of isolated cracks with some statistical treatment to account for the random distribution of the microcracks. The model presented here combines the effective bulk modulus expression given by Budiansky and O’Connell [4] with the fragment size equation given by Grady [5]. As such, it is not fundamentally different from the existing model by Kipp and Grady [3]. However, the use of two independent theories reduces the dependence of the model parameters on a single set of material data. This extends the range of applicability of the current model. Since dynamic material data are expensive and difficult to obtain, this represents a major improvement in the ability to model dynamic fracture of rock materials. Detailed equations for the damage model have been derived and are given in the text. These equations have been implemented into the finite element code DYNA2D [6]. All of the numerical algorithms are vectorized for vector processing computers and the model is computationally efficient. An application of the model to analyze an oil-shale blasting experiment is presented. For the range of geometry and loading conditions of the experiment, the numerical predictions compare favorably with field observations. There is presently very little published on the application of continuum damage models. This is combined with the lack of experimental data available to compare to numerical predictions. The main purpose of this work was to address this issue by presenting a complete description of the damage constitutive model with its numerical implementation and to compare numerical predictions to field data.

2. Constitutive

equations and assumptions

Predicting the behavior of rock under dynamic loading requires a constitutive model which is capable of treating the brittle failure of the material under tensile loading as well as the crushing type of failure associated with compressive loading of the material. The constitutive model presented here utilizes a damage growth model to predict the rock response under volumetric tensile loading. Under a volumetric compressive loading, the material is treated within the context of J2 flow theory as an elastic/perfectly plastic material. The basic thrust of this work was focused on the modeling of brittle failure, hence, a simple plasticity model was used in compression. A more complex model for compressive behavior would be desirable and will be considered in future work.

L.M. Taylor et al., Microcrack-induced

damage accumulation

303

2.1. Damage model The response of a brittle material to tensile loading is strongly affected by the strain rate. At very low strain rates a few large fractures dominate the fragmentation process and only a few large fragments are formed. At the other extreme, under high strain rates numerous small fractures simultaneously grow and interact, resulting in many small fragments. Also, the strain rate has a strong effect on the failure stress. A change of three decades in the strain rate, from 10 to 103, can result in a change in the failure stress by approximately one order of magnitude. The model described here explicitly includes these strain rate effects. The fundamental assumption of the damage model is that the material is permeated by an array of randomly distributed microcracks which grow and interact with one another under tensile loading. Hence, the model does not try to treat each individual crack, rather it treats the growth of the microcracks as an internal state variable which determines the accumulation of damage in the material. This damage is assumed to degrade the material stiffness of the rock and results in a strain softening behavior. Budiansky and O’Connell [4] used the self-consistent method to determine the elastic moduli of a cracked solid. Their results showed that for a random array of penny-shaped cracks in an isotropic elastic medium the effective bulk modulus of the cracked solid is

(1) where K and v are the original elastic constants for the intact or undamaged material, and the barred quantities represent the corresponding degraded constants of the damaged material. The crack-density parameter, Cd, is given by 45 (Y - V)(2 - F) cd = 16 (1 - F2)[10Y- F(l+ 3V)] If the crack density Cd is known, (2) can be used to determine the degraded Poisson’s ratio V. Kipp and Grady [3] used a Weibull distribution to determine the number of flaws per unit volume N that are active at a given pressure level P according to N = k(PJ3K)” ,

(3)

where k and m are material constants which must be determined as described below. It is important to note that the major advantage of this model over other similar damage models is the relative ease in determining these constants for different materials. Based on kinetic energy considerations, Grady [5] derived an expression for the nominal fragment diameter for dynamic fragmentation in a brittle material as d =

(;=-I)“,

In (4) p is the mass density,

C is the uniaxial wavespeed

(d/E/p),

and KIc is the fracture

304

L.M. Taylor et al., Microcrack-induced damage accumulation

toughness of the material. Also, d,, is the volumetric strain rate experienced by the material at fracture. This expression has been shown to provide a reasonable description of the fragment size data in oil shale, as well as for other brittle materials [5]. Since the fragmentation process can be considered as the interaction of active cracks under applied loads, a correlation between the nominal fragment size and the characteristic flaw dimension in the material is expected. A reasonable approach, which is followed here, is to assume that the characteristic flaw size, a, is directly proportional to the radius of the nominal fragment (i.e. a -id). Now, the crack-density parameter represents the volume fraction of material made up of voids. While the theory is not capable of providing the crack opening displacements, it is clear that the characteristic volume of a void is proportional to the characteristic length cubed. Hence, the crack-density parameter is given by Cd =

N+la3)

(5)

)

where p is the unknown proportionality constant. Combining (3)-(5) and absorbing the proportionality constants for the fragment size and volume into the constant k gives an expression for the crack density, cd=--

5

k

2(3K)”

&

-

( pc >

2pmd_2 Inax*

(6)

Now, (2) and (6) represent two independent expressions for C,. Using (6) to evaluate the crack density in terms of the known pressure and strain rate history, (2) can be solved for the current degraded Poisson’s ratio V. This value, along with the crack density, can then be used in (1) to determine the degraded bulk modulus. At this point it is convenient to define a single damage parameter D, which varies from zero, initially, to a maximum of one. The damage is defined as the volume fraction of material that has lost its load-carrying capability. Inspection of (1) suggests that the damage be defined bY D = gfl(Y)C,, (7) where the function fi(fi) is defined as fl(fi)

=

(1- F’) (1- 2Y) *

(8)

There is a one-to-one correspondence between the damage and the crack density. For a given crack density and original Poisson’s ratio, (2) represents a cubic equation for the degraded or damaged Poisson’s ratio V. For a crack density of zero (undamaged material), (2) yields a value of the original Poisson’s ratio for 5. As C, +&, Y+ 0, and the damage, as defined in (7) approaches one. Hence, the effective bulk modulus given by (1) vanishes as the crack density approaches a critical value of &. It seems plausible that, as the microcracks grow and intersect to the point that slightly more than one-half of the material volume is composed of cracks, the material stiffness is completely lost.

L.M. Taylor et al., Microcrack-induced

damage accumulation

305

The usual finear elastic strain relations for stress can be written in terms of the volumetric and deviatoric parts as P=3Ks,,

(9)

SQ= 2Geij,

(10)

where e, is the deviatoric part of the strain tensor, E, is the mean volumetric strain, S, is the deviatoric part of the stress tensor, P is the mean or bulk pressure, and G is the shear modulus of the material. Using the definition of damage given by (7), the stress-strain relations of (9) and (10) can be written in terms of the effective moduli and the single damage parameter, P = 3K(l-

D)E, ,

(11)

Sij = 2G(l-

LI)e#.

(12)

Here, the shear modulus is degraded in the same manner as the bulk modulus. 2.2. Elastic/perfectly plastic model The e~astic/pe~ectly plastic model used here in compression follows the des~~ption given by Krieg and Krieg [7]. It is described here for completeness. The total deviatoric strain rate, C;j,is assumed to consist of elastic and plastic parts, tz and kc’, respectively, such that &jz+?_!+&.

(13)

The deviatoric stress rate is assumed to be related to the elastic deviatoric strain rate by S..11=: 2Gi??819

(14)

where G is the original elastic shear modulus (not degraded by any tensile damage). The bulk response in compression is assumed to be purely elastic and is given by (9). The von Mises equivalent stress CFis defined in terms of the deviatoric stress as @=

V~SjiSji.

(15)

The yield stress, cr,, is assumed constant and the von Mises yield function is written as 9

=a-cr,.

3. Numerical impIementation The constitutive

(16) of the constitutive

model

model described here was implemented

into the DYNA2D [6] program, an

306

L.M. Taylor et al., Microcrack-induced damage accumulation

explicit finite element code for the transient, large-deformation analysis of two-dimensional solids (plane-strain and axisymmetric). The model has been vectorized for use on computers with vector processing, and the algorithms described here have been implemented in a fully vectorized form. In the following sections a complete description of the numerical algorithms required to implement the constitutive models is presented. The numerical finite element method has been discussed by Zienkiewicz [8] for solving solid mechanics problems and will not be described here. 3.1. Damage algorithm One of the most important requirements of any numerical algorithm is that it can be efficiently implemented. Unfortunately, for given values of C, and v (2) represents a cubic in V and the determination of 5 is a nontrivial numerical calculation to find the roots, of the equation. As a simplification (2) has been approximated with a linear analytic function for fi in terms of v and C,, v= V(l_$Cd).

(17)

Fig. 1 shows a plot of fi as a function of crack density for values of Poisson’s ratio between zero and a half. The solid line represents the results of using (2) to determine 6 and the dashed line shows the linearized result of (17). For realistic values of Poisson’s ratio for rocks (generally less than 0.3), the error in using the linearized equation is less than 5 percent; the maximum error for any value of Poisson’s ratio is 6.61 percent. Furthermore, the error in the

0.5

Ia p

0.4

2

a

_* g 0.3 3 B +

0.2

z a w 5 0.1 s 0.0 0.25

0.00 CRACK

DENSITY

, Cd

Fig. 1. Effect of crack density on Poisson’s ratio.

L.M. Taylor et al., Microcrack-induced

damage accumulation

307

degraded bulk modulus as determined from (1) is even less. Fig. 2 shows the bulk modulus as a function of crack density. The finite element algorithms require the constitutive equations to be written in an incremental or rate form. Taking the rate of (7), (ll), and (12) gives B = F-fi(fi)C, + +Ldfi(fi) ) 3K(l-

O).& - 3&B,

Sij = 2G(l-

D)$ - 2G&‘,

i, =

(18) (19) (20)

where (21) From (17),

(22)

(23) then (18) can be written as D =

~lf~(fi)-$Vf2(~)cd]cd.

(24)

n K

0.25

0.25

0.00 CRACK

0.50

DENSITY

Fig. 2. Effective bulk modulus as a function of crack density.

L.M. Taylor et al., Microcrack-induced damage accumulation

308

Finally, taking the rate of (6) gives

c=5 d

km

( > &

.2 (3K)“_l

*pi

pc

i-2

” Inax*

(25)

Equations (19), (20), (24), and (25) represent seven coupled ordinary differential equations to be integrated. Integrating these rate equations while maintaining the high computational speed available on a vector processing computer was one of the primary goals of the numerical algorithm. A fourth-order Runge-Kutta procedure [9] was used to integrate the differential equations. This integration scheme was used by Taylor [lo] for integrating the Kipp and Grady [3] damage model in the DYNA2D code. Those interested in the details of coding such a scheme for a vector processing computer should see the FORTRAN listing in [lo, appendix]. The Runge-Kutta integration operator requires four function evaluations-one at the beginning of the time step, two at the center of the time step, and one at the end of the time step. The strain rates are assumed to be constant over the time step. Hence, the strains E, and &iiare known at any time within the time step. The seven rate equations can be written in vector notation as

where y is a seven-dimensional

vector,

y’ = [D, cd, P, sll, &2, s33&I .

(27)

If At is the time step size and the subscript 0 indicates quantities evaluated at the beginning of the time step, then the four function evaluations can be written as g,

=f(to, YoPt

g2

=f(t,

+

fAt,

~0

+ $g,)At7

(29)

g, =f(t,

+

IAt,

Y,

+ ig,)At 9

(30)

g,

=f(t, + At,

7

Y,

+ g,)At.

The increment in the dependent AY = t
(31)

variables is then given by

2g, + 2g3+ a) .

(32)

Since only function evaluations are required, this integration method is easily vectorizable on a computer with vector processing. The computer program integrates the seven rate constitutive equations for 128 integration points at a time. The constitutive routine first calculates the stresses at the integration points as if they were elastic/perfectly plastic, irrespective of whether or not the material point is in tension or compression. This is done using the elastic-plastic algorithm described in the next section. The routine then integrates the damage

L.M. Taylor et al., Microcrack-induced

damage accumulation

309

equations as if the material points were in tension. Finally, it picks the correct stress state from one of the two calculations, based on the sign of the volumetric strain. This may seem inefficient, but, on a vector processing computer, this allows the calculations to be performed in a vector operation. There are no logical tests or subroutine calls, and all 128 integration points are processed in one vector operation. The default time step in DYNA2D is calculated internally by the code and is a function of the smallest element diagonal. The equations presented here for rate of damage accumulation represent a stiff system of equations and cannot be integrated accurately using this internally computed time step. Experience with the model has shown that if the default time step is reduced by a factor of six, the calculation always remains stable. Consequently, all the results shown below use this smaller time step size. 3.2. Elastic-plastic algorithm The elastic/perfectly plastic algorithm used here for the response of the material in compression is the usual elastic predictor/radial return algorithm found in most nonlinear finite difference or finite element codes. The volumetric response is totally elastic and calculates the increment in pressure as AP = 3K&At.

(33)

All of the plasticity takes place through the deviatoric response. An elastic trial stress is first calculated, S; = S; + 2Gi,At,

(34)

then the algorithm checks the yield function (16). If negative, the material is elastic and the stress calculation is finished. If positive, the trial stress must be modified to bring the state of stress back to the yield surface. This is accomplished by a radial return to the yield surface, with the final deviatoric stress given by (35) This approach to the integration accurate for most problems [7]. 4. ~termin~tion

of the equations

of plasticity has been shown to be quite

of the damage material constants

As discussed above, evaluation. of the material constants for the damage model should be readily accomplished with laboratory data that can be obtained relatively easily, One straightforward approach for determining the damage parameters k and m follows by assuming that Poisson’s ratio is zero. The effect of this assumption appears only in the value of k. Fortunately, the model is not particularly sensitive to variations in k. Changing k by 25 percent only affects the peak stresses in the example problems in the next section by 5 percent. Defining E as P/3K and ignoring Poisson’s ratio, the damage under constant strain rate loading

310

L.M. Taylor et al., Microcrack-induced

damage accumulation

can be written as D = AkL”i-‘,

(36)

where (37) At a constant strain rate, the mean stress is related to the average strain by P(t) = 3K(l-

D)&(t),

(38)

P(t) = 3K(l-

AkC’-2t”)it.

(39)

or

The time at which the stress is a maximum in (39) is given by

. 1-l=[(m+l)Ak]-‘/“$-“‘1” Using (40) in (39) gives P(t,,,) = 3K(m + l)-(m+1)‘m(Ak)-1’m~2’m.

(41)

The logarithm of (41) provides InP(t,)=

C,+&ln m

6.

(42)

Two points on the fracture stress versus strain rate curve can be used with (42) to determine m. Using (41), an expression for k can be determined by 3Km

m+l

PC*,>

m

In---

ln(m + 1) - In(AE-2). I

(43)

If laboratory data for fracture stress versus strain rate are not available, it is possible to generate this data using an expression derived by Kipp et al. [ll],

where N, is a shape factor (1.12 for penny-shaped cracks) and C, is the shear wave velocity of the material. Equation (44) has been shown to be a reasonable approximation for a number of rock types.

L.M. Taylor et al., Microcrack-induced

damage accumulation

311

5. Numerical simulations The first simulation on the model was on a simple one-dimensional, constant volumetric strain rate test. The material used is a nominal 30 gallon-per-ton (GPT) oil shale with a Young’s modulus of 10.8 GPa, Poisson’s ratio of 0.2, density of 2270 kg/m3, and a yield stress in compression of 200 MPa. Kipp and Grady [3] reported fracture stress versus strain rate data for this material, and this shown in Fig. 3. Using these data in (42) and (43) results in values of 7.0 for m and 5.116 x 10Z l/m3 for k. Fig. 4 shows the pressure or mean stress and damage as a function of time for a strain rate of lOOO/sec. Fig. 5 shows the pressure versus time, and Fig. 6 shows pressure versus volumetric strain for three decades of strain rate. In the formulation of the numerical model, when the material unloads volumetrically, it does so along the original (undamaged) bulk modulus. During such an unloading event the damage remains constant (damage is irreversible). Fig. 7 shows pressure versus strain during such an unloading excursion. The nominal strain rate in this calculation was lOOO/sec., and the dashed line in the figure shows the pressure if no unloading occurred. The material reloads along the damaged bulk modulus. Fig. 8 shows damage versus strain during the unloading process. Kipp and Grady [3] used a value of damage of 0.2 to identify the extent of cratering in the blasting experiments they analyzed. This value was chosen because the 0.2 damage contour provided a reasonable match to the crater which was excavated in the field. It is interesting to note that in the one-dimensional calculations the peak stress before strain softening generally occurs at a damage level of approximately 0.2. Next, the model was applied to a blasting calculation to simulate the SB-1 blasting experiments at the Anvil Points Mine. In this cratering test a single column of high explosive (HE) was bottom detonated. The HE was Iregel 1175U. The HE column was 2.5 m in height, 0.081 m in radius, and was centered at 3.75m below the ground surface. The axisymmetric mesh used in the analysis is shown in Fig. 9 along with the layering of the different oil-shale grades in gallons per ton (GPT). Table 1 shows the material property data FRACTURE STRESS DEPENDENCE ON STRAIN RATE 1

100

1

OIL SHALE

0 GAS GUN 0 CAPACITOR a HOPKINSON EAR 0 TENSILE EAR

10

I

100

I 10’

I

I

103 102 STRAIN RATE (s-l)

I 104

Fig. 3. Fracture stress dependence on strain rate.

L.M. Taylor et al., Microcrack-induced damage accumulation

312

- 0.9 60

50 3 ii -

40

ii 2 2

-

0.6

-

0.5

-

0.4

-

0.3

-

0.2

-

0.1

8 r” d

30

E

20 k=

x 102*l/m3

5.116

10

0 0

2

4

6 TIME

Fig. 4. One-dimensional 70

I

10

6

12

14

(psec)

constant

strain rate pressure

I

I

and damage

1

I

I

P= 2270 m= 7

i,=103

K9/m3

22 k=5.115 x 10 l/m3 KIC= 1 MPa Viii

50

\

z

I

I

K= 6 GPa v= 0.2

60

histories.

% 40 z .z v) w 30

1f=10*

E

ti

0 0

20

40

60 TIME

Fig. 5. One-dimensional

constant

60

100

120

140

(psec)

strain rate pressure

histories

for three strain rates.

L.M. Taylor et al,, Microcrack-induced

313

damage accumulation

70

60

50 3 c

40

2 2 v) g n

30 K=6

GPa

u=o.2

20

f’=2270

Kg/m3

m=7 k=5.115

10

0

x lo**

K~c=l

MPa diE

0.005

0.006

L

0.000

0.00 1

0.002

0.003

0.004

I l/m3

0.007

STRAIN

Fig. 6. One-dimensional

constant I

I

strain rate pressure I

strain curves for three strain rates.

I

I

I

I

I

60

P= 2270 K=6

kg/m3

GPa

u=o.2 k= 5.116

x lo**

l/m3

m=7 Krc= 1 MPa \ms

0.000

0.002

0.004

0.006

0.006

0.010

0.012

0.014

0.016

0.016

STRAIN

Fig. 7. One-dimensional

pressure

strain curve during a loading,

unloading,

and reloading

cycle.

314

L.M. Taylor et al., Microcrack-induced I

I

I

damage accumulation

I

1

I

I

I

0.6

0.7

0.6 Y 3

0.5

2 0.4

0.3

L’=

2270

K= Y=

6GPa 0.2

k=

5.116

m=

7

kr/m3

x 1022Pa

KIC = 1 MPa X% 0.2

0.1

0.0 0.000

0.002

0.004

0.006

0.000

0.010

0.012

0.014

0.016

0.016

STRAIN

Fig. 8. One-dimensional

damage

strain curve during a loading,

unloading,

and reloading SHALE

GRADE

(OPT)

1

cycle.

24.0 -21.2 -11.012.2 21.2 ‘1422.0 *15.2

l ,*

17.2

‘loI11.4 .‘la.2 8.7 -10.1 z11.2 s!!

32.6

2s.s 22.0 G21.2 t-32.2

-6.0

24.0

32.2

21.2

0.0 Fig. 9. Finite element

1.5

3.0

4.5

grid and stratigraphy

6.0

7.5

9.0

10.5

12.0

used in the analysis of the SB-1 oil-shale

blasting

experiment.

L.M. Taylor et al., Microcrack-induced

damage accumulation

315

Table 1 Material properties of different oil-shale grades used in the analysis of the SB-I oil-shale blasting experiment Grade

P

m

V

t&cc)

W’V 8.7 10.1 11.4 13.5 15.8 17.3 21.2 24.0 25.5 29.0

2.523 2.493 2.465 2.420 2.370 2.338 2.261 2.216 2.192 2.136

30.23 28.82 27.51 25.39 23.07 21.56 17.83 15.46 14.20 11.24

0.259 0.260 0.261 0.264 0.266 0.267 0.271 0.274 0.276 0.279

102 103 103 104 106 107 106 102 100 95

38.6

2.008

7.03

0.289

78

51.1

1.870

4.34

0.301

70

4.34 X 1oz7 2.%X10Z 2.04 x lp 1.08 x 1e7 5.01 x 10% 2.92x 10% 6.44X1025 2.09 x lO= 1.06x 1025 1.61x 10% 3.70x lo= 7.26x 10”’

7 7 7 7 7 7 7 7 7 7

5.12 5.20 5.27 5.38 5.51 5.59 5.80 5.95 6.05 6.38

x x x x x x x x x x

16 Id lo5 los los 16 16 ld 10’

7

10’ 7.26 x 10’

7

8.96 x lo5

used for the 12 grades of oil shale present. The HE was modeled using the JWL equation of state [12] in DYNA2D with density detonation

1.25 gmfcc, velocity

Chapman-Jouguet

A B RI R2 0 & V,

pressure

6178 m/see, 11.4 GPa, 47.6 GPa, 0.529 GPa, 3.5, 0.9, 1.3, 4.5 GPa, 1.0.

The experiment was analyzed twice. First, the rock was assumed to be isotropic with material properties corresponding to the 21.2 GPT oil shale. The volume of rock having a damage of at least 0.2 is shown in Fig. 10. Next, the analysis was repeated using the layered stratigraphy. Fig. 11 shows the resulting damage and the effect the layering has on the volume of rock that was damaged to the 0.2 level. Although the radial extent of the two damaged regions appears to differ only slightly, the volume of damaged rock is much smaller in the layered case. Since it is the volume of rock that can be fragmented and heaved during a blasting event that is of interest, the effect of the stratigraphy on the analysis appears important. The total CPU time required was 600 seconds for each case.

316

L.M. Taylor et al., Microcrack-induced

damage accumulation ,

SURFACE

ND SE-2

CRATER

, PROFILES

ISOTROPIC ROCK (Axisymmetric Case)

0.0

Extent

of damage

exceeding

EXPLOSIVE

m

FRACTURED

5.0

2.5

RADIAL Fig. 10.

m

0.2 assuming

DISTANCE

a homogeneous

7.5

ROCK

10.0

12.5

(meters) and isotropic

rock mass in the SB-1 simulation.

A stress gauge was grouted in place at a radial distance of 3.0 m and a depth of 2.5 m. Fig. 12 shows the measured radial stress at that point along with the numerical predictions for the two different cases (layered and isotropic). The arrival time and peak pressures are predicted in both cases quite accurately. In Fig. 13 the boundary of the 0.2 damaged region is shown again for the layered analysis. The average profile of the excavated crater is shown. This excavation was performed by using a backhoe and removing all ‘loose’ material. It represents a somewhat subjective measure of cratering since it is dependent upon what the backhoe operator considered loose material. However, it tends to lend credence to the numerical results. Note that the damage model predicts a zone of undamaged rock centered approximately 2.5 m below the surface. If the backhoe operator encountered a similar region during his excavations, he would have concluded that he had reached the bottom of the crater. The numerical model predicts a significant quantity of fragmented and loose rock under this undamaged zone. Also, viewed from inside the crater, fingers of loose fragmented rock were observed extending back along the layering of the oil shale as is seen in the numerical predictions. A core was taken along line AA shown in Fig. 13. At the point where the core was taken, the actual excavated crater did not extend quite as far radially. Hence, the core was begun on a point outside the crater, and line AA shows only the direction of coring. Fig. 14 shows the

L.M. Taylor et al., Microcrack-induced

317

damage accumulation

.;:. _ ____ :.. ..:.* ::.. Pan.‘... .‘..

b

t4J 7.5 n

LAYERED ROCK (Axisymmetric Case)

10.0

12.5 0.0

of damage

EXPLOSIVE

m

FRACTURED

ROCK

I

I

I

1

2.5

5.0

7.5

10.0

RADIAL Fig. 11. Extent

m

exceeding

DISTANCE

12.5

(meters)

0.2 using the site stratigraphy

in the SB-1 simulation.

1.2 1.0 0.8 0.6

- -.-

0.4

-

SE-1 MEASUREMENT DYNAPD LAYERED MODEL DYNA2D ISOTROPIC MODEL

0.2 0.0 -0.2 0.0

Fig. 12. Comparison

0.5

1.0

1.5

2.0

ELAPSED

TIME

of numerical

predictions

2.5

3.0

(milliseconds) of radial stress with experimental

data at r = 3.0 m, z = 2.5 m.

318

L.M. Taylor et al., Microcrack-induced damage accumulation b

LAYERED ROCK (Axisymmctric Case)

10.0

12.5

0.0

m

EXPLOSIVE

m

FRACTURED

ROCK

I

I

I

I

2.5

5.0

7.5

lo.0

12.5

RADIAL DISTANCE (meters) Fig. 13. Crater profile and core location for the SB-1 oil-shale blasting experiment.

LEGEND 0 INTACT (Odk0.2) DFRACTURED (O.%D<0.4) aWELL FRACTURED (0.4
5.0

CORE

DYNAPD

6.0 Fig. 14. Comparison

of core sample with model predictions.

L.M. Taylor et ai., Microcrack-induced damage accumulation

319

core log along line AA. This core log was compiled by inspecting the core material and characterizing the damaged rock according to the legend shown in the figure. Also shown is a simulated core log, compiled by assigning numerical values of damage to the levels shown in the legend. These levels are shown in parentheses in the legend. The model prediction for the core is then just a one-dimensional representation of a contour plot. The model predictions correlate very well with the core logs, especially where the core log shows that there were, in fact, zones of intact and undamaged rock.

6. Conclusions

A damage model for materials exhibiting strain-rate-dependent fracture behavior has been developed in the present study. Damage is defined to be the volume fraction of material that has been stress relieved by multiple microcrack interaction and growth. Damage accumulation is reflected in the continuous degradation of the material moduli during the loading process. Using Budiansky and O’Connell’s effective moduli approach for a cracked solid, methods for determining the crack-density and characteristic crack-dimension parameters were developed. Because of the statistical nature of the distribution of microcracks, a Weibull distribution was assumed in the model for the crack-density parameter. Constants for the Weibull distribution representation can be determined from laboratory data characterizing fracture stress as a function of the strain rate. The characteristic crack dimension was expressed by the fragment size equation given by Grady [5] which takes into consideration the effect of kinetic energy in a fragmentation event. In this manner, the determination of the crack-density parameter is independent of that for the characteristic crack dimension. This reduces the amount of dynamic data required to determine the material constants required in the model. Judging from the scarcity of available dynamic data, this represents a major advantage in the application of the model. The constitutive model has been implemented into the DYNALZD f6] finite element code. Detailed derivations of the rate equations and their numerical algorithms for code implementation were presented. All of the coding was vectorized, and the model runs very efficiently on vector processing computers. The model has been applied to treat oil-shale blasting problems. Results indicate the model reproduces, qualitatively, the brittle behavior of rock under blasting conditions. The model was shown to apply over a wide range of strain rates. Moreover, experience with the numerical implementation has shown the algorithms to be extremely robust and easy to use. Hence, routine blasting calculations can be performed to aid the engineer in designing blastiiig operations. Future work on the model will include a more sophisticated model of soil behavior in compression probably in the form of a Drucker-Prager type model. Also, the model may be extended to three dimensions so that the interactions of multiple explosive charges can be investigated.

References [l] D. Krajcinovic,

Continuum

damage mechanics, Appl. Mech. Rev. 37(l) (1984) l-6.

320

L.M. Taylor et al., Microcrack-induced damage accumulation

[2] L.G. Margolin and T.F. Adams, Numerical simulation of fracture, in: Issues in Rock Mechanics, Proceedings 23rd Symposium on Rock Mechanics (Society of Mining Engineers, 1982) 637-644. [3] M.E. Kipp and D.E. Grady, Numerical studies of rock fragmentation, SAND79-1582, Sandia National Laboratories, Albuquerque, NM, 1978. [4] B. Budiansky and R.J. O’Connell, Elastic moduli of a cracked solid, Internat. J. Solid and Structures 12 (1976) 81-97. [5] D. Grady, The mechanics of fracture under high-rate stress loading, in: Z.P. Bazant, ed., Preprints of the William Prager Symposium on Mechanics of Geomaterials: Rocks, Concretes and Soils, Northwestern University, Evanston, IL, 1983. [6] J.O. Hallquist, Users manual for DYNA2D-an explicit two-dimensional hydrodynamic finite element code with interactive rezoning, UCID-18756, Rev 1, Lawrence Livermore Laboratory, Livermore, CA, 1982. [7] R.D. Krieg and D.B. Krieg, Accuracies of numerical solution methods for the elastic-perfectly plastic model, ASME J. Pressure Vessel Technol. 99 (1979) 510-515. [8] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, New York, third ed., 1977). [9] C.R. Wylie, Advanced Engineering Mathematics (McGraw-Hill, New York, fourth ed., 1975). [lo] L.M. Taylor, Implementation of a vectorized damage constitutive model in the DYNA2D finite element program, SAND83-2660, Sandia National Laboratories, Albuquerque, NM, 1984. [ll] M.E. Kipp, D.E. Grady and E.-P. Chen, Strain-rate dependent fracture initiation, Internat. J. Fracture 16 (5) (1980). [12] B.M. Dobratz, LLNL explosives handbook, properties of chemical explosives and explosive simulants, UCRL-52997, Lawrence Livermore Laboratory, Livermore, CA, 1981.