Numerical modelling
Contents

Introduction
Rock mass is a natural geological material, therefore the physical and engineering properties of it have to be established rather than defined through a manufacturing process like steel. Rock mass is a massively Discontinuous, Inhomogeneous, Anisotropic and NonElastic (DIANE) material^{[1]}. Due to dynamic movements of the upper crust of the earth, the rock mass is under continuous dynamic stresses; this happens through diversity of earth related movements like tectonic movements, earthquakes, glaciation cycles, subsidence, etc. The rock mass is a fractured porous material which contains fluids like water, oil, and gas which they can significantly affect the in situ stress fields under different temperature and pressure conditions.
Rock mechanics projects are getting more complex, larger, and also more demanding for modelling requirements. Therefore, designing a rock engineering model requires comprehensive knowledge of the geomaterial, engineering properties, and parameters of the rock mass. However, the quality and quantity of the supporting data for rock engineering projects can never be complete due to the inherent variability of the rock mass property. This is the reason that the rock mechanics modelling for rock engineering projects are both science and art. In most cases the model does not need to be complete but to be adequate for the purpose of the study ^{[2]}.This can be achieved only on the basis of scientific knowledge and also empirical judgments supported by accumulated experiences through longterm practices.
Rock mechanics modelling
Rock engineering projects can be categorized into surface or underground rock engineering. These projects can be for civil engineering purposes or for mining, petroleum or environmental engineering. The ability of the rock mechanics model of predicting the behavior of the rock mass is achieved through the variety of modelling methods. As discussed before, the rock mass property is exclusively depended on the location and the site under study, so a computer model based on the specific site characteristics of the rock mass should be adopted in order to have a model with most similar behavior to the reality.
There is a wide spectrum of modelling approaches as rock engineering models are designed for different purposes. These modelling approaches can be presented in different ways. Hudson in 2001 developed eight approaches based on four methods and two levels, (Figure 1) ^{[3]}.
As illustrated in Figure 1 the four columns represent the four main modelling methods:
Method A: Engineering design based on the previous empirical design and experiences,
Method B: Engineering design based on simplified analytic solutions e.g. Kirsch equations for stress analysis around circle excavations,
Method C: Basic numerical modelling which captures most relevant rock mass mechanisms, and
Method D: Extensive numerical methods which involves all rock mass mechanisms.
Rock mass characterization can be obtained in the first stage getting the benefit of site investigation. One important point in rock engineering is the lack of data and that is the reason of applying empirical approaches i.e. classification systems ^{[3]}.
Numerical methods in rock mechanics
In order to get the benefit of numerical modelling for rock engineering purposes, in some cases the model can be represented using finite number of welldefined components. The behavior of each component then can be calculated using mathematical solutions. Therefore the global behavior of the rock mass can be defined using the behavior of each component as well as the interrelation between individual components (elements) ^{[2]}. These models are also called discrete models and the solution is typically straightforward.
In some other problems, the problem domain is subdivided into infinite components and the problem can be handled using the mathematical assumption of an infinitesimal element. This can be dealt with using differential equations to describe the behavior of the system at the field under study. Such systems are called continuous and have infinite degrees of freedom ^{[2]}.
The rock mass which mostly consists of intact rock, joints, faults and dykes is recognized as a discrete system. There is not proper closedform solution for this kind of environment and numerical methods must be applied for solving practical rock engineering problems. Due to the differences in the material assumptions, different numerical methods based on continuous and discrete systems have been developed and introduced. The most commonly used numerical methods for rock engineering practices are:
1The Finite Difference Method (FDM)
2The Finite Element Method (FEM)
3The Boundary Element Method (BEM)
4Discrete Element Method (DEM)
From which the FDM, FEM and BEM are considered as Continuum methods while the DEM is a Discontinuum method.
Characterization of rock masses for numerical methods
For each of the above mentioned numerical methods, the rock mass properties should be specified as well as boundary and initial conditions. To achieve high accuracy, the main components of the rock i.e. the fractures, inhomogeneity, anisotropy and inelasticity should be indicated which leads to a more extensive model. One important parameter in rock engineering modelling is the scale effect. It is an additional problem especially in a fractured rock mass ^{[4]}.
Some of the other problems in regard to rock characterization are as follow:
1Difficulty of characterization of the in situ rock stress in the modelling region,
2Deviation of rock properties in large scale in the field from lab test results, (scale effect),
3Estimation of rock properties using empirical characterization techniques.
These problems indicate that the whole issue of rock characterization should be realized and adjusted for using different methods of numerical modelling. The use of different numerical models requires different types of rock mass characterization. Therefore the success of a numerical model is in capturing the rock reality characterization related to the specific method that the model has been constructed under its governing rules ^{[2]}. Figure 2 illustrates the different interpretation and discretization of numerical modelling methods for a specific fractured rock mass.
Finite Difference Method
Among the different numerical methods, the FDM is the oldest numerical method to obtain approximate solutions to Partial Derivatives Equations (PDEs) in engineering. It is a mathematical expression of the form f(x+xa)f(x+xb). If a finite difference is divided by xb xa, one gets a difference quotient. In simple words, finite difference method relies on discretizing a function on a grid (Figure 3). The approximation of derivatives by finite differences plays a central role in Finite Difference Methods for numerical solutions, especially boundary value problems ^{[2]}.
FDM is generally based on a Taylor series representation of a function. The accuracy in this concept depends on the order of approximation. In other word, the accuracy depends on the number of terms used in the Taylor series representation of the function. The method is widely used in rock engineering practices especially in seismic modelling in which a discrete version of the wave equation can be derived where the wave field is propagated starting from the source location (initial conditions) ^{[5]}. Accuracy of the derivative calculation for a given order depends on grid spacing. Having a small grid size leads to a high accuracy but it causes longer time and memory for computation.
FDM generally deals with replacing of the partial derivatives of the objective function (for example displacement) by differences defined over certain intervals in the coordinate directions, Δx, Δy, and Δz ^{[5]}. ^{[2]}. Therefore a system of algebraic simultaneous equations of the predefined objective function can be introduced to each node in the grid system of the domain of the interest ^{[6]}. (Figure 4).
For an elastic solid in 2D, the FDM equation of equilibrium at a point like (i,j) is given as:
where ak and bk (elastic properties of the solid) are functions of the grid intervals Δx and Δy and F_x^(i,j) and F_y^(i,j) are body forces in the point (i,j).
Boundary conditions also should be defined at the boundary nodes of the domain. Taking into account the boundary conditions, solution of the simultaneous algebraic system equation will assure the production of the required values for the objective function at all nodes.
FLAC
Introduced by ITASCA consulting group, FLAC is a finite difference numerical modelling code for advanced geotechnical analysis of soil, rock or other materials that may undergo plastic flow when they are yielded in 2D (Figure 5). In FLAC materials are represented by elements or zones which form a grid that can be adjusted by the user to fit the shape of the modelling object ^{[7]}.
Simulation of faults, joints or boundaries is possible with FLAC due to large strain simulation of continua. The software provides an interactive modelling environment for the user.
FLAC3D
FLAC3D is a numerical modeling code for advanced geotechnical analysis in 3D. It is used in analysis, testing and design by geotechnical, civil and mining engineers. The program is designed to enhance any kind of geotechnical engineering practice where continuum analysis is required (Figure 6).
Finite Element Method
FEM is one of the popular numerical methods in engineering including rock mechanics. This popularity in the rock engineering comes from FEM’s ability in dealing with complex material like rock mass. The method is applied for solving problems that are described by partial differential equations or that can be formulated as functional minimization ^{[8]}. Like the FDM, this is a method to represent the continuum. The domain of the study consists of finite elements, (e.g. triangle elements with three nodes in 2D and brick elements with eight nodes in 3D) (Figure 7).
There are three main steps in the FEM analysis:
1discretization of the domain,
2local approximation,
3solution of the global matrix equation.
Using the networks of elements, displacement in any point within the element (u_i^e) can be obtained using:
Where N_ij are often called the shape functions (or a matrix of interpolation) and u_i^j is the displacement at the corresponding nodes. Using the interpolation matrix (shape function) the original partial derivatives functions of the problem can be replaced by:
Where vector {u_i^e } is the nodal value vector of the displacement and vector 〖{f〗_i^e} is the initial state of stress and the body forces. In rock engineering related problems, 〖[K〗_ij^e] is the stiffness matrix given by:
Where matrix [D_i ] is the elasticity matrix and the matrix [B_i ] is the geometry matrix which is determined by the relation between strain and the displacement ^{[2]}.
Because of many special challenges that rock engineering projects encounter  fractures and discontinuities, heterogeneity, anisotropy, nonliner material, scale, and time effects  FEM has been applied to rock mechanics due to its ability in dealing with such problems ^{[10]}.
Phase2
Phase2 is a 2 dimensional FEM code for geoengineering applications, it can be used for excavation design, slope stability analysis, probabilistic analysis, consolidation, ground water seepage and dynamic analysis capabilities ^{[11]}. It also offers a wide range of support system modeling options. Material models for rock and soil include MohrCoulomb, HoekBrown and CamClay. Analysis features for modeling jointed rock allows the user to automatically generate discrete joint or fracture networks according to a variety of statistical models (Figure 8).
ABAQUS
Abaqus FEA (formerly ABAQUS) provides a complete and flexible solution for understanding the behavior of materials using the Finite Element Method. It also takes the advantage of the latest highperformance parallel computing environment which allows the user to include details in the designing model. It also allows the user to minimize assumptions while reducing turnaround time for desired results ^{[12]} (Figure 9).
Comparison of FDM and FEM for rock mechanics applications
In rock mechanics, there are different commercial codes available for FDM and FEM, so it is important to find the code that is more relevant and appropriate to the domain of the application. Also in most cases the choice depends on the availability of the code, its cost and the user’s personal preference. However each method has its own advantages and disadvantages. Table 1 shows different features of each modeling method.
Boundary Element Method
The boundary integration equation method has been introduced to stress analysis in solids in the late 60s by Rizzo and Cruse ^{[14]}. Development of BEM application in rock engineering practices have been first done in 1973 by Crouch and Fairhurst. The method is widely used for general stress and deformation analysis for underground excavations, soil structure interaction, ground water flow and fracturing processes.
The basic idea of the boundary element method is the approximation of the solution to a PDE by solving the PDE on the boundary and then use it to find the solution inside the domain ^{[15]}. Although it sounds like a strange idea, but BEM is recognized as a powerful tool for finding solution. The main difference between boundary element methods with finite element or finite difference methods is on discretizing the model. Unlike FDM and FEM, in the boundary element methods only the boundary needs to be discretized. In this case the problem will much easier especially in 2D where the boundary is just a curve.
In practical engineering problems, the stress and displacement at a point inside a domain is of the interest. But as mentioned before, unlike the FEM, in BEM the displacement and stress values at the given point, P, is evaluated by:
Where the parameters, Tij (tarction vector), Uij (displacement vector), and Bi (body force vector) should be revaluated to the new position of the source point of P in the domain.
Examine2D
Examine2D is a plane strain boundary element program for calculation of stresses and displacements around underground and surface excavation in rock ^{[16]}. The program is interactive, easy to use, ideal for quick parametric analysis, and preliminary design. Examine2D provides an integrated graphical environment for better visualization for the users. It is a CAD based program which allows for point and click geometry input and edit (Figure 10).
Map3D
Map3D is an integrated three dimensional CAD based program using boundary element numerical method for stability analysis ^{[17]}. It is suitable for modelling of rock and soil engineering problems. Models that can be designed in Map3D include underground excavations, rock slopes, open pits, tunnels, fractures and surface infrastructure loads. The software can simulate yielding zones of different moduli (e.g. stiff dykes or soft ore zones) (Figure 11).
Discrete Element Method
Discrete element method (DEM) or distinct element method is a numerical method which is used for computation of the motion and effect of a large number of small particles ^{[18]} . Unlike the other numerical methods, DEM is introduced as a mesh free approach to geoengineering projects including rock mechanics and mining sciences. The grid based continuum methods can provide a quite acceptable approximation of a fractured rock mass when the scale of the model is sufficiently large enough, however there are a large group of problems in which the behavior of individual rock joints is of the interest ^{[19]}. Using conventional mesh based continuum methods it is possible to take into account the discontinuities or slide lines, however, they still have limitations when new contacts form between blocks.
The distinct element method is an alternative approach which can directly approximate the structure of the blocks of a jointed rock mass using arbitrary polyhedral. Therefore the model in DEM has the ability of detecting new contacts between block surfaces resulted from block motions.
The followings are the formulation steps of DEM in rock mechanics applications:
1.Identification of block system topology
2.Importing of block deformation properties using FEM
3.Introducing an algorithm for contact detection
4.Obtaining constitutive models for the rock blocks as well as the fractures
5.Integration of the equation of the motion (Figure 12).
UDEC
UDEC (Universal Distinct Element Code) is a distinct element tool for modeling of jointed and blocky material. It is a twodimensional numerical code which simulates the quasistatic or dynamic response for the loading of a jointed structure represented as blocks ^{[20]} . UDEC utilizes a discontinuous approach by assembling of discrete blocks. Getting the benefit of the intersections between blocks which are considered as boundary conditions for the individual blocks, large displacements along the discontinuities and even block rotation can be simulated. Because of its explicit solution approach, UDEC can model nonlinear behavior and complex models.
The properties of discontinuities and blocks are assigned separately in UDEC. The blocks can also be assigned as rigid or deformable; in case of the latter, the blocks are defined as finite difference zones which each zone is behaving based on predefined linear or nonlinear stressstrain law. The motion of blocks along the discontinuities are also controlled by either linear or nonlinear forcedisplacement relations in normal and shear directions (Figure 13).
3DEC
3dimensional Distinct Element Code (3DEC) is a numerical code for advanced 3D geotechnical analysis of rock, soil and structural support ^{[22]}. It simulates the response of a jointed rock mass to a static or dynamic loading. Similar to UDEC, the model blocks can be assigned as rigid or deformable and the discontinuities are considered as boundary conditions between blocks. Material models include elastic, anisotropic, MohrCoulomb, Drucker Prager, bilinear plasticity, strain softening, creep, and user defined (Figure 14).
Advantages and disadvantages of different numerical methods for rock mechanics
Every numerical method has its inherent constraints based on the specific mathematical approach to finding the solution. This leads the user to first analyse the problem and find out which approach is the best to model the reality of the problem. Table 2 shows the advantages and disadvantages of each numerical method for rock mechanics problems.
Reference
 ↑ Harrison, J.P. and Hudson, J.A.(2000) Engineering rock mechanics. Part 2.
 ↑ ^{2.00} ^{2.01} ^{2.02} ^{2.03} ^{2.04} ^{2.05} ^{2.06} ^{2.07} ^{2.08} ^{2.09} ^{2.10} Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering.
 ↑ ^{3.0} ^{3.1} ^{3.2} Hudson, J.A. (2001) "Rock engineering case histories: key factors, mechanisms and problems, Proceedings of the ISRM Regional Symposium EUROC2001."
 ↑ Dacunha, A.P. (1990) scale effects in rock masses.
 ↑ ^{5.0} ^{5.1} Yang, L. and Mirnal, K. (2009) Advanced finitedifference methods for seismic modeling.
 ↑ Wheel, M.A. (1996) A geometrically versatile finite volume formulation.
 ↑ ^{7.0} ^{7.1} ^{7.2} ITASCA. FLAC 5.0 User's guide, technical brochure.
 ↑ Nikishkof, G.P. (2004).Introduction to the finite element method.
 ↑ http://en.wikipedia.org/wiki/Finite_element_method.
 ↑ Pande, G.N. and Beer, G. and William, J.R. (1990) Numerical methods in rock mechanics. NewYork : Wiley.
 ↑ https://www.rocscience.com/products/3
 ↑ ^{12.0} ^{12.1} http://www.3ds.com/productsservices/simulia/portfolio/abaqus/abaqusportfolio/
 ↑ McKinnon, S. (2014) Stability analysis in mine design, Depatment of mining engineering, Queen's University.
 ↑ Cruse, T.A. (1973) Application of the boundary integral equation.
 ↑ Costabel, M. (1986) Principles of Boundary Element Methods.
 ↑ ^{16.0} ^{16.1} Rocscience (2014) Quick start tutorial of Examine2D.
 ↑ ^{17.0} ^{17.1} http://www.map3d.com/example_grid_one.php
 ↑ http://en.wikipedia.org/wiki/Discrete_element_method
 ↑ Morris, J. and Glenn, S. and Blair, F. (2001) The distinct element methodapplication to structures in jointed rock. U.S. Department of Energy.
 ↑ http://www.itascacg.com/software/udec
 ↑ UDEC 4.0 User's guide.
 ↑ ^{22.0} ^{22.1} http://www.itascacg.com/software/3dec
 ↑ Hoek, E. and Grabinsky, M.W. and Diederichs, M.S. (1991) Numerical Modelling for Underground Excavation Design.