ENGINNERING CALCULATIONS AND PRODUCTION EFFICIENCY

N~ERiCAL D~ERGROD~

MODELING OF A FRACTURED ROCK MASS AROUN~ WORKINGS D~C 624.19.001.24

O. N. Zoiotov

When excavating underground workings, the development of zones of inelastic deformations accompanied by such complex rheoiogicai phenomena as plasticity, creep, etc., is observed in the contour region of the enclosing mass. The intensity of this process is largely determined by the structure and degree of fracturing of the rocks. Mathematical modeling of the aforementioned process for various rock media on the basis of universal numerical methods is the most prompt and economical method of obtaining the necessary predictive evaluation of the changes in the natural state of the mass caused by driving a working. The method of initial stresses (MiS) [i] is finding wide use when solving such problems. in particular, in works [2-4] is described a numerical scheme of realization of this method on the basis of mathematical models of rock media proposed in them. Both isotropic [2, 3] rocks and rocks transversely anisotropic [4] in strength properties are examined. Strength anisotropy is due both to weakening of the mass by oriented joint systems - discrete anisotropy - and to a change in the properties in the plane of anisotropy described by a certain continuous function - continuous anisotropy ("feather" joints). Underlying the strength models is the Mohr criterion. The lines of limit states (LLS) of this criterion (pattern of rock strength) are approximated by a piecewise regular function (Lm, Fig. i) defined by the vector of the parameters {Op, Oc, Rmp, Rmc } [2]. Models of discrete anisotropic media realize the idea of the generalized Mohr criterion proposed by G. N. Kuznetsov [5], according to which the strength of the surface of joints of each of the systems is characterized in the criterion by an additional LLS (Lt), defined by an analogous set of parameters. The numerical scheme of MIS described in [i, 2] is developed below in the direction of modeling the step character of load application and consideration of changes in the strength properties of the medium. A. As is known, nonlinear problems do not have uniqueness of solution. Therefore, to obtain the most reliable result possible when realizing iterative schemes, it is recommended [i] to use the method of small increments. This method is usually expressed in a step character of assigning the load. The problem of assigning the value of the load increment in the next step is decisive in this case. in the indicated computational scheme the load is represented in the form of the vector of the point forces {F}. Since in the MiS the zeroth approximation {&60} F for each step {&F} is found by solving the linear problem with a constant matrix of the solving system of equations, the sum of the zeroth approximations will be equal to {60} F. For this reason it makes sense to construct the step process at once relative to the vector {60} F , i.e., the vector of the elastic solution of the problem. Below is given an algorithm for the step process, which is a superstructure on the numerical scheme of MIS described in [i, 2]. As the step of realization of the vector {60} F in the next k stage (k = i, 2 . . . . ) we take the vector {&6k}, equal to {a

j=

(1)

where {60}F, k is the unrealized (remaining) part of the vector; calculated by the formula =

(,

{60} F of the k stage is

-

(2)

for k = i {G}F, I = {8o}r,

h k is the specific part (0 ~ h k ~ i) of the vector {60}F, k realized at the k stage. Translated from Gidrotekhnicheskoe

178

0018-8220/90/2403-0178512.50

Stroitei'stvo,

No. 3, pp. 18-22, March, 1990.

9 1990 Plenum Publishing Corporation

I

Fig. i.

Approximation of the line of limit states of the Mohr criterion.

L~

~An

a

Fig. 2. Diagrams illustrating determination of the value of the load step on characteristic segments of the LLS: a) (-~, OEC]; b) ]OEC , OEp[; c) [OEp, Op*]. The value of h k is found by means of the Mohr criterion determined by line Lm~ (Fig. i) obtained from L m by transformation of the vector of the parameters {op, Oc, Rmp, Rmc} of the latter in proportion to a certain given value (for example, to the va• of the error in estimating these parameters). For each e-th finite element (FE) is found a value of hke, for which the Mohr circle rhk (Fig. 2), determined by the vector of the stresses {o}k e, calculated by the formula

{,,};

=

{:},I_~

+

