Computational methods for fracture in rock: a review and recent advances

Abstract

We present an overview of the most popular state-of-the-art computational methods available for modelling fracture in rock. The summarized numerical methods can be classified into three categories: Continuum Based Methods, Discrete Crack Approaches, and Block-Based Methods. We will not only provide an extensive review of those methods which can be found elsewhere but particularly address their potential in modelling fracture in rock mechanics and geotechnical engineering. In this context, we will discuss their key applications, assumptions, and limitations. Furthermore, we also address ‘general’ difficulties that may arise for simulating fracture in rock and fractured rock. This review will conclude with some final remarks and future challenges.

Keywords

numerical modelling / method development / rock mechanics / fractured rock / rock fracturing

Introduction

The importance of rock mass behavior is important in civil and underground engineering as most of the construction (e.g., dams, roads, etc.) are built on rock or soil foundation. Studying the rock behavior is challenging due to uncertainties and limited knowledge in the input parameters such as pre-existing faults and micro-cracks inside the rock mass, just to name a few. Although many experimental methods have been developed, a complete scanning over rock masses is not yet possible. Proper laboratory or in situ examination are costly or in some cases infeasible. Numerical modeling of rock reduces costs and complements experimental testing. In last decades, an enormous amount of work has been reported on modeling fracture. This paper provides an overview of the most popular computational methods for fracture in rock (in the absent of fluids). Rock is an anisotropic, heterogeneous material and commonly modeled in the framework of continuum mechanics which cannot capture fracture and the associated fractured rock kinematics.

Computational methods for fracture can be classified into three categories, namely Continuum Based Methods, Discrete Crack Approaches, and Block-Based Methods. Brittle fracture in continuum mechanics can be modeled by Linear Elastic Fracture Mechanic (LEFM). LEFM is only applicable to problems where the failure process zone (l) is less than one percent of the smallest critical dimension (D) of the structure. When the failure process zone compared to the smallest critical dimension of the structure is between one and twenty percent [¹,²], suggest to employ a (nonlinear) quasi-brittle fracture approach which accounts for the energy dissipation at postlocalization. Most applications in rock require a quasi-brittle approach to fracture. In the next section, continuum based methods for fracture are summarized. Section 3 presents discrete crack methods while Section 4 is devoted to all the methods that divide the domain into blocks and study their kinematics. The article closes with future approaches and some final remarks.

Continuum based methods

Quasi-brittle fracture (e.g., rock) refers to a stress drop after the yield stress (softening behavior) in the stress-strain diagram. This softening behavior is a sign of the reduced effective cross-sectional area due to initiation, coalescence, and growth of micro cracks, see e.g., Ref. [³]. for more details. Rock fracture is mostly quasi-brittle.

The underlying boundary value problem can be solved either in the strong form or in the weak form. Due to difficulties in imposing natural boundary conditions, most of the computational approaches for quasi-brittle fracture in rock are based on the weak form. Subsequently, the focus will be on computational methods based on the weak form or hybrid methods such as the Finite Volume-Finite Difference Method. The overview will begin with continuous approaches for fracture which smear out the crack over a finite width followed by discrete fracture approaches mainly in the context of Galerkin methods.

Combined Finite Volume, Finite Difference Method (FDM-FVM)

The Finite Volume Method (FVM) discretizes the integral form of the PDEs over the problem domain with a finite number of control volumes. Gauss divergence theory is employed in order to transform the domain integrals to boundary integrals. A set of algebraic equations can replace the integral form of PDE calculating the objective function of PDE at centre or vertexes of control volumes. One can combine the Finite Difference Method (FDM) – which is based on the strong form – with the FVM (FDM-FVM) by using a finite volume formulation for discretizing the problem domain and then solve the equilibrium equations using the neighbourhood control volume displacements in an FDM scheme [⁴–⁶].

The FDM-FVM cannot easily capture strong discontinuities and it is quite well-known that FDM-FVM crack propagation schemes are highly sensitive to meshing. Nonetheless, compared to a pure FDM, FDM-FVM has advantages with respect to meshing and applying boundary conditions [⁵–⁸]. Also, material nonlinearities can be modeled [⁹] more reliably. Different material properties in control volumes can be employed for heterogeneous material [⁶], geomechanical simulations [¹⁰], and fracture problems in rock [¹⁰–¹⁴].

Finite Element Method (FEM)

The robustness, reliability, and possibility to deal with complex geometries make the FEM one of the most popular methods for applications in geomechanics [¹⁵–²⁰]. While FEM is not well-suited to deal with discrete fracture, continuous approaches to fracture can ideally be implemented in the FEM framework. Note that the contributions described subsequently are not restricted to the finite element method and can theoretically be used also in other discretization techniques such as meshless methods.

Smeared crack models [²¹] which has also been used to model fracture in rock [²²–²⁷] horizontally scale the stress-strain curve with a characteristic element length ensuring the correct amount of energy dissipation after the material becomes unstable. They are easy to implement but are not true localisation limiters. The material instability can be prevented by using regularization techniques such as higher order continuum modes, nonlocal models [²⁸], gradient based models [²⁹–³¹] (also known as strongly [²⁹] or weakly [³] nonlocal models), polar models (e.g. Cosserat continuum [³²]), viscous models [³³,³⁴] and cohesive zone models [³⁵]; the latter models are better suited in combination with a discrete crack approach. The implementation of Cosserat continua for rock for instance is described by Mühlhaus et al. [³⁶–³⁸]. The enhanced review of Rabczuk [³⁹] provides a detailed overview, and therefore a thorough discussion is omitted in this manuscript. What all of the mentioned approaches above have in common is, they introduce a characteristic length and therefore cannot capture the discontinuity in the displacement field [⁴⁰–⁴⁷]. The ‘crack width’ is captured with several elements which makes the approach computationally expensive. If the global response is the main interest, continuous approaches for fracture are to our experience as well suited as discrete crack approaches. However, there are several interesting applications in geotechnical engineering such as hydraulic fracturing or compressed air energy storage [⁴⁸] which require modeling the fluid flow through discrete cracks. This cannot be accomplished with continuous approaches to fracture, and the fluid flow is commonly modeled by linking the permeability to the ‘damage’ of the material. When more fine-scale features need to be captured around the crack tip such as in heterogeneous materials, it might be feasible to employ multiscale methods for fracture as proposed in Refs. [⁴⁹–⁵⁵].

Discrete Crack Approaches

Element erosion and remeshing

One of the simplest approaches is to zero the stresses of the failed element(s) [⁵⁶]. However, such methods are quite sensitive with respect to the underlying discretization. Interelement-separation Methods are a better alternative [⁵⁷–⁵⁹] but they are also sensitive to the discretization and allow only crack propagation along predefined element boundaries unless remeshing techniques are exploited [⁶⁰–⁶³]. Very efficient remeshing techniques based on edge rotation have recently proposed by Areias et al. [⁶⁴–⁶⁹]; also in combination with continuous approaches to fracture [⁷⁰–⁷²].

Interelement-separation methods are often combined with cohesive zone models which ensure the well-posedness of the boundary value problem and account for the energy dissipation at postlocalization. They were frequently applied to dynamic fracture of various materials including rock; see, e.g., the contribution in Ref. [⁷³]. One advantage of remeshing techniques – compared to enriched FE formulations described in the next sections – is that they can be based on existing element which are well tested. If enrichment techniques are employed, the performance of the elements still need to be proven and it is often not simple to deal with problems with constraints, e.g., in plasticity or thin shell analysis. Also, the combination with cohesive elements is simpler in remeshing techniques when the crack is not located inside an element. On the other hand, remeshing in 3D might be complicated though much progress on meshing complex geometries have been made in the past decades. Such approaches have also frequently been developed and employed to large deformation problems including dynamic fracture and fragmentation.

Embedded Element Methods (EFEM)

