See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/284354574

A preliminary evaluation of an enhanced FDEM code as a tool to

simulate hydraulic fracturing in jointed rock masses

Chapter · May 2014

DOI: 10.1201/b16955-248

CITATIONS READS

11 1,288

5 authors, including:

Andrea Lisjak Omid Mahabadi

Geomechanica Inc. Geomechanica Inc.
55 PUBLICATIONS 2,246 CITATIONS 45 PUBLICATIONS 1,828 CITATIONS

SEE PROFILE SEE PROFILE

Tim Vietor
Nationale Genossenschaft für die Lagerung radioaktiver Abfälle
86 PUBLICATIONS 1,395 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

The Use of Fiber Optics in order to Determine the Behaviour of Ground Support Elements View project

Analysis and prediction of the excavation damaged zone (EDZ) in Opalinus Clay View project

All content following this page was uploaded by Andrea Lisjak on 10 December 2015.

The user has requested enhancement of the downloaded file.

 Rock Engineering and Rock Mechanics: Structures in and on
Rock Masses – Alejano, Perucho, Olalla & Jiménez (Eds)
© 2014 Taylor & Francis Group, London, 978-1-138-00149-7

A preliminary evaluation of an enhanced FDEM code as a tool to

simulate hydraulic fracturing in jointed rock masses

A. Lisjak, O.K. Mahabadi & P. Kaifosh

Geomechanica Inc., Toronto, ON, Canada

T. Vietor
National Cooperative for the Disposal of Radioactive Waste (NAGRA), Wettingen, Switzerland

G. Grasselli
Department of Civil Engineering, University of Toronto, Toronto, ON, Canada

ABSTRACT: Experimental evidence based on microseismic data clearly show that fluid-pressure-driven frac-
tures interact with preexisting discontinuities, thus highlighting the strong influence of rock mass structures on
hydraulic fracture development. However, these mechanisms are not well accounted for by analytical models
and conventional continuum-based numerical approaches. The purpose of this paper is to assess the use of
an alternative hybrid finite-discrete element (FDEM) code, enhanced with hydraulic fracturing capabilities, to
model pressure-driven fracturing in jointed rock masses. The proposed approach is first validated by comparing
the emergent pressure response and fracture patterns simulated under homogeneous isotropic conditions with
available analytical solutions. Then, the effect of rock mass discontinuities on the hydraulic fracturing process
is investigated for several joint geometries.

1 INTRODUCTION (FEM) in describing the continuum behaviour of the

intact material, with the ability of the Discrete Element
In rock engineering, the term hydraulic fracturing (HF) Method (DEM) to capture the interaction between
refers to the injection of fluid into a sealed-off bore- fractured discrete blocks. Furthermore, FDEM cap-
hole interval to induce and propagate fractures into tures the transition from continuum to discontinuum
the surrounding rock mass. Typical applications of this by explicitly considering fracture and fragmentation
technique include the estimation of the stress state at processes. In this work, homogeneous isotropic mod-
depth as well as the stimulation of unconventional els were first considered in order to validate the numer-
reservoirs. Conventional analyses of HF operations ical results with respect to analytical solutions. Then,
are based upon the solution of the stress distribu- the numerical simulations investigated if and how a
tion around a circular borehole in a homogeneous, preexisting fracture may influence the fracture trajec-
isotropic, elastic material subjected to far-field com- tory and the injection pressure response. Input data and
pressive stresses (e.g., Hubbert and Willis (1957)). stress conditions from the in-situ stress measurement
Considering these simple and idealistic assumptions, campaigns in the geothermal borehole Schlattingen-1
the accuracy of analyses carried out by this approach (Switzerland) (Klee, 2012) and the borehole BDS-5
may be debatable. In particular, the rock mass is often (Derriere, Mont Terri, Switzerland) (Klee, 2011) were
characterized by the presence of preexisting disconti- used as reference.
nuities (e.g., fractures, faults) or layering (e.g., bedding
planes) having variable orientations with respect to the
principal stress directions. In this case, the fundamen- 2 FUNDAMENTAL PRINCIPLES OF THE
tal assumptions of the analytical approaches, namely FDEM-HF CODE
continuity, homogeneity and isotropy, are clearly not
satisfied. The FDEM-HF simulation tool consists of three main
In this context, the scope of this study was to eval- computational modules exchanging information at
uate a hybrid finite-discrete element (FDEM) code every time step:
enhanced with HF capabilities as a tool to simu- 1. the FDEM geomechanics solver (Munjiza, 2004;
late fluid-pressure-driven fracturing in jointed rock Mahabadi et al., 2012) captures the mechanical
masses. FDEM is a numerical technique which com- response of the rock mass (i.e., deformation and
bines the advantages of the Finite Element Method fracturing);

