Modelling fracturing process of geomaterial using Lattice Element Method

Fracturing of geomaterial is involved in many geological processes and engineering applications. However, modelling of fracturing process is considered challenging owing to the heterogeneity of geomaterial. In this paper, a simple three-dimensional discontinuum method Lattice Element Method (LEM) is introduced to simulate the fracturing process. Geomaterial is modelled as interconnected 1D spring elements. Fracturing is modelled by simply removing lattice element exceeding a specified threshold related to the critical energy release rate of rock. Mesh dependency phenomena can be manipulated by introducing disorder in model which also incorporates heterogeneity in model. An in-house C++ code using a parallel conjugate gradient solver has been developed which is capable to handle large scale model composed of millions of lattices. Three simulations of fracturing process of a geomaterial with a pre-existing penny shape crack under uniaxial tension are presented. A simple discretisation of domain into 1D springs and the use of efficient solver enable LEM to model the heterogeneity of geomaterial by including large amount of rock features such as faults and joints inferred from geophysical surveys. This can shed light on explaining the complicated fracture patterns observed in brittle geomaterial.


INTRODUCTION
Fracturing of geomaterial is a fundamental geological process governing the formation of tectonic plates on earth crust to fissures in stiff overconsolidated clay. It also relates to earthquakes, volcanic eruptions and groundwater flows in geomaterial. In engineering applications, understanding of how geomaterial fractures is crucial in mining, tunnelling in rock and grouting in soil. Its application can also be extended to exploitation of natural resources such as Enhanced Geothermal System and unconventional resources involving hydraulic fracturing.
Geomaterial is considered as a challenging material for modelling because of its heterogeneity arised from discontinuities in different scales. Modelling fracturing process is another challenge. There are different attempts to simulate fracturing of geomaterial using continuum based or discontinuum based approaches. Finite Element Method (FEM) (Tang 1997), Boundary Element Method (BEM) (Thomas 1993) or more recently Extended Finite Element Method (XFEM) (Gordeliy and Peirce 2013) are continuum approaches in geomaterial fracturing simulation whereas Discontinuous Deformation Analysis (DDA) (Lin et al. 1996) and Discrete Element Method (DEM) / Bonded Particle Model (BPM) (Potyondy and Cundall 2004) are some examples of discontinuum approaches.
A lot of literature studies on fracturing of geomaterial in laboratory scales. However, little literature covers numerical modelling on fracturing process in field scales or on highly fracture rock mass which are relevant to many engineering applications. Modelling discontinuities in different scales explicitly leads to large models for numerical simulation. Also, modelling complex fracture propagation such as non-planar fracture surface and fracture branching and coalescence cannot be modelled easily in continuum based approaches. Therefore, for large scale problems, research mainly focuses on ground characterisation from various monitoring techniques and data analysis using geostatistical methods. This paper presents a simple numerical method, Lattice Element Method (LEM), which is suitable for large scale three-dimensional fracture simulation in geomaterial. It simplifies heterogeneous geomate- rial into 1D lattice elements. Fracture initiation and growth are simply modelled by removing lattices. Complex fracture evolution and coalescence can be simulation in a stable manner.

