Numerical Methods in Rock Mechanics




1 Introduction



The complexity associated in the discipline of rock mechanics necessitates the use of modern numerical methods. With the rapid advancements in computer technology, numerical methods provide extremely powerful tools for analysis and design of engineering systems with complex factors that was not possible or very difficult with the use of the conventional methods, often based on closed form analytical solutions. It is only rarely that analytical solutions can be found to rock mechanics problems of practical concern. This may be due to the problem in describing the boundary conditions by simple mathematical functions, inhomogeneity of the problem domain, non-linear nature of the governing partial differential equations, or the constitutive relation for the rock mass are non-linear or otherwise insufficiently simple mathematically (E. T. Brown, 1987). In all these cases, approximate solutions may be found by using computer-based numerical methods.







There have been significant advances in the use of the computational methods in rock mechanics in the last three decades. Numerical methods and computing techniques have become tools for formulating conceptual models and mathematical theories integrating diverse information about geology, physics, construction techniques, economy, the environment and their interactions. This has greatly enhanced the development of modern rock mechanics from the traditional empirical art of rock deformability and strength estimation and support design to the rationalism of modern mechanics.







2 Computational Methods



Generally, numerical methods of solving boundary value problems can be divided into two categories – differential methods and integral methods. In differential methods, the problem domain is discretised into a set of sub-domains or elements. This method requires that physical or mathematical approximations be made throughout a bounded region. Solution procedure is based on numerical approximations of the governing equations, i.e. the differential equations of equilibrium, the strain displacement relations and stress-strain equations, as in classical finite difference methods. Alternatively, this procedure may exploit approximations to the connectivity of the elements and continuity of displacements and stresses between elements (as in finite element method).







In the integral methods, only the problem boundary is defined and discretised. Numerical solutions use analytical solution for simple singular problems in such a way as to satisfy approximately for each element the boundary conditions in terms of imposed tractions and displacements. Integral methods effectively provide a unit reduction in the dimensional order of a given problem since only the problem boundary is defined and discretised. This reduces the size of the system of equations to be solved and offers significant advantages in computational efficiency over differential methods. This makes integral methods useful for solving 3D problems in elastostatics which are sometimes of concern in engineering rock mechanics (Hocking et al., 1976; Watson & Cowling, 1985). Integral methods model far-field boundary conditions correctly, restrict discretisation errors to the problem boundary and ensure continuous variation of stress and displacement throughout the material. They are best suited for homogeneous materials and linear material behaviour.







3 Numerical methods of modelling rock masses



Rock mass is largely discontinuous, anisotropic and inhomogeneous in natural geological state. Difficulties arise in numerical modelling due to such complex and non-homogeneous geological conditions of rock mass. The complex combination of its constituents and its long history of formation make rock masses a difficult material for mathematical representation via numerical modelling.







Rock mass can be generally classified into three groups, i.e., (a) continuous, (b) discontinuous and (c) pseudo continuous groups. Type (a) refers to intact rock mass, type (b) represents jointed rock mass and type (c) is for highly fractured or weathered rock mass. The behaviour of type (a) rock mass can be analysed by means of model based on continuum mechanics, while a discontinuous model such as those proposed by Cundall (1971) may be used for analyzing the type (b) rock mass where joint elements in the finite element analysis are also useful. Discontinuous model similar to that of type (b) can be used for type (c) rock mass. However, it is almost impossible to explore all the joint systems and it seems that this type of rock mass behaves just like a continuous body in a global sense. Therefore, a continuum mechanics model can be used with the effect of discontinuities adequately considered in the model. This is achieved by homogenization technique where equivalent continuum properties of the rock mass is derived based on the geometry of the contained fracture systems and physical properties of the intact rock matrix and the fractures.







The concepts of continuum and discontinuum are not absolute but relative. The choice of continuum or discrete (discontinuum) methods depends on many problem-specific factors, and mainly on the problem scale and fracture system geometry. This is particularly true for rock mechanics problems. There are no absolute advantages of one method over another. Some of the disadvantages of each type can be avoided by using the continuum-discrete models, termed hybrid models.







3.1 Numerical techniques for rock mechanics



Rock mechanics problems reduce to the solution of boundary value problems having mixed boundary conditions. Singularities in the stress and strain field occur due to shape and mixed boundary conditions (Liebowitz et al., 1995). Material and geometrical non-linearities may further complicate the formulation and makes prediction of convergence very difficult. The advances in computer technology with rapid growth in algorithm methods make it possible in solving many complex problems in rock mechanics that do not have closed form analytical solutions. Commercial numerical analysis codes have become relatively user friendly in recent times.







A number of numerical techniques have been applied to problems in rock mechanics including finite element, finite difference and boundary integral methods. A brief outline on different numerical techniques in rock mechanics will be discussed below. The aim is to identify the essential principles of each method rather than to provide comprehensive account of the methods.







3.1.1 Finite element method (FEM)



