Whole failure process analysis for jointed rock masses based on coupling method of DDA and FEM

(1)

Whole Failure Process Analysis for Jointed Rock

Mass Based on Coupling Method of DDA and

FEM

Huaizhi SU a, Zhiping WEN b and Meng YANG c a

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, China

b

Dept. of Computer Engineering, Nanjing Institute of Technology, China c

College of Water Conservancy and Hydropower Engineering, Hohai University, China

Abstract. The elastic-plastic mechanical behaviour is a typical characteristic of rock mass. The load action will bring on the local destruction, large deformation, even whole failure of rock mass with the discontinuous mediums (e.g. joint, crack and fault). It is a coupling process of the continuous deformation and the discontinuous deformation. The discontinuous deformation analysis (DDA) and finite element method (FEM) are combined to build the elastic-plastic mechanical model. The rock block is divided into the finite element meshes. FEM is used to solve the displacement field and the stress field inside the block. The contacts between the deformable blocks are simulated DDA method. The parametric variational principle is derived to analyze the elastic-plastic problem with above coupling model. The theoretical calculating formulae are obtained from the variational principle. The governing equations of mechanical model are established. The proposed method coupling DDA and FEM is used to implement the simulation and analysis for the deformation process of jointed rock masses around one underground cavern. It is easy to simulate the whole process from plastic to elastic yielding failure, and to the large deformation under the condition of plastic flow or instability.

Keywords. jointed rock mass, failure process, discontinuous deformation analysis, finite element method, coupling method

1. Introduction

It is necessary and important to explore the macro-mechanical properties, mechanical parameters and failure law of rock mass in water conservancy, railways and civil construction etc. Rock structure is regarded as a continuous body or a discontinuous body, and its deformation mechanism and the influencing factors have been studied with numerical analysis methods. If a single continuous or discontinuous analysis method is used to solve the above problem under complex geological condition, its result often has great limitation because some assumptions must be made (Bobet et al. 2009). It has been a trend to couple or combine the different numerical methods to implement the numerical simulation study for jointed rock mass. There exist two kinds of pattern for coupling DDA and FEM, namely, the external coupling pattern and the internal coupling pattern (Cheng and Zhang, 2002). If the external coupling pattern is adopted,

different models are selected in the different domains, and it is satisfied to keep the compatible conditions of force and displacement on the interface between adjacent domains. The internal coupling pattern means that the finite element meshes are set in DDA block to solve the displacement field and stress distribution in the block, while the contacts between the blocks satisfy the block kinematical theory of DDA.

Nonlinear property of rock mass is a very important mechanical characteristic of rock mass. The coupling method of DDA and FEM is rarely used to implement the physical nonlinear analysis of materials, which seriously hinder its application and promotion in engineering practice. In this paper, the coupling method of DDA and FEM is introduced to implement the elastic-plastic analysis of rock mass. According to the characteristics and displacement mode of the coupling model of DDA and FEM, the variation principle of coupling analysis is studied from the nonlinear theoretical system. The

This article is published online with Open Access by IOS Press and distributed under the terms of the Creative Commons Attribution Non-Commercial License.

(2)

incremental theoretical formula of the elastic-plastic coupling method based on variational principle is derived in detail, and the control equation is established. Lastly the proposed method is used in a special example to verify its calculation accuracy and efficiency for solving the continuous and discontinuous deformation of the rock mass elastic-plastic problem.

2. Coupling Model of DDA and FEM

In order to overcome the defects of the description ability of complex deformation in DDA method as a result of the assumption of constant strain coupled DDA and FEM method together, the block element is adopted to simulate the block system which is cutting into concrete shape by the structural interface. Blocks combine into a whole through the contacts of the contacts of the crack structural interface, while the shell meshes are divided by the finite element inside the block (Figure 1) .The stress field and the strain field in the block is described by FEM. DDA method is selected to simulate the contact interface among the block. As a result, this model combines the advantages of DDA and FEM method. It can reflect the continuous or discontinuous locations, and can describe the nonlinear characteristics of the materials and the discontinuous properties of the geometry.

Finite element mesh