In FEM, the displacement jump can be captured by enrichment of displacement field inside an element. Enrichment of displacement field is accomplished by including additional parameter(s) which are inherent to the element and therefore can be condensed at the element level. The name ‘Embedded Element’ is associated with the fact that the enrichment zone is completely embedded inside one element. The enrichment of the strain field for obtaining the kinematic of an element with one weak discontinuity [⁵⁹], two weak discontinuities [⁷⁴], and one strong discontinuity [⁷⁵] are the origins of the Embedded Elements Method (EFEM). The implementation of the EFEM has several advantages: 1) No remeshing is necessarily for capturing the crack growth; 2) The enrichment parameters are included at the element level. Therefore, its implementation in an existing FEM code straightforward. 3) Some approaches do not require crack path continuity [⁷⁶,⁷⁷] which facilitates tracing the crack path.

Since the element rigid body motions are not guaranteed, and elements are sensitivity to the direction of the discontinuity, the piece-wise constant crack opening of EFEM introduces some inaccuracies including stress locking, violation of traction continuity or an inaccurate representation of the displacement field [⁷⁸]. Recent developments aim to avoid these shortcomings [⁷⁹–⁸¹]. EFEM has also been employed for problems in geomechanics [⁸²] including complex failure mechanisms and multiscale analysis of rock [⁸³,⁸⁴]. While they are frequently applied to small strain theories, there are only a few contributions to large deformation problems and also coupled problems such as fluid-flow through discrete cracks.

Extended finite element method (XFEM)

The extended finite element method (XFEM) [⁸⁵,⁸⁶] introduces additional nodal parameters utilizing the local partition of unity [⁸⁷] and allows crack propagation without remeshing similar to EFEM. However, in contrast to EFEM, the additional degrees of freedom cannot be condensed at the element level. The principle idea of XFEM is shown in Fig. 1 [³⁹]. While the early applications for XFEM were focused on fracture, it has later on been applied to numerous areas such as biofilm growth, two-phase flow, fluid-structure interaction, inverse problems and optimization to name a few [⁸⁸–⁹²]. Applications to rock have been reported in [⁹³–⁹⁶]. XFEM also shows potential to model interesting coupled problems in rock mechanics coupling fluid flow through the propagation of cracks in porous media [⁹⁷]. However, due to lack of reliable criteria for crack branching and interaction, reliable predictions of complex fracture patterns in such problems remains a key challenge. In general, XFEM seems to be well suited for problems involving a few to a moderate number of cracks.

A special case of XFEM is the so-called phantom node method [⁹⁸]. The phantom node method employs overlapping elements in order to capture the strong discontinuity. This drastically facilitates the implementation of the method. On the other hand, it restricts its extension to other applications than fracture. Moreover, the crack is requested to propagate through an entire element each time step though there are some approaches in 2D which allow the crack to end inside an element [⁹⁹,¹⁰⁰] which is important for instance for heterogeneous materials. However, the extension of such crack tip elements to 3D seems cumbersome. XFEM, as well as the phantom node method, has been implemented into ABAQUS. An interesting approach is to combine isogeometric analysis [¹⁰¹] and XFEM as presented for instance in Refs. [¹⁰²,¹⁰³]. Such an approach allows for modelling cracks in so-called “exact geometries”, e.g., conic cross sections.

Numerical manifold method (NMM)

The numerical manifold method (NMM), also known as the finite cover method (FCM), is based on the concept of finite covers and the partition of unity (PU) and shows some similarities to the phantom node method. The method is formulated using a topological manifold and a differentiable manifold with two independent meshes, namely, a mathematical mesh and a physical mesh, see Fig. 2. It provides a unified framework to analyze continuum and discontinuum without changing predefined mesh in a discretized way. The NMM has been applied in the modeling of fluid structure interaction as well as in rock mechanics including the analysis of block system, jointed rock and fractured body, showing particular advantages over other PU based methods. Unlike other PU methods, the degrees of freedoms in the NMM are associated with the physical covers, rather than the nodes, which allow it to be naturally adapted to the changing geometries in analyzing complex discontinuum such as multiple intersecting cracks and branched cracks. It provides a natural solution by using the cover system to accommodate problem domains with strong or weak discontinuities, and unifies the continuum analysis and discontinuum analysis [¹⁰⁴–¹⁰⁶] the NMM has also attracted much interest from researchers in solid mechanics in recent years, e.g., simplex integration scheme in NMM [¹⁰⁷,¹⁰⁸]; the imposition of boundary conditions [¹⁰⁹–¹¹¹]; treatment of material discontinuities or arbitrary discontinuities [¹¹²]; modeling discontinuities or crack growth with NMM [¹¹³–¹¹⁵]; development of high accuracy manifold element [¹¹⁶]; development of 3D numerical manifold method (NMM) [¹¹⁷]; application of NMM to geotechnical or fluids problems [¹¹²].

In the NMM, the mathematical mesh is independently defined from the physical cover, and regular mesh type such as rectangle can be used. The shape of the mesh does not need to coincide with the problem domain. It is allowed to exceed over the size of the domain as shown in Fig. 2. The mathematical mesh is further partitioned by the physical lines, e.g., crack or joints into physical covers. The interpolation/approximation over the physical elements formed by the overlapping of the physical covers is the multiplication of the finite element interpolation over the mathematical mesh and the cover functions over the physical covers. The physical cover is defined by the intersection of mathematical cover and physical domain. The mathematical cover is the same as the physical cover if there are no material discontinuities or interfaces. Displacements are introduced at the physical cover and are represented as physical cover functions u^{{{\displaystyle L}}_i}\left(\mathit{X}\right), where i is the index of physical cover. For example, the physical cover associated with node 10 is shown in Fig. 2. The intersection of the physical cover is physical element. Within each element, displacement is approximated from the finite cover functions of its composing nodes. Finally, the finite cover approximation is obtained as

(1)u^e\left(\mathit{X}\right)={\displaystyle {\sum }_{i=\mathrm{1}}^{{{\displaystyle n}}_e}{{\displaystyle \omega }}_i}{u}^{{{\displaystyle L}}_i}\left(\mathit{X}\right),i=\mathrm{1},...,{{\displaystyle n}}_e,

and {{\displaystyle \omega }}_i is the weight function related with the ith physical cover and {{\displaystyle n}}_e is the number of associated physical covers. The weight function in Eq. (1) should be non-negative and conforms to partition of unity, i.e., {\displaystyle {\Sigma }_{i=\mathrm{1}}^{{{\displaystyle n}}_e}{{\displaystyle \omega }}_i}=\mathrm{1}. Here, linear weight functions in 2D, which are exactly the same the standard FEM, are used. For an arbitrary physical cover e shown in Fig. 2 the displacement is approximated by nodes 9, 10, 13, 14. The NMM can be regarded as a special form of the PU methods. The cover functions constructed by the NMM are free to take different orders of expansion for specific parts of the problem, and different types of analytical expansions in subdomains are possible. The displacements along the interfaces of subdomains are naturally consistent.

Apart from the theoretical advantages, the NMM has several issues that practically limit its wider applications. Firstly, when the order of the basis function increases, the number of unknowns increases significantly. Besides, the imposition of the boundary conditions is not straightforward. Secondly, the cover generation is not a simple task in the NMM, but it is an essential part of the NMM. Existing literature of the NMM have referred to various aspects of the finite cover generation such as the numbering of physical covers. However, there are very limited works that describe the automatic physical cover generation of the NMM in a systematic way or giving a general principle of cover numbering. Indeed, the absence of a general theory and algorithms has become the bottleneck of the NMM especially for 3D applications of the NMM.

Meshfree methods (MMs)