[o1: h,~ {,~

(3)

where {r e is the increment of elastic deformations corresponding to {A60}F,k, will be the limit with respect to the criterion Lm*. The minimum of the values of hke is taken for h k. The step process ends with fulfillment of the condition h~=1V It~olr.,,
(4)

where a is a certain small value. Formula (3) can be written differently

{oh = {o} ~ - , + h ~ {o0} ~,, w h e r e {o0} k i s t h e v e c t o r index "e" is omitted.

of stresses

corresponding

to the

(

vector

(a0} k.

5

)

Here and below the

179

TABLE i.

Value of the Calcu- characterof fated the computational procases

t e s s (NNN)

Maximum Total num- Number movement of bet of walls of the Iofiterati~MIS ~176 working, cm

4,20 9,20 9,70 10,23 9,66

0 1 II III

lOl

Fig.

3.

63 153 385 199

?

Number of cycles of correction of the properties of the FE

7 3

Scheme of correction of the parameters of the strength properties of finite elements beyond the limit state.

The expression for the limit radius rhk has the form r~k =

~ . k - 1 - - %.k-1 + hk ~ , k - - :,o.k 2

2

+

(6)

+ (~xy+ h ~ , ~ ) 2 . eerzorming elementary

transformation

on ( 6 ) , r~,k

=

we o b t a i n

2 hl~2_{_ 2dhk + r~ , k l rok - -

(7)

where r0,k, rh, k-i are the radii of the Mohr circles respectively for vectors {o0} k and {C}k_ l calculated by the known formula

'

a-

, + ~"'

2

~

~,~.k-1

~

k - - %o,k

2

+ ~y.k-~ ~y,,~.

(s) (9)

9_

According to formula ( 5 ) , a change i n the values of h from 0 to i corresponds to t r a n s f e r of the i n i t i a l Mohr c i r c l e w i t h radius rh,k_ l to the p o s i t i o n of the c i r c l e r~, k. in t h i s case, a c e r t a i n c h a r a c t e r i s t i c p o i n t A of the intermediate c i r c l e s moves in the general case

along a certain curve L h (Fig. 2) intersecting line Lm* for circles r~, k beyond the limits. The unknown circle is the circle rh, k located in the region of allowable values bounded by line Lm*, which touches this line. Since circle rh, k can be located in various characteristic segments of the investigated strength model (-~, OEC ] . ]OEC , oEp[, [OEp , ap*], the solving equation for h k should be obtained for each of them. For the segment (--A, OEC] we have (Fig. 2a)

r~,~=RMo.

(i0)

Substituting (i0) into (7), we obtain the solving equation of this segment 2 h k2 + 2dhk-F r2h.k_l - - Rue = O. r a,~

(11)

For the segment ]ogc, OEp[ we have (Fig. 2b) u

180

(12)

~ . v 84

/////////////////////////,

Fig. 4.

Calculation scheme of underground working: the calculated stages.

i and 2) numbers of

3

\

a

/L . ~ ' ~

b

Fig. 5. Results of calculations: a) zones of limit state of the mass; b) deformation of the contour of the working; 1-5) Nos. of the calculated cases. where rm, k is calculated by the formula r~,~=C~* cos q~*--oE sin q~*.

(13)

Substituting (12) into (7) with consideration of (13) and performing the necessary transformations, we obtain the solving equation for the main inclined (Fig. 2b) segment of line Lm* (,.2 " 2 + 2 ( d + rra,k-1 ~y, sin %~)hk O,k -=~oSin~ . ~M)hk

+r~k_1 ,

--r~

M,k--1

=0.

(14)

Finally, for segment [OEp , Op*] we have (Fig. 2c)

181

(15)

r h , k = Cp-- =Eh,k~ Rup,

w h e r e OEh ' k i s t h e c o o r d i n a t e

of the center

of the circle

rh, k calculated

by t h e f o r m u l a

=Eh,k = =E.k--] + h k : E O , k .

(i6)

Here ~ k ~ ~ k a r e t h e c o o r d i n a t e s o f t h e c e n t e r s o f t h e c i r c l e s r h , k , r k _ l and r 0 , k, r e s p e c t i v e l y . Substituting ( i 5 ) w i t h c o n s i d e r a t i o n o f ( i 6 ) i n t o ( 7 ) , we o b t a i n f o r the segment under consideration 2 ( ~ , k - o~,o)h~ + 2(e + o~0~p)h~ + ~o.~A2~ = o,

(i7)

where ~'EP -----6p ~ aE,k_ I.

(18)

it is easy to see that the solving equations (ii), (i4), and (17) can be described in the same way

Ah2.--~Bh--~-C=O. where f o r t h e Mohr c i r c l e s

beyond the l i m i t

( i 9)

we h a v e

( 20 )

A-JcB-~-C>O. The solution of Eq. (19) in the general case is determined by the known formulas h, =

-- B-+- ] / D

.

-- B - - I / D 2A

(21)

(22)

where

D=B~--4AC.

(23)

investigating the solution of Eq. (19) for each of the characteristic segments of Lm*, we can show that the unknown value of h is determined when A # 0 by expression (21) and when A = 0 by the formula h =--(C/B).

( 24 )

At the same time, the values of h calculated by formula (21) satisfy the inequality