Block Three-node element x

y

Figure 1. Finite element meshes in block element.

3. Coupling Analysis of Elastic-plastic

Problems

3.1. Variational Principle of Elastic-plastic Mechanics

(1) Basic equations

The stress and strain state is related to the history in plastic mechanics, hence, the incremental method is used to solve the elastic-plastic problems. For the given time t, assume that ui,

ij

H

and

V

ij represent the deformation and stress

variables of every point on the deformable body in balance; N is the intensive parameter; the body force increment i is dF_i, the surface force increment on the boundary *_p is dp_i, the displacement increment on the boundary *_u is

i

du , the total boundary is * * *p u. All state

variables should satisfy the following basic control equations (Sciarra, 2008).

1) Equilibrium equation , 0 ij j i d

V

dF (in ) (1) 2) Coordinate equation , , 2d

H

ij dui jduj i (in ) (2) 3) Boundary conditions ij j i d

V

n dp (on *p) (3) i i du du (on *u) (4) 4) Constitutive Relations

Stress-strain relationship, yield condition and flow rule:

e p ij ijkl kl ijkl kl kl dV D dH D dH dH (5)

, p,

0 ij ij f V H N d (6) p ij ij g dH O V w w (7)

where is ratio factor of plastic flow; Dijkl is

elastic coefficient tensor; f is yield function;

g g V is plastic potential function.

(2) Principle of minimum potential energy The increment functional of elastic-plastic potential energy for the can be written as (Zhang and Wang, 2006)

(3)

p dp ij ij ij i i i i A d R d d dF du d dp du d H O H : : * ª º 3 _¬ _¼ : : *

³

(8)

where dui is an independent function; is the

control function without the participation of variation, its physical meaning is the flow parameter; R_ij D_ijklw wg V_kl;

1

2

ij ij ijkl kl

A dH dH D dH

is the strain potential energy density in the incremental theory (Sun et al. 2000).

The minimum potential energy principle for elastic-plastic analysis can be described as follows. For the current configuration of the mechanical model on coupling computation of continuous and discontinuous deformation, under the action of given force increment, and in all possible displacement increment field that satisfy the compatible equation as Eq.(2) and displacement boundary condition as Eq. (4), the real solution makes the functional as Eq. (8) reach the general minimum value under the control of constraint condition (Eqs.(5) and (6)). The whole control equation of elastic-plastic analysis with the coupling mechanics model of continuous and discontinuous deformation can be obtained.

(3) Variational principle for elastic-plastic coupling analysis

Assume that the studied structure (namely block system) is composed of N blocks, each block is divided into several finite elements. According to Eq. (8), the potential energy functional of the system can be written as:

* 1 m m m p m m N ij ij ij dp _m _m m _i _i _i _i A d R d d dF du d dp du d H O H : : * ª_¬ º_¼ : ½ ° ° 3 _® _¾ : * ° ° ¯ ¿

³

¦

³

(9)

where m denotes the block number.

Assume that the m-th block is composed of q finite elements. Then the displacement increment mode can be expressed as:

^

`

1 1 2 , , , q m m b r b r r du x y T x y D du x y ½ ° ° _ª _º _' ® ¾ _¬ _¼ ° ° ¯ ¿

¦

(10)

The functional 3*dp is only one more term,

which is related with elastic-plastic deformation, than the functional in the elastic analysis. In view of parametric variable without the participation of variation, the first-order variation for dp

3 can do with the elastic mechanics. The following equation can be obtained by taking the first-order variation of functional 3*dp as Eq.(9).

> @

^

`

> @

^

`

* 1 1 1 m m m p m N m dp m ij ij m _ij N i mi m m q m N i mi m m q m A R d d d dF T D d dp T D d G O G H H G G : : * ª _w º ½ ° _« _» ° 3 _® :_¾ «w » ° _¬ _¼ ° ¯ ¿ ½ ° ° _® ' :_¾ ° ° ¯ ¿ ½ ° ° ® ' *¾ ° ° ¯ ¿

¦ ³

¦

³

¦

³

¦

(11)

The derivation process as follows is done by integration by parts, Green theorem and the given conditions. Based on

1

2

m m m

ij ij ijkl kl

A dH dH D dH , the first term on

right hand in Eq.(11) can be obtained as

^

`