Finite element is the most widely employed numerical method for rock mechanics and rock engineering. It does not require detailed programming experience to make efficient use of the finite element approach to problem solving in rock mechanics. However, familiarity with the fundamentals of the technique and with practical guidelines for generating reliable results is essential not only for the preparation of the program input, but also for recognition of faulty output.







Much of the FEM developments have been specifically oriented towards rock mechanics since early 1960s. This is mainly due to its flexibility for treatment of heterogeneity, non-linear deformability, in situ stresses, complex boundary conditions and gravity in early 1960s when traditional FDM with regular grids could not satisfy these essential requirements for rock mechanics problems. The well-known Goodman joint element in rock mechanics literature has been widely implemented in FE codes and improved by development of joint elements subsequently to address the limitation of the former and applied to many practical rock engineering problems.







The formulation of FEM is based on variational statement of the governing physics. FEM analysis constitutes three steps mainly, domain discretisation, local approximation and assemblage and solution of global matrix equation. This method involves the representation of continuum as an assembly of elements which are connected at discrete points called nodes. The problem domain is divided into discrete elements of various shapes, e.g. triangles and quadrilaterals in two-dimension cases and tetrahedrons and bricks in three dimensions. All forces are assumed to be transmitted through the body by the forces that are set up at the nodes. Expressions for these nodal forces, which are essentially equivalent to forces acting between elements, are required to be established. Continuum problem is analyzed in terms of sets of nodal forces and displacements for the problem domain.







The displacement components within the finite elements are expressed in terms of nodal displacements. Derivation of these displacements describes strain in the element. The stiffness of the medium to this induced strain determines stress in the element. Total stress within an element can be found out by superimposition of initial and induced stresses. The matrix of each element describes the response characteristics of the elements. These coefficient matrices are based on minimization of total potential energy. The elemental stiffness matrices are assembled to give the global stiffness matrix which is related to global force and displacement. As the number of elements in a problem domain tends to infinity, this is equivalent to solving differential equation.







FEM suffers limitations when simulating fracture problems mainly due to the limitation of small element size, continuous remeshing with fracture growth, conformable fracture path and element edges. However to overcome this limitation, discontinuous shape functions (Wan, 1990) are used for implicit simulation of fracture initiation and growth through bifurcation theory. Disadvantage of this method is that considerable time is required in preparing input data for a typical problem. This is particularly crucial in 3D problems and has led to the development of sophisticated mesh generation programs which eliminate much of the tedium involved in data preparation (Sheppard, 1988). FEM is computationally expensive. A large number of simultaneous equations must be solved to obtain a solution. If the problem is non-linear, the computation time increases enormously because the sets of simulataneous equations must be solved a number of times.







Inspite of all the above disadvantages, FEM is widely used due to its generality and flexibility to handle material heterogeneity, non-linearity and boundary conditions, with many well developed commercial codes with large capacity in terms of computing power, material complexity and user friendliness. It is one of several well-developed techniques that can provide useful information for engineering surface and underground excavations in rock. Static as well as dynamic analyses in two and three dimensions are possible. Most attractive advantage of FEM is the capability for direct inclusion of geological information in an analysis. Geometrical complexities, directional rock properties and various lithological units associated with surface topography, fault zones, igneous intrusions, existing excavations can be readily accommodated in FE approach. There are many finite element programs available with different degree of sophistication, ease of use and with considerable variation in cost.







3.1.2 Finite difference method (FDM)



This is one of the oldest numerical techniques used for the solution of sets of differential equations, given initial values and/or boundary values (e.g. Desai and Christian, 1977). The difference equations for a triangle are derived from the generalized form of Gauss’ divergence theorem (e.g. Malvern, 1969). Differential equations are solved by dividing the domain into connected series of discrete points called nodes. These nodes are the sampling points for the solution and are linked using finite difference operators to the governing equations. It is not necessary to combine the element matrices into a large global stiffness matrix as in the FE model. Instead, the FD method regenerates finite difference equations at each step. Derivatives of governing equations are replaced directly by algebraic expressions written in terms of field variables, e.g. stress or displacement, at discrete points in space (nodes).







The Finite Difference method allows one to follow a complicated loading path and highly non-linear behaviour without requiring the complex iterative procedure of a standard implicit code. Finite difference method can be used to discretize both time and space. It also provides easy error estimation techniques. It is particularly suitable for large, non-linear problem which may involve collapse or progressive failure.







Finite difference method is difficult to use for irregular shape domain or for problems involving singularities, because the fine meshing required near the singularity cannot be easily reduced for the rest of the domain. The conventional FD method with regular grid systems does suffer from shortcomings, most of all in its inflexibility in dealing with fractures, complex boundary conditions and material heterogeneity. This makes the standard FD method generally unsuitable for modelling practical rock mechanics problems. However, with the use of irregular meshes (triangular grid or Voronoi grid systems), which leads to Control Volume or Finite Volume techniques, significant progress has been made (Jing and Hudson, 2002). FLAC is the most well known computer code for stress analysis for engineering problems using FVM/FDM approach. Explicit representation of fractures is not easy in FD method as they require continuity of the functions between the neighbouring grid points. In addition, it is not possible to have special fracture elements as in FE method. However, this is the most popular numerical methods in rock engineering with applications covering from all aspects of rock mechanics, e.g., slope stability, underground openings, coupled hydro-mechanical etc.