LATTICE ELEMENT METHOD
LEM is also referred as lattice models or Lattice Spring Model in literature. LEM is often used to simulate fracturing process of heterogeneous medium in micro-or meso-scales. Concrete (Schlangen and Garboczi 1997), cemented granular material (Topin et al. 2007), cellular material (Wang and Stronge 1999), composite material (Snyder et al. 1992) and biomaterial (Hansen et al. 1996) are some of the application areas of LEM. Recently, Zhao et al. (2012) applied LEM on simulation of rock fracturing. To generate a LEM model, the domain is first discretised into nodes. A lattice network is obtained by connecting nodes by 1D lattice elements. A 3D problem is simplified into a 1D lattice network in 3D space. Each node represents a sub-domain called cell. Forces between cells can only be transmitted through lattices, which are the basic elements in LEM. It models the interaction between two cells on their shared surfaces. Therefore, each surface represents surface between two adjacent cells (nodes).
There are only two ways of specifying mechanical properties of a material, by specifying lattice network geometry and by specifying constitutive relationship of lattice. For regular lattice network, the Poission's ratio υ of lattice network depends on the choice of lattice element. If Hookeans spring is used, the Possion's ratio is cannot be chosen freely. However, Possion's ratio can be specified in LEM by using other types of lattice element.
One of the major advantages of LEM is the ease of fracture simulation, even for complicated fracture network. Fracturing is modelled by simply removing or degrading lattices that meets the failure criteria. Fractures are possible along surface between two cells. This is an analogy of preferred fracture path along existing discontinuities of geomaterial.

METHODOLOGY
An in-house C++ code LEM3D is under development. The code can generate large-scale 3D disordered lattice network and can simulate fracturing process according to a specified failure criteria. A highly efficient parallel solver is developed to handle millions of DOFs so fracture simulation involving thousands of steps can be finished within a reasonable time. One of the important issues in generation of lattice models is mesh dependency. For a regular lattice where only several fracture orientations are allowed, fracture path is biased along one of the available orientations (see Figure 2). Such bias can be eliminated by introducing disorder in the mesh. Local heterogeneity of rock is modelled at the same time. Anisotropic material can be modelled by favouring or suppressing certain lattice orientation.

Disordered Network
In this paper, an isotropic geomaterial is studied in which lattice orientation is uniformly distributed such that fracture path is not biased in any orientations. To achieve this in LEM3D, the domain is filled by nodes with randomly generated coordinates. Then, the subdomain (cell) represented by each node is determined by Voronoi tessellation. A lattice network is formed by searching for two nodes whose cells share a common surface and connect them by lattice.
There are two additional criteria imposed in generating disordered network. The first criterion is im-posed on generation of nodes to avoid generating lattice of very small length. The separation between any two nodes cannot smaller than a specified value l s,min . The second criterion is imposed on generation of lattice to avoid lattice representing a small surface area and formation of tiny fracture surface much smaller than the resolution we required. This is done by rejecting lattices whose area is less than a specified value A s,min .

Spring stiffness
The simplest elastic-brittle Hookean's spring is chosen as lattice element in this paper. Only two parameters, spring constant k s and failure force F st are required to define a lattice. The scaling rule of k s is provided by where A s is the surface area represented by lattice, l s is the lattice length and E s is a proportionality constant related to Young's Modulus of geomaterial.
It should be noted that E s is not the macroscopic Young's modulus and needs to be calibrated prior to simulation.

Heterogeneity
The heterogeneity of rock is automatically included in LEM model when disorder is introduced to avoid mesh dependency issue in fracturing simulation. There are several input parameters of LEM controlling the heterogeneity of geomaterial. The number of node N num indicates the heterogeneity of geomaterial. For the same domain, fewer nodes means higher degree of heterogeneity because of wider variation of lattice length and fracture area. The degree of heterogeneity can also be increased by specifying smaller l s,min and A s,min .

Failure Criterion in relation to fracture mechanics
For elastic-brittle lattice, it breaks when its spring force F s exceeding a threshold F st . Such lattice is removed from calculation in subsequence load step. The determination of such threshold can be related to Linear Elastic Fracture Mechanics (LEFM). According to Griffith (1921), total energy of the material-crack system should be unchanged when fracture propagates. Energy is required to extend a crack which is proportional to new fracture surface created. Such energy is provided by release of strain energy of material and is stored in surface of crack called surface energy. For an elastic-brittle spring, all the strain energy is released when it breaks, hence Left hand side of (2) is strain energy stored in lattice under elongation e. G Ic is called critical energy release rate and the subscript I denotes model I (tensile) failure. From constitutive relationship of Hookean's spring F s = k s e, the lattice force threshold F st is expressed as Putting (1) into (3), Hence, F st is proportional to surface area A s and inversely proportional to square root of lattice length l s . Denote microscopic tensile strength to be σ ts = F st /A s , (4) becomes The length scale l s is now expressed into two material parameters, tensile strength σ ts and critical energy release rate G Ic . In other words, the lattice length has to be chosen to match both parameter for a given material.

