## 7 Discrete Numerical Methods

In this section we will present a brief description of different element-based numerical methods that can be used in the modelling of fracking problems. The boundary element method (BEM) has been used in brittle anisotropic problems including crack propagation. The extended finite element method (X-FEM) has been developed recently and is also a good choice for fracture mechanics problems, and can be easily applied in cohesive models. Meshless methods are becoming popular in fracture mechanics problems. The discrete element method (DEM) is particularly used in problems with rock materials. The phase-field method and the configurational force method are also reviewed in this section.

#### Authors:

Gabriel Hattori^{1}, Jon Trevelyan^{1}, Charles E. Augarde^{1}, William M. Coombs^{1}, Andrew C. Aplin^{2}

^{1}School of Engineering and Computing Sciences, Durham University, South Road, Durham DH1 3LE, UK. ^{2}Science Labs, Department of Earth Sciences, Durham University, Durham DH1 3LE, UK

Received: 2 September 2015, Accepted: 6 January 2016, Published online: 29 January 2016

©The Author(s) 2016.

### 7.1 Boundary Element Method (BEM)

The boundary element method has first appeared in the work of Cruse and Rizzo [52] for elasticity problems, but it was effectively named as BEM in the work of Brebbia and Domínguez [41] and represented a series of advances in comparison to the existent domain discretisation methods as the finite element method (FEM) and the finite differences method (FDM) [109]:

- Accurate mathematical representations of the underlying physics are employed, resulting in the ability of the BEM to provide highly accurate solutions;
- The problem is defined only at the boundaries, which gives a reduction of dimensionality in the mesh (linear for 2D problems and surface for 3D problems), therefore resulting in a reduced set of linear equations to be solved;
- In spite of the boundary-only meshing, results at any internal point in the domain can be calculated once the boundary problem has been solved;
- There is a great advantage in certain classes of problem that can be characterised by either (1) infinite (or semi-infinite) domains, or (2) discontinuous solution spaces. These advantages have resulted in the BEM gaining popularity for acoustic scattering, fracture mechanics, re-entry corners and other stress intensity problems, where domain discretisation methods have poorer convergence.

However, there are some drawbacks which may deterred FEM users from migrating to the BEM:

- The system of equations is both non-symmetric and fully populated, which may lead to longer computing times (compared to FEM for example), especially in 3D problems. In this case, techniques such as the fast multipole method [228] have been introduced to accelerate the solution in large-scale problems;
- A Fundamental Solution (FS) or Green’s function, describing the behaviour of a point load in an infinite medium of the material properties is required as part of the kernel of the method. This can make the use of BEM infeasible in problems where a FS is not available;
- calculation of the FS must be computationally efficient, which makes explicit FS formulations very desirable in this sense. Dynamic problems usually have implicit formulations, see [60, 227, 277] for instance, where the FS is expressed in a integral form by means of the Radon transform;
- The BEM formulation requires the evaluation of weakly singular, strongly singular and sometimes hypersingular integrals which must be carefully treated. This can be done through a variety of methods, including singularity substraction, e.g. [100], or analytical regularisation, e.g. [91];
- Non-linear problems (e.g. material non-linearities) are difficult to model;

The constitutive equations are given as

with C_{ijkl} and σ_{ij} denoting the elastic stiffness and the mechanical stresses, respectively, andwhere u_{i}

are the elastic displacements. The Einstein summation notation applies in Eqs. (24) and (25).

The elastic tractions p_{ij} are given by

with n=(n_{1},n_{2},n_{3}) being the outward unit normal to the boundary.

The time-harmonic equilibrium equations in the absence of body forces can be written as

where t is the time and ρ is the mass density of the material.

From Fig. 6, let Ω be a cracked domain with boundary Γ, which can be decomposed into two boundaries, an external boundary Γc and an internal crack Γ_{crack}=Γ+∪Γ− represented by two geometrically coincident crack surfaces.