3.1.3 Boundary element method (BEM)



Rock mass is predominantly very large and for practical purposes can be assumed to be of infinite extent. Because of its volume discretisation the Finite element is not very well suited for problems with a low ratio of boundary surface to volume since a large number of elements are required to model the response of the domain. Boundary element method is particularly attractive for such analyses in rock mechanics where the surface of the excavation has to be discretised. The amount of input data required to describe a problem is greatly reduced and the influence of infinite rock mass is automatically considered in the rock mass.







The Boundary Element method requires the discretisation of the domain and, if necessary the boundaries between the regions with different properties. For two-dimensional situations line elements at the boundary represent the problem, while for three-dimensional problems, surface elements are required. Thus, the dimensionality of the problem is reduced by one. This is particularly attractive as the amount of data required to describe the problem is greatly reduced as compared to FEM. Due to the BEMs advantage in reducing model dimensions, 3D application using the displacement discontinuity method for stress analysis has become efficient (Kuriyama and Mizuta, 1993; Kuriyama et al., 1995; Cayol and Cornet, 1997).







However, the computational algorithm is not so straightforward. In BEM, a system of simultaneous equations in terms of unknowns associated with nodes of the surface elements is solved. Boundary integral equation method solves linear boundary value problems with known green function solutions. Green’s function solution and the governing differential system are used to formulate boundary value problems as an equivalent surface integral. BEM can be direct or indirect depending upon the different mathematical approach, but for either of them to be practical it is necessary to be able to compute economically a function of two points in space known as the fundamental or basic singular solution. Both of these two methods in their simplest form fail if two surfaces in space such as crack are modelled. This is of consequences not only for modelling cracks but slots which are created by the mining of tabular ore bodies where the distances between the two parallel surfaces are taken negligible compared to other dimensions (Watson, 1993). Discontinuity displacement method (DDM) is developed for modelling such problems. DDM was initially developed as an indirect method of BEM, but it can be also derived as a direct method.







Although BEM method has been used in complex rock mechanics problems involving non-linear constitutive equations and number of materials, this method of analysis is particularly efficient in homogeneous, linear elastic problems in three dimensions. Advantage of this method is reduced in complex non-linear material laws with sets of materials because the surface needs be discretised wherever there is change of material properties, and hence the preparation of input data becomes more complicated. It is necessary to perform integrations over the volume of the rock mass for such an inhomogeneous rock mass or one with nonlinear behaviour and so some of the advantage over other methods are lost. Matrices of the equations in this method are not symmetric and banded as in FEM. Though the number of equations to be solved is less, the computation time is not reduced by the same amount.







Though BEM is most effective when the rock mass exhibits a linear elastic behaviour, there is no restriction on the complexity of boundary conditions which can be accommodated. Three basic approaches are used in the modelling of the discontinuity boundaries in BEM, mainly (i) interfacing two or more boundary element regions at the location of the discontinuities; (ii) using displacement discontinuities in combination with boundary elements; and (iii) coupling the BEM with FEM. BEM are readily applicable to analysis of stress in hard rocks at the present stage of development but usually not in coal measures and weak rock masses. BEM has enabled a significant practical problems of 3D analysis of stress in hard rock mass possible in simple desktop computers. Although BEM is never expected to replace the FEM and FDM, there is clearly a type of problems where they are the most effective solution.







Some notable applications of BEM in rock mechanics include:



· stress analysis of underground excavations with and without fractures,



· simulation of mining in faulted rock,



· dynamic problems,



· back analysis of in situ stress and elastic properties,



· borehole tests for permeability measurements.







3.1.4 Distinct element method (DEM)



DEM is one of the most rapidly developing areas of computational mechanics with a broad variety of applications in rock mechanics. DEM is relatively new and many think of it as “not yet proven” numerical technique for analysis and design in rock mechanics. Formulation and development of DEM have progressed over a long period of time since first studied by Cundall (1971). This method has been extended to other areas of research even though it was originally developed for rock mechanics applications. This method was initially developed for 2D representation of jointed rock mass, but now has been extended to applications in particle flow research (Walton, 1980), studies of granular material (Cundall and Strack, 1983), and crack development in rock and concrete (Plesha and Aifantis, 1983; Lorig and Cundall, 1987). UDEC and 3DEC are the most popular computer codes used to perform static as well as dynamic analysis. 3DEC has been primarily used to study rockbursting phenomena in deep underground mines.







