OF THE NUMERICAL METHODS FOR THE MODELLING OF ROCK MECHANICS PROBLEMS

Subject review The numerical methods have their origin in the early 1960s and even at that time it was noted that numerical methods can be successfully applied in various engineering and scientific fields, including the rock mechanics. Moreover, the rapid development of computers was a necessary background for solving computationally more demanding problems and the development process of the methods in general. Thus, we have many different methods presently, which can be separated into two main branches: continuum and discontinuum-based numerical methods. Some problems require the strengths of both main approaches which brought the hybrid continuum/discontinuum methods. The first goal of this paper is to present the state of the art numerical methods and approaches for solving the rock mechanics problems, as well as to give the brief explanation about the theoretical background of each method. The second goal is to emphasise the area of applicability of the methods in rock mechanics.


Introduction
Rock material is a natural geological material consisting of different minerals.It is discontinuous, anisotropic, inhomogeneous, inelastic and contains numerous randomly oriented zones of initiation of potential failure.These zones can be initial cracks, defects, cavities or other natural flaws.Rocks have been used (by humans) from the early beginnings of the human race in many different ways especially for the building purposes.We were constantly facing its unpredictable nature trying to build roads, tunnels, underground excavations, mining shafts or avoid the sudden rock falls and slopes.From the early days there was a need for predicting the behaviour of this material and up to present times, we are still trying to fully understand it.
Because of the irregular nature of rock material, the prediction of its behaviour has always been a challenging task.Rock mechanics developed mostly for the design of rock engineering structures and today we have a wide spectrum of different modelling approaches for many various rock mechanics problems.That includes the approaches based on the previous experiences, simplified mathematical models that can be solved analytically (e.g.Bishop Method of slices for slope stability) or general classification systems for rocks.The latter found large number of applications in various types of engineering projects like designing and construction of excavations in rock.There are several rock mass classification systems including RMR, Q and GSI classifications which are used to estimate the rock mass properties.The lack of information is a common fact in rock mechanics and engineering design so empirical approaches (including the classification systems) are also widely used.The development of computers made a significant contribution to the field called computational mechanics which found its wide use in rock mechanics.Namely, the numerical methods for approximately solving the partial differential equations are used today as the main approach in engineering design and research of rock mechanics.In this paper we are presenting the state of the art of numerical methods in terms of historical background as well as their strengths and weaknesses with the focus on deterministic approach.
Numerical model for the analysis of rock mass problems has to include the possibility of description of the rock behaviour by relevant material model and fracture mechanism, the presence of the pre-existing cracks, the pre-existing stress state, inhomogeneity and anisotropy as well as a time dependent behaviour caused by creep and plastic deformation.The additional challenge is coupling of the hydraulic and thermal process in analysing the mechanical behaviour of the rock.The choice of the appropriate model and numerical methods depends on the problem which has to be solved.
The numerical methods can be classified into three main categories: continuum, discontinuum and hybrid continuum/discontinuum methods.
Namely, the continuum concept implies that the domain of interest cannot be separated and the continuity between the points must be preserved in order to establish the derivatives.Thus the continuity between the elements in numerical approach must be preserved as well.Contrary to the analytical solution of differential equations where the solution is known in each point, the numerical solution is calculated in the finite number of pre-defined nodes which reduces the system complexity.The discontinuum approach on the other side treats the separate elements as discrete ones which are individually continuous and they interact between each other.In the discontinuous methods, the rigid body motion (usually with large movements) is the main case of interest, while in continuum-based methods the deformation of the system is the main objective.Thus, the sliding of the rock block on the pre-existing discontinuity would rather be computed with the discrete methods, while the choice of continuum method would be better used in calculating the deformation of the rock mass above the excavation of the tunnel.In addition there also exist the hybrid continuum/discontinuum methods that use the best properties of both approaches.Such an example is to have the discrete rock block which moves and also deforms under the external loading.The Finite Difference Method (FDM) is one of the oldest widely-applied numerical methods for solving the partial differential equations that found their application in the rock mechanics.The general principle of the method is replacing the partial derivatives of the function by the finite differences defined over a certain interval in the coordinate directions.More precisely, the domain needs to be partitioned into grid of nodes, among which the finite differences are defined.The different choices can be made for finite difference integration schemes: explicit or forward, implicit or backward and Crank-Nicolson or central difference scheme are largely used.It is worth noting the FDM can be used for solving the time dependent problems.
As a result of meshing the domain, a system of algebraic equations with unknowns related to the predefined nodes will arise.Each algebraic equation connected to its corresponding grid node is expressed as combination of function values at its own node, as well as at the surrounding nodes.After introducing the boundary conditions, the system of algebraic equations is finally solved (usually by direct or iterative methods) which produces the values of unknowns at each node leading to the approximate solution of the partial differential equation with an error made because of the difference in partial derivatives and finite differences.Such a direct kind of discretization, together with no practical need to use interpolation functions as in other methods like FEM or BEM, position this method as the most direct and intuitive technique for solving the partial differential equations.The advantage of the method is also in possibility of the simulation of complex nonlinear material behaviour without iterative solutions.
The standard FDM uses regular grid, such as rectangular which is also the most important shortcoming of the method.Thus, representing the irregular geometries together with dealing with heterogeneous nature of rocks and complex boundary conditions seems like a significant disadvantage of this method for simulating rock behaviour in practical rock mechanics tasks.Another important disadvantage lies in continuity requirement of the proposed equations which is not suitable for dealing with fractures.However, the general FDM method [1,2] evolved eventually to deal with these shortcomings, mostly with irregular grids, which include irregular quadrilateral, triangular and Voronoi grids.
Further development of this approach continued with the Finite Volume Method (FVM) [3].The FVM is also a method for solving partial differential equations, not directly as in FDM, but in integral sense.This leads to the formulation of finite volumes, which represent the volume surrounding each node in mesh.The basic principle is to replace the integrals with algebraic functions of the unknowns in the nodes, taking into account the initial and boundary conditions which lead to the set of algebraic equations.
The main advantages of the FVM are the possibilities of using the irregular unstructured meshes such as triangles, arbitrary quadrilaterals or Voronoi cells and considering the material heterogeneities [4].Namely, it is possible to join different material properties to different finite volumes.
The continuity requirement between the neighbouring nodes still makes the fracture propagation impossible to include, which is a main disadvantage of FDM/FVM.The analysis of the fracturing processes in FDM/FVM models can be conducted through the material failure at the nodes or cell centres, but in this way it is not possible to simulate the true fracture propagation.Despite this lacking, the FVM is still one of the most popular numerical methods in rock mechanics with applications in slope stability, tectonic process, rock mass characterization and especially in coupled hydromechanical problems.The overview of the FVM applications in rock mechanics can be seen in [5].The latest improvements of the FDM/FVM can be seen in [6][7][8][9].One of the most well-known computer programs for solving a variety of rock mechanics problems by FDM/FVM is FLAC/FLAC3D [10].

