of Rock

Memoirs of the Faculty of Engineering,Okayama UniversitY,Vo1.28, No.1,pp.l29-138, November 1993

Numerical Simulation of Rock Toughness Testing

Takeo TANIGUCHI*, Sohichi HIROSE*, Finn OUCHTERLONY**, Kohji NAKAGAWA***, Akihiko MIYAJI**** and Yasufumi FUKUOKA*****

(Received September 30 , 1993)

SYNOPSIS

The testing method of rock toughness is proposed by the international society of rock mechanics (ISRM), but the results may be influenced by the test pieces, and the details of the crack propagation and the stress intensity factors are not clarified through the testing. Also the experimental test requires tedious works for the preparation of test specimen and economical responsibility. The present study aims to simulate numerically the rock toughness testing which is proposed by ISRM. For this purpose, the authors propose a numerical method which can simulate the experimental testing, and they show the propriety of the proposed method by comparing the results with the experimental and other numerical methods. At the same time, they clarify the details of crack propagation behaviors in rocks, and show the change of the stress intensity factors. The proposed method is based on the displacement-type finite element method, and several techniques are introduced to obtain accurate solution of the mechanical behavior near the crack-tip area.

1. INTRODUCTION

ISRM (International Society for Rock Mechanics) proposed a standard testing method of the fracture toughness of rocks, and the specimen for the testing is called Chevron bend specimen.

Using this specimen, fracture toughness of various rocks has been tested. The requirement of the numerical simulation is, however, increasing, since the results of the fracture toughness depend on the specimen, and the structural experiment requires tedious works to prepare the specimen and long period, and it also requires economical responsibility.

*

**

***

****

*****

Department of Engineering Science, Okayama University Swedish Blasting Laboratory, Sweden

Department of Civil Engineering, Yamaguchi University Japan Development and Construction Company Ltd.

Kawasaki Heavy Industry Company Ltd.

129

serious problems to be resolved for the development of the numerical simulation of the fracture toughness testing. The first problem is how to obtain accurate stress distribution which appears near the crack-tip area, The second problem is the geometrical modelling of the three dimensional specimen with complex boundary geometry. The last problem is that the anisotropic property must be taken into consideration at the analysis.

The finite element method has been used to solve problems mentioned above[l]. Using the displacement-type FEM, stable solutions can be expected to show the crack opening deformation, and the introduction of the solution into the expressions proposed by Ingraffea can directly lead to the estimation of the stress intensity factors[2]. For getting accurate stress distribution (or more precisely the crack opening displacement), there exist singular isoparametric elements which can sufficiently express the characteristics of the crack. Then, at the geometrical modelling of the specimen, we introduce the elements at the crack-tip domain.

The material used for the numerical analysis is rock, and, therefore, the anisotropy must be introduced to express the constitutive law of the material. In this investigation the transversely isotropic body is dealt instead of the general anisotropy, since the number of independent param- eters is too many in case of the general anisotropy.

Summarizing above items, the authors show the numerical simulation of the rock toughness testing of ISRM, which is based on the displacement-type finite element method, and the propriety of the method is investigated by comparing the result with the experiments and other numerical results,

2. ESTIMATION OF STRESS INTENSITY FACTORS

It is known that accurate solutions using the finite element method can be obtained if an appropriate type of element is prepared to express the stress distribution at the crack-tip and sufficiently fine meshes are set in the domain. For the former, we introduce a I5-node singular isoparametric element which is obtained from a 20-node isoparametric element by condensating one rectangular face of the hexahedron and by moving mid-nodes on edges adjacent to the condensated surface to quarter-points. As a result, the singular isoparametric element becomes a triangular prism, and this element is placed as to envelope the crack-tip as shown in Fig. 1. The residual domain is modelled using 20-node isoparametric elements.

The stress intensity factors can be estimated using the nodal displacements facing to the crack- tip. Incase of the isotropic body, the stress intensity factors, K l , KII andKIllfor Mode I, II and III are evaluated using the nodal displacements as follows[2],

(1)

Numerical Simulation of Rock Toughness Testing

z,w

Figure 1: Modelling at crack-tip area.

where matrices [B] and {A} are expressed as

131

[ ^~

[B] = 0

o

o

4(I-v²)

E

o

o ] o ,

4(1i

^{v ) (2)}

{A} =

2VB - Vc

+

^{2VE - VF}

+

^{VD - 2VB'}

+

^{VC' - 2VE'}

+

^{VF' - VD'}

+17(

^-4VB

+

^Vc

+

^{4VE - VF}

+

4VB' - VC' - 4VE'

+

^vF')/2

+1]2(VF

+

^{Vc -} 2VD - VF' - VC'

+

^2VD,)/2

2UB - Uc

+

^{2UE - UF}

+

^{UD - 2UB'}

+

^{uc' - 2UE'}

+

^{UF' - UD'}

+17(

^-4UB

+

^Uc

+

^{4UE - UF}

+

4UB' - UC' - 4UE'

+

^{UF' )/2}

+172

(UF

+

^{Uc - 2UD - Up - Uc'}

+

^{2UD' )/2 (3)}

2WB - Wc

+

^{2WE - WF}

+

^{WD - 2WB'}

+

^{WC' - 2WE'}

+

^{WF' - WD'}

+1]( -4WB

+

^Wc

+

^{4WE - WF}

+

^{4WB' - WC' - 4WE'}

+

wF')/2

+1]2(WF

+

Wc - 2WD - Wp - Wc'

+

2WD,)/2

In above equations, E and v show the Young's modulus and the Poisson's ratio, anduo" Va and

W a show the displacements at the node a(a= B, C, ... ,see Fig. 1 ) alongx, y and z axes of the local coordinates. 1]is equal to -(2z/L2

+

1), and L2 expresses the length of the element along the crack-front.

Now, we rewrite the above equations which can be applicable to anisotropic material with the stress-strain relationship expressed by fi = aijo"j or O"i = Cijfj (i,j = 1,2, ... ,6). Assume that all axes of the coordinates coincide with the axis of the anisotropic material. The stress intensity factors for anisotropic materials can, then, be evaluated using eq.(l) on the condition that eq.(2) is replaced by

(4)

5

Figure 2: Chevron bend specimen.

Here,⁵¹ and⁵² in eq.(4)are the solutions with positive imaginary parts fo~the following equation;

Two parameters, Pj and qj,are given by

an⁵2^j

+

^{a12 - a165j,}

a12^Sj

+

^{a2dSj - a26'}

3. NUMERICAL RESULTS FOR ISOTROPIC MATERIAL

(5)

(6) (7)

The geometry of the Chevron bend specimen and its testing method are illustrated in Fig. 2.

As obvious from the geometry of the specimen, the specimen has two symmetric surfaces. We may, therefore, treat only a quarter of the specimen as its numerical model. Fig. 3 shows the finite element model for a quarter of the Chevron bend specimen. Finer meshes are located near the domain of the crack-tip, and rather coarse meshes for the residual.

Five stages of crack propagation are considered for the estimation of the stress intensity factors, and the finite element meshes facing the crack surface are illustrated in Fig. 4. In these five cases, we assume the crack-front B (see Fig. 2) always lies along the 3-axis, and locations of the crack-front are given as follows;

Case 1 : a - ^ao= 0.05D Case 3 :a - ao= 0.15D Case5 : a - ao

=

0.25D

Case2 : a - ao= O.lOD Case 4: a - ao=0.20D

Following numerals are used for determining the Chevron bend specimen.

D

=

100mm, S

=

333mm, a

=

l5mm, 2(}

=

90^0, t

=

Omm

The results of the stress intensity factors are shown in Fig. 4, in which the stress intensity factors are evaluated using eq.(l). Note that Ain the figure show the non-dimensional stress intensity

Numerical Simulation of Rock Toughness Testing

L1= 1.5 L2= 1. 5 L3= I.5 L4= 4.5 L5= 9. 0 L6=IH.5

L6 (mm)

Figure 3: Finite element model for a quarter of the specimen.

Table 1: Average stress intensity factors

A.

Case 1 Case 2 Case 3 Case 4 Case 5 Case 5' a - ao 0.05D 0.10D 0.15D 0.20D 0.25D 0.25D

<1 > Ouchterlony[l] 10.548 10.031 9.889 10.305 10.935 ^-

<2 > With singular elements 10.191 9.691 9.640 9.956 10.433 10.436

<3> W

10

singular elements 9.632 9.172 9.076 9.410 9.839 -

133

factor defined by A= Kr/[F

1

D1.^5]. Symbols in Fig. 4 are defined as followings.

<

1 >

0:

The results by Ouchterlony[l]

<2 > . : The results of FE analysis with the singular isoparametric elements

<3 > .A.: The results of FE analysis without the singular isoparametric elements

Due to the geometry of the crack front, the distribution of the stress intensity factor along the crack front shows the highest value at the corner where the crack front bends. The analysis using the 15-node singular isoparametric elements shows good coincidence with the result by Ouchterlony comparing to the value obtained by the FE analysis without the singular elements. Table 1 shows the average value

A

of the non-dimensional stress intensity factors along the crack front B. The rows, < 1 >, < 2 > and<3 >, show the average values of the symbols,

0, •

and.A. in Fig. 4. In order to find the influence of the meshes, a different mesh pattern is prepared as shown in Fig. 4(f) and used for the analysis. The result is shown in the column "Case 5' " of Table 1. Comparing the results of Case 5 and Case 5', no difference between them is found. Hence we may conclude that the models using the 15-node singular isoparametric elements as shown in (a) through (e) of Fig. 4 are valid for our numerical investigations. Fig. 5 shows the comparison of our result with the ones by other investigators. From the figure we find that the stress intensity factor becomes the smallest for aiD

=

0.3, i.e. a - ao

=

0.15D. This characteristic is originally due to the Chevron bend model.

Ai---~

10 9q~9:,,9,

"..,,""

o

A1---,

10 r f r P.o.,

.

" " " ..^"..

rl--l--+--1°

Figure 4: Finite element meshes and calculated stress intensity factors. (a) - (e) correspond to cases 1 - 5, respectively. (f) shows finite element meshes for case 5'.

Numerical Simulation of Rock Toughness Testing 135

AI I 0 1I

.

I I

1:1

Th

Oucbterlony (hV ) II

\1 ^HFSD

Y.

\I

^\.\ A'

^{1/ 11/ 1}

::\

^{\ \\} Gerstle (BEH) ^{A',vCFD; , ,6/ 'LPO}

\,~TO\ ~ ^~'Gentl' (BEH)

\', " .... - -_ ... .g...;,' ,a',-.- ^_~'T

CTODlFEMI _&

&

I aiD

15 14 13 12

11 10

9 8

.15 .2 .3 .4 .5

Figure 5: Comparison of non-dimensional average stress intensity factors. < 1

> 0,

< 2

> .,

<

3

> .... :

results shown in Table 1; other results : reference [1].

4. NUMERICAL RESULTS FOR TRANSVERSELY ISOTROPIC MATERIAL

The aim of this section is to investigate the stress intensity factors of the Chevron bend model

III Case 3 under the condition of the transversely isotropy. The analysis is achieved using the displacement-type finite element method using the 15-node singular isoparametric elements along the crack front, and the stress intensity factors are evaluated using the formula proposed by Saouma and Sikiotis[3].

·In case of the isotropic body, we may treat only two parameters, i.e., the Poisson's ratio and the Young's modulus. Six parameters are, however, required for the transversely isotropic mate- rials. We firstly examine the number of parameters and decrease it as small as possible, since all combinations of selected parameters are too many for numerical tests.

Assume that the relation between strain and stress is expressed as Ei

=

^aijO'j. Then, most of 2-D anisotropic body cases can be analyzed using following two parameters[4];

\ _ an

A - ,

an (8)

The above equations can be rewritten using engineering quantities of Young's modulus, Poisson's ratio and shearing modulus as following;

VE1E2

P= 2G

12

-

V

^{V 12V 21' (9)}

On the other hand, following characteristics are observed from experiments of various rocks[5];

1. For the Young's modulus, we find the tendency ofE1= E2

>

E3 •

2. For the Poisson's ratio, we find V23

=

^V13 > ^V12' The Poisson's ratio has the value of 0.06

<

v

<

0.45.

Here it is supposed that the isotropic face of transversely isotropic material coinsides with the 1-2 plane.

Taking these results into consideration, we select following three parameters for an 3D anisotropic body in this investigation.

I. Young's modulus; E1

=

E2

=

,-lE3 ( ,

<

1),

11. Poisson's ratio; ^V23== ^{V13= ,-IV12,}

111. Shearing modulus; I/GI2

=

2(1

+

VI2)/E1 and I/G23

=

I/G31 ~ 1/E2

+

1/E3

+

2^{V 23/}E2

where the shearing moduli,G23and G_{31 ,} are approximate values using the equations by Lekhnitskii [6].

In case that the axes of the transversely isotropic body coincide with the Cartesian coordinate axes, we have to consider following three cases:

a) The isotropic face coincides with the 1-2 plane (E1

=

E2

=

,-IE3 ,

,<

^1),

b) The isotropic face coincides with the 2-3 plane (E2

=

E3

=

,E1 , ,

>

1),

c) The isotropic face coincides with the 3-1 plane (E3= E1= A-¹E2 , A= EdEI

<

1) All necessary material constants for the cases of b) and c) are obtained by recurring subscripts of the parameters for the case of a). Note that in case c), the parameter Adefined by A= Ed E1is introduced instead of ,.

Summarizing above considerations, independent material constants are one Young's modulus, , ( or A ) and one Poisson's ratio. We should notice that only two parameters, , ( or A ) and one Poisson's ratio, can be used as independent variables for our investigation, since we treat only the linear problems in this investigation and therefore the absolute value of Young's modulus has no influence to non-dimensional stress intensity factors.

Fig. 6 (a) to (c) show the average non-dimensional stress intensity factors

A

due to the various combination of the Poisson's ratio and the ratio of the Young's modulus for above three cases of a) to c), respectively. The results are normalized by the average non-dimensional stress intensity factor

A

i •o for the isoparametric body with the Poisson's ratio of 0.2. The logarithmic scale is used for the abscissa. The Poisson's ratio is changed between 0.1 and 0.4. Henceforth, we consider mainly on the influence ofE2and E3 to

A

under the condition of the constant E_{1 •}

In case of a) : Fig. 6 (a) shows that

A

increases in accordance with the increase ofE_{3 •} Also, the tendency of the increase is almost linear.

In case of c) : Fig. 6 (c) shows that

A

decreases in accordance with the increase of A. This phenomenon may be explained as following; the Young's modulus E₂₁ which represents the elastic constant in a normal direction to the crack face, b~comes large in accordance with the increasing ofA, and hence the crack opening becomes small.

Numerical Simulation of Rock Toughness Testing 137

0.91.

1123=^V13 (a)

1.01 ^{o· 0.4 c}

A 6: 0.3 ^A

0.2 c

1.

..

~

Aioo 0: 0.1 ^A

a

c

0.99 t _(c)

cH

i 1.03

0.98 3

1.02 c

0.97

.,

1.01 Ae cA 3

V12=^V32 ₈ ^C

A 0: 0.4 ^A

(b)

~

^{Aioo 1. 6: 0.3}

^•

1.01 _II .: 0.2 ⁰

A ^{A 0.99 0: 0.1}

Aioo

1.

•

0 0.98

0.99

i

1131=1I2l fil

o· 0.4

&

0.98 _{6: 0.3}

..

^0.2

0.97 ^{0: 0.1}

0.96

Figure 6: Average stress intensity factors

A

vs. , ( or>.)for various Poisson's ratios. (a) through (c) correspond to cases a) - c) in the text.

In case of b) : SinceE2= E3= , E1 and E1 is constant, both ofE2 and E3 vary in accordance with the change of ,. It indicates that the case of b) includes the effect as shown in both cases of a) and c). Fig. 6 (b) shows that

A

decreases according to the increasing of" but the decreasing is not linear as shown in the figures (a) and (c). The reason is that both of the decreasing tendency of c) and the increasing tendency of a) occur at the same time, since both ofE2and E3 increase in accordance with the increasing of,. We may conclude that the influence ofE2 to

A

is larger than that ofE3 ,since the figure shows the decreasing tendency as a whole.

As the material characteristic changes from the anisotropy to the isotropy, the influence of the Poisson's ratio to the stress intensity factor becomes large, but still smaller comparing to the effect of the ratio of the Young's moduli to the stress intensity factors.

The most important fact observed from Fig. 6 is that the change of the stress intensity factors due to the parameters mentioned above is limited within several percents of the values in the isotropic case. Then, according to the numerical experiments shown in this investigation, it can be said that the anisotropic characteristics do not give large influence to the stress intensity factors.

But, at present we can't conclude that anisotropic effects never appear to stress intensity factors, since above results are obtained through a few numerical experiments.

Through this investigation the authors proposed a numerical simulation method of the standard testing method, which is proposed by ISRM, and they showed several results of the fundamental investigation on the crack propagation and the stress intensity factors. The numerical experiments could not clarify the evident characteristics of the transversely isotropy of the material. The transversely isotropic characteristic is introduced for the simplification of the analysis, but there still remain a number of uncertain factors. It indicates that the proposed theory is one of the approximations of the anisotropic characteristics. And, only a few test problems could be solved through this study.

In further investigation it is expected that the mutual relation among material properties of the Poisson's ratio, the Young's ratio and so on should be clarified under the assumption of the transversely isotropy and also the general anisotropy, and that more appropriate parameters should be selected. In this investigation, only the Chevron bend model is considered. Another model called the Short Rod Model will also be numerically treated in the successive investigation.

ACKNOWLEDGMENT

This work has been supported in part by a grant from Electric Technology Research Foundation of Chugoku.

REFERENCES

[lJ Ouchterlony, F.: A FEM calibration of a core bend specimen with Chevron edge notch for fracture toughness measurements, Report DS 1987:11, SveDeFo, Stockholm, Sweden, 1987.

[2J Ingraffea, A. R. and Manu, C.: Stress-intensity factor computation in three dimensions with quarter-point elements, Int. J. Num. Meth. Eng., VoLl5, pp.1427-1445, 1980.

[3J Saouma, V. E. and Sikiotis, E. S.: Stress intensity factors in anisotropic bodies using singular isoparametric elements, Eng. Fract. Mech., Vol.25, pp.115-121, 1986.

[4J Suo, Z., Bao, G., Fan, B. and Wang, T. C.: Orthotropy rescaling and implications for fracture in composites, Int. J. Solids Structures, Vol.28, pp.235-248, 1991.

[5J Ouchterlony, F.: private communication, 1992..

[6J Lekhnitskii, S. G.: Theory of Elasticity of an Anisotropic Body, Mir Publ., Moscow" 1981.