This method treats domain as a discontinuum rather than a continuum in contrast to FEM and BEM. In DEM, rock mass is treated as an assemblage of rigid or deformable discrete blocks/particles. The contact displacements at the interfaces of a stressed assembly of blocks are identified and continuously updated throughout the deformation process and represented by appropriate constitutive models. The elements interact with one another through the forces developed at contact points. The equations of equilibrium are repeatedly solved until the laws of contacts and boundary conditions are satisfied.







Discrete element technique method is capable of analyzing multiple interacting deformable continuous, discontinuous or fracturing bodies undergoing large displacements and rotations. Dynamic equilibrium equation is solved for each body subjected to boundary interaction forces. There is no restriction on where one element may make contact with another, and nodes may interact with nodes or nodes with element faces. Forces generated between the contacting elements can be made to obey various interacting laws depending upon the physical nature of simulation. For example, interaction relationships for rock joints may include cohesion, dilation, damage to asperities, and stress dependent friction. Since elements rapidly change neighbouring elements upon fracturing or motion, automatic algorithm is used to compute connectivity or interaction of element to element. The governing dynamic equilibrium for each discrete element can be written in the general for as













where is the displacement, the superscript dots refer to differentiation with time, is the mass matrix, is the damping matrix, is the stiffness matrix, and is the applied load.







1.3.14.1 Discontinuous deformation analysis (DDA)



Implicit DEM is represented mainly by the Discontinuous deformation analysis (DDA) approach. DDA method has emerged as an attractive model for geomechanical problems because of its advantages over continuum based methods or the explicit DEM formulations. DDA has two advantages over the explicit DEM, mainly this method has relatively larger time steps and closed-form integrations for the stiffness matrices of elements. It is also capable to handle three-dimensional block system analysis and use of high order elements with more comprehensive representation of the fractures (Zhang and Lu, 1998). This code has been developed with application focusing on tunneling, caverns, fracturing and fragmentation process of geological and structural materials and earthquake methods (Lin et al., 1996; Yeung and Leong, 1997; Jing, 1998; Ohnishi and Cheng, 1999; Hsiung and Shi, 2001).







1.3.1.4.2 Key block approach



Key block approach similar to DEM, but without considerations of block deformation and motion was initiated by Warburton (Warburton, 1983; 1993) and Goodman and Shi (Goodman and Shi, 1985). This method identifies the key blocks (without stress or deformation analysis) formed by intersecting fractures and excavated free surfaces in the rock mass which have the potential for sliding and rotation in certain direction. This method is thus utilized for the analysis of the stability of rock structures which are characterized by rock blocks and fracture systems. This theory and its associated code developed has wide applications in rock engineering mainly in hard rocks in stability analysis, support design for slopes and underground excavations in fractured rocks.







DEM work for granular materials for geomechanics and civil engineering applications are widely reported in literature. The most well known code for this field are the PFC codes for both two dimensional and three dimensional problems, and the DMC code. Due to its conceptual attraction in the explicit representation of fractures, the DEM has wide application in rock engineering. Literature review shows a wide range of applicability of the methods as below:



Rock dynamics

Underground Works

Rock Slopes

Rock fracture

Laboratory test simulations and constitutive mode development

Hard rock reinforcement

Reservoir simulation

Fluid injection

Nuclear waste repository design and performance assessment

Stress-flow coupling

Intraplate earthquakes

Well and borehole stability

Rock permeability characterization

Acoustic emission in rock

Derivation of equivalent properties of fractured rocks





Although DEM is a general, flexible and powerful tool for analysing discontinuous rock mass, there are drawbacks to its use as standard analysis method. In addition to the difficulty associated with obtaining reliable data on location, orientation and persistence of the discontinuities, there is usually a lack of information on material behaviour at contacts as well as how to define damping of the system. DEM requires a considerable computation time to solve even simple problems. It remains as a qualitative tool and extremely useful in deformation and failure of blocky rock mass and provides insight into failure mechanisms.







3.1.5 Lattice models



Lattice models usually assume a linear-elastic material constitutive relation and, for this reason, they are conceptually very simple. This model has been applied in simulating fracture initiation and propagation in rocks, and to study the physics of rocks and nonlinear dynamics of earthquakes (Mora and Place, 1993; van Mier, 1995).







This model is similar to the discrete element method proposed by Cundall and Strack (Cundall, 1979). This model consists of interacting particles that are arranged in mesh of regular elements such as triangular elements and linked with massless springs whose stiffness and strength are based on the medium to be modelled (Mora and Place, 1993). The particle mass is derived from the density of the material which can be generated randomly to represent inhomogeneity of the medium. This method is based on molecular dynamics principles to model interacting particles by numerically solving their equations of motion. This technique is similar to that of the DEM for particle systems except it represents the continuum behaviour of the medium by assembly of particles and springs rather than as a direct discrete medium.







