Finite element analysis of blast-induced fracture propagation in hard rocks

Computers and Structures 182 (2017) 1–13

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Finite element analysis of blast-induced fracture propagation in hard rocks Marko Bendezu, Celso Romanel ⇑, Deane Roehl Department of Civil Engineering, Pontifical Catholic University of Rio de Janeiro (PUC-Rio), Rio de Janeiro, RJ, Brazil

a r t i c l e

i n f o

Article history: Received 17 May 2016 Accepted 16 November 2016

Keywords: Rock blasting Fracture propagation Extended finite element method Cohesive zone model

a b s t r a c t This paper presents a numerical analysis based on the finite element method to simulate blast-induced hard rock fracture propagation. Three different approaches are compared: the extended finite element method, with a cohesive zone model to represent the growth of fractures; the conventional finite element method using a remeshing technique and based on the linear fracture mechanics; the element deletion method to simulate a rock fragmentation process. The rock mass is a sound granite that remains linear elastic right up the breakage. Two numerical examples are presented in order to discuss the advantages and limitations of each approach. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Rock blasting is generally carried out by drilling into a rock mass, charging the blastholes and firing the igniters located in cylindrical charges. The detonation waves propagate with high velocities (between 2000 and 7000 m/s depending on the type of explosive), which justifies the assumption that the overpressure of the explosive gases begins to act on all points of the blasthole walls simultaneously. The energy released in the explosion is converted into two main forms that are responsible for rock fracturing, creating new cracks and widening the already existing ones: blastinduced stress waves (dynamic load) and the overpressure of the explosive gases (quasi-static load). There is no clear indication on how much energy is converted into stress wave energy, how much is available as high-pressure gases and how much is lost to other sources such as temperature increase, air blast and smoke. The energy partition will depend on the type of explosive, with TNT (trinitrotoluene) and like explosives classified as high in stress wave and low in gas production, while ANFO (ammonium nitrate/fuel oil) and others are considered high in gas production and low in stress wave energy [1]. The blast-induced P waves that travel out into the rock mass provoke sudden increases in the normal compressive stresses along the radial direction and normal tensile stresses in the tangential direction. Crushing of the rock adjacent to the blast hole ⇑ Corresponding author at: Department of Civil Engineering, Pontifical Catholic University of Rio de Janeiro (PUC-Rio), Rua Mq. S. Vicente 225 – Gavea, 22451-900 Rio de Janeiro, RJ, Brazil. E-mail addresses: [email protected], [email protected] (C. Romanel). http://dx.doi.org/10.1016/j.compstruc.2016.11.006 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

occurs whenever the resulting compressive stresses exceed the rock compressive strength and a system of radial fractures, emanating from the drill hole, arises as consequence of the high tangential tensile stresses. Since the problem geometry is generally bounded by a free surface (a soil-air interface), compressive P waves reflect back from the free surface as tensile P waves, providing an important contribution to the rock fracturing, in addition to S waves generated by a modal conversion phenomenon [2]. These reflected waves, by their turn, will give rise to new S and P waves due to a multiple reflection mechanism that involves both the free surface and the surfaces of the growing fractures. The region of dense radial fracturing around the blast hole extends to a distance that depends on the amplitudes of the stress waves, the mechanical properties of the rock and the characteristics of the explosive such as detonation velocity, type of contact between the charge and the blasthole walls. Reports in the literature [3–5] suggest that when ANFO explosive is used this region generally extends between 4 and 8 blasthole radiuses away. There is not practical interest in using high-intensity explosives that will generate very high pressures on the borehole walls. The crushing and the radial fracturing of the nearby rock will consume energy, but will only contribute with a very small volume of excavation, besides producing unnecessary damage to the rock surfaces and affecting the strength and stability of the remaining material. Beyond the region of severe damage, some main fractures continue to grow around the blasthole perimeter. In the literature, there is no agreement about the number of dominant fractures but some numerical results [6,7] reported the existence of 8–12 dominant fractures around the blasthole when TNT explosive is employed. The number of dominant fractures depends on the char-

2

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