Due to the absence of a mesh, Meshfree Methods (MMs) are also suitable to model arbitrary crack growth. There are different approaches to model fracture in meshfree methods which also depends on the application. For quasi-static fracture problems with a small to moderate number of cracks, the visibility, transparency or diffraction method [¹¹⁸] can be exploited though partition of unity approaches similarly to XFEM have been exploited in meshfree methods as well [⁹⁹,¹¹⁹–¹²⁷]. For dynamic fracture and fragmentation, i.e., applications with a huge number of cracks, simple particle splitting approaches can be employed. A very simple and efficient method to deal with fracture in meshfree methods is the cracking particles method (CPM) [¹²⁸–¹³⁰]. The CPM describes the crack as a set of crack segments. It requires no track cracking procedures and complex fracture patterns such as crack branching and coalescence is a natural outcome of the simulation. Hence, the CPM is particularly suited for modeling fracture in rock with extensive micro-cracking, see e.g., the contribution in Ref. [¹²⁵]. It has also frequently been developed for large deformation problems and finite strain theory. However, due to lack of crack path continuity, fluid flow through discrete cracks cannot be captured straight forward in such approaches. An interesting CPM method which enforces crack path continuity in 2D has been proposed by Ai and Augarde [¹³¹] and could be a good pathway to address such problems. Another approach by Rabzuck et al. [¹³²] is capable of simulating 3D fluid-structure interactions. Other applications of MMs to fracture in rock have been done for instance [¹³³–¹³⁵].

Peridynamics

Peridynamics (PD) is a nonlocal continuum mechanics theory proposed by Silling [¹³⁶]. One of the key advantages of PD over methods such as XFEM is that the crack is a natural outcome of the formulation. This is achieved by rewriting the equation of motion and substituting the divergence of the Cauchy stress tensor with an integral expression which contains the so-called bond forces. The first version proposed in 2000 is called bond based peridynamic (BB-PD). It allows the symmetric interactions between material points. BB-PD has several restrictions. For instance, it allows only material models with a Poisson ratio of 0.33 in 2D and 0.25 in 3D. Furthermore, the incorporation of arbitrary material models is quite cumbersome to impossible as the physical interpretation of the bond force vector is complicated, especially with increasing complexity of the constitutive behavior. Therefore, the so-called state-based peridynamic (SB-PD) has been develop [¹³⁷]. In SB-PD, the bond force is related to continuum mechanics stresses and hence allows for the incorporation of more complex material models. Also, the Poisson ratio limitation has been removed. On the other hand, several papers have shown extreme similarities between SB-PD and specific meshfree methods (SPH and RKPM) [¹³⁸–¹⁴⁰]. Hence, all difficulties inherent to meshfree methods such as instabilities due to nodal integration should also apply to SB-PD.

SB-PD can be further categorized into ordinary state based and non-ordinary state based PD. All PD formulations requires a uniform discretization which leads to high computational cost. It has been shown [¹⁴¹] that the introduction of refined areas leads to artificial wave reflections and high sensitivity in the crack path. The dual horizon peridynamic (DH-PD) [¹⁴²,¹⁴³] drastically alleviates this issue by distinguishing between active and reactive bond forces. DH-PD can be applied to BB-PD as well as SB-PD. An illustrative figure briefly highlighting the key idea of the different PD formulation is shown in Fig. 3. PD has mainly been applied for dynamic fracture and fragmentation including fragmentation of rock [¹⁴⁴–¹⁴⁶]. Although, SB-PD facilitates defining some of the material behavior, defining shear behavior is still problematic [¹⁴⁷]. Also, the extension to fluid flow through discrete cracks seem more cumbersome, particularly when the fluid is to be modeled by a PD approach. All PD formulations requires a uniform discretization which leads to high computational cost [¹⁴⁸].

Boundary element method (BEM)

Another popular method for fracture is the boundary element method (BEM) which reduces the problem by one dimension (see Fig. 4). As the name indicates, in the BEM only the (crack) boundary is discretized and therefore simplifies meshing and remeshing during crack propagation. However, BEM requires an exact solution of the differential equation inside the domains. Thus, the applicability of BEM is limited to problems where Green’s functions can be computed. It cannot easily be applied to heterogeneous materials or materials with complex nonlinear constitutive models. A method that can overcome this restriction is the scaled boundary element method [⁹⁸]. A good introduction to BEM is given in the textbook by Brebbia and Walker [¹⁴⁹]. While the BEM has been applied to fracture problems, see e.g., the contributions in [¹⁵⁰,¹⁵¹], its application to rock mechanics is scarce due to mentioned restrictions above [¹⁵²,¹⁵³].

Block-based methods

Block based methods (BBMs) model the material or structure through interaction/contact of an assembly of blocks. The idea behind this class of methods is to bound of an assembly of blocks to each other using contacts. BBMs allow for simulating large displacement of fractured rock and studying such systems [¹⁵⁴,¹⁵⁵]. Subsequently, we review four of the most popular BBMs, namely Discrete Element Method (DEM), Key Block Theory, Discontinues Deformation Analyze (DDA) and Lattice Network Model. Most of these methods are employed to mechanical problems, and their implementation require the following five steps:

1) Defining geometry configuration of the rock blocks in 2D or 3D.

2) Computation of the rock blocks deformations.

3) Contact detection between rock blocks.

4) Calculating the contact forces regarding the contact constitutive equations or law.

5) Determining the blocks positions and update their defined value in the first step.

Each method may have different sequences (e.g., DDA) due to their different approaches to the solution of the system. Some methods also consider the blocks to be rigid while others allow for deformable blocks. However, all methods require efficient contact detection models, and the material response is mainly governed by the interaction of the blocks.

Discrete element method (DEM)

The ‘Discrete element method’ (DEM) has been originated by Cundall [¹⁵⁶,¹⁵⁷]. In the DEM, the material behavior is described by the interaction of rigid and/or deformable bodies of arbitrary shape. Hence, fracture can be modeled by the DEM quite easily. The method seems most suitable for granular materials or to simulate the movement of fractured rock masses. In our opinion, it is not well suited to model fracture for any other material or application. One key difficulty of the DEM is the choice of the contact laws and the calibration of its material parameters. Reliable predictions of a variety of mechanical properties (for one realization) remains a major challenge and procedures to uniquely determine the material parameters for the discrete elements in order to capture a variety of macroscopic properties are missing. There is not a single manuscript which predicts macroscopic properties such as Young’s modulus, Poisson’s ratio, density, multi-axial (static and dynamic) strength of materials, fracture toughness and macroscopic measured yield and failure surfaces in 3D for complex material behavior based on a single realization. Most contributions report only a few macroscopic material parameters and a thorough sensitivity analysis quantifying the influence of the material and geometry input parameters is missing. However, it is expected that the results obtained by DEM will be quite sensitive not only with respect to the geometry and sizes of the discrete elements but also their statistical distribution. Even for granular materials, it is often unfeasible to account for small grains, so they are commonly neglected. Simulation of large 3D structures requires an efficient implementation. One of the key issue for such an effective implementation are fast contact detection algorithms. There are three types of contacts in 2D (vertex-to-vertex, vertex-to-edge, and edge-to-edge) and six in 3D (face-to-face, face-to-edge, face-to-vertex, vertex-to-edge, edge-to-edge, and vertex-to-vertex). Most of the contact detection algorithms focus on contact of convex bodies and the literature on general (concave and convex) bodies [¹⁵⁸] is scarce. A good overview on contact detection algorithms is given for instance in the textbook by Wriggers [¹⁵⁹]. Furthermore, nearly all DEM formulations employ an explicit time integration scheme which requires an upper limit on the time step. For some applications this increases the simulation time as well. A comprehensive report of DEM can be found in Refs. [¹⁶⁰,¹⁶¹].

As mentioned above, discrete elements can also be deformable. In principle, each discrete element can be discretized with finite elements which would also allow fracture inside the discrete elements. In this case, fracture criteria as employed in finite elements are needed. On the other hand, on the cost of accuracy simpler criteria are often exploited when the discrete element is assumed to be under constant strains. A linear displacement field is commonly assumed inside a discrete element since higher order elements or carved edge elements increase the complexity of the contact detection algorithm. Note that, the deformability at DEM leads to extra computations and decreases the critical time steps compared to rigid Rigid-Blocks Cluster Methods.

The contact laws govern the macroscopic behavior of the material. The contact is split into normal and tangential behavior. Linear, as well as nonlinear, contact laws have been developed in the past; many of those for granular materials [¹⁶²,¹⁶³] where the DEM has a more profound ‘physical basis.’ DEM solves the equation of motion with an explicit time integration scheme [¹⁶⁴]:

(2)m{{\displaystyle \ddot{u}}}_i^t+c{{\displaystyle \dot{u}}}_i^t={{\displaystyle F}}_i^t,

(3)I{\dot{\omega }}^t+c{\omega }^t=M^t,

where t presents time; i, direction in global system; u, displacement of element centre of gravity; w, angular velocity; c, viscose damping; m, mass; I, mass moment of inertia; F, resultant of applied force on element; and M presents applied momentum on the element centre of gravity. The extension to multi-physics problems (such as thermo-mechanical problems) which are important in some applications in geotechnical engineering seems difficult. Nonetheless, DEM is a very popular method to model fracture in fractured rock masses and rock mechanics including jointed rock masses [¹⁶⁵], hard rock reinforcement [¹⁶⁶], tunnelling [¹⁶⁷], underground excavations and mining [¹⁶⁸], well and borehole stability [¹⁶⁹], reservoir simulations [¹⁷⁰], earthquakes [¹⁷¹] and even coupled problems such as fluid injection [¹⁷²].

Key block theory

Key blocks (keystones) are blocks that have the potential to rotate or slide without any geometrical obstruction. Key block theory identifies the critical rock block which will lose stability first, and its motion will trigger potential sequential motion of other blocks. It was firstly proposed during the 1970s by Shi [¹⁷³]. The key block theory is a static geometry and topology analysis method in which no stress or deformation analyses is possible [¹⁷⁴–¹⁷⁶]. The method was developed with an aim to find the blocks that will become unstable. The method provides a practical stability analysis in hard rock, but since the blocks are assumed to be rigid, the method does not present a reliable result for soft rock stability. The key block theory is faster than other Block-Based methods. Including water effects [¹⁷⁷], the probabilistic predictions and Monte Carlo simulations for developing the fracture patterns [¹⁷⁸,¹⁷⁹], finite block size effect [¹⁸⁰], linear programming [¹⁸¹] and secondary blocks [¹⁸²] to the method, made the key block theory as an effective way for slope and tunnel surround rock stability analysis in hard rock conditions [¹⁸³,¹⁸⁴]. However, the theory cannot be used to model any advanced material behavior or detailed fracture pattern.

There are two different ways to implement the key block theory, namely stereographic projection method and vector methods. Conventional key block theory employs stereographic projection method to transfer three-dimensional objects to two-dimensional shapes. A reference sphere with centre at the coordinate system’s origin is generated. For a joint plane or free plane of a block, the dip direction and dip angle are ascertained by geological survey. The relative position relation of the plane and the block can be found by defining the outer normal vector. Warburton (1981) [¹⁸⁵] presented a vector method for rock. Goodman and Shi (1985) [¹⁷⁵] elaborated this theory systematically and provided several calculation examples using this method. Vector method describes block morphology with vector equations and matrices. This method provides an analytical solution using matrix operations based on vector analysis.

Discontinues deformation analyze (DDA)

Rocks found in nature contains a large amount of joints and fissures. These discontinuities, depending on their density, spatial distribution, and extensions, cut the rock into different formations and different block sizes of blocks. For heavily jointed rock with a dense distribution of joints, the rock is regarded continuum and can be modeled by using reduced material parameters compared to the intact rock. For rock with few dominant directions of discontinuities and relative large spacing between the discontinuities, large blocks can be formed. To analyze the kinematics, stability, and motion of rock blocks according to their shape, topology and spatial locations, will be of primary importance in rock engineering, e.g., for tunnel and slope.

For analyzing the motion and kinematics of rock mass, the discontinuous deformation analysis (DDA) method was originally proposed in the 1980s by Shi and Goodman [¹⁷⁶,¹⁸⁶] in 2D and later in 3D [¹⁸⁷]. It can be regarded as an extension of the key block theory proposed to the kinetics of rock blocks. DDA is a discontinuum based numerical method, especially applicable to solve large displacement and deformation problems in a block system. Accordingly, DDA has been applied to comprehensive joint surface representation [¹⁸⁸] and full internal discretization of blocks [¹⁸⁹,¹⁹⁰]. Extensions to fracture and fragmentation [¹⁹¹], coupled flow-stress analysis [¹⁹²], 3D block system [¹⁹³] and higher order displacement field to allow non-constant strains in each rock block [¹⁹⁴,¹⁹⁵] were proposed later.

In DDA, each rock block is treated as an elastic body with constant strain through the entire block. The motion and deformation of the block are represented by the geometric centroid of the block (assuming the rock block comprises single material). The condition of linear elasticity, uniform stresses and strains are assumed for the whole process of simulation. The movement and deformation of the block are defined by 12 independent deformation variables in the deformation matrix. The movement is described by 6 degrees of freedom in three-dimensional space, namely 3 translation degrees of freedom, and 3 rotational degrees of freedom. The deformation of the block is described by 6 strain components, namely 3 normal and 3 shear strain components.

The key challenge in large displacement modeling using DDA is the efficiency and robustness of the contact algorithm. The contact is decomposed into normal and tangential directions [¹⁹⁶,¹⁹⁷]. Contact detection includes the searching of neighbour blocks and detection of contact patterns, i.e., to find candidate neighbouring blocks and identify contact groups of vertices, edges and facets (geometrical elements) in a multi-block system. Contact relations between discontinuous objects are complicated especially in 3D. It requires high computational cost in contact detection of large scale problem. For 3D particle geometries, contact detection can take up to 80% of the total analysis time [¹⁹⁸]. Thus, the computational speed and stability of 3D DDA strongly depend on the efficiency and robustness of the contact searching and definition algorithm, which affects the simultaneous equilibrium equations solved in each step of open-close iteration (OCI) [¹⁹⁹]. A quite efficient contact detection algorithm has recently been proposed by Wu et al. [²⁰⁰].

To enable realistic applications of DDA to rock engineering, model generation in 3D remains a primary challenge. Since the major part of joint lies inside the rock mass and is usually invisible, it is difficult to obtain the joint information directly though outcrop though boreholes can assist in providing some information with a limited number of samples. Apart from the uncertainties of the spatial distribution of discontinuities and joint surface (interface) parameters, the complex intersecting relations of joint surfaces in 3D result in various sizes and shape of blocks. Various efforts have been made to find and determine the shape of the blocks from the site including photogeometry technique and 3D surface reconstruction technique. In the photogeometry technique, photos are taking at the boring face of a tunnel or outcrop of slopes [¹⁵⁸]. From the photo obtained, processed and analyzed, then location, dip direction and dip angle of joints as well as other geometric information of tunnels can be identified. The identified information is transferred as the input to generate 3-dimensional discontinuous models for DDA analysis. The models can be analyzed fast using key blocks theory and DDA to evaluate the stability condition of a tunnel.

While DDA and DEM share many similarities such as the contact detection, their key difference can be summarized in 4 points:

1) In contrast to DEM, DDA minimizes the potential energy of the whole system by solving a set of system equations similar to FEM.

2) DDA calculates the displacement directly from the stiffness of the system while DEM finds the forces and stresses from relative positioning of blocks.

3) DDA allows no block overlapping.

4) DDA is also known as an implicit form of DEM (explicit) because of the solution algorithm and geometry topology.

On the other hand, the application area of DDA is similar to the application area of DEM.

Lattice network model

The network of concentrated masses that are connected with the massless springs (element edges) is the domain of interest in this model. The generation of the masses and spring stiffness gives the possibility to define stochastic, inhomogeneous, and statistical inhomogeneous medium. This fact defers the lattice methods from other BBMs and the continuum methods. A dynamic local equation of motion will be solved over the domain resulting the deformation of the springs. The integration of the rock material through hydro-fracture procedure is the main application of lattice method in rock fracturing and fractured rock interaction [²⁰¹–²⁰⁴].

Conclusion and future challenges

A literature review on numerical methods for fracture and their application to rock fracturing and fractured rock has been provided. The methods were classified in three categories:

• Continuum Based Methods,

• Discrete Crack Approaches, and