The finite element method (FEM)
The Finite Element Method (FEM) originates from the early 1960s [11,12].It was developed as an alternative to the finite difference method for the numerical solution of stress concentration in continuum mechanics and was the first numerical method which was able to account for material heterogeneities, nonlinearities, complex geometries and boundary conditions.Due to this, FEM immediately became the most applied numerical method in rock mechanics, especially because the FDM at that time was limited only to regular grids.The rapid application of the method started in the late 1970s and early 1980s with the early works of Zienkiewicz and Bathe [13,14].Following these works, many rock mechanics problems at that time were solved with FEM [15÷17].The FEM was developing further during the years, and today it is still the most applied numerical method for many advanced rock and soil mechanics simulations of the non-linear, time-dependent and anisotropic behaviour.
The FEM is the numerical method for finding approximate solutions to the boundary value problems for differential equations.The general principle is to divide the domain of the problem into smaller sub-domains called finite elements, do the local approximation inside each finite element, perform the finite element assembly and find the solution of the global matrix equation.More precisely, the unknown function (usually displacement function) needs to be approximated with a trial function of the nodal values in a polynomial form, where the numerical integration is performed in each element in Gaussian quadrature points.After the finite element assembly is performed, the algebraic system of equations on a global level is obtained.The FEM derived as a special case of Galerkin method where trial functions are presented globally, contrary to local approximation in FEM.
One of the mostly used advantages of FEM is a possibility of representing heterogeneous rocks, where it is possible to assign different material properties to different finite elements.Presently, there are many various shapes of finite elements with different number of nodes for 1D, 2D and 3D cases.Special case of elements called the 'infinite elements' was developed to simulate the far-field domain in geotechnical applications [18÷20].
Because of the continuum assumptions, the usual FEM methods have restrictions in efficient application of the failure analysis, cracking and damage induced discontinuities or singularities [21,22].Since rock is a discontinuous material and FEM is a continuum method, there have been many attempts to improve it in order to simulate the fracture propagation and other discontinuous effects with it.From the early experiments on rock specimens, it was observed that the experimentally obtained stress-strain curves up to failure are not linear.One of the earliest models that could approximately simulate the stress-strain nonlinear curve due to crack opening was smeared-crack models [23].These models were perfectly brittle in the beginning, although the rock has some residual load-carrying capacity after reaching its strength resulting with softening behaviour.Another attempt was an indirect representation of the discontinuities with the 'joint' elements where their influence on physical behaviour is considered through constitutive laws of the discontinuities as equivalent continuum.However, no real detachment is possible with them.The 'joint' elements are also limited to small displacements, while large movements across the discontinuities (e.g.sliding on the discontinuities) are not possible due to the continuum assumptions.The first element of this kind was 'Goodman joint element' developed especially for rock applications [24,25].Later improvements of the joint elements can be found in literature [26÷29].
The ultimate load computation where structural elements are subjected to progressive failure which leads to the collapse of the structure has been the topic of interest of many authors recently since this is one of the crucial types of failure mechanisms.The key difficulty in failure analysis is correct and mesh-independent representation of the post-peak softening behaviour related to crack propagation.Also, when trying to simulate the fracture growth with FEM, it is essential to have small size of elements and perform the continuous re-meshing as the crack propagates.The 'enhanced' FEM methods have been evolving to overcome these difficulties, and a couple of new methods derived to help alleviate shortcomings of the standard FEM.On one side there is the finite element method with embedded discontinuities (ED-FEM), representing cracks truly in each element (e.g.see [30÷38]).On the other side there is extended finite element method (X-FEM) where cracks are represented globally (e.g.[39÷42]).The ED-FEM and X-FEM methods are equivalent in their capabilities to handle the most demanding kinematics incorporating both strong and weak discontinuities (e.g.see [43]).The strong discontinuities serve for truly simulating the crack propagation, while the weak discontinuities enable the heterogeneous representation of material within element.General approach is to enhance the standard kinematics of the finite elements with additional discontinuous functions to simulate discontinuous behaviour.Thus, the enhanced kinematics in the strain and displacement fields serves respectively for dealing with heterogeneities and localization phenomena.The ED-FEM relies upon the incompatible mode method and Hu-Washizu mixed variational formulation [44,45].The most recent contribution to ED-FEM approach in rock mechanics was made in the following works [46÷49] where the crack propagation due to modes I, II and III, as well as their combination, provides successful simulation of complex failure mechanisms that occur in rocks.This kind of representation is also suitable for connecting the scales where the natural crack growth from the micro to macro scale can be simulated.Namely, the cracks form as a result of accumulation of micro-cracks leading to forming larger macro-cracks.
Another example of enhanced finite elements is a Generalized Finite Element Method (G-FEM) [50,51].The G-FEM uses the local functions which are usually analytical solutions to specific problems.These additional functions are not necessarily polynomials, but are bonded together with the standard space through the partition of unity principle [52].The advantage of G-FEM is that complex geometries can be represented easily because the mesh is independent of the geometry of the domain of the problem.The fractures can be simulated by additional functions and their corresponding additional nodes surrounding the crack.The other strength of this approach is the possibility to represent voids, problems with microscales and boundary layers, which is important for rock mechanics engineering.
Apart from FEM being suitable for representing heterogeneous materials and using the unstructured and irregular meshes, it also proved to be the appropriate tool for representing various non-linear and inelastic types of behaviour.These observations carried the challenges up to present time and many different models have been made for representation of material hardening and softening.Among them damage [53,54] and plasticity rock models [55÷57] are the most used frameworks for simulation of the non-linear and inelastic behaviour of material.The FEM is also suitable for representing the geometric non-linearities, contact mechanisms, fluidstructure interaction, connecting the scales from nano and micro scale to large macro scales, etc.All of these reasons position the FEM as the mostly used numerical method in rock mechanics.The generalized and enhanced FEM makes the method even more attractive because of its capabilities of fracture simulations without re-meshing.One of the most notable commercial FEM software for geotechnical applications including rock mechanics is PLAXIS [58].