Adaptive load step
In each load step, lattice forces are calculated to check whether failure criterion is met. In theory, load step should be chosen small enough such that only one lattice fails within one load step. However, this requires a large number of load step unless a small model is used. To reduce computation time, load step is chosen large enough to allow multiple lattices to be failed within one load step. However, it should be small enough to avoid too much lattice breakage at one load step which may not capture the progressive damage of geomaterial.
Since the system is linearly elastic given no lattice breakage within a load step, the load step can be adaptively determined after calculation of each load step rather than adopting a small constant load increment. First, a small initial load F 0 is applied on the model. After lattice force calculation, the most critical lattice is identified which has the highest load to capacity ratio p max where F s and F st are lattice axial force and lattice failure force respectively. p max can also be the average of most critical n max lattices, p n,max . The new force F i+1 depends on the value p n,max as follow  When there is no lattice breakage (i.e. p n,max < 1.0), loading increases and the parameter α > 1.0 controls number of lattice to be broken in next load step. For load step that involves lattice removal (i.e. p n,max ≥ 1.0), the applied force decreases to capture the post peak behaviour and to avoid too much lattice breakage within a load step. The parameter β < 1.0 indicates the possibility of lattice breakage in next load step.

Computation
Given the large DOFs and load steps involved, an efficient linear algebra solver is required to perform simulation in reasonable time. In LEM3D, a solver using Preconditioned Conjugate Gradient (PCG) method is developed. The major merit of PCG is the significant reduction of storage as only several vectors are required to be stored. To speed up using multi-core CPUs, the solver uses OpenMP directives for parallel computing. A LEM model with 2 millions DOFs and 4 millions lattices only takes less than half a minute in solving one load step using 8-core Intel Xeron E5-2670 (2.6GHz) CPU. The solver will be modified for GPU computation for even greater speed up.

SIMULATIONS
Three simulations have been carried out by LEM3D. The dimensions of model are all 100mx100mx100m. All the samples are subjected to uniaxial tension along z-axis. The models are only constrained along z-direction on bottom face. All other faces are free. Three different models (Spare, Medium and Dense models) of isotropic rock with a hypothetical penny shape crack were performed using LEM3D. All three models use same set of input parameters as listed on Table 1 except the number of nodes N n in the model. A 50m diameter penny shape crack perpendicular to z-axis is placed at the centre (0,0,50) of model.

From diffusive breaking to connected fractures
The evolution of fractures in three LEM models are illustrated in Figure 3. At the beginning of simulation, all three models show diffusive breaking of lattice clustering in the vicinity of the edge of penny shape crack. The most critical lattice may not be exactly adjacent to main crack. Instead, the most critical lattice are the result of the combination of factors including lattice orientation, lattice parameters k s and F st , connectivity in addition to its vicinity to main crack. The variation of first three factors contribute the local heterogeneity and induces local variation of lattice forces.
At the beginning, broken lattices are largely unconnected and the geomaterial behaves roughly elastic macroscopically. This is because broken lattices are the weaker ones and their loading can be shared to adjacent lattices which are stronger. Formation of fractures is under controls because the stronger ones have enough spare capacity to take up additional loading without breaking.
The diffusive breaking of lattice can be regarded as fracture process zone ahead of crack tip. Before the crack is extended, the process zone is softened by the formation of microcracks.
When loading increases and more fractures formed, fractures start to be connected and joins the main crack. Breaking and coalescence of fractures lead to overstress of adjacent lattices and further breaking of lattice. Fracture growth starts to become unstable.