^

`

1 1 1 m m m m N m ij ij m m _ij N m m ij ijkl kl m N m ij kl m A R d d d d D d d R d d O G H H H G H O G H : : : ª _w º ½ ° _« _» _:° ® _« _» ¾ w ° _¬ _¼ ° ¯ ¿ : :

¦ ³

(12)

In the elastic-plastic domain,

ijkl ij ij ij

D d

H O

R d

V

. But in the elastic domain,

ijkl ij ij

D d

H

d

V

. So Eq.(12) can be rewritten as

^

`

1 1 m m m N m ij ij m m _ij N m ij kl m A R d d d d d d O G H H V G H : : ª _w º ½ ° _« _» _:° ® _« _» ¾ w ° _¬ _¼ ° ¯ ¿ :

¦ ³

(13)

(4)

^

`

_{> @}

_^

_`

> @

^

`

1 1 _, m m p m N m ij kl m m ij j mi m N q m m _{ij j} _mi _m q m d d d d n T D d d T D d V G H V G V G : * : : ' *½ ° ° ° ° ® ¾ ' : ° ° ° ° ¯ ¿

¦ ³

¦

³

¦

_¦

³

(14)

Eq.(11) is transformed into

> @

^

`

> @

^

`

* , 1 1 m m p N dp ij j i mi m m q m N m ij j i mi m m q m d dF T D d d n dp T D d G V G V G : * ½ ° ° 3 ® ' :¾ ° ° ¯ ¿ ½ ° ° _® ' *_¾ ° ° ¯ ¿

¦

³

¦

³

¦

(15) Due to the infinitesimal and arbitrariness of

^

G 'Dm

`

in m : and on m p * , the equilibrium equation and the boundary condition can be derived from * 0 dp G3 .

> @

, 0 ij j mi mi dV T dF T in _{: (16)} m m

_{> @}

₀ ij j mi i mi dV n T dp T on * (17) _pm 3.2. Global Control Equation Establishment Structural

:

is discrete into N blocks, and the block m-th is divided into q finite elements. The relationship of strain increment and general displacement increment can be expressed as:

m m b

d

H

B 'D (18)

The stress-strain relationship can be described as m m bp bp b D D B D V H ' ' ' (19)

The following equation can obtained from the above analysis

T T * 1 1 2 N m m m m m dp bp m D K D D F ½ 3 ® ' ' ' ' ¾ ¯ ¿

¦

(20) where T m m bp bp K B D Bd : :

³

(21)

T T , , m m p m m m T F T x y dFd T x y dpd : ª º ' _¬ _¼ : ª º _¬ _¼ *

³

(22) m bp

K is coefficient matrix of DDA; 'F m is external load increment vector of DDA.

The overall potential energy of the system can be obtained by integrating of Eq.(20).

* 1 T T 2 dp D Kbp D D F 3 ' ' ' ' (23) 1 N m bp bp m K

¦

K , 1 N m m F F '

¦

' (24)

where 'D, Kbp, 'F are the vector of overall

general displacement increment, the coefficient matrix and external load increment vector.

The governing equation, namely

bp

K ' 'D F, can be obtained by the extreme variation condition. It is the overall control equation of elastic-plastic coupling analysis. 3.3. Algorithm Implementing Elastic-plastic Model within DDA Method

An established contact condition between adjacent blocks needs to be checked and revised following yielding of the material during elastic-plastic iterations. Therefore, the open-close iterations of the contact algorithm need to be repeated during the elastic-plastic iterations that are carried out to achieve equilibrium. The solution illustrating the interaction between the contact iterations and elastic-plastic iterations is given herein as follows

1. Assemble and integrate generalized stiffness

> @

K and force

> @

F .

2. Detect contact, add or subtract contact penalties, and solve for deformation

> @

D . 3. For element i in blocks, update stress, and

compute yield function f , if f ! then 0 compute new stress state and update force vector.