1427
 are interspersed across the edges of all triangular
element pairs (Fig. 1). Since no adaptive remeshing
is performed as the simulation progresses, potential
fracture trajectories are restricted to the existing mesh
topology. Therefore, to minimize the bias induced
on the model response, sufficiently refined unstruc-
tured meshes should be used. In the crack element
Figure 1. Representation of solid material in FDEM. model, the bonding stresses transferred by the mate-
rial are functions of the displacement discontinuity
across the crack elements according to the cohesive
laws illustrated in Figure 2.

2.1.3 Rock joint model

Rock discontinuities (i.e., either preexisting or newly
created fractures) are treated by a rock joint model,
computing the contact forces between all pairs of tri-
angular elements that overlap in space. Two types of
forces are applied to the elements of each contacting
pair: (i) repulsive forces and (ii) frictional forces. The
repulsive forces are calculated using a penalty func-
tion. Contacting pairs tend to penetrate into each other,
generating distributed contact forces, which depend on
the shape and size of the overlap between the two bod-
ies and the value of a stiffness parameter, namely the
normal penalty coefficient, pn . In the tangential direc-
Figure 2. Constitutive behaviour at the contact between tion, the frictional forces between contacting couples
bonded triangular element pairs defined in terms of nor- are calculated using Coulomb’s law of friction. Rock
mal and shear bonding stress, σ and τ, versus crack relative joints can be directly incorporated into the geome-
displacement, o and s (i.e., opening and sliding). chanics solver by aligning the mesh topology with
preexisting cracks at the time of mesh generation. In
this study, rock joints were modeled as purely fric-
2. the cavity volume calculator dynamically tracks the tional surfaces. If at any time during the simulation
evolution of wet fractures within the model and a joint intersects a fluid-driven fracture, the result-
computes variations in cavity volume due to rock ing fluid percolation is accounted by having the fluid
elastic deformation and fracturing as well as fluid pressure to be applied to the entire newly-connected
compressibility; discontinuity.
3. the pump model computes the fluid pressure based
on the injection flow rate and cavity volume.
2.2 Fluid injection and pressure-driven
fracturing
2.1 Geomechanics solver
Fluid injection and pressure-driven fracture propaga-
In the geomechanics solver the solid rock is discretized tion are captured by a simplified approach based on
as a mesh consisting of nodes and triangular elements the principle of mass conservation for a compressible
(Fig. 1). An explicit second-order finite-difference fluid injected into a deformable solid. The model is
integration scheme is applied to solve the equations hydro-mechanically coupled exclusively in the sense
of motion for the discretized system and to update the that variations in cavity volume, due to either rock
nodal coordinates at each simulation time step. elastic deformation or fracturing, affect the pressure
of the compressible fluid, which, in turn, affects rock
2.1.1 Elastic deformation deformation and failure. On the other hand, the actual
The elastic deformation of the solid rock is modeled fluid flow is neglected.
according to the continuum theory of linear elastic-
ity using constant-strain triangular elements (Munjiza, 2.2.1 Cavity volume calculator
2004). Under the assumptions of isotropic behaviour, Prior to the start of the HF simulation, the boundary
the elastic response is fully characterized by the of the initial pressurized surface (typically coincident
Young’s modulus, E, and Poisson’s ratio, ν. with the borehole perimeter) is specified by labeling
the element edges adjacent to the fluid (Fig. 3). At each
2.1.2 Rock fracture model subsequent time step, the cavity volume is computed
Rock fracturing is modeled using a cohesive-zone based on the fracture topology and interconnectivity.
approach. With this technique, crack nucleation and As the simulation progresses and new fractures nucle-
growth is captured by dedicated four-noded interface ate, new wet edges are dynamically tracked based on
elements (referred herein to as crack elements), that their connectivity with the initial pressurized surface.

