- Email: [email protected]
Efficient numerical modelling of faulted rock using the boundary element method
Int. J. Rock Mech. Min. Sci. & Geomech, Ahstr. Vol. 31, No. 5, pp. 485-506, 1994
Pergamon
0148-9062(93)E0030-R
Copyright '~i'; 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0148-9062/94 $7.00 + 0.00
Efficient Numerical Modelling of Faulted Rock Using the Boundary Element Method G. BEERt B. A. POULSEN~ This paper presents an efficient method of modelling faulted rock using the multiregion boundary element method. The boundary element method is ideally suited for modelling problems in rock mechanics because of the requirement to descretize surfaces only. This contrasts with the volume discretization needed for the finite element method, finite difference method and the distinct element method. The disadvantage of the boundary element method is that the rock mass is assumed to be isotropic and elastic. However, it is possible to model the non-linear behaviour on fault planes by using the concept of multiple regions and by connecting these regions via non-linear springs or Goodman type joint elements. The disadvantage of this method has been that joint stiffnesses have to be specified and that a great number of iterations are needed because the stiffness coefficient of a spring may change very rapidly from a very high value (for example when a joint is in compression) to a very low value (when a joint is in tension). This paper describes a novel method for the treatment of faults with the boundary element method which eliminates the need to specify joint stiffness parameters and which is able to model significant fault movement with very few iterations. The numerical scheme makes very efficient use of computer resources because the non-linear iterations only involve degrees of freedom at the interfaces. Examples of applications of the new model to practical problems in mining and geological engineering are given.
INTRODUCTION The boundary element method is ideally suited for analysing problems in rock mechanics because only surfaces have to be discretised into boundary elements. Figure 1, for example, shows the boundary element mesh used for analysing mine excavations. With the conventional boundary element method only elastic and homogeneous domains can be analysed. Practical problems in rock mechanics do not comply with these ideal conditions and this has been one of the main reasons why the method has not found wider acceptance in rock mechanics. In reality, the rock mass is non-homogeneous and elasto-viscoplastic and has joints and faults. In general, if non-linear material behaviour is to be modelled with boundary elements then a volume discretization is needed. However,
the influence of discrete joints and faults present in the rock mass can be considered with the boundary element method in two ways: (1) by using displacement discontinuity boundary elements. These are special boundary elements which have two opposing surfaces; (2) by using the multiregion concept in the conventional boundary element method. The first approach is explained in detail in [1]. The present paper will concentrate on the second approach. THE MULTIREGION BOUNDARY ELEMENT METHOD
This method of dealing with piecewise homogeneous domains and joints and faults has been used by various authors [2-4]. Basically the approach involves two or more boundary element regions which are connected to
tlnstitut fiir Baustatik, University of Technology, Graz, Austria. :[:CSIRO Division of Exploration and Mining, P.O. Box 883, Kenmore, Qld 4069, Australia. 485
486
BEER and POULSEN:
BEM MODELLING ( ) F I.AI.I.Tt!I) ROCK
rrMVUE /
,/
:
'~\
iI
!/, /
/
SCALE:
j
12.255
Fig. 1. Boundary element mesh of a mining excavation.
each other. Some kind of joint behaviour is then assumed at the interfaces between regions. A boundary element region is a region which is bounded by a boundary element mesh and, in the case of an infinite and semi-infinite region, by the surface of a
INFINITE
sphere or half sphere with a radius of infinity. Figure 2 shows an example of a finite, infinite and semi-infinite region. Each region can be assigned elastic material properties, E, v etc. For a region one may establish a relation
C) C3 FINITE
SEMI-INFINITE
Fig. 2, Boundary element regions in two dimensions.
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
487
Fig. 3. Boundary element regions with a fault at the interface.
(a)
(d)
I
kP t",, ~.
(b)
(e)
(c)
(f) Fig. 4. Methods 1 and 2 for non-linear interface behaviour.
488
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK
sT
ts
n
Fig. 5. Local coordinate system at the interface.
between the tractions (boundary stresses) and the displacements. Using Betti's law we may write
c(e)u(e)
- .Is T(P, O)u(O) dS
+ ;sU(P, Q)t(Q)dS
(1)
where e(P) is the free term, t(Q) and u(Q) are vectors containing the values of traction and displacement at boundary point Q and T(P, Q) and U(P, Q) are the fundamental solutions for the traction and displacement at point Q due to orthogonal unit loads at P in the coordinate directions. Using isoparametric boundary elements [5] the displacement and traction can be approximated by J
u(~)
= ~ N:(~)u/ j=l
Substituting equations (2) and (3) into (1) the following matrix expression relating the displacements {u} and tractions {t} at all the nodes of the boundary of a region can be obtained using the point collocation method [5]
(2)
[Al{u} = [Bl{t}
One may now combine two or more regions to analyse problems with joints or faults. Figure 3 shows an example of the combination of two semi-infinite regions to analyse an excavation with one fault. For this example we can rewrite equation (4) for region I
[A]' {u}' = [B]' {t}'
(5)
[A] '~{ul H = [B]" {t}"
(6)
and for region II
J
t(~) = ~ Nj(~)t/
(3)
l-I
where uj and tj are values of u and t at n o d e j and N:(~) are shape functions of the intrinsic coordinates ¢ and r/.
At the nodes of each region, except for the interface nodes, either tractions or the displacements are known. On the interface nodes both displacements and
i!
tf = - r . V 3
Its
>
-tnj
tan ej + cj
(a)
(4)
(b)
Fig. 6. Conditions for a joint slip satisfied 2-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (shear link broken) for next time step.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK ~rl q~; +
489
c3
IR1 = F s.
lsl
tR2 = F s.
ts2
i73
ts 1
(a)
(b)
Fig. 7. Conditions for a joint slip satisfied 3-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (shear link broken) for next time step. tractions are unknown. Let us rearrange equation (5) as follows:
[A f, -- Bi]t{::} ' = [n f, - Ai]'{:f, } !
(7)
Where u r and t r are vectors containing the displacements and tractions of nodes not on the interface and u i and f are corresponding vectors for the nodes on the interface. Assuming that all tractions t f (i.e. the loads applied to the free surface) are known we can write:
~fuf] '
[.$.] ~ ~ = {b} ! - [Ai]'{ui} ! ( ti )
(8)
where {b}'= [W]'{tr} '
and
f~ '
ej
faro "~I =~foj
+
[ui 1 K i {ui} '
{t i }1 = {ti0 }1 + [Ki]l{u i }!
{u'} !
=
{¢}!!
(] 1)
and provided no tractions are applied at the interface (9)
{t'}! + {f}" = 0
(12)
~
" "
"
(a)
(10)
For the case where the fault/joint does not slip or dilate the conditions on the interface are:
[.~.]'= [Af, - B i ] '.
By treating each column of [Ai]l as the right-hand side, the system of equations (8) can be solved to give
{u
The first part of the solution represents the displacement of the "free" nodes u f° and the traction of the interface nodes f0 assuming the displacements at the interface are zero. Each column of the second part of the solution represents the corresponding displacement and traction values for unit displacements of the interface nodes. Thus [K i] represents a "stiffness" matrix which relates the traction of the interface with the displacements at the interface. For BE region I we can write the following relationship between the interface displacements and interface tractions
R
__tn
(b)
Fig. 8. Conditions for a joint separation satisfied 2-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (all links broken) for next time step.
490
BEER and POULSEN:
BEM M O D E L L I N G OF FAUL'I-EI) RO('K
tsR1 = -tsl
ls 1
(a)
(b)
Fig. 9. Conditions for a joint separation satisfied 3-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (all links broken) for next time step.
which yields the following system of equations for the displacements at the interface ([K']l + [Ki]II){u i } = {ti° }1 ---]- {t'°} u
(13)
The interface stresses f may be obtained by substituting the interface displacements u ~ into equation (10). The displacements at the free nodes u f can be obtained by substituting into equation (9). For the case where the joint slips or dilates equations (11) and (12) no longer hold and an iterative scheme has to be employed. This is discussed in detail in the following section. MODELLING OF JOINTS AND FAULTS There are basically two main approaches for the modelling of slip and dilation on fault planes: (1) By connecting the nodes at the interface by non-linear springs. Initially a high stiffness is assumed for the spring. This stiffness is reduced in the direction normal to the interface if the yield condition for joint opening is reached and in the direction tangential to the interface if the yield condition for joint slip is reached; (2) By connecting the nodes at the interface by rigid links in the direction perpendicular and tangential to the interface. If the yield condition in tension is reached then the link is broken. If the joint yields in shear then the link in the tangential direction is broken. Figure 4 illustrates the two approaches. In method 1 the pair of interface nodes is connected by springs. Initially the springs are assigned high values of stiffness in the normal direction, k,p and shear direction, k~ [Fig. 4(a)]. If the yield condition for joint slip is reached then the shear stiffness is reduced to a small but non-zero value k~ for any subsequent load step [Fig. 4(b)]. If the
yield condition for tension is reached then both the shear and normal stiffness are reduced [Fig. 4(c)]. In method 2 the nodes are connected by rigid links which prevent shear and separation [Fig. 4(d)]. If the yield condition in shear is reached the link which prevents slip is removed [Fig. 4(e)]. If the yield condition in tension has been reached then all links are broken [Fig. 4(f)]. The advantage of method 2 is two-fold: (1) there is no need to specify initial and residual values of joint stiffness, quantities which are not easily obtained in the field; (2) convergence of the method is significantly faster since we do not have to deal with large differences in stiffness between iterations. For these reasons method 2 has been used in the present work and this will be explained in more detail here.
n
Asl Fig. 10. Sign convention tbr joint displacement.
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D ROCK
491
'4
!
!
I
I
I
I
URS /
URS ..,
J
1.1111
ORS
ORS
QRS
\ \ -
\
ilCm:.
\ \ \
PLS $o¢q.
Ul
I; -
ors I
\ PLS
I I
~. _ _
\ HVS
~'Dll,
"~
440~
\
\
i STZ
x,,, -
STZ
\1
400 LLr~L
~°L I
t O
PLS I
I
[]
'~
!
on1
:
Ji
\~1
. |,,_,~ !_,~s ± : - . , _ _ , ~ : . ~ _ _ ~ _ ~ ,--..~
10950 N Og,dldJ~4k~l9~1
Fig. I I. Stope section, Hellyer mine, Australia.
492
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
LOCAL INTERFACE COORDINATE SYSTEM To be able to implement the method, the stresses at the interface have to be expressed in a local system of coordinates with one direction pointing normal to the joint and the other one or two directions being perpendicular to this direction (Fig. 5). The transformation from global traction/displacement ti, ui to local tractions/displacements t~, u~ at a node./is given by
u~ = T iuj
( ! 4)
t~ = Tjtj
(15)
Iij =
Us
X,=
Is
-n>,
u~, Us 2
t~=
.~),
~/:
'~"2 ~
"'2 ~
$2:
(19)
It'} = [T]{¢}' = [T]{f° }' + [T][Kq'[T]T{u' }'
{t' }' = {t o' }' + [K']'{u' }J
(21)
{t o' }' = [T] {f0},
(22)
[K' ]~ = [T] [K']' [T] T
(23)
where
The matrix [T] is given, for example, for four interface nodes, by
ni~
(17)
t,, ts z
08)
Ii
T2 0 0
T3
0
0
(24) T4
The transformation matrix Tj is computed with the outward normal n/to the node i. For three dimensional
ITE dDARY IENTS FAULT
91-400
(20)
or
For 3-D problems we have (Fig. 5)
u;=
.~/,
Equation (I0) can now be rewritten in terms of local components as:
Where for 2-D problems
{.n}
1l ll: ]
H ~
T,=
-D E
o 1 = 16.3MPa = 9.9MPa
:66~ \
~MPa
°
/
//ff -
b = 159 ° d=l°
< r T ~
//
X/
Fig. 12. Boundary element mesh of the Hellyer mine. Also shown is the in situ stress tensor.
~-~
BEER and FQULSEN:
BEM MODELLING
Fig. 13. Contours
OF FAULTED
of normal stress on Jack fault.
Fig. 14. Shear stress along strike.
ROCK
493
494
BEER
and POULSEN:
BEM MODELLING
Fig. 15. Shear stress up/down Fig. 16. Contours
of slip along
OF FAULTED
dip.
strike after two iterations
ROCK
BEER and POULSEN:
BEM MODELLING
OF FAULTED
Fig. 17. Contours of slip up/down after two iterations.
ROCK
495
496
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D R O C K
Fig. 18. Stress vectors at 360 level, no slip on Jack fault.
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D R O C K
problems vectors st and s2 are computed using the following vector x product (Fig. 5). sl = n x y
(25)
s2 = sl x n
(26)
where y is a unit vector in the y direction, i.e.
I f n points exactly in the y direction then equation (25) is replaced by s, = n x x
(28)
497
For smooth interfaces without kinks there is a unique value for nj at each node. For interfaces with kinks average values of n.j. can be assumed. If there are initial stresses tr0 present then we have to work out the corresponding initial tractions in the local interface coordinate system. To convert the initial stress to tractions at the interface we have t.,.0 = trx.onx + zx:..on,. + rx:.on-
(30)
t:,o = a:.on,. + zxy.onx + z:y.on..
(31)
lz, 0
(32)
=
tTz,on: --I- "Cxz.onx --l- Tyz,0ny
The tractions in the directions normal and tangential to the interface are given by
where x is a unit vector in the x direction t~ = Tt0
(33)
These initial interface tractions have to be added to the vector of traction due to a rigid interface {to'} ~.
[
10.000 11.gTg
Fig. 19. Stress vectors at 360 level, slip on Jack fault after two iterations.
lniiiai Discovery atx,y,z-0.0,0
Fig. 20. Structure
1
0
1
250
Fig. 21. Structure
contours
of Granite,
Ultramafic
contact
in plan view. Elevation
given in metres.
1
5OQm
contours
of the Lucky
and Golden 49x
faults in plan view. Elevation
given in metres.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK O N S E T OF SLIP AND SEPARATION--YIELD CONDITIONS
The condition for the onset of slip and separation (yield conditions) may now be written in terms of the values of tnj and ts,j, t~2/at interface node j as
Stip Fs = I tsjl + t,j tan q~j- CJ= 0
Separation
F,=t.j- ~=0
499
ITERATION (TIME STEP) PROCEDURE An iterative (time step) scheme has to be used for the solution of problems which involve slip and separation. For the first iteration (time step = 0) we assume rigid links between all the nodes on the interface. For the two region problem in Fig. 3 we can write
{u'},=o = {u'}' = {u' }"
(35)
and assuming no external forces are applied at the joint
{t'},=0 = {t'}l + {t'} t' = 0 (36) where the tension positive sign convention is used and ~bj is the angle of friction, cj is the cohesion and ~ is the This leads to the following system of equations tensile strength at node j. For three dimensional prob(37) ([K']' + [K']"){u'},:o = {t'°}' + {t'°} '' lems This system can be solved to give the displacements Itsjl - ~ lJ dr- t s2j 2 (34) {u' }, = 0 in the joint plane. The tractions at the joint plane
FEHVLIE
I
I / I
I •
I I
/ /
/ I
SCALE: I~1
Fig. 22. Boundary element mesh used for the geological modelling example. RMMS M ,5--H
201.507
500
BEER and POULSEN:
BEM MODELLING OF FAULTI:D RO(K
in the local (joint) coordinate system can be obtained by substituting {u'},=0 for {u'} ~ in equation (21)
{t'}',=0 = {t'°}~ + [Wl{u' I, =,,
(38)
where {t '° }~ includes the initial joint tractions mentioned above. Once we know these tractions we can check the yield conditions and determine if the joint slips or opens. The action taken for the next iteration depends on the value or F, or Ft- If both F, and F~ at an interface node are zero then no action is taken for that node. If F, is greater than or equal to, zero then the condition for joint slip is satisfied. If F, is greater than zero then at one stage during the increment the yield condition has been reached and the tractions at the end of the time step violate the yield condition. The component of the shear traction which is above F, = 0 is also called "excessive" shear traction. The procedure for joint yielding in shear is as follows (Figs 6 and 7): (1) the link(s) in the direction perpendicular to the vector normal to the joint plane are broken, (2) any "excessive" shear tractions are applied to the node in the opposite direction.
For yielding in tension the procedure is as tk~llows ( Figs 8 and 9): (I) all links are broken: (2) "excessive" shear and normal tractions are applied in the opposite direction. PROGRAMMING CONSIDERATIONS
The implementation of the algorithm for connecting and disconnecting degrees of freedom is explained here in more detail. The iteration time steps involve the displacements on joint planes only via the "region stiffness matrices" [Ki] ', [Ki] u, etc. The system of equations is solved using a frontal solution method [6]. In this method each variable is assigned a "destination" which determines where its stiffness coefficients are to be assembled. If two elements are connected by rigid links at a node then the stiffness coefficients for that node are assembled in the same location (that is they are added). If links are broken then the stiffness coefficients corresponding to the direction of the broken link (degree of freedom) are simply assembled in different locations. In the software implementation the allocation of "'destination" is handled as follows: each node is assigned a "restraint code" for each degree of freedom which is set to I for each connected degree of freedom
FEHVUE
2 / /
(a)
(b)
J
(c)
(d)
Fig. 23. The four boundary element regions of the geological modelling example.
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
and to 2 for each disconnected degree of freedom. Initially all "restraint codes" are set to 1. If a link is broken then the "restraint code" for the corresponding df. is set to 2. JOINT OPENING/CLOSING, SLIP/STICK MONITOR During the iteration (time stepping) it is important to keep track of joints which are closing after being open and joints where the direction of slip is reversed. This is done by a device known as a opening/closing, slip/stick monitor. The joint accumulated relative displacements (i.e. slip and opening) at time step t at node j can be computed from the absolute displacements at the interface by
Opening Anj(t) = u,i(t " ) -- ui,j(t)
(39)
501
Slip I As,i( t )
= u'~',j(t)-
As2j(t)
=
u'~,j( t )
(40)
Slip 2
(41)
.','~j(t) - u~j(t)
Figure 10 shows the sign convention used for the relative joint displacement. The monitor is implemented in the software in such a way that a rigid link is reestablished for the case where the joint closes again or the slip is reversed. Joint closing is detected when Ajj(t) changes sign and slip reversal occurs if either Aslj(t) or As2/(t) changes sign. APPLICATIONS
The numerical method presented above has been implemented in program BEFE [8] and extensively tested. Comparisons were made with theoretical solutions or solutions obtained with other methods or
FEMVUE
DISP~. SCAkE: I
I
4.545
SCALE: I
Fig. 24. View from above, note the slip of Lucky and Golden faults.
I
201 .Mg
502
BEER and P O U L S E N :
BEM M O D E L L I N G O F F A U L T E D R O C K
software. Test examples are reported in [7] and will not be repeated here. Instead we will show applications in mining and geological engineering.
MINE EXCAVATIONS IN FAULTED ROCK
This application relates to the modelling of the excavation of ore by open stoping and caving at the Hellyer Mine in Tasmania, Australia. The region modelled is in the northern area of the mine. Figure 11 shows a typical section through the caved stope. The geological feature modelled is the Jack fault, which is sub-parallel to the caved stope and limits it in the upper part. Figure 12 shows a perspective view of the boundary element mesh for modelling the caved stope, the 87A stope, the 400-91 stope and the Jack fault. The mesh was prepared from horizontal and vertical sections using the preprocessor F E M C A D [9]. Four node boundary elements with a linear shpae function were used. Although
the Jack fault consists of two surfaces with boundary elements on each side of the fault, only one surface had to be specified with the program automatically generating the other side. Although there are about four different materials only one material was assumed lbr the purposes of the analysis. The elastic properties of this material are: E = 30,000 MPa v =0.3 These values were derived from uniaxiai compression tests on drill core. The properties of the Jack fault were assumed to be Cohesion = 0.0 MPa Angle of friction = 20 '-~ Dilation angle = O These properties were derived from observation of the fault surface where it is exposed, and from experience.
FEMVHE
J
DISPL. SCALE: 1 I 4.545
SCALE:
I
Fig. 25. View from below, note the slip of the faults.
I
164.141
BEER
and POULSEN:
Fig. 26. Normal
BEM MODELLING
stress distribution Fig. 28. Contours
on lithological
OF FAULTED
contact
of slip on lithological
(Fig. 2% Overleaf)
ROCK
after two iterations. contact.
503
504
BEER
and POULSEN:
Fig. 27. Contours
BEM MODELLING
of slip up/down
OF FAULTED
dip on Lucky
and Golden
ROCK
faults.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK Virgin stress measurement results from only one site were available and these have been used. The stress field was assumed to vary linearly from the surface and have the following values:
At surface:
all stress zero
At depth of 280m: a~ = 16.3 M P a a2 = 9.9 M P a tr3 = 3.8 MPa
Figure 22 shows the boundary element mesh used for the analysis. The mesh consists of four boundary element regions which are shown in Fig. 23. Note that the mesh does not have any internal elements, only surface elements thus making mesh generation easier. The following material properties are used: Rock mass
bearing = 159 ° bearing = 66 ° bearing = 249 °
dip = 1° dip = 71 ° dip = 19 °
These stress tensors are shown in Fig. 12. The following figures show some results of the analysis. Figure 13 shows contours of the stress normal to Jack fault. One can clearly see the destressed zone adjacent to the excavation and the stress concentration at the corners where the Jack fault surface is exposed. Figures 14 and 15 show contours of shear stress along strike and up dip. The contours of slip along strike after 2 iterations (time steps) are shown in Fig. 16. The m a x i m u m values of slip of 15 m m occurs to the south of the caved area. Figure 17 shows contours of slip up/down dip after two iterations. The maximum value of slip ( 1 4 m m ) occurs where the Jack fault surface is exposed to the caved area. The predicted directions and values of slip were found to agree well with observations although actual magnitudes were found to be larger than those computed. Figures 18 and 19 show the effect of fault movement on the stress distribution at the 360 level. Figure 18 shows the stress tensors before any movement on the fault occurs, whereas Fig. 19 shows the stress tensors after the Jack fault has slipped by 15 mm. The main effect is to turn the stress tensors so that they are more perpendicular to the fault in the area to the south of the caved zone. Consequently the 87 Pillar, which was in a stress shadow before, has become more highly stressed.
Geological modelling This application relates to gaining a better understanding of the development of ore deposits. Ore deposits usually occur along contacts between dissimilar materials. In the particular case presented here, a gold deposit was found along the contact between the basalt and ultramafic rock. Figure 20 shows structure contours of the top of the granite/ultramafic lithological contact whereas Fig. 21 shows the position of the Lucky and Golden faults. The analysis was performed to test a theory of geologists about what occurred in the region during the Precambrian period (the inferred age of the deposit is 830 million years.) The geologists are fairly certain that the region was shortened in an approximate east/west direction. It is assumed that the deposit was formed at approx. 5 km depth below the surface. On this basis it was suggested that uniform stress boundaries should be applied to a volume or rock of 2000 x 2000 x 1000 m which contains the lithological contact and the Lucky and Golden faults.
505
Lithological contact Fault zones
E = v= ~b = c= q5 = c=
75 GPa 0.25 3ff ~ 5 MPa 20' 0 MPa
The model was restrained on one surface and loaded on the opposing surface. The boundary stresses were as follows: a~ = tr2 = tr~ =
145 MPa 132 MPa 50 M P a
horizontal, east west horizontal, north south vertical
Figure 24 shows a view from the top of the displaced shape of the mesh after two iterations. Figure 25 shows a view from below the mesh. It can be clearly seen that slip occurs on the Lucky and Golden faults. Figure 26 shows the distribution of stress normal to the lithological contact after slip has occurred. Figure 27 shows contours of slip up/down dip on the Lucky and Golden faults. The values are in metres. Finally, Fig. 28 shows the contours of the slip on the lithologicai contact. The red area seems to be near the place where the actual deposit was found. SUMMARY AND CONCLUSIONS An efficient numerical model has been presented which is able to simulate large amounts of fault movement in a rock mass with very few iterations. The advantages of this numerical model are: (1) A boundary discretization only is needed to analyse the problem. This essentially reduces the dimensionality of the problem. For example, it would have been extremely difficult to analyse the Hellyer problem with a reasonable amount of human and computer resources using a domain type method, i.e. finite element method, finite difference method or distinct element method. An analysis of the second example problem presented in this paper had been attempted unsuccessfully using a finite difference code previously. (2) There is no need to specify stiffness coefficients for the fault planes. Very often these are difficult to obtain in the field and unless the fault has infill the model is not really required to be able to simulate elastic joint deformation in compression. (3) Very few iterations are needed to obtain large values of slip and separation. For example slip in the order of 1 m was modelled after only two iterations in the geological example.
506
BEER and POULSEN:
BEM MODELLING OF FAUI.TED R O ( K
(4) The model is very efficient in terms of computer resources as long as few faults are modelled. As the number of faults/joints increases the model becomes less efficient, and for a large number of faults a distinct element code such as 3DEC would probably be a more efficient model to apply. The model is applicable to a large range of problems involving joints/faults two of which have been presented here. It is also applicable to the modelling of the initiation and propagation of cracks along material interfaces [10].
Acknowledgements--The authors would like to thank Dan Sullivan for drafting many of the figures used in this paper. Thanks also to Aberfoyle Resources, Hellyer Division for permitting the presentation of the first example of the application. Greg Marshall and Kevin Rossengren have assisted in this analysis. The data for the second example have been supplied by Bill Power of the CSIRO Rock mechanics Research Center, Perth. Accepted for publication 9 November 1993.
REFERENCES I. Crouch S. L. and Starfield A, M. Boundary Eh'ment ,~,lcth,~L~ i, Solid Mechanics, Allen & Uwin, London (1983). 2. Cruse T. A. and Myers G. J. Three-dimensional fracture mechanics analysis. J. Struct. Dir, .Ira. Soc. ('it. En~,r~ 103, 309 320 (1977). 3. Gerstle W. H., lngraffea A. R. and Perucchio R. Threedimensional fatigue crack propagation analysis using the houndary element method, Int. J. Fatigue 10, 189 192 (1988). 4. Crotty J. M. and Wardle L. J. Boundary integral analysis of piecewise homogeneous media with structural discontinuities. Int. J. Rock Mech. Mh~. Sci. & Geomech Ahstr 22, 419 427 (1985). 5. Beer G. and Watson J. O. h~troduction to [~Tnite and Boundary Element MethodsJbr Engineers. Wiley, New York (1992). 6. Irons B. and Ahmad S. Techniques q/'Finite Elements. Horwood, Chitchester (1980). 7. Beer G. An efficient numerical method for modelling initiation and propagation of cracks along material interfaces. Int. J. Numer. Method Engng 36, 3579 3594 (1993). 8. BEFE--Users and Reference Manual. CSS International, Carnerigasse 10/45, Graz, Austria. 9. Beer G. and Mertz W. A user-friendly interface for computer aided analysis and design in mining, hu. J. Rock. Mech. Sci. & Geomech. Abstr. 27, 541-552 (1990). 10. Beer G., Bunker K. and Poulsen B. A. Modelling of crack propagation due to a mechanical rock excavator. In preparation.
Pergamon
0148-9062(93)E0030-R
Copyright '~i'; 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0148-9062/94 $7.00 + 0.00
Efficient Numerical Modelling of Faulted Rock Using the Boundary Element Method G. BEERt B. A. POULSEN~ This paper presents an efficient method of modelling faulted rock using the multiregion boundary element method. The boundary element method is ideally suited for modelling problems in rock mechanics because of the requirement to descretize surfaces only. This contrasts with the volume discretization needed for the finite element method, finite difference method and the distinct element method. The disadvantage of the boundary element method is that the rock mass is assumed to be isotropic and elastic. However, it is possible to model the non-linear behaviour on fault planes by using the concept of multiple regions and by connecting these regions via non-linear springs or Goodman type joint elements. The disadvantage of this method has been that joint stiffnesses have to be specified and that a great number of iterations are needed because the stiffness coefficient of a spring may change very rapidly from a very high value (for example when a joint is in compression) to a very low value (when a joint is in tension). This paper describes a novel method for the treatment of faults with the boundary element method which eliminates the need to specify joint stiffness parameters and which is able to model significant fault movement with very few iterations. The numerical scheme makes very efficient use of computer resources because the non-linear iterations only involve degrees of freedom at the interfaces. Examples of applications of the new model to practical problems in mining and geological engineering are given.
INTRODUCTION The boundary element method is ideally suited for analysing problems in rock mechanics because only surfaces have to be discretised into boundary elements. Figure 1, for example, shows the boundary element mesh used for analysing mine excavations. With the conventional boundary element method only elastic and homogeneous domains can be analysed. Practical problems in rock mechanics do not comply with these ideal conditions and this has been one of the main reasons why the method has not found wider acceptance in rock mechanics. In reality, the rock mass is non-homogeneous and elasto-viscoplastic and has joints and faults. In general, if non-linear material behaviour is to be modelled with boundary elements then a volume discretization is needed. However,
the influence of discrete joints and faults present in the rock mass can be considered with the boundary element method in two ways: (1) by using displacement discontinuity boundary elements. These are special boundary elements which have two opposing surfaces; (2) by using the multiregion concept in the conventional boundary element method. The first approach is explained in detail in [1]. The present paper will concentrate on the second approach. THE MULTIREGION BOUNDARY ELEMENT METHOD
This method of dealing with piecewise homogeneous domains and joints and faults has been used by various authors [2-4]. Basically the approach involves two or more boundary element regions which are connected to
tlnstitut fiir Baustatik, University of Technology, Graz, Austria. :[:CSIRO Division of Exploration and Mining, P.O. Box 883, Kenmore, Qld 4069, Australia. 485
486
BEER and POULSEN:
BEM MODELLING ( ) F I.AI.I.Tt!I) ROCK
rrMVUE /
,/
:
'~\
iI
!/, /
/
SCALE:
j
12.255
Fig. 1. Boundary element mesh of a mining excavation.
each other. Some kind of joint behaviour is then assumed at the interfaces between regions. A boundary element region is a region which is bounded by a boundary element mesh and, in the case of an infinite and semi-infinite region, by the surface of a
INFINITE
sphere or half sphere with a radius of infinity. Figure 2 shows an example of a finite, infinite and semi-infinite region. Each region can be assigned elastic material properties, E, v etc. For a region one may establish a relation
C) C3 FINITE
SEMI-INFINITE
Fig. 2, Boundary element regions in two dimensions.
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
487
Fig. 3. Boundary element regions with a fault at the interface.
(a)
(d)
I
kP t",, ~.
(b)
(e)
(c)
(f) Fig. 4. Methods 1 and 2 for non-linear interface behaviour.
488
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK
sT
ts
n
Fig. 5. Local coordinate system at the interface.
between the tractions (boundary stresses) and the displacements. Using Betti's law we may write
c(e)u(e)
- .Is T(P, O)u(O) dS
+ ;sU(P, Q)t(Q)dS
(1)
where e(P) is the free term, t(Q) and u(Q) are vectors containing the values of traction and displacement at boundary point Q and T(P, Q) and U(P, Q) are the fundamental solutions for the traction and displacement at point Q due to orthogonal unit loads at P in the coordinate directions. Using isoparametric boundary elements [5] the displacement and traction can be approximated by J
u(~)
= ~ N:(~)u/ j=l
Substituting equations (2) and (3) into (1) the following matrix expression relating the displacements {u} and tractions {t} at all the nodes of the boundary of a region can be obtained using the point collocation method [5]
(2)
[Al{u} = [Bl{t}
One may now combine two or more regions to analyse problems with joints or faults. Figure 3 shows an example of the combination of two semi-infinite regions to analyse an excavation with one fault. For this example we can rewrite equation (4) for region I
[A]' {u}' = [B]' {t}'
(5)
[A] '~{ul H = [B]" {t}"
(6)
and for region II
J
t(~) = ~ Nj(~)t/
(3)
l-I
where uj and tj are values of u and t at n o d e j and N:(~) are shape functions of the intrinsic coordinates ¢ and r/.
At the nodes of each region, except for the interface nodes, either tractions or the displacements are known. On the interface nodes both displacements and
i!
tf = - r . V 3
Its
>
-tnj
tan ej + cj
(a)
(4)
(b)
Fig. 6. Conditions for a joint slip satisfied 2-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (shear link broken) for next time step.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK ~rl q~; +
489
c3
IR1 = F s.
lsl
tR2 = F s.
ts2
i73
ts 1
(a)
(b)
Fig. 7. Conditions for a joint slip satisfied 3-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (shear link broken) for next time step. tractions are unknown. Let us rearrange equation (5) as follows:
[A f, -- Bi]t{::} ' = [n f, - Ai]'{:f, } !
(7)
Where u r and t r are vectors containing the displacements and tractions of nodes not on the interface and u i and f are corresponding vectors for the nodes on the interface. Assuming that all tractions t f (i.e. the loads applied to the free surface) are known we can write:
~fuf] '
[.$.] ~ ~ = {b} ! - [Ai]'{ui} ! ( ti )
(8)
where {b}'= [W]'{tr} '
and
f~ '
ej
faro "~I =~foj
+
[ui 1 K i {ui} '
{t i }1 = {ti0 }1 + [Ki]l{u i }!
{u'} !
=
{¢}!!
(] 1)
and provided no tractions are applied at the interface (9)
{t'}! + {f}" = 0
(12)
~
" "
"
(a)
(10)
For the case where the fault/joint does not slip or dilate the conditions on the interface are:
[.~.]'= [Af, - B i ] '.
By treating each column of [Ai]l as the right-hand side, the system of equations (8) can be solved to give
{u
The first part of the solution represents the displacement of the "free" nodes u f° and the traction of the interface nodes f0 assuming the displacements at the interface are zero. Each column of the second part of the solution represents the corresponding displacement and traction values for unit displacements of the interface nodes. Thus [K i] represents a "stiffness" matrix which relates the traction of the interface with the displacements at the interface. For BE region I we can write the following relationship between the interface displacements and interface tractions
R
__tn
(b)
Fig. 8. Conditions for a joint separation satisfied 2-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (all links broken) for next time step.
490
BEER and POULSEN:
BEM M O D E L L I N G OF FAUL'I-EI) RO('K
tsR1 = -tsl
ls 1
(a)
(b)
Fig. 9. Conditions for a joint separation satisfied 3-D problem: (a) shows tractions and state of rigid links at the end of time step 0 and (b) shows applied tractions and state of rigid links (all links broken) for next time step.
which yields the following system of equations for the displacements at the interface ([K']l + [Ki]II){u i } = {ti° }1 ---]- {t'°} u
(13)
The interface stresses f may be obtained by substituting the interface displacements u ~ into equation (10). The displacements at the free nodes u f can be obtained by substituting into equation (9). For the case where the joint slips or dilates equations (11) and (12) no longer hold and an iterative scheme has to be employed. This is discussed in detail in the following section. MODELLING OF JOINTS AND FAULTS There are basically two main approaches for the modelling of slip and dilation on fault planes: (1) By connecting the nodes at the interface by non-linear springs. Initially a high stiffness is assumed for the spring. This stiffness is reduced in the direction normal to the interface if the yield condition for joint opening is reached and in the direction tangential to the interface if the yield condition for joint slip is reached; (2) By connecting the nodes at the interface by rigid links in the direction perpendicular and tangential to the interface. If the yield condition in tension is reached then the link is broken. If the joint yields in shear then the link in the tangential direction is broken. Figure 4 illustrates the two approaches. In method 1 the pair of interface nodes is connected by springs. Initially the springs are assigned high values of stiffness in the normal direction, k,p and shear direction, k~ [Fig. 4(a)]. If the yield condition for joint slip is reached then the shear stiffness is reduced to a small but non-zero value k~ for any subsequent load step [Fig. 4(b)]. If the
yield condition for tension is reached then both the shear and normal stiffness are reduced [Fig. 4(c)]. In method 2 the nodes are connected by rigid links which prevent shear and separation [Fig. 4(d)]. If the yield condition in shear is reached the link which prevents slip is removed [Fig. 4(e)]. If the yield condition in tension has been reached then all links are broken [Fig. 4(f)]. The advantage of method 2 is two-fold: (1) there is no need to specify initial and residual values of joint stiffness, quantities which are not easily obtained in the field; (2) convergence of the method is significantly faster since we do not have to deal with large differences in stiffness between iterations. For these reasons method 2 has been used in the present work and this will be explained in more detail here.
n
Asl Fig. 10. Sign convention tbr joint displacement.
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D ROCK
491
'4
!
!
I
I
I
I
URS /
URS ..,
J
1.1111
ORS
ORS
QRS
\ \ -
\
ilCm:.
\ \ \
PLS $o¢q.
Ul
I; -
ors I
\ PLS
I I
~. _ _
\ HVS
~'Dll,
"~
440~
\
\
i STZ
x,,, -
STZ
\1
400 LLr~L
~°L I
t O
PLS I
I
[]
'~
!
on1
:
Ji
\~1
. |,,_,~ !_,~s ± : - . , _ _ , ~ : . ~ _ _ ~ _ ~ ,--..~
10950 N Og,dldJ~4k~l9~1
Fig. I I. Stope section, Hellyer mine, Australia.
492
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
LOCAL INTERFACE COORDINATE SYSTEM To be able to implement the method, the stresses at the interface have to be expressed in a local system of coordinates with one direction pointing normal to the joint and the other one or two directions being perpendicular to this direction (Fig. 5). The transformation from global traction/displacement ti, ui to local tractions/displacements t~, u~ at a node./is given by
u~ = T iuj
( ! 4)
t~ = Tjtj
(15)
Iij =
Us
X,=
Is
-n>,
u~, Us 2
t~=
.~),
~/:
'~"2 ~
"'2 ~
$2:
(19)
It'} = [T]{¢}' = [T]{f° }' + [T][Kq'[T]T{u' }'
{t' }' = {t o' }' + [K']'{u' }J
(21)
{t o' }' = [T] {f0},
(22)
[K' ]~ = [T] [K']' [T] T
(23)
where
The matrix [T] is given, for example, for four interface nodes, by
ni~
(17)
t,, ts z
08)
Ii
T2 0 0
T3
0
0
(24) T4
The transformation matrix Tj is computed with the outward normal n/to the node i. For three dimensional
ITE dDARY IENTS FAULT
91-400
(20)
or
For 3-D problems we have (Fig. 5)
u;=
.~/,
Equation (I0) can now be rewritten in terms of local components as:
Where for 2-D problems
{.n}
1l ll: ]
H ~
T,=
-D E
o 1 = 16.3MPa = 9.9MPa
:66~ \
~MPa
°
/
//ff -
b = 159 ° d=l°
< r T ~
//
X/
Fig. 12. Boundary element mesh of the Hellyer mine. Also shown is the in situ stress tensor.
~-~
BEER and FQULSEN:
BEM MODELLING
Fig. 13. Contours
OF FAULTED
of normal stress on Jack fault.
Fig. 14. Shear stress along strike.
ROCK
493
494
BEER
and POULSEN:
BEM MODELLING
Fig. 15. Shear stress up/down Fig. 16. Contours
of slip along
OF FAULTED
dip.
strike after two iterations
ROCK
BEER and POULSEN:
BEM MODELLING
OF FAULTED
Fig. 17. Contours of slip up/down after two iterations.
ROCK
495
496
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D R O C K
Fig. 18. Stress vectors at 360 level, no slip on Jack fault.
BEER and POULSEN:
BEM M O D E L L I N G OF F A U L T E D R O C K
problems vectors st and s2 are computed using the following vector x product (Fig. 5). sl = n x y
(25)
s2 = sl x n
(26)
where y is a unit vector in the y direction, i.e.
I f n points exactly in the y direction then equation (25) is replaced by s, = n x x
(28)
497
For smooth interfaces without kinks there is a unique value for nj at each node. For interfaces with kinks average values of n.j. can be assumed. If there are initial stresses tr0 present then we have to work out the corresponding initial tractions in the local interface coordinate system. To convert the initial stress to tractions at the interface we have t.,.0 = trx.onx + zx:..on,. + rx:.on-
(30)
t:,o = a:.on,. + zxy.onx + z:y.on..
(31)
lz, 0
(32)
=
tTz,on: --I- "Cxz.onx --l- Tyz,0ny
The tractions in the directions normal and tangential to the interface are given by
where x is a unit vector in the x direction t~ = Tt0
(33)
These initial interface tractions have to be added to the vector of traction due to a rigid interface {to'} ~.
[
10.000 11.gTg
Fig. 19. Stress vectors at 360 level, slip on Jack fault after two iterations.
lniiiai Discovery atx,y,z-0.0,0
Fig. 20. Structure
1
0
1
250
Fig. 21. Structure
contours
of Granite,
Ultramafic
contact
in plan view. Elevation
given in metres.
1
5OQm
contours
of the Lucky
and Golden 49x
faults in plan view. Elevation
given in metres.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK O N S E T OF SLIP AND SEPARATION--YIELD CONDITIONS
The condition for the onset of slip and separation (yield conditions) may now be written in terms of the values of tnj and ts,j, t~2/at interface node j as
Stip Fs = I tsjl + t,j tan q~j- CJ= 0
Separation
F,=t.j- ~=0
499
ITERATION (TIME STEP) PROCEDURE An iterative (time step) scheme has to be used for the solution of problems which involve slip and separation. For the first iteration (time step = 0) we assume rigid links between all the nodes on the interface. For the two region problem in Fig. 3 we can write
{u'},=o = {u'}' = {u' }"
(35)
and assuming no external forces are applied at the joint
{t'},=0 = {t'}l + {t'} t' = 0 (36) where the tension positive sign convention is used and ~bj is the angle of friction, cj is the cohesion and ~ is the This leads to the following system of equations tensile strength at node j. For three dimensional prob(37) ([K']' + [K']"){u'},:o = {t'°}' + {t'°} '' lems This system can be solved to give the displacements Itsjl - ~ lJ dr- t s2j 2 (34) {u' }, = 0 in the joint plane. The tractions at the joint plane
FEHVLIE
I
I / I
I •
I I
/ /
/ I
SCALE: I~1
Fig. 22. Boundary element mesh used for the geological modelling example. RMMS M ,5--H
201.507
500
BEER and POULSEN:
BEM MODELLING OF FAULTI:D RO(K
in the local (joint) coordinate system can be obtained by substituting {u'},=0 for {u'} ~ in equation (21)
{t'}',=0 = {t'°}~ + [Wl{u' I, =,,
(38)
where {t '° }~ includes the initial joint tractions mentioned above. Once we know these tractions we can check the yield conditions and determine if the joint slips or opens. The action taken for the next iteration depends on the value or F, or Ft- If both F, and F~ at an interface node are zero then no action is taken for that node. If F, is greater than or equal to, zero then the condition for joint slip is satisfied. If F, is greater than zero then at one stage during the increment the yield condition has been reached and the tractions at the end of the time step violate the yield condition. The component of the shear traction which is above F, = 0 is also called "excessive" shear traction. The procedure for joint yielding in shear is as follows (Figs 6 and 7): (1) the link(s) in the direction perpendicular to the vector normal to the joint plane are broken, (2) any "excessive" shear tractions are applied to the node in the opposite direction.
For yielding in tension the procedure is as tk~llows ( Figs 8 and 9): (I) all links are broken: (2) "excessive" shear and normal tractions are applied in the opposite direction. PROGRAMMING CONSIDERATIONS
The implementation of the algorithm for connecting and disconnecting degrees of freedom is explained here in more detail. The iteration time steps involve the displacements on joint planes only via the "region stiffness matrices" [Ki] ', [Ki] u, etc. The system of equations is solved using a frontal solution method [6]. In this method each variable is assigned a "destination" which determines where its stiffness coefficients are to be assembled. If two elements are connected by rigid links at a node then the stiffness coefficients for that node are assembled in the same location (that is they are added). If links are broken then the stiffness coefficients corresponding to the direction of the broken link (degree of freedom) are simply assembled in different locations. In the software implementation the allocation of "'destination" is handled as follows: each node is assigned a "restraint code" for each degree of freedom which is set to I for each connected degree of freedom
FEHVUE
2 / /
(a)
(b)
J
(c)
(d)
Fig. 23. The four boundary element regions of the geological modelling example.
BEER and POULSEN:
BEM MODELLING OF FAULTED ROCK
and to 2 for each disconnected degree of freedom. Initially all "restraint codes" are set to 1. If a link is broken then the "restraint code" for the corresponding df. is set to 2. JOINT OPENING/CLOSING, SLIP/STICK MONITOR During the iteration (time stepping) it is important to keep track of joints which are closing after being open and joints where the direction of slip is reversed. This is done by a device known as a opening/closing, slip/stick monitor. The joint accumulated relative displacements (i.e. slip and opening) at time step t at node j can be computed from the absolute displacements at the interface by
Opening Anj(t) = u,i(t " ) -- ui,j(t)
(39)
501
Slip I As,i( t )
= u'~',j(t)-
As2j(t)
=
u'~,j( t )
(40)
Slip 2
(41)
.','~j(t) - u~j(t)
Figure 10 shows the sign convention used for the relative joint displacement. The monitor is implemented in the software in such a way that a rigid link is reestablished for the case where the joint closes again or the slip is reversed. Joint closing is detected when Ajj(t) changes sign and slip reversal occurs if either Aslj(t) or As2/(t) changes sign. APPLICATIONS
The numerical method presented above has been implemented in program BEFE [8] and extensively tested. Comparisons were made with theoretical solutions or solutions obtained with other methods or
FEMVUE
DISP~. SCAkE: I
I
4.545
SCALE: I
Fig. 24. View from above, note the slip of Lucky and Golden faults.
I
201 .Mg
502
BEER and P O U L S E N :
BEM M O D E L L I N G O F F A U L T E D R O C K
software. Test examples are reported in [7] and will not be repeated here. Instead we will show applications in mining and geological engineering.
MINE EXCAVATIONS IN FAULTED ROCK
This application relates to the modelling of the excavation of ore by open stoping and caving at the Hellyer Mine in Tasmania, Australia. The region modelled is in the northern area of the mine. Figure 11 shows a typical section through the caved stope. The geological feature modelled is the Jack fault, which is sub-parallel to the caved stope and limits it in the upper part. Figure 12 shows a perspective view of the boundary element mesh for modelling the caved stope, the 87A stope, the 400-91 stope and the Jack fault. The mesh was prepared from horizontal and vertical sections using the preprocessor F E M C A D [9]. Four node boundary elements with a linear shpae function were used. Although
the Jack fault consists of two surfaces with boundary elements on each side of the fault, only one surface had to be specified with the program automatically generating the other side. Although there are about four different materials only one material was assumed lbr the purposes of the analysis. The elastic properties of this material are: E = 30,000 MPa v =0.3 These values were derived from uniaxiai compression tests on drill core. The properties of the Jack fault were assumed to be Cohesion = 0.0 MPa Angle of friction = 20 '-~ Dilation angle = O These properties were derived from observation of the fault surface where it is exposed, and from experience.
FEMVHE
J
DISPL. SCALE: 1 I 4.545
SCALE:
I
Fig. 25. View from below, note the slip of the faults.
I
164.141
BEER
and POULSEN:
Fig. 26. Normal
BEM MODELLING
stress distribution Fig. 28. Contours
on lithological
OF FAULTED
contact
of slip on lithological
(Fig. 2% Overleaf)
ROCK
after two iterations. contact.
503
504
BEER
and POULSEN:
Fig. 27. Contours
BEM MODELLING
of slip up/down
OF FAULTED
dip on Lucky
and Golden
ROCK
faults.
BEER and POULSEN: BEM MODELLING OF FAULTED ROCK Virgin stress measurement results from only one site were available and these have been used. The stress field was assumed to vary linearly from the surface and have the following values:
At surface:
all stress zero
At depth of 280m: a~ = 16.3 M P a a2 = 9.9 M P a tr3 = 3.8 MPa
Figure 22 shows the boundary element mesh used for the analysis. The mesh consists of four boundary element regions which are shown in Fig. 23. Note that the mesh does not have any internal elements, only surface elements thus making mesh generation easier. The following material properties are used: Rock mass
bearing = 159 ° bearing = 66 ° bearing = 249 °
dip = 1° dip = 71 ° dip = 19 °
These stress tensors are shown in Fig. 12. The following figures show some results of the analysis. Figure 13 shows contours of the stress normal to Jack fault. One can clearly see the destressed zone adjacent to the excavation and the stress concentration at the corners where the Jack fault surface is exposed. Figures 14 and 15 show contours of shear stress along strike and up dip. The contours of slip along strike after 2 iterations (time steps) are shown in Fig. 16. The m a x i m u m values of slip of 15 m m occurs to the south of the caved area. Figure 17 shows contours of slip up/down dip after two iterations. The maximum value of slip ( 1 4 m m ) occurs where the Jack fault surface is exposed to the caved area. The predicted directions and values of slip were found to agree well with observations although actual magnitudes were found to be larger than those computed. Figures 18 and 19 show the effect of fault movement on the stress distribution at the 360 level. Figure 18 shows the stress tensors before any movement on the fault occurs, whereas Fig. 19 shows the stress tensors after the Jack fault has slipped by 15 mm. The main effect is to turn the stress tensors so that they are more perpendicular to the fault in the area to the south of the caved zone. Consequently the 87 Pillar, which was in a stress shadow before, has become more highly stressed.
Geological modelling This application relates to gaining a better understanding of the development of ore deposits. Ore deposits usually occur along contacts between dissimilar materials. In the particular case presented here, a gold deposit was found along the contact between the basalt and ultramafic rock. Figure 20 shows structure contours of the top of the granite/ultramafic lithological contact whereas Fig. 21 shows the position of the Lucky and Golden faults. The analysis was performed to test a theory of geologists about what occurred in the region during the Precambrian period (the inferred age of the deposit is 830 million years.) The geologists are fairly certain that the region was shortened in an approximate east/west direction. It is assumed that the deposit was formed at approx. 5 km depth below the surface. On this basis it was suggested that uniform stress boundaries should be applied to a volume or rock of 2000 x 2000 x 1000 m which contains the lithological contact and the Lucky and Golden faults.
505
Lithological contact Fault zones
E = v= ~b = c= q5 = c=
75 GPa 0.25 3ff ~ 5 MPa 20' 0 MPa
The model was restrained on one surface and loaded on the opposing surface. The boundary stresses were as follows: a~ = tr2 = tr~ =
145 MPa 132 MPa 50 M P a
horizontal, east west horizontal, north south vertical
Figure 24 shows a view from the top of the displaced shape of the mesh after two iterations. Figure 25 shows a view from below the mesh. It can be clearly seen that slip occurs on the Lucky and Golden faults. Figure 26 shows the distribution of stress normal to the lithological contact after slip has occurred. Figure 27 shows contours of slip up/down dip on the Lucky and Golden faults. The values are in metres. Finally, Fig. 28 shows the contours of the slip on the lithologicai contact. The red area seems to be near the place where the actual deposit was found. SUMMARY AND CONCLUSIONS An efficient numerical model has been presented which is able to simulate large amounts of fault movement in a rock mass with very few iterations. The advantages of this numerical model are: (1) A boundary discretization only is needed to analyse the problem. This essentially reduces the dimensionality of the problem. For example, it would have been extremely difficult to analyse the Hellyer problem with a reasonable amount of human and computer resources using a domain type method, i.e. finite element method, finite difference method or distinct element method. An analysis of the second example problem presented in this paper had been attempted unsuccessfully using a finite difference code previously. (2) There is no need to specify stiffness coefficients for the fault planes. Very often these are difficult to obtain in the field and unless the fault has infill the model is not really required to be able to simulate elastic joint deformation in compression. (3) Very few iterations are needed to obtain large values of slip and separation. For example slip in the order of 1 m was modelled after only two iterations in the geological example.
506
BEER and POULSEN:
BEM MODELLING OF FAUI.TED R O ( K
(4) The model is very efficient in terms of computer resources as long as few faults are modelled. As the number of faults/joints increases the model becomes less efficient, and for a large number of faults a distinct element code such as 3DEC would probably be a more efficient model to apply. The model is applicable to a large range of problems involving joints/faults two of which have been presented here. It is also applicable to the modelling of the initiation and propagation of cracks along material interfaces [10].
Acknowledgements--The authors would like to thank Dan Sullivan for drafting many of the figures used in this paper. Thanks also to Aberfoyle Resources, Hellyer Division for permitting the presentation of the first example of the application. Greg Marshall and Kevin Rossengren have assisted in this analysis. The data for the second example have been supplied by Bill Power of the CSIRO Rock mechanics Research Center, Perth. Accepted for publication 9 November 1993.
REFERENCES I. Crouch S. L. and Starfield A, M. Boundary Eh'ment ,~,lcth,~L~ i, Solid Mechanics, Allen & Uwin, London (1983). 2. Cruse T. A. and Myers G. J. Three-dimensional fracture mechanics analysis. J. Struct. Dir, .Ira. Soc. ('it. En~,r~ 103, 309 320 (1977). 3. Gerstle W. H., lngraffea A. R. and Perucchio R. Threedimensional fatigue crack propagation analysis using the houndary element method, Int. J. Fatigue 10, 189 192 (1988). 4. Crotty J. M. and Wardle L. J. Boundary integral analysis of piecewise homogeneous media with structural discontinuities. Int. J. Rock Mech. Mh~. Sci. & Geomech Ahstr 22, 419 427 (1985). 5. Beer G. and Watson J. O. h~troduction to [~Tnite and Boundary Element MethodsJbr Engineers. Wiley, New York (1992). 6. Irons B. and Ahmad S. Techniques q/'Finite Elements. Horwood, Chitchester (1980). 7. Beer G. An efficient numerical method for modelling initiation and propagation of cracks along material interfaces. Int. J. Numer. Method Engng 36, 3579 3594 (1993). 8. BEFE--Users and Reference Manual. CSS International, Carnerigasse 10/45, Graz, Austria. 9. Beer G. and Mertz W. A user-friendly interface for computer aided analysis and design in mining, hu. J. Rock. Mech. Sci. & Geomech. Abstr. 27, 541-552 (1990). 10. Beer G., Bunker K. and Poulsen B. A. Modelling of crack propagation due to a mechanical rock excavator. In preparation.