4. Repeat Step 3 (i.e. elastic-plastic iteration) until converge, and then repeat Step 2 until

(5)

no-tension and no-penetration.

5. Update the positions of blocks and elements, update stresses of elements.

There are many yield criterions available for deformation and stress analysis of geomedia (Priest, 2010). The simplest and generally accepted yield criterion is the Mohr-Coulomb criterion which involves only the maximum and minimum principle stresses. However, the Mohr-Coulomb yield function, which geometrically is a right hexagonal pyramid in the principle stress space, is not smooth, it is numerically difficult to use for the plastic deformation calculation by the associated flow rule theory. The Drucker-Prager yield function, which is an approximation of the Mohr-Coulomb yield function (a smooth right circular cone) in the principle stress space, was used herein for the plastic deformation analysis.

1 2

f

D

I J k (25)

where D and k are material constants; I₁ is the first invariant of stress tensor; J2 is the second invariant of stress deviator.

4. Example Analysis

4.1. Analyzed Object

As an example, a level transportation roadway of underground is studied. The cavity depth is about 350m, cavity height is 9m, and the width is 7m. The internal semi-circular excavation is adopted. There are two non-penetrating faults on the top of cavity, as shown in Figure 2. The main mechanical properties of rock mass and structural surface around the cavity are investigated and listed in Table 1. The mode of changing the initial stress is adopted as the loading method to simulate the variation of the depth of the underground project. The rock mass with unit thickness is analyzed. The whole calculation domain is as follows: horizontal(x direction)

u

longitudinal(y direction) = 70

u

90m. The numerical control parameters were: normal

penalty stiffness 11

5 10

n

k

u

N m

, number of

time steps N=1000, time step size ' t 0.01s; and maximum allowable displacement per time step (i.e. displacement allowed ratio) g=0.02.

The shear spring penalty stiffness can by estimated by k_s k_n

2 1

v

, where v is the Poisson’s Ratio of the block.

The underground excavation is simulated under different depths by changing the initial stress. This initial stress is loaded before excavating, which is coinciding with the actual excavation situation. Take unit thickness rock to analysis, whole calculation domain for: horizontal (X)

u

longitudinal (Y) = 70

u

90m. The computation mesh and boundary conditions are shown in Figure 2.

Figure 2. Underground transportation roadway (Left) and its computational model, boundary conditions (Right).

4.2. Calculation Results Analysis

Figure 3 presents the distribution of plastic domain in surrounding rock and deformation vector graph around cavern. In order to compare with a single method, the DDA method is adopted to simulate above domain by using the same of calculation parameters. The results on the displacement of the feature points around cavity calculated by the DDA method and coupling method are shown in Table 2. Using the parameters in Table 1, the displacement of Point A tends to converge to 10mm within 300 time steps. In order to quantitatively study the relevant rule between the displacement of surrounding rock mass and its influence factors, the displacement varieties of the feature points at the top and bottom of the cavity are calculated and analyzed. During the analysis for a chosen parameter, the other parameters are fixed with the values given in Table 1. In practice, a reasonable range for each of the parameter values can be selected according to the results of field tests and engineering judgment. The discontinuities are assumed to be clean planar faults without surface roughness or infilling. The dilation angles of joints are set to zero herein. According to the numerical calculations, the relation curves between the displacement of two feature points and the variety of each parameter are drawn in Figure 4 and Figure 5. The

(6)

following conclusions can be obtained by analyzing above calculated results.

(1) It is shown in Figure 3(a) that the plastic domains distribute at the top of the cavity and near the faults. According to the calculation process, the elements near excavation edge enter firstly into the plastic, and then toward to the top of cavity and along fault direction to be deeper, which is in line with the actual situation.

(2) Figure 3(b) shows that there are displacements tending into the cavern around the cavern. It reflects bottom swelling caused by the excavation, lateral pressure and the top subsidence. There have deformation dislocation in the top of cavity due to the fault.

(3) It is seen from Table 2 that the calculated displacement of the coupling method is less than that of the DDA method. It is caused that the DDA method is a numerical analysis method based on blocky theory, and its calculated displacement is larger than the corresponding ones of FEM.