1428
Figure 3. Conceptual diagram of the hydraulic fracture
propagation process in a FDEM-HF model. For illustration
purposes, only selected triangular elements are shown.
Figure 4. Geometry and boundary conditions of the HF
models. D is the diameter of the borehole.
The cavity volume calculated at each time step is fed
into the pump model.
Table 1. Rock input parameters.
2.2.2 Pump model
Parameter Schla BDS
The pump model is used to determine the fluid pres-
sure, P, to be applied to all edges marked as wet. In Bulk density, ρ (kg/m3 ) 2500 2500
general, the fluid pressure depends on the input flow Young’s modulus, E (GPa) 30 0.25
rate, Q, and responds to variations in cavity volume, Poisson’s ratio, ν 0.27 0.25
V . At every time step, the fluid mass, m, is integrated Damping coefficient, µ (kg/m · s · 105 ) 10.4 5.6
from the flow rate (specified as a time-varying bound- Tensile strength, ft (MPa) 4 3.5
ary condition). Then the pressure, P, is calculated Cohesion, c (MPa) 20 20
from this mass and the cavity volume as Mode I fracture energy, GIc (N/m) 0 0 or 1
Mode II fracture energy, GIIC (N/m) 50 50
Friction angle, φ (◦ ) 38 38
Normal contact penalty, pn (GPa · m) 300 200
Tangential contact penalty, pt (GPa/m) 30 20
Fracture penalty, pf (GPa) 150 100
where β is the bulk modulus of the fluid, which has
pressure P0 at density ρ0 . Finally, based on the length
and orientation of each wet element edge, the fluid
pressure is converted into equivalent nodal forces.
graded towards the external boundaries, where an ele-
ment size equal to 0.2 m was used for both models.
All models were meshed using the Delaunay triangu-
3 HOMOGENEOUS ISOTROPIC lation scheme available in the program Gmsh v. 2.8.3
MODELS (Geuzaine and Remacle, 2009). The input parameters
adopted for the host rock are reported in Table 1.
3.1 Model description
3.1.1 Geometry, mesh and input parameters 3.1.2 In-situ stress field, borehole excavation and
The model geometry for the FDEM-HF analyses con- fluid injection
sisted of a circular cross-section of a borehole placed The models were preliminarily initialized with the
at the centre of a 5 m × 5 m square domain (Fig. 4). in-situ stress fields reported in Table 2. The Schla-a
The borehole diameter, D, was equal to 0.16 and condition corresponds to the stress field measured in
0.10 m for the Schlattingen-1 (Schla) and BDS-5 the Schlattingen-1 borehole at a depth of 1,122 m while
(BDS) models, respectively (Klee, 2011, 2012). To the BDS-a stresses refer to the values measured in the
maximize the model resolution in the borehole near- BDS-5 borehole at a depth of 156 m. The Schla-b and
field while keeping the run times within practical BDS-b are the associated isotropic stress models. Sub-
limits, a mesh refinement zone, with average element sequently, a zero-displacement condition was applied
h = D/27 was adopted around the borehole boundary. to the far-field boundary of the model and the borehole
From the refinement zone the element size was then excavated.

1429
Table 2. In-situ stress conditions for the model.

Model σH (MPa) σh (MPa)

