Get Complete Project Material File(s) Now! »

## Finite Difference and Finite Volume Method (FDM, FVM)

The Finite Difference Method (FDM) is one of the oldest widely-applied numerical methods for solving the partial differential equations that found its application in the rock me-chanics field. The general principle of the method is replacing the partial derivatives of the function by the finite differences defined over the 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. It’s worth noting that 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 pre-defined 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 grids, 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 [Perrone and Kao, 1975, Brighi et al., 1998] evolved eventually to deal with this shortcomings, mostly with irregular grids, which include irregular quadrilateral, triangular and Voronoi grids.

The further development of this approach continued with the Finite Volume Method (FVM) [Wheel, 1996]. The FVM is also a method for solving the partial differential equations, not directly as in FDM, but in integral sense. The method is based on 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 arbitrary triangles, quadrilaterals or Voronoi cells and considering the material heterogeneities [Fallah et al., 2000]. 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 hydro-mechanical problems. The latest improvements of the FDM/FVM can be seen in [Benito et al., 2001, Onate et al., 1994, Lahrmann, 1992, Jasak and Weller, 2000].

### Finite Element Method (FEM)

The Finite Element Method (FEM) originates from the early 1960s and the works of [Clough, 1960, Argyris, 1960]. 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, non linearities, 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 when many rock mechanics problems at that time were solved with FEM [Naylor et al., 1981]. The FEM was evolving fast 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 general principle of FEM 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 functions 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. One of the mostly used advantages of FEM is a possibility of representing heterogeneous rocks, where it is possible to assign the 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 [Zienkiewicz et al., 1983].

Because of the continuum assumptions, the FEM method has restrictions in efficient application of the failure analysis, cracking and damage induced discontinuities or singularities [Ibrahimbegovic, 2009]. 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 model.

This model was 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 [Goodman et al., 1968].

#### Boundary Element Method (BEM)

The basic principle of BEM is 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 interpolation can be used for calculating the solutions inside the domain. The main advantage of the BEM is reduction of model dimensions by 1. The approximation of the solution at the boundary elements is performed using the shape functions similarly to FEM. 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. The BEM has its origins in the early 1980s [Crouch and Starfield, 1983]. This method has found applications in underground excavations, dynamic rock problems, analysis of in situ stress and borehole drilling. 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 simulated in the works of [Blandford et al., 1981, Mi and Aliabadi, 1992], using domain subdivision into subdomains or by Dual BEM. The notable variants of BEM are also Galerkin BEM (GBEM) [Bonnet et al., 1998] and Dual-reciprocity BEM (DRBEM) [Kontoni, 1992].

**Discrete Element Method (DEM)**

The DEM started to develop in the field of the rock mechanics applications due to its requirements for the modelling of discontinuous behaviour [Cundall, 1971]. 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. The basic concept is to threat 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.

The constant evolution of DEM led to many different numerical approaches 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 moves and interacts between each other including the fragmentation process etc. All of DEM based methods have the similar basic approach, while the difference 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. The overview of the contact models can be found in [Wriggers, 2006].

The two main approaches of DEM exist according to the numerical integration scheme that can be explicit or implicit. 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 [Cundall, 1988]. 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, nuclear waste disposal, rock slopes, acoustic emission in rock, boreholes, laboratory test simulations etc.

The most well-known method of implicit DEM approach is Discontinuous Deformation Analysis (DDA) [Shi, 1992]. 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. The interaction between the independent blocks is performed by contact algorithms. It is an attractive method for solving rock mechanics problems where the rock blocks interact and deform between each other.

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, tunnelling, hydraulic fracturing.

Another approach similar to DEM is the Discrete lattice network methods used to simulate rock fracture initiation and propagation [Ostoja-Starzewski, 2002]. 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. When the cracks propagate through the domain, cohesive links are usually removed and remeshing is performed to simulate the evolving discontinuities.

**Aims, scopes and methodology**

The main scientific goal of this thesis is to provide the enhanced predictive models for localized failure in rock mechanics, by taking into account material heterogeneities and pre-existing defects. Localized failure represents the typical type of rock failure where full set of failure modes is needed: modes I and II in 2D setting and modes I, II and III in 3D setting. Mode I stands for tensile opening leading to development of macrocracks, while modes II and III represent shear sliding along the cracks. However, the true rock failure consists of two dissipative mechanisms: fracture process zone with microcracks and localized failure mechanism with macro-cracks. The aim of the thesis is to combine these two mechanisms within the single model, which depends on fine scale heterogeneities.

The novel approach is proposed here, by considering the rock material as a assembly of grains held together by cohesive forces. For this purpose, the discrete lattice model approach is chosen for rock representation, while cohesive links are enhanced with additional discontinuous properties. Such an enhancement is enabled through the embedded discontinuity framework (ED-FEM).

Another goal pertains to understanding the role of the fluid-flow to the development of micro-cracks and macro-cracks during the localized failure process. This leads to a complex multiphysics problem of volumetric fluid structure-interaction where all of the ingredients of the rock mechanics model need to be combined with micro-porosity of rock which allows for flow through the porous medium.

**Table of contents :**

Contents

List of Figures

List of Tables

**1 Introduction **

1 Motivation

2 Overview of the numerical methods applied to rock mechanics

2.1 Finite Difference and Finite Volume Method (FDM, FVM)

2.2 Finite Element Method (FEM)

2.3 Meshless methods

2.4 Boundary Element Method (BEM)

2.5 Discrete Element Method (DEM)

2.6 Discrete Fracture Network (DFN)

2.7 Hybrid methods

3 Aims, scopes and methodology

4 Outline

**2 2D rock mechanics model **

1 Model description

2 2D model formulation

2.1 Enhanced kinematics

2.2 Equilibrium equations

2.3 Constitutive model

2.4 Computational procedure

2.5 Anisotropic model with multisurface criterion

3 Numerical results

3.1 Validation test on beam with embedded discontinuities

3.2 Preparation of 2D plain strain rock specimens

3.3 Numerical tests on heterogeneous specimens

4 Final comments on the presented 2D rock mechanics model

Rock mechanics and failure phenomena

**3 3D rock mechanics model **

1 Model description

2 3D model formulation

2.1 Kinematics of strong disontinuity

2.2 The finite element approximation

2.3 Virtual work

2.4 Constitutive model

2.5 Computational procedure

2.6 The global solution procedure

2.7 The failure criteria

3 Numerical examples

3.1 Embedded discontinuity beam test

3.2 The construction of rock specimens

3.3 Rock specimen under uniaxial (unconfined) compression and tension test

3.4 Influence of pre-existing defects

4 Final comments on the presented 3D rock mechanics model

**4 Influence of specimen shape deviations on uniaxial compressive strength **

1 Research motivation

2 Preparation of specimens with targeted shape deviations

3 Numerical model

4 Numerical results

5 Final comments on research of shape deviations influence to the UCS

**5 2D model for failure of fluid-saturated rock medium **

1 Model introduction

2 The porous media formulation

3 The discrete lattice model description

3.1 The discrete lattice mechanical and fluid flow formulations

3.2 The strong discontinuities in poro-elastoplastic solid

4 The enhanced finite element formulation

4.1 The finite element interpolations

4.2 The enhanced weak form

4.3 The finite element equations of a coupled poroplastic problem .

4.4 The operator split algorithm

5 Numerical simulations

5.1 Uncoupled fluid flow across the lattice

5.2 Uniaxial coupled poroelastic problem

5.3 Drained compression test of the poro-plastic sample with the localized failure

6 Final comments on the presented failure model of fluid-saturated rock medium

**6 Conclusions and future perspectives **

1 Conclusions

2 Future perspectives

**Bibliography **