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 Non-Elastic (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 long-term 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].

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 well-defined 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 inter-relation 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 sub-divided 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 closed-form 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:

1-The Finite Difference Method (FDM)

2-The Finite Element Method (FEM)

3-The Boundary Element Method (BEM)

4-Discrete 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:

1-Difficulty of characterization of the in situ rock stress in the modelling region,

2-Deviation of rock properties in large scale in the field from lab test results, (scale effect),

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.

Figure 2- Interpretation of a rock mass shown in (a), by FDM or FEM shown in (b), BEM shown in (c), and DEM shown in (d).[2].

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 pre-defined 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 2-D, 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 2-D and brick elements with eight nodes in 3-D) (Figure 7).

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, non-liner 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 geo-engineering 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 Mohr-Coulomb, Hoek-Brown and Cam-Clay. 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).

Figure 8- Stress analysis around mine access shafts using Phase2

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 high-performance 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).

Figure 11-Modelling of slip on the fault intersecting a tabular mine [17].

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:

UDEC

UDEC (Universal Distinct Element Code) is a distinct element tool for modeling of jointed and blocky material. It is a two-dimensional numerical code which simulates the quasi-static 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 non-linear 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 non-linear stress-strain law. The motion of blocks along the discontinuities are also controlled by either linear or non-linear force-displacement relations in normal and shear directions (Figure 13).

3DEC

3-dimensional 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, Mohr-Coulomb, 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.

Table 2-Advantages and disadvantages of different numerical methods [23].