- Email: [email protected]
Numerical simulation of progressive rock failure and associated seismicity
Int. J. Rock Mech. Min. Sci. Vol. 34, No. 2, pp. 249-261, 1997
Pergamon PII: S0145-9062(96)00039-3
© 1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0148-9062/97 $17.00 + 0.00
Numerical Simulation of Progressive Rock Failure and Associated Seismicity CHUN'AN TANGt In this paper, a numerical approach to modeling progressive failure leading to collapse in rock and associated seismicity is reported. In the first part, a newly developed numerical code, RFPA 2~ (Rock Failure Process Analysis), is introduced. The program allows modeling of the observed evolution of damage and associated seismic events due to progressive failure leading to collapse in brittle rock. There are three features distinguishing the approach from conventional numerical methods such as FEM: (1) by introducing heterogeneity of rock parameters into the model, RFPA eD can simulate non-linear behavior in rock using a linear method; (2) by introducing elastic modulus reduction for failed elements, RFPA 2~ can process discontinuum mechanics problems by a continuum mechanics method; and (3) by recording the event-rate of failed elements, the seismicities associated with the progressive failure in rock can be simulated. In the second part of the paper, the applications of RFPA eo in simulating geological process and in solving mining design problems are illustrated. Numerical simulation of a fault initiation process indicates that some of the important phenomena, such as coalescence of microfracture, the nucleation and growth of crack clusters, fault initiation and development, elastic rebound, dilatation, uplift and seismic behavior, etc. can be simulated with this numerical code. The simulation of progressive failure leading to collapse in underground openings demonstrates the capacity of this code in solving mining design problems. © 1997 Elsevier Science Ltd. All rights reserved.
demonstrates non-linear behavior. It also prevents us from approaching seismicity problems since the seismicity-prone material is by no means homogeneous. One of the most important factors affecting the progressive failure is heterogeneity. When rock is subjected to a stress field, cracks may nucleate, propagate, interact and coalesce. During fracturing, the heterogeneity plays a marked influence in determining the fracture paths and the resulting fracture patterns. The influence of heterogeneity is pronounced on the progressive failure process [4]. Even core specimens obtained from a seemingly homogeneous block of rock show variability both in deformation and strength properties. Thus, the distributive character of the heterogeneity plays a crucial role in determining the evolution of fractures. Therefore, a more reasonable numerical model for the rock or rock mass should be able to demonstrate the progressive failure due to heterogeneity, which results in non-linear behavior [5]. This may only succeed via a statistical approach. The recent work by Basista, Kracinovic, Blair and Cook, and Kim and Yao, etc. has shed new light on the understanding of the progressive
INTRODUCTION Non-linear and discontinuous numerical methods have become an important tool in modeling geological or geotechnical engineering processes for the analysis of stresses and deformations in geological or rock engineering structures such as faults, underground openings, slopes and foundations [1-3]. Although many numerical methods, such as finiteelement, boundary-element, finite-difference and discrete-element methods, can do well in simulating non-linear behavior in rock deformation, most of them are not a physical modeling of the non-linear behavior of brittle rock. In the non-linear FEM, for example, the non-linear behavior has been assigned by introducing a non-linear constitutive law to elements which are considered to be homogeneous material. This undoubtedly results in a non-linear outcome. Though this method has gained a sufficient degree of functionality, and may also incorporate compaction laws, it may not always help in our understanding of why rock material ~'Center for Rockbursts and Induced Seismicity Research, Northeastern University, Shenyang, 110006 PR China. 249
250
TANG:
SIMULATION OF PROGRESSIVE FAILURE
failure phenomenon and the damage mechanisms [4, 6, 7]. Some of these authors have developed or used "lattice" models, where random strengths and/or sizes are assigned to the individual members of a lattice. The overall response of the lattice due to both the elastic deformation and progressive breaking of individual members can be established. Even though the individual members in the lattice do not explicitly model actual micromechanics in rock, the overall response of the models is very "rock-like". However, the post-peak (strain-weakening) response is non-deterministic and represents a challenge to the field of damage mechanics [6]. Rock deformation is associated with very complicated progressive failures, as characterized by coalescence of microfracture, the nucleation and growth of cluster, fault initiation and development, dilatation, and seismic behavior, etc. which, as a global manifestation, demonstrates post-peak or strain-weakeningt behavior. It is known from laboratory experiments that near the peak stress, microfractures coalesce and bifurcate to form discrete macroscopic shear faults or splitting fractures. As pointed out by Kemeny [6], this level of sophistication has generally been lacking in most of the micromechanical damage models to date, and correctly incorporating these phenomena has been a challenge for numerical modelers in geotechnical engineering for the past several decades. In this paper, a new approach to modeling progressive failure in rock is reported. The numerical method, RFPA 2D (Rock Failure Process Analysis Code), developed recently by Tang [8], allows one to model the observed evolution of damage and induced seismicity due to the progressive failure leading to collapse in brittle rock or rock mass. A brief description of the development of RFPA 2Dis presented in the next section, followed by two examples of modeling rock failure as a geological process and in underground openings for a crown pillar collapse are presented, respectively.
include rock samples with the laboratory test set-up, geological structure, underground openings, pillars, slopes, etc. The program is used to construct models, analyze and then display the failure process with information including stress, strain, displacement and micro-seismic events. When using continuum mechanics and the linear method for numerically processing non-linear and discontinuum mechanics problems in rock failure, some problems are encountered, which require specialized approaches. These problems and the appropriate solutions are discussed in more detail in the following sub-sections.
Heterogeneity The most important hypothesis reflected in the RFPA 2° code is that heterogeneity in rock strength causes progressive failure behavior. Hence, strength heterogeneity is introduced into the model by defining an element strength (ao) distribution. In the current version of RFPA 2D, a Weibull distribution, defined by equation (1) is used, but other probability density functions can also be used depending on the user's understanding of the material property.
~- m(O'c~ m-I (a c)m e"°"
(1)
qg(a~) ao\ao ]
where the scale parameter o-0 is relative to all element strengths. The parameter m defines the shape of the density function which defines the degree .of material homogeneity, and is referred to as the homogeneity index (as shown in Fig. 1). The same density function is used to generate variation in elastic constants.
Non-linear behavior
RFPA 2Dis based on the idea that the heterogeneity is the source of non-linearity [8]. It will be shown in this paper that the global non-linear behavior observed in brittle rock can be simulated with brittle-elastic elements if the heterogeneity is considered. Therefore, it is not B R I E F D E S C R I P T I O N O F R F P A zD necessary to introduce more complicated constitutive The demand for new tools which may contribute to laws for every element, because the deformation improved understanding of seismicities induced in rock behavior of a single element contributes little to the failure, has initiated the development of the RFPA 2D macro-behavior of rock. Furthermore, the author code [8]. Since its initiation, significant progress has been made, especially in its capacity to deal with practical N problems including geological process and mining ~ m=1,5 250 induced seismicity. This section gives a brief introducm=3 tion and reports on the recent developments. 200 i ~ m=8 RFPA 2D, based on the linear finite element method, is 150 a comprehensive two-dimensional rock failure analysis package. The program allows simulation of the 100 ] progressive failure of rock leading to collapse via a 50 simple approximation, eliminating the numerical complexities of non-linear, discontinuum codes. Models can O,
1"At normal temperature, rock never becomes soft. The so-called strain-softening in rock is in fact strain-weakening due to progressive failure under loading.
0
2°°
°
60°
60°
Fig. 1. Distribution of element strength.
oo
TANG: SIMULATION OF PROGRESSIVE FAILURE
251
Fig. 2. Linear constitutive laws and strength distribution (homogeneity index m = 3).
considers that using a non-linear plastic constitutive law for describing the microfracture behavior of rock is not consistent with the fact that the deformation of brittle rock emits AE, which is considered to be caused by elastic energy release. This means that the A E sources must be in an elastic-brittle behavior. In R F P A 2D, linear constitutive laws have been introduced for all elements, and elements are given different strength and elastic constant parameters depending on the heterogeneity of the rock materials (Fig. 2). When the pre-described strength is reached, the element is considered to have lost most of its strength and released most of its elastic energy. This makes the model very simple but, as will be shown, successful in simulating non-linear rock behavior and the associated induced seismic events. The stress-strain curves for three homogeneity indices are shown in Fig. 3(a). It can be seen that materials with
higher m values represent more homogeneous materials, whereas those with lower m values are more heterogeneous. It is clear that the models predict non-linear stress-strain behavior similar to that observed in laboratory tests. This shows that macro-non-linear behavior of brittle rock can be modeled by heterogeneous elements with micro-linear (elastic-brittle) behavior, which is more reasonable than the conventional method using micro-non-linearity to model macro-non-linearity. Seismic events
Nowadays, one method of observing damage during rock deformation experiments is by monitoring the acoustic emissions (AE) or seismic events produced during deformation. I f we assume a unique association between a single AE event and micro-crack forming
200. 150, ~l 100, 50, O,
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
500 400 ¢,J 300
< 2o0 100
c(%) Fig. 3. Stress-strain and AE count-strain curves for materials with different homogeneity properties of element strengths. RMMS 34/2--F
252
TANG: SIMULATION OF PROGRESSIVE FAILURE
Fig. 4. AE source distribution modeled with RFPA2D.
event, then the possibility of using this technique to assess indirectly the damage evolution exists. Cox and Meredith [9] have analyzed catalogues of AE events recorded during compression tests in rock in terms of the information they give about the accumulated state of damage in a material, and combine this measured damage state with a model for the weakening behavior of cracked solids, which shows that reasonable predictions of the mechanical behavior are possible. Based on this background knowledge, it seems reasonable to assume that the number of seismic events or AE should be proportional to the number of failed elements. This may be questioned on the ground that there may be some fracture which is not dynamic or for other reasons does not release most of its elastic energy. However, as pointed out by Cox and Meredith [9], this is not a problem w i t h the materials assumed to be
elastic-brittle, since most of the energy in an element will be released as elastic energy when it fails. Thus, by recording the counts of failed elements, the seismicities associated with the progressive failure can be simulated in R F P A 20 that allows elements to fail when overstressed. Figure 3(b) shows the AE count-strain curves for three homogeneity indices. It is clear that there exists a good relation between stress-strain curves and AE count-strain curves. These results are consistent with findings presented by Cox and Meredith [9], and Khair et al. [10]. Also, the energy released in the failure of elements can be calculated approximately according to the element strengths, and then, the magnitude of the energy release can be calculated. Figure 4 shows an example of the locations of AE sources and their magnitude when a rock sample is loaded in compression. The white circles represent the seismicities of the current loading step and
Fig. 5. Different friction considerations between RFPA2Dand conventional codes.
TANG:
SIMULATION OF PROGRESSIVE FAILURE
Fig. 6. Numerical simulation of progressive failure leading to collapse.
253
TANG:
SIMULATION OF PROGRESSIVE F A I L U R E
the black circles record the seismicities of the previous steps.
Discontinuum When an element fails, a discontinuity has been introduced. As mentioned by Plischke et al. [11], when the conventional numerical simulation reaches such a state, the calculation will be terminated or the user should interrupt the calculation and modify the model by introducing a joint or special gap elements into the model. However, this kind of artificial discontinuity usually has a too simple geometry. An alternate method proposed here is to run the same type of meshed model, but to vary the assumed parameters for the failed elements describing its new mechanical properties. In this way, there will be no difference between an unfailed element and a failed element, except that the latter has a weak and soft mechanical property. When failed elements connect to form a fault, the fault will have an unpredictable and more realistic geometry depending on the failed elements connected to it. Thus, in R F P A 2D, the newly formed fault can be modeled by the self-connecting of failed neighboring elements during the progressive failure process, and no sliding interface or special "gap" element is needed. For the pre-existed crack or fault, RFPA m also uses the continuum method, i.e. elements are used. For a crack, the elements have a very weak and soft mechanical property like an "air element". For the fault, the gouge materials can be taken into account.
255
to collapse of shear fracture, in which fault nucleation and associated seismicity is demonstrated. The rock properties are chosen so as to simulate brittle failure behavior near the earth's surface. The scale parameter for strength, a0, is assumed to be 2 0 0 M P a and homogeneity index m -- 2; for elastic constants, E0 = 60 GPa and m = 30; and v = 0.25. The ratio of compression to tensile strength, s = a~/at, is assumed to be 3. The size of the mesh was 100 × 60 with 6000 elements. The simulations are limited to a two-dimensional problem in plane strain conditions. In order to study the influence of flaws on the initiation of a fault, a weaker ellipse area is assigned in which the mean strength of elements is 40% lower than the surrounding areas and the elastic constant is 10% higher. Once the mechanical properties are assigned to the elements, the model is compressed axially and horizontally in a displacement control mode to simulate a constant far-field displacement rate. The ultimate displacement is divided into 50 steps. As the numerical simulation progresses, the positions of elements change, and so do the stresses and strains. These dynamic changes caused by the increasing displacement at the boundary are analyzed to help understand the mechanisms of progressive failure development and induced seismicity. Simulation results are discussed in the following sections.
Progressive failure leading to collapse Friction In conventional numerical models, if a joint or fault is modeled, the friction coefficient must be taken into account, but the geometry of these flows is not "joint-like" or "fault-like" since their surfaces are not really rough-like, therefore, it is not a physical model of the joint or fault geometry [as shown in Fig. 6(a)]. In RFPA w, however, a joint or fault is modeled more physically by a self-induced mechanical closure and slide on a real geometrical basis (roughness), rather than calculating its friction force. As mentioned above, when failed elements connect to form a fault, the fault will have an unpredictable geometry, depending on the failed elements connected to it [as shown in Fig. 5(b)]. The unfailed elements along the two sides of the fault may close and touch each other at some points since the fault elements are very soft due to the lower elastic modulus. Then, the touching elements will resist the movement of elements towards each other, which is the real mechanism of friction. Thus, although no friction coefficient has been introduced to the model, RFPA 2D can model the real friction behavior on a more physical interlock basis. APPLICATION IN MODELING FAULT NUCLEATION The first presented application is based on a simple geological model to simulate progressive failure leading
Figure 6 shows the progressive failure process obtained in the simulation. The relative stresses are presented in these plots (Fig. 6), which are calculated by the following equation: o's_ a, O'c
O-c
Sa3
(2)
O'c
where as/at represents the severity of element or its proximity to failure. The brightness of the gray shading in the plots indicates the stress levels (high = white, low = black). The following stages are observed during the progressive failure process: (1) During the initial loading phase (low stress level), a few fractures are localized and relatively sparse (see steps 12-27 in Fig. 6). The stress-strain behavior is nearly elastic during this stage [Fig. 7(a)]; (2) Then, a population of randomly located, non-interacting fractures is observed. These sequences involve a limited number of isolated fractures because of the wide distribution of weak elements. The principal characteristic of this stage of deformation is diffuse, heterogeneous microfracturing of the entire rock mass. The stress-strain curve becomes non-linear in this stage; (3) As microfracture damage accumulates, fractures become clustered involving more elements, leading to fracture interaction. As more damage accumulates, a few
256
TANG: SIMULATION OF PROGRESSIVE FAILURE development of a process zone that advances along the fault by propagating in-plane (step 36). More elements fail at the advancing fault tip. It is interesting to find that, due to the propagating of the fault, the fracture development in the left area becomes sparse, while in the bottom right area the failure only occurs along the fault. It should be noted that the self-induced fault geometry is much more complex than artificially introduced faults in a conventional numerical method. It is also worth noting that during the whole failure process, stress redistribution plays an important role. Dilation, uplift and tilt
Many solids show an increase in volume during deformation. In rock, the dilatant expansion takes the form of new cracks which open up between and through the grains. Much research effort has gone into the study of microcrack dilatancy, i.e. the form of dilatant cracking that one sees in a small laboratory sample under compressive stress. It has been found that at about half the peak strength of rock, dilatancy begins to increase clearly. Although it is not measured quantitatively, the dilation of the rock mass during progressive failure can be identified from Fig. 6. As observed in laboratory fracture experiments, the simulation shows that dilation occurs within the stressed area surrounding an impending rupture zone and develops at an accelerating rate at the surface area, especially in the left upper surface area (steps 29-36). Like the lateral strain and volumetric strain accompanying the axial strain in the specimen loaded in a laboratory compression test, the vertical movements of ground surface, such as uplift (or dilatancy) and tilt, have often been reported in leveling surveys [12, 13]. Figure 6 shows realistic pictures about the dilation, uplift and tilt. If we take the model shown in this figure as a tectonic plate, then after the nucleation of the fault, it becomes two tectonic plates and forms a converging boundary where plates move toward each other. Due to the converging and the slipping along the fault, uplift Fig. 7. Stress-strain, AE count-strain and the.frequency-magnitude and tilt occur at the ground surface (steps 36-42). relation. Seismic events and source locations
large clusters emerge before main fracture nucleation (step 28). The stress-strain curve shows strain-weakening behavior, with a small stress drop [Fig. 7(a)]; (4) Eventually the number of failed elements increases drastically as the central area, where the elements have lower strengths than the surrounding areas, becomes mechanically unstable, i.e. the elements in the weaker zone suddenly collapse forming a fault (step 29). The newly formed fault is composed of a large number of small fractures. As a result, the fault produces its own stress field which dominates further fracture growth. At this stage, elastic rebound in the surrounding areas of fault (step 29 in Fig. 6) and a large stress drop [Fig. 7(a)] are observed; (5) The stress field of the fault concentrates around it, producing further microfracture damage around and ahead of its tips. This stress concentration leads to the
Two steps for AE source locations are presented in Fig. 8. The results related to AE behavior provide a number of important clues regarding the progressive failure leading to brittle faulting. The pre-nucleation of AE locations indicates that due to the heterogeneity in a rock mass, AE source locations remain diffuse before the main shocks [Fig. 8(a)]. During the main shock, a large number of seismic events occur along a narrow zone [Fig. 8(b)], i.e. the rock fails with a "self-organizing" mechanism--from "disorder to order". The AE behavior observed in these simulations is comparable to the experimental results described by Lockner et al. [14]. The curve of event rate with time is shown in Fig. 7(b), and the frequency-magnitude relation is also shown in Fig. 7(c). It is worth noting that every sudden increase of event rate in Fig. 7(b) results in a sudden stress drop in Fig. 8(a), as in the case of laboratory tests [9, 10]. This
TANG: SIMULATION OF PROGRESSIVE FAILURE
257
forepeak part of the relation, however, it shows a difference between the numerical model and the G u t e n b e r ~ R i c h t e r linear relation. This difference has been found by many seismologists in their observations or experiments [15]. APPLICATION IN MODELING PILLAR COLLAPSE AND ASSOCIATED SEISMICITY
Fig. 8. Seismic events and source location.
quantitative approach can be an improvement on the more traditional methods of AE analysis in correlating the AE activity with the weakening of the material. Although the relation between AE frequency and magnitude is not the focus of this paper, it is still noteworthy that the relation between the frequency of events and magnitude is approximately linear in the post-peak part of the curve shown in Fig. 7(c). This conclusion is significant because it does verify the empirical Gutenberg-Richter relation and in turn provides validity for our numerical model. For the
When designing a tunnel in a civil engineering project, elastic analysis is generally accepted as a tool to help in analyzing the stability of a permanent opening. However, in designing a mining underground opening, it is common for designers to be faced with distinctive problems related to non-linear behavior of the rock mass in which the stresses may reach magnitudes large enough to cause failure of rock surrounding the openings. Especially in burst prone ground, it is more difficult to obtain solutions to these problems than in non-burst prone ground. Since ground control problems in mining are always associated with progressive failure due to material heterogeneity, it follows that any attempt to predict ground conditions with conventional elastic analysis or even plastic analysis techniques without considering the influence of heterogeneity is useful if only to indicate potential failure zones or, at best, to predict yielding caused by geometrical reasons. In this section, progressive failure leading to a crown pillar collapse is simulated. The model set up is illustrated in Fig. 9. Although the geometry of this model comes from the INCO's Creighton Mine, Sudbury, Canada [16], it is not intended in this paper to model exactly according to the situation in the mine. The purpose of this example is for demonstrating the potential of RFPA 2Din solving mining design problems. There is a fault located along the footwall contact. The upper excavation represents the current mining extending from hangingwall to footwall. Between the upper
Fig. 9. Numerical model of a cross-section perpendicular to a strike.
258
TANG: SIMULATION OF PROGRESSIVE FAILURE
excavation and the lower excavation is the crown pillar. Assume that due to the mining process surrounding the modeled area there is a large horizontal far-field displacement rate compared with the vertical far-field displacement. Therefore, the vertical confined boundary condition is selected for the simulation. The horizontal boundary is compressed in a displacement control mode. The rock properties are chosen so as to simulate a brittle rock. The scale parameter for strength, a0, is assumed to be 60 MPa and homogeneity index m = 4; for elastic constants, E0 = 65 GPa; and v = 0.25. The ratio of compression to tensile strength, s = ~/a~, is assumed to be 3. Due to the lack of material properties for the fault gouge, the strength and elastic constant of
the gouge material are arbitrarily chosen as 90% lower than the strength and elastic constant of rock, respectively. There is no difference for elements between the fault and the surrounding rock except that lower strength and elastic constants are set for the fault. The size of the mesh was 100 x 70 with 7000 elements. The simulations are limited to a two-dimensional problem in plane strain conditions. The simulated result is shown in Fig. 10. The seismic event rate for comparison is also given in Fig. 10. It can be observed that, although mathematically R F P A 2D is completely a simple static continuum elastic code, the outcome of this numerical model shows that the code has the ability to reveal the evolutionary nature of the
ii~
~W
...............~;~i~Nii~ii~
ii!!ii~i~;: ............
Fig. lO(a).
7
....
~
~¸~
° ~ ~.~N~
~
~~ ' ~ ~ ' ~..................... ; ~ ~ ~iii!~i iil
TANG: SIMULATION OF PROGRESSIVE FAILURE ........
....
......
............. ....
Step:
-~
i
II
259
5:
: : : :::'~
,
;_;
'
Fig. 10. Numerical simulation of progressive faiiure leading to collapse of crown pillar and associated seismicity. fracture phenomenon from microfracture scale to global failure. It can be seen that the whole failure process can be divided into four main stages: (1) The fault slips over most of its length with the bottom end connected to the upper left corner of the bottom excavation, inducing many seismic events along the fault. Some seismic events are also detected in the central area of the crown pillar. This corresponds to steps 1 - 3 i n Fig. 10; (2) A large triangular block forms in the roof of the upper excavation (step 4), and the block attempts to move into the opening (steps 5-50). In the bottom part of the excavation, however, more fractures develop in the floor leading to dilation of the ground surface, and
then upheaval phenomena are observed in the following steps (steps 17 50); (3) The fault slipping becomes drastic resulting in the sudden bending fracture of the crown pillar in the end near the hangingwaU (step 17), and then the crown pillar becomes one long block rotated in the manner shown in the following steps. The modeling of this example again demonstrates the advantage of RFPA 2D which can model the progressive failure process step by step, which is important in the numerical simulation of failure mechanisms in rock materials since progressive failure will result in a stress redistribution which may be very different from the original stress distribution.
260
TANG: SIMULATION OF PROGRESSIVE FAILURE
Fig. 11. Numerical simulation of crown pillar without supports, with weakened supports and with strong supports.
Another advantage of RFPA 2D is that it can be used to help examine the feasibility of a support design of underground openings in burst prone ground. For example, from the above simulation, three areas are found to be the locations with intense micro-seismic activity (as shown in Fig. 10). In order to find the way to prevent these failures, support methods are modeled by installing supports in those event-prone locations to keep the crown pillar stable. The simulation result for three supports installed in the upper and lower parts of
the excavation is shown in Fig. 11 (b). It is found that the supports are not strong enough and failure occurs again. Even though, the fault slip is prevented. Figure 1 l(c) shows a result for stronger supports in which the roof and floor in both excavations are confined completely. Although failure still occurs in the roof of the upper excavation and in the floor of the lower excavation, the crown pillar becomes less failed .and is stable during the whole loading process, as shown in Fig. 1 l(c). Comparing with the weaker supports shown in
TANG: SIMULATION OF PROGRESSIVE FAILURE Fig. 1 l(b), it is o b s e r v e d t h a t the stability function o f the s u p p o r t system p l a y s an i m p o r t a n t role as s h o w n in Fig. 1 l(c). CONCLUSION M o s t o f the n u m e r i c a l s i m u l a t i o n o f r o c k failure has involved stress a n d strain d e t e r m i n a t i o n s . Very little w o r k has been d o n e in n u m e r i c a l m o d e l i n g o f seismic activities i n d u c e d in progressive failure leading to collapse. Since the stability c o n t r o l p r o b l e m s in m i n i n g are often associated with r o c k mass d a m a g e a n d seismic activities, it follows t h a t a n y a t t e m p t to predict u n s t a b l e failure b a s e d on stress o r strain d e t e r m i n a t i o n s m u s t be limited to where the m e a s u r e m e n t s are taken. It is seen in this p a p e r t h a t the a m o u n t o f d a m a g e occurring to the r o c k m a s s d u r i n g failure, a n d hence the possibility o f r u p t u r e , can be p r e d i c t e d better b y m e t h o d s which a c c o m m o d a t e the progressive failure a n d a s s o c i a t e d seismic activities. T h e c o m b i n e d use o f failure event m o d e l i n g a n d event source l o c a t i o n d e t e r m i n a t i o n s (or as an a p p r o x i m a t i o n , AE), in c o n j u n c t i o n with s t r e s s - s t r a i n simulation, will p r o v i d e a d d i t i o n a l insight into the p r o b l e m . T h e successful a p p l i c a t i o n s in this p a p e r have s h o w n t h a t the R F P A 2D is in principle c a p a b l e o f m o d e l i n g geological processes a n d r o c k engineering p r o b l e m s . The r e p r o d u c i b i l i t y o f p h e n o m e n a , such as fault initiation a n d d e v e l o p m e n t , elastic r e b o u n d , uplift o r seismic activities, etc. offers an attractive s u p p l e m e n t to physical m o d e l tests in b o t h geological a n d g e o m e c h a n i c a l p r o b l e m s . This refers especially to m i n i n g in b u r s t - p r o n e g r o u n d , tunnel projects a n d the u n d e r g r o u n d storage o f nuclear waste at great depth. C u r r e n t d e v e l o p m e n t efforts focus on a t h r e e - d i m e n sional model. Acknowledgements--The work was funded by China National Natural
Science Foundation and China National Education Commission. This work was initiated in 1985 under the supervision of Professor X. H. Xu. The academic visit to Professor J. A. Hudson in 1991-1992 resulted in the formation of the theoretical frame work. The discussion with Professor C. Fairhurst in Kingston, Canada, 1993, gave me important encouragement to continue my work in rock failure
261
instability, and the recent visit to Professor P. K. Kaiser made me feel specially interested in the problem of mining induced seismicity. The author would like to take this opportunity to give my thanks to all of them.
Accepted for publication 21 June 1996.
REFERENCES 1. Sterpi D., Cividini A. and Donelli M. Numerical analysis of a shallow excavation in strain softening rock. Proceedings o f 8th Int. Congress on Rock Mechanics, pp. 545-549. Tokyo (1995). 2. Potyondy D., Cundall P. and Lee C. Modeling of notch formation in the URL mine-by tunnel: Phase I--Micromechanical experiments, Report to URL by Itasca (1995). 3. Shi G. H. Block system modeling by discontinuous deformation analysis. Southampton U.K. and Boston U.S.A. (1993). 4. Kim Kunsoo and Yao C. Y. Effects of micromechanical property variation on fracture processes in simple tension. In Rock Mechanics (Edited by Daemen and Schultz), pp. 471-476. Balkema, Rotterdam (1995). 5. Tang C. A., Hudson J. A. and Xu X. H. Rock failure instability and related aspects of earthquake mechanisms. China Coal Industry Publishing House, Beijing (1993). 6. Kemeny J. M. Meeting reports. Geotech. Geological Engng 9, 83-85 (1991). 7. Blair S. C. and Cook N. G. W. Statistic model for rock fracture. In Rock Mechanics (Edited by Tillerson and Wawersik), pp. 729-735. Balkema, Rotterdam (1992). 8. Tang C. A. Numerical simulation of rock failure process. Proceedings of the 2nd Youth Symposium on Rock Mechanics in China, Chengdu (1995).
9. Cox S. J. D. and Meredith P. G. Microcrack formation and material softening in rock measured by monitoring AE. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 30, 11-24 (1993). 10. Khair A. W. et al. Acoustic emission pattern. Third Conference on Acoustic Emission/Microseismic Activity in Geological Structures and Materials. The Pennsylvania State University, U.S.A. (1981).
11. Plischke B., Dallmann R. and Chedmail J. F. Numerical simulation of geological processes. Proceedings of the 7th Int. Congress on Rock Mechanics, pp. 1909-1912. Paris (1991). 12. Rikitake T. Earthquake Forecasting and Warning. Centre for Academic Publications Japan/Tokyo (1982). 13. Mogi K. Earthquake Prediction. Academic Press, Tokyo (1985). 14. Lockner D. A., Moore D. E. and Reches Z. E. Microcrack interaction leading to shear fracture. In Rock Mechanics (Edited by Tillerson and Wawersike), pp. 807-815. Balkema, Rotterdam (1992). 15. Latousakis J. and Drakopoulos J. A modified formula for frequency-magnitude distribution. PAGEOPH. 125, 753-764 (1987). 16. Wiles T. D. Non-linear modeling of mining at great depth: a case history at Creighton mine. Int. Symposium on Rock Mechanics and Rock Physics at Great Depth. Pau, France (1989).
Pergamon PII: S0145-9062(96)00039-3
© 1997 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0148-9062/97 $17.00 + 0.00
Numerical Simulation of Progressive Rock Failure and Associated Seismicity CHUN'AN TANGt In this paper, a numerical approach to modeling progressive failure leading to collapse in rock and associated seismicity is reported. In the first part, a newly developed numerical code, RFPA 2~ (Rock Failure Process Analysis), is introduced. The program allows modeling of the observed evolution of damage and associated seismic events due to progressive failure leading to collapse in brittle rock. There are three features distinguishing the approach from conventional numerical methods such as FEM: (1) by introducing heterogeneity of rock parameters into the model, RFPA eD can simulate non-linear behavior in rock using a linear method; (2) by introducing elastic modulus reduction for failed elements, RFPA 2~ can process discontinuum mechanics problems by a continuum mechanics method; and (3) by recording the event-rate of failed elements, the seismicities associated with the progressive failure in rock can be simulated. In the second part of the paper, the applications of RFPA eo in simulating geological process and in solving mining design problems are illustrated. Numerical simulation of a fault initiation process indicates that some of the important phenomena, such as coalescence of microfracture, the nucleation and growth of crack clusters, fault initiation and development, elastic rebound, dilatation, uplift and seismic behavior, etc. can be simulated with this numerical code. The simulation of progressive failure leading to collapse in underground openings demonstrates the capacity of this code in solving mining design problems. © 1997 Elsevier Science Ltd. All rights reserved.
demonstrates non-linear behavior. It also prevents us from approaching seismicity problems since the seismicity-prone material is by no means homogeneous. One of the most important factors affecting the progressive failure is heterogeneity. When rock is subjected to a stress field, cracks may nucleate, propagate, interact and coalesce. During fracturing, the heterogeneity plays a marked influence in determining the fracture paths and the resulting fracture patterns. The influence of heterogeneity is pronounced on the progressive failure process [4]. Even core specimens obtained from a seemingly homogeneous block of rock show variability both in deformation and strength properties. Thus, the distributive character of the heterogeneity plays a crucial role in determining the evolution of fractures. Therefore, a more reasonable numerical model for the rock or rock mass should be able to demonstrate the progressive failure due to heterogeneity, which results in non-linear behavior [5]. This may only succeed via a statistical approach. The recent work by Basista, Kracinovic, Blair and Cook, and Kim and Yao, etc. has shed new light on the understanding of the progressive
INTRODUCTION Non-linear and discontinuous numerical methods have become an important tool in modeling geological or geotechnical engineering processes for the analysis of stresses and deformations in geological or rock engineering structures such as faults, underground openings, slopes and foundations [1-3]. Although many numerical methods, such as finiteelement, boundary-element, finite-difference and discrete-element methods, can do well in simulating non-linear behavior in rock deformation, most of them are not a physical modeling of the non-linear behavior of brittle rock. In the non-linear FEM, for example, the non-linear behavior has been assigned by introducing a non-linear constitutive law to elements which are considered to be homogeneous material. This undoubtedly results in a non-linear outcome. Though this method has gained a sufficient degree of functionality, and may also incorporate compaction laws, it may not always help in our understanding of why rock material ~'Center for Rockbursts and Induced Seismicity Research, Northeastern University, Shenyang, 110006 PR China. 249
250
TANG:
SIMULATION OF PROGRESSIVE FAILURE
failure phenomenon and the damage mechanisms [4, 6, 7]. Some of these authors have developed or used "lattice" models, where random strengths and/or sizes are assigned to the individual members of a lattice. The overall response of the lattice due to both the elastic deformation and progressive breaking of individual members can be established. Even though the individual members in the lattice do not explicitly model actual micromechanics in rock, the overall response of the models is very "rock-like". However, the post-peak (strain-weakening) response is non-deterministic and represents a challenge to the field of damage mechanics [6]. Rock deformation is associated with very complicated progressive failures, as characterized by coalescence of microfracture, the nucleation and growth of cluster, fault initiation and development, dilatation, and seismic behavior, etc. which, as a global manifestation, demonstrates post-peak or strain-weakeningt behavior. It is known from laboratory experiments that near the peak stress, microfractures coalesce and bifurcate to form discrete macroscopic shear faults or splitting fractures. As pointed out by Kemeny [6], this level of sophistication has generally been lacking in most of the micromechanical damage models to date, and correctly incorporating these phenomena has been a challenge for numerical modelers in geotechnical engineering for the past several decades. In this paper, a new approach to modeling progressive failure in rock is reported. The numerical method, RFPA 2D (Rock Failure Process Analysis Code), developed recently by Tang [8], allows one to model the observed evolution of damage and induced seismicity due to the progressive failure leading to collapse in brittle rock or rock mass. A brief description of the development of RFPA 2Dis presented in the next section, followed by two examples of modeling rock failure as a geological process and in underground openings for a crown pillar collapse are presented, respectively.
include rock samples with the laboratory test set-up, geological structure, underground openings, pillars, slopes, etc. The program is used to construct models, analyze and then display the failure process with information including stress, strain, displacement and micro-seismic events. When using continuum mechanics and the linear method for numerically processing non-linear and discontinuum mechanics problems in rock failure, some problems are encountered, which require specialized approaches. These problems and the appropriate solutions are discussed in more detail in the following sub-sections.
Heterogeneity The most important hypothesis reflected in the RFPA 2° code is that heterogeneity in rock strength causes progressive failure behavior. Hence, strength heterogeneity is introduced into the model by defining an element strength (ao) distribution. In the current version of RFPA 2D, a Weibull distribution, defined by equation (1) is used, but other probability density functions can also be used depending on the user's understanding of the material property.
~- m(O'c~ m-I (a c)m e"°"
(1)
qg(a~) ao\ao ]
where the scale parameter o-0 is relative to all element strengths. The parameter m defines the shape of the density function which defines the degree .of material homogeneity, and is referred to as the homogeneity index (as shown in Fig. 1). The same density function is used to generate variation in elastic constants.
Non-linear behavior
RFPA 2Dis based on the idea that the heterogeneity is the source of non-linearity [8]. It will be shown in this paper that the global non-linear behavior observed in brittle rock can be simulated with brittle-elastic elements if the heterogeneity is considered. Therefore, it is not B R I E F D E S C R I P T I O N O F R F P A zD necessary to introduce more complicated constitutive The demand for new tools which may contribute to laws for every element, because the deformation improved understanding of seismicities induced in rock behavior of a single element contributes little to the failure, has initiated the development of the RFPA 2D macro-behavior of rock. Furthermore, the author code [8]. Since its initiation, significant progress has been made, especially in its capacity to deal with practical N problems including geological process and mining ~ m=1,5 250 induced seismicity. This section gives a brief introducm=3 tion and reports on the recent developments. 200 i ~ m=8 RFPA 2D, based on the linear finite element method, is 150 a comprehensive two-dimensional rock failure analysis package. The program allows simulation of the 100 ] progressive failure of rock leading to collapse via a 50 simple approximation, eliminating the numerical complexities of non-linear, discontinuum codes. Models can O,
1"At normal temperature, rock never becomes soft. The so-called strain-softening in rock is in fact strain-weakening due to progressive failure under loading.
0
2°°
°
60°
60°
Fig. 1. Distribution of element strength.
oo
TANG: SIMULATION OF PROGRESSIVE FAILURE
251
Fig. 2. Linear constitutive laws and strength distribution (homogeneity index m = 3).
considers that using a non-linear plastic constitutive law for describing the microfracture behavior of rock is not consistent with the fact that the deformation of brittle rock emits AE, which is considered to be caused by elastic energy release. This means that the A E sources must be in an elastic-brittle behavior. In R F P A 2D, linear constitutive laws have been introduced for all elements, and elements are given different strength and elastic constant parameters depending on the heterogeneity of the rock materials (Fig. 2). When the pre-described strength is reached, the element is considered to have lost most of its strength and released most of its elastic energy. This makes the model very simple but, as will be shown, successful in simulating non-linear rock behavior and the associated induced seismic events. The stress-strain curves for three homogeneity indices are shown in Fig. 3(a). It can be seen that materials with
higher m values represent more homogeneous materials, whereas those with lower m values are more heterogeneous. It is clear that the models predict non-linear stress-strain behavior similar to that observed in laboratory tests. This shows that macro-non-linear behavior of brittle rock can be modeled by heterogeneous elements with micro-linear (elastic-brittle) behavior, which is more reasonable than the conventional method using micro-non-linearity to model macro-non-linearity. Seismic events
Nowadays, one method of observing damage during rock deformation experiments is by monitoring the acoustic emissions (AE) or seismic events produced during deformation. I f we assume a unique association between a single AE event and micro-crack forming
200. 150, ~l 100, 50, O,
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
0.1
0.2
0.3
0.4
0.5
500 400 ¢,J 300
< 2o0 100
c(%) Fig. 3. Stress-strain and AE count-strain curves for materials with different homogeneity properties of element strengths. RMMS 34/2--F
252
TANG: SIMULATION OF PROGRESSIVE FAILURE
Fig. 4. AE source distribution modeled with RFPA2D.
event, then the possibility of using this technique to assess indirectly the damage evolution exists. Cox and Meredith [9] have analyzed catalogues of AE events recorded during compression tests in rock in terms of the information they give about the accumulated state of damage in a material, and combine this measured damage state with a model for the weakening behavior of cracked solids, which shows that reasonable predictions of the mechanical behavior are possible. Based on this background knowledge, it seems reasonable to assume that the number of seismic events or AE should be proportional to the number of failed elements. This may be questioned on the ground that there may be some fracture which is not dynamic or for other reasons does not release most of its elastic energy. However, as pointed out by Cox and Meredith [9], this is not a problem w i t h the materials assumed to be
elastic-brittle, since most of the energy in an element will be released as elastic energy when it fails. Thus, by recording the counts of failed elements, the seismicities associated with the progressive failure can be simulated in R F P A 20 that allows elements to fail when overstressed. Figure 3(b) shows the AE count-strain curves for three homogeneity indices. It is clear that there exists a good relation between stress-strain curves and AE count-strain curves. These results are consistent with findings presented by Cox and Meredith [9], and Khair et al. [10]. Also, the energy released in the failure of elements can be calculated approximately according to the element strengths, and then, the magnitude of the energy release can be calculated. Figure 4 shows an example of the locations of AE sources and their magnitude when a rock sample is loaded in compression. The white circles represent the seismicities of the current loading step and
Fig. 5. Different friction considerations between RFPA2Dand conventional codes.
TANG:
SIMULATION OF PROGRESSIVE FAILURE
Fig. 6. Numerical simulation of progressive failure leading to collapse.
253
TANG:
SIMULATION OF PROGRESSIVE F A I L U R E
the black circles record the seismicities of the previous steps.
Discontinuum When an element fails, a discontinuity has been introduced. As mentioned by Plischke et al. [11], when the conventional numerical simulation reaches such a state, the calculation will be terminated or the user should interrupt the calculation and modify the model by introducing a joint or special gap elements into the model. However, this kind of artificial discontinuity usually has a too simple geometry. An alternate method proposed here is to run the same type of meshed model, but to vary the assumed parameters for the failed elements describing its new mechanical properties. In this way, there will be no difference between an unfailed element and a failed element, except that the latter has a weak and soft mechanical property. When failed elements connect to form a fault, the fault will have an unpredictable and more realistic geometry depending on the failed elements connected to it. Thus, in R F P A 2D, the newly formed fault can be modeled by the self-connecting of failed neighboring elements during the progressive failure process, and no sliding interface or special "gap" element is needed. For the pre-existed crack or fault, RFPA m also uses the continuum method, i.e. elements are used. For a crack, the elements have a very weak and soft mechanical property like an "air element". For the fault, the gouge materials can be taken into account.
255
to collapse of shear fracture, in which fault nucleation and associated seismicity is demonstrated. The rock properties are chosen so as to simulate brittle failure behavior near the earth's surface. The scale parameter for strength, a0, is assumed to be 2 0 0 M P a and homogeneity index m -- 2; for elastic constants, E0 = 60 GPa and m = 30; and v = 0.25. The ratio of compression to tensile strength, s = a~/at, is assumed to be 3. The size of the mesh was 100 × 60 with 6000 elements. The simulations are limited to a two-dimensional problem in plane strain conditions. In order to study the influence of flaws on the initiation of a fault, a weaker ellipse area is assigned in which the mean strength of elements is 40% lower than the surrounding areas and the elastic constant is 10% higher. Once the mechanical properties are assigned to the elements, the model is compressed axially and horizontally in a displacement control mode to simulate a constant far-field displacement rate. The ultimate displacement is divided into 50 steps. As the numerical simulation progresses, the positions of elements change, and so do the stresses and strains. These dynamic changes caused by the increasing displacement at the boundary are analyzed to help understand the mechanisms of progressive failure development and induced seismicity. Simulation results are discussed in the following sections.
Progressive failure leading to collapse Friction In conventional numerical models, if a joint or fault is modeled, the friction coefficient must be taken into account, but the geometry of these flows is not "joint-like" or "fault-like" since their surfaces are not really rough-like, therefore, it is not a physical model of the joint or fault geometry [as shown in Fig. 6(a)]. In RFPA w, however, a joint or fault is modeled more physically by a self-induced mechanical closure and slide on a real geometrical basis (roughness), rather than calculating its friction force. As mentioned above, when failed elements connect to form a fault, the fault will have an unpredictable geometry, depending on the failed elements connected to it [as shown in Fig. 5(b)]. The unfailed elements along the two sides of the fault may close and touch each other at some points since the fault elements are very soft due to the lower elastic modulus. Then, the touching elements will resist the movement of elements towards each other, which is the real mechanism of friction. Thus, although no friction coefficient has been introduced to the model, RFPA 2D can model the real friction behavior on a more physical interlock basis. APPLICATION IN MODELING FAULT NUCLEATION The first presented application is based on a simple geological model to simulate progressive failure leading
Figure 6 shows the progressive failure process obtained in the simulation. The relative stresses are presented in these plots (Fig. 6), which are calculated by the following equation: o's_ a, O'c
O-c
Sa3
(2)
O'c
where as/at represents the severity of element or its proximity to failure. The brightness of the gray shading in the plots indicates the stress levels (high = white, low = black). The following stages are observed during the progressive failure process: (1) During the initial loading phase (low stress level), a few fractures are localized and relatively sparse (see steps 12-27 in Fig. 6). The stress-strain behavior is nearly elastic during this stage [Fig. 7(a)]; (2) Then, a population of randomly located, non-interacting fractures is observed. These sequences involve a limited number of isolated fractures because of the wide distribution of weak elements. The principal characteristic of this stage of deformation is diffuse, heterogeneous microfracturing of the entire rock mass. The stress-strain curve becomes non-linear in this stage; (3) As microfracture damage accumulates, fractures become clustered involving more elements, leading to fracture interaction. As more damage accumulates, a few
256
TANG: SIMULATION OF PROGRESSIVE FAILURE development of a process zone that advances along the fault by propagating in-plane (step 36). More elements fail at the advancing fault tip. It is interesting to find that, due to the propagating of the fault, the fracture development in the left area becomes sparse, while in the bottom right area the failure only occurs along the fault. It should be noted that the self-induced fault geometry is much more complex than artificially introduced faults in a conventional numerical method. It is also worth noting that during the whole failure process, stress redistribution plays an important role. Dilation, uplift and tilt
Many solids show an increase in volume during deformation. In rock, the dilatant expansion takes the form of new cracks which open up between and through the grains. Much research effort has gone into the study of microcrack dilatancy, i.e. the form of dilatant cracking that one sees in a small laboratory sample under compressive stress. It has been found that at about half the peak strength of rock, dilatancy begins to increase clearly. Although it is not measured quantitatively, the dilation of the rock mass during progressive failure can be identified from Fig. 6. As observed in laboratory fracture experiments, the simulation shows that dilation occurs within the stressed area surrounding an impending rupture zone and develops at an accelerating rate at the surface area, especially in the left upper surface area (steps 29-36). Like the lateral strain and volumetric strain accompanying the axial strain in the specimen loaded in a laboratory compression test, the vertical movements of ground surface, such as uplift (or dilatancy) and tilt, have often been reported in leveling surveys [12, 13]. Figure 6 shows realistic pictures about the dilation, uplift and tilt. If we take the model shown in this figure as a tectonic plate, then after the nucleation of the fault, it becomes two tectonic plates and forms a converging boundary where plates move toward each other. Due to the converging and the slipping along the fault, uplift Fig. 7. Stress-strain, AE count-strain and the.frequency-magnitude and tilt occur at the ground surface (steps 36-42). relation. Seismic events and source locations
large clusters emerge before main fracture nucleation (step 28). The stress-strain curve shows strain-weakening behavior, with a small stress drop [Fig. 7(a)]; (4) Eventually the number of failed elements increases drastically as the central area, where the elements have lower strengths than the surrounding areas, becomes mechanically unstable, i.e. the elements in the weaker zone suddenly collapse forming a fault (step 29). The newly formed fault is composed of a large number of small fractures. As a result, the fault produces its own stress field which dominates further fracture growth. At this stage, elastic rebound in the surrounding areas of fault (step 29 in Fig. 6) and a large stress drop [Fig. 7(a)] are observed; (5) The stress field of the fault concentrates around it, producing further microfracture damage around and ahead of its tips. This stress concentration leads to the
Two steps for AE source locations are presented in Fig. 8. The results related to AE behavior provide a number of important clues regarding the progressive failure leading to brittle faulting. The pre-nucleation of AE locations indicates that due to the heterogeneity in a rock mass, AE source locations remain diffuse before the main shocks [Fig. 8(a)]. During the main shock, a large number of seismic events occur along a narrow zone [Fig. 8(b)], i.e. the rock fails with a "self-organizing" mechanism--from "disorder to order". The AE behavior observed in these simulations is comparable to the experimental results described by Lockner et al. [14]. The curve of event rate with time is shown in Fig. 7(b), and the frequency-magnitude relation is also shown in Fig. 7(c). It is worth noting that every sudden increase of event rate in Fig. 7(b) results in a sudden stress drop in Fig. 8(a), as in the case of laboratory tests [9, 10]. This
TANG: SIMULATION OF PROGRESSIVE FAILURE
257
forepeak part of the relation, however, it shows a difference between the numerical model and the G u t e n b e r ~ R i c h t e r linear relation. This difference has been found by many seismologists in their observations or experiments [15]. APPLICATION IN MODELING PILLAR COLLAPSE AND ASSOCIATED SEISMICITY
Fig. 8. Seismic events and source location.
quantitative approach can be an improvement on the more traditional methods of AE analysis in correlating the AE activity with the weakening of the material. Although the relation between AE frequency and magnitude is not the focus of this paper, it is still noteworthy that the relation between the frequency of events and magnitude is approximately linear in the post-peak part of the curve shown in Fig. 7(c). This conclusion is significant because it does verify the empirical Gutenberg-Richter relation and in turn provides validity for our numerical model. For the
When designing a tunnel in a civil engineering project, elastic analysis is generally accepted as a tool to help in analyzing the stability of a permanent opening. However, in designing a mining underground opening, it is common for designers to be faced with distinctive problems related to non-linear behavior of the rock mass in which the stresses may reach magnitudes large enough to cause failure of rock surrounding the openings. Especially in burst prone ground, it is more difficult to obtain solutions to these problems than in non-burst prone ground. Since ground control problems in mining are always associated with progressive failure due to material heterogeneity, it follows that any attempt to predict ground conditions with conventional elastic analysis or even plastic analysis techniques without considering the influence of heterogeneity is useful if only to indicate potential failure zones or, at best, to predict yielding caused by geometrical reasons. In this section, progressive failure leading to a crown pillar collapse is simulated. The model set up is illustrated in Fig. 9. Although the geometry of this model comes from the INCO's Creighton Mine, Sudbury, Canada [16], it is not intended in this paper to model exactly according to the situation in the mine. The purpose of this example is for demonstrating the potential of RFPA 2Din solving mining design problems. There is a fault located along the footwall contact. The upper excavation represents the current mining extending from hangingwall to footwall. Between the upper
Fig. 9. Numerical model of a cross-section perpendicular to a strike.
258
TANG: SIMULATION OF PROGRESSIVE FAILURE
excavation and the lower excavation is the crown pillar. Assume that due to the mining process surrounding the modeled area there is a large horizontal far-field displacement rate compared with the vertical far-field displacement. Therefore, the vertical confined boundary condition is selected for the simulation. The horizontal boundary is compressed in a displacement control mode. The rock properties are chosen so as to simulate a brittle rock. The scale parameter for strength, a0, is assumed to be 60 MPa and homogeneity index m = 4; for elastic constants, E0 = 65 GPa; and v = 0.25. The ratio of compression to tensile strength, s = ~/a~, is assumed to be 3. Due to the lack of material properties for the fault gouge, the strength and elastic constant of
the gouge material are arbitrarily chosen as 90% lower than the strength and elastic constant of rock, respectively. There is no difference for elements between the fault and the surrounding rock except that lower strength and elastic constants are set for the fault. The size of the mesh was 100 x 70 with 7000 elements. The simulations are limited to a two-dimensional problem in plane strain conditions. The simulated result is shown in Fig. 10. The seismic event rate for comparison is also given in Fig. 10. It can be observed that, although mathematically R F P A 2D is completely a simple static continuum elastic code, the outcome of this numerical model shows that the code has the ability to reveal the evolutionary nature of the
ii~
~W
...............~;~i~Nii~ii~
ii!!ii~i~;: ............
Fig. lO(a).
7
....
~
~¸~
° ~ ~.~N~
~
~~ ' ~ ~ ' ~..................... ; ~ ~ ~iii!~i iil
TANG: SIMULATION OF PROGRESSIVE FAILURE ........
....
......
............. ....
Step:
-~
i
II
259
5:
: : : :::'~
,
;_;
'
Fig. 10. Numerical simulation of progressive faiiure leading to collapse of crown pillar and associated seismicity. fracture phenomenon from microfracture scale to global failure. It can be seen that the whole failure process can be divided into four main stages: (1) The fault slips over most of its length with the bottom end connected to the upper left corner of the bottom excavation, inducing many seismic events along the fault. Some seismic events are also detected in the central area of the crown pillar. This corresponds to steps 1 - 3 i n Fig. 10; (2) A large triangular block forms in the roof of the upper excavation (step 4), and the block attempts to move into the opening (steps 5-50). In the bottom part of the excavation, however, more fractures develop in the floor leading to dilation of the ground surface, and
then upheaval phenomena are observed in the following steps (steps 17 50); (3) The fault slipping becomes drastic resulting in the sudden bending fracture of the crown pillar in the end near the hangingwaU (step 17), and then the crown pillar becomes one long block rotated in the manner shown in the following steps. The modeling of this example again demonstrates the advantage of RFPA 2D which can model the progressive failure process step by step, which is important in the numerical simulation of failure mechanisms in rock materials since progressive failure will result in a stress redistribution which may be very different from the original stress distribution.
260
TANG: SIMULATION OF PROGRESSIVE FAILURE
Fig. 11. Numerical simulation of crown pillar without supports, with weakened supports and with strong supports.
Another advantage of RFPA 2D is that it can be used to help examine the feasibility of a support design of underground openings in burst prone ground. For example, from the above simulation, three areas are found to be the locations with intense micro-seismic activity (as shown in Fig. 10). In order to find the way to prevent these failures, support methods are modeled by installing supports in those event-prone locations to keep the crown pillar stable. The simulation result for three supports installed in the upper and lower parts of
the excavation is shown in Fig. 11 (b). It is found that the supports are not strong enough and failure occurs again. Even though, the fault slip is prevented. Figure 1 l(c) shows a result for stronger supports in which the roof and floor in both excavations are confined completely. Although failure still occurs in the roof of the upper excavation and in the floor of the lower excavation, the crown pillar becomes less failed .and is stable during the whole loading process, as shown in Fig. 1 l(c). Comparing with the weaker supports shown in
TANG: SIMULATION OF PROGRESSIVE FAILURE Fig. 1 l(b), it is o b s e r v e d t h a t the stability function o f the s u p p o r t system p l a y s an i m p o r t a n t role as s h o w n in Fig. 1 l(c). CONCLUSION M o s t o f the n u m e r i c a l s i m u l a t i o n o f r o c k failure has involved stress a n d strain d e t e r m i n a t i o n s . Very little w o r k has been d o n e in n u m e r i c a l m o d e l i n g o f seismic activities i n d u c e d in progressive failure leading to collapse. Since the stability c o n t r o l p r o b l e m s in m i n i n g are often associated with r o c k mass d a m a g e a n d seismic activities, it follows t h a t a n y a t t e m p t to predict u n s t a b l e failure b a s e d on stress o r strain d e t e r m i n a t i o n s m u s t be limited to where the m e a s u r e m e n t s are taken. It is seen in this p a p e r t h a t the a m o u n t o f d a m a g e occurring to the r o c k m a s s d u r i n g failure, a n d hence the possibility o f r u p t u r e , can be p r e d i c t e d better b y m e t h o d s which a c c o m m o d a t e the progressive failure a n d a s s o c i a t e d seismic activities. T h e c o m b i n e d use o f failure event m o d e l i n g a n d event source l o c a t i o n d e t e r m i n a t i o n s (or as an a p p r o x i m a t i o n , AE), in c o n j u n c t i o n with s t r e s s - s t r a i n simulation, will p r o v i d e a d d i t i o n a l insight into the p r o b l e m . T h e successful a p p l i c a t i o n s in this p a p e r have s h o w n t h a t the R F P A 2D is in principle c a p a b l e o f m o d e l i n g geological processes a n d r o c k engineering p r o b l e m s . The r e p r o d u c i b i l i t y o f p h e n o m e n a , such as fault initiation a n d d e v e l o p m e n t , elastic r e b o u n d , uplift o r seismic activities, etc. offers an attractive s u p p l e m e n t to physical m o d e l tests in b o t h geological a n d g e o m e c h a n i c a l p r o b l e m s . This refers especially to m i n i n g in b u r s t - p r o n e g r o u n d , tunnel projects a n d the u n d e r g r o u n d storage o f nuclear waste at great depth. C u r r e n t d e v e l o p m e n t efforts focus on a t h r e e - d i m e n sional model. Acknowledgements--The work was funded by China National Natural
Science Foundation and China National Education Commission. This work was initiated in 1985 under the supervision of Professor X. H. Xu. The academic visit to Professor J. A. Hudson in 1991-1992 resulted in the formation of the theoretical frame work. The discussion with Professor C. Fairhurst in Kingston, Canada, 1993, gave me important encouragement to continue my work in rock failure
261
instability, and the recent visit to Professor P. K. Kaiser made me feel specially interested in the problem of mining induced seismicity. The author would like to take this opportunity to give my thanks to all of them.
Accepted for publication 21 June 1996.
REFERENCES 1. Sterpi D., Cividini A. and Donelli M. Numerical analysis of a shallow excavation in strain softening rock. Proceedings o f 8th Int. Congress on Rock Mechanics, pp. 545-549. Tokyo (1995). 2. Potyondy D., Cundall P. and Lee C. Modeling of notch formation in the URL mine-by tunnel: Phase I--Micromechanical experiments, Report to URL by Itasca (1995). 3. Shi G. H. Block system modeling by discontinuous deformation analysis. Southampton U.K. and Boston U.S.A. (1993). 4. Kim Kunsoo and Yao C. Y. Effects of micromechanical property variation on fracture processes in simple tension. In Rock Mechanics (Edited by Daemen and Schultz), pp. 471-476. Balkema, Rotterdam (1995). 5. Tang C. A., Hudson J. A. and Xu X. H. Rock failure instability and related aspects of earthquake mechanisms. China Coal Industry Publishing House, Beijing (1993). 6. Kemeny J. M. Meeting reports. Geotech. Geological Engng 9, 83-85 (1991). 7. Blair S. C. and Cook N. G. W. Statistic model for rock fracture. In Rock Mechanics (Edited by Tillerson and Wawersik), pp. 729-735. Balkema, Rotterdam (1992). 8. Tang C. A. Numerical simulation of rock failure process. Proceedings of the 2nd Youth Symposium on Rock Mechanics in China, Chengdu (1995).
9. Cox S. J. D. and Meredith P. G. Microcrack formation and material softening in rock measured by monitoring AE. Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. 30, 11-24 (1993). 10. Khair A. W. et al. Acoustic emission pattern. Third Conference on Acoustic Emission/Microseismic Activity in Geological Structures and Materials. The Pennsylvania State University, U.S.A. (1981).
11. Plischke B., Dallmann R. and Chedmail J. F. Numerical simulation of geological processes. Proceedings of the 7th Int. Congress on Rock Mechanics, pp. 1909-1912. Paris (1991). 12. Rikitake T. Earthquake Forecasting and Warning. Centre for Academic Publications Japan/Tokyo (1982). 13. Mogi K. Earthquake Prediction. Academic Press, Tokyo (1985). 14. Lockner D. A., Moore D. E. and Reches Z. E. Microcrack interaction leading to shear fracture. In Rock Mechanics (Edited by Tillerson and Wawersike), pp. 807-815. Balkema, Rotterdam (1992). 15. Latousakis J. and Drakopoulos J. A modified formula for frequency-magnitude distribution. PAGEOPH. 125, 753-764 (1987). 16. Wiles T. D. Non-linear modeling of mining at great depth: a case history at Creighton mine. Int. Symposium on Rock Mechanics and Rock Physics at Great Depth. Pau, France (1989).