Lattice model formulated by Brandtzaeg (Brandtzaeg, 1927) and modified by Reinius (Reinius, 1956) and Baker (Baker, 1959) treats material as a network of brittle breaking bar or beam elements. Crack growth is simulated by removing elements that exceed their tensile strength. This model seems to predict some of the fracture behaviour in compression as well other tests (Burt and Dougill, 1977; Schlangen and van Mier, 1992; Schlangen, 1993). Material heterogeneity in the model is introduced by assigning different properties to the elements in the lattice, or by varying bar stiffness.







Beam lattice model for simulating fracture processes in concrete was used by van Mier (van Mier, 1995). The fracture law in this method becomes very simple when the microstructure of the material is included in the lattice model. A regular or random triangular framework of beam elements discretises the continuum. The elastic stiffness and Poisson’s ratio of the complete lattice is ensured to be equal to that of the continuum by adjusting the size of the beams. Heterogeneities in the material is introduced by overlaying mesh on top of a computer-generated particle distribution and different mechanical properties are assigned to the beams falling in each phase. Fracture process is reproduced by sequential removal of elements which are assumed to have linear elastic behaviour until failure. The element with maximum value of the ratio between the effective stress and the tensile strength is removed at each step. It is assumed that only mode I fracture occurs locally. The coefficient is used to adjust the contribution of the bending moment in a beam. Numerical uniaxial tensile experiments conducted showed that the model can reproduce the fracture processes observed in real physical experiments.







3.1.6 Discrete fracture network method (DFN)



DFN is a discrete model that has a wide applications for fluid flow and transport processes through a system of connected fractures where it is difficult to derive equivalent continuum flow and transport properties of fractured rock (Yu et al., 1999; Zimmerman and Bodvarsson, 1996). DFN uses FEM and BEM mesh, pipe models and channel lattice models. The pipe model and the channel lattice model provide simpler representations of the fracture system geometry. DFN approach is based on the stochastic simulation of fracture systems. The fractal concept has been applied to the DFN approach in order to consider the scale dependence of the fracture system geometry and for up-scaling the permeability properties (Barton and Larsen, 1985; Chiles, 1988; Barton, 1992).







FRANCMAN/MAFIC (Dershowitz et al., 1993) and NAPSAC (Stratford, 1990; Herbert, 1994; 1996) are some of the DFN formulations and computer codes with many applications for rock engineering. There are some limitations of this method. This method is highly dependent on the interpretation of the in situ fracture system geometry which can only be roughly estimated and cannot be validated practically. However, the DFN model has wide applications for problems of fractured rock. Below are some of the examples of developments and applications of the DFN approach.







Development of multiphase fluid flow

Hot-dry-rock reservoir simulations

Characterisation of permeability fractured rocks

Water effects on underground excavations and rock slopes.





3.1.7 Hybrid methods



Hybrid methods are used in rock engineering for stress/deformation and flow problems in fractured rocks. Hybrid methods combine FEM/BEM, DEM/FEM and DEM/BEM to take advantages of the strength of each method while avoiding many of its disadvantages. FEM/DEM is used for non-linear or fractured near fields where explicit representation of the fracture is needed. BEM is used for simulating far field rocks as equivalent elastic continuum. Hence, the hybrid of these methods provide numerical technique for effective representation of the effects of the far field to the near field rocks.







Zienkiewicz (1977) first proposed to couple FEM/BEM to avoid shortcomings of each of the models when used separately and optimise on computation time, efficiency and accuracy. This optimised approach of analysis has been found to give good results in geotechnical investigations conducted by Brebbia and Georgiou (1979), Beer and Meek (1981, 1983 and 1986), Brady and Wassyng (1981), Yeung and Brady (1982), Ohkami et al. (1985), Varadrajan et al. (1985 and 1987), Swoboda et al. (1987 and 1989) and Xiao and Carter (1993). In this method, region close to point of interest is discretised into finite elements while the far field is discretised using boundary element formulation. In rock mechanics, this method has been mainly used for simulating the mechanical behaviour of underground excavation (Varadarajan et al., 1985; Ohkami et al., 1985; Gioda and Carini, 1985; Swoboda et al., 1987; van Estrorff and Firuziaan, 2000).







In the hybrid BEM/DEM, the BEM region which surrounds DEM is represented by super block having contacts with smaller blocks along the interface with DEM. UDEC/3DEC are the computer codes that are based on this technique for stress/deformation analysis in rock mechanics. Literature review show a development of hybrid discrete-continuum models for coupled hydro-mechanical analysis of fractured rocks using combinations of DEM, DFN and BEM approaches (Wei, 1992; Wei and Hudson, 1998). In hybrid DEM/FEM model, DEM region consists of rigid blocks and the FEM region can have non-linear material behaviour. A hybrid beam-BEM model is based on the same principal as of hybrid BEM/FEM model and used to simulate the support behaviour of underground openings.







Hybrid methods have many advantages in the field of rock mechanics. However, care should be taken to ensure continuity or compatibility conditions in the interface between regions of different models. This is particularly important when different material assumptions are made.







