Modelling of primary fragmentation in block caving mines using a ﬁnite-element based fracture mechanics approach

The growth of fractures around an undercut of a block cave is simulated. A ﬁnite element based approach is used, in which fractures are represented as non-planar 3D surfaces that grow in response to boundary stresses and interaction. A new mesh is recreated at each step to compute the displacement ﬁeld. Stress intensity factors are computed around fracture tips using a technique that computes the interaction integral over a virtual disk. Fracture geometry is updated using Paris and Scho¨llmann propagation laws, and a geometric fracture pattern ensues from the simulation. The growth of fractures is examined in the lower 20 m of a mine at 812 m depth. The growth of 30, 60 and 90 fractures is examined. Realistic extraction schedules for over 100 draw points control the rate of mass extraction. The effect of rock bridges as overburden stress shields is investigated. Bridges are modelled by constraining the vertical displacement of the top boundary. This case is compared to a Neumann-type overburden stress boundary condition in which the overburden is felt throughout the top of the cave. In both cases, fractures grow to form a dome shape above and around the cave during extraction. For the case of a ﬁxed top boundary, fracture growth is observed away from the cave, while in the direct overburden stress case, fractures tend to grow close to the cave. Over-arching fractures concentric to the undercut continue to grow as the cave progresses.


Introduction
Underground mining methods extract large volumes of rock directly from the subsurface by undercutting rocks that contain minerals of interest, and relying on blasting and gravity-driven fragmentation to extract ore at low cost (Brown 2003). Block caving, in situ leaching, and in situ fragmentation, all aim to decimate the rock whilst still underground so as to reduce the environmental footprint of the mining activity, whilst maintaining low production costs. Fracture growth during tunnelling, caving and blasting is a key controlling parameter in designing extraction procedures, in order to attain specific fragment size distributions and avoid excessive fines production (Rance et al. 2007). Numerical modelling can be useful in determining how fractures grow under specific mine scenarios, and may conduce to large economic benefits to the industry if energy and mass conservative methods are systematically applied as part of subsurface mining procedures. In addition, high-accuracy numerical methods can be instrumental in defining safety limits and identifying mine-specific risk scenarios.
Numerous numerical methods have been developed to model fracture and fragmentation. These are discussed here in the context of subsurface mining, and their suitability to model primary fragmentation. In particular, a main focus of the discussion is finite element derived techniques that model multiple fractures explicitly and account for variable matrix material heterogeneities.
The finite element method (FEM) is suitable for the analysis of fracture interaction during growth, and for combination with flow and poro-elastic analysis for hydraulic fracture growth. FEM requires the generation of meshes to discretise the solution domain. The extended finite element method (XFEM) is a specialisation of FEM, in which elements have enriched shape functions which model a partition of unity function, allowing for discontinuities to exist within each element sub-domain (Moës et al. 1999). The main advantage of XFEM is that it can represent discontinuities within a pre-existing mesh, and can model discrete fracture propagation without re-meshing. The enrichment of the shape functions can be used to model fracture growth and friction, as well as fluid flow through fractured media. However, by avoiding re-meshing, XFEM cannot ensure that the mesh is sufficiently refined ahead of the crack tip as growth progresses, leading to large errors in stress and energy estimations around the tips. It is ideally suited for the problem of moving interfaces, and for fracture analysis cases for which fracture paths can be predicted. Furthermore, enrichment functions become unstable for fracture interaction, and when fractures share nodes.
In contrast, the discrete element method (DEM) models the movement of volumetric objects and their interaction based on the solution of dynamical equations, using a penalty-based damping to resolve collisions. Agglomerates of smaller-scale spherical or polygonal/polyhedral elements represent bodies. Extensions, such as Discontinuous Deformation Analysis (DDA), solve for displacements instead of forces (Shi 1988). DEM is simple to implement, and efficient to run in both 2D and 3D. It can be combined with FEM to model non-linear deformation of the bodies (FEMDEM; Munjiza et al. 1995), usually using Drucker-Prager based plastic or de-bonding pseudofracture models. The standard DEM method utilises arbitrary dampening penalty values that affect the ability of the method to conserve energy. Furthermore, this method requires the definition of micro-properties using laboratory calibration tests that are rock type and scale dependent (Lisjak and Grasselli 2014), and are specific to a fixed mesh size-and therefore cannot be readily translated to the field scale without additional empirical assumptions. De-bonding fracturing is arbitrary and is not based on rigorous fracture growth principles. In contrast, a DEM impulse approach uses directly measurable material properties, such as Young's modulus and Poisson's ratio, and ensures energy conservation during fragmentation by relying on impulses instead of penalties to resolve contact (Tang et al. 2014).
The present paper investigates primary fragmentation around a growing cave, by modelling fracture propagation during caving. As opposed to XFEM or DEM approaches, the current paper directly applies FEM to simulate the growth of multiple fractures. The finite element method is combined with a geometrydriven approach in which fracture geometry evolves independently from its discretisation. This method provides three key ingredients to running multiple fracture growth and interaction simulations: (1) it allows fractures to be represented discretely, as subdimensional entities; (2) it allows for the definition of geologically realistic, site-specific matrix heterogeneities; (3) it uses input properties that are measurable and well-understood rock properties, such as Young's modulus and Poisson's ratio; (4) by applying adaptive re-meshing, the mesh is always suited for the current geometry; (5) it can be directly validated against analytical solutions for stress intensity factors as well as against laboratory experiments for single and multiple fractures under tension, compression and shear loading. Furthermore, this approach represents fracture and matrix domains explicitly, and allows modelling coupled fluid flow, as well as thermal and poro-elastic effects.

