View source for Numerical modelling

==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<ref name="Harrison"> Harrison, J.P. and Hudson, J.A.(2000) Engineering rock mechanics. Part 2. </ref>. 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 <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.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) <ref name="Hudson"> Hudson, J.A. (2001) "Rock engineering case histories: key factors, mechanisms and problems, Proceedings of the ISRM Regional Symposium EUROC2001." </ref>.



[[File:Figure 1.png|thumb|center|800px| Figure 1: Eight different approaches to rock engineering design and rock mechanics modelling <ref name="Hudson"> Hudson, J.A. (2001) "Rock engineering case histories: key factors, mechanisms and problems, Proceedings of the ISRM Regional Symposium EUROC2001". </ref>]]


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 <ref name="Hudson"> Hudson, J.A. (2001) "Rock engineering case histories: key factors, mechanisms and problems, Proceedings of the ISRM Regional Symposium EUROC2001" </ref>.

==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) <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>. 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 <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.

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 <ref name="Dacunha"> Dacunha, A.P. (1990) scale effects in rock masses. </ref>.

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),

'''3-'''Estimation 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 <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>. Figure 2 illustrates the different interpretation and discretization of numerical modelling methods for a specific fractured rock mass.


[[File:Interpretation_of_a_rock_mass.png|thumb|center|800px|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).<ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.]]

===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 <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.



[[File:Grid_concept_in_discretizing_of_a_function.png|thumb|center|300px|Figure 3-Grid concept in discretizing of a function <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.]]



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) <ref name="Yang"> Yang, L. and Mirnal, K. (2009) Advanced finite-difference methods for seismic modeling. </ref>. 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 <ref name="Yang"> Yang, L. and Mirnal, K. (2009) Advanced finite-difference methods for seismic modeling. </ref>. <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>. 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 <ref name="Wheel"> Wheel, M.A. (1996)  A geometrically versatile finite volume formulation. </ref>. (Figure 4).



[[File:Regular_quadrilateral_grid_for_FDM.png|thumb|center|700px|Figure 4-Regular quadrilateral grid for FDM <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.]]



For an elastic solid in 2-D, the FDM equation of equilibrium at a point like (i,j) is given as:

[[File:Formule_1-AH_(1).PNG|center|700px]]


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 <ref name="ITASCA"> ITASCA. FLAC 5.0 User's guide, technical brochure. </ref>. 

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.



[[File:Slope_stability_analysis_using_FLAC.png|thumb|center|700px|Figure 5-Slope stability analysis using FLAC-Plot of displacement vectors and x-velocity contours <ref name="ITASCA"> ITASCA. FLAC 5.0 User's guide, technical brochure. </ref>.]]

====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).



[[File:Simple_bench_in_FLAC_3D.png|thumb|center|1000px|Figure 6-A simple bench slope in FLAC3D <ref name="ITASCA"> ITASCA. FLAC 5.0 User's guide, technical brochure. </ref>. ]]

===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 <ref name="Nikishkof"> Nikishkof, G.P. (2004).Introduction to the finite element method. </ref>. 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).



[[File:Domain_discretization_FEM.png|thumb|center|600px|Figure 7-2-D Domain discretization using triangle elements <ref name="wiki"> http://en.wikipedia.org/wiki/Finite_element_method. </ref>. ]]


There are three main steps in the FEM analysis:

	''1-discretization of the domain,''

	''2-local approximation,''

	''3-solution of the global matrix equation.''



Using the networks of elements, displacement in any point within the element (u_i^e) can be obtained using:


[[File:Formule_2-AH_(2).PNG|center|200px]]



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:


[[File:Formule_3-AH_(1).PNG|center|350px]]


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:

[[File:Formule_4-AH.PNG|center|450px]]




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 <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.

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 <ref name="Pande"> Pande, G.N. and Beer, G. and William, J.R. (1990) Numerical methods in rock mechanics. NewYork : Wiley. </ref>. 

====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 <ref name="rocscience"> https://www.rocscience.com/products/3 </ref>. 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).



[[File:Phase_2_FEM.png|thumb|center|800px|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 <ref name="Abaqus"> http://www.3ds.com/products-services/simulia/portfolio/abaqus/abaqus-portfolio/ </ref> (Figure 9).



[[File:Abaqus.png|thumb|center|700px|Figure 9-Stress analysis around a hole in a plate <ref name="Abaqus"> http://www.3ds.com/products-services/simulia/portfolio/abaqus/abaqus-portfolio/ </ref>.]]

===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.


[[File:Table_1-FDM_vs_FEM.PNG|thumb|center|800px|Table 1-Comparison of FDM and FEM <ref name="steve"> McKinnon, S. (2014) Stability analysis in mine design, Depatment of mining engineering, Queen's University. </ref>]]

===Boundary Element Method===

The boundary integration equation method has been introduced to stress analysis in solids in the late 60s by Rizzo and Cruse <ref name="Cruse"> Cruse, T.A. (1973) Application of the boundary integral equation. </ref>. 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 <ref name="Costabel"> Costabel, M. (1986) Principles of Boundary Element Methods. </ref>.  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:


[[File:Formule_5-AH.png|center|600px]]


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 <ref name="Rocscience"> Rocscience (2014) Quick start tutorial of Examine2D. </ref>. 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).

[[File:Examine2D.png|thumb|center|800px|Figure 10-Displacements around a horseshoe shaped tunnel <ref name="Rocscience"> Rocscience (2014) Quick start tutorial of Examine2D. </ref>.]]

====Map3D====

Map3D is an integrated three dimensional CAD based program using boundary element numerical method for stability analysis <ref name="Map3D"> http://www.map3d.com/example_grid_one.php </ref>. 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).

[[File:Map_3D.png|thumb|center|800px|Figure 11-Modelling of slip on the fault intersecting a tabular mine <ref name="Map3D"> http://www.map3d.com/example_grid_one.php </ref>.]]

===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 <ref name="DEM"> http://en.wikipedia.org/wiki/Discrete_element_method </ref> . 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 <ref name="Morris"> Morris, J. and Glenn, S. and Blair, F. (2001) The distinct element method-application to structures in jointed rock. U.S. Department of Energy. </ref>. 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 [http://minewiki.engineering.queensu.ca/mediawiki/index.php/Numerical_modelling_-_constitutive_models constitutive models] for the rock blocks as well as the fractures

5.Integration of the equation of the motion (Figure 12).


[[File:DEM_flowshit.png|thumb|center|400px|Figure 12-DEM Computation Cycle <ref name="Jing"> Jing, L. (2003) A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. </ref>.]]

<br>

====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 <ref name="Itasca Udec"> http://www.itascacg.com/software/udec </ref> . 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).


[[File:UDEC.png|thumb|center|800px|Figure 13-Principal stress distribution around tunnels in UDEC <ref name="UDEC"> UDEC 4.0 User's guide. </ref>.]]

====3DEC====

3-dimensional Distinct Element Code (3DEC) is a numerical code for advanced 3D geotechnical analysis of rock, soil and structural support <ref name="3DEC"> http://www.itascacg.com/software/3dec </ref>. 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).

[[File:3DEC.png|thumb|center|800px|Figure 14- Stability analysis of an open pit using 3DEC <ref name="3DEC"> http://www.itascacg.com/software/3dec </ref>.]]

==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.


[[File:Table_2-Ah.PNG|thumb|center|800px|Table 2-Advantages and disadvantages of different numerical methods <ref name="Hoek"> Hoek, E. and Grabinsky, M.W. and Diederichs, M.S. (1991) Numerical Modelling for Underground Excavation Design. </ref>.
]]

== Reference ==
<references />