3.1.8 Neural Networks



This is a new emerging concept where mechanisms occurring in rock mass in reality is not mapped directly but rather the rock mass is represented indirectly by a system of connected nodes. These nodes do not necessarily have any physical interpretation nor the input-output values. These network models provide descriptive and predictive capabilities and have been used extensively for rock parameter identification and engineering activities. This method tries to “mimick” the perception of human brain in the neural network so that the programmes can incorporate judgments based on empirical methods and experience. For this to happen, the model has to be trained with a large set of parameters.







There have been applications of neural networks to rock mechanics and rock engineering problems, some of them are outlined below:



· Rock mass classification



· Rock mass properties



· Rock fracture analysis



· Earthquake information analysis



· Rock slope displacement



· Stress-strain curves for intact rock



· Tunnel support



Limitation of this method is there is lack of theoretical basis for verification and validation of the techniques and their outcomes. The model may not reliable estimate outside its range of training parameters and this method has not yet provided an alternative to conventional modelling.







4 Concluding Remarks



There have been significant advances in computational methods over the last decade, specifically in numerical methods in solving rock mechanics problems. Formulation of conceptual models and mathematical theories integrating diverse information about geology, physics, construction techniques, economy, the environment and their interactions have been possible due to the development of numerical methods and computing techniques today. This has led to the development of modern rock mechanics from the traditional empirical art of rock deformability and strength estimation and support design to the rationalism of modern mechanics.







Due to the inherent nature of rock mass containing discontinuity, fractures and inhomogeneity, numerical modelling has become more challenging. Success of numerical modelling for rock mechanics can entirely depend upon the quality of the characterisation of the fracture system geometry, physical behaviour of the individual fractures and their interaction. Today’s numerical modelling capability can handle very large scale and complex equation systems, but still there are limitations in the quantitative representation of the physics of fractured rocks.







It is not possible to completely validate numerical models by experiments in rock mechanics due to the assumptions in mathematical models and complexities like fracture in rock mass. However, numerical models can be calibrated against laboratory and in situ experiments and the output of the results used to successfully analyse practical problems. This needs a combined scientific and engineering support is needed for applying numerical methods to rock mechanics and rock engineering.







5 References



Baker, A. L. L. (1959). “An Analysis of Deformation and Failure Characteristics of Concrete.” Mag. Concrete Res. 11: 119-128.



Barton, C. C. (1988). Fractal Analysis of the Scaling and Spatial Clustering of Fractures in Rock. Proc. 1988 GSA annual meeting on fractals and their applications to geology.



Barton, C. C. and E. Larsen (1985). Fractal Geometry of Two Dimensional Fracture Networks at Yucca Mountain, Southwestern Nevada. Proc. Int. Symp. on Rock Joints, Bjorkliden, Sweden.



Beer, G. (1983). “Finite Element, Boundary Element and Coupled Analysis with Applications in Geomechanics.” International Journal for Numerical Methods in Engineering 19: 567-580.



Beer, G. (1986). Implementation of Combined Boundary Element-Finite Element Analysis with Applications in Geomechanics. Devlopment. In Boundary Element Method. Banarjeeand Watson. London, Elsevier Applied Sciences. 4: 191-225.



Beer, G. and J. L. Meek (1981). Coupled Finite Element-Boundary Element Analysis of Infinite Domain Problems in Geomechanics. Numerical Methods for Coupled Problems. Bettessand Lewis. Swansea, U.K., Pineridge Press,: 605-629.



Brady, B. G. and A. Wassyng (1981). “A Coupled Finite Element-Boundary Element Method of Stress Analysis.” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 23(4): 475-485.



Brandtzaeg, A. (1927). “Failure of Material Composed of Non Isotropic Elements.” Det. Kgl. Norske. Videnskabers Selskabs 12: Grandheim.



Brebbia, C. A. and P. Georgiou (1979). “Combination of Boundary and Finite Elements in Elastostatics.” Appl. Math. Modelling 3: 212-220.



Brown, E. T. and J. Bray (1987). Analytical and Computational Methods in Engineering Rock Mechanics. London ; Boston, Allen & Unwin.



Burt, N. J. and J. W. Dougill (1977). “Progressive Failure in a Model Heterogeneous Medium.” J. Engng Mech. 103: 365-376.



Cayol, Y. and F. H. Cornet (1997). “3d Mixed Boundary Elements for Elastostatic Deformation Field Analysis.” Int. J. Rock Mech. Min. Sci. 34(2): 275-287.



Chiles, J. P. (1988). “Fractal and Geostatistical Method for Modelling a Fracture Network.” Math Geol. 20(6): 631-654.



Cundall, P. A. (1971). A Computer Model for Simulating Progressive Large Scale Movements in Blocy Rock Systems. Proc. Sympo. Int. Soc. Rock Mech.



Cundall, P. A. and O. D. L. Strack (1979). “A Discrete Numerical Model for Granular Assemblies.” Geotechnique 29: 47-65.