**Fig. 6** Elastic body with a crack

The Dual BEM formulation for time-harmonic loading relies on two boundary integral equations (BIEs), one with respect to the displacements at a point ξ of the domain Ω

and a BIE with respect to the generalised tractions

which follows from the differentiation of the displacement BIE and further substitution into the constitutive laws equation (for details see [90]). N_{r} stands for the outward unit normal to the boundary at the collocation point ξ,_{cij} is the free term that comes from the Cauchy Principal Value integration of the strongly singular kernels p^{∗}_{ij},u^{∗}_{ij} and p^{∗}_{ij} are the displacement and traction FS and d^{∗}r_{ij} and s^{∗}r_{ij} follow from derivation and substitution into Hooke’s law of u^{∗}_{ij} and p^{∗}_{ij}, respectively.

In most cases, the cracks are considered to be free of mechanical tractions. These boundaries conditions can be summarised as

where the ‘+’ and ‘−’ superscripts represents the upper and lower crack surfaces, respectively. Eqs. (28) and (29) can be redefined in terms of the crack tip opening displacement (Δu_{J}=u^{+}J−u^{−}J) in function of the crack-free boundary Γ_{c} and one of the crack surfaces, say Γ_{+}

In this latter equation, the free term has been set to unity due to the additional singularity arising from the coincidence of the two crack surfaces. The inconvenience of this approach is that the BEM formulation will now involve integrals including both strong singularities which require special treatment. Numerous hypersingular approaches have been developed, in particular to anisotropic materials under static [90, 91, 150, 282] and time-harmonic [6, 93, 94, 226, 232, 283, 294] loadings.

The use of a hypersingular formulation does not limit at all the crack shape, being valid for curved and branched cracks, for example. However, it is commonplace to make use of discontinuous boundary elements to ensure that all collocation points lie on the smooth surface within the body of an element; this is required to satisfy the Hölder continuity requirement of the hypersingular BIE.

As stated previously, the Stress Intensity Factors (SIF) are the measure of the stress amplification at the crack tip. They are used extensively when estimating the structural life in a number of applications, from civil engineering structures to aerospace devices. Therefore, a precise calculation of this parameter is essential. The principal difficulty, faced throughout the development of BEM and FEM approaches for modelling LEFM problems, is the use of these discrete techniques to capture the singular stress solution. Traditional finite element piecewise polynomial shape functions are ineffective. We now describe some common approaches to obtain the SIFs:

1- Quarter-point: Developed by Henshell and Shaw [119] and Barsoum [19] for finite elements, it consists in moving the mid-side node of a quadratic boundary element from the centre to 1/4 of the element length from the crack tip. It was shown that the mapping between the element in real space and in the space of the intrinsic coordinates automatically captures the asymptotic displacement behaviour of 1/r√ present in the vicinity of the crack tip (refer to [231] for further explanations).

2- J-integral: Proposed by Rice [224], a path independent integral (assuming a non-curved crack) is used to evaluate the energy release rate due to the presence of the crack,

where n_{1} is the component of the outward unit normal vector in the x_{1} direction, u_{i} are the displacement and t_{i} are the tractions. The term W=1/2σ_{ij}ε_{ij }is the strain energy density.

3- Interaction integral: the J-integral can be decomposed into 3 parts [110, 176]

^{(1)}is the J-integral of the so-called principal state, which represents the energy release rate of the material; J

^{(2)}is the J-integral of the auxiliary state, which depends on the displacements around the crack tip; M

^{(1,2)}is the interaction integral containing terms of the principal and auxiliary state, and is defined as

where *A* is the area inside the contour Γ*j* surrounding the crack tip, and *W*^{(}^{1}^{,}^{2}^{)} is given as

Let us remark that the indices (1) and (2) correspond to the principal and auxiliary states, respectively.