0
(25)

An analysis at the next k-th step of the stress state of each e-th FE for the purpose of determining the value of hke for it is carried out according to the following scheme. i. The position of the centers OE,k- I, SEH ' k (Fig. 2) and radii rh,k_1, r~, k of the corresponding Mohr circles are found, if circle r~, k is in the region bounded by line Lm*, then hke is assumed equal to i and the analysis ends on this, otherwise see paragraph 2. 2. if the points OE,k_ I and OEH~k belong to the same segment of Lm* then the value of hke is calculated by the appropriate rormuias and the analysis ends on this, otherwise see par. 3. 3. If OE,k_ I < SEC and OEH ' k ~]~ ~ the value of hk e is calculated by the formulas corresponding to the segment ]oEC, apE[ and the analysis ends on this (the case ZE,k__1< see A ~ OEH,kE [OEp, Gp*] is not examined). 182

4. if OE, k-l ~]OEC , OEp [ proceed to par. 5, otherwise, if OEE,k < OEC , then the value of hke is calculated by the formulas corresponding to the segment -~, OEC], otherwise, i.e., [OEH,k > OEp] by the formulas of the segment [oEp, Op~]. The anaiysis ends on this. 5. The value of hk e is calculated by the formulas of the segment ]OEC , OEp[ and the analysis ends on this, i.e., the case oE,k--I E[OEp, %*] AOEH,k
,.c,}~.

=

"/

(26)

in the continuous anisotropic model [3] the vectors of the strength characteristics the principal axes of anisotropy are corrected.

of

in the computational scheme described above the procedure under consideration is represented in the form of the next superstructure on the classical scheme of MiS enveloping the step process of realization of the initial approximation (see par. A). The algorithm for correcting the strength properties of the FE is constructed in the following way. After completing the iterative procedure of MiS for the next load step, the stress state of each e-th FE is estimated by a generalized strength criterion, if in this case the condition k~y = (rM/r) > l + ~

(27)

is not satisfied for a certain line L t (criterlon), then the correction of its approximating parameters described above is performed, if in the next inspection of the discrete model at least one FE is found for which condition (27) is not satisfied, a new cycle of the iterative procedure of MIS realizing the discrepancy in the solution {40} obtained for the corrected strength models begins. After its completion the indicated analysis is again repeated until condition (27) is satisfied for all FEs. 183

During performance of the described procedure, information about all "disturbed" FEs, i.e., those elements whose strength properties were changed, is preserved. As a result, after completing the calculation it is possible, in addition to the zones of the beyondlimit state of the mass (kmy > i + e) at the current stage, to estimate also the extent of the "disturbed" zones in which condition (27) was not satisfied at least one of the preceding stages, load steps, i.e., to obtain a complete picture of the development of irreversible deformations in the mass, including the destressed zone. C. As an example of realization of the above-described computational scheme, we will give the solution for a single working without support, being excavated to full section. The rock mass was regarded as homogeneous, isotropic with the following values of the physical and mechanical properties: modulus of elasticity E = 2"103 MPa, Poisson ratio ~ = 0.22, unit weight of the ground 7 = 27.5 kN/m 3, c = 0.06 MPa, ~ = 48 ~ . The origin of the field of natural stresses was assumed in accordance with Dinnik's hypothesis (gravitational origin). A two-stage structurally nonlinear [2] calculation scheme, shown in Fig. 4, was introduced into consideration, in the first stage the calculated origin of the field of natural stresses was modeled; the bulk load 7 was assigned, the pressure of the overlying rock stratum was replaced by a load q = 3.3 MPa distributed along the upper boundary of the calculated region, in the second stage excavation of the section of the working to the full depth was modeled. The discrete model of the finite element method included in the first stage 2400 plane barycentric [7] quadrangular finite elements (2499 points). A series of calculations was carried out with a different character of the computational process. The latter was defined by a binary number character NN~, the corresponding digits of which assigned the regime of turning on (the value of the digit i) and off (the value of the digit 0) the computational procedures of the step-iterative scheme described above. The low-order digit of NN~ determines the type of model: 0 is an elastic model (in this case N ~ = 0); 1 is a nonlinear model with realization of the classical MIS scheme; the middle is the character of the load step procedure (see par. A); the high-order digit is the character of the procedure of correction of the strength properties of the FE (see par. B). Table 1 gives the parameters of the computational process of each of the calculated cases and the values of the maximum movements of the wails of the working. Figure 5 shows the results of the calculations in the form of zones of the limit state of the mass (kmy < i) and deformed states of the contour of the working. in comparison with the elastic model, in the nonlinear model the maximum values of the movements increased by more than 2 times and the configuration of zones with kmy <. i changed substantially. The maximum difference in the variants of realization of the nonlinear model is observed on comparing calculated cases 2 and 4 (Table i): the maximum movements of the wails increased by 10%, the configuration of the zone with kmy < i changed - its depth increased somewhat in several sections. An analysis of the results of calculations and parameters of the computational process of solving nonlinear problems permits the conclusion that in the investigated case we can limit ourselves to the use of only the classical scheme of the MIS with an accuracy sufficient for design practice. The question about the expediency of realizing one or another computational scheme in the general case within the scope of the algorithm described above, especially when computer time is limited, should be solved with consideration of the values of the parameters of the model on the basis of a numerical experiment or estimate calculations. LITERATURE CITED i.