• Block-Based Methods.

Continuous methods and discrete crack methods seem to be most suitable for modeling fracture in rock while they seem less suitable for the interaction of fractured rock which plays a key role for instance in avalanches. Continuous approaches to fracture introduce a characteristic length scale and spread the crack width over a certain domain. They are good candidates when global responses are of interest as they occur in large rock system. They seem less suited for applications where detailed information around a crack tip is of interest.

BBMs allow simulations of rock fracturing as well as the interaction of fractured rock with ease. Most BBMs model the macroscopic behavior through interaction of rigid or deformable blocks. Deformable BBMs commonly have assumed quite simple kinematics in each block; commonly constant strains though there are also more complex models as discussed in Section 4 (Block-based methods). The macroscopic response is therefore highly sensitive with respect to the shapes, sizes, and distribution of the blocks as well as the contact laws between them. Calibrating these parameters reliably and uniquely remains one major challenge of those BBMs.

References

[1] Carpinteri A. Post-peak and post-bifurcation analysis of cohesive crack propagation. Engineering Fracture Mechanics, 1989, 32(2): 265–278

[2] Planas J, Elices M. Nonlinear fracture of cohesive materials. In: Current Trends in Concrete Fracture Research, Springer, 1991, 139–157

[3] Bažant Z P, Jirásek M. Nonlocal integral formulations of plasticity and damage: survey of progress. Journal of Engineering Mechanics, 2002, 128(11): 1119–1149

[4] Wheel M. A geometrically versatile finite volume formulation for plane elastostatic stress analysis. Journal of Strain Analysis for Engineering Design, 1996, 31(2): 111–116

[5] Selmin V. The node-centred finite volume approach: bridge between finite differences and finite elements. Computer Methods in Applied Mechanics and Engineering, 1993, 102(1): 107–138

[6] Fallah N, Bailey C, Cross M, Taylor G. Comparison of finite element and finite volume methods application in geometrically nonlinear stress analysis. Applied Mathematical Modelling, 2000, 24(7): 439–455

[7] Bailey C, Cross M. A finite volume procedure to solve elastic solid mechanics problems in three dimensions on an unstructured mesh. International Journal for Numerical Methods in Engineering, 1995, 38(10): 1757–1776

[8] Mishev I D. Finite volume methods on voronoi meshes. Numerical Methods for Partial Differential Equations, 1998, 14(2): 193–212

[9] Fryer Y, Bailey C, Cross M, Lai C H. A control volume procedure for solving the elastic stress-strain equations on an unstructured mesh. Applied Mathematical Modelling, 1991, 15(11–12): 639–645

[10] Detournay C, Hart R. FLAC and numerical modelling in geomechanics. In: Proceedings of the International FLAC symposium on Numerical Modelling in Geomechanics, Minneapolis. Rotterdam: Balkema, 1999

[11] Fang Z. A local degradation approach to the numerical analysis of brittle fractures in heterogeneous rocks. Dissertation for PhD degree. Imperial College London (University of London), 2001

[12] Martino S, Prestininzi A, Scarascia Mugnozza G. Mechanisms of deep seated gravitational deformations: parameters from laboratory testing for analogical and numerical modeling. In: Proc. Eurock, 2001, 137–142

[13] Kourdey A, Alheib M, Piguet J, Korini T. Evaluation of slope stability by numerical methods. The 17th International Mining Congress and Exhibition of Turkey, 2001, 705–710

[14] Marmo B A, Wilson C J L. A verification procedure for the use of FLAC to study glacial dynamics and the implementation of an anisotropic flow law. In: Särkkä P, Eloranta P, eds. Rock Mechanics—A Challenge for Society. Lisse: Swetz and Zeitlinger, 2001

[15] Jing L, Hudson J. Numerical methods in rock mechanics. International Journal of Rock Mechanics and Mining Sciences, 2002, 39(4): 409–427

[16] Wittke W, Sykes R. Rock Mechanics. Springer Berlin, 1990

[17] Peng W. The damage mechanics model for jointed rock mass and its nonlinear FEM analysis. Chinese Journal of Rock Mechanics and Engineering, 1988, 7(3): 193–202 (in Chinese)

[18] Zheng Y R, Zhao S Y. Application of strength reduction FEM in soil and rock slope. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(19): 3381 (in Chinese)

[19] Zhao S, Zheng Y, Deng W. Stability analysis on jointed rock slope by strength reduction FEM. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(2): 254–260

[20] Cai M, Horii H. A constitutive model and FEM analysis of jointed rock masses. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1993, 30(4): 351–359

[21] Bazănt Z P, Cedolin L. Fracture mechanics of reinforced concrete. Journal of the Engineering Mechanics Division, 1980, 106(6): 1287–1306

[22] Goodman R E, Taylor R L, Brekke T L. A model for the mechanics of jointed rocks. Journal of Soil Mechanics & Foundations Division, 1968, 94(3): 637–660

[23] Zienkiewicz O C, Best B, Dullage C, Stagg K G. Analysis of nonlinear problems in rock mechanics with particular reference to jointed rock systems. In: Proceedings of International Society of Rock Mechanics, 1970

[24] Ghaboussi J, Wilson E, Isenberg J. Finite element for rock joints and interfaces. Journal of Soil Mechanics & Foundations Division, 1973, 99(10): 849–862

[25] Desai C, Zaman M, Lightner J, Siriwardane H. Thin-layer element for interfaces and joints. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(1): 19–43

[26] Goodman R E. Methods of geological engineering in discontinuous rocks. New York: West Publishing, 1976

[27] Katona M G. A simple contact–friction interface element with applications to buried culverts. International Journal for Numerical and Analytical Methods in Geomechanics, 1983, 7(3): 371–384

[28] Bažant Z P. Why continuum damage is nonlocal: micromechanics arguments. Journal of Engineering Mechanics, 1991, 117(5): 1070–1087

[29] Peerlings R H J, de Borst R, Brekelmans W A M, de Vree J H P. Gradient enhanced damage for quasi-brittle materials. International Journal for Numerical Methods in Engineering, 1996, 39(19): 3391–3403

[30] Peerlings R H J, de Borst R, Brekelmans W, Geers M. Localisation issues in local and nonlocal continuum approaches to fracture. European Journal of Mechanics. A, Solids, 2002, 21(2): 175–189

[31] de Borst R, Pamin J, Geers M G. On coupled gradient-dependent plasticity and damage theories with a view to localization analysis. European Journal of Mechanics. A, Solids, 1999, 18(6): 939–962

[32] Pasternak E, Dyskin A, Mühlhaus H B. Cracks of higher modes in Cosserat continua. International Journal of Fracture, 2006, 140(1–4): 189–199

[33] Etse G, Willam K. Failure analysis of elastoviscoplastic material models. Journal of Engineering Mechanics, 1999, 125(1): 60–69

[34] Rabczuk T, Eibl J. Simulation of high velocity concrete fragmentation using SPH/MLSPH. International Journal for Numerical Methods in Engineering, 2003, 56(10): 1421–1444

[35] Hillerborg A, Modéer M, Petersson P E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6(6): 773–781

[36] Miihlhaus H B, Triantafyllidis T. Surface waves in a layered half-space with bending stiffness. Developments in Geotechnical Engineering, 1987, 44: 277–290

[37] Mühlhaus H B. Application of Cosserat theory in numerical solutions of limit load problems. Archive of Applied Mechanics, 1989, 59(2): 124–137

[38] Vardoulakis I, Mühlhaus H. Local rock surface instabilities. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1986, 23: 379–383

[39] Rabczuk T. Computational methods for fracture in brittle and quasi-brittle solids: state of the art review and future perspectives. ISRN Applied Mathematics, 2013, 2013: 332–369

[40] de Borst R. Fracture in quasi-brittle materials: a review of continuum damage-based approaches. Engineering Fracture Mechanics, 2002, 69(2): 95–112

[41] Jirásek M, Zimmermann T. Analysis of rotating crack model. Journal of Engineering Mechanics, 1998, 124(8): 842–851