The quarter-point approach allows a direct extrapolation of the SIF by using the crack opening displacement. The J-integral is more cumbersome numerically since the displacements and tractions at the closed path integral are part of the BEM domain and have to be evaluated first; however it is more accurate than the direct extrapolation.

Chen [47] has analysed mixed mode SIFs of anisotropic cracks in rocks with a definition of the J-integral for anisotropic materials and the relative displacements at the crack tip. In Ke et al. [133], the authors have suggested a methodology to obtain the fracture toughness of anisotropic rocks through experimental measurements of the elastic parameters and further comparison with a BEM code. In another work, Ke et al. [132] have proposed a crack propagation model for transversely isotropic rocks. Let us remark that all the previously mentioned works have used the Lekhnitskii formalism [145] in order to model the anisotropy of the material. The Lekhnitskii formalism is a polynomial analogy form of the matricial Stroh formalism.

Crack propagation problems have also been studied under the BEM framework. Portela et al. [214] used the maximum stress criterion as crack growth criteria in a dual BEM. Quasi-static 3D crack growth is analysed in [169].

Cohesive models have also been developed with the BEM: Oliveira and Leonel [195, 196] have proposed a cohesive crack growth model, where the zone ahead of the crack tip is modelled as a fictitious crack model. This formulation gives rise to a volume integral, which must be regularised. The cohesive stresses are dependent on the crack tip opening displacement.

Yang and Ravi-Chandar [286] have proposed a cohesive model where the single-domain dual integral equations are used as an artifice to avoid the mathematical degeneration of the formulation imposed by the crack. In this case, the domain is divided in two sub-domains, where the crack is in the fictional domain division. Moreover, the cohesive zone is modelled as an elastic spring connecting both crack faces. Normal and tangential crack tip opening displacements are considered, and the crack growth is obtained from successive iterations of the non-linear system of equations, where the stiffness of the cohesive zone is taken into account.

Saleh and Aliabadi [233, 234, 235] and Aliabadi et al. [7] have studied the crack propagation problem in concrete using a fictitious crack tip zone. The cohesive zone is modelled with additional boundary elements at the fictitious crack tip that satisfy a softening cohesive law. A major drawback of this methodology is that the crack growth path has to be known a priori.

#### 7.1.1 Fast Multipole Method (FMM)

The linear system formed in the BEM framework is much smaller than its equivalent with FEM formulation. However, the resulting matrix is full, not sparse like the FEM stiffness matrix, and this considerably increases the computational time required to solve a large problems. In 1985, Rokhlin [228] developed a method to reduce the complexity of solving the system of equations to O(n) instead of O(n^{3}), where n is the number of unknowns. This technique was named the Fast Multipole Method (FMM), and generally involves using an iterative solver (such as GMRES [230]) to solve the linear system

which comes from the discretisation of Eqs. (31) and (32).

The Green’s functions in the BIEs can be expanded as follows

where **x**_{e} is an expansion point near **x** obtained through Taylor series expansion, for instance. The original integral containing u^{∗}_{ij} can be rewritten as

where Γ_{a} is a boundary away from ξ. This change allows the collocation point ξ to be independent of the observation point x due to the introduction of a new point **x**_{e}. Equation (39) has to be evaluated only once for different collocation points.

The FMM applied in BEM can be described by the following steps [150]:

- Discretise the boundary Γ;
- Determine a tree structure of the elements. For example, in a 2D domain, define a square containing the entire boundary and call this square the cell of level 0. Then, divide the square into 4 equal cells and call them level 1. Repeat until each cell contains a predetermined number of elements (in Fig. 7, each cell has one element). Cells with no children cell are called leafs. For 3D cases, the same principle applies using cubic cells instead of square cells;
- Compute the moments on all cells for all levels l≥2 and trace the tree structure (shown in Fig. 8). The moment is the term from Eq. (39) that is independent from the collocation point. The moment of parent cells is calculated from the summation of the moments of its 4 children cells;
- Compute the local expansion coefficients on all cells starting from level 2 and tracing the tree structure downward to all leaves. The local expansion of the cell C is the sum of the contributions from the cell in the interaction list of the cell and the far cells. The interaction list is composed by all the cells from the level l that do not share any common vertices with other cells at the same level, but their parent cells do share at least one common vertex at level l−1. Cells are said to be far cells of C if their parent cells are not adjacent to the parent cell of C;
- Compute the integrals from element in leaf cell C and its adjacent cells as in standard BEM. The cells in the interaction list and the far cells are calculated using the local expansion; Obtain the solution of
**Ax**=**b**. The iterative solver updates the unknown solution of**x**and goes to step 3 to evaluate the next matrix vector product**Ax**until the solution converges within a given tolerance.

**Fig. 7** Hierarchical tree structure

**Fig. 8** Hierarchical quad-tree structure

The FMM has been used in 3D fracture mechanics problems as can be seen in [192, 290], and some recent works on GPU can be found in [101, 108, 278]. The FMM is largely detailed in [149].

#### 7.1.2 Adaptive Cross Approximation (ACA)

The Adaptive Cross Approximation (ACA) approach uses a different technique in order to reduce the complexity of the BEM with respect to the storage and operations. ACA uses the concept of hierarchical matrices introduced by Hackbusch [107], where a geometrically motivated partitioning into sub-blocks takes place, and each sub-block is classified as either admissible or inadmissible according to the separation of the node clusters within them.

The main idea is that admissible blocks are approximated by low-rank approximants formed as a series of outer products of row and column vectors, greatly accelerating the evaluation of the matrix vector product that lies within each iteration of an iterative solver. While the FMM deals with the analytical decomposition of the integral kernels, ACA can evaluate only some original matrix entries, or use a full pivoted form where all terms of matrix are calculated, to get an almost optimal approximation. The approximation of matrix A∈C^{t×S} is given by

where k is a low-rank compared to t and s. It is important to remark that the low-rank representation can only be found when the generating kernel function in the computational domain of **A** is asymptotically smooth. It has been shown in [20] that elliptic operators with constant coefficients have this property.

In hierarchical matrices, the near and far fields have to be separated. The index sets I for row and J for columns so that elements far away will have indices with a large offset.

By means of a distance based hierarchical subdivision of I and J cluster trees T_{I} and T_{J} are created. In each step of this procedure, a new level of son clusters is inserted into the cluster trees. A son cluster is not further subdivided and is said to be a leaf if its size reaches a prescribed minimal size b_{min}. Usually one of two different approaches is considered. First, a subdivision based on bounding boxes splits the domain into axis-parallel boxes which contain the son clusters. Alternatively, a subdivision based on principal component analysis splits the domain into well-balanced son clusters leading to a minimal cluster tree depth.

Now, the hierarchical (H)-matrix structure is defined by the block cluster tree T_{IJ}=T_{I}×T_{J} using the following admissibility criterion: min(diam(t),diam(s))≤ηdist(t,s), with the clusters t⊂T_{I},s⊂T_{J}, and the admissibility parameter 0<η<1. The diameter of the clusters t and s, and their distance, are obtained as

A block b is said to be admissible if it satisfies this admissibility criterion. Otherwise, the admissibility is recursively verified for each son cluster, until the block becomes admissible or reaches the minimum size.

Finally, the ACA method idea is to split the matrix **A**∈C^{t×s} into **A**=**S**_{k}+**R**_{k}, where **S**_{k} is the rank k approximation and **R**_{k} is the residuum which has to be minimised. We now present the ACA method itself:

1-Define k=0 where **S**_{0}=0 and **R**_{0}=**A** and the first scalar pivot to be found is γ1=(R0)−1_{ij}, and i, j are the row and column indices of the actual approximation step;

2-For each step υ, obtain