acteristics of the explosive, represented by the shape and duration of the pressure pulse applied on the borehole walls [8]. The fracturing of rock masses can be numerically simulated by the conventional finite element method (FEM) but changes in the fracture topology with time will require remeshing of the entire domain, a severe limitation for situations of fracture propagation involving complex geometries. The use of adaptive remeshing schemes may also be awkward to implement because of the large computational burden and the consideration of multiple fractures will make remeshing almost intractable and inapplicable. The extended finite element method [9] (XFEM) is an alternative approach that is particularly useful for the approximation of solutions with severe non-smooth characteristics in small parts of the domain. The finite element space is enriched with discontinuous functions that give a greater accuracy and computational efficiency to the numerical solution when compared to the conventional FEM applications. Furthermore, the mesh is not required to match the geometry of the fractures since they are not physically represented. XFEM is a very attractive and effective way to simulate the propagation of fractures along an arbitrary, solution-dependent path that eliminates the need for remeshing. This paper investigates the dynamic fracturing of a sound rock mass, admitted as an isotropic and homogeneous material that remains linear elastic right up the moment of breakage. The behavior of hard rocks containing a large percentage of quartz closely corresponds to this material behavior, such as granite. Three methods are compared in this research: the extended finite element method (XFEM), the conventional finite element method (FEM) using a remeshing technique and the element deletion method that simulates the evolution of a rock fragmentation process. Different types of numerical methods may be also useful to investigate the dynamic fracture propagation in solids. Among others, can be mentioned the molecular dynamic simulation [10], meshfree methods [11], the discrete particle method [12] and a version of the lattice model where a network of truss elements is used in the zone where the fracture process is expected to take place and a FEM element mesh may be considered in the surrounding areas [13]. The blasting effects on the rock immediately around the drill hole (crushed zone) are ignored and fracture propagation is considered due to the blast-induced stress waves only. The initial stresses existing near the blasthole are also neglected since the stress increments generated by the detonation are admitted much larger than the original stress state. The elastodynamic problem is formulated under plane strain condition, an assumption that is approximately satisfied only at the plane normal to the longitudinal axis of the cylindrical blasthole passing through its mid-length; a true three-dimensional simulation would have been very complex, lengthy and expensive for the time being. 2. Extended finite element method The main concept in XFEM is to enrich the usual finite element space with additional degrees of freedom that allow fractures to open and increase the accuracy of the approximation near the fracture tip. The functions in the XFEM space have the general form [14]:

uðxÞ ¼

X

ui ui ðxÞ þ

i

þ

X bj uj ðxÞHðxÞ j

X k

! 4 X uk ðxÞ c‘k F ‘ ðrðxÞ; hðxÞÞ

ð1Þ

‘¼1

where fui g are the conventional interpolation functions, H(x) the Heaviside function associated to the current fracture geometry,

fbi g the enrichment degrees of freedom associated with fracture separation away from the tip, fc‘k g the enrichment degrees of freedom associated with near-tip displacement and fF ‘ ðr; hg the asymptotic branch functions expressed in polar coordinates, from the fracture tip, as: fF ‘ ðr; hg ¼

pffiffiffi pffiffiffi pffiffiffi pffiffiffi  r sinðh=2Þ; r cosðh=2Þ; r sinðh=2Þ sinðhÞ; r cosðh=2ÞsinðhÞ ð2Þ

The sum over i in Eq. (1) is taken considering all the nodes in the mesh while the set of nodes j over which the second sum is performed contains all the nodes belonging to an element entirely cut by the fracture. The set of nodes k is built in such a way that it contains all nodes located within a certain distance from the fracture tip [15,16]. XFEM requires additional degrees of freedom (enrichment) in order to allow fractures to open. In a region separated by a fracture into two pieces away from its tip, the enrichment is provided by a Heaviside function defined to be 1 on one side of the fracture and 1 on the other side. This technique is easy to implement in the case of a single straight fracture but it is more complicated as the geometry of the fracture becomes irregular. It has been shown [17,18] that the discontinuity in the displacement field with the Heaviside function is equivalent to the addition of an extra element on top of an existing one, a method that has become increasingly popular and is known as the phantom node method. Consider an element with nodes n1, n2, n3, n4 (Fig. 1). The fracture Cc divides the element domain into two complementary subdomains, XA and XB. In the phantom node method, a discontinuity in the element displacement field is constructed by adding phan~2 , n ~3 , n ~ 4 ) on top of the existing nodes. ~1 , n tom nodes (here labeled n The finite element is replaced by two new elements, referred to as elements A and B. The connectivity of the overlapping elements is ~1 , n ~ 2 , n3, n4] for the new element A and [n1, n2, n ~3 , n ~ 4 ] for the new [n element B. The elements do not shares nodes, and therefore have independent displacement fields. Both elements are only partially active; the active part of element A is denoted by XA and the active part of element B by XB. The displacement uðxÞ of a point with coordinates x is computed with the standard finite element interpolation functions fuðxÞg and the integrations are carried out considering either the subdomains XA or XB depending on the location of the point with respect to the fracture. Closure of the fracture tip is enforced when no phantom nodes are added on the element boundary that contains the tip; thus, the overlapping elements do share nodes and their displacement fields are not completely independent. Since the interpolation functions associated to the enriched elements are the same as for the intact elements, the phantom node method is easy to be implemented, being available in several commercial FEM solvers [20,21]. While XFEM is suitable for additional enrichment of the displacement field around the fracture tip to capture the singular field (third term of Eq. (1)), the method of phantom nodes is only applicable to model cohesive fractures, where the singularity in the stress field is removed due to the presence of a cohesive traction. The discontinuity generally grows elementwise, with the fracture tip located at an element boundary, although variations in the phantom element method may admit the tip inside the finite element [22]. In the cohesive zone model [23–25] the process region is admitted as an extension of the fracture length up to a point called fictitious fracture tip (or mathematical crack tip) in which a specific constitutive law is considered in order to relate stress decreases with increase in fracture opening. The real fracture tip (or physical fracture tip) is the point on the fracture surface where there is no stress and the normal opening is bigger than the critical opening.

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

3

Fig. 1. Connectivity and active parts of two overlapping elements in the phantom node method [19].

Given the importance of the process region in fracture mechanics, it is crucial to consider a model that takes into account the specific behavior of material within this very small region. In finite element analysis, a difficulty is that the energy dissipation through fracture surfaces converges to zero by mesh refinement, a phenomenon denoted as spurious mesh sensitivity, physically not accepted. A possible solution may be achieved by specifying the energy dissipated over fracture surfaces, and one of the simplest approaches to model the mechanical behavior around a fracture tip in brittle materials is the cohesive zone model in which the fracture characteristics are embedded in a surface traction vs fracture opening constitutive relationship. An inherent length scale is introduced into the model and no fracture criterion (K-dominant field) is required; outside the cohesive zone the material behavior is admitted linear elastic. The area under the surface traction - fracture opening curve represents fracture energy and in the cohesive zone model the energy release rate is specified as a material parameter. The required parameters, cohesive strength and critical separation displacement, must be determined from careful fracture experiments. The onset of damage refers to start of degradation of the response at a point. A maximum principal stress criterion has been considered in this research to predict damage, which is assumed to initiate whenever the following interaction function reaches unity:

f ¼

hrmax i ¼1 T max

ð3Þ

where Tmax is the tensile strength of the rock and rmax the maximum principal stress; the Maculay brackets are used to indicate that a compressive stress state does not initiate damage, i.e. hrmax i = 0 if rmax < 0. The damage evolution describes the rate at which the toughness of the material is degraded once the corresponding initiation criterion is reached. For mixed mode failure, the most widely used criterion to predict fracture propagation is the power law criterion [26]:



GI GIC

g

 þ

GII GIIC

g ¼1

ð4Þ

where g is a material parameter, GI and GII are the energy release rates in the normal (mode I) and shear (mode II) directions, respectively, GIC and GIIC denote the corresponding critical fracture energy rates required to cause failure. The mixed mode energy release rate can be expressed as GT = GI + GII which, for fracture propagation, should be equal to the critical energy release rate GC,

,

GC ¼ GI þ GII ¼ 1

m1 GIC

g

 þ

where m1 = GI/GT and m2 = GII/GT.

m2 GIIC

g 1=g ð5Þ

3. Conventional finite element method In the linear elastic fracture mechanics the stress components in the fracture tip region can be described by single-term parameters, known as stress intensity factors, whose analytical expressions are available for some simple cases. Fracture advances can be monitored by the attainment of a critical value KC (quasistatic loading) or KD (dynamic loading), considered to be a material property. For the mixed mode I–II, the literature reports several criteria for condition of unstable fracture growth under static loading [27], among them the following elliptical criterion:



KI K IC

2

 þ

K II K IIC

2 ¼1

ð6Þ

The orientation of fracture propagation is determined as the direction hm of maximum tensile tangential stresses,

2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3  2  2 1 KI 1 KI 4 hm ¼ 2 arctan  þ 85 4 K II 4 K II

ð7Þ

For a dynamic loading condition, it has been observed that both the tensile strength of the rock and the strain energy release rate are significantly greater than those measured quasi-statically. The authors used in this paper a criteria for mixed mode I-II (Eq. (6)), strictly developed for quasi-static loads, but replacing the critical stress intensity factors KIC and KIIC by the corresponding dynamic values KID and KIID. The stress intensity factors are numerically determined using triangular quadratic elements with the mid-nodes moved to the quarter points adjacent to the fracture tip (quarter-point singular elements) [28,29]. It can be shown [30] that for triangular elements the radial strain component along all rays emanating from the pffiffiffi node at the fracture tip possesses a 1 r singularity. Eight triangular quarter-point elements with a 45° internal angle at the common vertex are connected around the fracture tip, as shown in Fig. 2. The stress intensity factors for mixed mode I-II are obtained by the modified crack-closure method [31]. The strain-energy release rate equations for mode I (GI) and mode II (GII) are given by:

1 F yi ½t11 ðvm  vm0 Þ þ t 12 ðvl  vl0 Þ 2da 1  F yj ½t21 ðvm  vm0 Þ þ t 22 ðvl  vl0 Þ 2da

GI ¼ 

1 F xi ½t11 ðum  um0 Þ þ t 12 ðul  ul0 Þ 2da 1 F xj ½t21 ðum  um0 Þ þ t 22 ðul  ul0 Þ  2da

ð8aÞ

GII ¼ 

ð8bÞ

4

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

Fig. 2. Finite element rosette at fracture tip (left) and vertical nodal forces in quarter-point triangular element (right).

4. Element deletion method The element deletion method is one of the simplest approaches to simulate blast-induced failure of brittle rocks within the framework of the conventional FEM. There is no need to represent strong discontinuities in displacement fields since the discontinuity is created by gradual element deletion from the mesh. Elements are not physically eliminated, but their contributions are ignored by prescribing that the stress state within deleted elements is close to zero. In order to avoid numerical problems related to strong local loss of equilibrium, this reduction is made in several relaxation steps and the modulus of elasticity of the deleted element is approximated to zero in the last relaxation step. The Rankine failure criterion detects the initiation of the fracture, when the maximum principal stress exceeds the tensile strength of the rock. Although this initiation is based on mode I only, the analysis can subsequently incorporate both modes I and II through a smeared crack model [35], in which the nucleation of one or more cracks is translated into a deterioration of the current stiffness and strength at the element integration points. Song et al. [36] reported that results obtained by the element deletion method are extremely sensitive to the structure and refinement of the finite element mesh, being more useful in studies on rock fragmentation than in true analysis of fracture propagation. Due to the simplicity of computational implementation, the

element deletion method is available in some general-purpose finite element programs such as LS-Dyna and Abaqus. 5. Dynamic model The fast conversion of chemical energy of the explosive generates a detonation pressure pD ranging from 0.5 to 50 GPa that, in consequence, creates a pressure pulse on the blasthole walls whose characteristics depend on the mechanical properties of the rock, the amount of produced gas, temperature, detonation velocity and the contact between the explosive and the walls, among other factors. It is generally agreed that the maximum pressure p0 applied on the blasthole walls is something less than the detonation pressure calculated from chemical considerations [37,38]. An appropriate pulse pressure can be approximated by either applying pD over an expanded cavity radius [39] or estimating p0 from empirical correlation with pD [40,5]. The pressure pulse has two distinct features depending on the explosive type: high-order or low-order detonations. A highorder detonation yields a pressure pulse showing a rapid growth to a well-defined peak, followed by a sharp pulse fall, while the low-order detonation produces a pressure pulse of smaller peak intensity and a longer attenuation in time [4]. High-order detonation is provoked by TNT and emulsion explosives while the loworder detonation is generated by ANFO explosives. Closed-form solutions for propagation of elastic waves from a pressurized spherical cavity generally use the following pressure decay function [41,42]

 pðtÞ ¼ po  eat  ebt

ð9Þ

where po is the maximum pressure applied by the pulse, t is time and a and b are positive frequency-dependent decay constants. 3000 2500 2000

p0 [MPa]

where Fxi and Fyj are the nodal forces at the ith node in the x and y directions, respectively, um and vm the displacement components of the mth node in the x and y directions, respectively, t11 = 6  3p/2, t12 = 6p  20, t21 = 1/2 and t22 = 1. Fig. 2 also shows the nodal forces in the y-direction. In the conventional FEM the geometry of the finite element mesh must vary continually in time as the dominant fractures propagate into the rock mass. In this research, an auto adaptive algorithm is used for generation of finite element meshes based on the ‘‘quadtree” scheme combined with a boundary triangulation technique [32]. Assuming that rock fracturing occurs under constant velocity, the fracture growth at a given time interval can be computed without difficulty. The velocity of fracture propagation must be lower than the propagation of R waves [33], being generally taken less than half of the S wave velocity [34]. As the problem geometry is modified by the dynamic fracturing, the values of nodal displacements, nodal velocities and nodal accelerations must be mapped from the old configuration at time t  Dt to the current one at time t. The mapping technique is easy to implement but requires a considerable computational effort for the identification of the corresponding nodes and elements in the two finite element models.

1500 1000 500 0

0

1000

2000

3000

4000

5000

6000

t [ s] Fig. 3. Pressure pulse applied on the blasthole walls with maximum pressure p0 = 1000 MPa and p0 = 2500 MPa when t0 = 450 ls.

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

5

Fig. 4. Finite element mesh in XFEM analysis.

Fig. 5. XFEM model at time t = 15.7 ls: (row a) distribution of major principal stresses (Pa) and finite element mesh around crack 1 and (row b) status of enriched elements; STATUSXFEM = 1 when element is completely cracked and STATUSXFEM = 0 when it is still intact [16].

6

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

Similar expressions have been applied to model long columns of explosives as sequential detonations of stacked spherical charges [43,44]. This research uses a modified form of Eq. (9) [8],

pðtÞ ¼ p0  nðeat  ebt Þ

ð10aÞ

n ¼ 1=ðeat0  ebt0 Þ

ð10bÞ

t 0 ¼ ð1=ðb=aÞÞ  logðb=aÞ

ð10cÞ

where n is a variable to represent the rising and decaying phases of the pressure pulse and t0 is the time required to achieve the maximum pressure p0.

As fractures advance into the rock mass, their surfaces tend to open and to close as consequence of the complex dynamic interaction among all dominant fractures. The condition of nonpenetration must be guaranteed by both surfaces of any given fracture at every time step. This impenetrability constraint may be approximately satisfied with a penalty formulation [45,46] in which the contact condition is enforced by imposing equal and opposite forces on the fractures surfaces along the contact region. The compressive behavior between the surfaces does not interact with the cohesive behavior, which contributes to the stress state only when the fracture is open. Another aspect in finite element analysis requiring careful control is the choice of the element size because large elements are

Fig. 6. Conventional FEM analysis at time t = 15.7 ls with distribution of major principal stresses (MPa) and finite element mesh around crack 1.

Fig. 7. Finite element mesh in the vicinity of the blasthole when neglecting penetration control: (a) XFEM at time t = 145.6 ls and (b) conventional FEM at time t = 523 ls.

Fig. 8. Finite element mesh in the vicinity of the blasthole when imposing penetration control: (a) XFEM at time t = 523 ls and (b) conventional FEM at time t = 523 ls.

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

7

Fig. 9. Horizontal displacements in the vicinity of the blasthole with no penetration control: (a) XFEM at time t = 145.6 ls and (b) conventional FEM at time t = 523 ls.

Fig. 10. Horizontal displacements in the vicinity of the blasthole with penetration control: (a) XFEM at time t = 523 ls and (b) conventional FEM at time t = 523 ls.

Fig. 11. Principal stress distribution and fracture propagation considering penetration control in XFEM analysis.

unable to transmit motions at high frequencies. From several researches published in the literature [47–49] it is suggested that the element size for effective transmission should not be more than 1/10 of the shear wavelength at the maximum frequency (cutoff frequency). Material damping is admitted null since the rock mass is considered as an elastic body. In order to take into account the radiation condition, avoiding the reflection of the stress waves striking at the boundaries of the finite size model, viscous boundary conditions [48] have been

implemented in the conventional FEM analysis while in the XFEM analysis infinite elements have been used. 6. Numerical examples 6.1. Single blasthole This example analyzes the dynamic rock fracturing around a blasthole of radius a0 = 2.54 cm located 40 cm from a free surface.

8

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

Fig. 12. Fracture propagation pattern in a glass plate [51].

Eight initiation cracks of length a0/3 are admitted equally spaced along the blasthole perimeter. The granite properties are Young’s modulus E = 60 GPa, Poisson’s ratio m = 0.25, mass density q = 2800 kg/m3, critical stress intensity factors KID = 1.65 MPa m0.5 and KIID = 1.03 MPa m0.5, energy release rates GIC = 42.5 Pa m and GIIC = 16.6 Pa m, tensile strength Tmax = 3 MPa. The fracture propagation speed is 1464 m/s, corresponding to half the shear wave velocity. The cohesive zone model uses the power law criterion (Eq. (4)) with g = 1.

The pressure pulse on the blasthole walls due to detonation is shown as the tallest curve in Fig. 3, with p0 = 2500 MPa, t0 = 450 ls and the decay constants a and b, experimentally determined for a granite rock [4], with a decay ratio b/a = 3.1. Fig. 4 shows the finite element mesh for the XFEM analysis [20] formed by 7465 four-node finite elements and 59 infinite elements. For the conventional FEM analysis, using a Fortran code written for this research, the mesh consisted of triangular quadratic elements (6-node elements), in a variable number depending on the successive meshes. In the blasthole detail, the starting location of pre-existing cracks, inserted as enriched elements, is identified by black elements. The status of an enriched element (STATUSXFEM) is labeled 1 when it is completely cracked or 0 when still intact. At time t = 15.7 ls after detonation, Fig. 5a presents the distribution of the major principal stresses around the blasthole (positive sign for tensile stresses) and in the neighborhood of crack 1. The maximum tensile stresses act near the fracture surfaces where the enriched elements, as shown in Fig. 5b, are already complete or partially cracked. The stresses on the fracture surfaces decrease with increasing opening of the fracture, although not suddenly to zero, given the concept of fictitious fracture tip assumed in the cohesive zone model. It can be also observed that the fracture does not follow the element boundaries and there is no need to refine the mesh around the fracture tip. Fig. 6 presents a comparison with the conventional FEM analysis [50] showing the degree of mesh refinement around the fracture tip required to capture the singular asymptotic stress fields adequately. The maximum tensile stresses are developed at the fracture tips in the conventional FEM analysis (linear facture

Fig. 13. Fragmentation pattern with time in the element deletion method.

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

9

Fig. 14. Fragmentation pattern at time t = 524 ls with three different meshes in the element deletion method.

4

3

0.4 m

blast hole

5 6

7

12

2

13

1 8

14 0.4 m

11 blast hole

15

10 9 16

free surface Fig. 15. Two blastholes detonated simultaneously in a sound granite rock mass.

mechanics) while in the XFEM analysis they occur along the fracture surfaces (cohesive zone model). In both methods (XFEM, conventional FEM) the simulation has been ran twice: firstly, penetration control was imposed through a penalty formulation, while in the second run this constraint was simply disregarded to permit a comparison between both approaches. Under the condition of no penetration control in the conventional FEM analysis the surfaces of fracture 5 overlap themselves at time t = 523 ls (Fig. 7), causing a program interruption, while in the XFEM analysis this interruption occurred at an earlier time (t = 145.6 ls) when surfaces of fractures 5 and 7 begin to superimpose. Fig. 8 shows the finite element meshes near the blasthole at time t = 523 ls with the imposition of penetration control, resulting that fractures propagate without numerical difficulty in both finite element models. The corresponding horizontal displacement field is shown in Figs. 9 and 10 while Fig. 11 illustrates the distribution of principal stresses in the XFEM analysis with penetration control. Fig. 10 also shows that as the fractures advance, the blocks between fractures 4 and 5 (Fig. 5) and between fractures 5 and 6 tend to be expelled from the rock mass. Fractures 4 and 6 do not intersect the free surface and grow following a path almost parallel to the free boundary. A similar behavior (Fig. 12) was experimentally observed by Porter [51] pressurizing a borehole with oil in glass plates to level near a critical value, followed by a sudden increase in pressure using an explosive wire.

The element deletion method has been also used to simulate the same problem, incorporating a brittle Rankine failure-type material model. The finite element mesh consisted of 11,935 four-node elements and 82 infinite elements, without any preexisting crack around the blasthole. The pressure of the detonation gases is also indirectly considered besides the stress waves generated by the transient pulse (Fig. 3) applied on the borehole walls. Fig. 13 presents a ring of deleted elements whose radius rapidly expands from 0.066 m at t = 100 ls after explosion to 0.492 m at t = 524 ls, when numerical simulation has been finished. In fact, it is not possible to observe the formation of any dominant fracture because this method rather represents a rock fragmentation process than a fracture propagation analysis. It is also clear that results depend strongly on the type and number of finite elements, as shown in Fig. 14 considering three different finite element meshes at time t = 524 ls: (a) 11,935 four-node elements; (b) 34,874 fournode elements; and (c) 31,557 three-node elements. Despite the large number of elements, the computer runtime is low when compared to XFEM runtime (7465 four-node elements) since the element deletion method is based on an explicit time integration scheme. 6.2. Two blastholes The second example analyzes the blast-induced rock fracturing around two blastholes of radius a0 = 2.54 cm drilled into a sound granite rock. Eight initiation cracks of length a0/3 are introduced, equally spaced along their perimeters (Fig. 15). The rock properties are the same already mentioned in the previous example. The pressure pulse on the blasthole walls provoked by the simultaneous explosions is graphically shown in Fig. 3 with a pressure peak p0 = 1000 MPa. The mesh for the XFEM analysis is formed by 37,040 four-node finite elements and 428 infinite elements to ensure the radiation condition. From Fig. 16 it can be observed that at time t = 60.7 ls fractures 1 and 13, parallel to the free surface, intersect themselves and at time t = 138.6 ls the fractures 7 and 15 reach the free surface. At time t = 450 ls after detonation, when the pressure pulse attains its peak p0 = 1000 MPa, some dominant fractures deviate from the radial direction, such as fractures 6 and 16, following a path almost parallel to the free surface. Fractures 2, 12, 8 and 14 stop growing, with their tips no longer experimenting a strong tensile stress field to foster propagation, as can been seen in the right col-

10

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

Fig. 16. Fracture patterns and corresponding tensile stress distributions at times t = 60.7 ls (top row), t = 138.6 ls (middle row) and t = 450 ls (bottom row) from XFEM analysis.

umn of Fig. 16 that shows the tensile stress fields at three different instants. This problem has been re-analyzed through the conventional FEM. Fig. 17 indicates quite similar fracture patterns obtained with XFEM at the final stage of analysis t = 538 ls. At this time, the analysis with the conventional FEM was interrupted due to difficulties of the remeshing algorithm to create the singular element rosettes (Fig. 2), what happens whenever fractures tends to intersect each other. Another numerical simulation has been carried out using XFEM, this time without any pre-existing cracks along the blasthole

perimeters. This approach is possible since the cohesive fracture model is adopted, which allows fracture initiation when the cohesive response in an enriched element begins to degrade, according to the chosen maximum principal stress criterion (Eq. (3)). Fig. 18 represents the fracture patterns at time t = 90.7 ls and t = 404 ls, being possible to observe that many small cracks are initially developed, but just some of them predominate with time. The number of 8–12 predominant fractures [6,7] seems to be confirmed, although fractures are not evenly distributed along the perimeters of the holes and their number is greater in the second hole. Although this approach based on the absence of pre-

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

11

Fig. 17. Fracture patterns from XFEM (left) and conventional FEM (right) analysis at time t = 538 ls.

Fig. 18. Fracture patterns in XFEM analysis without pre-existing cracks at times t = 90.7 ls and t = 404 ls.

existing cracks is possible, it brings some inconveniences. The computer runtime is much higher when compared to the analysis with pre-existing cracks and the probability of numerical difficulties that may arise due to fracture intersections is also greater. Finally, for the same problem, the element deletion method has been applied and the numerical results are presented in Fig. 19, showing the evolution in time of rock fragmentation. It can be observed that around the blastholes there is a very fragmented rock ring, without any definition of dominant fractures. 7. Conclusion This paper dealt with the application of three different approaches, based on the finite element method, to investigate blast-induced fracturing of hard rocks. In the conventional FEM the finite element mesh must vary continually in time as the dominant fractures propagate into the

rock mass. This requirement constitutes its main disadvantage since as the problem geometry is modified the values of nodal displacements, nodal velocities and nodal accelerations must be mapped from the old mesh to the current one, requiring a considerable computational effort for the identification of the corresponding positions of nodes and elements between the two finite element models. The element deletion method is a simple approach with no need to represent strong discontinuities in the displacement fields since the discontinuity is created by gradual element elimination from the mesh. The Rankine failure criterion detects the initiation of the fracture, when the maximum principal stress exceeds the tensile strength of the rock. Its main disadvantages are the extreme dependence of the numerical results on the structure and refinement of the finite element mesh and also due to the fact that the element deletion method represents rather a rock fragmentation event than a fracture propagation process.

12

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13

Fig. 19. Evolution in time of rock fragmentation with the element deletion method.

The extended finite element method XFEM is an alternative approach that is particularly useful for the approximation of solutions with severe non-smooth characteristics in small parts of the domain such as in fracture propagation problems. The finite element space is enriched with discontinuous functions that give a greater accuracy and computational efficiency to numerical solutions when compared to the conventional FEM applications. The mesh is not required to match the geometry of the fractures since they are not physically represented. In the phantom node method, used in this research, a finite element is replaced by two new elements constructed by adding new nodes (phantom nodes) on top of the existing nodes and the singularity in the stress field is removed by considering a cohesive zone model near the fracture tip (process region). One disadvantage of XFEM application is that discontinuity is generally admitted growing elementwise, with the fracture tip always located on an element boundary. If a coarse mesh is adopted then fracture propagation and traction distribution within the cohesive zone may not be properly captured. Moreover, in dynamic problems the mesh refinement depends on the wavelength at the maximum frequency, since large elements are unable to transmit motion at high frequencies. In the conventional FEM analysis, the degree of mesh refinement around the fracture tip is also high in order to represent the singular asymptotic stress fields adequately. An alternative approach would be to replace the linear elastic behavior by an elastoplastic constitutive stress-strain law to avoid stress singularities in the fracture tip region. In this research, a brittle material behavior was assumed for hard rocks containing a large percentage of quartz. This assumption made also possible the qualitative comparison of the conventional finite element results with the experimental results obtained by Porter [51] considering fracture propagation in glass, as shown in Fig. 12.

References [1] Fourney WL. Mechanisms of rock fragmentation by blasting. In: Hudson JA, editor. Comprehensive rock engineering; principles, practice & projects. Excavation, support and monitoring, vol. 4. England: Pergamon Press Ltd.; 1993. p. 39–69. [2] Achenbach JD. Wave propagation in elastic solids. 1st ed. North Holland; 1987. [3] Siskind DE, Fumanti RR. Blast-produced fractures in Lithonia granite. report of investigation 7901. U.S. Bureau of Mines; 1974. [4] Aimone CT. Three-dimensional wave propagation model of full-scale rock fragmentation [Ph.D. Dissertation]. Northwestern University; 1982. [5] Jimeno CL, Jimeno EL, Carcedo FJA. Drilling and blasting of rocks. Taylor & Francis; 1995. [6] Ghosh A, Daemen JJK. Rock fragmentation in bench blasting – a numerical study. In: Daemen JJK, Schultz RA, editors. Proceedings of the 35th U.S. symposium on rock mechanics. p. 553–8. [7] Song J, Kim K. Blasting induced fracturing and stress field evolution at fracture tips. In: Daemen JJK, Schultz RA, editors. Proceedings of the 35th U.S. symposium on rock mechanics. p. 547–52. [8] Cho SH, Kaneko K. Influence of the applied pressure waveform on the dynamics fracture processes in rock. Int J Rock Mech Min Sci 2004;41 (5):771–84. [9] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [10] Furuya Y, Noguchi H. A combined method of molecular dynamics with micromechanics improved by moving the molecular dynamics region successively in the simulation of elastic-plastic crack propagation. Int J Fract 1998;94:17–31. [11] Belytschko T, Lu YY, Gu L. Crack propagation by element-free Galerkin methods. Eng Fract Mech 1995;51:295–315. [12] Rabczuk T, Song JH, Belytschko T. Simulations of instability in dynamic fracture by the cracking particles method. Eng Fract Mech 2009;76:730–41. [13] Schlangen E, Garboczi FJ. Fracture simulations of concrete using lattice models: computational aspects. Eng Fract Mech 1997;57:319–32. [14] Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. [15] Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. Int J Numer Meth Eng 2005;64:354–81. [16] Béchet E, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. Int J Numer Meth Eng 2005;64:1033–56. [17] Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Meth Appl Mech Eng 2004;193:3523–40.

M. Bendezu et al. / Computers and Structures 182 (2017) 1–13 [18] Song JH, Areias PMA, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. Int J Numer Meth Eng 2006;67:868–93. [19] Van Der Meer FP, Sluys LJ. A phantom node formulation with mixed mode cohesive law for splitting laminates. Int J Fract 2009;158(2):107–24. [20] Abaqus 6.12 user’s manuals. Pawtucket (RI): Abaqus Inc.; 2014. [21] Prévost J. Dynaflow – a nonlinear transient finite element analysis program. Princeton University; 1983. [22] Rabczuk T, Zi G, Gerstenberger A, Wall WA. A new crack tip element for the phantom node method with arbitrary cohesive cracks. Int J Numer Meth Eng 2008;75(5):577–99. [23] Barenblatt GI. The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses. J Appl Math Mech 1959;23:22–36. [24] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. Adv Appl Mech 1962;7:55–129. [25] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [26] Wu EM, Reuter Jr RC. Crack extension in fiberglass reinforced plastics. T. & AM report 275. University of Illinois; 1965. [27] Whittaker BN, Singh RN, Sun G. Rock fracture mechanics: principles, design e applications. 1st ed. Elsevier Science Ltd.; 1992. [28] Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. Int J Numer Meth Eng 1976;10:25–37. [29] Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. Int J Numer Meth Eng 1975;9:495–507. [30] Owen DRJ, Fawkes AJ. Engineering fracture mechanics. Swansea (UK): Pineridge Press Ltd.; 1983. [31] Raju IS. Calculation of strain-energy release rates with higher order and singular finite elements. Eng Fract Mech 1987;28(3):251–74. [32] Araujo TDP, Cavalcante Neto JB, Carvalho MTM, Bittencourt TN, Martha LF. Adaptive simulation of fracture processes based on spatial enumeration techniques. Int J Rock Mech Min Sci 1997;34:188.e1–188.e14. [33] Kostrov BV. On the crack propagation with variable velocity. Int J Fract 1975;11(1):47–56. [34] Grady DE, Kipp ME. The micromechanics of impact fracture of rock. Int J Rock Mech Min Sci 1979;5(16):293–302. [35] Saharan MR, Mitri HS. Numerical procedure for dynamic simulation of discrete fractures due to blasting. Rock Mech Rock Eng 2008;41(5):641–70.

13

[36] Song JH, Wang H, Belytschko T. A comparative study on finite element methods for dynamic fracture. Comput Mech 2008;42:239–50. [37] Kutter HK, Fairhurst C. On the fracture process in blasting. Int J Rock Mech Min Sci 1971;8:181–202. [38] Laroque GE, Favreau RF. Blasting research at the mines branch. In: Proceedings twelfth symposium of rock mechanics. Rolla (USA): University of Missouri; 1970. p. 341–51. [39] Kutter HK. The interaction between stress wave and gas pressure in the fracture process of an underground explosion in rock, with particular application to presplitting [Ph.D. Dissertation]. University of Minnesota; 1967. [40] Persson P-A, Holmberg R, Lee J. Rock blasting and explosives engineering. CRC Press; 1994. [41] Duvall WI. Strain-wave in rock near explosions. Geophysics 1953;18 (2):310–23. [42] Favreau RR. Generation of strain waves in rock by an explosion in a spherical cavity. J Geophys Res 1969;74(17):4267–80. [43] Chua KM, Aimone C, Majtabai N. A numerical study of the effectiveness of mechanical rock bolts in an underground opening excavated by blasting. In: Proceedings of the 33th U.S. symposium on rock mechanics, New Mexico. p. 285–93. [44] Dowding CH, Aimone CT. Multiple blast-hole stresses and measured fragmentation. Rock Mech Rock Eng 1985;18(1):17–36. [45] Papadopoulos P, Jones RE, Solberg JM. A novel finite element formulation for frictionless contact problems. Int J Numer Meth Eng 1995;28:2603–17. [46] Saracibar CA. A new frictional time integration algorithm for large-slip multibody frictional contact problems. Comput Meth Appl Mech Eng 1997;142:303–34. [47] Celep Z, Bazant Z. Spurious reflection of elastic waves due to gradually changing finite element size. Int J Numer Meth Eng 1983;19:631–46. [48] Kuhlemeyer RL, Lysmer J. Finite element method accuracy for wave propagation problems. J Soil Mech Found Div, ASCE 1973;99(5):421–7. [49] Mullen R, Belytschko T. Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation. Int J Numer Meth Eng 1982;18:11–29. [50] Lima ADR, Romanel C, Roehlm DM. Dynamic fracturing of hard rocks by blasting. In: Proceedings of the 8th NUMOG. p. 345–51. [51] Porter DD. A role of the borehole pressure in blasting: the formation of cracks [Ph.D. Dissertation]. University of Minnesota; 1970.