Cundall, P. A. and O. D. L. Strack (1983). Modelling of Microscopic Mechanisms in Granular Material. In Mechanics of Granular Materials: New Models and Constitutive Relations. J. T. Jenkinsand M. Satake. Amsterdam, Elsevier: 137-149.



Dershowitz, W. S., G. Lee, J. Geier, S. Hitchcock and P. la Pointe (1993). User Documentation: Francman Discrete Feature Data Analysis, Geometric Modelling and Exploration Simulations. Seattle, Golder Associates.



Desai, C. S. and J. T. Christian (1977). Numerical Methods in Geomechanics, New York: McGraw-Hill.



Gioda, G. and A. Carini (1985). “A Combined Boundary Element-Finite Element Analysis of Lined Openings.” Rock Mechanics and Rock Engineering 18: 293-302.



Goodman, R. E. and G. Shi (1985). Block Theory and Its Application to Rock Engineering., Prentice-Hall: Englewood Cliffs, NJ.



Herbert, A. W. (1994). Napsac (Release 3.0) Summary Document. AEA D&R 0271 Release 3.0, A. E. A. Technology, Harwell, U.K.



Herbert, A. W. (1996). Modelling Approaches for Discrete Fracture Network Flow Analysis. Coupled Thermo-Hydro-Mechanical Processes of Fractured Media-Mathematical and Experimental Studies. O. Stephansson, L. Jingand C. F. Tsang. Amsterdam, Elsevier: 213-229.



Hocking, G., E. T. Brown and J. O. Watson (1976). Three Dimensional Elastic Stress Analysis of Underground Openings by the Boundary Integral Equation Methods. Proc. 3rd Symp. Applns. Solid Mechs., Toronto, University of Toronto Press.



Hsiung, S. M. and G. Shi (2001). Simulation of Earthquake Effects on Underground Excavations Using Discontinuous Deformation Analysis (Dda). Rock Mechanics in the National Interest. Elworth, Tinucciand Heasley, Swets & Zeitlinger Lisse: 1413-1420.



Jing, L. (1998). “Formulation of Discontinuous Deformation Analysis (Dda) - an Implicit Discrete Element Model for Block Systems.” Eng. Geol. 49: 371-381.



Jing, L. (2003). “A Review of Techniques, Advances and Outstanding Issues in Numerical Modelling for Rock Mechanics and Rock Engineering.” Int. J. Rock Mech. Min. Sci. 40: 283-353.



Jing, L. and J. A. Hudson (2002). “Numerical Methods in Rock Mechanics.” Int. J. Rock Mech. Min. Sci. 39: 409-427.



Kuriyama, K. and Y. Mizuta (1993). “Three Dimensional Elastic Analysis by the Displacement Discontinuity Method with Boundary Division into Triangle Leaf Elements.” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 30(2): 111-123.



Kuriyama, K., Y. Mizuta, H. Mozumi and T. Watanabe (1995). “Three Dimensional Elastic Analysis by the Boundary Element Method with Analytical Integrations over Triangle Leaf Elements.” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 32(1): 77-83.



Liebowitz, H., J. S. Sandhu, J. D. Lee and F. C. M. Menandro (1995). “Computational Fracture Mechanics: Research and Application.” Engineering Fracture Mechanics 50(5/6): 653-670.



Lin, C. T., B. Amadei, J. Jung and J. Dwyer (1996). “Extensions of Discontinuous Deformation Analysis for Jointed Rock Masses.” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 33(7): 671-694.



Lorig, L. J. and P. A. Cundall (1987). Modelling of Reinforced Concrete Using the Distinct Element Method. In Fracture of Concrete and Rock. S. P. Shahand S. E. Swarty. SEM, Bethel, CT: 459-471.



Malvern, L. E. (1969). "Introduction". In Mechanics of a Continuous Medium, Englewood Cliffs, New Jersey: Prentice Hall.



Mora, P. and D. Place (1993). “A Lattice Solid Model for the Nonlinear Dynamics of Earthquakes.” International Journal of Modern Physics 4(6): 1059-1074.



Ohkami, T., Y. Mitsui and T. Kusama (1985). “Coupled Boundary Element/Finite Element Analysis in Geomechanics Including Body Forces.” Computers and Geotechnics 5: 263-278.



Ohnishi, Y. and G. Chen (1999). “Simulation of Rock Mass Failure with Discontinuous Deformation Analysis.” J. Soc. Mater. Sci. Japan 48(4): 329-333.



Plesha, M. E. and E. C. Aifantis (1983). On the Modelling of Rocks with Microstructure. In Rock Mechanics - Theory - Experiment - Practice. C. C. Mathewson. New York, Association of Engineering Geologists: 27-35.



Reinius, E. (1956). “A Theory of Deformation and Failure of Concrete.” Mag. Concrete Res. 8: 157-160.



Schlangen, E. (1993). “Experimental and Numerical Analysis of Fracture Process in Concrete.” 28: 1-111.