Schla-a 20.8 17.1

Schla-b 17.1 17.1
BDS-a 6.8 4.6
BDS-b 4.6 4.6

Table 3. Comparison of theoretical and simulated break-

down pressure values.

Pb theoretical Pb simulated
Model (MPa) (MPa)

Schla-a 34.5 34.9

Schla-b 38.2 37.4
Figure 5. Injection pressure as a function of time for the
BDS-a 11.0 11.0
Schla and BDS models with GIc = 0 N/m.
BDS-b 12.7 12.6

was recorded in response to an increase of cavity vol-

ume induced by crack growth. Finally, the injection
In the last stage of the simulation, the fluid was pressure value tended to reach a steady-state value,
injected into the borehole at a constant volumetric Psi , roughly corresponding to the minimum principal
flow rate Q = 50 l/s. Although this injection rate is in-situ stress.
significantly faster than that typically used in actual It is noteworthy that an actual comparison between
experiments (i.e., 1 to 5 l/min), an injection rate sensi- simulated and analytical pressure values is only mean-
tivity analysis indicated that the simulated breakdown ingful if the Mode I fracture energy was assumed
and shut-in pressure values approached constant val- equal to zero. That is, the post-peak strain softening
ues for flow rates smaller than approximately 100 l/s. behaviour of the crack elements (Fig. 2) is neglected
Initial fluid pressure, fluid density and bulk modu- and a crack nucleates as soon as the tensile strength,
lus were set equal to 0 MPa, 1,000 kg/m3 and 2.2 GPa, ft , is reached (i.e., ft = T ), thus resulting in an elastic-
respectively. brittle response. On the other hand, if a non-nil GIc
value is assumed, the tensile strength of the rock effec-
tively increases together with the Pb and Psi values.
3.2 Results and discussion For the case of GIc = 1 N/m, Pb showed increases of
As summarized in Table 3, a good agreement was 0.7 and 0.4 MPa for the BDS-a and BDS-b models,
obtained between the simulated values of breakdown respectively. Under anisotropic in-situ stress condi-
pressure, Pb , and the analytical values proposed by tions, the hydraulic fractures started to develop, as
Hubbert and Willis (1957): expected, approximately at the intersection of the σH
direction with the borehole boundary (Fig. 6). Crack
branching could be observed early in the simulation
in close proximity to the borehole boundary. In con-
trast, under isotropic in-situ stress conditions, a radial
where T is the tensile strength of the rock. Notice that
pattern of multiple fractures initiated and propagated
in the original formulation an additional term, Pw , cor-
from the borehole.
responding to the rock pore pressure is added to the
right member of Equation (2). In this case, the term
was neglected as no poromechanical formulation is
currently available in FDEM code. 4 MODELS WITH JOINTS
For each case, the simulated Pb value was derived
from the plots of the injection pressure, P, as func- The effect of preexisting discontinuities on the HF pro-
tion of the simulation time, T , depicted in Fig. 5. The cess and injection pressure response was investigated
simulated injection curves correctly reproduced the for the BDS model for the geometrical configurations
typical behaviour observed in-situ. Initially, a linear depicted in Figure 7. In general, in the presence of
increase of injection pressure was captured owing to preexisting discontinuities intersecting the pressurized
the coupled elastic deformation of the fluid inside the borehole, a sensible reduction of breakdown pres-
borehole and the surrounding solid rock. The appear- sure was observed with respect to the corresponding
ance of the first fluid-pressure-induced tension crack homogenous and isotropic model (Fig. 8). Also, preex-
was associated with a maximum in the injection pres- isting discontinuities sensibly influenced the simulated
sure curve (i.e., Pb ). Subsequently, a drop in pressure fracture trajectories, as depicted in Figure 9.

1430
 Figure 8. Injection pressure as function of time for the
BDS models with GIc = 1 N/m.

discontinuities due to the stress conditions more favor-