O. Zienkiewicz,

Finite Element Method in Engineering

[Russian translation], Mir, Moscow

(1975). 2. 3.

4. 5.

184

O. N. Zoiotov and v. ~. ~senofontov, "Calculation of subsurface hydraulic structures by the finite element method in a nonlinear formulation," Gidrotekh. Stroit., No. 12 (1983). O. N. Zolotov and V. K. Ksenofontov, "Realization of the Mohr criterion in numerical solution of problems of the mechanics of underground structures," Tr. Gidroproekta, No. 115 (1986). O. N. Zolotov, "Strength of anisotropic rock media under conditions of plane strain," Tr. Gidroproekta, No. 115 (1986). ~. N. ~uznetsov, "Graphic methods of estimating limit states of a fractured mass around mine workings," in: Current Problems of Rock Mechanics [in Russian], Nauka, Leningrad (1972).

6. 7.

K. V. Ruppeneit and Yu. M. Liberman, introduction to Rock Mechanics [in Russian], Gosgortekhizdat (1960). B. V. Fradkin, "An improvement of the finite element method in solving three- and twodimensional problems of elasticity theory," Tr. VODGE0, No. 44 (1974).

ENGINEERING MODEL OF DEFORMATION AND FAiLD~E u~ ROCKS AND CONCRETE D~DER COMPRESSION M. G. Zertsaiov Editor's Note: The selection of a model of the behavior of rocks and concrete under loads is exceptionally important in calculations of the stress-strain state and bearing capacity of hydraulic structures and their foundations. The phenomenoiogicai models now being used in engineering practice do not reflect to the necessary extent the physical processes accompanying deformation and failure of the indicated materials. The three articles of M. G. Zertsaiov, united under the general title "Engineering model of deformation and failure of rocks and concrete under compression," present the model developed by the author, which reflects and takes into account the real structure of materials and their behavior under loads. The first article examines deformation of rocks and concretes from the start of loading to the appearance of the first microcracks. The second and third articles, to be published in the next issues of the journal, examine respectively deformation of the materials from the start of microcracking to fracture (peak strength) and in a "beyond-limit" state. A virtue of the works is that, despite the complexity of the problem, the author obtained results expressed by rather simple relations and permitting their direct use in practice in engineering calculations of concrete structures and rock foundations. Deformation before the Start of Microcracking. Power engineering is related to the design and construction of high concrete dams and underground hydraulic structures, including large-span underground workings of the machine halls of hydroelectric stations and tunnels. As a rule, the construction sites of the indicated objects are located in areas with complex engineering-geological conditions caused by the presence of fracturing, stratification, structural inhomogeneity, and anisotropy of the rock mass. Numerical methods, mainly the finite element method (FEM), are presently used when calculating the interaction of the rock mass with structures, in this case, the stress-strain state of the rock mass-structure system in most cases is determined with consideration of the fact that the rock is modeled by a homogeneous isotropic (or anisotropic) elastic or eiastopiastic quasi-solid medium. The use of such models in calculations by virtue of the characteristics of the rock mass indicated above is a rather rough approximation, and in most cases is quite unacceptable, since a real rock mass, as a rule, is weakened both by individual large joints and by joint systems of various orders. With consideration of this, the development of the FEM should be directed at the development in the calculation scheme the disturbance of the more fully reflecting the behavior of joints and

of calculations of rock masses on the basis of elements making it possible to reproduce continuity of the rock mass and of models of the materials filling them.

At the same time, it is necessary to take into account that the behavior of a rock mass is substantially affected also by those blocks and rock separates into which it was broken by the existing disturbances of continuity, in the existing calculation practice, the behavior of such blocks and separates is modeled, as a rule, by a linearly elastic relation. Translated from Gidrotekhnicheskoe

Stroitei'stvo,

0018-8220/90/2403-0185512.50

No. 3, pp. 22-27, March, 1990.

9 1990 Plenum Publishing Corporation

185