Schlangen, E. and J. G. M. van Mier (1992). “Simple Lattice Model for Numerical Simulation of Fracture of Concrete Materials and Structures.” Mater. Structures 25: 534-542.



Sheppard, M. S. (1988). “Approaches to the Automatic Generation and Control of Finite Element Meshes.” Appl. Mech. Rev. 41: 169-185.



Stratford, R. G., A. W. Herbert and C. P. Jackson (1990). A Parameter Study of the Influence of Aperture Variation on Fracture Flow and the Consequences in a Fracture Network. Rock Joints. N. Bartonand O. Stephansson, Rotterdam: Balkema: 413-422.



Swoboda, G., W. Mertz and G. Beer (1987). “Rheological Analysis of Tunnel Excavation by Means of Coupled Finite Element (Fem)-Boundary Element (Bem) Analysis.” Int. J. Numer. Anal. Meth. Geomech. 11: 115-129.



Swoboda, G., W. Mertz and A. Schmid (1989). Three Dimensional Numerical Models to Simulate Tunnel Excavation. Proc. 3rd Int. Symp. on Numerical Models in Geomechanics (NUMOG III), Elsevier Applied Science, London and New York.



van Mier, J. G. M., E. Schlangen and A. Vervuurt (1995). Lattice Type Fracture Models for Concrete. Continuum Models for Materials with Microstructure. H. B. Muhlhaus, John Wiley & Sons: 341-377.



Varadarjan, A., K. G. Sharma and R. B. Singh (1985). “Some Aspects of Coupled Fem-Bem Analysis of Underground Openings.” Int. J. Numer. Anal. Methods Geomech. 9: 557-571.



Varadarjan, A., K. G. Sharma and R. B. Singh (1987). “Elasto-Plastic Analysis of an Underground Opening by Fem and Coupled Fem-Bem.” Int. J. Numer. Anal. Methods Geomech. 11: 475-487.



von Estorff, O. and M. Firuziaan (2000). “Coupled Bem/Fem Approach for Non-Linear Soil/Structure Interaction.” Eng. Anal. Boundary Elements 24: 715-725.



Walton, O. R. (1980). Particle Dynamic Modelling of Geological Materials. Livermore, CA, Lawrence Livermore National Laboratory.



Wan, R. C. (1990). The Numerical Modelling of Shear Bands in Geological Materials. Edmonton, University of Alberta.



Warburton, P. M. (1983). Application of a New Computer Model for Reconstructing Blocky Rock Geometry, Analysing Single Rock Stability and Identifying Keystones. Proc. 5th Int. Cong. ISRM, Melbourne.



Warburton, P. M. (1993). Some Modern Developments in Block Theory for Rock Engineering. Comprehensive Rock Engineering : Principles, Practice, and Projects. J. A. Hudson. Oxford ; New York, Pergamon Press. 3: 293-315.



Watson, J. O. and R. Cowling (1985). Applications of Three-Dimensional Boundary Element Method to Modelling of Large Mining Excavations at Depth. Proc. 5th Int. Symp. Numerical methods in geomechanics., Rotterdam: Balkema.



Wei, L. (1992). Numerical Studies of the Hydromechanical Behaviour of Jointed Rocks., Imperial College of Science and Technology, University of London.



Wei, L. and J. A. Hudson (1988). “A Hybrid Discrete-Continuum Approach to Model Hydro-Mechanical Behaviour of Jointed Rocks.” Eng. Geol. 49: 317-325.



Xiao, B. and J. P. Carter (1993). “Boundary Element Analysis of Anisotropic Rock Masses.” Engg. Anal with Boundary elements 11: 293-303.



Yeung, D. and B. G. Brady (1982). A Hybrid Quadratic Isoparametric Finite-Boundary Element Code for Underground Excavation Analysis. Proc. 23rd U.S. Symp. on Rock Mechanics, University of California, Berleley.



Yeung, M. R. and L. L. Loeng (1997). “Effects of Joint Attributes on Tunnel Stability.” Int. J. Rock Mech. Min. Sci. 34(3/4): Paper No. 348.



Yu, Q., M. Tanaka and Y. Ohnishi (1999). An Inverse Method for the Model of Water Flow in Discrete Fracture Network. Proc. 34th Janan National Conference on Geotechnical Engineering, Tokyo.



Zhang, X. and M. W. Lu (1998). “Block-Interfaces Model for Non-Linear Numerical Simulations of Rock Structures.” Int. J. Rock Mech. Min. Sci. 35(7): 983-990.



Zienkiewicz, O. C., D. W. Kelly and P. Bettess (1977). “The Coupling of the Finite Element Method and Boundary Solution Procedures.” Int. J. Numer. Methods Eng. 11: 355-375.



Zimmerman, R. W. and G. S. Bodvarsson (1996). “Effective Transmissivity of Two-Dimensional Fracture Networks.” Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 33(4): 433-436.