The meshless methods
The application of the FEM in practical engineering problems with complex geometries, material properties and boundary conditions requires mesh generation which sometimes, especially in 3D problems, can be equally as demanding as solving the problem.Another weakness of the FEM is locking phenomenon which is often reflected in a numerical instability due to the mesh distortion.Both of these problems can be avoided by meshless methods which developed in a way that elements connecting the nodes are not required.Instead, the trial functions are no longer as in standard FEM, but generated from the neighbouring nodes within a domain of influence.More precisely, it is only necessary to generate the nodes across the domain without defining fixed element topology.The interpolation functions obtained in this way are no longer polynomial functions and it is more difficult to integrate them compared to standard FEM where Gauss integration points inside the elements are used for numerical integration.This makes the meshless methods computationally more demanding, but with an advantage of not using the standard meshes generators and easier representation of more complex geometries.Many meshless methods have been developed, but some of the mostly used ones are: Smoothed particle hydrodynamics [59], Diffuse element method [60], Element-free Galerkin method [61], Reproducing kernel particle methods [62], Moving least-squares reproducing kernel method [63], hp-cloud method [64], the method of finite spheres [65], Finite point method [66] etc.For the specific purpose in rock mechanics, where the treatment of the fractures is of essential influence, the meshless approach has significant potential which is shown by Zhang et al [67].