[42] Jirásek M, Zimmermann T. Rotating crack model with transition to scalar damage. Journal of Engineering Mechanics, 1998, 124(3): 277–284

[43] Rabczuk T, Akkermann J, Eibl J. A numerical model for reinforced concrete structures. International Journal of Solids and Structures, 2005, 42(5–6): 1327–1354

[44] Ohmenhäuser F, Weihe S, Kröplin B. Algorithmic implementation of a generalized cohesive crack model. Computational Materials Science, 1999, 16(1): 294–306

[45] Carpinteri A, Chiaia B, Cornetti P. A scale-invariant cohesive crack model for quasi-brittle materials. Engineering Fracture Mechanics, 2002, 69(2): 207–217

[46] François M, Royer-Carfagni G. Structured deformation of damaged continua with cohesive-frictional sliding rough fractures. European Journal of Mechanics. A, Solids, 2005, 24(4): 644–660

[47] de Borst R, Remmers J J, Needleman A. Mesh-independent discrete numerical representations of cohesive-zone models. Engineering Fracture Mechanics, 2006, 73(2): 160–177

[48] Zhuang X, Huang R, Liang C, Rabczuk T. A coupled thermo-hydro-mechanical model of jointed hard rock for compressed air energy storage. Mathematical Problems in Engineering, 2014, 179169

[49] Silani M, Talebi H, Hamouda A M, Rabczuk T. Nonlocal damage modelling in clay/epoxy nanocomposites using a multiscale approach. Journal of Computational Science, 2016, 15: 18–23

[50] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modeling of three-dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80: 82–92

[51] Silani M, Ziaei-Rad S, Talebi H, Rabczuk T. A semi-concurrent multiscale approach for modeling damage in nanocomposites. Theoretical and Applied Fracture Mechanics, 2014, 74: 30–38

[52] Talebi H, Silani M, Bordas S P, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[53] Budarapu P R, Gracie R, Bordas S P, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[54] Budarapu P R, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[55] Talebi H, Silani M, Bordas S P, Kerfriden P, Rabczuk T. Molecular dynamics/XFEM coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture. International Journal for Multiscale Computational Engineering, 2013, 11(6): 527–541

[56] Belytschko T, Lin J I. A three-dimensional impact-penetration algorithm with erosion. Computers & Structures, 1987, 25(1): 95–104

[57] Camacho G T, Ortiz M. Computational modelling of impact damage in brittle materials. International Journal of Solids and Structures, 1996, 33(20): 2899–2938

[58] Xu X P, Needleman A. Void nucleation by inclusion debonding in a crystal matrix. Modelling and Simulation in Materials Science and Engineering, 1993, 1(2): 111–132

[59] Ortiz M, Leroy Y, Needleman A. A finite element method for localized failure analysis. Computer Methods in Applied Mechanics and Engineering, 1987, 61(2): 189–214

[60] Pandolfi A, Krysl P, Ortiz M. Finite element simulation of ring expansion and fragmentation: the capturing of length and time scales through cohesive models of fracture. International Journal of Fracture, 1999, 95(1–4): 279–297

[61] Pandolfi A, Guduru P, Ortiz M, Rosakis A. Three dimensional cohesive-element analysis and experiments of dynamic fracture in c300 steel. International Journal of Solids and Structures, 2000, 37(27): 3733–3760

[62] Zhou F, Molinari J F. Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency. International Journal for Numerical Methods in Engineering, 2004, 59(1): 1–24

[63] Falk M L, Needleman A, Rice J R. A critical evaluation of cohesive zone models of dynamic fracture. Journal de Physique. IV, 2001, 11(PR5): Pr5-43–Pr5-50

[64] Areias P, Reinoso J, Camanho P, Rabczuk T. A constitutive-based element-by-element crack propagation algorithm with local mesh refinement. Computational Mechanics, 2015, 56(2): 291–315

[65] Areias P, Rabczuk T, Camanho P. Finite strain fracture of 2d problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[66] Areias P, Rabczuk T, Dias-da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[67] Areias P, Rabczuk T, Camanho P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[68] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[69] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[70] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh refinement. Computational Mechanics, 2016, 58(6): 1003–1018

[71] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 322–350

[72] Areias P, Msekh M, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[73] Gens A, Carol I, Alonso E. Rock joints: FEM implementation and applications. Studies in Applied Mechanics, 1995, 42: 395–420

[74] Belytschko T, Fish J, Engelmann B E. A finite element with embedded localization zones. Computer Methods in Applied Mechanics and Engineering, 1988, 70(1): 59–89

[75] Dvorkin E N, Cuitio A M, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. International Journal for Numerical Methods in Engineering, 1990, 30(3): 541–564

[76] Feist C, Hofstetter G. Three-dimensional fracture simulations based on the SDA. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(2): 189–212

[77] Sancho J M, Planas J, Fathy A M, Galvez J C, Cendon D A. Three-dimensional simulation of concrete fracture using embedded crack elements without enforcing crack path continuity. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(2): 173–187

[78] Jirásek M. Comparative study on finite elements with embedded discontinuities. Computer Methods in Applied Mechanics and Engineering, 2000, 188(1): 307–330

[79] Linder C, Zhang X. Three-dimensional finite elements with embedded strong discontinuities to model failure in electromechanical coupled materials. Computer Methods in Applied Mechanics and Engineering, 2014, 273: 143–160

[80] Linder C, Armero F. Finite elements with embedded strong discontinuities for the modeling of failure in solids. International Journal for Numerical Methods in Engineering, 2007, 72(12): 1391–1433

[81] Oliver J, Huespe A, Blanco S, Linero D. Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach. Computer Methods in Applied Mechanics and Engineering, 2006, 195(52): 7093–7114

[82] Foster C, Borja R, Regueiro R. Embedded strong discontinuity finite elements for fractured geomaterials with variable friction. International Journal for Numerical Methods in Engineering, 2007, 72(5): 549–581

[83] Nikolic M, Ibrahimbegovic A. Rock mechanics model capable of representing initial heterogeneities and full set of 3d failure mechanisms. Computer Methods in Applied Mechanics and Engineering, 2015, 290: 209–227

[84] Saksala T. Rate-dependent embedded discontinuity approach incorporating heterogeneity for numerical modeling of rock fracture. Rock Mechanics and Rock Engineering, 2015, 48(4): 1605–1622

[85] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[86] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131–150

[87] Melenk J M, Babuška I. The partition of unity finite element method: basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1): 289–314

[88] Rabinovich D, Givoli D, Vigdergauz S. Crack identification by arrival timeusing XFEM and a genetic algorithm. International Journal for Numerical Methods in Engineering, 2009, 77(3): 337–359

[89] Béchet E, Scherzer M, Kuna M. Application of the x-FEM to the fracture of piezoelectric materials. International Journal for Numerical Methods in Engineering, 2009, 77(11): 1535–1565

[90] Mayer U M, Gerstenberger A, Wall W A. Interface handling for three-dimensional higher-order XFEM-computations in fluid–structure interaction. International Journal for Numerical Methods in Engineering, 2009, 79(7): 846–869

[91] Verhoosel C V, Remmers J J, Gutiérrez M A. A partition of unity-based multiscale approach for modelling fracture in piezoelectric ceramics. International Journal for Numerical Methods in Engineering, 2010, 82(8): 966–994

[92] Nanthakumar S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[93] Zheng A X, Luo X Q, Shen H. Numerical simulation and analysis of deformation and failure of jointed rock slopes by extended finite element method. Rock and Soil Mechanics, 2013, 34(8): 2371–2376 (in Chinese)

[94] Wan L L, Yü T T. Pre-processing of extended finite element method for discontinuous rock masses. Rock and Soil Mechanics, 2011, 32: 772–778 (in Chinese)

[95] Goodarzi M, Mohammadi S, Jafari A. Numerical analysis of rock fracturing by gas pressure using the extended finite element method. Petroleum Science, 2015, 12(2): 304–315

[96] Zhuang X, Chun J, Zhu H. A comparative study on unfilled and filled crack propagation for rock-like brittle material. Theoretical and Applied Fracture Mechanics, 2014, 72: 110–120

