International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Journal Review Article
A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering$
L. Jing*
Division of Engineering Geology, Royal Institute of Technology, Technikenringen 72, Stockholm S-100 44, Sweden Accepted 20 January 2003
Abstract The purpose of this review paper is to present the techniques, advances, problems and likely future developments in numerical modelling for rock mechanics. Such modelling is essential for studying the fundamental processes occurring in rocks and for rock engineering design. The review begins by explaining the special nature of rock masses and the consequential dif?culties when attempting to model their inherent characteristics of discontinuousness, anisotropy, inhomogeneity and inelasticity. The rock engineering design backdrop to the review is also presented. The different types of numerical models are outlined in Section 2, together with a discussion on how to obtain the necessary parameters for the models. There is also discussion on the value that is obtained from the modelling, especially the enhanced understanding of those mechanisms initiated by engineering perturbations. In Section 3, the largest section, states-of-the-art and advances associated with the main methods are presented in detail. In many cases, for the model to adequately represent the rock reality, it is necessary to incorporate couplings between the thermal, hydraulic and mechanical processes. The physical processes and the equations characterizing the coupled behaviour are included in Section 4, with an illustrative example and discussion on the likely future development of coupled models. Finally, in Section 5, the advances and outstanding issues in the subject are listed and in Section 6 there are speci?c recommendations concerning quality control, enhancing con?dence in the models, and the potential future developments. r 2003 Elsevier Science Ltd. All rights reserved.
Keywords: Review; Rock mechanics; Numerical modelling; Design; Coupled processes; Outstanding issues
Contents 1. Introduction . . . . . . . . . . . . 1.1. Special nature of rock masses 1.2. Rock mechanics modelling for 1.3. Scope of this review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . rock engineering design and construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284 285 286 287 287 287 291 291 293 293 293
2.
Numerical methods in rock mechanics . . . . . . . . . . . . . . . . . 2.1. Numerical methods for modelling continuous and discontinuous 2.2. Characterization of rock masses for numerical methods . . . . . 2.3. Enhanced understanding provided by numerical methods . . . .
. . . . . . . rock masses . . . . . . . . . . . . . .
3.
Numerical techniques for rock mechanics: states-of-the-art . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Finite Difference Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1. Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
$ This is the second of a series of Journal Review Articles commissioned by the Editor. The series consists of articles reviewing signi?cant or topical subjects, or subjects requiring expert explanation. This Review is a signi?cantly expanded version of the ‘‘Numerical Methods in Rock Mechanics’’ CivilZone Review paper which was published in Vol. 39, No. 4, June 2002, pp. 409–427, and is longer than usual papers in order to do justice to the subject. Also, an enhanced referencing system has been used here with the references being provided in two ways: ?rstly, by author and date, so that this information is contained directly in the text; and, secondly, by bracketted numbers, following the standard Journal format. *Tel.: +46-8-790-6808; fax: +46-8-790-6810. E-mail addresses: lanru@kth.se (L. Jing).
1365-1609/03/$ - see front matter r 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S1365-1609(03)00013-3
284
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
3.2.
3.3.
3.4.
3.5.
3.6.
3.7. 3.8.
3.1.2. Finite volume approach of FDM and its application to stress analysis . . . . 3.1.3. Fracture and non-linear analyses with FDM/FVM . . . . . . . . . . . . . . Finite Element Method and related methods . . . . . . . . . . . . . . . . . . . . . . 3.2.1. Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2. Fracture analysis with the FEM . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3. Meshless (meshfree) methods . . . . . . . . . . . . . . . . . . . . . . . . . Boundary Element Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1. Basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2. Fracture analysis with BEM . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3. Alternative formulations associated with BEM . . . . . . . . . . . . . . . . Basic features of the Discrete Element Method (DEM) . . . . . . . . . . . . . . . . 3.4.1. Explicit DEM—Distinct Element Method: block systems . . . . . . . . . . . 3.4.2. Implicit DEM—Discontinuous Deformation Analysis method: block systems 3.4.3. Key Block theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.4. DEM formulations for particle systems . . . . . . . . . . . . . . . . . . . . 3.4.5. Dynamic lattice network models . . . . . . . . . . . . . . . . . . . . . . . . Discrete Fracture Network method . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1. DFN-the basic concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2. Stochastic simulations of fracture networks . . . . . . . . . . . . . . . . . . 3.5.3. Solution of the ?ow ?elds within fractures . . . . . . . . . . . . . . . . . . 3.5.4. Issues of importance and dif?culty . . . . . . . . . . . . . . . . . . . . . . 3.5.5. Alternative formulations—percolation theory . . . . . . . . . . . . . . . . . Hybrid models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1. Hybrid FEM/BEM models . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2. Hybrid DEM/BEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.3. Alternative hybrid models . . . . . . . . . . . . . . . . . . . . . . . . . . . Neural networks/empirical techniques . . . . . . . . . . . . . . . . . . . . . . . . . Constitutive models of rocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.1. Classical constitutive models of rocks . . . . . . . . . . . . . . . . . . . . . 3.8.2. Failure criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.3. Time effects and viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.4. Size effects and homogenization . . . . . . . . . . . . . . . . . . . . . . . . 3.8.5. Damage mechanics models . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.6. Rock fracture models . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
294 294 295 295 296 298 301 301 303 304 307 308 312 314 314 314 315 315 315 316 316 317 318 318 318 318 319 320 320 320 321 321 322 323 325 327 327 328 329 329 331 332 333 333 333 333 333 333 333 335 335
4. 5.
Coupled thermo-hydro-mechanical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inverse solution methods and applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. Displacement-based back analysis for rock engineering . . . . . . . . . . . . . . . . . . . . . . . 5.2. Pressure-based inverse solution for groundwater ?ow and reservoir analysis . . . . . . . . . . . . Advances and outstanding issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. Advances in numerical modelling in rock mechanics . . . . . . . . . . . . . . . . . . . . . . . . . 6.2. Issues of special importance and dif?culty in numerical modelling for rock mechanics . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1. Main numerical modelling methods . . . . . . . . . . . . 7.1.1. Finite Element Method . . . . . . . . . . . . . . 7.1.2. Boundary Element Method . . . . . . . . . . . . 7.1.3. Finite Difference Method/Finite Volume Method . 7.1.4. Discrete element method (DEM) . . . . . . . . . 7.1.5. Discrete Fracture Network model . . . . . . . . . 7.2. General comments regarding numerical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.
7.
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1. Introduction The purpose of this Journal Review Article is to present the techniques, advances, problems and likely future developments in numerical modelling for rock
mechanics. In this Section, the review is prefaced by noting the special nature and idiosyncracies of rock masses—and hence some of the dif?culties associated with capturing the rock reality in the numerical models. The utility of numerical modelling in providing
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
285
understanding for rock engineering design and construction is also explained. Finally, we note here the scope and content of the Review with its emphasis on summarizing trends and providing an extensive literature source. 1.1. Special nature of rock masses The reason for the general dif?culty in modelling rock masses, by whatever numerical method, is that rock is a natural geological material, and so the physical or engineering properties have to be established, rather than de?ned through a manufacturing process. The rock mass is largely Discontinuous, Anisotropic, Inhomogeneous and Not-Elastic (DIANE), (Harrison and Hudson, 2000) [1]. Rock masses are under stress and continuously loaded by dynamic movements of the upper crust of the Earth, such as tectonic movements, earthquakes, land uplifting/subsidence, glaciation cycles and tides. A rock mass is also a fractured porous medium containing ?uids in either liquid or gas phases, e.g. water, oil, natural gases and air, under complex in situ conditions of stresses, temperature and ?uid pressures. The complex combination of constituents and its long history of formation make rock masses a dif?cult material for mathematical representation via numerical modelling. In relation to the generally discontinuous nature of rock masses, the photograph of a blasted rock surface in Fig. 1 highlights the fact that rock masses contain through-going pre-existing fractures,1 as well as fractures introduced by the excavation process. Most of the fractures visible in Fig. 1 are pre-existing natural fractures. Although these rock fractures have occurred naturally through geological processes, their formation is governed by mechanical principles, as illustrated by the three main sets of fractures that, in this case, are mutually orthogonal and divide the rock mass into cuboids. The fractures are most often clustered in certain directions resulting from their geological modes and history of formation. One of the main tasks of numerical modelling in rock mechanics is to be able to characterize such mechanical discontinuities in a computer model—either explicitly or implicitly—the socalled ‘material conceptualization’. Additionally, the interaction between the rock mass and the engineering structure has to be incorporated in the modelling procedure for design, so that consequences of the construction process have also to be characterized. To adequately represent the rock mass in computational models, capturing such fracturing and the complete DIANE nature of the rock mass, plus the
The word ‘fracture’ is used in this Article to indicate natural breaks in the rock continuum, e.g. faults, joints, bedding planes, ?ssures. Thus, the term ‘fracture’ is used here as a synonym for ‘discontinuity’.
1
Fig. 1. Surface of a blasted rock mass, illustrating that pre-existing fractures can divide the rock mass into discrete blocks, and that the interaction between the rock mass and the engineering processes also needs to be modelled for the engineer to have a predictive capability for design purposes. Note the ‘half-barrels’ of the blasting boreholes.
consequences of engineering, it is necessary to be able to include the following features during model conceptualization:
*
*
*
*
*
the relevant physical processes and their mathematical representations by partial differential equations (PDEs), especially when coupled thermal, hydraulic and mechanical processes need to be considered simultaneously; the relevant mechanisms and constitutive laws with the associated variables and parameters; the pre-existing state of rock stress (the rock mass being already under stress); the pre-existing state of temperature and water pressure (the rock mass is porous, fractured, and heated by a natural geothermal heat gradient or manmade heat sources) the presence of natural fractures (the rock mass is discontinuous);
286
*
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
*
*
*
*
variations in properties at different locations (the rock mass is inhomogeneous); variations of properties in different directions (the rock mass is anisotropic); time/rate-dependent behaviour (the rock mass is not elastic and may undergo creep or plastic deformation); variations of properties at different scales (the rock mass is scale-dependent); the effects resulting from the engineering perturbations (the geometry is altered).
rest on a scienti?c foundation but require empirical judgements supported by accumulated experiences through long-term practices. This is the case because the quantity and quality of the supporting data for rock engineering design and analysis can never be complete, even though they can be perfectly de?ned in models. 1.2. Rock mechanics modelling for rock engineering design and construction Some form of predictive capability is necessary in order to coherently design an engineered structure, whether it be on the rock mass surface or within the underground rock mass, and whether it be for civil engineering addressed in this CivilZone review or for mining, petroleum or environmental engineering. The predictive capability is achieved through a variety of modelling methods. Even if one simply adopts the same design as a previously constructed structure, the rock mass condition is generally site-speci?c and one should use a computer model adopted for the speci?c site conditions to ensure that the rock mass is likely to behave in similar fashion. As rock mechanics modelling has developed for the design of rock engineering structures with widely different purposes, and because different modelling methods have been developed, we now have a wide spectrum of modelling approaches. These can be presented in different ways: the categorization into eight approaches based on four methods and two levels, as illustrated in Fig. 2, is from (Hudson, 2001) [2].
The extent to which these features can actually be incorporated into a computer model will depend on the physical processes involved and the modelling technique used; hence, both the modelling and any subsequent rock engineering design will contain subjective judgements. Rock engineering projects are becoming larger and more demanding in terms of the modelling requirements, one of which, for example, may be to include coupled thermo-hydro-mechanical (THM) behaviour into the model. A truly fully coupled model (including extra processes, such as chemistry) requires complete knowledge of the geometrical and physical properties and parameters of the fractured rock masses. Thus, the challenge is to know how to develop an adequate model. The model does not have to be complete and perfect: it only has to be adequate for the purpose. For these reasons, rock mechanics modelling and rock engineering design are both a science and an art. They
Objective
Method A
Use of pre-existing standard methods Site Investigation Precedent type analyses and modifications
Method B
Method C
Basic numerical methods, FEM, BEM, DEM, hybrid
Method D
Extended numerical methods, fully-coupled models
Analytical methods, stress-based
Lev el 1 1:1 mappi ng
Rock mass classification, RMR, Q, GSI
Database expert systems, & other systems approaches
Integrated systems approaches, internet-based
Level 2 Not 1:1 mapping
Design based on forward analysis
Design based on back analysis
Construction
Fig. 2. Four basic methods, two levels and hence eight different approaches to rock mechanics modelling and rock engineering design, from Hudson [2].
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
287
The modelling and design work starts with the objective, the top box in Fig. 2. Then there are the eight modelling and design methods in the main central box. The four columns represent the four main modelling methods:
*
* *
*
Method A: Design based on previous design experiences, Method B: Design based on simpli?ed models, Method C: Design based on modelling which attempts to capture most relevant mechanisms, and Method D: Design based on ‘all-encompassing’ modelling.
There are two rows in the large central box in Fig. 2. The top row, Level 1, includes methods in which there is an attempt to achieve one-to-one mechanism mapping in the model. In other words, a mechanism which is thought to be occurring in the rock reality and which is to be included in the model is modelled directly, such as explicit stress–strain relations. Conversely, the lower row, Level 2, includes methods in which such mechanism mapping is not direct. The consequences of, for example, the constitutive models and associated parameters may well be contained within the four modelling and design methods in Level 2, but one cannot explicitly identify the relation within the methodologies, e.g. in the rock mass classi?cation techniques. Some supporting rock mass characterization parameters will be obtained from site investigation, the lefthand box. Then the rock engineering design and construction proceeds, with a feedback loop to the modelling from construction. An important point is that in rock mechanics and engineering design, having insuf?cient data is a way of life, rather than a simple local dif?culty, and that is why the empirical approaches (i.e. classi?cation systems) have been developed and are still required. Therefore, we will also be discussing the subject of parameter representability associated with sample size, representative elemental volume (REV), homogenization/ upscaling, because these are fundamental problems associated with modelling, and are relevant to the ABCD method categories in Fig. 2. 1.3. Scope of this review The use of computers makes signi?cant contributions to all the eight modelling and design methods in Fig. 2; however, the speci?c numerical methods and approaches that are being reviewed here are used directly in Methods 1C and 1D. Also, there is concentration on the actual numerical methods (rather than computing per se or design per se) and discussion on the rock mass characterization issues related to the numerical methods. Highlighted are the techniques, advances,
coupled mechanisms, technical auditing and the ability to present the content of the modelling, the outstanding issues, and the future of this type of modelling. In short, highlighted is the special contribution that numerical models are currently making to rock mechanics. Because the focus of this Review is on the modelling concepts, the associated special features of modelling rock fractures, the main development milestones, typical application requirements, development trends, and outstanding issues of importance and dif?culty, special attention is paid to Section 3 for alternative formulations in each of the modelling methods, noting the potentials for rock mechanics problems. It is hoped that this treatment will provide readers with a comprehensive presentation of the state-of-the-art of numerical analysis in rock mechanics in general, and civil engineering applications in particular—in terms of historical background, presents status and likely future trends.
2. Numerical methods in rock mechanics Before considering the details and advances in the speci?c numerical modelling methods (presented in Section 3), an introduction is provided here to the methods and there is discussion on the continuum vs. discrete approaches. Also considered is the characterization of rock masses which is necessary to provide input to the numerical models, and there is illustration of how enhanced understanding is obtained through the use of such models. 2.1. Numerical methods for modelling continuous and discontinuous rock masses In numerical modelling of engineering problems, some problems can be represented by an adequate model using a ?nite number of well-de?ned components. The behaviour of such components is either well known, or can be independently treated mathematically. The global behaviour of the system can be determined through well-de?ned inter-relations between the individual components (elements). One typical example of such discrete systems is a beam structure. Such problems are termed discrete and the discrete representation and solution of such systems by numerical methods are usually straightforward. In other problems, the de?nition of such independent components may require an in?nite sub-division of the problem domain, and the problem can only be treated using the mathematical assumption of an in?nitesimal element, implying in theory an in?nite number of components. This usually leads to differential equations to describe the system behaviour at the ?eld points. Such
288
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
systems are termed continuous and have in?nite degrees of freedom. To solve such a continuous problem by numerical methods using digital computers, the problem domain is usually subdivided into a ?nite number of sub-domains (elements) whose behaviour is approximated by simpler mathematical descriptions with ?nite degrees of freedom. These sub-domains must satisfy both the governing differential equations of the problem and the continuity condition at their interfaces with adjacent elements. This is the socalled discretization of a continuum. It is an approximation of a continuous system with in?nite degrees of freedom by a discrete system with ?nite degrees of freedom. The continuity referred to here is a macroscopic concept. The continuum assumption implies that at all points in a problem domain, the materials cannot be torn open or broken into pieces. All material points originally in the neighbourhood of a certain point in the problem domain remain in the same neighbourhood throughout the deformation or transport process. Of course, at the microscopic scale, all materials are discrete systems. However, representing the microscopic components individually is intractable mathematically and unnecessary in practice. The individual components (elements) of a discrete system are usually treated as continuous. Their properties may either be obtained from laboratory tests if the components are indeed continuous and macroscopically homogeneous, such as elastic beam structures, or be mathematically derived from homogenization processes if the components themselves are heterogeneous or/and fractured, such as the fractured rock masses we are considering here. The concepts of continuum and discontinuum are therefore not absolute but relative and problem-speci?c, depending especially on the problem scales. This is particularly true for rock mechanics problems. For example, a block of rock isolated by large fractures zones may be treated as one of many block components in a computer model, but the block itself may contain a large number of smaller fractures that cannot be explicitly represented if the problem is to be tractable. Homogenization is then needed to derive the equivalent continuum properties of such blocks, which are then functions of the geometry of the contained fracture systems and physical properties of the intact rock matrix and the fractures. The fractured rock mass comprising the Earth’s upper crust is a discrete system. Closed-form solutions do not exist for such geometries and numerical methods must be used for solving practical problems. Due to the differences in the underlying material assumptions, different numerical methods have been developed for continuous and discrete systems.
The most commonly applied numerical methods for rock mechanics problems are: Continuum methods
* * *
the Finite Difference Method (FDM), the Finite Element Method (FEM), and the Boundary Element Method (BEM). Discontinuum methods
* *
Discrete Element Method (DEM), Discrete Fracture Network (DFN) methods. Hybrid continuum/discontinuum models
* * * *
Hybrid FEM/BEM, Hybrid DEM/DEM, Hybrid FEM/DEM, and Other hybrid models.
The FDM is a direct approximation of the governing PDEs by replacing partial derivatives with differences at regular or irregular grids imposed over problem domains, thus transferring the original PDEs into a system of algebraic equations in terms of unknowns at grid points. The solution of the system equation is obtained after imposing the necessary initial and boundary conditions. This method is the oldest member in the family of numerical methods, one that is widely applied and is the basis of the explicit approach of the DEMs. The FEM requires the division of the problem domain into a collection of sub-domains (elements) of smaller sizes and standard shapes (triangle, quadrilateral, tetrahedral, etc.) with ?xed number of nodes at the vertices and/or on the sides—the discretization. Trial functions, usually polynomial, are used to approximate the behaviour of PDEs at the element level and generate the local algebraic equations representing the behaviour of the elements. The local elemental equations are then assembled, according to the topologic relations between the nodes and elements, into a global system of algebraic equations whose solution then produces the required information in the solution domain, after imposing the properly de?ned initial and boundary conditions. The FEM is perhaps the most widely applied numerical method in engineering today because its ?exibility in handling material heterogeneity, non-linearity and boundary conditions, with many well developed and veri?ed commercial codes with large capacities in terms of computing power, material complexity and userfriendliness. (It is also the basis of the implicit approach of the DEM.) Due to the interior discretization, the FDM and FEM cannot simulate in?nitely large domains (as sometimes presented in rock engineering problems, such as half-plane or half-space problems) and the ef?ciency of the FDM and FEM will decrease
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
289
with too high a number of degrees of freedom, which are in general proportional to the numbers of nodes. The BEM, on the other hand, requires discretization at the boundary of the solution domains only, thus reducing the problem dimensions by one and greatly simplifying the input requirements. The information required in the solution domain is separately calculated from the information on the boundary, which is obtained by solution of a boundary integral equation, instead of direct solution of the PDEs, as in the FDM and FEM. It enjoys greater accuracy over the FDM and FEM at the same level of discretization and is also the most ef?cient technique for fracture propagation analysis. It is also best suited for simulating in?nitely large domains due to the use of fundamental solutions of the PDEs in such domains. The DEM for modelling a discontinuum is relatively new compared with the three methods described above and focuses mostly on applications in the ?elds of fractured or particulate geological media. The essence of the DEM is to represent the fractured medium as assemblages of blocks formed by connected fractures in the problem domain, and solve the equations of motion of these blocks through continuous detection and treatment of contacts between the blocks. The blocks can be rigid or be deformable with FDM or FEM discretizations. Large displacements caused by rigid body motion of individual blocks, including block rotation, fracture opening and complete detachments is straightforward in the DEM, but impossible in the FDM, FEM or BEM.
Fig. 3 illustrates the discretization concepts of the FDM/FEM, BEM and DEM for fractured rocks. An alternative DEM for ?uid ?ow in fractured rock masses is the DFN method that simulates ?uid ?ow through connected fracture networks, with the matrix permeability either ignored or approximated by simple means. The stress and deformation of the fractures are generally ignored as well. This method is conceptually attractive for simulating ?uid ?ow in fractured rocks when the permeability of the rock matrix is low compared to that of the fractures, and has wide applications in groundwater ?ow for civil engineering, reservoir simulation in petroleum engineering and heat energy extraction in geothermal engineering. An important difference between the continuum and discrete methods is the treatment of displacement compatibility conditions. In the continuum methods, the displacement compatibility must be enforced between internal elements, which is automatic in the cases of the FDM and BEM, but for the FEM it is maintained by keeping constant element-node connectivity topology and consistent orders of the trial (shape) functions between the neighbouring elements. However, displacement compatibility is not required between blocks in the DEM, and is replaced by the contact conditions between blocks with specially developed constitutive models for point contacts or fractures. The complete decoupling of rigid body motion mode and continuous deformation mode of individual blocks is usually adopted in DEM through the co-rotation scheme. The rigid body motion does not produce strains
joints
faults
(a)
joint element
(b)
region 1 block region 4 region 3 element of displacement discontinuity Regularized discontinuity
(d)
region 2
block
(c)
Fig. 3. Representation of a fractured rock mass shown in (a), by FDM or FEM shown in (b), BEM shown in (c), and DEM shown in (d).
290
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Persistent discontinuities
continuum
(a)
continuum Sets of discontinuities
(b)
(c)
(d)
Fig. 4. Suitability of different numerical methods for an excavation in a rock mass: (a) continuum method; (b) either continuum with fracture elements or discrete method; (c) discrete method; and (d) continuum method with equivalent properties.
inside the blocks, but it does produce displacements of blocks, often of large scale. In the continuum approach, the rigid body motion mode of deformation is generally not included because it does not produce strains in the elements. Therefore, a continuous system re?ects mainly the ‘‘material deformation’’ of the system and the discrete system re?ects mainly the ‘‘member (unit, or component) movement’’ of the system. The choice of continuum or discrete methods depends on many problem-speci?c factors, but mainly on the problem scale and fracture system geometry. Fig. 4 illustrates the alternative choices for different fracture circumstances in rock mechanics problems. Continuum approaches should be used for rock masses with no fractures or with many fractures, the behaviour of the latter being established through equivalent properties established by a homogenization process (Fig. 4a and d). The continuum approach can be used if only a few fractures are present and no fracture opening and no complete block detachment is possible (Fig. 4b). The discrete approach is most suitable for moderately fractured rock masses where the number of fractures too large for continuum-withfracture-elements approach, or where large-scale displacements of individual blocks are possible (Fig. 4c). Modelling fractured rocks demands high performance numerical methods and computer codes, especially regarding fracture representations, material heterogeneity and non-linearity, coupling with ?uid ?ow and heat transfer and scale effects. It is often unnecessarily restrictive to use only one method, even less one code, to provide adequate representations for the most signi?cant features and processes: hybrid models or multiple process codes are often used in combination in practice. There are no absolute advantages of one method over another, as is explained further in the later part of this
Continuum for the far-field
excavation
Discontinuum for the near field
Boundary elements on the outer boundary
Boundary elements on the interface
Fig. 5. Hybrid model for a rock mass containing an excavation—using the DEM for the near-?eld region close to the excavation and the BEM for the far-?eld region.
review. However, some of the disadvantages inherent in one type can be avoided by combined continuumdiscrete models, termed hybrid models. In 1984, Lorig and Brady [3] presented an early computational scheme in which the far-?eld rock is modelled as a transversely isotropic continuum using the BEM and the near-?eld rock as a set of discrete element blocks de?ned by rock fractures. This type of hybrid BEM-DEM is shown in Fig. 5. The complex rock mass behaviour caused by fractures and matrix non-linearity in the near-?eld of the excavation can be ef?ciently handled by the DEM or FEM, surrounded by a BEM representation of the far?eld region with linear material behaviour without fractures. The basis for such simple representation of the far-?eld is the fact that the gradients of variation of the physical variables, such as stress, displacement or ?ow, decrease rapidly with distance from the excavation. Therefore, if the interface between the near-?eld
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
291
FEM/DEM region and far-?eld BEM region is far enough from the excavation, the BEM representation will provide an accurate enough representation to model the effects of the far-?eld on the near-?eld. The ‘model’ and the ‘computer’ now provide essential support for rock mechanics analyses and understanding. The numerical methods and computing techniques assist in formulating conceptual models and mathematical developments to integrate and unify diverse geological, mechanical, hydraulic and thermal phenomena, whose interactions would not otherwise be revealed otherwise. Moreover, such developments and progress in computer methods for rock engineering will continue because they are mainly stimulated by the prospect that they will provide the information that cannot be obtained by experiments, because conducting large-scale in situ experiments is most often not possible. In fact, full veri?cation of computer models by experiments in rock mechanics is not possible: this is because the complete geometry and properties of the fractured rock mass components will never be completely known, and so veri?cations/validations can only be partial. Working with uncertainty and variability (about processes, properties, parameters, loading conditions and histories, initial and boundary conditions, etc.) becomes a way of life in rock engineering, requiring clari?cation of source information, understanding the signi?cance of assumptions, studying propagation paths relating to the assumptions and their mathematical treatment. Clearly de?ned mathematical approaches do exist to describe, analyse, and model uncertainties and error propagation, but their application in mathematical and computer models for rock engineering is still dif?cult— simply because we do not have a reference point for making judgements, except for broad empirical judgements. Modelling errors may be found as a result of the failure of rock engineering structures or accidents, but conceptual failures and modelling mistakes may be hidden under the thick blanket of the operational success of structures. Model reliability and credibility are always relative, subjective and case-dependent. This current lack of a rigorous treatment of uncertainty in rock engineering may well be a major reason why many practising engineers, and even researchers, remain uncertain about the validity and hence applicability of mathematical models and computer methods. 2.2. Characterization of rock masses for numerical methods For the different types of numerical modelling methods described in Section 2.1, the modelling is linked to generic or speci?c rock masses by the boundary and initial conditions and the rock properties. For example, the elastic modelling of a tunnel excava-
tion at a speci?c location requires a knowledge of the in situ rock stress state and the elastic properties of the rock. If the modelling is to incorporate the main components of the rock reality—the fractures, inhomogeneity, anisotropy and inelasticity, including failure—a more extensive model and a more extensive rock mass characterization are required. The scale effect is a related and additional problem, especially where fractures are affecting the rock mass properties (da Cunha, 1990, 1993; Amadei, 2000) [4,5,6]. Some of the rock characterization problems are as follows:
*
*
*
*
*
the in situ rock stress is not easy to characterize over the region to be modelled; rock properties measured in the laboratory may not represent the values on a larger scale; rock properties cannot be measured directly on a large scale; rock properties may have to be estimated from empirical characterization techniques; the uncertainty in the rock property estimates is not easy to quantify.
These problems do not mean that we cannot supply the necessary rock characterization parameters but they do mean that the whole issue of rock characterization in relation to numerical methods must be carefully considered. Needless to say, the use of different numerical models will require different types of rock property characterization. Thus, the question of whether the numerical modelling is successful in capturing the rock reality relates to both the type of numerical model and the associated rock property characterization. Other connected issues are:
*
*
*
*
is it necessary to explicitly represent the fractures or can equivalent properties be used, i.e. discontinuum vs. continuum models? to what extent is it necessary to simulate all the operating mechanisms, i.e. to use a 1:1 mechanism mapping approach, cf. Level 1 vs. Level 2 in Fig. 2? how can the combined numerical modelling technique and rock characterization method be calibrated? how can the rock characterization method be technically audited to provide some guidance on whether it is an adequate procedure? These issues are discussed in Section 5.
2.3. Enhanced understanding provided by numerical methods The purpose of numerical modelling in rock mechanics is not only to provide speci?c values of, say, stress and strain, at speci?c points but is also to enhance
292
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Tunnel Roof 140 Initial Stress State: σ3 = 27.0 MPA σ2 = 40.5 MPA σ1 = 54.0 MPA 100 80 60 40 20 σ3 0 -10 -7.5 -5 -2.5 0 roof
tunnel face advancement wall tunnel face advancement
(MPa)
σ1
120 2-D Kirsch Sol, n: σθθ = 135 MPA
σ2
Stress Vector Orientation
2.5
5
7.5
10
(m)
D = -5 m
D=5m
Fig. 6. Variation in the 3-D stress state along the roof line of an advancing tunnel face, from Eberhardt (2001) [7].
our understanding of the processes involved, particularly the changes that result from the perturbations introduced by the engineering. For example, we can use different numerical models based on linear elasticity to estimate the rock stress changes that occurs as a result of tunnelling. We can answer the question ‘What is the maximal compressive stress in the wall of a tunnel after excavation?’ However, enhanced understanding comes from a numerical demonstration of the progressive change in stresses throughout the full stress path from the original natural rock stress state to the ?nal disturbed stress state. The illustration in Fig. 6 shows the variation in the three-dimensional (3-D) stress state ahead of an advancing tunnel face and provides much more information than just the before and after values Eberhardt (2001) [7]. Also, the effects of engineering actions, such as the time when support is introduced, can be studied more coherently. Similarly, enhanced understanding comes from studying the development of all the processes that occur in rock masses as a result of engineering actions. Some of the advantages of numerical modelling in this context are the ability to:
*
Fig. 7. Use of the Particle Flow Code (PFC) to model progressive failure and acoustic emission in the uniaxial compression test (from Hazzard and Young, 2001 [8]).
*
*
develop qualitative understanding through quantitative evaluation, establish to what extent 1:1 mapping of mechanisms and properties is necessary.
* *
study rock mechanics processes from beginning to end, conduct sensitivity analyses rapidly, conduct numerical experiments for design and construction options,
An example of studying the development of a fundamental failure mechanism with a numerical code is given in Fig. 7 where the development of fracturing during the uniaxial compression of a rock specimen is shown (Hazzard and Young, 2000, [8]). This is a good example of a numerical model being able to illustrate the
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
293
progressive failure mechanism. This approach can be applied to modelling the rock mass response to engineering actions in different ?eld circumstances, leading to enhanced understanding and hence enabling the engineer to design more coherently.
ui;j ? b1 ui?1;j ? b2 ui;j?1 ? b3 ui;j?1 ? b4 ui?1;j y y y y y
i;j ? b5 ui?1;j?1 ? b6 Fy ; y
?1?
3. Numerical techniques for rock mechanics: states-of-the-art 3.1. Finite Difference Methods 3.1.1. Basic concepts The FDM is the oldest numerical method to obtain approximate solutions to PDEs in engineering, especially in ?uid dynamics, heat transfer and solid mechanics. The basic concept of FDM is to replace the partial derivatives of the objective function (e.g. displacement) by differences de?ned over certain spatial intervals in the coordinate directions, Dx; Dy; Dz; which yields a system of algebraic simultaneous equations of the objective functions at a grid (mesh) of nodes over the domain of interest (Fig. 8a) (Wheel, 1996 [9]). Solution of the simultaneous algebraic system equations, incorporating boundary conditions de?ned at boundary nodes, will then produce the required values of the objective function at all nodes, which satisfy both the governing PDFs and speci?ed boundary conditions. The conventional FDM utilizes a regular grid of nodes, such as a rectangular grid as shown in Fig. 8a. Using a standard FDM scheme, the so-called 5-point difference scheme (Fig. 8b), the resultant FDM equation at grid node ?i; j? will be expressed as combinations of function values at its four surrounding nodes. For a Navier equation of equilibrium for elastic solids in 2-D, the FDM equation of equilibrium at point ?i; j? is given as ui;j ? a1 ui?1;j ? a2 ui;j?1 ? a3 ui;j?1 ? a4 ui?1;j x x x x x
i;j ? a5 ui?1;j?1 ? a6 Fx ; x
where coef?cients ak and bk ?k ? 1; 2; y; 6? are functions of the grid intervals Dx and Dy; and the elastic i;j i;j properties of the solids, and Fx and Fy are the body forces lumped at point ?i; j?; respectively. Assembly of similar equations at all grid points will yield a global system of algebraic equations whose solution can be obtained by direct or iterative methods. FDM schemes can also be applied in the time domain with properly chosen time steps, Dt; so that function values at time t can be inferred from values at t ? Dt: The fundamental nature of FDM is the direct discretization of the governing PDEs by replacing the partial derivatives with differences de?ned at neighbouring grid points. The grid system is only a convenient way of generating objective function values at sampling points with small enough intervals between them, so that errors thus introduced are small enough to be acceptable. No local trial (or interpolation) functions are employed to approximate the PDE in the neighbourhoods of the sampling points, as is done in FEM and BEM. It is therefore the most direct and intuitive technique for the solution of the PDEs. The conventional FDM with regular grid systems does suffer from shortcomings, most of all in its in?exibility in dealing with fractures, complex boundary conditions and material inhomogeneity. This makes the standard FDM generally unsuitable for modelling practical rock mechanics problems. However, signi?cant progress has been made in the FDM so that irregular meshes, such as quadrilateral grids (Perrone and Kao, 1975, [10]) and the Voronoi grids (Brighi et al., 1998, [11]) can also be used. Although such irregular meshes can enhance the applicability of the FDM for rock mechanics problems, however, the most signi?cant improvement comes from the so-called Control Volume or Finite Volume approaches.
?y
Cell 4 l
Cell 3 k
Cell 2
Cell 1
i,j+1 i-1,j i,j i+1,j i,j-1
P
Cell 5 i Cell 6 Cell 7 j Cell 8
(a)
?x
(b)
Fig. 8. (a) Regular quadrilateral grid for the FDM and (b) irregular quadrilateral grid for the FVM (after Wheel, 1995 [9]).
294
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
3.1.2. Finite volume approach of FDM and its application to stress analysis The Finite Volume Method (FVM) is also a direct approximation of the PDEs, but in an integral sense. An elastostatics problem with body O; is divided into a ?nite number, N; of internal contiguous cells of arbitrary polyhedral (or polygonal in 2-D cases) shape, called Control Volumes (CV) Ok with boundary Gk of unit outward normal vector nk ; k ? 1; 2; y; N: The boundi ary Gk of CV Ok is comprised of a number, M k ; polygonal side (faces or line segments), Gp ; p ? k k 1; 2; y; M k ; Gk ? ,M Gp : Assuming isotropic, linear p?1 k elasticity and using Gauss’ divergence theorem, the Navier–Cauchy equation of equilibrium in terms of stress can be rewritten in terms of displacement as " k # Z N M X XZ tk dG ? fi dO i
k?1 p?1
?
" k N M X XZ
k?1 p?1 Gp k
Gp k
Ok
# sk np ij j
k dG ? Fx ? 0;
?2?
where Fik ? rgi V k is the body force vector of the CV of volume V k lumped at its centre, r is the material density and gi is the body force intensity vector, such as gravity acceleration. The task is to formulate the integrals into algebraic functions of the displacements at nodes de?ning the boundary sides Gp of Ok ; which vary with different grid k schemes. For an unstructured quadrilateral grid system (Fig. 8b), a typical cell P ?CV?; with its centre at node P; has four sides ?ij; jk; kl; li? and four nodes ?i; j; k; l?; surrounded by eight neighbouring cells with centre nodes I; J; y; O: The integral terms in Eq. (2) for the cell P are written in terms of displacement variables at the centres of cells [9], written as X X K A p up ? Ar ur ? Bp up ? Br ur ? Fx ? 0; x x y y C p up ? y X
r r
Cr ur ? Dp up ? y x
X
r
r
K Dr ur ? Fy ? 0; x
?3?
as ?exible as FEM in handling material inhomogeneity and mesh generation. As a branch of the FDM, the FVM can overcome the in?exibility of the grid generation and boundary conditions in the traditional FDM with unstructured grids of arbitrary shape. It has similarities with the FEM and is also regarded as a bridge between FDM and FEM, as pointed out in Selim (1993) and Fallah et al. (2000) [12,13]. A FVM model can be readily constructed using a standard FEM mesh, as shown in Bailey and Cross (1995) [14]. Similar examples of FVM for non-linear stress analysis with elasto-plastic and visco-plastic material models is given in Fryer et al. (1991) [15]. With proper formulations, such as static or dynamic relaxation techniques, no global system of equations in matrix form needs to be formed and solved in the FDM/ FVM approach. The formation and solution of the equations are localized, which is more ef?cient for memory and storage handling in the computer implementation. This also provides the additional advantage of more straightforward simulation of complex constitutive material behaviour, such as plasticity and damage, without iterative solutions of predictor–corrector mapping schemes that must be used in other numerical methods using global matrix equation systems, as in the FEM or BEM. The FDM/FVM approaches are therefore specially suited to simulate non-linear behaviour of solid materials. The reason is its special advantage of no-matrix-equation-solving formulation and data structure, so that integration of non-linear constitutive equations is a straightforward computer implementation step, rather than iterative prediction-mapping integration loops required in FEM. This is one of the main attractiveness of FVM such as demonstrated in Winkins (1963) and Taylor et al. (1995) [16,17]. At present, the most well-known computer codes for stress analysis for non-linear rock engineering problems using the FVM/FDM approach is perhaps the FLAC code group (ITASCA, 1993) [18], with a vertex scheme of triangle and/or quadrilateral grids. 3.1.3. Fracture and non-linear analyses with FDM/FVM Explicit representation of fractures is not easy in FDM/FVM because the ?nite difference schemes in FDM and interpolations in FVM require continuity of the functions between the neighbouring grid points. During the early development of FVM approaches, it is possible to represent weakness zones of certain thickness as collections of cells of different materials, which are not permitted to have openings or to be detached from their neighbouring cells. However, it is possible today to have special ‘‘fracture elements’’ in FVM models as in FEM, such as reported in Granet et al. (2001) and Caillabet et al. (2000) [19,20] for ?uid ?ow in deformable
where coef?cients Ap ; Ar ; Bp ; Br ; Cp ; Cr ; Dp ; Dr are functions of the cell geometry and the elastic properties of the solids, with r ? 1; 2; y; 8 running through the eight surrounding cells. This formulation of FVM with displacement variables at cell centres is called the cell-centred scheme of the FVM. If, on the other hand, the nodal displacement variables are kept as the system unknowns and the displacements at the cell centres are replaced by a combination of nodal displacements de?ning the cells, the scheme is called the vertex-centred scheme of the FVM. It is also possible to consider different material properties in different cells in the FVM, in similar ways as in the FEM. The FDM/FVM approach is therefore
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
295
porous media. On the other hand, the FDM/FVM models have been used to study the mechanisms of macroscopic fracturing processes, such as shear-band formation in the laboratory testing of rock and soil samples (Fang, 2000; Martino et al., 2002) [21,22], slope stability (Kourdey et al., 2001) [23] and glacial dynamics (Marmo and Wilson, 2001) [24]. This is achieved as a process of material failure or damage propagation at the grid points or cell centres, without creating fracture surfaces in the models. Another important improvement of FVM over the classical FDM is the use of unstructured meshes, such as triangles, arbitrary quadrilaterals, or Vorinoi grids (Mishev, 1998) [25]. This advantage, plus the ?exibility of the FVM approach in material models and boundary condition enforcement, ensures that the FDM/FVM is still one of the most popular numerical methods in rock engineering, with applications covering almost all aspects of rock mechanics, e.g. slope stability, underground openings, coupled hydromechanical or THM processes, rock mass characterization, tectonic process, and glacial dynamics. The most comprehensive coverage in this regard can be seen in Detournay and Hart (1999) [26]. The late developments in fundamentals and computer formulations can be seen in Benito et al. (2001), Onate et al. (1994), Lahrmann ? ! ! (1992), Demirdmic and Muzaferija (1994) Demirdmic et al. (2000), Jasak and Weller (2000) and Cocchi (2000) [27–33]. 3.2. Finite Element Method and related methods 3.2.1. Basic concepts Although the concept of domain discretization can be traced back to Courant (1943), and Prager and Synge (1947) [34,35], the ground-breaking work in FEM development is described in Turner et al. (1956) [36] when triangle elements were ?rst invented for structural analysis (Clough, 1960, [37]) when the term FEM was ?rst used for plane stress problems, and in Argyris (1960) [38] presenting the matrix method for structural analysis, and describing the duality of force and displacement transformations and the virtual work principle. The method was rapidly adopted and promoted in many scienti?c and engineering ?elds, as illustrated by the text books of Zienkiewicz (1977) and Bathe (1982) [39,40]. Indeed, the FEM has been the most popular numerical method in engineering sciences, including rock mechanics and rock engineering. Its popularity is largely due to its ?exibility in handling material inhomogeneity and anisotropy, complex boundary conditions and dynamic problems, together with moderate ef?ciency in dealing with complex constitutive models and fractures, i.e. the DIANE features. All these merits were very appealing to researchers and practising
engineers alike during early development in the 1960s and 1970s when the main numerical method in engineering analysis was the FDM with regular grids. Since then, the FEM method has been extended in many directions. Basically, three steps are required to complete an FEM analysis: domain discretization, local approximation, and assemblage and solution of the global matrix equation. The domain discretization involves dividing the domain into a ?nite number of internal contiguous elements of regular shapes de?ned by a ?xed number of nodes (e.g., triangle elements with three nodes in 2-D and brick elements with eight nodes in 3-D). A basic assumption in the FEM is that the unknown function, ue i over each element, can be approximated through a trial function of its nodal values of the system unknowns, uji ; in a polynomial form. The trial function must satisfy the governing PDF and is given by ue ? i
M X j?1
Nij uji ;
?4?
where the Nij are often called the shape functions (or interpolation functions) de?ned in intrinsic coordinates in order to use Gaussian quadrature integration, and M is the order of the elements. Using the shape functions, the original PDF of the problem is replaced by an algebraic system of equations written
N X i?1 e ?Kij ?fue g ? j N X ?fie ? i?1
or
Ku ? F;
?5?
e where matrix ?Kij ? is the coef?cient matrix, vector fue g is j the nodal value vector of the unknown variables, and vector ffie g is comprised of contributions from body force terms and initial/boundary conditions. e For elasticity problems, the matrix ?Kij ? is called the element stiffness matrix given by Z e ?Kij ? ? ??Bi ??Ni ??T ?Di ??Bj ? dO; ?6? Oi
where matrix ?Di ? is the elasticity matrix and matrix ?Bi ? is the geometry matrix determined by the relation between the displacement and strain. The global stiffness matrix K is banded and symmetric because the matrices ?Di ? are symmetric. Material inhomogeneity in FEM is most straightforwardly incorporated by assigning different material properties to different elements (or regions). To enforce the displacement compatibility condition, the order of shape functions along a common edge shared by two elements must be the same, so that no displacement discontinuity occurs along and across the edge. ‘‘In?nite elements’’ have also been developed in FEM to consider the effects of an in?nite far-?eld domain on the near-?eld behaviour, most notably the ‘‘in?nite domain elements’’ of Beer and Meck (1981) [41] and the
296
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
‘‘mapped in?nite elements’’ of Zienkiewicz et al. (1983) [42], with focus on geo-mechanical applications. The original concept was proposed by Bettess (1977) [43] for ?uid mechanics problems. An in?nite element formulation with body force terms was given recently by Cheng (1996) [44] with the emphasis also on geotechnical problems. The mapped in?nite elements are simply implemented using special shape functions that project boundary nodes at in?nite distances in one or two directions, where the displacements are either zero or have prescribed values. Additional nodes are needed at the imaginary in?nite locations. The in?nite domain element technique does not require additional in?nite nodes, but requires a ‘‘decay function’’ to describe the manner in which the displacements vary from mesh boundary to in?nity. The shape functions used in the in?nite element formulations are singular at the ‘‘in?nite’’ nodes. Because rock mechanics is one of the most stimulating ?elds for development of numerical methods—with many special challenges, such as fractures, property heterogeneity and anisotropy, material and geometrical non-linearity, and scale and time effects—much FEM development work and application has been speci?cally oriented towards rock mechanics problems, as illustrated in the publications of Owen and Hinton (1980), Naylor et al. (1981), Pande et al. (1990), Wittke (1990), and Beer and Watson (1992) [45–49]. The FEM has been the most widely applied numerical methods for rock mechanics problems in civil engineering because it was the ?rst numerical method with enough ?exibility for treatment of material heterogeneity, non-linear deformability (mainly plasticity), complex boundary conditions, in situ stresses and gravity. A typical recent development is given in Tang et al. (1998) [50] for simulating fracturing processes in inhomogeneous rocks with FEM. Also, the method appeared in the late 1960s and early 1970s, when the traditional FDM with regular grids could not satisfy these essential requirements for rock mechanics problems. It out-performed the conventional FDM because of these advantages.
l
thickness = 0
3.2.2. Fracture analysis with the FEM Representation of rock fractures in the FEM has been motivated by rock mechanics needs since the late 1960s, with the most notably contributions from Goodman et al. (1968), Goodman (1976), Zienkiewicz et al. (1970), Ghaboussi et al. (1973), Katona (1983), Desai et al. (1984) [51–56]. Assuming that the contact stresses and relative displacements along and across the rock fractures of a theoretical zero thickness (Fig. 9a) follow a linear relation with constant normal and shear stiffness, Kn and Ks ; Goodman et al. (1968) [51] proposed a ‘joint element’ which can be readily incorporated into an FEM process, with its local equilibrium equation given by kG uG ? f G ; ?7?
where the matrix kG is a symmetric matrix with its entries de?ned by the normal and shear stiffness, the element’s length and its orientation to the global coordinate system, respectively. The vector uG ? ?uix ; uiy ; ujx ; ujy ; uk ; uk ; ulx ; uly ?T is the nodal displacement x y vector of the four nodes (i; j; k and l) de?ning the joint element (Fig. 9b) and vector f G : The above formulation, the well-known ‘Goodman joint element’ in rock mechanics literature, has been widely implemented in FEM codes and applied to many practical rock engineering problems. Also, it has been extended to consider peak and post-peak behaviour in the shear direction. However, its formulation is based on continuum assumptions—so that large-scale opening, sliding, and complete detachment of elements are not permitted. The displacements of a joint element are of the same order of magnitude as its neighbouring continuum elements, allowing the displacement compatibility condition to be kept along and across the joint elements. Because of the zero thickness of the joint element, numerical ill-conditioning may arise due to large aspect ratios (the ratio of length to thickness) of joint elements. Zienkiewicz et al. (1970) [53] proposed a six-node fracture element with two additional nodes in the
k j t i (c) m l j k
n
s
n
i
l/2
(a) n
l/2
thickness =t
i
j
s m
p l n i
ζ
η o k j ξ
l
(b)
(d)
Fig. 9. Fracture elements in FEM by (a) Goodman et al. (1968) [51], (b) Ghaboussi et al. (1973) [54], (c) Zienkiewicz et al. (1970) [53] and (d) Buczkowski and Kleiber (1997) [60].
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
297
middle section of the element, and a small thickness (Fig. 9c). The elements can, therefore, be curved. The formulation may be seen as a ‘degenerate’ ordinary solid element of narrow thickness, and is subject to numerical ill-conditioning when the aspect ratio is too large. Using the relative displacements between the two opposite surfaces of fractures as the independent system unknowns, Ghaboussi et al. (1973) [54] proposed an FEM joint element based on the theory of plasticity (Fig. 9b). The use of the relative displacement components across and along the fractures of ?nite thickness reduces the number of unknowns of the fracture elements by half, de?ned at two nodes instead of four nodes as in Goodman’s joint elements. A ?nite thickness t is also used. The normal and shear strain components of the element are de?ned as the corresponding ratios of relative normal and shear displacements over the fracture thickness. An elasto-plastic relation between the normal and shear stresses and the normal and shear strains of the fracture element is formulated and can be implemented in the usual manner for continuum FEM analysis. This formulation is more robust in terms of numerical ill-conditioning as compared with those proposed in Goodman et al. (1968) and Zienkienwicz et al. (1970) [51,53], due to the use of the relative displacements. The ‘thin-layer’ elements developed by Desai et al. (1984) [56] are also based on a continuum assumption; these are a solid element with a specially developed constitutive model for contact and frictional sliding. The fracture element formulation in FEM has also been developed with interface element models in contact mechanics, using the FEM approach, instead of the continuum solid element approximation as mentioned above. Katona (1983) [55] developed an FEM interface element model de?ned by mating pairs of nodes, without using the normal and shear stiffness parameters, and with three states—sticking, slipping and opening—based on the Coulomb friction law. A similar approach was further discussed in Wang and Yuan (1997) [57]. In Gens et al. (1989, 1995) [58,59] 3-D FEM interface models simulating the behaviour of rock fractures were developed using the theory of plasticity. Based on the same principles, recent work by Buczkowski and Kleiber (1997) [60] considered orthotropic friction for contact interface elements in the FEM based on the theory of plasticity. The FEM interface models described above present signi?cant improvements over the early joint element models through a more systematic consideration of the kinetic and thermodynamic constraints, but they are still limited to the small displacement assumptions in the FEM with the consequence that large-scale movements across and along fracture elements are not possible. Despite these efforts, the treatment of fractures and
fracture growth remains the most important limiting factor in the application of the FEM for rock mechanics problems, especially when large number of fractures needs to be represented explicitly. The FEM suffers from the fact that the global stiffness matrix tends to be ill-conditioned when many fracture elements are incorporated. Block rotations, complete detachment and large-scale fracture opening cannot be treated because the general continuum assumption in FEM formulations requires that fracture elements cannot be torn apart. When simulating the process of fracture growth, the FEM is handicapped by the requirement of small element size, continuous re-meshing with fracture growth, and conformable fracture path and element edges. This overall shortcoming makes the FEM less ef?cient in dealing with fracture problems than its BEM counterparts. However, special algorithms have been developed in an attempt to overcome this disadvantage, e.g. using discontinuous shape functions (Wan, 1990) [61] for implicit simulation of fracture initiation and growth through bifurcation theory. A special class of FEM, often called ‘enriched FEM’, has been especially developed for fracture analysis with minimal or no re-meshing, as reported in Belytschko and Black (1999), Belytschko et al. (2001), Daux et al. (2000), Duarte et al. (2000, 2001), Dolbow et al. (2000), Jirasek and Zimmermann (2001a, b), Mo. s et al. (1999) and Sukumar et al. (2000) [62–71]. e The basic concept is direct representation of the objective function (such as displacements) with arbitrary discontinuities and discontinuous derivatives in FEM, but without need for the FEM meshes to conform to the fractures and no need for re-meshing for fracture growth. The treatment of fractures is at the element level. The surfaces of the fractures are de?ned by assigned distance functions so that their representation requires only nodal function values, represented by an additional degree of freedom in the trial functions, a jump function along the fracture and a crack tip function at the tips. The motions of the fractures are simulated using the level sets technique (Stolarska et al., 2001) [72]. Fig. 10a illustrates the non-conformal fracture-mesh relation in this technique with an arbitrary fracture intersecting a regular mesh, where circled nodes are ‘enriched’ with additional jump functions and squared nodes are ‘enriched’ with additional crack tip functions. On the other hand, the regular elements intersected by the fractures are changed into general polygons (Fig. 10b) and quadrature of the weak form in elements requires that these polygons be subdivided into standard FEM elements (such as triangles, Fig. 10c for numerical integration (Mo. s et al., 1999; Belytschko et al., 2001) e [70,71]. In Belytschko et al. (2001) [63], the enriched
298
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
(a)
(b)
(c)
Fig. 10. (a) Representation of an arbitrary fracture in a regular FEM mesh with nodes enriched by jump functions (circled nodes) or crack tip functions (squared nodes) (after Mo. s et al., 1999) [70]; (b) Details of rectangular elements intersected by a fracture, thus forming polygons; and e (c) triangularization of polygons into triangle elements for quadrature integration (after Belytschko et al., 2001) [63].
C3
1 a)
(a)
C2 b)
(b)
C1
0
(c)
Fig. 11. (a) A regular non-conformal mesh of GFEM with ‘enriched’ nodes surrounding a fracture; (b) general coverings in the manifold method; and (c) manifold method coverings with a standard triangular FEM mesh and shape functions (after Chen et al., 1998 [78]).
method was applied to a tunnel stability analysis with fractures simulated as displacement discontinuities. The ‘enriched’ FEM with jump functions and crack tip functions has improved the FEM’s capacity in fracture analysis. Coupled with the FEM’s advantage in dealing with material heterogeneity and non-linearity, this makes the ‘enriched’ FEM suitable for non-linear fracture analysis. One such example is the so-called ‘generalized ?nite element method’ (GFEM), which was developed based on the partition of unity principle (Duarte et al., 2000; Strouboulis et al., 2000, 2001) [73– 75]. The mesh in GFEM is independent of the geometry of the domain of interest and therefore can be regular regardless of the object geometry. Fractures can be simulated by their surrounding nodes ‘enriched’ by jump functions and crack tip functions (Fig. 11a). This greatly simpli?es the meshing tasks. The GFEM is in many ways similar to the so-called ‘manifold’ method except for the treatment of fractures and discrete blocks (Shi, 1991, 1992; Chen et al., 1998) [76–78]. The manifold method uses the truncated discontinuous shape functions to simulate the fractures and treat the continuum bodies, fractured bodies and assemblage of discrete blocks in a uni?ed form, and is a natural bridge between the continuum and discrete representations.
The method is formulated using a node-star covering system for constructing the trial functions (Fig. 11b). A node is associated with a covering—star, which can be a standard FEM mesh (Fig. 11c) or generated using leastsquare kernel techniques with general shapes. The integration, however, is performed analytically using Simplex integration techniques. Like GFEM, the manifold method can also have meshes independent of the domain geometry, and therefore the meshing task is greatly simpli?ed and simulation of the fracturing process does not need remeshing. The technique has been extended for applications to rock mechanics problems with large deformations and crack propagation (Wang et al., 19971a, 1997b) [79,80]. Most of the publications are included in the series of proceedings of the ICADD2 symposia (Li et al., 1995; Salami and Banks, 1996; Ohnishi, 1997; Amadei, 1999) [81–84]. 3.2.3. Meshless (meshfree) methods Besides the fracture analysis problem, the traditional FEM suffers from other shortcomings, especially the ‘locking’ effects and high demand for mesh generation.
2 ICADD (acronym for International Conference on Analysis of Discontinuous Deformation).
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
299
There are two types of ‘locking’ effects in FEM: numerical locking and element locking. ‘Numerical locking’ is the phenomenon by which numerical approximation deteriorates near some limiting values of material properties (Arnold, 1981; Babu$ka s and Suri, 1992; Suri, 1996) [85–87], or special geometry limits. Typical examples are the Poisson’s ratio for elasticity problems when the value of Poisson’s ratio is near 0.5, and shear locking for shells and plates when the thickness of the shells and plates is reduced to a small amount (Bucalen and Bathe, 1995) [88]. Proper h-, p- or hp-convergence measures need to be taken to avoid such locking effects and ensure solution convergence (Babu$ka and Suri, 1990; Chilton and Suri, 1997) s [89,90]. The ‘element locking’ in FEM is the numerical instability and can be caused by local mesh distortions, such as large aspect ratios of elements under highly concentrated loads, especially in dynamic large deformation analysis. The mesh generation is a demanding task in applying FEM for practical problems with complex interior structures and exterior boundaries. The meshes must be detailed enough to ensure proper representation of the problem geometry, solution convergence and result accuracy; yet they should also enable the computations to be completed in an economically reasonable time. These con?icting issues are problem-speci?c and must be balanced carefully between resolution and resources. The problem is critical when dealing with 3-D problems with complex geometry. Although commercial software is available for fully or semi-automatic generation of FEM meshes (as preprocessor) and for results evaluation and presentation (as post-processor), it is still up to the engineers and analysts as code users to address the issues of mesh resolution, result accuracy and computing resources. Mesh generation usually takes a much longer time than the actual calculations—because it depends almost entirely on judgement and experience, rather than theoretical guidance. A number of trial-and-error cycles are often needed to settle the issues properly. The most common method for improving the solution convergence and avoiding locking effects is by successively increasing mesh resolution, i.e. increasing the number of elements. This is called the h-convergence approach. A different approach, called p-convergence, is to increase the order of the trial (shape) functions, i.e. increasing node numbers per element while maintaining a constant element numbers. The combination of the two approaches is the hp-convergence approach in FEM, and is built into many commercial FEM codes. The aim is to avoid the ‘hourglass’ phenomenon due mainly to too low an order of trial functions, and achieve a robust solution convergence rate (Oden, 1990; Babu$ka and Suri, 1990) [91,89]. s
Signi?cant progress has been made in the last decade in the ‘meshless’ (or ‘meshfree’, ‘element-free’) method, which is closely related to FEM. In this approach, the trial functions are no longer standard but generated from neighbouring nodes within a domain of in?uence by different approximations, such as the least-square technique (Belytschko et al., 1996) [92]. The requirement for mesh generation is only generation and distribution of discrete nodes, without ?xed element-node topological relations as in the FEM. A large number of different meshless formulations have been developed over the years. Most notable among them are:
*
*
*
*
*
*
*
*
*
*
*
‘‘Smooth particle hydrodynamics’’ method (SPH) (Monaghna, 1988; Randles and Libersky, 1996) [93,94]; ‘‘Diffuse element method’’ by Nayroles et al. (1992) [95]; ‘‘Element-free Galerkin method’’ (EFG) (Belytschko et al., 1994) [96]; ‘‘Reproducing kernel particle methods’’ (RKPM) (Liu et al., 1995, 1996; Chen et al., 1996) [97–99]; ‘‘Moving least squares reproducing kernel method’’ (MLSRK) (Liu et al., 1997) [100]; ‘‘hp-cloud method’’ (Duarte and Oden, 1996; Liszka et al., 1996) [101,102]; ‘‘Partition of unity method’’ (PUM) (Melenk and Babu$ka, 1996) [103]; s ‘‘Local Petrov–Galerkin’’ (MLPG) and ‘‘local boundary integral equation’’ (LBIE) methods (Atluri and Zhu, 1998; Atluri et al., 1999) [104,105]; ‘‘Method of ?nite spheres’’ (De and Bathe, 2000) [106]; ‘‘Finite point method’’ (Onate et al., 1996; Sulsky and ? Schreyer, 1996) [107,108]; ‘‘Natural element method’’ (NEM) (Sukumar et al., 1998) based on a Voronoi tessellation of a set of nodes [109].
Naturally, the main advantage of the meshless approaches is the much reduced demand for meshing compared with standard FEM and FDM/FVM for both continuous and fractured bodies. They can still be viewed as classes of weighted residual techniques in computational mechanics like FEM and BEM, and are performed with three key operations: interpolation using trial (shape) functions; integration to derive governing algebraic equations; and solution of the ?nal system equations. There are basically three interpolation techniques: wavelets, moving least-square functions, and the partition of unity or hp-clouds. The moving leastsquare techniques are special cases of the partition of unity. The interpolation functions produced using these techniques are non-polynomial functions and this makes the integration of the weak form more demanding compared with standard FEM.
300
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Another shortcoming of many meshless approaches is the dif?culty in enforcement of essential boundary conditions. The Kronecker delta function property in the FEM/BEM shape functions ensures that the essential boundary functions are met ef?ciently. However, in many meshless methods, the generated interpolation functions do not have the Kronecher delta property at nodes, and special techniques must be applied to overcome this dif?culty, such as coupling with FEM at boundaries (Krongauz and Belytschko, 1996) [110], Lagrange multipliers (Belytschko et al., 1994) [96], penalty factors (Atluri and Zhu, 1998) [104], etc. The early meshless formulations, such as diffuse elements, EFG, hp-clouds, PUM, RKPM, and MLSPK are not true meshless techniques because a background mesh (cell) structure is still needed for integration, although not for interpolation. An example of such a cell–node structure in EFG is given in Fig. 12a. The discretization is therefore not adequately ‘free’. The ?nite point method uses a weighted least-square interpolation (therefore no element structure needed) and point collocation (thus by-passing the elements for integration), and therefore is a true meshless formulation. However, care needs to be taken in the choice of collocation points to ensure correct solutions of the derived equations. Similar formulations are also reported by Zhang et al. (2001) [111]. The MLPG and LBIE methods are true meshless techniques without using background meshes for either interpolation or integration. The trial functions are generated using moving least squares, partition of unity or Shepard functions with the shapes of circles, rectangles or ellipse at the bases. The test functions are the same as the trial functions for the MLPG and the fundamental solutions for the LBIE, respectively, over their respective support and neighbouring nodes of in?uences (Fig. 12b). The integration is performed over
sub-domains de?ning the supports for the test functions (Ote in Fig. 12b) or their intersections. A panel factor is used in MLPG to ensure the satisfaction of the essential boundary conditions. The ?nite sphere method (De and Bathe, 2000) [106] may be viewed as a special class of MLPG with circular support domain shapes. The difference between the MLPG and the method of ?nite spheres is the numerical integration scheme. The method of ?nite spheres uses Gaussian product rules on 2-D annuli and annular sectors with a larger number of integration points compared with the FEM. A similar development by Atluri and Li (2001) [112], called the ?nite cloud method, uses the combined point collocation and ?xed least-square kernel technique for constructing interpolation functions. The meshless approach greatly simpli?es the task of mesh generation in FEM since no ?xed elements are required and largely eliminates the element locking effect, with the cost of more computational effort in generating numerically the trial functions over the selected node clusters. From the pure computing performance point of view, it has not yet outperformed the FEM techniques, but it has potential for civil engineering problems in general, and rock mechanics applications in particular, due to its ?exibility in treatment of fractures, as reported by Zhang et al. (2000) [113] for analysis of jointed rock masses with block-interface models, and Belytschko et al. (2000) [114] for fracture growth in concrete. Its development was stimulated largely for simulating the mechanics of fractured rocks—because of the latter’s unusual complexity in geometry and behaviour of the hosted fractures. A contact-detection algorithm using the meshless technique was also reported by Li et al. (2001) [115] that may pave the way to extending the meshless technique to discrete block system modelling. The concept was also extended to the BEM (see Section 3.3).
Shadow Cell Node i Crack
Support of a node and its nodes of influence for trial functions
? ite
j ? tr
?? te ∩Γ = Γs
Domain of influence
Node j Quadrature point
(a)
Node
?? te =Γs
(b)
Fig. 12. (a) The cell-and-node structure of the EFG meshless method (Belytschko et al., 1994) [96]; and (b) node and support structure in MLPG and LBIE where symbols Oite and Ojtr indicate the support domains of nodes i and j; respectively, and @Ote indicates the boundaries of the support domains.
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
301
3.3. Boundary Element Methods 3.3.1. Basic concepts Unlike the FEM and FDM methods, the BEM approach initially seeks a weak solution at the global level through an integral statement, based on Betti’s reciprocal theorem and Somigliana’s identity. For a linear elasticity problem with domain O; boundary G of unit outward normal vector ni ; and constant body force fi ; for example, the integral statement is written as cij uj ? Z
G
Substitution of Eqs. (10) into (9) and for Z Z t? Nj dG; Uij ? u? Nj dG; Tij ? ij ij Bi ? Z
Gk Gk
fj
Gk
qu? ij qn
dG
?11?
Eq. (8) can be written in matrix form as ?Tij ?l; k?? fuj ?k?g ? ?Uij ?l; k?? ftj ?k?g ? fBi ?k?g;
?2N?2N? ?2N?1? ?2N?2N? ?2N?1? ?2N?1?
?12?
t? uj dG ? ij
Z
G
u? tj dG ? ij
Z
G
qu? ij qn
fj dG;
?8?
where uj and tj are the displacement and traction vectors on the boundary G; the terms u? and t? are called ij ij displacement and traction kernels. The term cij is called the free term determined by the local geometry of the boundary surfaces, cij ? 1 when the ?eld point is inside the domain O: The solution of the integral Eq. (8) requires the following steps: (1) Discretization of the boundary G with a ?nite number of boundary elements. For 2-D problems, the elements are 1-D line segments which may have one node at the centre of the element (constant element), two nodes at the two ends of the line segment (linear elements) or three nodes with two end nodes and one central node (quadratic elements). Let N denote the total number of boundary elements. The boundary integral equation then is re-arranged into a sum of local integrals over all elements cij uj ?
N XZ k?1 Gk
where i; j ? 1; 2 for 2-D and 1, 2, 3 for 3-D problems, respectively, l; k ? 1; 2; y; N; and Z Tij ?l; k? ? cij dlk ? t? Nj dG: ?13? ij
Gk
(3) Evaluation of the integrals Tij ; Uij and Bi with point collocation method by setting the source point P at all boundary nodes successively. Closed-form solutions exist only for some particular cases (see, for examples, Fratantonio and Rencis, 2000; Carini et al., 1999 [116,117]), and numerical integration using Gaussian quadrature is often used. Note that singularity occurs in the above integrals when the source and ?eld points are located on the same elements, and special integration schemes need to be used to evaluate them in a Cauchy Principal Value sense. (4) Incorporation of boundary conditions and solution. Incorporation of the boundary conditions into the matrix Eq. (12) will lead to ?nal matrix equation ?A?fxg ? fbg; ?14?
t ? uj ij
dG ?
N XZ k?1 Gk
u? t j ij
dG fj dG: ?9?
?
N XZ k?1 Gk
qu? ij qn
(2) Approximation of the solution of functions locally at boundary elements by (trial) shape functions, in a similar way to that used for FEM. The difference is that only 1-D shape functions with intrinsic coordinate ?1pxp1 is needed for 2-D BEM problems, and 2-D shape functions with two intrinsic coordinates ?1pxp1 and ?1pZp1 are needed for 3-D problems. The displacement and traction functions within each element are then expressed as the sum of their nodal values of the element nodes: ui ?
m X k?1
N k uk ; i
ti ?
m X k?1
N k tk ; i
?10?
where m is the element order (=1, 2 or 3 for 2-D problems, for example), and uk and tk are the nodal i i displacement and traction values at node k; respectively.
where the global matrix ?A? is a mixture of Tij and Uij ; the unknown vector fxg is a composite of both unknown displacements and unknown boundary tractions, and the known vector fbg is the sum of the body force vector fBi g and the products of Tij with known displacements and Uij with known tractions, respectively. The resultant Eq. (14) is usually fully populated and asymmetric, leading to fewer choices for ef?cient equation solvers, compared with the sparse and symmetric matrices encountered in the FEM. The solution of Eq. (14) will yield the values of unknown displacements and tractions at boundary nodes. Therefore all boundary values of displacements and tractions are obtained. (5) Evaluation of displacements and stresses inside the domain. For practical problems, it is often the stresses and displacements at some points inside the domain of interest that have special signi?cance. Unlike the FEM in which the desired data are automatically produced at all interior and boundary nodes, whether some of them are needed or not, in BEM the displacement and stress values at any interior point, P; must be evaluated
302
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
separately by ui ?P? ? ?
M X k?1 M X l?1
# % Tij uk ? j
M X k?1
# % Uij tk ? j
M X k?1
# Bk ;
?15?
sij ?P? ? ?
Skij ulk ? %
M X l?1
% Dkij tlk ;
?16?
# # # where kernels Tij ; Uij ; Skij ; Dkij and Bi must be reevaluated according to the new position of the source point inside the domain (closed-form formulae for them are available in many text books on the BEM), usually without singularities unless the point is very close to %j boundary, and uk and tk are known or calculated %j displacement and traction vectors at all boundary nodes. The boundary integral equation method was used for the ?rst time by Jaswon (1963) and Symm (1963) [118,119] for solving potential problems. The breakthrough for stress analysis in solids came with the work of Rizzo (1967), Cruse and Rizzo (1968) and Cruse (1973), and Cruse (1978) [120–123] for fracture mechanics applications, based on Betti’s reciprocal theorem (Betti, 1872) [124] and Somigliana’s identity in elasticity theory (Somigliana, 1885) [125]. The basic equations can also be derived using the weighted residual principle, as presented in Brebbia and Dominguez (1977) [126]. The introduction of isoparametric elements using different orders of shape functions in the same fashion as that in FEM, by Lachat and Watson (1976) and Watson (1979) [127,128], greatly enhanced the BEM’s applicability for stress analysis problems. All these works have established the status of BEM as an ef?cient numerical method for solving general engineering mechanics problems. The most notable original developments of BEM application in the ?eld of rock mechanics may be attributed to Crouch and Fairhurst (1973), Brady and Bray (1978) and Crouch and Star?eld (1983) [129–131], quickly followed by many as reported in the rock mechanics and BEM works (Hoek and Brown, 1982; Brebbia, 1987; Pande et al., 1990; Beer and Watson, 1992) [132,133,47,49] for general stress and deformation analysis for underground excavations, soil-structure interactions, groundwater ?ow and fracturing processes, and a large number of journal and conference publications. Notable examples are the work for stress/deformation analysis of underground excavations with or without faults (Venturini and Brebbia, 1983; Beer and Pousen, 1995a, b; Kayupov and Kuriyagama, 1996; Cerrolaza and Garcia, 1997; Pan et al., 1998; Shou, 2000) [133–139], dynamic problems (Tian, 1990; Siebrits and Crouch, 1993; Birgisson and Crouch, 1998) [140–142], in situ stress and elastic property interpretation (Wang and Ma, 1986; Jing 1987) [143,144], and borehole tests for permeability measurements (Lafhaj and Shahrour, 2000) [145]. Since the early 80s, an important develop-
mental thrust concerns BEM formulations for coupled thermo-mechanical and hydro-mechanical processes, such as the work reported in Pan and Maier (1997), Elzein (2000) and Ghassemi et al. (2001) [146–148]. Due to the BEM’s advantage in reducing model dimensions, 3-D applications are also reported, especially using DDM for stress and deformation analysis, such as Kuriyama and Mizuta (1993), Kuriyama et al. (1995) and Cayol and Cornet (1997) [149–151]. The main advantage of the BEM is the reduction of the computational model dimension by one, with much simpler mesh generation and therefore input data preparation, compared with full domain discretization methods such as the FEM and FDM. Using the same level of discretization, the BEM is often more accurate than the FEM and FDM, due to its direct integral formulation. In addition, solutions inside the domain are continuous, unlike the pointwise discontinuous solutions obtained by the FEM and FDM groups. The solution domains of BEM can be divided into several sub-domains with different material properties, and this will often reduce the calculation time as well. The method is also suitable for considering in?nite domains (full or half space/plane), due to its use of the fundamental solutions. However, in general, the BEM is not as ef?cient as the FEM in dealing with material heterogeneity, because it cannot have as many sub-domains as elements in the FEM. The BEM is also not as ef?cient as the FEM in simulating non-linear material behaviour, such as plasticity and damage evolution processes, because domain integrals are often presented in these problems. The BEM is more suitable for solving problems of fracturing in homogeneous and linearly elastic bodies. The BEM formulation described above is called the direct formulation in which the displacements and tractions in the equations have clear physical meanings, are the basic unknowns of the boundary integral equations which are explicitly described on the problem boundary, and can be directly obtained by the solution of the integral equations. In the indirect formulation, on the other hand, the basic unknowns have no physical meanings and are just ?ctitious source densities related to the physical variables such as displacements and tractions. The typical indirect BEMs are the Displacement Discontinuity Method (DDM) by Crouch (1976) [152] for 2-D problems and Weaver (1977) [153] for 3-D problems and the Fictitious Stress Method by Crouch and Star?eld (1983) [131]. The basic concept of the indirect approach is to place the ?nite domain of interest into an imaginary in?nitely large domain (full or halfplane or spaces) to derive the boundary integral equations relating the physical variables, such as displacements and tractions, to ?ctitious source densities, such as ?ctitious load (stress) or displacement
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
303
discontinuity. Dunbar (1985) [154] showed the equivalence between the direct and indirect BEM approaches. 3.3.2. Fracture analysis with BEM To apply standard direct BEM for fracture analysis, the fractures must be assumed to have two opposite surfaces, except at the apex of the fracture tip where special singular tip elements must be used. Denote Gc as the path of the fractures in the domain O with its two opposite surfaces represented by G? and G? ; respecc c tively, Somigliana’s identity (when the ?eld point is on the boundary) can be written as Z Z uj ? cij Duj ? t? uj dG ? t? Duj dG ij ij ? Z
G G
u? t j ij
dG ?
Z
Gc
Gc
u? Dtj ij
dG ?
Z
qu? ij fj dG; G qn
?17?
where Dui and Dti are the displacement and traction jumps across the two opposite surfaces of the fractures. Because of the very small thickness of fractures, the two nodes at the opposite surfaces of a fracture will in fact occupy the same coordinates. This will naturally lead to singular global stiffness matrices if the same boundary conditions (or unknowns) are speci?ed at the two opposite fracture surfaces. Also, any set of equal and opposite tractions on the fracture surfaces will lead to the same equation since Dti ? 0: In addition, the displacement difference Duj becomes an additional unknown on Gc besides uj : To overcome these dif?culties, two techniques were proposed. One was to divide the problem domain into multiple sub-domains with fractures along their interfaces (Fig. 13a), by Blandford et al. (1981) [155]. This way, the stiffness matrices contributed by opposite surfaces of the same fracture will belong to different sub-domain stiffness matrices; thus, the singularity of the global matrix is avoided. This technique, however, requires the knowledge of fracturing paths (used for deciding sub-regions) and growth rate (for deciding element sizes), which is determined by the solution of the problem itself, before the problem solution, and may not
be applicable for many practical problems without symmetry in geometry and boundary conditions. The second technique is the Dual Boundary Element Method (DBEM). The essence of this technique is to apply displacement boundary equations at one surface of a fracture element and traction boundary equations at its opposite surface, although the two opposing surfaces occupy practically the same space in the model. The general mixed mode fracture analysis can be performed naturally in a single domain (Fig. 13b). The term DBEM was ?rst presented in Portela (1992) and Portela et al. (1992, 1993) [156–158], and was extended to 3-D crack growth problems by Mi and Aliabadi (1992, 1994) [159,160]. However, the original concept of using two independent boundary integral equations for fracture analysis, one displacement equation and another its normal derivative, was developed ?rst by Watson (1979) [128]. Special crack tip elements, such as developed in Yamada et al. (1979) and Aliabadi and Rooke (1991) [161,162], are used at the fracture tips to account for the stress and displacement singularity. The DDM has been widely applied to simulate fracturing processes in fracture mechanics in general and in rock fracture propagation problems in particular due to the advantage that the fractures can be represented by single fracture elements without need for separate representation of their two opposite surfaces, as should be done in the direct BEM solutions. It was developed by Crouch and Star?eld (1983) [131] with open fractures, and was extended to fractures with contact and friction by Wen and Wang (1991) [163] and Shen (1991) [164] for mechanical and rock engineering analyses, respectively. The ?ctitious unknowns are the displacement discontinuity Dui acting on the boundary G of a ?nite body of domain O; inserted in an in?nitely large half (or full) space. The displacement and traction is given by Z Z ? u? Duj dG; ui ? cij Duj ? uij Duj dG ? % ij G Gc Z Z %ij ti ? cij Duj ? t? Duj dG ? t? Duj dG ?18? ij
G Gc
?2
Γc Γ Γc+
_
Γc Γ Γ
Γc
?1
?
(b) (c)
?
(a)
Fig. 13. Illustrative meshes for fracture analysis with BEM: (a) sub-domain, direct BEM; (b) single domain, dual BEM; and (c) single domain DDM.
304
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
%ij for ?eld points on G: The kernels u? and t? are the % ij fundamental solutions due to unit discontinuity per unit length (for 2-D) or area (for 3-D), see Appendix A in Wen (1996) [165] for details. By writing displacement and traction equations in Eq. (18) at boundary and fracture nodes with corresponding known displacements and tractions, respectively, a determinate matrix equation can be produced and its solution will yield all discontinuities Dui at all boundary and fracture nodes, which can then be used to produce displacements or stresses at required points inside domain O: For fracturing analysis using numerical methods, there are two main tasks, besides the establishment of the system equations. They are the evaluation of the stress intensity factors (SIF) and simulating fracture growth, based on determination of fracture deformation modes (I, II and III) and different criteria for fracture growth, such as maximum tensile strength or energy release rate. The details can be obtained in Aliabadi and Rooke (1991) [162], Wen (1996) [165] and Mi (1996) [166], for both displacement discontinuity and DBEM formulations applied for fracture growth analysis. A comprehensive coverage of the subject is given in Alliabadi (1999) [167], together with topics of rock fragmentation using DFEM and transport using DFN. Analysing fracturing processes using BEM is challenging, especially for rock mechanics problems. On the one hand, what happens exactly at the fracture tips in rocks still remains to be adequately understood, with the additional complexities caused by the microscopic heterogeneity and non-linearity at the fracture tip scale, especially regarding the fracture growth rate. On the other hand, complex numerical manipulations are still needed for re-meshing following the fracture growth process so that the tip elements are added to where new fracture tips are predicted, and updating of system equations following the re-meshing, although the task is much less cumbersome than that required for domain discretization methods such as standard FEM. Due to the above dif?culties, fracture growth analyses in rock mechanics have not been widely applied, and were mostly performed with 2-D BEM using indirect formulations such as DDM, considering often a small number of isolated, non-intersecting fractures in usually, small 2-D models concerning local failure mechanisms, such as borehole breakout and mechanical breakage (Tan et al., 1998) [168]. A 3-D DDM code POLY3D is developed at Stanford University, which is able to consider a number of non-intersecting fractures in 3-D (La Pointe et al., 1999) [169]. Usually, fracture growth is ignored completely in rock mechanics applications, and only stress and deformations of large-scale fractures are included, such as faults or fracture zones, using the multi-region formulations (Crotty and Wardle, 1985) [170] and in Beer and Pousen (1995a, b) [134,135]. This
can be justi?ed by the fact that the scale of fracture growth, under the normal loading conditions encountered in most civil engineering projects, is small and omitting this factor will not cause major misinterpretations of the behaviour of rock masses under consideration. The assumption, however, will not be true for other problems—such as borehole stability and rock spalling, where the former is dominated by in situ stresses, and the latter is controlled by the local stresses, local fracture geometry and dynamic fracturing processes. The effects of ?uid pressures and heat gradients on fracturing process in rock, either static or dynamic, is not properly understood, even less analysed by numerical methods, due to the added complexity in the physics, equations and numerical solution procedures. However, considering their potential signi?cance for the performance and safety of many environment-oriented civil engineering projects in fractured rocks, the need for such studies is clear. 3.3.3. Alternative formulations associated with BEM The standard BEM, dual BEM and DDM as presented above have a common feature: the ?nal coef?cient matrices of the equations are fully populated and asymmetric, due to the traditional nodal collocation technique. This makes the storage of the global coef?cient matrix and solution of the ?nal equation system less ef?cient, compared with FEM. They suffer also from another drawback: the need for special treatment of sharp corners on the boundary surfaces (curves) or at the fracture intersections, because of the change of directions of the outward unit normal vectors at the corners. This causes an inadequate number of equations in the ?nal system compared with that of the degrees of freedom. Arti?cial corner smoothing, additional nodes or special corner elements are usually the techniques applied to solve this particular dif?culty. In a special formulation of BEM, called the Galerkin BEM, or Galerkin Boundary Element Method (GBEM), these two shortcomings are removed automatically as the consequences of the Galerkin formulation. 3.3.3.1. Galerkin Boundary Element Method. The GBEM produces a symmetric coef?cient matrix by multiplying the traditional boundary integral by a weighted trail function and integrates it with respect to the source point on the boundary for a second time, in a Galerkin sense of weighted residual formulation. For an elasticity problem, the Somigliana identity becomes a double integral equation Z Z Z wi cij uj dG ? wi t? uj dG dG ij
G
?
Z Z
G G
G
G
wi u? tj dG dG ? ij
Z wi
G
qu? ij qn
fj dG;
?19?
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
305
where wi is the weight function. Discretizing the boundary into elements and choosing the weight functions the same as the trial (shape) functions fNu g for displacement and fNt g for tractions, respectively, in vector form and using the Galerkin approximation technique, the boundary integral equation in Eq. (19) becomes a matrix equation of the form 9 8 9 2 38 ?C uu ? ??C ut ? ?C uc ? > ftg > > ff u g > = < = < 6 7 ? ?ff t g ; 4 ??C tu ? ?C tt ? ??C tc ? 5 fug > > > > ; : ; : ?C cu ? ??C ct ? ?C cc ? ff c g fDug ?20? where the vectors ftg; fugand fDug are the unknown traction and displacement vectors at the boundary and displacement discontinuity vector along the fractures inside the domain, respectively. The right-hand vectors ff u g; ff t gand ff c g are obtained from known boundary displacement and traction conditions on the boundary and loading condition along the fractures, respectively. The coef?cient sub-matrices ?C mn ? (m; n ? u; t; and c) are contributions from double integrals over elements with speci?ed displacements (Gu ), traction (Gt ), or along the fracture paths (Gc ), respectively, written as (Carini et al., 1999) [117] ?C uu ? ? Z Z
Gu Gu
The GBEM is an attractive approach due to the symmetry of its ?nal system equation, which paves the way for the variational formulation of BEM for solving non-linear problems. It is also ?exible in choosing different formulations, and the traditional BEM can be seen as a special case. A recent review on GBEM is given by Bonnet et al. (1998) [171] in which the theoretical foundation and applications in different mechanics ?elds are summarized. However, integration of these strong and hypersingular kernels are more dif?cult than the common singular kernels in the traditional BEM. Analytical integration techniques were also proposed for elements with low order shape functions (see Salvadori, 2001) [172], and singularity subtraction techniques were also proposed to ease the task of integrations (Michael and Barbone, 1998) [173]. Applications of GBEM into rock mechanics problems has commenced (Wang et al., 2001) [174]. 3.3.3.2. Boundary Contour Method. The Boundary Contour Method (BCM) involves rearranging the standard BEM integral Eq. (8) so that the difference of the two integrals appearing on the right-hand side of Eq. (8) can be represented by a vector function Fi ? u? tj ? t? uj which is divergence free, i.e. r ? F ? 0; except ij ij at the point of singularity (Nagarajan et al., 1994, 1996). This property of Fi ensures the existence of a vector function Vi so that F ? r ? V and the BEM equation can be rewritten as Z Z I Fi dG ? Vi dl; ?22? cij uj ? ?u? tj ? t? uj ? dG ? ij ij
G G qG
fNt ?Q?gT ?G uu ?P; Q??fNt ?Q?g dG?Q? dG?P? ; ?21a?
?C ui ? ?
Z Z
Gi
fNt ?Q?gT ?G ut ?P; Q??fNt ?Q?g dG?Q? dG?P? ;
Gu
?i ? t; c?
?21b? ?C ij ? ?
Z Z
Gi
fNu ?Q?gT ?G tt ?P; Q??fNu ?Q?g dG?Q? dG?P?
Gj
?i; j ? t; c?;
?21c? ff u g ? ff t g ? ff c g ? Z Z Z
Gu
% fNt gffgT dGu ; fNu gff%t gT dGt ;
Gt
fNu gff%c gT dGc ;
Gc
?21d?
where elements in the matrices ?G mn ? (m; n ? u; t; and c) are kernel functions derived from fundamental solutions of the problems. Due to the extra integration, the singularity of these kernels increases, and the kernels ?C ut ? and ?C tu ? are called strongly singular kernels whose integration must be evaluated in a Cauchy principal value sense and kernel ?C pp ? is called hypersingular kernels which must be evaluated in the Hadamard ?nite parts.
where qG is the boundary of the boundary G of the domain O of interest. Therefore, if the vector potential function Vi can be obtained, the dimension of the computational model can be reduced further by one, i.e. it is only necessary to evaluate values of Vi at nodes for 2-D BEM without integration, and integrals along the closed contours of psuedo-2-D elements for 3-D BEM, thus explaining the name BCM. Since the vector function Fi contains the unknown ?elds of displacements and tractions, special shape functions must be chosen for them to determine the potential function Vi : For 2-D linear elasticity problems, a ?ve node quadratic element with two traction nodes and three displacement nodes with equal intervals in between (Fig. 14a) is used to completely de?ne the shape functions. Three different kinds of 3-D elements were proposed for 3-D elasticity problems, and the simplest one is a four-node triangle element with three displacement nodes and one traction node at the centre (Fig. 14b). The values of vector function Vi will then be determined based on the shape functions (in closedform for 2-D problems and using numerical integrations for 3-D problems), which are then used to determine
306
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Displacement node (a)
Traction node (b)
Fig. 14. (a) Boundary elements used for 2-D BCM and (b) 3-D BCM (after Nagarajan et al., 1996) [176].
formulation. This problem will also occur when considering initial stress/strain effects, and non-linear material behaviour such as plastic deformation. The traditional technique to deal with such domain integrals is the division of the domain into a number of internal cells, which will seriously compromise the advantages of BEM’s ‘‘boundary only’’ discretization. Different techniques have been developed over the years to overcome this dif?culty (see for example Brebbia et al., 1984; Partridge et al., 1992) [189,190], as listed below.
*
nodal values of boundary displacements and tractions. Details for the numerical implementation of BCM can be seen in Nagarajan et al. (1994, 1996), Phan et al. (1997), Mukherjee et al. (1997) and Zhou et al. (1999) [175–179]. The method has also been combined with Galerkin approximation, to become the so-called Galerkin Boundary Contour Method (GBCM) (Zhou et al., 1998; Novati and Springhetti, 1999, [180,181]), and generated a symmetric ?nal coef?cient matrix when all elements are straight segments for 2-D problems. Similar developments, combined with meshless approaches, are also reported in Chati et al. (2001) [182]. The BCM approach is attractive mainly because of its further reduction of computational model dimensions and simpli?cation of the integration. The savings in preprocessing of the simulations are clear. Treatment of fractures and material non-homogeneity has not been studied in BCM; these may limit its applications to rock mechanics problems considering the present state-ofthe-art. 3.3.3.3. Boundary Node Method. Another recent development associated with the BEM is the Boundary Node Method (BNM) (Mukherjee and Mukherjee, 1997; Chati et al., 1999; Chati and Mukherjee, 2000; Kothnur et al., 1999) [183–186]. The method is a combination of traditional BEM with a meshless technique using the moving least squares for establishing trial functions without an explicit mesh of boundary elements. It further simpli?es the mesh generation tasks of the BEM at the cost of increasing computational operations for establishing the trial functions. Its applications concentrate on shape sensitivity analysis at present and solution of potential problems (Gu and Liu, 2002; Gowrishankar and Mukherjee, 2002) [187,188], but can be extended to general geomechanics problems, especially groundwater ?ow and stress/deformation analysis. 3.3.3.4. Dual Reciprocity Boundary Element Method (DRBEM). When source terms are presented in BEM formulations, such as gravitational body forces, heat sources, sinks/source terms in ?ow problems, thermal stress ?elds, etc., domain integrals often appear in the
*
*
*
Analytical integration of domain integrals, which is applicable to limited cases of simple geometry and boundary conditions. Fourier expansion of integrand functions, which is limited by domain geometry and boundary conditions, with additional computational cost of calculating the expansion coef?cients. Galerkin vector technique based on Green’s identity and higher order fundamental solutions, which can be readily applied to transform the domain integrals into boundary ones when the source terms are simple, such as constant body forces or heat sources. When the source terms are complex functions of space and time, the higher fundamental solutions may be dif?cult to obtain and the technique cannot be effectively used. An extension of this technique using multiple higher order fundamental solutions, instead of just one as in the Galerkin vector method, is called the Multiple Reciprocity Method in the literature. The Dual Reciprocity Method (DRM), which is a more generalized method closely related to the Galerkin vector and multiple reciprocity techniques (MDR) for constructing particular solutions suitable for non-linear and time-dependent problems. The source terms can be more generalized functions of space and time. Global interpolation functions, or Radial basis functions, are often used in converting the domain integrals into boundary ones (Cheng et al., 1994) [191].
The numerical treatment of the domain integrals, when the initial domain ?elds must be considered, is a signi?cant subject in BEM, since it relates to some important rock mechanics problems as mentioned above. The DRM approach appears to be the most widely applied technique so far since it has a uni?ed approach to treat different source terms and initial ?elds. A revisit of the technique for thermo-elasticity and elastic body force problems is presented in Cheng et al. (2001) [192] and an application for groundwater ?ow problem is presented in El Harrouni et al. (1997) [193]. Treatment of initial stress and strain ?elds related to plastic deformation is presented in Ochiai and Kobayashi (1999, 2001) and Gao (2002) for both 2-D and 3-D elasto-plastic problems [194–196]. Despite the
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
307
efforts, the numerical ef?ciency of the DRM is still a subject of debate and has signi?cant in?uence on the suitability and ef?ciency of BEM for problems of material non-homogeneity and non-linearity. The BEM methods as a whole is less ef?cient in treatment of material non-homogeneity and non-linearity compared with FEM and FVM, despite the fact that a limited number of multiple sub-regions of different material properties can be ef?ciently handled in BEM. However, it is much more ef?cient in simulating fracturing process (initiation, growth and coalescence) in elastic solids. It is therefore not surprising that BEM is most often applied for stress analysis problems, or fracturing process simulation, of CHILE continua, and is used for far-?eld representations in hybrid models. 3.4. Basic features of the Discrete Element Method (DEM) Rock mechanics is one of the disciplines from which the DEM originated (Burman, 1971; Cundall, 1971, 1974; Chappel, 1972, 1974; Byrne, 1974) [197–202]. The other engineering branches that stimulated the development of the DEM are structural analysis and multi-body systems. The theoretical foundation of the method is the formulation and solution of equations of motion of rigid and/or deformable bodies using implicit (based on FEM discretization) and explicit (using FVM discretization) formulations. The method has a broad variety of applications in rock mechanics, soil mechanics, structural analysis, granular materials, material processing, ?uid mechanics, multi-body systems, robot simulation, computer animation, etc. It is one of most rapidly developing areas of computational mechanics. The key concept of DEM is that the domain of interest is treated as an assemblage of rigid or deformable blocks/ particles/bodies and the contacts among them need to be identi?ed and continuously updated during the entire deformation/motion process, and represented by proper constitutive models. This fundamental conception leads naturally to three central issues: (i) identi?cation of block or particle system topology based on the fracture system geometry, or particle shape assumptions within the domain of interest; (ii) formulation and solution of equations of motion of the block (particle) system; (iii) detection and updating of varying contacts between the blocks (particles) as the consequences of motions and deformations of the discrete system. The basic difference between DEM and continuumbased methods is that the contact patterns between components of the system are continuously changing with the deformation process for the former, but are ?xed for the latter.
To formulate a DEM method to simulate the mechanical processes in rock mechanics applications, the following problems must be solved: (1) space sub-division and identi?cation of block system topology; (2) representation of block deformation (rigid or deformable, using FVM or FEM); (3) developing an algorithm for contact detection (penalty function, Lagrange multiplier, or augmented Lagrange multiplier); (4) obtaining constitutive equations for the rock blocks and fractures; (5) integration of the equations of motion of the blocks/particles (dynamic relaxation; time-marching FVM). The block system identi?cation depends on the reconstruction of the fracture system in situ according to usually very limited data from borehole logging or surface/underground mapping of fracture systems at generally very limited exposure areas. The usual procedure is to establish probabilistic density functions (PDFs) of the fracture parameters (orientation, frequency, size/trace length, etc.) and then use the random number ?eld technique to re-generate a number of realizations of synthetic fracture systems that share the same statistical geometric properties of the sampled fracture population from the logging and mapping. The reliability of such stochastic models of fracture system is therefore dependent on the quality of the logging and mapping operations that, in turn, depends on the quantity of data (therefore areas for fracture mapping and numbers/length of boreholes) available. It is evident that reliability of such fracture system models are largely unknown or the level of uncertainty is very high, due simply to the fact that the real fracture systems are hidden inside rock masses and will never be fully accessed by measurements. Regarding this dif?culty, a large number of such realizations need to be generated so that they, collectively, provide a much improved representation of the stochastic nature of the fracture system. This is called Monte Carlo simulation, and will naturally cause signi?cant increase of DEM computational cost. When the fracture systems are available, either regenerated or assumed, the next step is to construct the block systems de?ned by the fracture system. The task is trivial if the fractures are in?nitely long (or large) and follow regular patterns of distributions (such as constant spacing and ?xed orientations). However, for random fracture systems of ?nite fracture size, special techniques of combinatorial topology is needed to construct the block systems according to the fracture system geometry (see Lin et al., 1987; Lin and Fairhurst, 1988; Lin, 1992; Jing and Stephansson, 1994a, b; Jing, 2000; Lu, 2002) [203–209].
308
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
For rigid block analysis, an explicit time-marching scheme is used to solve the dynamic equations of motion of the rigid block system, based on a dynamic or static relaxation scheme, or an FDM approach in the time domain. For deformable block systems, the solution strategies are different for the treatment of block deformability. One is explicit solution with ?nite volume discretization of the block interiors, without the need for solving large-scale matrix equations. The other is an implicit solution with ?nite element discretization of the block interiors, which leads to a matrix equation representing the deformability of the block systems, similar to that of the FEM. The most representative explicit DEM methods is the Distinct Element Method created by Cundall (1980, 1988) [210,211] with the computer codes UDEC and 3DEC for 2- and 3-D problems of rock mechanics (ITSACA, 1992, 1994) [212,213]. Other developments were made in parallel with the distinct element approach and used the name ‘discrete element methods’, such as in Taylor (1983), Williams et al. (1985), Williams and Mustoe (1987), Williams (1988), Willaims and Pentland (1992), Mustoe (1992), Hocking (1977, 1992), Williams and O’Connor (1995) [214–222]. Another approach, socalled ‘Block-Spring Model’ (BSM) is essentially a version of DEM with rigid blocks linked by springs and applied for structural (Kawai, 1977a, b; Kawai et al., 1978) [223–225] and rock engineering problems (Wang and Garga, 1993; Wang et al., 1997; Li and Vance, 1999; Hu, 1997; and Li and Wang, 1998) [226–230]. However, the approach and codes by the Distinct Element Methods appears to be the main direction of application in rock mechanics problems, even the term ‘discrete element methods’ is more universally adopted. The implicit DEM was represented mainly by the Discontinuous Deformation Analysis (DDA) approach, originated by Shi (1988) [231] and further developed by Shyu (1993) and Chang (1994) [232,233] for stress/ deformation analysis, and Kim et al. (2000) and Jing et al. (2001) for coupled stress-?ow problems [234,235]. The method uses standard FEM meshes over blocks and the contacts are treated using the penalty method. Similar approaches were also developed by Ghaboussi (1988), Barbosa and Ghaboussi (1989, 1990) [236–238]. The technique uses four-noded blocks as the standard element and is called the Discrete Finite Element Method. Another similar development, called the combined ?nite-DEM (Munjiza et al., 1995, 1999; Munjiza and Andrews, 2000) [239–241], considers not only the block deformation but also fracturing and fragmentation of the rocks. However, in terms of development and application, the DDA approach occupies the front position. DDA has two advantages over the explicit DEM: permission for relatively larger time steps and closed-form integrations for the stiffness matrices of elements. An existing FEM code can also be
readily transformed into a DDA code while keeping all the advantageous features of the FEM. 3.4.1. Explicit DEM—Distinct Element Method: block systems The Distinct Element Method was originated in the early 70s by a landmark paper on the progressive movements of rock masses as 2-D rigid block assemblages (Cundall, 1971) [198]. The work was extended later into a code, RBM, written in machine language for a NOVA mini-computer (Cundall, 1974) [199]. The method and the RBM code later progressed, ?rstly by approximating the deformation of blocks of complex 2-D geometry by a constant strain tensor, with the code translated into the FORTRAN language and called SDEM for block systems (Cundall and Marti, 1979) [242]. A separate version of the SDEM code, called CRACK, was created to consider fracturing, cracking and splitting of intact blocks under loading, based on a tensile failure criterion. The representation of ‘‘simply deformable blocks’’ causes incompatibility between the complex block geometry and constant strain tensor, and the dif?culty was overcome later by using full internal discretization of blocks by ?nite volume meshes of triangle elements, leading to early versions of the code UDEC (Cundall, 1980; Cundall and Hart, 1985) [210,243], which has a BEM function representing the far-?eld (Lemos, 1987) [244]. Extension to 3-D problems was developed by Cundall (1988) [211] and Hart et al. (1988) [245], leading to the code 3DEC. The technique of the explicit DEM is presented comprehensively in Cundall and Hart (1992), Hart (1993) and Curran and Ofoegbu (1993) [246–248]. The principle of simulating large-scale deformations of elasto-plastic materials using ?nite difference/volume schemes developed in Wilkins (1963) [16] and the dynamic relaxation principles (Southwell, 1935, 1940, 1956) [249–251] are the mathematical basis. The contact detection and updating is performed based on the ‘‘contact overlap’’ concept. The method and codes were then developed further by coupling heat conduction and viscous ?uid ?ow through fractures (treated as interfaces between block boundaries). 3.4.1.1. Block discretization. Blocks are represented as convex polyhedra in 3-D with each face a planar convex polygon having a ?nite number of rectilinear edges. Their 2-D counterparts are general polygons with a ?nite number of straight edges (Fig. 15). The 2-D polygons can be either convex or concave, but the 3-D polyhedral must be convex. These blocks are formed by fractures which are represented in the problem domain either individually (for larger-scale fractures) or by a fracture sets generator (for smaller-scale fracture sets) using random distributions—based on site or modelling requirement data—of dip angles, dip directions, spacing
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
309
triangle element tetrahedral element a tetrahedral element
(a)
(b)
(c)
Fig. 15. Discretization of blocks by: (a) constant strain triangles; (b) constant strain tetrahedral; and (c) a typical tetrahedral element.
Table 1 Types of contacts for polygons and polyhedral Block shapes Arbitrary polygons (convex or concave) (2-D block) Convex polyhedral (3-D block) Contact types vertex-to-vertex, vertex-to-edge, edge-to-edge vertex-to-vertex, vertex-to-edge, vertex-to-face, edge-to-edge, edge-to-face, face-to-face
and apertures of the sets. The vertices (corners), edges and faces of individual blocks and their connection relations are identi?ed during the block generation process. The deformable blocks are further divided into a ?nite number of constant strain triangles in 2-D or tetrahedra in 3-D. These triangles or tetrahedra form a mesh of the FVM (zones). Rectangular element meshes can also be used for 2-D problems when the problem geometry is favourable. 3.4.1.2. Representation of deformation. An explicit, large strain Lagrangian formulation for the constant strain elements is used to represent the element deformations. The displacement ?eld of each element varies linearly and the faces or edges of the elements remain as planar surface or straight line segments. Higher order elements may also be used, but curved boundary surfaces (or edges) may be obtained, which may in turn complicate the contact-detection algorithm. Based on Gauss’ theorem to convert volume (area) integrals into surface (line) integrals, the increments of element strain can be written Deij E
N Dt X m ??v ?nj 7?vm ?ni ?DS k ; j 2 k?1 i
?23?
where DSk is the area (or length) of the kth boundary face (or edge) with unit normal nk ; and vm is the mean i i value of velocity over DSk : The summation extends over the N faces (or edges) of an element (zone). The sign ‘‘+’’ is used if i ? j; otherwise, the sign ‘‘?’’ is used. Dt is the time step. The stress increments are obtained by invoking the constitutive equations for the block materials.
3.4.1.3. Representation of contacts. Kinematically, block contacts are determined by the smallest distance between two blocks, pre-set in the codes or models. When this distance is within a prescribed threshold, a potential contact between these two blocks is numerically established. The contact-detection algorithm in the Distinct Element Method programs determines the contact type (different touching patterns between vertices, edges and faces), the maximum gap (if two blocks do not touch but are separated by a gap close to the pre-set tolerance), and the unit normal vector de?ning the tangential plane on which sliding can take place. Table 1 lists all types of contacts. Mechanically, the interaction between two contacting blocks is characterized by a stiffness (spring) in the normal direction and a stiffness and friction angle (spring-slip surface series) in the tangential directions with respect to the fracture surface (contact plane, see Fig. 16a). Interaction forces developed at contact points are determined as linear or non-linear functions of the deformations of springs and slip surfaces (i.e., the relative movements of blocks at contact points) and resolved into normal and tangential components, depending the constitutive models of the contacts (point contacts or edge/face contacts). The concept of contact ‘overlap’, though physically inadmissible in block kinematics—because blocks should not interpenetrate each other—may be accepted as a mathematical means to represent the deformability of the contacts. However, it does present a numerical shortcoming that is dif?cult to overcome when the normal forces or stresses at contact points are large. In this case, even with high normal stiffness, the ‘overlap’ may be too excessive to be acceptable and the calculation has to be stopped to implement some remedial measure (for example, to increase the normal
310
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
Kn
Φ
Kt
l
l
(a)
Old position New position C1 ?Un ?U t Ft Fn l1
(b)
For deformable blocks, the equations of motion are written for grid points—the vertices of internal difference elements. The central difference scheme is similar to the ?rst equation in Eq. (24) with only one modi?cation to the resultant out-of-balance force, fi fi ? fic ?
C3 C4
N X k?1
sij ?nk DS k ?; j
?26?
C2
l2 l3 l4
(c)
(d)
Fig. 16. Mechanical representation of contacts in the 2-D DEM.
stiffness) and start again. It also presents a problem for ?uid ?ow calculation in which the apertures of fractures may become negative if ‘overlap’ occurs at contact points. The mathematical representation of the contact ‘overlap’ is thus not fully compatible with physical reality. 3.4.1.4. Numerical integration of the equations of motion. An explicit central difference scheme is applied in the Distinct Element Method to integrate the equations of motion of the block system, as opposed to the implicit approach utilized in other continuumbased numerical methods. The unknown variables (contact forces or stresses) on the block boundary or in the internal elements are determined locally at each time step from the known variables on the boundaries, in the elements and their immediate neighbours. There is no need to set up and solve a matrix form of the equations of motion. The non-linearity in the material behaviour (of the fractures or intact blocks) can be handled in a straightforward manner. The equations of motion for a rigid block, in terms of translational and rotational velocities, are written P fi ?t?Dt=2? ?t?Dt=2? vi ? bi Dt; ? vi ? m P Mi ?t?Dt=2? ?t?Dt=2? Dt; ?24? ? oi ? oi I where m is the block mass, I is the moment of inertia, bi are the volume force components of the block and Mi are the components of the resultant moment. The displacement at the next time step is then given by u?t?Dt? ? u?t? ? vi i i
?t?Dt=2?
where fic is the resultant contact force if the grid point is on the boundary of the block. The symbol N denotes the number of difference elements connected by this grid point. At each time step, the kinematic quantities (velocities, displacements and accelerations) are ?rst calculated and the contact forces or stresses, as well as the internal stresses of the elements, are then obtained via constitutive relations for contacts. In the general calculation procedure, two basic tasks are performed in turn. The kinematic quantities are updated ?rst, followed by invoking the constitutive relations to provide the corresponding forces and stresses, see Fig. 17. 3.4.1.5. Applications and remarks. Due mainly to its conceptual attractions in the explicit representation of fractures, the DEM, especially the Distinct Element Method, has been enjoying wide application in rock engineering. A large quantity of associated publications has been published, especially in conference proceedings: it is not practical to list these even at a moderate level for this review. Therefore, a few representative references, mainly in international journals, are given here to show the wide range of the applicability of the methods:
*
*
*
*
*
*
Dt; Dt; ?25?
*
y?t?Dt? ? y?t? ? oi i i
?t?Dt=2?
where yi is the angular displacement of the block.
Tunnelling, underground excavations and mining: Barton (1991), Jing and Stephansson (1991), Nordlund et al. (1995), Chryssanthakis et al. (1997), Hanssen et al. (1993), Kochen and Andrade (1997), McNearny and Abel (1993), Souley et al. (1997a, b), So?anos and Kapenis (1998) and Lorig et al. (1995) [252–262]; Rock dynamics: Zhao et al. (1999) and Cai and Zhao (2000) [263,264]; Nuclear waste repository design and performance assessment: Chan et al. (1995), Hansson et al. (1995), Jing et al. (1995, 1997) [265–268]; Reservoir simulations: Gutierrez and Makurat (1997) [269]; Fluid injection: Harper and Last (1989, 1990a, b) [270–272]; Rock slopes, caving and gravity ?ow of particle systems: Zhu et al. (1999) [273]; Laboratory test simulations and constitutive model development for hard rocks: Jing et al. (1993, 1994), Lanaro et al. (1997) [274–276];
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
311
Kn Constitutive relations for all contacts Kt detail F n := F n - K n ?Un F t := F t - K Rigid blocks
t
?Ut ?Un Ft Fn
?U t Deformable blocks
F t := min{ ?Fn, abs(F t )}sgn(Ft )
All blocks Fc i Equations of motion Xi M
All blocks Fi
c
element
grid point
At block centroids F i = ΣF i M =Σ (e x F ) d Ui = Fi /m dt 2 2 d θ = M/I dt 2 etc. . . .
2 ij i j c
At elements (zones): ?ε = 1 2 ij d dU i + d dU j dx( dt ) dx( dt )
j i ij
?t
σ = C(σ , ? ε , ...) ij ij
at grid points: F i = Fie + F
2 c i
d U i = F /m i 2 dt etc. .......
t := t + ?t Back to
Fig. 17. Calculation cycles in the Distinct Element Method (after Hart, 1993) [247].
*
* * *
*
*
Stress-?ow coupling: Makurat et al. (1995) and Liao and Hencher (1997) [277,278]; Hard rock reinforcement: Lorig (1985) [279]; Intraplate earthquake: Jing (1990) [280]; Well and borehole stability: Rawlings et al. (1993) and Santarelli et al. (1992) [281,282]; Acoustic emission in rock: Hazzard and Young (2000) [283]; Derivation of equivalent hydro-mechanical properties of fractured rocks: Zhnag et al. (1996), Mas-Ivars et al. (2001), Min et al. (2001) [284–286].
A recent book by Sharma et al. (2001) [287] includes reference to a collection of DEM application papers for various aspects of rock engineering. The applications
concentrate on hard rock problems and have increasing focus on coupled hydro-mechanical behaviour—because of the dominating effects of the rock fractures on these aspects, and so where the explicit representation of fractures is necessary. For the softer and weaker rocks, equivalent continuum models are more applicable because there is less difference between the deformability of the fractures and the rock matrix. Despite the advantages of DEM, lack of knowledge of the geometry of the rock fractures limits its more general applications. In general, the geometry of fracture systems in rock masses cannot be known and can only be roughly estimated. The adequacy of the DEM results in capturing the rock reality are therefore highly dependent on the interpretation of the in situ fracture
312
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
system geometry—which cannot be even moderately validated in practice. Of course, the same problem applies also to the continuum models, such as the FEM or FDM, but the requirement for explicit fracture geometry representation in the DEM highlights the limitation and makes it more acute. Monte Carlo fracture simulation may help to reduce the level of uncertainty, albeit with increased computation. A prime subject for research, therefore, is increased quality of rock fracture system characterization with more advanced and affordable means, possibly using geophysical exploration techniques. 3.4.2. Implicit DEM—Discontinuous Deformation Analysis method: block systems DDA originated from a back analysis algorithm for determining a best ?t to a deformed con?guration of a block system from measured displacements and deformations (Shi and Goodman, 1985) [288]. It was later further developed to perform complete deformation analysis of a block system (Shi, 1988) [231]. The early formulation used a simple representation of block motion and deformation, with six basic variables (three rigid body motion and three constant strain components) and is not suitable for irregularly shaped blocks. The major improvements come from full internal discretization of blocks by triangular or fournoded FEM elements (Shyu, 1993; Chang, 1994) [232,233], as also demonstrated in Jing (1998) [289]. These improvements make the DDA method more suitable for arbitrarily shaped deformable blocks. Improvements of the method have been made for frictional contacts (Jing, 1993) [290], user-friendly code structure and environment (Chen et al., 1996; Chen, 1998) [291,292], rigid block systems (Koo and Chern, 1998) [293], fractured rock masses (Lin et al., 1996) [294] and coupled ?ow-stress analysis (Kim et al., 2000; Jing et al., 2001) [234,235] with ?uid conducted only in fractures. Numerous other extensions and improvements have been implemented over the years in the late 90s’, e.g. recently in Doolin and Sitar (2001) [295], with the bulk of the publications appearing in a series of ICADD conferences (Li et al., 1995; Salami and Banks, 1996; Ohnishi, 1997; Amadei, 1999) [81–84]. 3.4.2.1. Basic concepts of DDA. By the second law of thermodynamics, a mechanical system under loading (external and/or internal) must move or deform in a direction which produces the minimum total energy of the whole system. For a block system, the total energy consists of the potential energy due to different mechanisms like external loads, block deformation, system constraints, kinetic and strain energy of the blocks and the dissipated irreversible energy. The minimization of the system energy will produce an
equation of motion for the block system, the same as that used in the FEM. For a system of N blocks, each having mi nodes (i ? 1; 2; y; N), the total number of nodes is m1 ? m2 ? ? ? ? ? mN ? M; and each node has two orthogonal displacement variables, u and v: Assuming, without losing generality, that nodes are numbered sequentially blockwise, the minimization will yield (2M ? 2M) simultaneous equations, written symbolically as 38 9 2 k11 k12 k13 y k1N > d1 > > > > > 7> > 6 > > 6 k21 k22 k23 y k2N 7> d2 > 7< = 6 6 k31 k32 k33 y k3N 7 d3 7> > 6 7> > 6 ^ ^ ^ ^ 5> ^ > > > 4 ^ > > > > : ; kN1 kN2 kN3 y kNN dN 8 9 > f1 > > > > > > > > f2 > > > < = ? f3 or ?K?fDg ? fFg; ?27? > > > > > ^ > > > > > > > : ; fN where diagonal sub-matrices kij is a ?2mi ? 2mi ? matrix representing the sum of contributing sub-matrices for the ith block of mi nodes. Vector di is a ?2mi ? 1? vector of displacement variables of the ith block and vector f i is a ?2mi ? 1? vector of resultant general forces acting on the ith block. The off-diagonal sub-matrices kij ?iaj? represent the sum of contributing sub-matrices of contacts between blocks i and j and other inter-block actions like bolting. The matrix ?K? can also be called the global ‘‘stiffness matrix’’. For the three-block system illustrated in Fig. 18, the matrix structure of the equation system by DDA is shown in Fig. 19. Compared with the explicit approach of the DEM, the DDA method has four basic advantages over the explicit DEM: (i) The equilibrium condition is automatically satis?ed for quasi-static problems without using excessive iteration cycles. (ii) The length of the time step can be larger, and without inducing numerical instability. (iii) Closed-form integrations for the element and block stiffness matrices can be performed without the need for Gaussian quadrature techniques. (iv) It is easy to convert an existing FEM code into a DDA code and include many mature FEM techniques without inheriting the limitations of the ordinary FEM, such as small deformation, continuous material geometry, and reduced ef?ciency for dynamic analysis. However, matrix equations are produced and need to be solved, using the same FEM technique.
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
(9)
313
(3) (1)
(1) 1 2
(4) (5)
3
(6) (14)
(3) 4 5
(7) (13)
(8)
2, 3 Contact connectivity 2, 6 6, 2 5, 8
(1), (2), (3) (2), (4), (3) (5), (6), (9) (6), (8), (9) (6), (7), (8) (10), (11), (14) (11), (12), (14) (12), (13), (14)
(2)
(10)
6 7
(11)
8
(2)
(12)
(1) - Block number 1 - Element number
1 2 3 Element connectivity 4 5 6 7 8 (1) - Node number
Fig. 18. Blocks and FEM element discretizations by DDA—an example of three blocks.
Assembly the element deformation stiffness submatrices element 1 element 2 block 1
+
element 6 element 7
=
element 8 block 2
+
+
=
element 3
element 4
element 5
block 3
+
+
(a)
=
Assembly the contact stiffness submatrices between element elements 2 & 3 elements 2 & 6 elements 5 & 8
+
+
(b)
=
Assembly the global stiffness matrix
+
=
- one entry - sum of 2 entries - sum of 3 entries - sum of 4 entries - sum of 5 entries
global matrix
All deformation stiffness submatrices
All contact stiffness submatrices
(c)
Fig. 19. Matrix structure of the DDA method for the three-block system in Fig. 18.
314
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
The DDA method has emerged as an attractive model for geo-mechanical problems because its advantages cannot be replaced by continuum-based methods or explicit DEM formulations. It was also extended to handle 3-D block system analysis (Shi, 2001) [296] and use of higher order elements (Hsiung, 2001) [297], plus more comprehensive representation of the fractures (Zhang and Lu, 1998) [298]. The applications focus mainly on tunnelling, caverns, fracturing and fragmentation processes of geological and structural materials and earthquake effects (see, for example, Yeung and Leong, 1997; Hatzor and Benary, 1998; Ohnishi and Chen, 1999; Pearce et al., 2000; Hsiung and Shi, 2001) [299–303]. 3.4.3. Key Block theory The Key Block approach, initiated independently by Warburton (1983, 1993) [304,305] and Goodman and Shi (1985) [306], with a more rigorous topological treatment of block system geometry in the latter (see also Shi and Goodman, 1989, 1990 [307,308]), is a special method for stability analysis of rock structures dominated by the geometrical characteristics of the rock blocks and hence fracture systems. It does not perform any stress and deformation analysis, but identi?es the ‘‘key blocks’’ (or ‘‘keystones’’ in the terms of Warburton (1983) [304]), which are formed by intersecting fractures and excavated free surfaces in rocks and have the potential for sliding and rotation towards certain directions without geometrical constraints. This technique is a powerful and ef?cient tool for stability analysis and support design for slopes and underground excavations in fractured rocks—with large-scale movements of isolated blocks (such as wedges in slopes) as the major mode of failure and instability, rather than deformation and stresses of the intact rock matrix. It is therefore suitable for ‘hard rock engineering’, but is less suitable for soft rocks because the stress/deformation/failure of the rock matrix is equally important for the latter. Key Block Theory, or simply Block Theory in the Warburton approach, and the associated code development enjoys wide applications in rock engineering, with further development considering Monte Carlo simulations and probabilistic predictions (Stone and Young, 1994; Mauldon, 1993; Hatzor, 1993; Jakubowski and Tajdus, 1995; Kusznaul and Goodman, 1995) [309–313], water effects (Karaca and Goodman, 1993) [314], linear programming (Mauldon et al., 1997) [315], ?nite block size effect (Windsor, 1995) [316] and secondary blocks (Wibowo, 1997) [317]. Predictably, the major applications are in the ?eld of tunnel and slope stability analysis, such as reported by Chern and Wang (1993), Scott and Kottenstette (1993), Nishigaki and Miki (1995), Yow (1990), Boyle and Vogt (1995), Lee and Song (1998) and Lee and Park (2000) [318–324].
3.4.4. DEM formulations for particle systems Simulating the mechanical behaviour of granular materials is another important application area of the DEM, for both the Distinct Element and DDA approaches. The principle of the DEM technique for granular materials is basically the same as for the blocks, with the additional simpli?cation that particles are rigid and their shape can be regular (circular, elliptical in 2-D and ellipsoidal in 3-D) or irregular (generally polygonal in 2-D and polyhedral in 3-D). The contacts between the particles are represented by springs, and friction may also be considered. The seminal work of the DEM for granular materials for geomechanics and civil engineering application is the series of papers by Cundall and Strack (1979a–c, 1982) [325–328], which was based on an earlier work by Cundall et al. (1978) and Strack and Cundall (1978) [329,330]. The development and applications are mostly reported in a series of proceedings of symposia and conferences, such as in Jenkins and Satake (1983), ! Satake and Jenkins (1988), Biarez and Gourves (1989), Thornton (1993), Siriwardane and Zaman (1994) in the ?eld of micro-mechanics of granular media in general [331–335], and in Mustoe et al. (1989) and Williams and Mustoe (1993) in the ?eld of geomechanics in particular [336,337]. The simulations of particle systems using the DDA approach appear mostly in ICADD conference systems (Li et al., 1995; Salami and Banks, 1996; Ohnishi, 1997; Amadei, 1999) [81–84]. The most well-known codes in this ?eld are the PFC codes for both 2-D and 3-D problems (Itasca, 1995) [338], and the DMC code by Taylor and Preece (1989, 1990) [339,340]. The method has been widely applied to many different ?elds such as soil mechanics, the processing industry, non-metal material sciences and defence research. The following examples of publications highlight the wide applications of the DEM for particle systems in the ?eld of rock engineering:
*
*
*
* *
Fracturing and fragmentation processes of rock blasting: Preece (1990, 1994), Preece and Knudsen (1992), Preece et al. (1993), Preece and Scovira ! (1994), Donze et al. (1997), Lee et al. (1997), Lin and Ng (1994) [341–348]; Ground collapse and movements: Iwashita et al. (1988), Zhai et al. (1997) [349,350]; Hydraulic fracturing in rocks: Thallak et al. (1991), Huang and Kim (1993), Kim and Yao (1994) [351– 353]; Tunnelling: Kiyama et al. (1991)[354]; Rock fracture: Blair and Cook (1992) [355].
3.4.5. Dynamic lattice network models A numerical method closely related to the DEM model of granular material, often called the dynamic lattice network model, has also been applied for
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
315
simulating fracture initiation and propagation in rocks. The main concept of the lattice model is that the medium consists of a mesh of regular elements, such as triangle elements, with particles of lumped masses located at the vertices of the mesh. The mass particles are then connected by massless springs along the edges of the mesh. The dynamic motion of the medium is simulated by the equations of motions of the mass particles and the deformation of the springs, whose stiffness and strength are derived from those of the medium, which may have local variations and be generated randomly. The masses of the particles are derived from the density of the material, and can also be generated randomly, for representing statistical inhomogeneity of the medium. This technique is similar to that of the DEM for particle systems, with the difference that it represents the continuum behaviour of the medium through a deterministic/stochastic assembly of particles and springs, rather than as a direct discrete medium. The applications of the lattice model focus on the fracturing processes of intact rock materials during loading and hydraulic fracturing, such as reported in Paterson (1988), Song and Kim (1994a, 1994b, 1995), Muhlhaus et al. (1997), Schlangen and Van Mier (1994), . Li et al. (2000), Napier and Dede (1997) and Place and Mora (2000) [356–364]. 3.5. Discrete Fracture Network method 3.5.1. DFN-the basic concepts The DFN method is a special discrete model that considers ?uid ?ow and transport processes in fractured rock masses through a system of connected fractures. The technique was created in the early 1980s for both 2-D and 3-D problems (Long et al., 1982, 1985; Robinson, 1984; Andersson, 1984; Endo, 1984; Endo et al., 1984; Smith and Schwartz, 1984; Elsworth, 1986a, b; Dershowtz and Einstein, 1987; Andersson and Dverstop, 1987) [365–375], and has been continuously developed afterwards with many applications in civil, environmental and reservoir engineering and other geoscience ?elds. The effects of mechanical deformation and heat transfer in a rock mass on ?uid ?ow and transport are dif?cult to model by the DFN approach and are perhaps crudely approximated without explicit representation of the fractures. Thus, this method is most useful for the study of ?ow and transport in fractured media in which an equivalent continuum model is dif?cult to establish, and for the derivation of equivalent continuum ?ow and transport properties of fractured rocks (Yu et al., 1999; Zimmerman and Bodvasson, 1996) [376,377]. A large number of publications has reported progresses in journals and international symposia and conferences. Systematic presentations and evaluations of the method
have also appeared in books, such as Bear et al. (1993), Sahimi (1995), the US National Research Council (1996) and Adler and Thovert (1999) [378–381]. The DFN model is established on the understanding and representation of the two key factors: fracture system geometry and transmissivity of individual fractures. The former is based on stochastic simulations of fracture systems, using the PDFs of fracture parameters (density, orientation, size, aperture or transmissivity) formulated according to ?eld mapping results, in addition to the assumption about fracture shape (circular, elliptical or generally polygonal). Due to the fact that direct mapping can only be conducted at surface exposures of limited area, boreholes of limited length/depth and the walls of underground excavations (tunnels, caverns, shafts, etc.) of more limited exposure areas, and with both lower and upper cut-off limits for mapping, the reliability of fracture network information is dependent on the quality of mapping and representativeness of the sampling, and hence its adequacy is dif?cult to evaluate. Equally dif?cult is also the representation of the transmissivity of the fracture population, due to the fact that in situ and laboratory tests can only be performed at a limited number of fracture samples at restricted locations, and the effect of sample scale is dif?cult to determine. Despite the above hurdles, the DFN model enjoys wide applications for problems of fractured rocks, perhaps mainly due to the fact that it is a so far, irreplaceable tool for modelling ?uid ?ow and transport phenomena at the ‘near-?eld’ scale—the ‘near-?eld’ because the dominance of the fracture geometry at small and moderate scales makes the volume averaging principle used in continuum approximations sometimes unacceptable at such scales. Its applicability diminishes for ‘far-?eld’ problems at large scales when explicit representation of large numbers of fractures make the computational model less ef?cient, and the continuum model with equivalent properties become more attractive, similar to the DEM. There are many different DFN formulations and computer codes, but most notable are the approaches and codes FRACMAN/MAFIC (Dershowitz et al., 1993) [382] and NAPSAC (Stratford et al., 1990; Herbert, 1994, 1996; Wilcock, 1996) [383–386] with many applications for rock engineering projects over the years. Some special features of DFN are brie?y reviewed below. 3.5.2. Stochastic simulations of fracture networks The stochastic simulation of fracture systems is the geometric basis of the DFN approach and plays a crucial role in the performance and reliability of the DFN model, in the same way as the DEM. The key process is to create PDFs of fracture parameters relating to the densities, orientations and sizes, based on ?eld
316
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
mapping results using borehole logging data and scanline or window mapping techniques, and generate the realizations of the fractures systems according to these PDFs and assumptions about fracture shape (circular discs, ellipses or polygons), (Dershowitz, 1984; Billaux et al., 1989) [387–388]. A critical issue in this technique is the treatment of bias in estimation of the fracture densities and trace lengths from conventional straight scanline or rectangular window mappings. A notable recent development using circular windows (Mauldon, 1998; Mauldon et al., 2001) [389,390] provides an important step forward in this regard. 3.5.3. Solution of the ?ow ?elds within fractures Numerical techniques have been developed for the solution of ?ow ?elds for individual fracture elements using closed-form solutions, the ?nite element model, the boundary element model, the pipe model and the channel lattice model. Closed-form solutions exist, at present, only for planar, smooth fractures with parallel surfaces of regular shape (i.e. circular or rectangular discs) for steady-state ?ow (Long, 1983) [391] or for both steadystate and transient ?ow (Amdei and Illangaseekare, 1992) [392]. For fractures with general shapes, numerical solutions must be used. The FEM discretization technique is perhaps the most well-known techniques used in the DFN ?ow models and has been used in the DFN codes FRACMAN/MAFIC and NAPSAC. The basic concept is to impose an FEM mesh over the individual discs representing fractures in space (Fig. 20a) and solve the ?ow equations. The aperture or transmissivity ?eld within the fracture can be either constant or randomly distributed. Similarly, the BEM discretization can also be applied with the boundary elements de?ned only on the disc boundaries (Fig. 20b),
with the fracture intersections treated as internal boundaries in the BEM solution. The compatibility condition is imposed at the intersections of discs. See Elsworth (1986a, b) [372,373] and Robinson (1986) [393] for detailed formulations. The pipe model represents a fracture as a pipe of equivalent hydraulic conductivity starting at the disc centre and ending at the intersections with other fractures (Fig. 20c), based on the fracture transmissivity, size and shape distributions (Cacas, 1990) [394]. The channel lattice model represents the whole fracture by a network of regular pipe networks (Fig. 3.13d). The pipe model leads to a simpler representation of the fracture system geometry, but may have dif?culties to properly represent systems of a number of large fractures. The channel lattice model is more suitable for simulating the complex ?ow behaviour inside the fractures, such as the ‘‘channel ?ow’’ phenomena (Tsang and Tsang, 1987) [395], and is computationally less demanding than the FEM and BEM models since the solutions of the ?ow ?elds through the pipe elements are analytical. The fractal concept has also been applied to the DFN approach in order to consider the scale dependence of the fracture system geometry and for upscaling the permeability properties, using usually the full box dimensions or the Cantor dust model (Barton and ! Larsen, 1985; Chiles, 1988; Barton, 1992) [396–398]. Power law relations have been also found to exist for trace lengths of fractures and have been applied for representing fracture system connectivity (Renshaw, 1999) [399]. 3.5.4. Issues of importance and dif?culty The in?uence of the rock matrix on ?ow in rock fractures is usually not considered in the DFN models. However, the related effects also need to be estimated
FE elements
Intersection BE elements
Intersection
(a)
(b)
Equivalent pipes
Equivalent pipe networks
(c)
(d)
Fig. 20. Representation of rock fractures for the ?ow equation solution: (a) FEM; (b) BEM; (c) equivalent pipes; and (d) channel lattice model.
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
317
when the permeability of the rock matrix is large compared with that of the fractures, or the time scale of the problem is long enough so that matrix diffusion cannot be ignored. In such cases, a fracture system embedded in a highly porous matrix needs to be properly represented, considering the different time scales between the fracture ?ow and matrix diffusion processes, such as the FEM technique used by Sudicky and McLaren (1992) [400]. Dershowitz and Miler (1995) [401] reported a simpli?ed technique for considering matrix in?uence on ?ow in DFN models, using a probabilistic particle tracking technique. The stress/deformation processes of the fractures, and the effects of stress/deformation of rock matrix on the deformation/?ow of fractures, are usually also ignored in the DFN models, simply due to the complexity of the numerical techniques and extraordinary computational efforts needed. Although simple estimates concerning the effects of in situ stresses on fracture aperture variations have been used in DFN models, e.g. in the NAPSAC code (Wilcock, 1996) [386], proper representation of the coupled stress/?ow process in fractures is needed, especially for near-?eld problems for performance and safety assessment of radioactive waste repositories. Simulation of multiphase ?ow and transport by DFN models is also an important subject, with signi?cance for reservoir simulations (for oil/gas recovery and hot-dryrock thermal energy extraction projects) and nuclear waster repositories (heat decay and waster phase change, and radionuclide transport). The conventional technique is the dual permeability models which treats the fracture system and matrix as two superimposed conducting media for ?ow, heat transfer and transport, because the heat transfer process is largely dominated by matrix heat conduction—but heat convection by ?ow in fractures may also become important issues, e.g. for hotdry-rock reservoir simulations and performance of the near-?eld buffer material for nuclear waste isolation (Pruess and Wang, 1987; Slough et al., 1999; Hughes and Blunt, 2001) [402–404]. Proper solutions of the governing equations for mass and energy balances, with phase change relations (water evaporation and condensation, moisture ?ow driven by thermal gradients, etc.) need to be properly considered, with embedded DFN models in porous media, which is a subject still needed to be further investigated. The main challenge is the computational effort for considering the matrix– fracture interaction, with fractures dominating ?ow and matrix dominating heat transfer. Phase changes may occur in both. Because of its origin based on fracture mapping at limited exposure areas, the DFN models are based on, and at the same time are limited by, their two cornerstones: the largely unknown true fracture geometry; and the hydraulic properties of fractures—with, in
both cases, the attendant dif?culties in property measurement and evaluation. The upper and lower cut-off limits for trace lengths in mapping, and the effects of scale, roughness, and stress-dependency of fracture aperture and transmissivity, play a signi?cant role in the adequacy and reliability of the DFN models: it is often dif?cult to evaluate their signi?cance because of the lack of any independent checking mechanisms. An additional drawback is the high computational demands when large fracture numbers are required for large-scale problems. These drawbacks are similar to those of the DEM models and can only be partially overcome or qualitatively addressed with the aid of advanced mapping and measurement techniques. Therefore, the applications of DFN models are concentrated more on characterization of the permeability of fractured rocks and generic studies of fracture in?uences, and the design of rock engineering works for near-?eld problems, see for example as follows:
*
*
*
Hot-dry-rock reservoir simulations: Layton et al. (1992), Ezzedine and de Marsily (1993), Watanabe and Takahashi (1995), Kolditz (1995), WillisRichards and Wallroth (1995), Willis-Richards (1995) [405–410]; Characterization of permeability of fractured rocks: Dershowitz et al. (1992), Herbert and Layton (1995), !! Doe and Wallmann (1995), Barthelemy et al. (1996), Jing and Stephansson (1996), Margolin et al. (1998), Mazurek et al. (1998), Zhang and Sanderson (1999) [411–418]; Water effects on underground excavations and rock slopes: Rouleau and Gale (1987), Xu and Cojean (1990), He (1997) [419–421].
3.5.5. Alternative formulations—percolation theory Percolation theory is a counterpart of the lattice model for solid deformation and fracturing for ?uid diffusion process. The theory is based on a random lattice model of conductors (fractures) for deriving ?uid transport conditions (permeability), based on the connectivity through a geometrical sample of the fractures (Robinson, 1984; Hestir and Long, 1990; Berkowitz and Balberg, 1993; Sahimi, 1995) [367,422,423,379]. The theory provides a theoretical platform for understanding the geometric conditions for ?uid conduction in fractured media in terms of a percolation threshold, expressed as a critical probability. Application of this theory in rock mechanics and engineering problems seems to concentrate on characterization of ?ow properties of fractured rocks. Some example works in this ?eld are given bellow:
*
Critical behaviour of deformation and permeability of fractured rocks: Zhang and Sanderson (1998) [424];
318
*
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
*
*
Flow and transport in fracture networks: Mo et al. (1998) [425]; Microstructure and physical properties of rock: ! Guequen et al. (1997) [426]; Fluid, heat and mass transport in percolation clusters: Kimmich et al. (2001) [427].
affected by the symmetrized BEM equation. A possible step forward in this direction is to use the Galerkin double integration techniques in the BEM region so that the ?nal BEM stiffness matrix is automatically symmetric, and therefore can be directly inserted in the ?nal hybrid BEM/FEM matrix without errors caused by arti?cial ‘‘symmetrization’’. 3.6.2. Hybrid DEM/BEM The hybrid DEM/BEM model was implemented only for the explicit Distinct Element Method, in the code group of UDEC and 3DEC. The technique was created by Lorig and Brady (1982, 1984, 1986) [436,3,437], and was implemented into UDEC by Lemos (1987) [244]. The basic concept is to treat the BEM region (which surrounds the DEM region) as a ‘super’ block having contacts with smaller blocks along the interfaces with the DEM region (cf. Fig. 5), which can be treated in standard DEM contact representations. The key conditions are: (1) the kinematic continuity along the interfaces of the two regions during the time-marching process; and (2) the elastic properties of the two regions near the interface are similar. Condition 2 indicates that blocks in the DEM regions must be deformable, i.e. not be rigid blocks. In the case of mixed rigid and deformable block systems, special equations of motion need to be developed to handle such cases. Wei (1992) and Wei and Hudson (1998) [438,439] reported a development of hybrid discrete-continuum models for coupled hydro-mechanical analysis of fractured rocks, using combinations of DEM, DFN and BEM approaches. The near-?eld of a fractured rock mass is simulated using DEM and DFN models, using independent DFN and DEM codes, for representing the dominance of fractures on the near-?eld ?uid ?ow and stress/deformation of rock blocks, and the far-?eld is simulated by BEM codes for ?ow and stress/deformation in a continuum. The equations of ?ow and motion are not directly coupled, but are solved independently by separate DFN, DEM and BEM codes, and they are coupled through an internal linking algorithm with the time-marching process. 3.6.3. Alternative hybrid models Besides the above mainstream hybrid formulations, there are other coupling techniques which take advantage of different numerical methods. Pan and Reed (1991) [440] reported a hybrid DEM/FEM model, in which the DEM region consists of rigid blocks and the FEM region can have non-linear material behaviour. The algorithm places the FEM calculations into the DEM time-marching process. Since the blocks in DEM
3.6. Hybrid models Hybrid models are frequently used in rock engineering, basically for ?ow and stress/deformation problems of fractured rocks. The main types of hybrid models are the hybrid BEM/FEM, DEM/BEM models. The hybrid DEM/FEM models are also developed. The BEM is most commonly used for simulating far-?eld rocks as an equivalent elastic continuum, and the FEM and DEM for the non-linear or fractured near-?eld where explicit representation of fractures or non-linear mechanical behaviour, such as plasticity, is needed. This harmonizes the geometry of the required problem resolution with the numerical techniques available, thus providing an effective representation of the effects of the far-?eld to the near-?eld rocks. 3.6.1. Hybrid FEM/BEM models The hybrid FEM/BEM was ?rst proposed in Zienkiewicz et al. (1977) [428], then followed by Brady and Wassyng (1981) [429], and Beer (1983) [430] as a general stress analysis technique. In rock mechanics, it has been used mainly for simulating the mechanical behaviour of underground excavations, as reported in Varadarajan et al. (1985), Ohkami et al. (1985), Gioda and Carini (1985), Swoboda et al. (1987) and van Estorff and Firuziaan (2000) [431–435]. The coupling algorithms are also presented in detail in Beer and Watson (1992) [49]. The standard technique is to treat the BEM region as a ‘super’ element with an arti?cially ‘symmetrized’ stiffness matrix, using the least-square techniques, so that it can be easily inserted into the symmetric FEM stiffness matrix for the ?nal solution, which is easier to handle than the non-symmetric BEM stiffness matrix. However, such arti?cial ‘symmetrization’ introduces additional errors into the ?nal system equations. The coupling can also be performed in the opposite direction, i.e. treat the FEM region as a ‘super’ BEM element, and insert the corresponding FEM stiffness matrix into the ?nal BEM stiffness matrix; this leads to a asymmetric stiffness matrix for the ?nal equation, which needs additional computational efforts for solution. The hybrid BEM/FEM models are as ef?cient computationally as the FEM, with the additional advantage of being able to deal with the non-linear behaviour of materials in the FEM region, using the FEM’s advantages. However, this advantage may be
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
319
region are rigid and the FEM region is an elastic continuum, the kinematical continuity condition along the interface of the DEM and FEM regions may not be satis?ed. . A hybrid beam-BEM model was reported by Pottler and Swoboda (1986) [441] to simulate the support behaviour of underground openings, using the same principle as the hybrid BEM/FEM model. Sugawara et al. (1988) [442] reported a hybrid BEMcharacteristics method for non-linear analysis of rock cavern. 3.7. Neural networks/empirical techniques All the numerical modelling methods described so far are in the category of ‘1:1 mapping’—using the terms in Fig. 2, Section 1, which illustrates the eight basic methods of rock mechanics modelling and rock engineering design. The term ‘1:1 mapping’ refers to the attempt to model geometry and physical mechanisms directly, either speci?cally or through equivalent properties. A completely different numerical method, located in Box 2C in Fig. 2, uses neural networks: a ‘non-1:1 mapping’ method. The rock mass is represented indirectly by a system of connected nodes, but there is not necessarily any physical interpretation of the nodes, nor of their input and output values. The model purports to operate like the human brain. Another example is the use of empirical techniques, the rock rating systems (Q, RMR, RMI, GSI, etc.—see Box 2B in Fig. 2), for characterization of rock mass properties and construction design, without resorting to the solution of the basic balance equations and use of thermodynamically acceptable constitutive models. This design method is possible because of the fact that a true 1:1 mapping solution can never be achieved with 100% con?dence and reliability for problems in fractured rocks because a 100% characterization of the rock/ fracture system can never be achieved and validated. This is the reason for the development and success of non-1:1 mapping techniques like the rating systems in rock engineering. Through the use of empirical techniques, including empirical failure criteria, a rock engineer can state whether a tunnel or a cavern is safe or dangerous, without solving any of the equations in the BEM, FEM or DEM methods. Such a ‘non-1:1 mapping’ system has its advantages and disadvantages. The advantages are:
*
*
there is the possibility that the ‘perception’ we enjoy in the human brain may be mimicked in the neural network, so that the programs can incorporate judgements based on empirical methods and experiences.
The disadvantages are that
*
*
*
*
the procedure may be regarded as simply supercomplicated curve ?tting (because the program has to be ‘taught’), the model cannot reliably estimate outside its range of training parameters, critical mechanisms might be omitted in the model training, and there is a lack of theoretical basis for veri?cation and validation of the techniques and their outcomes.
Neural network models provide descriptive and predictive capabilities and, for this reason, have been applied through the range of rock parameter identi?cation and engineering activities. Recent published works on the application of neural networks to rock mechanics and rock engineering includes the following:
*
*
* *
*
*
* *
*
*
*
*
*
*
*
that the geometrical and physical constraints of the problem, which appear in governing equations and constitutive laws when the 1:1 mapping techniques are used, are no longer so dominating, different kinds of neural networks and empirical models can be applied to a problem, and
*
* *
Stress–strain curve for intact rock: Millar and Clarici (1994) [443]; Intact rock strength: Alvarez Grima and Babuska (1999); Singh et al. (2001) [444,445]; Fracture aperture: Kacewicz (1994) [446]; Shear behaviour of fractures: Lessard and Hadjigeorgiou (1999) [447]; Rock fracture analysis: Sirat and Talbot (2001) [448]; Rock mass properties: Qiao et al. (2000); Feng et al. (2000) [449,450]; Rock mechanics models: Feng and Seto (1999) [451]; Rock mass classi?cation: Sklavounos and Sakellariou (1995); Liu and Wang (1999) [452,453]; Displacements of rock slopes: Deng and Lee (2001) [454]; Tunnel boring machine performance: Alvarez Grima et al. (2000) [455]; Displacements/failure in tunnels: Sellner and Steindorfer (2000); Lee and Sterling (1992) [456,457]; Tunnel support: Leu (2001); Leu et al. (2001) [458,459]; Surface settlement due to tunnelling: Kim et al. (2001) [460] Earthquake information analysis: Feng et al. (1997) [461]; Rock engineering systems (RES) modelling: Millar and Hudson (1994); Yang and Zhang (1998) [462,463]; Rock engineering: Yi and Wanstedt (1998) [464]; Overview of the subject: Hudson and Hudson (1997) [465].
320
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
As evidenced by the list of highlighted references above, the neural network modelling approach has already been applied to the variety of subjects in rock mechanics and rock engineering. It is also evident that the method has signi?cant potential—because of its ‘non-1:1 mapping’ character and because it may be possible in the future for such networks to include creative ability, perception and judgement. However, the method has not yet provided an alternative to conventional modelling, and it may be a long time before it can be used in the comprehensive Box 2D mode envisaged in Fig. 2 and described in Feng and Hudson (2003) [466]. 3.8. Constitutive models of rocks The constitutive models of rocks, including those for both rock fractures and fractured rock masses, are one of the most important components of numerical solutions for practical rock engineering problems, and one of the most intensively and continuously investigated subjects in rock mechanics. The most recent developments in the area are brie?y introduced with some comments, supported by literature sources. To make the presentation clearer, the models are divided into ?ve groups according to their different formulation platforms and traditional application areas: classical constitutive models, failure criteria, time effects and viscosity, size effects and homogenization, damage mechanics models, and rock fracture models. 3.8.1. Classical constitutive models of rocks The classical constitutive models are the models based mostly on the theory of elasticity and plasticity, but with special considerations of fracture effects. The model of linear elasticity based on the generalized Hooke’s law is still by far the most widely adopted assumption for the mechanical behaviour of rocks, especially for hard rocks. When the CHILE assumption (Continuous, Homogeneous, Isotropic, Linear Elastic) is adopted, the constitutive law is simply characterized by two independent material properties, either the most commonly known, Young’s modulus (E) and Poisson’s ratio (n), or Lame’s two parameters, G and l: More sophisticated constitutive models of anisotropic elasticity can be derived in closed-form by considering alternative elastic symmetry conditions for intact rocks (such as transversely isotropic elasticity) or equivalent or effective continuum elastic rocks intersected by orthogonal sets of in?nitely large or ?nite fractures, see Singh (1973), Gerrard (1982), Fossum (1985), Wei and Hudson (1986), Yashinaka and Yamabe (1986), Wu (1988), Murakami and Hegemier (1989), Chen (1989), Singh (2000), and Wu and Wang (2001) [467–476]. Oda (1986) [477] developed a crack tensor approach to derive
the equivalent elastic compliance tensors of fractured rocks with randomly distributed ?nite fractures, and Li et al. (1998) developed effective stress–strain relations for rocks containing deformable micro-cracks under compression. In the above developments, interactions between the fractures were not considered. A comprehensive summary of the constitutive laws for geomaterials in oil and gas reservoir rocks is given by Papamichos (1999) [478]. Plasticity and elasto-plasticity models have been developed and widely applied to fractured rocks since the 1970s, based mainly on the classical theory of plasticity, with typical models using Mohr–Coulomb and Hoek–Brown failure criteria (Hoek, 1983; Hoek and Brown, 1982, 1997) [479,132,480] as the yield functions and plastic potentials. A comprehensive and detailed description of such plasticity models is given in Owen and Hinton (1980) [45], together with FEM implementations. Parallel and similar developments can also be seen in Zheng et al. (1986) [481] with a strain space formulation, Adhikary and Dyskin (1998) [482] for layered rocks, Sulem et al. (1999) [483] for elasto-plasticity models of a sandstone and Boulon and Alachaher (1995) [484] for a non-linear model based on generalized strain paths, suitable for FEM implementations. Strain-hardening and strain-softening are the two main features of plastic behaviour of rocks, with the latter more often observed under uniaxial compression test conditions. Related works are reported ! by Dragon and Mroz (1979), Nemat-Nasser (1983), ! Zienkiewicz and Mroz (1984), Read and Hegemier (1984), Gerrard and Pande (1985), Desai and Salami (1987), Rowshandel and Nemat-Nasser (1987), Kim and Lade (1988), and Sterpi (1999) [485–493]. Strain-localization (e.g. shear-banding) is a deformation phenomenon closely related to the constitutive models of rocks, and has been studied intensively using plasticity and damage mechanics models (see, for example, Rudnicki and Rice, 1975; Fang, 2001) [494,21]. 3.8.2. Failure criteria The failure criteria of rocks are important components of constitutive relations and are usually used as yield surfaces or/and plastic potential functions in a plasticity model. Besides the most well known and perhaps also the most widely used Mohr–Coulomb and Hoek–Brown criteria, different failure (or strength) criteria have been proposed for rock masses over the years. The dimensionless forms of the Mohr-type failure criteria were discussed in Pariseau (1994) [495]. Sheorey (1997), Mostyn and Douglas (2000) and Parry (2000) have provided a comprehensive reviews of the subject [496–498]. Some of the recent developments are listed
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
321
below with their contexts:
*
*
*
*
*
*
*
*
*
*
* *
Anisotropy of jointed rock mass strength: Amadei and Savage (1989) [499]; Time-dependent tensile strength of saturated granite: Sun and Hu (1997) [500]; Failure criteria for anisotropic rocks: Duvean et al. (1998) [501]; . Compressive failure of rocks: Gupta and Bergstrom (1997), Tharp (1997) [502,503]; Testing factors for establishing rock strength: Hawkins (1998) [504]; Shear failure envelope for modi?ed Hoek–Brown criterion: Kumar (1998) [505]; Effect of intermediate principal stress on strength of anisotropic rocks: Singh et al. (1998) [506]; Modi?ed Mohr–Coulomb failure criterion for layered rocks: Lai et al. (1999) [507]; Short- and long-term strength of isotropic rocks: Aubertin et al. (2000) [508]; Relation between failure criterion and deformability of rock: Hibino et al. (2000) [509]; Water effects on rock strength: Masuda (2001) [510]; Failure criterion for transversely isotropic rocks: Tien and Kuo (2001) [511].
plasticity or plasticity (leading to constitutive models of so-called visco-elasticity, and its plasticity counterparts, visco-elasto-plasticity/visco-plasticity). Comprehensive descriptions of the constitutive models using viscoplasticity were given in Valanis (1976) [514] and in Owen and Hinton (1980) [45]. The link to and coupling with ?uid ?ow are also important. Representative and recent works are listed below:
*
*
*
*
*
*
*
*
3.8.3. Time effects and viscosity Time effects are one of the most important and also perhaps one of the least understood aspects of the physical behaviour of rock masses. There are two main aspects: the effects caused by the rock (and fracture) viscosity and the effects caused by the dynamic loading conditions. The former concerns mainly the stationary behaviour of rock over long- or extremely long-terms, such as geological time periods, and the latter is just the opposite—dynamic and even violent behaviour over short durations, such as earthquake effects. In such cases, not only the magnitudes, directions and durations, but also the rate of change in loading parameters are important. The time and rate effects are therefore often discussed in combination. The effect of viscosity, also termed a rock rheology effect, is signi?cant in rock salt and other weak rocks, and is caused by two mechanisms: creep and relaxation. The former is the behaviour of increasing deformation (strain) under constant loading (stress); and the latter is the decreasing loading (stress) while the deformation (strain) state is kept constant. The fundamentals of the physics, experimental basis and the constitutive models for these two mechanisms of viscosity were described in detail in Jaeger and Cook (1969) [512] and Cristesu and Hunsche (1998) [513], with applications in tunnel/cavern failure and borehole closure in mining and petroleum engineering. The effect of viscosity has been considered in constitutive modelling in combination with other basic deformation mechanisms, such as elasticity and elasto-
Parameter identi?cation of visco-elastic materials using back-analysis: Ohkami and Ihcikawa (1997), Ohkami and Swoboda (1999), Yang et al. (2001) [515–517]; Viscosity and yield strength degradation of rock: ! Nawrocki and Mroz (1999) [518]; Visco-elasto-plastic behaviour with ?nite strains: Nedijar (2002a, b) [519,520]; Comparison of formats and algorithms of viscoplasticity models: Runesson et al. (1999) [521]; Combination with ?uid ?ow and rock porosity: Abousleiman et al. (1993, 1996), [522–523]; Material softening simulations: Diez et al. (2000) [524]; Rheology of fractured rocks: Patton and Fletcher (1998) [525]; Poroelasticity of rocks with anisotropic damage: Shao et al. (1997) [526].
3.8.4. Size effects and homogenization Size effect is a special feature of fractured rocks, mainly due to two factors caused by the existence of fractures of various sizes in rock masses. The ?rst is the fact that the fracture systems divide the rock mass into a large number of sub-domains or blocks, whose sizes and interactions dominant the overall behaviour of rock masses. The second is the fact that the physical behaviour of fractures themselves is dependent on the sizes of fractures, due to scale-dependence of surface roughness of the rock fractures (Fardin et al., 2001) [527]. The state-of-the-art of this subject for rocks is presented in two edited volumes by Da Cunha (1990, 1993) [4,5] and a recent comprehensive survey of size effects in the strength and behaviour of structures, including geotechnical structures, was given by Bamant (2000) [528]. A shortcoming of classical plasticity theory is the lack of an intrinsic length scale in the models that is otherwise needed to explain the size effect in such theories. Since the 1990s, however, a theory of gradient plasticity was developed, which may be applied to consider plastic behaviour and Strain-localization of geological materials with regular fracture patterns, such as Aifantis (1992), Zbib and Aifantis (1989, 1992) [529–531]. Frantziskonis et al. (2001) [532] proposed a new scale-dependent constitutive model for heterogeneous materials like concretes, using wavelet analysis
322
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
techniques. Application of this theory has not been found in rock mechanics publications. The most representative and contemporary methods to take the size effect into account for the physical behaviour of fractured rocks is the equivalent continuum approach established on the basis of the Representative Elementary Volume (REV), through analytical or numerical processes of homogenization and/or upscaling and based on assumed constitutive models. The properties established are often called effective or equivalent properties. Except for a few cases with speci?c fracture sets (often assumed to be in?nitely large in size), closed-form solutions for equivalent properties do not exist for generally fractured rocks, and numerical simulations are often used as the tools of derivation. One exception is the crack tensor theory proposed by Oda (1986) [477], in which random fracture populations can be readily incorporated into the analytical functions de?ning the compliance and permeability tensors. For homogenization and upscaling, much work has been undertaken in general solid mechanics ?elds, especially for composite materials. The general methodology is well known, see for example the most recent publication by Fraldi and Guarracino (2001) [533]. For the mechanical properties of fractured rocks, a number of closed-form solutions were obtained with simpli?ed fracture system geometry, as mentioned in Section 3.8.1. Homogenization and upscaling techniques have a long history of application in deriving equivalent hydraulic properties of fractured rocks, such as in Snow (1965) [534]. Two comprehensive reviews by Renard and ! ! de Marsily (1997) [535] and Wen and Gomez-Hernandez (1996) [536] summarized the state-of-the-art. The more recent works on the subject are reported in Lee et al. (1995), Li et al. (1995), Tran (1996), Hristopulos and Christakos (1997), Scheibe and Yabusaki (1998), Pozdniakov and Tsang (1999), Shahimi and Mehrabi (1999), Lunati et al. (2001) and Zhang and Sanderson (1999) [537–544,418]. Derivation of coupled hydro-mechanical effective properties of general porous media was reported in Kachanov et al. (2001) [545]. Similar works were also reported for heat transfer processes in porous media, see Quintard et al. (1997) [546]. For ?ow in fractured porous media like rock masses, new techniques were also created to separate the contributions from the rock matrix and from the fracture systems (which are also treated as equivalent continua), the so-called dual (double) porosity models, dual (double) permeability models and dual (double) continuum models, in order to simplify the complexity in the fracture–matrix interaction behaviour and to partially consider the size effects caused mainly by the fractures. Some recent works in this ?eld can be seen in Zimmerman et al. (1996), Bai (1997), Bai et al. (1999), Choi et al. (1997), Masters et al. (2000), McLaren et al.
(2000), Vogel et al. (2000), Zhang et al. (2000), and Landereau et al. (2001) [547–555]. In these approaches, although the matrix and fractures are numerically ‘separated’, the size effects in the fracture system still exist, and the interactions between the matrix continua and fracture ‘continua’ still need to be considered. It is more than likely that the two ‘continua’ will have different REVs. Proper numerical methods to treat the effect of such differences in size effects of co-existing continua remain to be developed. In the above approaches and models, the interactions between the ?nite fractures were not considered. It appears that an effective approach to derive equivalent elastic constitutive models of rocks with randomly distributed ?nite fractures is the numerical approaches of homogenization and upscaling, such as reported by Stietel et al. (1996) and Lee and Pande (1999) using FEM [556,557], and Mas-Ivas et al. (2001) and Min et al. (2001) using DEM [285,286]. The existence of the REV for fractured rocks is still in debate and non-REV approaches were also proposed in Pariseau (1999) [558]. The focus of debate is whether such REVs can exist physically considering the presence of a hierarchic structure of fracture sizes and widths (apertures) and their vastly different physical behaviour and properties. 3.8.5. Damage mechanics models Constitutive models of rocks have also been developed using continuum damage mechanics principles, proposed ?rst by Kachanov (1958) [559], based on scalar, vector or tensor representations of the void formation, micro-cracking or embedded fracture phenomena in rock under loading. This theory is very closely related to both continuum mechanics and fracture mechanics, and serves as a bridge connecting the two (see Oliver, 2000; Oliver et al., 2002) [560,561]. It has a certain parallelism in the formulation with the plasticity models, such as using damage evolution laws in place of ?ow rules. The restriction by the normality rule in plasticity theory is, however, absent in damage mechanics principles. The damage mechanics theory has also a certain advantage in simulating the Strainlocalization factors using continuum approaches and study of brittle-ductile deformation model transitions observed during testing rock samples. Comprehensive reviews on its development, characteristics, trends and weaknesses are given in Krajcinovic (2000) and de Borst (2002) [562,563]. The damage mechanics approach has been applied to study strength degradation and Strain-localization phenomena in rocks and to formulate damage related constitutive models of rock and rock like materials, such as Kawamoto et al. (1988) and Ichikawa et al. (1990) [564,565]. A large number of papers has been published on these subjects and cannot be summarized in even
L. Jing / International Journal of Rock Mechanics & Mining Sciences 40 (2003) 283–353
323
moderate detail. Below are some recent publications concerning rock damage:
*
*
*
*
*
*
*
Constitutive models and properties of rocks: Aubertin and Simon (1997), Grabinsky and Kamaleddine (1997), Shao (1998), Shao and Rudnicki (2000), Yang and Daemen (1997), Chazallon and Hicher (1998), Homand-Etienne et al. (1998), Basista and Gross (1998), Zhao (1998), Carmelier (1999), Chen (1999), Dragon et al. (2000), Jessell et al. (2001), Li et al. (2001), Brencich and Gambarotta (2001) [566–580]; Damage induced permeability change in rock: Souley et al. (2001) [581]; Strain-localization and failure predictions: Comi (1999, 2001), Peerlings et al. (2002) [582–584]; Scale effects of damage in cemented granular rocks: Pisarenko and Gland (2001) [585]; Effects of damage on explosions, blasting and fragmentation of rocks: Yang et al. (1996), Liu and Katsabanis (1997a, b), Rossmanith and Uenishi (1997), Ma et al. (1998), Wu et al. (1999) [586–590]; Disturbed state concept-based model: Desai and Zhang (1998) [591]; Effects of damage on excavations and boreholes: Sellers and Scheele (1996), Cerrolaza and Garcia (1997), Elata (1997); Nazimko et al. (1997a, b), Swoboda et al. (1998), Zhu and Li (2000), Zhu et al. (1999) [592–598,273].
3.8.6. Rock fracture models Constitutive models for rock fractures play an important role in almost every aspect of rock mechanics and rock engineering, especially in the ?elds of numerical modelling and characterization. It is therefore not surprising that the mechanics and models of rock fractures has become one of the main themes in almost every international or national conference on rock mechanics, rock physics and rock engineering. The most directly relevant volumes are those edited by Stephansson (1985), Barton and Stephansson (1990), Myer et al. (1995) and Rossmanith (1990, 1995, 1998) [599–604]. The subject has also become an inevitable part of many text and reference books and parts of many edited volumes, such as Chernyshev and Dearman (1991), Lee and Farmer (1993), Selvadurai and Boulon (1995), Hudson (1993), Indrarantna and Hague (2000) and Indrarantna and Ranjith (2001) [605–610]. Some early comprehensive reviews on the experimental aspects, formulations of constitutive models and strength envelopes of rock fractures are given in Jing and Stephansson (1995), Ohnishi et al. (1996) and Maksimovic (1996), respectively [611–613]. In this section, we will only present a limited number of fundamental developments of importance. The constitutive models of rock fractures are formulated with mainly two approaches: empirical and
theoretical. The primary variables are contact tractions and relative displacements (instead of stresses and strains in continuum models), and aperture and ?ow rates (?uxes). Implementation of constitutive models into continuum numerical methods such as FEM often leads to so-called interface elements, which may also cause numerical instabilities when zerothickness of interface elements is employed, such as discussed in Kaliakin and Li (1995), Day and Potts (1994), and Lee and Pande (1999) [614–615,557]. Implementation into DEMs is generally a more straightforward use of contact mechanics principles, but the prevention of inter-penetration of solid blocks must be applied, using methods like penalty function, Lagrangian multipliers or augmented Lagrangian multiplier techniques. Besides the classical Coulomb or Mohr–Coulomb friction laws, the most well-known constitutive laws developed for rock fractures with an empirical approach are G