able for tensile failure at a distance from the borehole.
As more fluid was pumped into the borehole, the
hydraulic fractures tended to re-align themselves in the
direction parallel to σH thus resulting in a characteristic
curved fracture trajectory. With respect to the homo-
geneous case, the breakdown pressure decreased from
11.7 MPa to 7.7 MPa. For the case of flaw geometry
F3, the injection pressure response was analogous of
the F1 case. As depicted in Fig. 16(b), fractures nucle-
ated at an angle from the flaw tip and then tended to
Figure 6. Simulated fracture trajectories and evolution of realign themselves in the direction parallel to σH , thus
σ1 contour for the BDS-a and BDS-b models (GIc = 1 N/m). resembling the wing-cracking phenomenon observed
during compression tests on flawed rock specimens
(e.g., Hoek and Bieniawski, 1965). For the case of flaw
geometry F4, no reduction of Pb was observed com-
pared with the homogeneous case as the two flaws were
not hydraulically connected to the pressurized bore-
hole.Also, the post-peak pressure response was similar
to the homogeneous case with the only exception that a
secondary(smaller) spike in injection pressure tended
to appear soon after the drop in pressure from the Pb
value.
Figure 7. Geometries used to investigate the influence of
rock mass fabric on the simulated HF process.

5 CONCLUDING REMARKS

For the case of flaw geometry F1, hydraulic frac- A preliminary evaluation of a hybrid continuum-
tures nucleated at the tip of the discontinuity and discontinuum numerical code to simulate hydraulic
propagated along a straight line for both in-situ stress fracture propagation was carried out. The FDEM code
fields. Since the flaw extended partially out of the was equipped with a computational module to sim-
borehole-induced stress concentration area, the injec- ulate fluid-pressure-driven fracturing. The approach
tion pressure Pb required to open a fracture (5.5 MPa) was first validated by reproducing injection pressure
became equal to the value of Psi (5.1 MPa), which responses as well as fracture patterns expected under
in turns approximated σh (4.8 MPa). Therefore, the homogeneous and isotropic conditions. Preliminary
pressure injection curve did not exhibit any drop in investigations on the effect of preexisting discontinu-
pressure but showed a bilinear behaviour (i.e., linear ities indicated that, unlike conventional analytical and
elastic increase followed by a constant pressure value). continuum-based numerical methods, FDEM has the
For the case of flaw geometry F2, hydraulic frac- capability to realistically account for the influence of
tures nucleated again at the tip of the preexisting rock mass fabric on the hydraulic fracture process.

1431
Figure 9. Simulated fracture trajectories and evolution of σ1 contour for the BDS models with joints (GIc = 1 N/m).

REFERENCES Terri). Technical Note 2011-45. Mont Terri Consortium.

St-Ursanne, Switzerland.
Geuzaine, C., Remacle, J.F., 2009. Gmsh: a three- Klee, G., 2012. Geothermal borehole Schlattingen 1.
dimensional finite element mesh generator with built-in Hydraulic fracturing stress measurements. Project Report
pre- and post-processing facilities. International Journal 12-08. NAGRA. Wettingen, Switzerland.
for Numerical Methods in Engineering 79: 1309–1331. Mahabadi, O.K., Lisjak, A., Grasselli, G., Munjiza, A., 2012.
Hoek, E., Bieniawski, Z., 1965. Brittle failure propagation Y-Geo: a new combined finite-discrete element numer-
in rock compression. International Journal of Fracture ical code for geomechanical applications. International
Mechanics 1: 137–155. Journal of Geomechanics 12: 676–688.
Hubbert, M.K., Willis, D.G., 1957. Mechanics of hydraulic Munjiza, A., 2004. The combined finite-discrete element
fracturing. Technical Report. US Geological Survey. method. John Wiley & Sons Ltd.
Klee, G., 2011. DS Experiment, Phase 15: Hydraulic fractur-
ing measurements in BDS-5 borehole at Derriere (Mont

1432

View publication stats