Spatial distribution of failure lattice
Initially, failure lattices distribute evenly around main penny shape crack. The density of fracture decreases rapidly away from main crack. Fractures grow in radial direction without apparent bias in certain directions. In Dense model, cracks are more localized along perimeter and failure lattice tends to diffusive in Sparse model.
Afterwards, fracture growth starts to be biased at locations where broken lattices are connected. Such bias is greater in Sparse model and appears at earlier stage.
The orientation of failure surface is shown in Figure 4. The orientation of failure surface deviated from horizontal locally because of local heterogeneity. The overall trend of horizontal fracture grow can be observed. The disorder introduced to LEM mesh successfully removes the mesh dependency observed in regular mesh.

From microscopic to macroscopic
The applied pressure-displacement curves of all three models are plotted in Figure 5. From the plot, brittle failure of rock are observed in all three LEM models. They all shows Class II behaviour (i.e. a snap back of curve after peak) as first observed in laboratory tests of brittle rock by Wawersik and Fairhurst (1970). They explained that the snap back portion of the curve represents unstable fracture propagation of geomaterial, or the fracture growth is 'self-sustaining' without any work done from external load. In other words, after the peak, the strain energy stored in geomaterial is sufficient to continue fracture growth until   collapse. Energy must be extracted from the system to capture the failure process. In LEM3D, this is done by unloading of model whenever breaking of lattice is detected.
In the simulation, the energy release rate G Ic is specified but the microscopic tensile strength σ ts is related to lattice length l s according to (5) where σ ts decrease with l s . From pressure-displacement curves in Figure 5, same trend can also be observed macroscopically.
The macroscopic 'brittleness' can be indicated by the amount of snap-back in the plot. Smallest degree in 'brittleness' is shown in Sparse model which provides highest degree of heterogeneity (due to larger variation of lattice length and area). On the other hand, Dense model provides higher resistance in tension in the expense of higher 'brittleness'.

CONCLUSIONS AND FUTURE WORK
This paper presents the principles of LEM and covers some details in LEM implementation such as generation of disordered model and scaling rule of lattice parameters from geometric information of Voronoi tessellation.
To demonstrate the possibility of LEM to model fracturing of geomaterial considering heterogeneity. Three simulations using an in-house C++ code, LEM3D, are presented. Three LEM models with DOFs up to 2 millions are used to simulate the fracturing process of brittle geomaterial with pre-existing penny shape crack under uniaxial tension. A transition from diffusive and even lattice breaking and localized crack coalescence is observed. The simulations can successfully reproduce the stress-strain curves Figure 4: Lattice network of Sparse Model at load step 2250. Color of lattice shows its force to strength ratio which varies from 0 (blue) to 1 (red). Broken lattices are shown in black. (a) shows a section along penny shape crack and (c) shows a cross section perpendicular to the crack. The isometric view of model is shown in (b). The fracture surface formed is non-planar due to heterogeneity (shown by curved trace of broken lattices along boundaries of model in (b) and two separated trace of lattices at bottom left corner of (a) which are connected underneath the cut plane).
Figure 5: Force-displacement plot of three LEM simulations. All three models demonstrate Class II behavior of brittle rock under laboratory test. The discrete data points corresponding to the snapshot shown in Figure 3 is indicated. observed in laboratory tests.
Although geomaterial is rarely subjected to purely uniaxial tension, the uniaxial compression and other loading conditions can be easily simulated by by specifying the compressive failure criteria of lattices.
More lattice models will be included such as elastic-softening spring to model fracturing process zone. A fluid model will be incorporated in future LEM3D code suitable for simulation of hydraulic fracturing in brittle and fractured geomaterial like shale. With the capacity of LEM to model large scale problem and the advances in parallel computing, LEM has great potential for practical applications by realistic simulation with input from large database from in-situ geophysics measurements such as microseismic monitoring.