[97] Réthoré J, Borst R, Abellan M A. A two-scale approach for fluid flow in fractured porous media. International Journal for Numerical Methods in Engineering, 2007, 71(7): 780–800

[98] Song C, Wolf J P. The scaled boundary finite-element methodalias consistent infinitesimal finite-element cell method for elastodynamics. Computer Methods in Applied Mechanics and Engineering, 1997, 147(3–4): 329–355

[99] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[100] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92: 242–256

[101] Hughes T J, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, nurbs, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39): 4135–4195

[102] Ghorashi S S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[103] Nguyen-Thanh N, Valizadeh N, Nguyen M, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff–love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[104] Shi G H. Manifold method of material analysis. Tech. rep., DTIC Document, 1992

[105] Shi G H. Modeling rock joints and blocks by manifold method. In: The 33th US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association, 1992

[106] Ma G, An X, He L. The numerical manifold method: a review. International Journal of Computational Methods, 2010, 7(1): 1–32

[107] Li S, Cheng Y, Wu Y F. Numerical manifold method based on the method of weighted residuals. Computational Mechanics, 2005, 35(6): 470–480

[108] Lin J S. A mesh-based partition of unity method for discontinuity modeling. Computer Methods in Applied Mechanics and Engineering, 2003, 192(11): 1515–1532

[109] Terada K, Asai M, Yamagishi M. Finite cover method for linear and non-linear analyses of heterogeneous solids. International Journal for Numerical Methods in Engineering, 2003, 58(9): 1321–1346

[110] Terada K, Kurumatani M. Performance assessment of generalized elements in the finite cover method. Finite Elements in Analysis and Design, 2004, 41(2): 111–132

[111] Terada K, Ishii T, Kyoya T, Kishino Y. Finite cover method for progressive failure with cohesive zone fracture in heterogeneous solids and structures. Computational Mechanics, 2007, 39(2): 191–210

[112] Zheng W, Zhuang X, Tannant D D, Cai Y, Nunoo S. Unified continuum/discontinuum modeling framework for slope stability assessment. Engineering Geology, 2014, 179: 90–101

[113] Kurumatani M, Terada K. Finite cover method with multi-cover layers for the analysis of evolving discontinuities in heterogeneous media. International Journal for Numerical Methods in Engineering, 2009, 79(1): 1–24

[114] Gao H, Cheng Y. A complex variable meshless manifold method for fracture problems. International Journal of Computational Methods, 2010, 7(1): 55–81

[115] Zhang H, Li L, An X, Ma G. Numerical analysis of 2D crack propagation problems using the numerical manifold method. Engineering Analysis with Boundary Elements, 2010, 34(1): 41–50

[116] Chen G, Ohnishi Y, Ito T. Development of high-order manifold method. International Journal for Numerical Methods in Engineering, 1998, 43(4): 685–712

[117] Jiang Q, Zhou C, Li D. A three-dimensional numerical manifold method based on tetrahedral meshes. Computers & Structures, 2009, 87(13): 880–889

[118] Belytschko T, Lu Y Y, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[119] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically nonlinear three-dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[120] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[121] Rabczuk T, Areias P, Belytschko T. A meshfree thin shell method for nonlinear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[122] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367–382

[123] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[124] Rabczuk T, Belytschko T. Application of particle methods to static fracture of reinforced concrete structures. International Journal of Fracture, 2006, 137(1–4): 19–49

[125] Rabczuk T, Areias P. A new approach for modelling slip lines in geological materials with cohesive models. International Journal for Numerical and Analytical Methods in Geomechanics, 2006, 30(11): 1159–1172

[126] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23): 1391–1411

[127] Mossaiby F, Bazrpach M, Shojaei A. Extending the method of exponential basis functions to problems with singularities. Engineering Computations, 2015, 32(2): 406–423

[128] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37): 2437–2455

[129] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[130] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29): 2777–2799

[131] Ai W, Augarde C E. An adaptive cracking particle method for 2D crack propagation. International Journal for Numerical Methods in Engineering, 2016, 108(13): 1626–1648

[132] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid–structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[133] Zhu H, Zhuang X, Cai Y, Ma G. High rock slope stability analysis using the enriched meshless shepard and least squares method. International Journal of Computational Methods, 2011, 8(2): 209–228

[134] Zhuang X, Augarde C, Mathisen K. Fracture modeling using meshless methods and level sets in 3D: framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[135] Zhuang X, Huang F, Zhu H. Modelling 2D joint propagation in rock using the meshless methods and level sets. Chinese Journal of Rock Mechanics and Engineering, 2012, 31: 21872196 (in Chinese)

[136] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[137] Silling S A, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. Journal of Elasticity, 2007, 88(2): 151–184

[138] Ganzenmüller G C, Hiermaier S, May M. On the similarity of meshless discretizations of peridynamics and smooth-particle hydrodynamics. Computers & Structures, 2015, 150: 71–78

[139] Bessa M, Foster J, Belytschko T, Liu W K. A meshfree unification: reproducing kernel peridynamics. Computational Mechanics, 2014, 53(6): 1251–1264

[140] Shojaei A, Mudric T, Zaccariotto M, Galvanetto U. A coupled meshless finite point/peridynamic method for 2D dynamic fracture analysis. International Journal of Mechanical Sciences, 2016, 119: 419–431

[141] Bobaru F, Yang M, Alves L F, Silling S A, Askari E, Xu J. Convergence, adaptive refinement, and scaling in 1D peridynamics. International Journal for Numerical Methods in Engineering, 2009, 77(6): 852–877

[142] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 108(12): 1451–1476

[143] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[144] Ha Y D, Lee J, Hong J W. Fracturing patterns of rock-like materials in compression captured with peridynamics. Engineering Fracture Mechanics, 2015, 144: 176–193

[145] Wang Y, Zhou X, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics. Engineering Fracture Mechanics, 2016, 163: 248–273

[146] Zhou X P, Wang Y T. Numerical simulation of crack propagation and coalescence in pre-cracked rock-like Brazilian disks using the non-ordinary state-based peridynamics. International Journal of Rock Mechanics and Mining Sciences, 2016, 89: 235–249

[147] Ren H, Zhuang X, Rabczuk T. A new peridynamic formulation with shear deformation for elastic solid. Journal of Micromechanics and Molecular Physics, 2016, 1(02): 1650009

[148] Mossaiby F, Shojaei A, Zaccariotto M, Galvanetto U. OpenCL implementation of a high performance 3D peridynamic model on graphics accelerators. Computers & Mathematics with Applications, 2017, 1856–1870

[149] Brebbia C A, Walker S. Boundary Element techniques in Engineering. Elsevier, 1980

[150] Mi Y, Aliabadi M. Three-dimensional crack growth simulation using BEM. Computers & Structures, 1994, 52(5): 871–878

[151] Simpson R N, Bordas S P, Trevelyan J, Rabczuk T. A two-dimensional isogeometric boundary element method for elastostatic analysis. Computer Methods in Applied Mechanics and Engineering, 2012, 209: 87–100

[152] Lafhaj Z, Shahrour I. Use of the boundary element method for the analysis of permeability tests in boreholes. Engineering Analysis with Boundary Elements, 2000, 24(9): 695–698

[153] Nguyen B, Tran H, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 252–275

[154] Cundall P A. A computer model for simulating progressive large-scale movements in blocky rock systems. In: Procedings of the Symposio of the International Society of Rock Mechanics, Nancy, 1971

[155] Cundall P A. Rational design of tunnel supports: a computer model for rock mass behavior using interactive graphics for the input and output of geometrical data. Tech. rep., DTIC Document, 1974, 1–195

[156] Cundall P A, Strack O D. A discrete numerical model for granular assemblies. Geotechnique, 1979, 29(1): 47–65

[157] Cundall P A. Formulation of a three-dimensional distinct element model part I. A scheme to detect and represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1988, 25: 107–116