(4) According to Figure 4(a), the overall displacement in the feature points is reducing with the increasing of the elastic modulus (E) of rock mass when E is in range of [5, 25GPa]. Decreasing the elastic modulus E (from 10GPa to 5GPa) resulted in a significant increase of 7.45mm (from 9.68mm to 17.13mm, an increase of 77%) in the displacement of Point A. Whereas the increase from 10GPa to 25GPa only produced a decrease of 0.62mm (about 6.41% ), namely from 9.68mm to 9.06mm. As E is large enough, the influence of the variation of E on the deformation is not obvious. This change law of displacement due to the elastic modulus can be ascribed to rock stiffness. Figure 4(b) shows that the displacement is only partly sensitive to variations of the value of rock mass the Poisson's ratio (ȣ), with only a difference of 1.59mm (increase rate is about 16.16%) between the upper and lower bounds of 0.1 and 0.45 for Point A. In general, the Poisson's ratio directly contributes to the stress field and deformation of side walls which was not considered herein. For deformation parameters, results (Figure 4) show that the displacement field is more sensitive to elastic modulus than Poisson's ratio.

(5) The following case can be seen from Figure 5(a). The displacement varies with selected scope of the parameters when the

cohesion(c) of structural interface is in the range of [0, 1MPa]. It is divided into two sections, namely, [0, 0.25] and [0.25, 1]. A rapid increase of 37.47mm of the displacement of point A was seen when the cohesion c was decreased (from 0.25MPa to 0.0). However, little effect of increasing cohesion (from 0.25MPa to 1.0MPa) on the decrease of the displacement. This can be interpreted by the low value of the friction angle (i.e. 25°), in which condition, shear deformation develops significantly under the combined effect of low friction angle and cohesion. Results (Figure 5(a)) show that increasing the cohesion of the major structural interfaces with engineering measurements can reduce its displacement.

(6) According to Figure 5(b), a similar degree of sensitivity to that of elastic modulus was seen in the results of varying the friction angle (

M

) between 15° and 45°. For models characterized by structural interfaces, a rapid increase of 68% was found when

M

was decreased from 25° to 15°due to the greater shear displacements driven by the high initial ground stress acting on the models. Whereas little change was seen when

M

increases from 30° to 45°. Figures 5(a) and 5(b) show that discontinuities have a significant influence on the mechanical behavior.

(7) The final input parameter tested with respect to displacement was the tension strength. It is shown from Figure 5(c) that it is not remarkable influence to the displacement when the tensile strength of the structural interface is in the range of [0, 2MPa].

(8) From the above analysis, the deformation and possible destruction forms of underground surrounding rock are no more than a strength problem of surrounding rock. Especially, when there are joint, fissures, faults etc. around the cavern, the integrity and uniformity of the rock is damaged in different degree, and its mechanical behavior becomes a complex comprehensive question. The rock mass displacement is mainly controlled by the mechanical properties of structural surface when the rock mass intensity is high. In the design of the supporting structure, it should be paid particular attention to increase the residual strength of surrounding rock, namely structural surface cohesion, tensile strength and internal friction angle, etc. The displacement of

(7)

surrounding rock is not uniform after the underground caverns excavation, the displacement destruction happen to one or a few part at first, which leads to the overall instability of underground cavity. The first damaged parts should be the key parts of the supporting structure.

(9) The coupling method of DDA and FEM can not only consider the discontinuous deformation characteristics of rock mass, also can describe accurately the block deformation filed. At the same time, the elastic-plastic properties and other nonlinear mechanical properties of the material can be considered fully as a result of the introduction of variation principle.

(a) (b)

Figure 3. Distribution of plastic zone around rock masses (a) and calculate deformation vector graph (b).

(a)

(b)

Figure 4. Displacements versus rock parameters: (a) Elastic modulus ~ Displacement, (b) Poisson's ratio ~ Displacement.

(a)

(b)

(c)

Figure 5. Displacements versus parameters of fault: (a) Cohesion ~ Displacement, (b) Friction angle ~ Displacement, (c) Tensile strength ~ Displacement.

Table 1. Parameters of rock mass and fault Materials Elastical modulus

(GPa) Poisson's ratio Density (kN/m3) Rock mass 10 0.24 25 Fault 㸫㸫㸫

Materials Tensile strength (MPa) Cohesion (MPa) Friction angle(deg) Rock mass 5.0 2.5 58 Fault 0 0.25 25

Table 2. Calculated displacement of typical points with DDA method and couple method

Method Displacement(mm)

A B C1 C2 D1 D2

Couple method 9.68 8.63 5.09 8.01 7.87 6.50 DDA method 12.33 10.85 8.12 10.55 9.94 9.26

(8)

5. Conclusions

Rock mass is a complicated multi-body system with high discontinuity, heterogeneity and anisotropy characteristics, etc. It shows the great nonlinear and non-elastic mechanical behavior. The structure and material deformation of actual rock engineering (e.g underground cavern, dams, etc) is usually a series of complicated gradual deformation and destructive process from accumulation of small deformation, damage, evolution, the formation and expansion of the macro-crack, holing through faults until the occurrence of the rigid body motion and other large deformation, large displacement of the scattered body. The above problems can be solved well by the elastic-plastic coupling methods of FEM and DDA.

(1) The variational principle is derived to analyze the elastic-plastic problem with the coupling method of FEM and DDA. The theoretical calculating formulae are obtained from the variational principle. It is an expansion to the coupling analysis theory of FEM and DDA, which has certain theoretical significance. The method proposed in this paper is used to simulate the excavation of the underground cavern. Its feasibility and superiority has been proven. The method can response the material nonlinear and geometry discontinuous characteristic of rock mass. It is easy to simulate the whole process from plastic to elastic yielding failure, and to the large deformation under the condition of plastic flow or instability. So it can be widely used in future in actual project.

(2) The quality of the results is completely dependent on the ability to acquire good data that accurately represents the field conditions. Discontinuity data are particularly important, such as strength data. Data must be measured and recorded for each major discontinuity, and ranges are useful in addition to average values. Therefore, the proposed method should be combined with the feedback analysis. According to the observed deformation data, the mechanics parameters of rock mass are obtained by the feedback analysis methods, and then the

deformation and stability of rock mass are forecast, which can provide the reliable basis for the engineering construction, operation and the design of reinforcement measure.

(3) This study also shows that the ideas and methods of the numerical analysis need to be improved constantly to better simulate the engineering mechanical properties of the practical rock material, and then better guide the engineering practice. An suitable numerical simulation technique need to be selected according to the characteristic, degree, timeliness, economical and reasonable factors for one actual rock mass problem.

Acknowledgements

This research has been partially supported by National Natural Science Foundation of China (SN: 51179066, 41323001, 51139001), Jiangsu Natural Science Foundation (SN: BK2012036), the Doctoral Program of Higher Education of China (SN: 20130094110010), Non-profit Industry Financial Program of MWR (SN: 201301061).

References

Bobet, A., Fakhimi, A., Johnson, S., Morris, J., Tonon, F., Ronald Yeung, M. (2009). Numerical models in discontinuous media: review of advances for rock mechanics applications. Geotech. Geoenviron. Eng. 135, 1547–1561.

Cheng,Y.M., Zhang, Y.H. (2002). Coupling of FEM and DDA methods. Int. J. Geomech. 2, 503–517.

Priest, S.D. (2010). Comparisons between selected three-dimensional yield criteria applied to rock. Rock Mech. Rock Eng. 43, 379–389.

Sun, D.H., Hu, Q., Xu, H. (2000). A neurocomputing model for the elastoplasticity. Comput. Meth. Appl. Mech. Eng. 182, 177–186.

Sciarra, F.M. DE. (2008). Variational formulations, convergence and stability properties in nonlocal elastoplasticity. Int. J. Solids Struct. 45, 2322–2354. Zhang, H.W., Wang, H. (2006). Parametric variational

principle based elastic-plastic analysis of heterogeneous materials with Voronoi finite element method. Appl. Math. Mech. 27, 904–912.