where the operators ()i and ()j indicate the i-th row and the j-th column vectors, respectively;

3-The next pivot γ_{υ+1 }is chosen to be the largest entry in modulus of the row (R_{υ})_{i} or the column (R_{υ})_{j}

4-The approximation stops when the following criterion holds

The main advantage comparing to the FMM method is that ACA can be implemented directly into an existing BEM code. Moreover, due to its inherently parallel data structure, parallel programming can be readily implemented, increasing the computational efficiency. However, the original matrix A will not be entirely recovered.

Note that it is not necessary to build the whole matrix beforehand. The respective matrix entries can be computed on demand [20]. Working on the matrix entries has the advantage that the rank of the approximation can be chosen adaptively while kernel approximation requires an a priori choice.

A few recent works on ACA implementation can be found in [81, 99]. Use of the method for problems in 3D elasticity can be found in [28, 158] and the application of ACA in crack problems was shown for the first time in [137].

### 7.2 Enriched Formulations

#### 7.2.1 eXtended Finite Element Method (X-FEM)

The motivation that lay behind the development of X-FEM was to eliminate some of the deficiencies of standard FEM for crack modelling, most importantly the requirement for highly refined meshing around the crack tips and the mandatory remeshing for crack growth problems. The partition of unity [15] is a general approach that allows the enrichment of finite element approximation spaces so that the FEM has better convergence properties. In X-FEM, the partition of unity method allows element enrichment such that degrees of freedom (dofs) are added to represent discontinuous behaviour. In this framework, the mesh is independent from the discontinuities, so that cracks may now pass through elements rather than being constrained to propagate along elment edges. This gives the FEM much more flexibiility to model crack growth without remeshing.

Two types of enrichment function are applied in the X-FEM: the Heaviside enrichment function, responsible for characterising the displacement discontinuity across the crack surfaces, and a set of crack tip enrichment functions (CTEFs), responsible for capturing the displacements asymptotically around the crack tip. This latter presents complex behaviour, varying for different constitutive laws (see [12, 79, 193], for some different CTEF). In this sense, it is similar to the FS, necessary in BEM formulations.

The displacement approximation **u**^{h}(**x**) with the partition of unity can be stated as [176]

where N_{i} is the standard finite element shape function associated with node i,**u**i is the vector of nodal dofs for classical finite elements, and **a**_{j} and **b**^{α}_{k} are the added set of degrees of freedom that are associated with enriched basis functions, associated with the Heaviside function H(x) and the CTEF Fα(x), respectively. N is the set of all nodes, NH is the set of all nodes lying on crack surfaces, and NCT is the set of all nodes belonging to elements touching a crack tip.

Since the CTEFs describe the displacements at the crack tip zone through the addition of several dofs, the stress concentration around the crack tip can be found more accurately with a significantly coarser mesh compared to the mesh used with standard FEM in a similar problem.

The presence of blending elements, which do not contain the crack but contain enriched nodes, is also important and has to be considered. These elements were analysed by Chessa et al. [48], and some studies have improved the convergence of blending elements (see [84], for instance). The X-FEM convergence rate can also be increased through the use of geometrical enrichment [142], where a number of elements around the crack tip receive the CTEF instead of a single element (this latter named topological enrichment).

**Fig. 9** Elastic body with a cohesive crack

Figure 9 illustrates an arbitrary elastic body with a cohesive crack. The governing equations for a cohesive crack model are given by [178]

where **Ω** is the domain, **f**^{b} is the body force vector, **f**^{t} is the external traction vector, σ is the stress tensor, α is the load factor which controls the load increments, **f**^{c} is the traction along the cohesive zone, and is a function of the crack opening Δu.

The discretisation of Eq. (50) yields

where **B** is the finite element strain-displacement matrix, b is the vector containing the body forces and **T**c(Δu) is the cohesive softening law relating the crack surface normal traction **f**_{c} to the crack opening Δu.

X-FEM has been widely used with cohesive models in the last few years. Some authors [45, 51, 175] have used a typical X-FEM formulation to model the cohesive crack, i.e., a Heaviside enrichment function is used to model the jump between the crack surfaces and a crack tip enrichment function is used to model the asymptotic behaviour at the crack tip.

Xiao and Karihaloo [284] have obtained the asymptotic displacement at the cohesive zone for isotropic materials based on the Williams expansion. The authors considered only the case where the crack is traction free and the crack is subject only to mode I. The obtained enrichment functions are

where κ is the Kolosov constant (for details refer to [284]), θ is the crack orientation with respect to the x1 axis, a_{11} is a real constant and comes from the Williams expansion. In this case, Eq. (57) receives a new crack tip enrichment term, as in the X-FEM formulation for linear elastic fracture mechanics (see [110, 175]). Zamani et al. [291] uses higher-order terms of the crack tip asymptotic field to obtain an enrichment function based on the Williams expansion.

This approach has provided good results for isotropic materials, but it may not be the same for anisotropic materials. An alternative approach is to model the crack with Heaviside elements only [139, 180, 262, 302]. Since there is no discontinuity at the crack tip, there are no SIFs at the crack tip, and therefore no crack tip enrichment function is required. The displacement field **u**(**x**) is given by

where N_{i} is the standard finite element shape function associated to node i and **a**_{j} is the additional set of degrees of freedom associated with the Heaviside enrichment function H, defined as +1 if it is evaluated above the crack or −1 if below the crack. The sets N and N^{H} denote the standard and enriched nodes, respectively.

The crack growth is modelled considering some rules, for example, if the level of stress at the crack tip is above the material tensile strength [178, 262].

In [139], a 2D cohesive model for an isotropic material was presented, where both fluid and porous material interact. The pressure inside a crack is also modelled. The Heaviside enrichment function is employed, as well as a pressure enrichment function, which allows the continuity of steep gradients without enforcing this condition. The crack propagation criteria depends on the stress state at the crack tip. The fluid behaviour can retard crack initiation and propagation. A local change of the flow can be seen immediately after crack propagation. The deformation around the crack causes fluid to flow mostly from the crack itself since the crack permeability is much higher than the medium permeability. This flow from the crack to the crack tip causes closing of the crack. However, a delamination test has shown that if the stiffness and permeability are high, the fluid does not influence crack growth.

More methods for crack propagation in X-FEM can be found in [151, 167, 182, 183, 225] for brittle fracture and [168, 179, 291] for cohesive cracks.

#### 7.2.2 Enriched BEM

The extended boundary element method (X-BEM) was first proposed by Simpson and Trevelyan [251] for fracture mechanics problems in isotropic materials. The main idea is to model the asymptotic behaviour of the displacements around the crack tips by introducing new degrees of freedom. The displacements **u**^{h}(**x**) are thus redefined as

where N and N^{CT} are the sets with non-enriched and enriched nodes, respectively, N_{i} is the standard Lagrangian shape function associated with node i,**u**_{i} is the vector of nodal degrees of freedom, and a^{α}_{k} represents the enriched basis functions which capture the asymptotic behaviour around the crack tips. In elastic materials, a^{α}_{k} is an 8-component vector for two-dimensional problems, since only two nodal variables (u_{1}, u_{2}) and four enrichment functions are needed to describe all the possible deformation states in the vicinity of the crack tip [110].

Hattori et al. [110] used the following anisotropic enrichment functions initially developed for the X-FEM

r is the distance between the crack tip and an arbitrary position, θ is the orientation measured from a coordinate system centred at the crack tip, and **A**,**B** and μ are obtained from the following eigenvalue problem

(no sum on m) with

Let us emphasise that the anisotropic enrichment functions can also be used for isotropic materials, since this is a degenerated case from anisotropic materials. For more details please refer to reference [110].