[158] Zhu H, Wu W, Zhuang X, Cai Y, Rabczuk T. Method for estimating normal contact parameters in collision modeling using discontinuous deformation analysis. International Journal of Geomechanics, 2016, 17(5): E4016011

[159] Wriggers P. Computational Contact Mechanics. Springer Science & Business Media, 2006

[160] Cundall P A, Hart R D. Numerical modelling of discontinue. Engineering Computations, 1992, 9(2): 101–113

[161] Curran J H, Ofoegbu G I. Modeling discontinuities in numerical analysis. Comprehensive Rock Engineering, 1993, 1: 443–468

[162] Walton O R. Force models for particle-dynamics simulations of granular materials. In: Mobile Particulate Systems. Springer Netherlands, 1995, 287: 367–380

[163] Luding S. About contact force-laws for cohesive frictional materials in 2D and 3D. In: Procedings of Behavior of granular media, 2006, 9: 137–147

[164] Jing L. A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering. International Journal of Rock Mechanics and Mining Sciences, 2003, 40(3): 283–353

[165] Huang D, Wang J, Liu S. A comprehensive study on the smooth joint model in DEM simulation of jointed rock masses. Granular Matter, 2015, 17(6): 775–791

[166] Lorig L. A simple numerical representation of fully bonded passive rock reinforcement for hard rocks. Computers and Geotechnics, 1985, 1(2): 79–97

[167] Kochen R, Andrade J C O. Predicted behavior of a subway station in weathered rock. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(3): 160–e1–160.e13

[168] Souley M, Hoxha D, Homand F. Distinct element modelling of an underground excavation using a continuum damage model. International Journal of Rock Mechanics and Mining Sciences, 1999, 35(4–5): 442–443

[169] Rawlings C, Barton N, Bandis S, Addis M, Gutierrez M. Laboratory and numerical discontinuum modeling of wellbore stability. Journal of Petroleum Technology, 1993, 45(11): 1086–1092

[170] Gutierrez M, Makurat A. Coupled HTM modelling of cold water injection in fractured hydro-carbon reservoirs. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(3): 113.e1–113.e15

[171] Jing L. Numerical modelling of jointed rock masses by distinct element method for two- and three-dimensional problems. Dissertation for PhD degree. Luleå University of Technology, Sweden

[172] Harper T, Last N. Response of fractured rock subject to fluid injection part III. Practical application. Tectonophysics, 1990, 172(1–2): 53–65

[173] Shi G H. Stereographic method for the stability analysis of the discontinuous rocks. Scientia Sinica, 1977, 3: 260–271

[174] Warburton P M. Some modern developments in block theory for rock engineering. Analysis and Design Methods: Comprehensive Rock Engineering: Principles, Practice and Projects 2, 2013: 293–315

[175] Goodman R E, Shi G H. Block Theory and Its Application to Rock Engineering. Prentice-Hall Englewood Cliffs, NJ, 1985

[176] Shi G H, Goodman R E. The key blocks of unrolled joint traces in developed maps of tunnel walls. International Journal for Numerical and Analytical Methods in Geomechanics, 1989, 13(2): 131–158

[177] Karaca M, Goodman R. The influence of water on the behaviour of a key block. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1993, 30: 1575–1578

[178] Jakubowski J, Tajdus A. The 3D Monte Carlo simulation of rigid block around a tunnel. In: Mechanics of Jointed and Faulted Rock. Rotterdam: Balkema, 1995, 551–6

[179] Kuszmaul J, Goodman R. An analytical model for estimating key block sizes in excavations in jointed rock masses. In: Fractured and Jointed Rock Masses. Rotterdam: Balkema, 1995,19–26

[180] Windsor C R. Block stability in jointed rock masses. In: Nedlands W A, ed. CSIRO Rock Mechanics Research Centre. Fractured and Jointed Rock Masses, Lake Tahoe, California. 1992, 65–72

[181] Mauldon M, Chou K, Wu Y. Linear programming analysis of key-block stability. In: Computer Methods and Advancements in Geomechanics, 1997, 1: 517–22

[182] Wibowo J L. Consideration of secondary blocks in key-block analysis. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(3): 333.e1–333.e2

[183] Song J J, Lee C I, Seto M. Stability analysis of rock blocks around a tunnel using a statistical joint modeling technique. Tunnelling and Underground Space Technology, 2001, 16(4): 341–351

[184] Lee I M, Park J K. Stability analysis of tunnel key-block: a case study. Tunnelling and Underground Space Technology, 2000, 15(4): 453–462

[185] Warburton P. Vector stability analysis of an arbitrary polyhedral rock block with any number of free faces. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1981, 18: 415–427

[186] Shi G H, Goodman R E. Two-dimensional discontinuous deformation analysis. International Journal for Numerical and Analytical Methods in Geomechanics, 1985, 9(6): 541–556

[187] Shi G H. Three-dimensional discontinuous deformation analyses. In: DC Rocks 2001, The 38th US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association, 2001

[188] Zhang X, Lu M. Block-interfaces model for non-linear numerical simulations of rock structures. International Journal of Rock Mechanics and Mining Sciences, 1998, 35(7): 983–990

[189] Shyu K. Nodal-based discontinuous deformation analysis. Dissertation for PhD degree. University of California, Berkeley, 1993

[190] Jing L. Formulation of discontinuous deformation analysis (DDA) an implicit discrete element model for block systems. Engineering Geology, 1998, 49(3): 371–381

[191] Lin C T, Amadei B, Jung J, Dwyer J. Extensions of discontinuous deformation analysis for jointed rock masses. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1996, 33: 671–694

[192] Kim Y I, Amadei B, Pan E. Modeling the effect of water, excavation sequence and rock reinforcement with discontinuous deformation analysis. International Journal of Rock Mechanics and Mining Sciences, 1999, 36(7): 949–970

[193] Jiang Q, Yeung M. A model of point-to-face contact for three-dimensional discontinuous deformation analysis. Rock Mechanics and Rock Engineering, 2004, 37(2): 95–116

[194] Hsiung S M. Discontinuous deformation analysis (DDA) with nth order polynomial displacement functions. In: DC Rocks 2001, The 38th US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association, 2001

[195] Koo C, Chern J. The development of DDA with third order displacement function. In: Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, Berkeley, CA. TSI Press: Albuquerque, 1996, 12–14

[196] Tonon F. Analysis of single rock blocks for general failure modes under conservative and non-conservative forces. International Journal for Numerical and Analytical Methods in Geomechanics, 2007, 31(14): 1567–1608

[197] Tonon F, Asadollahi P. Validation of general single rock block stability analysis (bs3d) for wedge failure. International Journal of Rock Mechanics and Mining Sciences, 2008, 45(4): 627–637

[198] Nezami E G, Hashash Y M, Zhao D, Ghaboussi J. A fast contact detection algorithm for 3-D discrete element method. Computers and Geotechnics, 2004, 31(7): 575–587

[199] Shi G H. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of deformable block structures. Engineering Computations, 1992, 9(2): 157–168

[200] Wu W, Zhu H, Zhuang X, Ma G, Cai Y. A multi-shell cover algorithm for contact detection in the three-dimensional discontinuous deformation analysis. Theoretical and Applied Fracture Mechanics, 2014, 72: 136–149

[201] Li H, Bai Y, Xia M, Ke F, Yin X. Damage localization as a possible precursor of earthquake rupture. Pure and Applied Geophysics, 2000, 157: 1929–1943

[202] Mühlhaus H, Sakaguchi H, Wei Y. Particle based modelling of dynamic fracture in jointed rock. In: Proceedings of the 9th international conference of the international association of computer methods and advances in geomechanics–IACMAG, 1997, 97: 207–216

[203] Napier J, Dede T. A comparison between random mesh schemes and explicit growth rules for rock fracture simulation. International Journal of Rock Mechanics and Mining Sciences, 1997, 34(3): 217.e1–217.e3

[204] Place D, Mora P. Numerical simulation of localisation phenomena in a fault zone. Pure and Applied Geophysics, 2000, 157, 11–12: 1821–1845

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature