Numerical simulation of 3-d failure process in heterogeneous rocks

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

NUMERICAL SIMULATION OF 3-D FAILURE PROCESS IN HETEROGENEOUS ROCKS Z.Z.Liang C.A.Tang H.X.Li T.Xu Y.B.Zhang Center for Rock Instability & Seismicity Research, Northeastern University, Shenyang, People’s Republic of China

Abstract: A numerical code RFPA3D based on linear elastic damage mechanics on mesoscopic scale and FEM is developed to simulate the complete progressive 3-D failure and macroscopic mechanical behaviours of rock-like materials. Rock properties such as elastic constants, peak strength, and Poisson ratio are randomly distributed to reflect the initial random distributed weakness in mesoscopic scale. By introducing heterogeneity of rock properties into the model, the RFPA3D code can simulate non-linear deformability of a quasi-brittle behaviour with an ideal brittle constitutive law. Each element in rock samples is considered to be homogenous and obeys a simple constituent law. A self-adaptive loading method is adopted to control the number of elements as less as possible that reach the failure threshold. A specimen was subjected to uniaxial tension. Complete stress-strain curves are obtained. Different from 2-D failure, 3-D samples show more ductile features. Crack initiation, propagation, coalescence and stress distribution all can be observed obviously during progressive failure process. Keywords: heterogeneity; numerical method; failure process; three-dimensional

1. INTRODUCTION It is well known that rock is one kind of highly heterogeneous material that contains many cracks and flaws on the micro-scale. Many traditional analytical models have to ignore some important factors influencing the failure behaviour of material. Heterogeneity is such an example for rock-like materials. However, the failure process analysis is much more complicated than stress analysis. Analytical models based on the hypothesis that rock is homogeneous uniform material are not much suitable. Many mechanical behaviours and phenomenon cannot be explained satisfactorily. Numerical simulation has become one of the most important methods to take into account of heterogeneities and provide an increased insight into material behaviours and structure behaviours. A numerical approach seems essential for fracture mechanisms because of the important interactions between fracture growth and the structural environment. As pointed by Van Mier [1977], with the aid of numerical tools, effects from the structural environment (loading system and loading conditions) can be studied in great detail. To construct heterogeneous models and corresponding analysis methods is one significant task of current numerical mathematics and material mechanics.

In the past decades many investigations have devoted to compressive failure mechanics, including the microfracture initiation, propagation and coalescence [Wawersik, 1970; Lockner, 1991]. The combinations of statistical theory with numerical models such as the lattice model [Cundall, 1998], the bonded particle model [Hart, 1998], and the RFPA model [Tang, 2000] based on FEM are found to be appropriate to simulate the progressive failure of brittle and heterogeneous materials such as rocks. In 1992 E. Schlangen and J. Mier modelled the brittle failure process of concrete-like materials by using a simplified grid model. G. Frantziskonis analyzed the influence of heterogeneity on the displacement field near the structure and obtained an analytical result. Using statistical theory [1997], S.C. Blair and N. Cook investigated the nonlinear mechanical behaviours and the meso-scale heterogeneities of stress field cased by geometry and size effects [1998]. In 1998 C.A. Tang and P.K. Kaiser simulated the progressive failure of rock by using RFPA2D(1998). Investigations mentioned above are all limited to two-dimensional conditions. As we know, many rock mechanical problems, except plane strain problems and plain stress problems, such as crack propagation and crack interaction are related to many directions and cannot be simplified to twodimension problems. Few presentations can be

1

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

found to explain the progressive failure of rock in three-dimensions and related nonlinear behaviours resulting from material heterogeneities. The rock specimen was subjected to numerical tests to investigate some characteristic features of the complete stress-strain curves and the phenomena observed during progressive fracture.

used. We define m as the homogeneity index of the rock. According to the Weibull distribution and the definition of homogeneity index, a larger m implies that more elements with the mechanical properties approximated to the mean value and a more homogeneous rock specimen. The distribution of homogeneity index in RFPA3D is the same as that in RFPA2D [Tang, 2000].

2. 2. BRIEF DESCRIPTION OF RFPA3D

2.2 Elastic damage model in RFPA3D

2.1 Heterogeneity in RFPA3D

m=2 models Figure 1. RFPA3D heterogeneity.

m=20 with

different

It is known that heterogeneous distribution in stress field and deformation field results from the heterogeneous distribution of material properties The RFPA3D code has been developed by considering the deformation of an elastic material containing an initial random distribution of microfeatures to simulate the progressive failure in a more visual way, including simulation of the failure process, failure induced further rupture and failure induced stress redistribution. As shown in Figure 1, all numerical specimens in RFPA3D are discretized into small elements with the same size. Mechanical properties in each element, such as elastic modulus, peak strength, passion ratio and weight, are the same. In other word, mechanical properties of each element at mesoscopic scale level are homogenous. For heterogeneity, the material properties for elements are randomly distributed throughout the specimen by following a Weibull distribution:

ϕ=

⎡ σ ⎤ m σ m−1 ( ) exp ⎢− ( ) m ⎥ σ0 σ0 ⎣ σ0 ⎦

(1)

σ

where σ is the element strength and 0 is the mean strength of the elements for the specimen. For an elastic modulus, E, the same distribution is

The mechanical properties throughout each element are the same and each element is homogenous. In RFPA3D, each element is has an elastic-brittle constitutive law during the loading process. Elastic damage constitutive relation for an element under uniaxial compressive stress and tensile stress is illustrated in Fig. 2. Before the stress of the element satisfies the strength criterion (such as the Coulomb criterion), the elastic modulus is a constant with the same vale as before loading. When the stress increases to a value leading to the failure of the element, in elastic damage mechanics the elastic modulus of the element may degrade gradually as damage progresses and the elastic modulus of the damaged element is defined as follows:

E = (1 − D ) E0

(2)

where D represents the damage variable, and E and E0 are the elastic modulus of the damaged and undamaged elements, respectively. It must be noted that the element and its damage are assumed to be isotropic and elastic, and therefore E; E0 and D are all scalar. In the presented paper, it was assumed that each element has a Mohr-Coulomb failure criterion envelope with a tensile cut-off on mesoscopic scale, and the failure of elements may be either in shear or in tension mode. The tensile failure mode has a higher priority than shear failure mode if the stress state of the element satisfies both the tensile failure criterion and shear failure criterion (Coulomb failure criterion). According to different failure modes, damage variable D can be described as following: If the element is subjected to uniaxial tensile stress, before the tensile stress (the minimum principal stress) of the element reaches its tensile strength σt, the element keeps linear elastic. When σ3>σt, the element fails and the elastic modulus changes to a small value and its strength falls to σrt, which we can call it residual tensile strength. When the tensile stress increases to a more large value,

2

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

σ 1 − kσ 3 ≥ σ c

the element loses his capability of loading, and the corresponding strain is σut. The evolution of damage variable D can be summarized as following:

⎧ 0 ⎪⎪1 − σ rt D=⎨ E ε 0 ⎪ ⎪⎩ 1

(6)

k=

(ε < ε t 0 ) (ε t 0 ≤ ε ≤ ε ut )

(7)

(ε > ε ut )

σ

(3) The variation of damage variable is obtained when the element is subjected to uniaxial tensile stress. In RFPA3D code, rock specimens are subjected to 3-D stress loading, according to Mazars investigation, we extend it from one dimensional damage model to a 3-D model. In Mazars study, he used an effective strain instead of

where c is the uniaxial compressive strength of the element, and K is the ratio of uniaxial compressive strength to uniaxial tensile strength..

φ is the friction angle. The damage variable under

3-D stress is described as: 0 ⎧ ⎪ 1−σ D=⎨ rc 1− ⎪⎩ E0ε 1

tensile stress. Effective strain ε can be defined as following:

ε=

ε1

2

+ ε2

2

+ ε3

1 + sin φ 1 − sin φ

ε1 ≥ ε c0 (8)

where

2

ε1 < ε c0

and

σ rc is the compressive residual strength

ε c 0 is can be described in Figure 3.

(4)

ε where ε 1 , ε 2 , and 3 are principal strains and

σ

Elastic

x

is a function and it can be defined as following:

σc

⎧x x ≥ 0 x =⎨ ⎩0 x < 0

Damaged

(5)

ε

ε ut

ε t0

σ rc

Further Damaged

ε r0

0

σ rt

σt

Figure 2. Elastic damage constitutive law for element in tensile failure mode. Shear damage variable caused by stress can also be obtained as tensile damage variable. MohrCoulomb failure criterion is used as the second failure criterion, that is:

ε1

Figure 3. Elastic damage constitutive law for element in shear failure mode.

2.3 loading procedure External loading is applied with displacement loading at each step. Since in each loading step the RFPA3D model can be regarded as a linear problem, a self-adaptive loading method is adopted, which can automatically determine the loading displacement in each loading step to control the number of elements as less as possible that reach the failure threshold. The external displacement loading is applied much in the same way self-serve control loading system works.

3

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

ni = f (d i ) (9)

λ=

ni Nη

3. NUMERICAL TESTS (10)

d i +1 =

di

λ (11)

We assume that ni is the number of failure elements in step i, di is the external displacement increment. N, η and λ are the total number of elements of the specimen, calculating precision coefficient, and increment coefficient. External displacement increment to be applied in the next step is decided according to the number of elements failed in the previous step. If the displacement applied is so large that it leads to a certain number of elements’ failure, the increment in the next step will decrease to control the calculating precision.

2.4 Implement of RFPA3D with computer tools

User Interface (Visual C++)

RFPA3D

3-D image display

FEM parallel computing

(OpenGL)

(FORTRAN 77, MPI)

is displayed by OpenGL (Open Graphic Library). Linear elastic calculating is implemented with a library package MPI (Message Process Interface) and Fortran 77 language.

Figure 4. Elastic damage constitutive law for element in shear failure mode. RFPA3D is composed of three parts: user interface, 3-D image display and FEM parallel computing. User interface, including pre-process and post-process, for example, data input and output, model establishment and final graphs and curves, is developed by Visual C++ language. 3-D image display can be seemed as a special postprocess of user interface. All 3-D images calculated

A user-friendly pre- and post-processor is integrated to generate the finite element mesh and prepare the input data. Different from other FEM software, RFPA3D code cannot only simulate the initiation of cracks, but also the crack propagation, coalescence and interaction between multi-cracks. The rock specimen containing enhanced regular grains was prepared to undertake uniaxial tension tests. The specimen had the size (200mm×200mm×200mm) and was divided into 128,000 elements (160×80×10). Homogeneity index of the numerical sample was 5. Initial displacement and coefficients of precision are 0.001 and 0.005 respectively. Complete stress-strain curves were obtained by using this self-serve control loading style. The stress-strain curves are similar to the typical curve of brittle rock observed in laboratory (Figure 5). It is found that although the constitutive law for the individual element in the numerical model is nearly brittle, a substantial non-linearity exists before the maximal stress. Even though each element in numerical specimen is homogeneous on meso-scale, the specimens are heterogeneous since mechanical property parameters of all elements are randomly distributed throughout the whole specimen. However, different from plane stress or plane strain tests simulated with RFPA2D, the curves are smoother and sharp stress decrease is not found in both curves. As shown in Figure 3, it is obvious that a clear post-peak region (strain softening) exists from the complete curves. When there is one element that satisfies with the failure criterion, the strength reduces to a small value (only 10% of the original strength). The fracture of the specimen is the accumulative results of elements at mesoscopic level. It should be noted that there is a big difference between the macro strength and the mean value of the elements since a number of elements with lower strength will fail before it reached to peak stress. In conventional experiments it is impossible to measure the microstructure stress and displacement distribution throughout the specimen. In order to study the crack propagation mode, many investigators had to use glass or colophony instead of rock or concert in laboratory. It is so difficult to

4

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

observe cracks in opaque materials. One of the most important advantage of RFPA3D is that we can obtain detailed information about the stress distribution induced stress redistribution (Figure 6 and Figure 7). In the second numerical test, cracks initiated from the corner of the prepared indentations and propagated along the boundary of the enhanced grains. As the loading displacement increased, two cracks propagated ahead to the centre of the numerical specimen. In RFPA3D, stress and displacement are illustrated by different colour grades. There are so many as 2048 colour grades you can select, so it is easy to find any parts where high stress exists. From the high stress parts, both localized zone and fracture nucleation could be observed before the macro failure of the specimens. Strain -0.01

-0.01

0.00

0.00

0.00 0

-20000 -30000

Stress (Mpa) Stress (Mpa)

-10000

-40000

Figure 5. Complete stress-strain curve of the specimen subjected to uniaxial tensile stress

Figure 7. Minimal principal stress of the specimen during progressive failure process.

Figure 6. Numerical specimen with enhanced grains subjected to uniaxial stress loading. Figure 8. Displacement of the specimen subjected to uniaxial tensile stress loading.

5

Paper 2A 02 — SINOROCK2004 Symposium Int. J. Rock Mech. Min. Sci. Vol. 41, No. 3, CD-ROM, © 2004 Elsevier Ltd.

4. CONCLUSIONS A numerical software package RFPA3D is developed to simulate complete progressive failure of rock-like brittle materials. The uniaxial compression test and uniaxial tension test are undertaken with RFPA3D in the presented paper. The following conclusions can be drawn from the RFPA3D model and numerical tests results. 1. Heterogeneity is introduced into the numerical specimen by following a statistical distribution function. 2. RFPA3D code can simulate non-linear deformability of a quasi-brittle behaviour with an ideal brittle constitutive law with a strength and elastic modulus reduction of weaken element. 3. RFPA3D implements a self-sever style loading with a self-adaptive loading method. 4. Complete failure process can be obtained. Crack initiation, propagation, coalescence and stress distribution all can be observed obviously during progressive failure process.

material sand structures [J]. Material Struct, 25:534-542. Tang CA, Liu H, Lee P.K.K, Tsui Y, Tham L.G. 2000. Numerical studies of the influence of microstructure on rock failure in uniaxial compression－Part : effect of heterogeneity. Int J Rock Mech Min Sci 37:555-569. Van Mier 1977. JGM. In: Fracture process of concrete. New York: CRC Press, p. 253. Wawersik WR, Fairhurst C. 1970. A study of brittle rock failure in laboratory compression experiments. Int J Rock Mech Min Sci 7:561-75. Tang CA, Kaiser PK. 1998. Numerical simulation of cumulative damage and seismic energy release during brittle rock failure － Part Ⅰ : fundamentals. Int J Rock Mech Min Sci 35(2):113-121.

Acknowledgements: The study presented in this paper is funded by the China National Natural Science Foundation (No.50174013, 50134040 and 50204003).

5. REFERENCES Blair SC, Cook N.G. 1998. Analysis of compressive fracture in rock using statistical techniques. Int J Rock Mech Min Sci 35(7):837848. Cundal PA. 1998. Formulation of a threedimensional distinct element model－ Part . Int J Rock Mech Min Sci 25:107-116. Frantziskonis G. 1997. Heterogeneous solids - Part I: Analytical and numerical-Dresults on boundary effects. Eur J Mech, A/Solids 16(3): 409-423. Hart R, Cundall PA, Lemos J. 1998. Formulation of a three-dimensional discinct element model － Part . Int J Rock Mech Min Sci 25:117-25. Lockner DA, Byerlee JD, Kuksenko V, Ponomarev A, Sidorin A. Quasi-static fault growth and shear fracture energy in granite. Natrure 350:39-42. Schlangen E, Mier JG. 1992. Simple lattice model for numerical simulation of fracture of concrete

6