Methodology
The method presented here utilises a mesh to compute the displacement field, and replaces that mesh as the geometry progresses Zimmerman 2011, 2013). Thus, the mesh can be substituted by a mesh-free point cloud based computation of the stress field, or by a boundary element based computation of the stress intensity factors. The method aims to remain geometrically focussed, with the sole aim of growing multiple interacting fractures in an accurate manner.

Constitutive model
The rock is assumed to behave in a linear elastic manner, and so the stress-strain relation is given by r ¼ D[e Àê, where e and r are the infinitesimal strain and stress vectors, D is the elasticity matrix, andê is the vector of initial strains. The finite element is applied to solve the deformation field, taking into account surface tractions and body forces.

Geometry
The model geometry, including boundaries and fractures, is defined using three-dimensional cubic Non-Uniform Rational B-Spline (NURBS) surfaces (Fig. 1). The model boundaries form a watertight domain. Multiple fractures grow within the model domain, in a planar or curved manner, depending on the effects of the boundary conditions and the interaction with other fractures. At each step, fracture geometry is updated by extending the NURBS representation of the fracture (Paluszny and Zimmerman 2013).

Discretisation
The volumetric domain is discretised using isoparametric quadratic tetrahedra. Fracture surfaces are discretised using isoparametric quadratic triangles. Side-and corner quarter-point tetrahedral elements are placed along the crack front. These are a type of quadratic elements in which the centre node is shifted towards the tip; they are placed in order to better capture the stress singularity at the fracture tip (Nejati et al. 2015a). Displacements are computed at the nodes, material properties are defined at Gaussian integration points, and stresses and strains are computed at the integration point locations. Tips are discretised into nodes and segments, and stress intensity factors are computed at each node along the tips.

Stress intensity factors and growth criteria
As the rock deforms, energy release rates are computed at the crack tips, after which energy-based methods are used to predict the onset and direction of fracture growth. These serve as local approximations of energy in each modality of growth (I, II, and III) that can be used to determine propagation. Specifically, the interaction integral (Yau et al. 1980;Nakamura and Parks 1988) is computed along discrete locations of the fracture tip in order to approximate model stress intensity factors, K I , K II , and K III (Fig. 2). Classic stress intensity factor computations require structured meshes to be defined around the tips. These are generally created by hand, and can be generated automatically into integrated hybrid meshes (Bremberg and Faleskog 2015). However, as the density of fractures increases, the complexity of the generation of such a structured mesh also increases. A virtual integration technique avoids the need of a structured mesh by using an artificial 'virtual' mesh to integrate around the tips, using a 'virtual' cylinder (Cervenka and Saouma 1997;Paluszny and Zimmerman 2011) or a 'virtual' disk (Nejati et al. 2015b). When aided by isoparametric quadratic elements (Daimon and Okada 2014) and with the quarter-point optimisation (Nejati et al. 2015a), this technique allows swift and accurate computation of stress intensity factors on unstructured tetrahedral meshes. The integration yields high quality solutions when the integration domain radius approximates the mesh size around the fracture (Nejati et al. 2015b). Thus, the use of quarter-point quadratic tetrahedra around the fracture front in combination with the virtual technique allows the use of coarse, unstructured tetrahedral meshes for growth computations. In this work, the disk virtual integration technique is used to compute the interaction integral [see details and equations in Nejati et al. (2015b)], which allows the independent computation of K I , K II , and K III along the fracture tips.
The Paris-law is used to govern fracture extension (Lazarus 2003; Zimmerman 2011) at every tip, and a 3D growth angle models the direction in which each tip extends (Schöllmann et al. 2002). The calculation of local Ks allows the modelling of growth as a function of both interaction and boundary conditions simultaneously (Fig. 3).
The direction of crack growth is not constrained to follow the pre-existing mesh; instead, a new mesh is created at each step to conform to the new fracture geometry. The quality and refinement of the mesh is adjusted to the new geometry. As a result, discrete fractures are allowed to grow, and coalesce, leading to the formation of primary fragmentation.

Implementation
The presented methods have been implemented into the Imperial College Geomechanics Toolkit (ICGT; Paluszny and Zimmerman 2011), which interacts with CSMP?? (Complex Systems Modelling Platform, C??), an object-oriented finite-element based library that is specialised in simulating complex multiphysics processes (Matthäi et al. 2001). It has been extensively validated, against both experimental and Fig. 2 Stress intensity factor J-and interaction integral computation techniques. Initially based on structured meshes, new methods allow for unstructured meshes, whilst reducing integration domain to a virtual disk field experiments, and can run on workstations as well as on high-performance computing systems.

Results and discussion
The growth of fractures embedded in a 1000 9 200 9 200 m cave is modelled. Fractures are initially disk-shaped, and propagate in response to the growth of the mined cave. The models contain between thirty and ninety fractures. The rock is assumed to contain many more fractures, and only the largest of those fractures are modelled here. It is assumed that the radii of the fractures of this scale follow a Gaussian distribution, with a mean of 10 m, and standard deviation of 0.1 m. Initially, fractures do not intersect each other or the boundaries. Growth is expected to be a function of the boundary conditions, and of the effect of the interaction of multiple fractures. A 200 9 200 9 20 m rectangular area undercuts the mine at the centre of the model, and the rock fragments are assumed to migrate through the undercut following a realistic draw point extraction schedule. The mine has 120 draw points that are uniformly distributed over the undercut area. The rock is assumed to form a muck pile, which is modelled as a continuum of lesser strength (one tenth of the original strength) directly below the cave. Details of the extraction schedule and the draw point have been presented by Paluszny et al. (2015). Fracture propagation is modelled in response to gravity-driven deformation of the cave (Fig. 4).
Material properties correspond to an underground mine. The rock is assumed to be primarily composed of Biotite Quartz Monzonite, with a Young's modulus of 61 GPa, an average volume-weighted UCS of 114 MPa, and a Poisson's ratio of 0.26. Tensile strength is taken as UCS/10, i.e., r t = 11.4 MPa, and K IC is taken as 1.4 MPa m 0.5 .
Two boundary conditions are modelled. In case A, the top surface is fixed in the z direction. This corresponds to a shielded overburden stress caused by the effect of local heterogeneities in the rock matrix, such as a local bridge. In case B, an overburden Fig. 3 Double notch validation tests show that fracture interaction is influenced by boundary conditions as well as spacing. In both cases compression is followed by an increase in tension. In a tension increase occurs over five steps while for b tension increase occurs over twenty steps   (Fig. 4). In both cases, all model sides are fixed in all directions except the vertical, and the bottom boundary (with the exception of the undercut) is fixed.
The two cases are not identical, as they have several key differences in addition to the boundary conditions.
(1) cases are made up of different distributions on initial fracture locations and orientations. Even though all datasets have between thirty and ninety initial diskshaped fractures, these have all slightly different radii and locations for each run. (2) Initial fractures are placed directly above the undercut, occupying 50% of the area across length of the model (x-axis, from left to right for the front views). (3) In addition, in the y and z directions, initial fractures are placed respecting a margin of 10 and 30% to the boundaries. (4) Fractures Fig. 7 Case B. Growth of 30 fractures at (1) 10% and (2) 30% margin away from the boundary Geomech. Geophys. Geo-energ. Geo-resour. (2017) 3:121-130 127 are grown with a beta propagation exponent of 2, and maximum extension length ranges between 10 and 20 m per month (see Zimmerman 2011, 2013).
3.1 Case A: shielded overburden Figure 5 shows sixty orthogonal fractures that propagate due to the growth of the cave. In general, a shielded overburden will lead to the growth of primary fractures close to the cave, in an arching manner. However, it also leads to massive growth away from the cave, directly influencing the shape of the ensuing cave, due to the intersection and coalescence of the interacting fractures. More instances of fracture patterns for case A have been discussed by Paluszny et al. (2015).

Case B: direct overburden
Fractures around an undercut grow in response to caving without any shielding of the overburden stresses. The modality also conduces to the systematic growth of fractures around the cave. Specifically, fractures grow around the undercut and pronounce the formation of the cave by growing in a dome shaped manner. Fractures interact and form a curved concentric discontinuity that further induces caving and accelerates the fragmentation of the rock immediately above the cave.  Figure 6 shows a rendered view of ninety fractures growing during extraction. Fractures are coloured differently to enhance their domains. In Fig. 7, four steps of growth of two different thirty-fracture sets are shown. On the left, fractures leave a 10% margin to the boundary, whereas on the right they leave a 30% margin. In both cases, the resulting pattern is comparable. Fractures closest to the extraction pattern grow first. The same behaviour is observed in Fig. 8, where sixty fractures grow around the undercut. The dome shape is even more pronounced, and interaction plays a role in creating a swarm of fractures that line the top region of the cave. Fractures undergo a complete reorientation caused by local stresses, and little influence of the initial orientation of the fractures is noticed. Fractures are concentric to the undercut, and grow along the regions of the rock that suffers the most differential displacement. Fractures influence the shape of the cave in the case of direct stress boundary conditions, as observed before for shielded boundary conditions.
All simulations run on a Dell Precision Workstation with two Intel Xeon CPUs X5460 @ 3.16 GHz each. Their duration varies between 2 and 6 h, depending on the amount of fractures and amount of extracted years (up to 2 years of extraction).

Conclusions
A finite element based method to model primary fragmentation driven by fracture growth has been applied to examine the influence of boundary conditions on fracture growth patterns during caving. The method represents fractures explicitly, and utilises stress intensity factors and energy-based criteria to model growth. The model grows limited numbers of fractures that can capture the overall behaviour of fracture growth during caving. Fractures are observed to influence the shape of the resulting cave. Results show that when overburden stresses are shielded, such as in the case of the formation of a bridge, fractures tend to grow away from the cave. When stresses are not shielded, fractures tend to grow closer to the cave. In both cases, fractures grow to form a dome shape around and above the undercut. These results can be useful when adjusting extraction schedules in response to cave growth.