The boundary element method (BEM)
The BEM is a numerical computational method for solving partial differential equations where the solution of weak form is obtained globally through an integral statement.The basic principle of BEM is using the given boundary conditions to fit boundary values into the integral equation.Thus the discretization is needed only at the boundary with a finite number of boundary elements.After finding the solution on a boundary, the integral equation can be used again for calculation of solutions directly inside the domain.The main advantage of the BEM is reduction of model dimensions by 1, so for 2D BEM problems, the boundary elements are 1D lines that can be constant, linear or quadratic.In the 3D case, the boundary elements are 2D elements.The approximation of the solution at the boundary elements is performed using the shape functions similarly to FEM.The numerical integration is usually performed by using Gaussian quadrature points.When applying the boundary conditions into the matrix equations obtained after the approximation stage, the final global matrix equation with unknowns at boundary is obtained.Contrary to standard FEM matrix equations, the global matrix in BEM is usually asymmetric.The BEM application in rock mechanics has its origins in the early 1980s [68,69].This method found applications in underground excavations [70÷72], dynamic rock problems [73], analysis of in situ stress [74] and borehole drilling [75].Besides reducing the model dimensions by 1, its strength is accuracy in finding the solution because of its direct integral formulation.The standard BEM was developed for continuous and linear elastic solutions inside domain which was a disadvantage of the method in the beginning.The other main disadvantage over FEM is not efficient dealing with heterogeneity because not complete domain is discretized with elements.The fracture propagation with BEM was possible with the later development of two new approaches.The first one is using the domain division into sub-domains and the pre-defined crack path [76].The second one was Dual BEM (DBEM) [77] with application of the displacement and traction boundary conditions at opposite side of the fracture.With these improvements, the BEM became a notable tool for simulating the fracturing process.The other variants of BEM are Galerkin BEM (GBEM) [78,79] and Dual-reciprocity BEM (DRBEM) [80,81].The list of most notable continuum methods is given in Tab. 2.

The discrete element method (DEM)
The DEM started to develop in the field of the rock mechanic applications due to its requirements for modelling the discontinuous behaviour [82,83].The method was primarily defined as the computational approach that can simulate finite displacements and rotations of discrete bodies including their detachment.The theoretical formulations are based upon the equations of motion of rigid or deformable bodies using implicit or explicit time integrations.The basic concept is to treat the domain of interest as an assembly of particles or blocks which are continuously interacting between each other.In the DEM approach, the contact between components of the system is constantly changing during the deformation process.
Later on, many different numerical methods based on DEM developed for various rock mechanics problems.The main strength of the approach was the fact that the real discontinuities could be simulated, as well as representing the rock blocks which move and interact between each other including the fragmentation process etc.All of DEM based methods have the similar basic approach, while the differences between them are the use of various shapes of discrete elements, the way of computing the contact forces between the discrete bodies, the way of recognizing the contact, the way of integration of equations of motion, etc.The contact between the discrete bodies is essential part of solving the task.Nowadays, the mostly used contact algorithms are penalty, Lagrangian multiplier and augmented Lagrangian types of contact.The overview of the contact models can be found in [84].
According to the numerical integration scheme, there exist the two main approaches of DEM: explicit and implicit time integration DEM.The application of the explicit DEM started in the 1970s and was used mostly for rigid block simulations.The most representative explicit DEM method is Distinct Element Method developed by Cundall [85,86].It was created to simulate the fracturing, cracking and splitting of the blocks under external loading.The method was applied in many various rock mechanic applications including tunnelling, underground works and rock dynamics [87÷91], nuclear waste [92,93], rock slopes [94], acoustic emission in rock [95], boreholes [96], laboratory test simulations [97], etc.It is worth noting that irregular geometry of rocks is usually not easy to incorporate into numerical simulations.Thus, one of the latest and modern approaches to capture the realistic geometry and use it with DEM software is LIDAR laser scanning technology [98,99].The overview of explicit DEM approaches in rock engineering can be found in [100,101].The most notable computer codes for two and three-dimensional problems in rock mechanics using distinct element method are UDEC and 3DEC [102].
The most well-known method of implicit DEM approach is Discontinuous Deformation Analysis (DDA) [103,104].It is similar to the finite element method for solving deformation of the bodies, but it also accounts for interaction of the independent blocks along discontinuities in fractured and jointed rock mass.DDA uses an implicit algorithm for simultaneous solution of the equations of equilibrium by minimizing the total potential energy of the blocky rock mass system.The interaction between the independent blocks is performed by equations of contact interpenetration and accounts for friction.DDA uses the large displacement formulations for discontinuous movements between blocks.In this method, the standard mesh generation is used similarly to FEM, while penalty method or Lagrange multiplier method are used for the contact.The early formulations were limited to only simple representation of regular block movements and deformations, while in the later developments irregular blocks could be used and meshed with the triangular and four-node finite elements [105,106].It is an attractive method for solving rock mechanics problems where the rock blocks interact and deform between each other.The list of applications is covered in [107÷109].
Besides using the block systems in rock mechanics, the particle based DEM is an appropriate method for simulating the granular materials such as soil.The principle of solving tasks with granular materials is the same as for blocks, but with the simpler contact mechanisms due to rigid regular or irregular particles that do not deform.The regular particles include circular, elliptical and ellipsoidal shapes, while irregular ones can be polygonal or polyhedral.Some of the applications of particle based DEM in rock mechanics are rock fracturing and fragmentation due to rock blasting [110÷112], tunnelling [113], hydraulic fracturing [114] and ground movements [115].
Another approach similar to DEM is the Dynamic lattice network method used to simulate rock fracture initiation and propagation.The general idea is to have rigid particles covering the domain of interest, with springs acting as cohesive links and connecting the particles.The springs usually do not possess the mass and have the characteristic stiffness properties of material.The mass of the particles depend on the density of the material.During the simulation process, the springs deform until they are completely broken leaving the particles to move and interact with other particles.The application examples in rock engineering are mostly used in fracturing processes of intact rock [116÷118].One example of lattice models is found in [119], where particles are represented with Voronoi cells with geometrically nonlinear Reissner beams acting as cohesive links between them.Such a Voronoi cell representation allows for simulating the fracturing of the brittle heterogeneous material such as rock and concrete.
Lately, many numerical methods that combine advantages of finite and discrete elements developed.One of the well-known methods which combine finite and discrete elements is Discrete Finite Element Method in which the four-node blocks are used [120÷122].The main advantage of the method is a possibility to simulate interaction between the blocks, as well as the deformation of blocks.The higher order shape functions are used here for considering non-linear deforming of the blocks.The method uses implicit time integration and it is similar to DDA.Another similar approach is Combined Finite Discrete Element Method (FEM/DEM) [123,124] which considers block deformation as well as the fragmentation processes.Contrary to Discrete Finite Element Method, the FEM/DEM uses explicit time integration scheme.The FEM/DEM method is used in rock mechanics applications [125,126] and simulation of the behaviour of stone structures [127].There were some developments of the method considering the reinforcement between discrete elements with application in concrete, but which can be easily used to simulate the rock slopes or tunnel excavations with anchor reinforcements [128].The list of the covered discontinuum methods is given in Tab. 3.

DFN
The flow equations are solved with: ▪Closed form solutions ▪FEM ▪BEM ▪Pipe model ▪Lattice models

Discrete fracture network method (DFN)
The DFN is a discrete model which started to develop in the late 1970s mostly for simulating fluid flow and transport processes through the connected rock fractures.The groundwater transport is dominated by the connected paths formed by fractures, cracks, voids, etc.Thus the central problem is to define the position, geometry and properties of discrete fractures.Once the discrete fractures are defined deterministically or stochastically, the DFN model can be solved to understand the transport issues and flow nature.The flow equations within the fractures can be performed using closed-form solutions [129], FEM [130], BEM [131], pipe model [132] or the lattice model [133].For the smooth and regular fractures close-form solutions can be used, while for the geometrically irregular shapes, discretization needs to be performed and solution is found numerically.The pipe models and channel lattice models are solved analytically and are computationally less demanding than FEM and BEM models.Most notable software incorporating DFN is FRACMAN [134], with applications in mining, slope and tunnel fracture geology and other geomechanical applications.

Conclusion
Rock mechanics is among many other fields which benefit from the development of computers and numerical methods.Due to this fact, the rock mechanics evolved from traditional empirical approaches and estimations of their properties to modern science field which uses advanced numerical and mathematical approaches.The irregular rock properties such as heterogeneity, as well as the fracturing process, were some of the main difficulties that influence the choice of method.The FDM is widely used because of its simplicity and the possibility of handling the non-linear behaviour.At its beginnings, FDM was limited to regular mesh and thus was not capable of simulating irregular geometries and complex boundary conditions.Thanks to the more recent developments, the general FDM overcame these difficulties and is presently capable to model complex and irregular geometries.The FVM was capable to represent complex geometries and material heterogeneities, but due to representing the material failure at the nodes or cell centres, it is not possible to create fracture propagation.The FEM has been mostly used among all of the numerical methods since its beginnings, because of its capabilities to represent material heterogeneities, nonlinear behaviour such as damage and plasticity, complex geometries and boundary conditions.The original FEM could simulate the propagation of discontinuities using the re-meshing process which is computationally demanding especially for 3D cases.The constant development led to the enhanced FEM methods such as X-FEM, ED-FEM, G-FEM and meshless methods which enable to represent the true crack propagation and other localized and discontinuous effects which happen regularly in the rock mass.The main strength of BEM is to represent fracturing in rocks due to the most recent formulations, to reduce the problem complexity from 3D to 2D, or 2D to 1D and solve the problem at boundary.This is appropriate in solving the large scale problems, where number of degrees of freedom drastically decreases.The DEM approach solves the equations of motion and allows de-bonding and detaching of elements, thus representing real discontinuities.This is suitable for problems with large number of fractures which are dominant in failure process.This method was primarily developed for the rock mechanics problems and presently it is widely used for many rock applications.
It is shown in this paper that presently many different methods and approaches exist.There is no generalized method which solves all of the rock mechanics problems, but rather we need to choose the method which is most appropriate for the problem being solved.

Table 1
Overview of the most notable numerical approaches and methods ▪Finite Difference Method FDM ▪Finite Volume Method FVM ▪Finite Element Method FEM ▪Meshless Methods ▪Boundary Element Methods BEM Discontinuum Methods: ▪Discrete Element Method DEM ▪Discrete Fracture Network Method DFN Hybrid Methods: ▪Discrete Finite Element Method ▪Combined Finite Discrete Element Method FEM/DEM The most notable continuum numerical methods are

Table 2
The most representative and widely used methods in the continuum approach

Table 3
The most representative and widely used methods in the discontinuum approach