University of Southern California Dissertations and Theses

An extended finite element method based modeling of hydraulic fracturing

(USC Thesis Other)

An extended finite element method based modeling of hydraulic fracturing

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

AN EXTENDED FINITE ELEMENT METHOD BASED MODELING OF
HYDRAULIC FRACTURING
by
Arman Khodabakhshnejad

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
in Partial Fulfillment of the
Requirement for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)

August 2016

Copyright 2016 Arman Khodabakhshnejad
i

Dedication

To my parents Mashalah and Sedigheh
My wife Rana
My sisters Azadeh, and Nooshin
Without whom none of my success would be possible
ii

Abstract
Maximizing hydrocarbon production from unconventional reservoirs requires the implementation
of best practices in order to create a connected network of hydraulic and natural fractures. The
desired result will be achieved if the operational hydraulic fracturing parameters are optimized.
Determining the appropriate number of fracturing stages needed for peak performance demands a
well calibrated model. Such a model allows for a more accurate prediction of fracture distribution.
To this end, the industry currently employs several forward simulation tools. However, owing to
the complexity of this problem, and the simplified assumptions made, existing methods often
provide a poor representation of reality, leading to inaccurate predictions. To address these
challenges, this work presents a novel, coupled approach to model hydraulically fractured
subsurface cracks, providing a dynamic evolution of the subsurface fracture network.
In this study, I introduce a new algorithm to model the initiation and propagation of
hydraulic fracturing. It integrates rock deformation and fluid flow using a coupling approach. The
algorithm enables us to accurately simulate fluid pressure inside the fracture and rock matrix and
to calculate the in situ stress distribution within the reservoir.
An in-house fluid flow simulator is developed based on a modified diffusivity equation
that takes the effect of fracture width changes into account. I calculate fluid pressure using Finite
Difference Method (FDM), and export it to the geo-mechanical simulator as a load. To calculate
the resulting stress distribution inside the model and to predict fracture initiation and propagation,
I apply the Extended Finite Element Method (XFEM). The workflow involves integrating the fluid
flow and rock deformation equations; the problem is solved by an iterative implicit scheme. I
iii

validate the performance of the algorithm, applying it to both a 2D laboratory scale and 3D field
scale modeled fracturing data.
The effect of numerical parameters on the fracture behavior in 2D modeling was examined;
the algorithm was used to investigate stress shadowing and the interaction between stimulated and
natural fractures. To increase the computation speed of the simulation, I apply uniform pressure
distribution to the entire fracture network. The algorithm later is used for performing sensitivity
analyses on a set of numerical parameters. Mesh resolution, plane stress and plane strain condition,
model size, and boundary condition are analyzed for a single fracture model. It became clear that
the aforementioned factors significantly impact the model response to the applied stresses. I
conduct sensitivity analysis for the resulting stress shadow as a function of the distance between
the fractures. In addition, the interaction between a hydraulic fracture and a natural fracture was
evaluated. The modeling results suggest that the hydraulic fracture can activate natural fractures
both remotely or after their intersection.
I also test the algorithm to assess its capabilities in predicting the hydraulic fracturing
parameters in a three dimensional model. The Griffith-Sneddon analytical solution of penny-
shaped fracture was used to investigate the accuracy of the model. I show that the algorithm can
accurately predict fracture width as well as fracture volume. Furthermore, the modeling results of
hydraulic fracturing propagation compare to the Geertsma-de Klerk analytical solution. The results
show some discrepancy between analytical and numerical solutions. Mesh size is the major factor
that can control the deviation of the numerical model results from those derived from the analytical
models.
iv

In field size reservoir modeling of a fracture, the algorithm reveals the mechanisms that
could not be modeled through conventional modeling of hydraulic fracturing. For example, the
approach accurately models pressure disturbance in the reservoir during fracturing.
Reservoir pressure modeling suggests that pressure builds up at the central elements of the
fracture, and consequently a high leak-off rate is expected. In contrast, pressure in the reservoir
and near the crack tip drops due to the creation of a suction zone ahead of fluid front and behind
the crack tip. Unlike conventional approaches, the results of this work demonstrate that my
algorithm successfully models hydraulic fracturing propagation in a 3D reservoir. Additionally,
my new approach is capable of modeling the interaction between hydraulic and natural fractures
in order to generate more reliable predictions.

v

Acknowledgement

First and foremost, I would like to express my sincere gratitude to my advisor, Prof. Aminzadeh
for his guidance and utmost support during my PhD studies at USC. His profound knowledge of
petroleum engineering and geoscience as well as his responsibility and commitment, afforded me
with an exceptional opportunity to nurture my abilities in different areas of science. Every
conversation with Prof. Aminzadeh provided me with new ideas and insight. I started my PhD as
a reservoir engineer and finishing it as a geomechanical engineer. He cultivated a constructive
environment in which new ideas could grow. He benefited me not only for my PhD but also for
my future career. I am incredibly proud to be Professor Aminzadeh’s student.
I would like to especially thank Prof. Iraj Ershaghi, Director of the Petroleum Engineering
Program, for his insightful advice and generous support during my PhD. I want to express my
appreciation to Prof. Iraj Ershaghi and to Prof. Charles Sammis for serving on my PhD and
qualifying examination committee. I am thankful for all your contributions, and the time and effort
you spent on reviewing my thesis. Thanks also to Prof. Roger Ghanem and Prof. Behnam Jafarpour
for serving as committee members for my qualifying exam.
I am grateful to Dr. Ahmed Ouenes, for whom I had the pleasure to work in my summer
internship in FracGeo. I have gained a wealth of experience and knowledge, and part of it is
presented in the thesis. I am also thankful for the valuable technical feedback from Dr. Yesser
HajNasser.
vi

The partial financial support was provided by the Reservoir Monitoring Consortium
(RMC) and Induced Seismicity Consortium (ISC) is greatly acknowledged. A specials thanks to
the members for providing me with valuable feedback.
There are many friends and colleagues who have helped me to accomplish my PhD, and
without whom I would not be here now. My especial thanks go to Mohammad Evazi, Shahram
Farhadi Nia, Mohammad Javahery, Hasan Shojaei, and Mohsen Taghavifar for being such
amazing friends. It is my privilege to have you as friends. I also want to mention that the
friendships of Hamidreza Jahangiri, Siavsh Hakim Elahi, Azerang Gol Mohammadi, Magdalene
Ante, Mehran Hosseini, Marjan Sherafati, Mehdi Shams, Yahya Taghavifar, Ahmed Bubashait,
and Robert Walker have been invaluable. Thank you for providing such a great environment.
Finally, I want to thank my parents, Mashalah Khodabakhshnejad, and Sedigheh
Daneshgar for their unconditional support and love throughout my entire life. My wife, Rana
Ghajar, has given endless love, patience, and support. I am very thankful; I could not have done
my PhD without her.

vii

Table of Contents

Abstract ........................................................................................................................................... ii
Table of Contents .......................................................................................................................... vii
List of Tables ................................................................................................................................. xi
List of Figures ............................................................................................................................... xii
Chapter 1 : INTRODUCTION........................................................................................................ 1
1.1 Background ............................................................................................................... 1
1.2 Objective of the Study .............................................................................................. 4
1.3 Outline....................................................................................................................... 5
Chapter 2 : BACKGROUND AND LITERATURE REVIEW ...................................................... 7
2.1 Introduction ............................................................................................................... 7
2.2 Rock Constitutive Models......................................................................................... 7
2.3 Fracture Models in Rocks ....................................................................................... 10
2.3.1 Energy Based Method ...................................................................................... 15
2.3.2 Stress Based Method ........................................................................................ 16
2.3.3 Displacement Based Method ............................................................................ 16
2.4 Cohesive Zone Model ............................................................................................. 17
2.5 Fracture Propagation in Rock ................................................................................. 22
2.6 Fracture Propagation in Hydraulic Fracturing: a Theoretical Approach ................ 30
2.6.1 Perkins-Kern-Nordgren (PKN) ........................................................................ 30
2.6.2 Khristianovitch-Zheltov-Geertsma-De Klerk (KGD) ...................................... 31
2.6.3 Radial Model .................................................................................................... 32
2.7 Fracture Propagation in Hydraulic Fracturing: A Numerical Approach ................ 33
viii

2.7.1 Finite Difference Method ................................................................................. 34
2.7.2 Discrete Element Method ................................................................................. 35
2.7.3 Boundary Element Method .............................................................................. 36
2.7.4 Finite Element Method ..................................................................................... 38
2.7.5 Extended Finite Element Method ..................................................................... 41
2.8 Concluding Remarks ............................................................................................... 45
Chapter 3 : DESCRIPTION OF HYDRAULIC FRACTURING SIMULATION ALGORITHM
....................................................................................................................................................... 47
3.1 Formulation of Fluid Flow Equation ...................................................................... 47
3.2 Numerical Analysis of Fluid Flow Equations ......................................................... 50
3.3 Stress Modeling ...................................................................................................... 52
3.4 Fracture Modelling.................................................................................................. 53
3.4.1 Rock Damage Model ........................................................................................ 53
3.4.2 Cohesive Segment Method and Phantom Nodes .............................................. 54
3.5 Algorithm Development and Numerical Simulation Scheme ................................ 56
3.6 Concluding Remarks ............................................................................................... 59
Chapter 4 : TWO DIMENSIONAL MODELING OF FRACTURING IN LABORATORY
SCALE .......................................................................................................................................... 61
4.1 Introduction ............................................................................................................. 61
4.2 Single Fracture Propagation .................................................................................... 62
4.2.1 Effect of Resolution .......................................................................................... 68
4.2.2 Effect of Plane Stress and Plane Strain Condition .......................................... 72
4.2.3 Effect of Type of Analysis ................................................................................. 74
4.2.4 Effect of Model Shape ...................................................................................... 79
4.2.5 Effect of Initial Fracture Location ................................................................... 83
4.3 Multiple Crack Propagation .................................................................................... 92
ix

4.3.1 Effect of Fracture Spacing In Rectangular Plate ............................................ 93
4.3.2 Effect of Fracture Spacing In Square Plate ..................................................... 97
4.3.3 Hydraulic and Natural Fracture Interaction ................................................. 103
4.4 Conclusion Remarks ............................................................................................. 112
Chapter 5 : THREE DIMENSIONAL SIMULATION OF HYDRAULIC FRACTURING IN
FIELD SCALE............................................................................................................................ 114
5.1 Introduction ........................................................................................................... 114
5.2 Description of Model ............................................................................................ 119
5.3 Parameter Sensitivity Analysis ............................................................................. 121
5.3.1 Sensitivity Analysis on Rock Strength and Stress Field ................................. 122
5.3.2 Sensitivity Analysis on Reservoir Size............................................................ 123
5.4 Field Scale Numerical Analysis of Hydraulic Fracturing ..................................... 126
5.5 Model Verification ................................................................................................ 129
5.5.1 Fracture Width ............................................................................................... 130
5.5.2 Fracture Volume ............................................................................................ 132
5.5.3 Fluid Pressure ................................................................................................ 133
5.5.4 Fracture Length ............................................................................................. 135
5.6 Discussion ............................................................................................................. 136
5.6.1 Fracture Length Analysis in Hydraulic Fracturing ....................................... 137
5.6.2 Reservoir Pressure and Leak-Off................................................................... 137
5.7 Conclusion Remarks ............................................................................................. 139
Chapter 6 : SUMMARY AND FUTURE RESEARCH DIRECTION ...................................... 141
6.1 Summary ............................................................................................................... 141
6.2 Recommendation and Future Work ...................................................................... 143
Appendix A Derivation of the Fluid Flow Equation .................................................................. 146
x

Appendix B Multiple Fracture Propagation under Tensile Load................................................ 152
References ................................................................................................................................... 162

xi

List of Tables
Table 2.1: Various Cohesive Zone Models and their advantages and drawbacks ........................ 20
Table 4.1: Fracture and plate characteristics, and boundary condition ......................................... 65
Table 5.1: Reservoir, geometrical, and damage properties of the model ................................... 121
Table 5.2: Hydraulic fracturing fluid properties ......................................................................... 126

xii

List of Figures
Figure 2.1: Constitutive response of different materials ............................................................... 10
Figure 2.2 Presence of the fracture in the large scale ................................................................... 12
Figure 2.3: Study of a problem in different scales (Buehler and Theodor, 2007). ....................... 12
Figure 2.4: Sample response to the load of lateral, volumetric, and axial forces (Nicksiar and
Martin, 2012). ............................................................................................................................... 13
Figure 2.5: Crack and its process zone in front of the tip. The area shows the degraded element.
....................................................................................................................................................... 18
Figure 2.6: The traction-separation law for a material with linear softening ............................... 19
Figure 2.7: Crack propagation modes. .......................................................................................... 23
Figure 2.8: Transient from micro scale to macro scale. Discrete medium is the middle stage and
fits for the uncemented rock materials (Rojek and Oñate, 2007). ................................................ 26
Figure 2.9: Damage evolution in rock in micro scale. a) Unloaded rock, b) Rock under load
lower than failure, c) Rock break down under a high compressive stress (DiGiovanni et al, 2007).
....................................................................................................................................................... 27
Figure 2.10: Fracture growth modeling in a material using Abaqus. ........................................... 28
Figure 2.11: Fracture network pattern in mineback experiments (Warpinski and Teufel, 1987). ....
....................................................................................................................................................... 29
Figure 2.12: Fracture network models to estimate fracture propagation in the reservoir (Cipolla et
al, 2008). ....................................................................................................................................... 29
Figure 2.13: Fracture schematic of PKN model (Economides and Nolte, 2000). ........................ 31
Figure 2.14: Fracture schematic of KGD model (Economides and Nolte, 2000). ........................ 32
Figure 2.15: Fracture schematic of penny-shape model (Economides and Nolte, 2000). ............ 32
Figure 2.16: Finite Difference Method implementation in the 2-D space .................................... 35
Figure 3.1: Fluid moves inside the crack and simultaneously transfers between crack and matrix
....................................................................................................................................................... 48
Figure 3.2: Fracture and matrix discretization for fluid flow modeling ....................................... 51
Figure 3.3: Workflow for the modeling of hydraulic fracturing ................................................... 57
Figure 4.1: Plate is subjected to two tensile loads acting from the top and bottom of the plate.
Initial fracture is placed at the center left part of the plate. .......................................................... 63
xiii

Figure 4.2: Plate is subjected to the two tensile loads acting from the top and bottom of the plate.
Initial fracture is placed at the center left part of the plate. .......................................................... 64
Figure 4.3: Plate is subjected to the two compressive loads acting from the top and bottom of the
plate. Initial fracture is placed at the left part of the plate. Pressure inside the fracture increases
and applies load on the fracture walls. .......................................................................................... 64
Figure 4.4: Stress field distribution in a low resolution rectangular plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by
factors of 10,000, 1,000, 100, 10, and 2, respectively. ................................................................. 69
Figure 4.5: Stress field distribution in a high resolution rectangular plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by
factors of 10,000, 1,000, 100, 10, and 2, respectively. ................................................................. 70
Figure 4.6: The Von Mises stress distribution for 100 by 50 meshing scheme. ........................... 72
Figure 4.7: Stress field distribution in a high resolution rectangular plate in a plane strain
condition. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The
width of fractures are magnified by factors 2. .............................................................................. 73
Figure 4.8: Stress field distribution in a high resolution rectangular plate for 10MPa/100s load
rate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width
of fractures are magnified by factors of 10,000, 1,000, 100, 10, and 2, respectively. .................. 75
Figure 4.9: Stress field distribution in a rectangular plate. The fracture is magnified by a factor of
2. Legend is modified for improved representation. ..................................................................... 76
Figure 4.10: Stress field distribution in a rectangular plate for load rate of 10 MPa/10 s. The
width of fracture is magnified by a factor of 2. ............................................................................ 78
Figure 4.11: Stress field distribution in a rectangular plate for load rate of 10 MPa/1 s. The width
of fracture is magnified by a factor of 2. ...................................................................................... 78
Figure 4.12: Stress field distribution in a square plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c)
10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by factors of 10,000,
1,000, 100, 10, and 2, respectively. .............................................................................................. 81
Figure 4.13: Von Mises stress distribution before fracture, (a) the x-axis represents the true
distance from fracture tip, (b) the x-axis represents the normalized distance from fracture tip. ......
....................................................................................................................................................... 83
xiv

Figure 4.14: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid
pressure is (a) 10
6
Pa, (b) 2×10
6
Pa, (c) 3×10
6
Pa, and (d) 3.308×10
6
Pa. The width of fractures
are magnified by factors of 10,000, 100, 100, and 20, respectively. Initial fracture is placed 50
cm from the top. ............................................................................................................................ 85
Figure 4.15: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid
pressure is (a) 10
6
Pa, (b) 2×10
6
Pa, (c) 3×10
6
Pa, and (d) 3.204×10
6
Pa. The width of fractures
are magnified by factors of 10,000, 100, 100, and 20, respectively. Initial fracture is placed 34
cm from the top. ............................................................................................................................ 87
Figure 4.16: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid
pressure is (a) 10
6
Pa, (b) 2×10
6
Pa, and (c) 2.844×10
6
Pa. The width of fractures are magnified
by factors of 100, 100, and 10, respectively. Initial fracture is placed 25 cm from the top. ........ 88
Figure 4.17: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid
pressure is (a) 10
6
Pa, (b) 2×10
6
Pa, and (c) 2.27×10
6
Pa. The width of fractures are magnified by
factors of 100, 100, and 10, respectively. Initial fracture is placed 17 cm from the top. ............. 90
Figure 4.18: Fracture initiation pressure in different fluid pressures ............................................ 91
Figure 4.19: Plate failure pressure for different fluid pressures ................................................... 91
Figure 4.20: Plate is subjected to the two compressive loads acting from the top and bottom of
the plate. Two fractures are placed at the plate and are separated by a distance of h. The fluid
pressure force the fractures to open and propagate. ...................................................................... 93
Figure 4.21: Stress field distribution in a rectangular plate. Distance between two fractures is 10
cm. Fluid pressure is (a) 2.677×10
6
Pa, (b) 3.628×10
6
Pa, (c) 3.785×10
6
Pa, (d) 3.878×10
6
Pa,
and (e) 3.879×10
6
Pa. The width of fractures are magnified by factor 10. ................................... 94
Figure 4.22: Stress field distribution in a rectangular plate. Distance between two fractures is 30
cm. Fluid pressure is (a) 2.183×10
6
Pa, (b) 3.174×10
6
Pa, (c) 3.614×10
6
Pa, (d) 3.723×10
6
Pa,
and (e) 3.735×10
6
Pa. The width of fractures are magnified by factor 10. ................................... 95
Figure 4.23: Stress field distribution in a rectangular plate. Distance between two fractures is 50
cm. Fluid pressure is (a) 2.130×10
6
Pa, (b) 3.146×10
6
Pa, (c) 3.3127×10
6
Pa, (d) 3.3128×10
6
Pa,
and (e) 3.3129×10
6
Pa. The width of fractures are magnified by factor 10. ................................. 97
Figure 4.24: Plate is subjected to the two compressive loads acting from the top and bottom of
the plate. Initial fractures are placed at the left part of the plate. Distance between the fractures is
shown by h. ................................................................................................................................... 98
xv

Figure 4.25: Stress field distribution in a square plate. Distance between two fractures is 5 cm.
Fluid pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.728×10
6
Pa, and (e)
3.744×10
6
Pa. The width of fractures at (a) to (d) are magnified by factors 10 and 2 for figure (e).
....................................................................................................................................................... 99
Figure 4.26: Stress field distribution in a square plate. Distance between two fractures is 10 cm.
Fluid pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.480×10
6
Pa, and (e)
3.734×10
6
Pa. The width of fractures are magnified by factors 10. ........................................... 101
Figure 4.27: Stress field distribution in a square plate. Distance between two fractures is 20 cm.
Fluid pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.551×10
6
Pa, and (e)
3.583×10
6
Pa. The width of fractures are magnified by factors 10. ........................................... 102
Figure 4.28: Fracture propagation in a plate in the direction of existing fracture ...................... 103
Figure 4.29: Stress field distribution in a rectangular plate. Distance between two fractures is
33.75 cm. Fluid pressure is (a) 2.34×10
6
Pa, (b) 2.35×10
6
Pa, (c) 2.37×10
6
Pa, (d) 2.381×10
6
Pa,
(e) 2.382×10
6
Pa, (f) 2.383×10
6
Pa, (g) 2.384×10
6
Pa, (h) 2.385×10
6
Pa, and (i) 2.390×10
6
Pa.
The width of fractures are magnified by factors 10. ................................................................... 104
Figure 4.30: Shear stress around the crack tips during hydraulic fracture growth. .................... 105
Figure 4.31: Plate is subjected to the four compressive loads acting from the four edges. Distance
between two fractures is shown by d. ......................................................................................... 107
Figure 4.32: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two
fractures is 25 cm. The width of fractures are magnified by factor 10. The pressure inside
fracture is 2.7 MPa. ..................................................................................................................... 108
Figure 4.33: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two
fractures is 37.5 cm. The width of fractures are magnified by factor 10. The pressure inside
fracture is 2.7 MPa. ..................................................................................................................... 110
Figure 4.34: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two
fractures is 62.5 cm. The width of fractures are magnified by factor 10. The pressure inside
fracture is 2.7 MPa. ..................................................................................................................... 111
Figure 5.1: Schematic of the reservoir constructed for modeling of hydraulic fracturing. ........ 119
Figure 5.2: Effect of rock tensile strength on fracture initiation pressure .................................. 122
Figure 5.3: Effect of reservoir width on fracture initiation pressure .......................................... 124
Figure 5.4: Effect of overlaying and underlying sand layers on fracture initiation pressure ...... 125
xvi

Figure 5.5: Initial fracture is place at the center of reservoir. The center of fracture and reservoir
coincide. ...................................................................................................................................... 126
Figure 5.6: Fracture propagation inside the reservoir after, (a) 9.30 s, (b) 19.48 s, (c) 26.43 s, (d)
27.17 s of fluid injection. Red zone represents the fracture. ....................................................... 128
Figure 5.7: Fracture propagation in time. Y-axis shows the equivalent radius of the fracture. .. 129
Figure 5.8: Analytical and numerical calculation of fracture width for a penny-shaped fracture.
..................................................................................................................................................... 131
Figure 5.9: Ratio of numerical to analytical solution in a penny-shaped fracture ...................... 132
Figure 5.10: Analytical and numerical representation of the fracture volume ........................... 133
Figure 5.11: Analytical and numerical calculation of fluid pressure .......................................... 134
Figure 5.12: Analytical and numerical calculation of the fracture radius in the early time of the
hydraulic fracturing ..................................................................................................................... 135
Figure 5.13: Pressure change in the reservoir elements adjacent to the fracture wall during
fracture propagation. The scale covers -400 psi to 100 psi change of the fluid pressure compared
to initial 3000 psi. ....................................................................................................................... 138

1

Chapter 1 : INTRODUCTION
1.1 Background
Injecting high pressure fluid into reservoir rock creates a complex network of connected fractures.
That operation is called hydraulic fracturing in the petroleum engineering industry, and is known
as fracking in environmental engineering. Various parameters influence fracture network
characteristics. These parameters can be part of an operational plan or they can be related to the
reservoir properties. Accurate characterization of the fracture network depends on understanding
reservoir properties, as well as the physics of the initiation and propagation of the fractures.
Reservoir properties can be obtained through reservoir characterization and monitoring tools that
are based upon the integration of geophysical, geological, and petro-physical data, as well as other
field measurements. The fundamental physics governing the creation and propagation of fractures,
and resulting from the stimulation process (induced fractures) could be determined through
experimental and theoretical study of fractures, rock, and fluid properties. Fracture mechanics is a
branch of mechanics and it is introduced to deal with such problems for different circumstances of
crack initiation and propagation. However, fracture study in hydraulic fracturing is usually
different from normal fracture mechanics problems, as many other factors contribute to the
fracturing process. Hydraulic fracturing progresses by the continuous hydraulic force provided by
the fluid injection. Therefore, the properties of fracture fluid should be incorporated in any
inclusive modeling of fracture network propagation. Kinematic, thermodynamic, and rheological
properties of fluid should be included in the hydraulic fracturing simulator in order to offer an
accurate predictive model.
Fractures should remain open when fluid pressure declines; towards that end, sand
particles, called proppant, are added to the fluid. Effects of the proppant in any post fracture
2

network modeling should be incorporated in order to calculate operational fractured volume. Other
mechanisms, such as proppant bridging, could limit the fluid flow in many fractures which would
impact the fracture network propagation.
An accurate hydraulic fracture network characterization and modeling depend on many
parameters. Incorporating all of these parameters in a single model is nearly impossible. Therefore,
a reduced model should be developed that takes into account the most influential components.
Many models have been developed that predict fracture growth in hydraulic fracturing
operations. Early attempts were mostly based on simplifying assumptions, such as growth of planar
fracture rather than arbitrary fracture structure (Khristianovitch et al, 1955; Geertsma and De
Klerk, 1969; Perkins and Kern, 1961; Nordgren, 1972; Daneshy, 1973; Hagoort et al, 1980). Since
computer simulation had not yet been established, most of those models were based on analytical
solutions of elasticity equations. Among them, PKN and GDK were widely used to predict fracture
growth. As advances grew in computer processing capacity, numerical models were developed in
order to capture different aspects of fracture growth. de Pater et al (2005) conducted an experiment
to model fracture growth in the rock. They also used Distinct Element Method (DEM) to model
the fracture network. Predictions were improved by adjusting the model parameters. Mofazzal
Hossain and Rahman (2008) used Boundary Element Method (BEM) to model fracture growth.
They assumed formation rock as an impermeable medium and therefore ignored the fluid transfer
between fracture and matrix. Weng et al (2011) solved a fully coupled fracture network
propagation in the presence of natural fractures. They assumed that the fractures were vertical and
ignored nonplanar behavior of fractures. They used criteria based models to include both induced
and natural fracture interaction. The effect of geological discontinuity in hydraulic fracturing
growth was studied by Zhang et al (2007). Finite Difference Method (FDM) was used to model
3

fluid flow, and six finite elements to model fracture growth behavior of the fracture tips. The model
was limited to one discontinuity, and the presence of a fracture network was disregarded. Riahi
and Damjanac (2013) studied the interaction of hydraulic and natural fractures using DEM, and
described pre-existing natural fractures using Discrete Fracture Network (DFN). They assumed
that natural fractures are elastic and deform during hydraulic fracture propagation. Olson (2008)
developed a simplified multiple hydraulic fracturing model to study hydraulic fracture propagation
in the presence of pre-existing natural fractures. His model was based on Displacement
Discontinuity (DD), Boundary Element Method (BEM) and was able to model rock deformation;
however, fluid flow was ignored to avoid undue complexity. Khodabakhshnejad and Aminzadeh
(2014) modeled fluid pressure change in hydraulic fracturing using a coupled approach of fluid
flow and rock deformation. According to the model, fluid pressure change in hydraulic fracturing
could be significant. Also, they reported that the 2D modeling might not be a reliable approach for
hydraulic fracturing simulation.
Reviewing the studies in the hydraulic fracturing modeling area reveals the need for
development of a new approach. Fluid flow in the fracture and matrix, fluid transfer between
matrix and fracture, stress field at the reservoir, heterogeneity of the formation rock, fracture
initiation and growth, and other important mechanisms and factors should be incorporated in a
comprehensive model to enable engineers to predict fracture propagation and to estimate operation
efficiency more accurately. It should be feasible to model fracture growth and behavior under
differing circumstances of operational and reservoir conditions. The model should be able to
address the following matters:
 Fracture and matrix fluid transfer
4

 Change of stress field in time
 Accurate rock model
 Reliable fracture model
 Presence of natural fractures
 Non planar fracture growth
 Include remote interaction of natural and induced fractures
 Numerical stability
Many studies addressed at least one of those items and have offered a detailed explanation
of each factor. However, real world hydraulic fracturing modeling needs a combination of all the
above-mentioned aspects. Therefore, a more thoroughly integrated approach would avail itself of
all the advances in different areas.
1.2 Objective of the Study
Many hydraulic fracturing models have been developed to tackle this problem from one specific
viewpoint. The complex nature of the reservoir, underdeveloped characterization tools, inaccurate
physical models, as well as low computational spaces necessitate the consideration of simplifying
assumptions to make the problem solvable. However, many of these assumptions result in the
creation of an unrealistic picture of the process. The lack of a comprehensive model that considers
the effects of different mechanisms such as leak-off, fracture growth, and stress shadow, motived
us to develop a procedure that encompasses all the important underlying processes in fracking. We
can build a platform to incorporate a representative rock model of the reservoir, if provided. Also,
5

we employ a fluid simulator, and fracture modeling software simultaneously. Each program is
developed separately, then tested and validated through the sensitivity analysis for a simple model.
The fluid simulator incorporates fundamental formulation of fluid flow in porous media to predict
the fluid flow transfer between the fracture and matrix. The fracture modeling software employs
eXtended Finite Element Method (XFEM) to assess fracture behavior for the simplified model.
We established a platform to integrate several numerical and analytical models that were
created for hydraulic fracturing modeling. MATLAB serves as the main platform to link Abaqus
and other MATLAB and Python functions in the program. As mentioned earlier, eXtended Finite
Element Method (XFEM) was developed in Abaqus for fracture modeling. We have used that
numerical modeling scheme to calculate fracture locations and parameters. Fluid flow in porous
media, and fracture is incorporated within the model using a coupling approach. The objective of
the project is to model fracture initiation and growth realistically. The model is capable of
incorporating the elasto-plastic behavior of rock. No prior assumption is made regarding the
fracture shape and size; fractures can grow in different directions and sizes. The coupling approach
enables us to model fluid and rock behavior simultaneously, and to determine the major factors in
hydraulic fracturing simulation. In addition, we have tested the effect of multi-fracturing
propagation in a 2-D plate, which could provide us with new insight into the hydraulic and natural
fracture interaction. Based on this model, a realistic geo-mechanical simulation program will be
created, which can be further modified to model fracture propagation in a realistic reservoir.
1.3 Outline
The reminder of this report is organized as follows. In Chapter Two, a review of the rock
constitutive and fracture models are presented. These models describe the rock behavior under the
stress. Recent advances in fracture zone characterization and modeling are discussed. Cohesive
6

Zone Model shed some lights on the fracture nucleation and propagation. A brief description of
the model will be provided. Different techniques are developed for hydraulic fracturing simulation.
The theory and numerical treatment of these models are discussed in this chapter. Chapter Three
explains the fundamental of the developed model. The fluid flow equations are derived and
discretized based on the Finite Difference Method. Constitutive equations for isotropic linear
elasticity are presented in this chapter. Also, the damage model that is used for process zone
characterization is explained. The algorithm to solve the equation is elaborated in this chapter.
Chapter Four describes the result of two-dimensional hydraulic fracturing simulation using the
simplified form of algorithm. The key assumption is the uniform pressure in fracture. The result
of the workflow is presented in Chapter Five. The model is described at the beginning of the
chapter and a few examples show the sensitivity of the model to different geometrical and
geomechanical parameters. Three dimensional simulation of the hydraulic fracturing is presented
to show the model capabilities. Chapter Six summarizes the thesis and discusses future work.

7

Chapter 2 : BACKGROUND AND LITERATURE REVIEW
2.1 Introduction
In this chapter I will provide the background needed to understand the importance of developing
a new hydraulic fracturing model. The problem in hydraulic fracturing simulation is twofold: (1)
to transcend the current Linear Elastic Fracture Mechanics (LEFM) method, and develop a more
realistic approach to rock behavior modeling based on experimental and field observation; and (2)
to utilize more sophisticated numerical tools to model fracture initiation and propagation in the
rock. The purpose of this research is to introduce a model capable of addressing the issues of rock
mechanics and numerical simulation.
2.2 Rock Constitutive Models
Rock behavior under differing circumstances is formulated using developed constitutive models.
These models are based on the physics of rocks -- their structure, stress condition, and so forth.
Classical constitutive law, which is based on Hooke’s law, has been widely adapted to explain
rock behavior under stress. It is developed to describe the relation between two physical quantities:
stress and strain. The relation is as follows:
𝜎 𝑖𝑗
= 𝜆 𝛿 𝑖𝑗
𝜖 𝑘𝑘
+ 𝜇 ( 𝜖 𝑖𝑗
+ 𝜖 𝑗𝑖
) , (2-1)
Where 𝜎 and 𝜖 are stress and strain respectively, 𝛿 is kronecker delta, and 𝜆 and 𝜇 are lamé
constants.
Advances in rock mechanics showed the deficiency of elasticity models to explain the
micro and macro scale physics in the rock. Therefore, more sophisticated laws have been suggested
that more accurately describe rock behavior (Papamichos, 1999). Two types of constitutive laws
have been developed to depict the physics of a rock: (1) a general model developed for different
8

rock and stress condition and, (2) a problem-specific constitutive law. In the oil and gas industry,
geo-mechanical experiments of specimens are not planned in the reservoir characterization
procedure. Therefore, general laws should be adapted in modeling -- which in itself may have
many disadvantages. Some are listed below:
 Effects of creep may be overestimated
 Without experimental analysis, there is no guarantee that the model matches the rock behavior
 Limited number of material properties are included in the model
In addition to those mentioned above, there are additional downsides when applying this
method. Nevertheless, using general constitutive law is inevitable.
Every constitutive model should have the following properties (Faruque, 1983):
 Yield condition
 Elastic part of stress-strain relationship
 Flow rule
 Hardening rule
 Strain softening
Yield condition is used to describe the initiation of inelastic behavior in the rock. In other
words, rock shows elasticity before a yield condition; when the stress exceeds this point,
inelasticity emerges. Therefore, every model which incorporates inelastic behavior should define
the plasticity relative to the yield condition.
The elastic part of the stress-strain relationship is obtained through experimental analysis
of the sample. When the sample no longer shows elastic behavior, it undergoes inelastic
deformation, which consists of two separate parts -- elastic deformation, and plastic deformation.
9

Elastic deformation at this range should be defined, and calculated into every model in which the
goal is to describe rock behavior.
Flow rule refers to the plasticity behavior during inelastic deformation. It is defined as the
fluid-like behavior of rock when the stress exceeds the elastic threshold. However, no unique
function can be defined, since the stress rate may change from one case to another. A basic value
can be used in case of infinitesimal increases of stress.
Hardening rule signifies the increase in stress level of the material. As mentioned, when
the stress exceeds the elastic limit, the material shows plastic behavior. However, plasticity does
not necessarily mean a decrease in stress levels. Hardening may occur because of microstructural
rearrangement in the material, and locking of crystals. Three kinds of hardening are defined:
isotropic, kinematic, and anisotropic.
Rock shows the softening behavior once the stress level exceeds a certain value. The stress
increases to peak failure, after which, with continuous deformation, the material degrades and
follows the softening regime (Lade and Kim, 1988). The softening curve can be exponential,
linear, or an experimental defined model.
Materials show different responses to the applied stress. Elastic materials only have the
elastic part of the constitutive relation. For the elastic-plastic materials, depending on the
morphology, structure, and rock type, the material may behave differently; a variety of responses
can be expected. In Figure 2.1, several stress-strain curves are shown, which can be attributed to
different materials (Ranjan et al, 2007).
10

Figure 2.1: Constitutive response of different materials

Based on the abovementioned parameters, different constitutive laws have been created to
describe rock behavior. Many of those laws are case-dependent, and the model is tuned with the
available experimental data (Papamichos, 1999; Schutjens et al, 2004; Sulem and Ouffroukh,
2006; Schutjens et al, 2001).
2.3 Fracture Models in Rocks
Fracture initiation has been the subject of many prior studies. Often, the models were based on
experimental observations of crack initiation in the rock. Those models were shown to be deficient
in real field applications. Scale dependency of fracture initiation is an important parameter that
cannot be accurately incorporated in experiments. However, some attempts were made to estimate
the value from mineback experiments (Warpinski and Teufel, 1987; Jeffrey and Settari, 1995;
Jeffrey et al, 1992; Jeffrey et al, 2009). In general, the two most important concepts in rock
Yeild Point
Linear Elastic
(a) Mild Steel (b) Non-linear elastic (c) Visco elastic
(d) Rigid plastic (e) Elasto plastic (f) Brittle, work-hardening, and
work-softening
Brittle
Work-hardening
work-softening
Velocity
Gradient
σ
σ
σ
σ σ σ
ϵ
ϵ
ϵ
ϵ
ϵ
τ
11

fractures are: scale dependency of the fracture growth, and non-linearity of rock properties. An
acceptable model should be able to address these issues and to offer solutions.
Fractures develop from micro to macro scale (Huang et al, 1993; Atkinson and Meredith,
1987; Zhao, 1998). The initial thinking related the propagation of macro scale fracture growth to
micro scale fracture formation. Small fractures coalesce and form a larger one. Large scale
fractures are those which contribute to fluid flow in the reservoir by acting as the conduits. The
degree to which large fractures contribute to the overall fluid flow varies. Matrix permeability,
fracture orientation, fracture density, and fracture width should all be estimated in order to
understand the benefits of fractures to the fluid conductivity of rock. Furthermore, rock matrix
should be characterized and combined with fracture properties to increase the accuracy of
calculations. Grain size, cementation, rock composition, temperature, and fluid in the pores may
alter rock behavior that has been loaded under different compression and tension forces.
In materials science, multiscale fracture initiation and evolution are investigated for
material characterization and damage studies. It begins with atomic scale and includes micro,
meso, and macro scales. In a macro scale study, fractures may extend for tens of meters in the
formation (Figure 2.2).
Based on the nature of the material and accuracy of the calculations, a scale is selected. In
engineering, problems such as fracture propagation in rocks during hydraulic fracturing, a macro
scale and continuum approach is taken. The science and engineering approach to material studies
in terms of scale is shown in Figure 2.3 (Buehler and Theodor, 2007).
12

Figure 2.2 Presence of the fracture in the large scale
(http://www.deuceofclubs.com/rv/cal230b.htm).

Figure 2.3: Study of a problem in different scales (Buehler and Theodor, 2007).

13

Based on the scale system, different techniques are utilized to investigate the problem. At
the atomic scale, statistical approaches are better suited to model the molecule interactions. Particle
based methods such as the Distinct Element Method is used to model system behavior. In the
engineering approach, continuum models are used to predict system response to the parameter
changes.
In addition to scale dependency of the solution, non-linearity of the system behavior may
cause deviation from conventional solutions and result in a misleading prediction.
Fracture initiation is a function of many parameters -- such as stress regime of the reservoir,
rock heterogeneity, rheology of the fluid, and existence and properties of natural fractures in the
formation (Busetti, 2009). For example, the effects of stress direction to crack initiation is shown
in Figure 2.4 (Nicksiar and Martin, 2012). Therefore, an inclusive model should take into
consideration the major characteristics of the rock and fluid so as to properly calculate the initiation
of the fracture, and its subsequence extension.

Figure 2.4: Sample response to the load of lateral, volumetric, and axial forces (Nicksiar and
Martin, 2012).

14

It is necessary to explain the crack initiation process before we elaborate on the developed
models in the field. As explained in the constitutive laws, the definition of the hardening rule is an
important part of the law for various materials. It is crucial for the rock since the process is slightly
different. Strain-softening characteristics which are observed in rock are totally dissimilar to many
other materials. Therefore, another type of model should be defined specifically for rocks.
A rock breaks down after high pressure fluid is injected inside the sample. Before the crack
initiation, rock experiences a certain deformation which is not necessarily categorized under the
elastic deformation laws. When the stress exceeds a threshold, rock starts to show plastic
deformation. In this state, a portion of rock near the stress zone begins to soften, and then carries
less stress. Meanwhile, the other portions start to unload elastically, which then creates a new stress
field. That phenomenon is called strain localization. The region under the influence of plasticity
has the same dimension as the crack. It can still transmit the stress, but its capacity to load stress
decreases if the strain increases. After the stress exceeds a threshold and the strain exceeds a certain
value, the softening region can no longer transmit the stress, and crack initiation begins. That is a
classic description of the crack initiation in rock (Boone and Ingraffea, 1989). To explain this
process mathematically, further theoretical study is needed.
Several models are suggested that describe fracture initiation and propagation. Each seems
suitable for the specific conditions for which that model was developed. Few of those models are
widely used in either industry or academia. Herein, we focus on methods which are broadly
accepted within the engineering community.
Two types of materials are usually considered in the field of fracture mechanics -- brittle
and ductile. In brittle materials, the fracture occurs suddenly, with little or no trace of plasticity
15

observed in the material. The process zone is small or negligible as compared to ductile materials.
Brittle materials could be modeled using Linear Elastic Fracture Mechanics (LEFM). In ductile
materials the effect of plasticity is not negligible; therefore that influence should be incorporated
into the formulation. Rock can be either brittle or ductile; the selection of the most suitable model
should be done carefully based on field evidence and experimental analysis. There are three main
approaches to deal with fracture in the material. They are:
 Energy based method
 Stress based method
 Displacement based method
2.3.1 Energy Based Method
This technique is based on the energy balance in the material. It has been proven that surface area
has more energy than the body. Therefore, creating a new surface in the form of a crack demands
energy consumption. A relationship can be developed if all the sources of input, output and saved
energy are identified. In other words -- the crack grows if the energy given to the crack surface is
sufficient to prevail over the atomic bonds in process zone. Griffith first proposed the idea of
energy balance and developed the formulation.
Energy absorbed in the crack growth is represented by surface energy; stored strain energy
(internal energy) is released when a crack grows. If and when the stress exceeds the fracture
toughness, the fracture starts to grow. However, the dynamic properties of crack growth are not
included in the model.
16

A parameter called energy release rate is defined to include the rate of potential energy
change by crack extension. Energy release rate at the fracture initiation, Gc, is a threshold value
that shows the starting point of fracture growth. Thus, a fracture starts to grow when G = Gc.
In plastic-elastic materials, a new behavior can be identified. These materials follow the
same path of elastic material during loading, but they deviate from the loading path when they are
unloaded. To model this behavior, J integral is defined to include nonlinear behavior of the
material. It is known that the nonlinear energy release rate J could be defined as a path independent
line integral. Because of the method used to define J, it can be seen as both an energy parameter
and a stress intensity parameter.
2.3.2 Stress Based Method
In addition to the fracture toughness, the stress intensity factor can be defined to characterize stress
distribution at the crack tip. If we assume that fracture initiation is related to the stress and stress
condition, the stress intensity factor represents a combination of those parameters as the criterion
for crack initiation. The crack begins to grow when the stress intensity factor exceeds the threshold
value.
As seen elsewhere (Shah, 1995), energy release rate and the stress intensity factor are
equivalent in a linear fracture mechanics approach. They can be related in an equation which will
be shown later.
2.3.3 Displacement Based Method
Also known as crack tip opening displacement (CTOD), the displacement based approach includes
the blunting of the crack tip. In certain conditions, when the material toughness is significant, the
initial assumption of the sharp crack tip is no longer valid. This phenomenon is a result of plastic
17

deformation of the material at the crack tip. Therefore, the infinite radius in the tip is substituted
with a finite radius that should be estimated. Several studies have investigated the phenomenon
and a model was developed which included it in the calculation. Two methods are outlined below:
2.3.3.1 Irwin model
In this model, Irwin (1957) extended the crack radius to include the plastic zone. He defined an
effective radius by the addition of both crack length and plastic zone. Plane strain and plane stress
formulations have been subsequently obtained.
2.3.3.2 Dugdale model
The strip yield model was developed by Dugdale (1960) to capture the extension of the plastic
zone in the deformation of materials. The existence of a strip of plastic zone is assumed in the
calculation. Based on that model, material is considered as the elastic–perfectly plastic and it is
subjected to the plane stress condition. Two linear elastic solutions can be obtained for this case
and superposition of these two solutions relaxes the assumption of singularity at the crack tip.
Dugdale (1960) and Barenblatt (1959) are co-founders of the cohesive approach, which
takes into account the process zone in front of the crack tip. Dugdale is considered to be a pioneer
of the cohesive zone model.
2.4 Cohesive Zone Model
As just stated, Barenblatt and Dugdale developed the theory of a cohesive zone. They included the
effect of the process zone in front of the crack tip. In their theory, fracture initiates in the
neighboring region of the crack. Two virtual surfaces exist that control the crack initiation and
propagation. The constitutive relation is defined for these two surfaces to characterize the crack
18

behavior and element degradation. The schematic of the crack and its process zone is shown in
Figure 2.5. The figure is modified from Backers (2005).

Figure 2.5: Crack and its process zone in front of the tip. The area shows the degraded element.

According to LEFM, there is a singularity at the fracture tip -- in a perfect elastic model,
the stress at the fracture tip is infinity. That is not an acceptable assumption and does not have a
physical meaning. Therefore, the assumption should be modified to realize the inelastic behavior
in the crack tip region. For metals, the region is represented by a plastic zone and Irwin (1957)
pioneered the research in this area. Hoagland (1973) suggested the existence of a microstructure
zone in the fracture region. The area is usually referred to as a process zone. He performed
experimental analyses on several rock samples, and used Scanning Electron Microscopy (SEM) to
observe the micro-fractures within the zone. Such deviation from LEFM dictates the need for a
new model. The Cohesive Zone Model (CZM) offered a solution to the traditional deficiency of
the LEFM, and paved the way for an advanced model which provided a more realistic
representation of rock behavior. The cohesive zone is defined for the region ahead of fracture tip

Traction Free Zone Inelastic Zone
Elastic Zone
Visible Fracture/Microfracture
Process zone
Effective Fracture Length
σ
x
19

and consists of arbitrary numbers of elements in the discretized model. It means that the number
of damaged elements determine the cohesive zone. Therefore, the method depends upon the
meshing scheme and element size. If the selected mesh size is sufficiently small, the effect of mesh
size could be minimized. A comprehensive study is done by Shet et al (2002) to identify different
CZM’s and their advantages and drawbacks. The result of the study is given in Table 2.1.
In CZM, based on traction-separation law, there are two surfaces in front of the fracture
tip. They are called cohesive surfaces, and the relation between traction and separation
characterizes the behavior of the surfaces. The relation is defined as follows:
𝜎 = 𝜎 𝑚𝑎𝑥
𝑓 (
𝛿 𝛿 𝑐 ) , (2-2)
Where 𝜎 is traction, 𝜎 𝑚𝑎𝑥
is maximum allowable traction, 𝛿 is separation across the surface, and
𝛿 𝑚𝑎𝑥
is characteristic opening. Figure 2.6 shows a sample of traction-separation relation curve.

Figure 2.6: The traction-separation law for a material with linear softening

𝜎 𝑑𝑛

𝛿 𝑑𝑛
1
1
20

Table 2.1: Various Cohesive Zone Models and their advantages and drawbacks
21

Table 2.1 (continued)
22

The traction is applied on two surfaces and increases the separation between them. Once
the stress exceeds maximum allowable traction, the material goes through the softening phase. The
displacement of, or separation between, the two surfaces increases until the separation at the tail
reaches a characteristic opening (Jin et al, 2006). Following that, the surfaces separate and cohesive
traction vanishes. The damaged zone in the material fails at this point and a fracture starts to grow.
The CZM is based on theoretical analysis of the fracture initiation and growth. It means
that no physical surfaces are created in the fracturing process. However, physical change occurs
in the fracture tip and micro-fractures nucleate and coalesce. The experimental analyses have been
performed, and acoustic emissions of creation of such small fractures is well correlated with the
development of the damage zone (Cox et al, 1993).
2.5 Fracture Propagation in Rock
The state of fracture propagation is defined as the condition where the crack grows and leaves the
zone of stress nucleation, and it is no longer associated with the zone (Lawn, 1993). Stress intensity
factors and energy models are usually applied to predict crack initiation and propagation time.
According to these models, when the stress intensity factor, or energy of the fracture, exceeds a
threshold, a crack starts to grow. The fracture growth can be static or dynamic. At this time, we
will not cover the dynamic fracture growth, or fatigue growth in rock. Therefore, the mechanisms
and processes are limited to quasi-static or static crack growth.
Fracture patterns in real materials are usually complex and it is impossible to find the exact
location of a fracture. Nevertheless, a model can be presented to simplify the fracture structure. In
the model, three basic modes of fractures are defined that distinguish fracture types based on the
23

relation between stress direction and fracture displacement. Modes of fracture are shown in Figure
2.7.

Figure 2.7: Crack propagation modes.

Three modes of cracks are as follows:
 Mode I: Perpendicular to the fracture and fracture front
 Mode II: Parallel to the fracture and perpendicular to fracture front
Mode III: Parallel to the fracture and fracture front
Mode I is the one most often reported in fracture mechanics. It is usually measured by the
P-wave recording tools at the surface. Mode II is detected by the shear wave receivers at the
surface. Mode III is not studied in geo-mechanics since tearing does not occur in the hydraulic
fracturing. In hydraulic fracturing, detection of shear wave and primary waves is common and can
be attributed to the presence of Modes I and II fractures in the rock. In fact, a geo-mechanical study
of rock breakdown shows mixed mode fracture growth in many scenarios (Al-Shayea, 2005).
24

Many factors should be incorporated in the calculations to predict crack propagation
trajectory in the rock. Uncertainty may emerge since no accurate method is available to estimate
all the parameters. These parameters could be determined through conducting experimental setup
and utilizing advanced characterization methods and tools, as well as accurate monitoring of the
operation. Stress field magnitude and orientation in the formation, rock heterogeneity, and natural
fractures should be determined to generate input data for the fracture propagation simulator. If so,
a physical model could be developed that would accommodate most of the mechanisms and
processes. Experimental study should be performed to characterize the rock and to ensure that
appropriate values are used to describe the rock. These experimental studies should aim for the
hydraulic fracturing in particular. Most of the studies are focused on tectonic force created
fractures in which the process varies from fluid driven fracture generation (Eichhubl and Atilla,
2003; Weinberg and Regenauer-Lieb, 2010; Lore et al, 2002; Eichhubl, 2004; Zhao, 1998). In the
tectonic induced fracture propagation, the fractures propagate by slow motion fracture elongation
through pores. In this process, fractures are formed by a coalescence of voids. It can also be
developed through high speed tectonic plate movements along the faults. In those cases, dynamic
fracture propagation is observed in the rock. In the dynamic fracture, process, rock properties, and
mechanisms are different from the low speed crack propagation. In addition to the speed of fracture
propagation, we understand that fractures from many scales propagate simultaneously. Micro-
fracture coalescence forms larger fractures which then merge to create a network. It is quite
different than the mechanism in hydraulic fracture propagation.
Many attempts are made to capture fracture propagation in the rock. Assumption of the
type of rock was a base from which to develop numerical models in these studies. Cementation of
the material determines bonding between particles in the formation. If the rock constitutes separate
25

granular particles, and no bonding is present, then distinct particle-based models can be utilized to
predict fracture initiation and propagation in the rock. However, this approach might not be valid
in most of the reservoir rocks because rock diagnosis, high pressure, and high temperature render
an environment that cements materials and creates the bonding. Therefore, a model should be
tailored that integrates the effect of rock properties in the calculations. The damage properties of
rock, as well as any natural fracture distribution, should be determined and then incorporated into
the model to offer a realistic description of the formation.
Finally, the results of micro scale modeling of rock should be integrated with the laboratory
tests at the macro scale. They may then be used for analysis in the continuum models and
simulators, such as Abaqus. The transient state from micro to discrete to macro is shown in Figure
2.8.
Multi-scale approaches are used to investigate the problems in the fracture mechanics.
Continuum models are commonly used in large scale systems, and are evaluated by differential
equations that are normally solved by numerical methods such as Finite Element (FEM) and Finite
Difference (FDM). Simplified forms of the solution can be obtained by applying analytical
techniques to solve the problem. These forms of solution are tradeoff between having higher
accuracy and having a faster solution method. To lessen the inaccuracy arising from utilizing
analytical methods, numerical methods are applied.
26

Figure 2.8: Transient from micro scale to macro scale. Discrete medium is the middle stage and
fits for the uncemented rock materials (Rojek and Oñate, 2007).

Multi-scale approaches are used to investigate the problems in the fracture mechanics.
Continuum models are commonly used in large scale systems, and are evaluated by differential
equations that are normally solved by numerical methods such as Finite Element Method (FEM)
and Finite Difference Method (FDM).
Analytical solutions are developed based on certain assumptions, which limit the
underlying physics incorporated in the model. Some of these methods will be discussed later. Most
of the time, Linear Elastic Fracture Mechanics (LEFM) is used to model rock behavior, although
this may not be the most accurate, as the rock may show plastic behavior. Other misleading
assumptions could be: homogeneity of the reservoir, plane strain model, planar fracture growth,
and homogeneous stress field. These assumptions are easily violated in real case studies;
nevertheless, they are widely used.
In reality, most of the reservoirs are heterogeneous in different scales. Geological features
such as faults and multi layers are present in the macro and mega scales. Also, in the micro scale,
27

grain size and composition, and natural fractures may vary at small intervals. Some features of the
micro scale should be refined before being incorporated into the engineering calculations.
Many experiments are conducted to study the effect of stress on the rock sample in the
micro scale. Result of DiGiovanni et al experiments is shown in Figure 2.9 (DiGiovanni et al,
2007).

Figure 2.9: Damage evolution in rock in micro scale. a) Unloaded rock, b) Rock under load
lower than failure, c) Rock break down under a high compressive stress (DiGiovanni et al, 2007).

In hydraulic fracturing operations, the initial fracture growth can be assumed as planar.
This initial fracture is a transient state between micro and macro scale; it transcends the micro
scale as the stress field after fracture propagation is different from the initial state. Also, it cannot
be captured accurately in continuum model since it still can be shaped by grain properties. Figure
2.10 shows fracture behavior at the early stages. Macro scale approaches can be taken when the
fracture size passes a certain threshold, depending on the type of reservoir rock.
a b c
28

Figure 2.10: Fracture growth modeling in a material using Abaqus.

Macro scale heterogeneity must be considered in the model to offer a realistic picture of
the reservoir behavior. Crack initiation and evolution mechanisms are different in the macro scale.
The effect of geological features is intensified in these types of models. In the reservoirs,
geological features may deviate hydraulic fractures from the predicted pattern and create complex
fracture networks. Mineback experiments have been performed seeking understanding of real
fracture behavior in rocks in the large scale (Kaiser et al, 2013; Warpinski and Teufel, 1987). As
seen in Figure 2.11, fractures propagate in different directions, which makes it impossible to
propose a comprehensive model to predict fracture pattern accurately. However, a few simplifying
assumptions are made to estimate the fracture network propagation. Several models are proposed
for this purpose, some of which are shown in Figure 2.12.
20 mm
29

Figure 2.11: Fracture network pattern in mineback experiments (Warpinski and Teufel, 1987).

Figure 2.12: Fracture network models to estimate fracture propagation in the reservoir (Cipolla et
al, 2008).

30

2.6 Fracture Propagation in Hydraulic Fracturing: a Theoretical Approach
Limited computational capacity of simulators and the lack of advanced characterization methods
have obliged researchers to seek simplistic models to simulate fracture propagation. Several
theoretical models are suggested which are applicable for certain conditions (Khristianovich and
Zheltov, 1955; Perkins and Kern, 1961; Geertsma and De Klerk, 1969; Nordgren, 1972; Abe et al,
1976; Settari and Cleary, 1986; Warpinski and Berry Smith, 1989). Most of the models are 2-D
and only a few are developed for the 3-D fracture modeling.
Three widely used fracture modeling methods and the applicability of each follows:
2.6.1 Perkins-Kern-Nordgren (PKN)
Perkins and Kern (1961) introduced a model to estimate fracture length and width as a function of
time. They assumed that the fracture is vertical and had a fixed height. They incorporated the effect
of stress field as the gradient between pressure inside the fracture and stress field at the boundary
of fracture. Also, the crack shape was assumed to be elliptical. Fluid loss was neglected in their
model. Subsequently, Nordgren modified the PK model to add the effect of leak-off to the
calculations (Figure 2.13).
In the PKN model, as mentioned, the height of fracture is fixed and the fracture propagates
in length. Since the width of horizontal fracture does not vary in horizontal direction, the state of
plane strain prevails in vertical planes perpendicular to the fracture plane. The assumption of plane
strain implies that strain presents only at the direction perpendicular to the fracture propagation.
This supposition is only applicable when the length of fracture is much larger than its height.
31

Figure 2.13: Fracture schematic of PKN model (Economides and Nolte, 2000).

2.6.2 Khristianovitch-Zheltov-Geertsma-De Klerk (KGD)
One of the initial attempts to model hydraulic fracturing was made by Khristianovitch and Zheltov
(1955). They assumed that fracture height was much greater than fracture length. Therefore,
fracture width is independent of the vertical position at any distance from the well. They
formulated a set of equations based on this assumption.
In order to solve the equations, flow rate at the fracture is assumed to be constant. Also,
they alleviated the effect of pressure in the fracture by assuming a constant pressure in the fracture,
except at the region near the crack tip.
This model was later modified by Geertsma and de Klerk (1969) in case of no fluid loss
(Figure 2.14).
32

Figure 2.14: Fracture schematic of KGD model (Economides and Nolte, 2000).

2.6.3 Radial Model
The early efforts to model penny-shaped or radial model fractures were made by Sneddon in a
general study of radial fracture growth (Sneddon, 1946). The method was adopted later by Abe et
al (1976) for the modeling of hydraulic fractures. Fluid pressure is assumed to be constant and
fluid loss is a linear function of any pressure difference across the fracture boundary (Figure 2.15).

Figure 2.15: Fracture schematic of penny-shape model (Economides and Nolte, 2000).

33

Theoretical models are developed based on some of the following simplifying assumptions:
 Fluid properties are constant (however, PKN was modified to incorporate power law
viscosity change)
 Problem is two-dimensional (a few three-dimensional models are developed with limited
applicability),
 Fracture is planar
 LEFM models fracture behavior
 Reservoir is homogeneous and isotropic, and has a constant stress field
 In case of leak-off, it is modeled by simple filtration theory
These assumptions lessen the reliability of the simulation. Therefore, advanced numerical
models should be used to offer a rigorous predictive model.
2.7 Fracture Propagation in Hydraulic Fracturing: A Numerical Approach
The application of numerical schemes in hydraulic fracturing has increased dramatically, due to
the abilities of computers to solve the Partial Differential Equation (PDE). The problems are
commonly formulated in a PDE form, and then simplified in order to find the solutions.
Computational development has enabled researchers to make problems more complex, and to use
numerical methods in the resolution. Several numerical approaches that are utilized to model
hydraulic fracturing will be discussed, including the merits of the most recent one - Extended Finite
Element Method (XFEM). We will detail how it builds upon the FEM to address the issues of
previously developed models.
34

2.7.1 Finite Difference Method
Finite Difference Method (FDM) is a widely used numerical method in engineering calculations.
To solve a PDE, a discrete form of derivatives is employed by using Taylor’s expansion to convert
the PDE to FD form. Implementation of FDM is relatively simple and it is fairly accurate when
compared to the other methods (Figure 2.16). However, its shortcomings are exposed when dealing
with irregular shapes. Although many efforts were made to develop an unstructured grid for FDM,
that grid is not as widely used as other numerical methods (Fernandez, 2004).
Azim and Abdelmoneim (2013) integrated Local Grid Refinement (LGR) with the FEM to
overcome the limitation of FEM in capturing discontinuity of the model in hydraulic fracturing.
However, the work did not address the hydraulic fracture propagation and tried only to model
fractures during production.
In a study conducted by Segura and Carol (2003) FEM and FDM were combined to
simulate hydraulic fracturing. They utilized the FEM to model solid deformation under stress and
took advantage of the flexibility of the FDM to model fluid flow within the fracture.
In addition, some meshless methods are practiced based on the FDM. For example,
Generalized Finite Difference Method (GFDM) was formulated to overcome deficiencies
concerning the conventional FDM (Dow et al, 1990).
35

Figure 2.16: Finite Difference Method implementation in the 2-D space

2.7.2 Discrete Element Method
Discrete Element Method (DEM) is not categorized as the continuum mechanic model. In this
method, the domain is divided into particles or blocks, depending on the nature of the problem.
Equation of motion is solved for the particles or blocks. Since the model is primarily developed
for rocks, it is most frequently applied to model granular rock behavior or soil mechanics (Cundall
and Hart, 1992). However, with just a few modifications, it can also be used to model rocks. Three
main features should be addressed when the DEM model is utilized (Jing, 2003):
1- Identification of contact surfaces between granular materials
2- Calculation of equation of motion under applied force
3- Topological characterization of the medium
36

The basic difference between the DEM and the continuum model is the nature of the
components that constitute the medium. In the continuum model, the discontinuity is introduced
to the system. However, for the DEM it is explicitly defined. DEM is more applicable when there
is a high deformation and rotation in the meso scale. In that case, grain movement under different
forces can be readily studied. However, at the macro scale, which is the subject of the hydraulic
fracturing process, more appropriate methods should be used. There are three issues that prevent
application of the DEM in the field scale hydraulic fracturing simulation:
1- To run a DEM, characterization of granular scale is necessary, which is unavailable for
the field scale hydraulic fracturing simulation in the reservoir.
2- The size of reservoir is considerably larger than the grain size; therefore, an exceptional
amount of computation time is needed to model the operation.
3- Reservoir grains are usually coated by cements, which cannot be handled accurately by
a DEM.
Other methods have been suggested to circumvent some of the limitations in the DEM to
model fracking (Deng et al, 2011; McClure et al, 2010). For example, Deng treated the grains and
cements separately to incorporate the effects of deformation and breakdown. A thorough review
of the DEM is available elsewhere (McClure et al, 2010).
2.7.3 Boundary Element Method
Boundary Element Method (BEM) was developed to solve linear PDEs numerically by integration
of the equation. In this method, the equation is changed in such a way as to solve the PDE based
on boundary of the system, and a solution is calculated in a manner that takes into account the
boundary condition. In the BEM, boundary of the domain is discretized; thus, more resolution can
be assigned for the boundary. Therefore, we can expect higher resolution in its domain when we
37

deal with the BEM. The BEM shows some deficiency when the domain is not homogeneous. In
this case, the method is not applicable most of the time. One of the main advantages of the BEM
over FEM is the computation time. In FEM, the entire domain should be discretized and mesh
assigned for every region of that domain. Hence, a large matrix is generated, which requires
extensive computation time.
The BEM can be useful in various areas of engineering. However, one should seek the
most applicable method to fit the model. Two of the most frequently used methods in LEFM to
capture model deformation under stress are the integral equation, and the displacement
discontinuity method (Sheibani and Olson, 2013).
Usually the following five steps are taken to calculate the solution using the BEM. (Jing,
2003):
1- Discretization of the boundary of the domain
2- Initialization of the problem by solution approximation
3- Calculation of obtained integrals from PDE
4- Solving the equation using assigned boundary condition values
5- Calculation of stress and displacement inside the domain
As discussed earlier, the BEM is an efficient and fast method in the homogeneous domains.
Basically, the problem is reduced by one dimension when BEM is applied. Therefore, computation
speed increases, but the solution accuracy decreases.
BEM is very efficient when the boundary is of interest. However, it shows inadequacy in
the calculation of internal field values (Anderson, 2005).
38

2.7.4 Finite Element Method
Finite Element Method (FEM) is one of the most widely used methods for numerical analysis of
engineering problems. In this method, a continuous domain is discretized and the function value
is estimated at the discretized nodes. In contrast to the BEM and FDM, the function value is
directly calculated in FEM.
As mentioned earlier, different types of meshes are created to divide the domain into
smaller areas, which enables us to increase the resolution and accuracy of the solution. Many types
of mesh generation methods are suggested for the various fields of engineering. They can generally
be divided into two groups, either two or three dimensional. Triangular and quadrilateral are
widely used in the 2D solid mechanics. In 3D, more choices are available, such as tetrahedral,
pyramid, hexahedron, and prism with the triangular base. Selection of one does not impose the
abandonment of another, since a combination of meshes can be used in hybrid mode. Advanced
software programs such as Abaqus utilizes more complicated mesh shapes which suit the
calculations.
Generally, the FEM is implemented in two steps to solve a problem numerically:
1- The domain is discretized to generate subdomains, called elements; their aggregation
makes the system of interest
2- Combining the equations for each node, and solving the global matrix to calculate the
function value at each node
There are many configurations for a discretized domain. Elements should be selected in such a
manner that the combination of all should span the entire domain, but they should not overlap.
Nodes should be selected carefully in order to capture changes in the domain. However, if the
39

number of nodes increase, the computation space, needed for the calculations, increases
exponentially. Therefore, an analysis should be performed to offer an optimized number of nodes,
and shapes of elements.
In the FEM, the function is written as the summation of values at each node. The function
is called shape function and it should be evaluated to obtain the solution for the problem.
Since this report aims to explain the advantages of XFEM over classical FEM, a quick
introduction to both methods is necessary. Here, we will review the FEM and then discuss the
XFEM in the next section.
Assume we have a Boundary Value Problem (BVP) PDE in the following form:
𝐿 ( 𝑢 ) = 𝑞 ( 𝑥 ) , (2-3)
Where u is function of x, and L is a differential operator. It can also be written in the following
form:
𝐿 ( 𝑢 )− 𝑞 ( 𝑥 ) = 0, (2-4)
If we multiply it with a weight function (v(x)) and integrate:
∫ ( 𝐿 ( 𝑢 )− 𝑞 ( 𝑥 ) ) 𝑣𝑑𝑥 1
0
= 0,
(2-5)
If we approximate this function using Lagrangian polynomials, we have:
𝑢 ( 𝑥 ) ≈ ∑ 𝑐 𝑖 ∅
𝑖 ( 𝑥 )
𝑛 𝑖 =1
,
𝑢 ( 𝑥 ) ≈ ∑ 𝑑 𝑖 𝜓 𝑖 ( 𝑥 )
𝑛 𝑖 =1
,
(2-6-1)
(2-6-2)
40

If approximation of u and v are substituted by the original functions, Equation 2.3 is no
longer valid and a residual value will appear. Therefore, we have:
∫ ( 𝐿 ( ∑ 𝑐 𝑖 ∅
𝑖 ( 𝑥 )
𝑛 𝑖 =1
)− 𝑞 ( 𝑥 ) ) ( ∑ 𝑑 𝑖 𝜓 𝑖 ( 𝑥 )
𝑛 𝑖 =1
) 𝑑𝑥 1
0
= 𝑟 ( 𝑥 ) ,
(2-7)
Equation 2.5 is the strong form of the integral function. The strong form should be
substituted by a weak form to be solved by the FEM. Here, we strive to address the formulation of
the problem in the FEM compared to XFEM. Therefore, we will not elaborate on the topic. Further
explanation of the formulation and solution approach is available elsewhere (Zienkiewicz, 1977).
As discussed previously, fracture can be detected in many areas of science and engineering.
Finite Element is used to model fracture in the orthotropic materials (Boone et al, 1987), fault
activation modeling (Dalguer et al, 2003), concrete materials (Hillerborg et al, 1976), tunnel
drilling (Saouma et al, 1981), and rocks (Saouma et al, 1984).
Salehi and Nygaard (2011) studied fracture initiation and growth using FEM in mud
induced fractures during drilling. They used the cohesive zone model, and incorporated stress
change around the fractures after fracture propagation. Settari et al developed a simulator, coupling
the FEM, to model the geo-mechanical properties of the reservoir and FDM to model fluid flow
inside the fractures (Settari, 1980; Settari and Cleary, 1984; Bachman et al, 2003; Ji et al, 2007).
The model consisted of a geo-mechanical section, reservoir fluid simulator, and fracture simulator.
The model was able to capture a 3-D planar fracture growth, as well as pore pressure changes.
Wangen (2013) took into account the ability of the FEM to model fluid flow, and developed a fully
FEM based simulator. He used Biot’s equation to couple fluid pressure and rock deformation. He
assigned a strength threshold to the nodes, and allowed fracture propagation when the stress
exceeded the threshold. Valencia et al (2005) studied hydraulic fracture growth in coalbed methane
41

reservoirs. It was shown that the 2-D modeling of hydraulic fracturing cannot give a realistic
picture of a fracture network. Li et al (2013) employed Mixed-FEM to model the planar fracture
growth in 3-D reservoir model by integration of rock stress change, turbulent flow, and interaction
of both natural and hydraulic fractures. It was demonstrated that accurate determination of
Young’s modulus should be the primary objective for a reliable simulation. Xu et al (2010)
employed FEM to model deformation of rock under stress, induced by fluid inside the fractures.
They used their model to determine the location of fractures through history matching. However,
the final solution was non-unique and inaccurate.
Finite Element is a suitable method to capture rock deformation. Different types of meshes
can be mixed to generate a robust meshing of the domain. It is based on a rigorous mathematical
formulation which captures the essence of fracture initiation and growth. However, a huge
computation expense would be incurred while calculating the fracture growth in the macro scale.
Therefore, a suitable version of FEM was developed to reduce some of the deficiencies associated
with the conventional FEM. Extended Finite Element Method (XFEM) was introduced to alleviate
some of those issues.
2.7.5 Extended Finite Element Method
While the Finite Element Method (FEM) was traditionally employed to capture fracture and crack
initiation and propagation, the re-meshing of an entire domain reduces the computation efficiency
of that method. Therefore, eXtended Finite Element (XFEM) was introduced to deal with that
issue. The XFEM retains the original capabilities of FEM, it can also be employed to model the
effect of cracks. To do so, extra functions have been added to the original shape functions to reflect
any effects of cracks in the element. The presence of cracks may influence the element in two
ways: (1) a fracture is discontinuity within the element and should be incorporated into the shape
42

function, and (2) if the fracture tip is located in the element, singularity arises which needs to be
taken into account. For those reasons, the two functions are added to the shape function. Few
concepts should be clarified before a discussion of XFEM.
2.7.5.1 Partition of Unity (PoU)
The concept of Partition of Unity (PoU) was introduced by Melenk and Babuška to enable the
classic FEM to handle discontinuities (Melenek and Babuška, 1996). In this method, the nodal
shape function is multiplied to the enrichment function (Dolbow et al, 2001) and is added to the
conventional FEM function to incorporate the effect of discontinuity. According to Babuška and
Melenk (1996), PoU can benefit FEM calculation in two aspects:
1- A prior knowledge about the problem can be incorporated into the solution approach
2- In light of its flexibility to adapt to any space, any higher degree shape function can be
constructed for the problem
The fundamental concept is to define a set of functions over space in such a way that the
subdomains form a partition of unity over the space; in other words, the sum becomes one. PoU is
important because it assigns the value of the function to one, and avoids the divergence that may
emerge due to the presence of the crack tip.
Partition of Unity is used for the smooth function to capture non-smooth behavior in the
domain. The function usually incorporates a non-smooth, and perhaps non-polynomial to model
discontinuity of the domain.
In hydraulic fracture problems, fractures usually span a wide area in the formation.
However, most of the elements in the reservoir remain without a fracture. Therefore, PoU is needed
to incorporate the fractures in the shape function only for small portion of the elements. Inclusion
43

of the fracture effect in the element through applying PoU saves some costly computation
expenses, since convergence to the solution does not require the domain re-meshing at each step.
2.7.5.2 Level Set Method
Osher and Sethian (1988) introduced the Level Set Method (LVS) to track front propagation in
fluid. However, it evolved to also capture front movement in solid materials. Stolarska et al (2001)
combined LSM with XFEM to increase the efficiency of technique. In FEM, re-meshing is
performed to keep the discontinuities out of the elements all the time. Therefore, no implicit
inclusion of cracks in the calculation is required. In contrast, XFEM deals with cracks inside the
element. LSM enables XFEM to capture crack attributes inside the element. It represents the crack
location, including crack tip location. Based on the location of the crack tip, the domain is divided
into 4 areas. Two signed functions are used to determine the sign of each area. Distance function
is also incorporated into the formulation to offer the qualitative value of closeness to the crack tip.
Normal level set function, and tangential level set function, serve as the basis on which to divide
the domain. An appropriate function may be assigned to describe the crack tip once it has been
located.
Now, formulation of XFEM can be discussed since fundamental knowledge needed to
understand the method is given.
Equation 2.3-7 described the general formulation of FEM to solve a PDE. Herein, no
general formulation is needed because we focus on the crack tip propagation under LEFM
assumption. Therefore, notion changes to suit the current problem.
44

𝑢 ( 𝒙 ) ≈ ∑ ∅
𝑖 𝑢 𝑖 ( 𝒙 )
𝑁 𝑖 =1
,
(2-8)
Equation 2.8 still describes shape function for FEM. Shape function is modified in XFEM
to consider effect of discontinuities in the element. A term is added to the estimated value of the
node for this purpose. It means:
𝑢 𝑖 ( 𝒙 ) 𝑢 𝑖 ( 𝒙 )+ ∑ Ψ
𝑖𝑗
𝑁 𝐸 𝑗 =1
𝛾 𝑗 ( 𝒙 )
After the transformation, we have:
𝑢 ( 𝒙 ) = ∑ ∅
𝑖 ( 𝑢 𝑖 ( 𝒙 )+ ∑ Ψ
𝑖𝑗
𝑁 𝐸 𝑗 =1
𝛾 𝑗 ( 𝒙 ) )
𝑁 𝑖 =1
,
(2-9)
If the terms be separated:
𝑢 ( 𝒙 ) = ∑ ∅
𝑖 𝑢 𝑖 ( 𝒙 )
𝑁 𝑖 =1
+ ∑ ∅
𝑖 ∑ Ψ
𝑖𝑗
𝑁 𝐸 𝑗 =1
𝛾 𝑗 ( 𝒙 )
𝑁 𝑖 =1
,
(2-10)

In Equation 2.10, two terms are added to estimate the function over the domain. The first term is
the classical FEM shape function which is used to estimate the smooth and continuous function.
The second term or term B is the enrichment term which is used to incorporate the fracture tip
behavior into the estimation.
Equation 2.10 improves upon the FEM for better handling of fracture. Moreover,
discontinuity of the element should be captured using an extra term. Discontinuity is represented
A
B
45

in form of jump across the domain. Every crack inside an element is modeled through adding an
extra term that has this jump. Heaviside function is used to describe the discontinuity in the
domain. It is also called unit step function and is used to represent the piecewise constant function
in the XFEM. Jump over the domain is captured by Heaviside function considering the interface
discontinuity. Heaviside function for the XFEM is defined as the following:
𝐻 ( 𝑥 ) = {
1 𝑥 > 𝑥 0
−1 𝑥 < 𝑥 0
,
(2-11)
Depending on the location of specific point in the domain, a value is assigned to the Heaviside
function.
When the Heaviside function is taken into account, shape function will be as the following:
𝑢 ( 𝒙 ) = ∑ ∅
𝑖 𝑢 𝑖 ( 𝒙 )
𝑁 𝑖 =1
+ ∑ ∅
𝑖 ∑ Ψ
𝑖𝑗
𝑁 𝐸 𝑗 =1
𝛾 𝑗 ( 𝒙 )
𝑁 𝑖 =1
+ ∑ Ω
𝑘 𝑢 𝑘 ( 𝒙 ) 𝐻 𝑁 𝑘 𝑘 =1
( 𝒙 ) ,
(2-12)
Equation 2.12 is the final form of XFEM shape function that is commonly utilized to represent the
crack behavior in the domain. Further explanation of the topic is available elsewhere (Pommier,
2013).
2.8 Concluding Remarks
This chapter described the behavior of rock under the stress, fracture initiation and growth in the
material, and the different approaches to model them. Depending on the type of rock, the response
of the material under stress may differ. To capture the complex behavior of materials, many rock
models have been developed. That experimental setup is needed, in order to offer an acceptable
characterization of a rock in a lab scale. Even with such a model, the generation of a valid field
model is not guaranteed. They should be translated into theoretical models and used for the
estimation of fracture initiation and propagation. The traditional LEFM could not be used in a
46

complex environment considering the simplistic approach. Therefore, the more advanced CZM
was developed to predict fracture initiation and propagation. The model needs to be exported in a
reliable mode to enable one to predict a fracture network. The current fracturing models show
deficiency in matching the field data, and could not be used to calculate the hydraulic fracturing
parameters in a complex geological environment accurately. Given the advances in numerical
modeling, incorporating a realistic rock model, and the introduction of XFEM, fracture initiation
and growth could be assimilated into the calculations, which would address several issues. The
method should then be integrated with a fluid flow simulator to create a comprehensive model for
hydraulic fracturing simulation.
47

Chapter 3 : DESCRIPTION OF HYDRAULIC FRACTURING
SIMULATION ALGORITHM
This chapter will present the algorithm that is proposed to model hydraulic fracturing in the
reservoir by coupling fluid flow and rock deformation in an iterative approach. Two different tools
were developed to model the processes separately. In the fluid flow modeling, a simulation code
was established in Matlab, and based on FDM; the formulation and numerical implementation will
be outlined forthwith.
Abaqus was utilized to predict rock deformation, and fracture initiation and propagation.
The fundamental equations of rock deformation will also be offered, as will the numerical
approach to problem solving.
3.1 Formulation of Fluid Flow Equation
A successful hydraulic fracturing operation depends upon creating conductive fluid paths that
allow the reservoir fluid to be produced. Hydraulic fluid contributes in two general ways: (1) it
acts as a conduit to transfer the energy generated at the surface to the rock, resulting in rock
breakdown and creation of a fracture network and, (2) it carries the proppants and deposits them
inside the created fractures in order to keep them open after the operation is complete. Therefore,
particular properties need to be considered to ensure the selection of the most appropriate hydraulic
fluid. In this study, the focus is on the first aspect of hydraulic fracturing fluid. The equations
describe only the fluid behavior associated with pad fluid injection.
Fluid flow in the reservoir is modeled using a comprehensive approach that incorporates
the fluid flow in the fracture, matrix, and fluid flow displacement at the boundary. Cross flow at
the boundary of fracture and matrix is taken into account in the model. However, it is assumed that
48

fluid does not violate the laminar flow in the fracture boundary; the fluid flow in the fracture can
be modeled through negligible boundary flow assumption. Fracture and matrix consist of two
different media in terms of fluid behavior and medium type. Fluid flow in porous media is modeled
through integration of conservation of mass and Darcy’s law, which is called diffusivity equation.
The diffusivity equation is given as the following:

∇
2
𝑝 =
𝜇𝑐𝜑
𝐾 ∂p
∂t
, (3-1)
Where p is fluid pressure, μ is fluid viscosity, c is compressibility factor, φ is porosity, K is
permeability of matrix, and t is time.

Figure 3.1: Fluid moves inside the crack and simultaneously transfers between crack and matrix

Diffusivity equation is utilized to model multiphase multi-dimensional fluid flow in the
porous media. However, an extra term is needed to describe the rock deformation under pressure
change. Therefore, the equation is modified to incorporate the effect of fracture opening in fluid
flow, as follows:

∇( 𝑤 3
∇𝑝 ) = 𝜇 (𝑐𝑤 ∂p
∂t
+
𝜕𝑤
𝜕𝑡
),
(3-2)
49

Where w is width of fracture. The parallel plate formulation for flow in a fracture was widely used
(Snow, 1965; Romm, 1966; Louis, 1969; Sharp and Maini, 1972; Wilson and Witherspoon, 1974;
Gale et al., 1974) and leads to a fracture permeability defined as

𝑘 =
𝑤 2
12
,
(3-3)
It is necessary to point out that the fluid flow from fracture to the reservoir does change over time.
According to the Carter’s model (Howard et al, 1957), it is a function of time and leak-off
coefficient. The equation is given as:
𝑢 𝐿 =
𝐶 𝐿 √𝑡 ,
(3-4)
Where 𝐶 𝐿 is leak-off coefficient, and t is time. The integrated form of the Carter equation is a
function of the fracture width (Economides et al, 1999). The model is developed based on a few
assumptions that might be violated under certain conditions. For example, fracture width might
not be constant in many cases. Also, other studies stated that the assumption of orthogonal flow to
the fracture-matrix boundary may be violated considering development of a leak-off zone inside
that matrix (Khodaverdian and McElfresh, 2000). Moreover, studies show that the leak-off
depends on various factors (Mayerhofer et al, 1995). At the beginning stages of fluid injection, the
filter cake is not yet formed on the fracture wall. Therefore, the effect of filter cake on the fluid
leak-off is negligible. The fluid flow into the reservoir is influenced primarily by filter cake once
an impermeable layer is formed on the fracture wall. Earlier, the effect of a pressure difference
between hydraulic fracture and the reservoir was a prominent factor in determination of the leak-
off. In this study, we assume there is not enough time for the filter cake to form in the fracture
wall. Therefore, the pressure difference between fracture and reservoir is a primary factor for
50

estimation of the leak-off. Based on this assumption, we will be able to create an integrated
approach to calculate pressure in the fracture and reservoir in a single function.
Traditionally, Darcy’s law was developed to describe the fluid flow in porous media. It is
essentially derived from Navier-Stokes by averaging procedure (Neuman, 1997). Here, a similar
technique is practiced for the fluid flow in the fracture to offer a uniform, consistent model and to
bridge the gap between fluid flow in the matrix and fluid flow inside the fracture. A fracture is
divided into cells, and conductivity of the cells is converted to the permeability using the following
relation:

𝐾 𝑓 ∝ ℎ
𝑓 𝑤 𝑓 ,
(3-5)
Where subscript f is fracture, h is height of fracture, and w is width of fracture.
Coupling of fluid flow in fracture and matrix is performed easily, since no complicated
procedure is incorporated for fluid flow modeling at the fracture and matrix boundary.
3.2 Numerical Analysis of Fluid Flow Equations
Finite Difference Method (FDM) is used to discretize the equations. Fluid flows in the fracture
and matrix; it also transfers from fracture to matrix and vice versa. Equations in each region need
to be discretized separately, considering the type of equations in fractures and in matrix.
Permeability of a fracture is defined in such a way that provides us with a tool to integrate the fluid
flow in porous and nonporous media. Therefore, a fracture is treated as a porous medium, and its
permeability is defined according to (3-3).
51

Figure 3.2: Fracture and matrix discretization for fluid flow modeling

Fracture width is narrow, hence, vertical flow can be neglected in the modeling. It means
that there is no pressure difference between two fracture walls.
As discussed previously, fluid flows inside the fracture and transfers from fracture to the
matrix and also flows inside the matrix. Three distinct regions can be recognized for fluid flow--
inside the fracture, fracture-matrix boundary, and inside matrix. Discretization of the matrix is
straightforward assuming the size of fracture block is small relative to the matrix block. Cubical
elements are used for meshing the space. However, it is more difficult for discretization of the
fracture and fracture-matrix boundary. Fracture-matrix region is discretized assuming that there is
no pressure change between fracture walls or Kf>>Km. Given that assumption, discretization of
the region is performed.
The fracture width variation is incorporated into the calculations. Therefore, the
neighboring elements do not have similar width and permeability. Here, fluid flow inside the
fracture is modeled using a method which incorporates the consequence of size variation in
i,j,k
i-1,j,k
i+1,j,k
i,j,k-1
i,j,k+1
i,j,k i-1,j,k i+1,j,k
i,j,k+1
i,j,k-1
52

modeling (Ertkin et al., 2001). It is based on the harmonic averaging of “Area×Permeability” for
adjacent elements and can be written as follows:

( 𝐾𝐴 )
𝑎𝑏
= 2
𝐾 𝑎 𝐴 𝑎 𝐾 𝑏 𝐴 𝑏 𝐾 𝑎 𝐴 𝑎 + 𝐾 𝑏 𝐴 𝑏 ,
(3-6)
Ka and Kb are the permeability of element a and b, and Aa and Ab are area of fracture sides for
element a and b at the interface of two elements. Area of the fracture is fracture width multiplied
by fracture length:

𝐴 = 𝑤𝑙 ,
(3-7)
Length of fracture is constant and fracture permeability is proportional to fracture width. After
substituting the variables and simplifying the equation, we have:
( 𝐾𝐴 )
𝑎𝑏
= 2𝐶 𝑊 𝑎 3
𝑊 𝑏 3
𝑊 𝑎 3
+ 𝑊 𝑏 3
, (3-8)
Where C is conversion factor from SI unit to field unit.
3.3 Stress Modeling
Many geomechanical constitutive equations are offered to model rock behavior before and after
damage (Homand-Etienne et al, 1998, Jing, 2003, Lu et al, 2014). These models are developed for
certain materials and conditions. Reservoir rock may show different behaviors under stress.
However, it does typically obey linear elasticity law, except at the crack tip region. Here, linear
elastic stress-strain constitutive model was employed to describe rock behavior under stress, but
before the yield point. The elastic stress-strain constitutive equation is as follows:
53

𝑒 𝑖𝑗
=
1 + 𝑣 𝐸 𝜎 𝑖𝑗
−
𝑣 𝐸 𝜎 𝑘𝑘
𝛿 𝑖𝑗
,
(3-9)
Where 𝑣 is Poisson’s ratio, E is Young’s modulus, e is strain, 𝜎 is stress, and 𝛿 is Kronecker delta,
and is 1 when i=j, otherwise it is 0. In our approach, it is assumed that crack propagates in the
equilibrium. Therefore, stress in the medium is always in the equilibrium state. Taking into account
this assumption, the partial differential equation is given as:
𝜎 𝑖𝑗 ,𝑗 + 𝑓 𝑖 = 0, (3-10)
Where f is body force. Stress equation is solved numerically to determine stress field at the crack
tip.
3.4 Fracture Modelling
Accurate prediction of the hydraulic fracturing process depends predominantly on utilizing reliable
fracture models that properly describe the fracture behavior under stress. Two major models should
be employed for this purpose. First, a damage model needs to be developed to characterize the
crack tip region, and explain the rock properties. Second, a fracture model that is capable of
explaining rock behavior, fracture initiation and propagation. In this study, we relied on both
models that were developed for rock strain-stress curve calculations. It is worth mentioning that
the models can be easily modified based on new laboratory and field property reports.
3.4.1 Rock Damage Model
Rock under stress undergoes a complex process of deformation and break. Accurate
characterization of the transition from the initial elastic behavior to crack propagation demands
comprehensive laboratory and field investigation. Many damage models are proposed to describe
rock behavior (Zhu and Tang, 2006; Golshani Aliakbar et al, 2007; Simo and Ju, 1987). They
54

mainly focused on the process zone in front of the crack tip. The reason for high number of studies
is that rock shows a wide range of stress-strain relationships while the crack is under stress.
Laboratory experiments and field analyses were conducted to determine the relationship. This
relationship is generally represented as a curve and it ordinarily starts with the linear elastic model,
and is then follows by a fracture initiation and propagation model.
The relation between damage initiation and evolution can be attributed to various
underlying mechanisms of crack nucleation and coalescence. In heterogeneous materials such as
reservoir rock, a complex pattern of fracture networks may form that depends on the spatial-
temporal loads on the rock. To model damage initiation and propagation, an accurate rock damage
model should be provided which is able to predict the fracture growth pattern. A comprehensive
study was performed by Busetti (2012) to characterize Berea Sandstone. We have employed the
model in our calculation to predict hydro-fracture growth. A detailed description of the model can
be found at Busetti (2012).
3.4.2 Cohesive Segment Method and Phantom Nodes
As previously stated, various models have been suggested to simulate crack growth onset, and
predict the fracture propagation pattern. Traditionally, LEFM was used to analyze stress regime
and fracture initiation and growth. However, it is understood that except for brittle materials with
low plasticity, LEFM is unsuitable for fracture initiation and propagation modeling in rocks.
Therefore, other realistic models were developed. One of the most recognized is the cohesive zone
model, in which the presence of a process zone is considered in the crack tip region. It is altered
by stress change in the fracture tip, and should be modeled to capture the underlying mechanisms
of the fracture process. In the cohesive zone, a nonlinear relation between surface displacement
and stress does exist and must be realized. The relation is able to propose the following function:
55

𝜎 = 𝜓 ( 𝛿 ) ,
(3-11)
Where 𝜓 is cohesive interaction function. Generally, a linear relation exists between stress and the
displacement of surfaces at low stress, and before the yield point. Once the displacement between
surfaces exceeds a threshold, softening occurs, which leads to materials failure at the
characteristics length. Modeling of the crack path in the cohesive zone approach needs the presence
of a predefined fracture path. In homogeneous materials, it is practical to do so, since the path of
the fracture can be estimated with above average accuracy. However, in heterogeneous materials
such as reservoir rock, a predefined fracture growth path may lead to erroneous results. To address
this issue, cohesive segment is employed, which is capable of modeling a fracture propagation
pattern without prior knowledge of the fracture path (Remmers et al, 2003). In that approach, a set
of degrees of freedom is added to the neighboring elements to calculate the crack growth inside
the element. If the step size is sufficiently small, crack propagation is calculated by constitutive
relation, and it is extended to the element boundary. The crack extends straight, but it can deviate
from the crack trajectory of the previous step. The phantom node method is introduced to enhance
the modeling of discontinuity in the fracture (Song et al, 2006). In this technique, an image domain
is introduced to the system to increase the degrees of freedom. Phantom nodes are defined in the
image domain. They match the real nodes if they are not part of the cracked element. When the
crack passes through the element, phantom and real nodes separate, and constitutive law is defined
for each node set. In this case, XFEM field is calculated by the summation of two displacement
fields which each contain both phantom and real nodes. Detailed discussions can be found
elsewhere (Song et al, 2006).
In our approach, the crack tip might be located inside the element. Therefore, the
asymptotic singularity of a fracture tip exists in form of enrichment function. To alleviate the issue,
56

it is assumed that crack jumps from boundary of element; hence, asymptotic function no longer
needed in the nodal shape function.
3.5 Algorithm Development and Numerical Simulation Scheme
A rigorous approach was developed by the incorporation of various rock, fluid, and rock-fluid
interaction parameters. Wide ranging factors are taken into account to provide a realistic picture
of hydraulic fracturing growth in the reservoir. Fractures are initiated and subsequently propagate
when stress at the crack exceeds the maximum principal stress criterion of the rock. However, the
stress criteria can be modified, based on experimental and field analysis of the reservoir.
Stress field in the reservoir is related to different factors. Prior to any fluid injection,
tectonic forces and overburden pressure determine the stress field. Once the injection begins, the
fluid exerts force against the fracture wall, and disturbs the stress field. The degree of stress field
distortion mainly depends on the fluid pressure and geo-mechanical properties of the reservoir
rock as well as initial state of stress field. Therefore, we see a complex interaction between fluid
flow in the fracture and rock matrix, geo-mechanical parameters of the rock, as well as the in situ
forces in the reservoir. Accurate prediction of hydraulic fracture propagation in the reservoir
demands the incorporation of all these phenomena into a single model. Towards that goal, a
workflow has been constructed that integrates the major mechanisms in hydraulic fracturing. It is
shown in Figure 3.3.
57

Figure 3.3: Workflow for the modeling of hydraulic fracturing

As seen from the Figure 3.3, the workflow offers an approach that takes into account the
complicated fluid flow, rock models and their interaction. Geo-mechanical modeling is performed
utilizing an extended finite element based numerical solver to simulate fracture initiation and
growth. Subsequently, fluid pressure is adjusted once the fractures are propagated in the reservoir.
The algorithm allows a 3D modeling of fracture growth in the reservoir. The model has not
been simplified by the plane theory of elasticity to avoid the limitations that are commonly
imposed in this type of modeling. Three dimensional models must be utilized to offer a
comprehensive simulator, as many features are spatially distributed in the reservoir and cannot be
simply integrated in a two dimensional model. Pre-existing fractures and multilayer reservoirs are
two major sources of spatial heterogeneity in the reservoir, and they can be incorporated into the
geo-mechanical model.
Many techniques exist for multi-physics, multiscale simulation of fracking. Some of those
methods are stable and solve the fluid flow and rock deformation equations in a fully coupled
58

scheme. Others are less stable and loosely couple fluid and geomechanical equations (Kim et al,
2009). The solution can be achieved through utilizing both approaches, but the degree of accuracy
and computation time would be different. Generally, an iterative approach can balance the
computation time and accuracy; hence, it is widely used in numerical simulations. In this section,
we discuss a procedure that may solve the hydraulic fracturing problem. An algorithm is developed
and is tested for the simulation of fracture initiation and growth. There are three major components
of the simulation -- fluid flow, rock deformation modeling, and a connecting rock and fluid model.
Fluid flow domain is meshed by using cubical elements. The equations are discretized by
means of a fully implicit time discretization, also known as backward Euler. However, the term
for fracture width change is treated explicitly. It is updated in each step by exercising the
calculations of a rock deformation modeling tool. LU decomposition was employed to solve the
system of linear equations.
Abaqus solves the rock deformation equations. A static approach is taken for this purpose.
One-step analysis is applied to model solid deformation when the stress is less than failure stress.
However, the program adjusts increment sizes automatically whenever the stress is close to the
failure stress and fracture initiation. Newton’s method was used to solve the system of nonlinear
equations. It is an iterative based method, and searches at the vicinity of the solution until it
converges. Also, Abaqus/Standard serves as a tool to resolve the geo-mechanical equations. It was
developed based on an implicit approach.
Flow equation, and solid deformation equations are solved on two separate domains. An
iterative approach is used to couple fluid flow and rock modeling tools. In each step, the rock
model from the previous step is used to model fluid flow. The calculated pressure field is exported
59

to the Abaqus to calculate deformation of solid model. The solid model is exported back to the
fluid flow model to calculate pressure field. The iteration continues until the error has become
minimized.
An iterative approach has been widely used to calculate a system of nonlinear equations. It
is common to use Jacobian matrix to increase the convergence speed. However, it is not applicable
to our problem. Calculation of Jacobian matrix is cumbersome and inefficient. Therefore, the
Secant method was used to increase the convergence speed.
Instability may occur in the calculations, especially during the early stages of the fracture
initiation. The fracture is small, so an incremental fluid change causes pressure to rise in the
fracture. Therefore, in each step, the fracture changes more than in the previous step, and it grows
exponentially. To make that correction, and ensure a stable solution, the fracture opening in each
iteration is adjusted to avoid an excessive change of fracture volume relative to the amount of
injected fluid. It should be noted that the control of fracture width change can only be applied to
those elements where fluid is injected. Fracture widths for the rest of elements are subsequently
adjusted.
3.6 Concluding Remarks
I have presented the fundamental equations of fluid flow, rock deformation, fracture initiation,
and propagation; their different aspects were explained in detail. Hydraulic fracturing is the
integration of various mechanisms. The interaction of which were discussed and a workflow was
shown which connects them, and all in a single platform. The model can handle various boundary
conditions. The benefits of this approach over previously suggested methods can be summarized
as:
60

- Real-time crack growth can be tracked in a realistic model
- Fluid flow model can model the leak-off term accurately
- XFEM enables us to solve the fracturing problem with a negligible impact on the accuracy
- A physical based approach was utilized to increase the speed of calculations
Integration of these benefits yields an inclusive method which addresses the conventional issues
of hydraulic fracturing simulation, and provides a platform for complicated analyses of the
problem using field data.
This method needs to be tested against the field data to assess its performance, and to more
fully understand the capabilities and limitations of the model. A theoretical analysis will be given
in the next chapter.

61

Chapter 4 : TWO DIMENSIONAL MODELING OF FRACTURING IN
LABORATORY SCALE
The goal of this chapter is to provide insight into fracture behavior in various structural and
boundary conditions. In the previous chapter, the fundamental equations of fluid flow and rock
deformation were discussed, and an algorithm was presented to model hydraulic fracturing. The
algorithm consists of two parts; fluid flow modeling and rock/fracturing simulation. In this chapter,
we disregard the modeling of fluid flow behavior inside the fracture, as it is filled with the fluid
all the time, and fluid pressure is constant throughout the fracture. The rock is downscaled in both
size and dimension to enable us to investigate various mechanisms of fracturing without
considering complicated relations of the fluid and rock parameters. We describe a set of numerical
analyses that have been performed to study crack initiation and propagation in a small scale model.
The analysis starts with an example of single fracture propagation. Sensitivity analysis is
performed on important factors of fracture initiation and propagation. Thereafter a model with
multiple fractures is designed to study stress shadowing, and hydraulic and natural fractures
propagation. For the former, fractures are placed in parallel and for the latter, fracture are
perpendicular. The model is capable of simulating the opening of natural fractures after crossing
of hydraulic fracture. The effect of model boundary and distance between two growing hydraulic
fractures is also studied.
4.1 Introduction
Hydraulic fracturing Simulator could be a useful tool if it takes into the consideration the
underlying phenomena. One of the major phenomena is the hydraulic fracture behavior in relation
to natural fractures. Hydraulic fractures interact with natural fractures and it could result in opening
62

of natural fracture and change in hydraulic fractures path. Many studies have been conducted to
investigate what effect natural fractures have on the hydraulic fracture propagation (Olson, 2008,
Renshaw and Pollard, 1995, Gu and Weng, 2010, Zhou et al, 2008, Raymond et al, 2015). Various
parameters influence the stress field around the fracture tip, which may alter the behavior near
natural fractures. An induced fracture may cross a natural fracture, or it may terminate in the
natural fracture. In the first scenario, a natural fracture can be reactivated and propagate, or it can
have no reaction to the hydraulic fracture, and remains closed. Also, shearing of a natural fracture
is plausible. A reliable model should be able to model hydraulic and natural fracture interaction
and predict their behavior when they are in vicinity of each other. We have used Abaqus and have
integrated it into our algorithm to be able to model interaction of hydraulic and natural fractures.
In addition, the model is capable of simulating the stress shadowing that is an important factor
during multistage hydraulic fracturing.
To be able to accurately model hydraulic fracturing, we need to identify sensitivity of the
model to the numerical parameters. Model size, boundary effect, model resolution, and solution
approach are major numerical factors that could affect the simulation results. Sensitivity analysis
on each factor in field scale is an assiduous task and demands high computational time. Therefore,
we have created small scale models with single fracture to investigate the effect of each factor on
the hydraulic fracturing modeling. Afterwards, the optimum model is selected for studying stress
shadowing and hydraulic and natural fractures interaction.
4.2 Single Fracture Propagation
In this section, we will study the effect of various parameters on a single fracture propagation in a
2-dimensional plate. Static and implicit dynamic analysis were performed for the model. We
considered both plane strain and plane stress conditions in our analysis. Effect of model resolution
63

as well as shape were studied in this section. Effect of the reservoir stress field on fracture
propagation is investigated subsequently.
Three models are constructed to study various parameters. Figure 4.1,
Figure 4.2, and Figure 4.3 show the schematics of the models with an initial fracture.

Figure 4.1: Plate is subjected to two tensile loads acting from the top and bottom of the plate.
Initial fracture is placed at the center left part of the plate.

No extra boundary condition is applied and all parts of the plates are allowed to move. The
characteristics of plates and fractures are given in Table 4.1.
12.5 cm
100 cm
50 cm
64

Figure 4.2: Plate is subjected to the two tensile loads acting from the top and bottom of the plate.
Initial fracture is placed at the center left part of the plate.

Figure 4.3: Plate is subjected to the two compressive loads acting from the top and bottom of the plate.
Initial fracture is placed at the left part of the plate. Pressure inside the fracture increases and applies load
on the fracture walls.

25 cm
100 cm
100 cm
25 cm
100 cm
100 cm
65

As seen from the Table 4.1, the geo-mechanical properties of the material is very close to
the sandstone. The plate shape is rectangular for Figure 4.1, with a ratio of 1:2 between its width
and height and is square with a ratio of 1:1 between its width and height for
Figure 4.2 and Figure 4.3. In the base model, the fracture is at the center, starting from the
left side of the plate. The initial size of the fracture is one-fourth of the total width of the plate.
Constant load is applied from both the top and bottom of the plate. Fracture propagation is based
on traction separation law. We assumed that the plates stay undamaged before stress exceeds the
maximum principal stress. Damage evolution is based on displacement condition.
Table 4.1: Fracture and plate characteristics, and boundary condition
Parameter Value
Plate size for Figure 4.1 (m×m) 1×0.5
Plate size for Figure 4.2 (m×m) 1×1
Initial fracture size for Figure 4.1 (m) 0.125
Initial fracture size for Figure 4.2 (m) 0.25
Young’s modulus (GPa) 20
Poisson’s ratio (fraction) 0.3
Damage initiation criterion Maximum principal stress
Maximum principal stress (MPa) 5×10
6

Damage evolution type Maximum displacement
Maximum displacement 0.01
Fracture criterion type Virtual Crack Closure Technique (VCCT)
Mixed mode behavior Benzeggagh and Kenane (BK)
Energy release rate in mode I, II, and III (J/m
2
) 30, 100, 1000
η 1

66

In this work, we tested the dynamic analysis of fracture propagation for modeling as well.
But that does not mean that the fracture propagation can be categorized under dynamic fracture
mechanics. Hydraulic fracture growth commonly follows the static fracture mechanics, and
dynamic refers to the time dependency of the process. The quasi static approach was selected to
describe the process. Abaqus requires definition of contact surfaces in the dynamic analysis of a
fracture growth. In such cases, we assume the validity of Virtual Crack Closure Technique (VCCT)
to model the behavior of two fracture surfaces while they are in contact with each other. It describes
how the two surfaces react to the mixed mode stresses. VCCT was developed to determine energy
release rate during delamination (crack growth) which depends on the mixed mode ratio; in such
cases, an explicit approach is needed to calculate fracture mode separation (Krueger, 2004).
Mixed mode behavior in VCCT was determined by an equation that was introduced by BK
(Benzeggagh and Kenane, 1996). The relationship is as follows:
𝐺 𝑇𝐶
= 𝐺 𝐼𝐶
+ ( 𝐺 𝐼𝐼𝐶
− 𝐺 𝐼𝐶
)(
𝐺 𝐼𝐼
𝐺 𝑇 )
𝑚 , for mode I and II, (4-1)
Where

𝐺 𝑇 = 𝐺 𝐼 + 𝐺 𝐼𝐼
,
(4-2)
And
𝐺 𝑇𝐶
= 𝐺 𝐼𝐶
+ ( 𝐺 𝐼𝐼𝐶
− 𝐺 𝐼𝐶
)(
𝐺 𝐼𝐼
+𝐺 𝐼𝐼𝐼
𝐺 𝑇 )
𝑚 , for mode I, II, and III, (4-3)
Where
67

𝐺 𝑇 = 𝐺 𝐼 + 𝐺 𝐼𝐼
+ 𝐺 𝐼𝐼𝐼
, (4-4)
And GTC is total critical strain energy release rate; and GIC, GIIC, GIIIC are Modes I, II and
III critical strain energy release rates, respectively.
It is assumed that the critical strain release rate in Mode I is lower than the other modes.
Therefore, Mode I of fracture propagation is dominant; and the energy release rate for the mode I
was calculated using the following equation:

𝐺 𝐼 =
𝐾 𝐼 2
𝐸 ′
,
(4-5)
Where
𝐸 ′
=
𝐸 1−𝑣 2
, for plane strain, (4-6)
And
𝐸 ′
= 𝐸 , for plane stress, (4-7)
Assuming stress intensity factor of 810 KPa/m
0.5
, critical strain energy release for Mode I
is 30 j/m
2
for plane strain condition.
Stress in the model is evaluated using Von Mises stress, which is a criterion that represents
the net energy stored in material due to distortion of the element, but in the form of stress.
Therefore, it is a proper criterion to report, and to compare various analyses. The definition is as
follows:
68

𝜎 𝑉𝑀
= √
1
2
(( 𝜎 𝑥𝑥
− 𝜎 𝑦𝑦
)
2
+ ( 𝜎 𝑥𝑥
− 𝜎 𝑧𝑧
)
2
+ ( 𝜎 𝑧𝑧
− 𝜎 𝑦𝑦
)
2
) + 3( 𝜏 𝑥𝑦
2
+ 𝜏 𝑦𝑧
2
+ 𝜏 𝑥𝑧
2
) ,
(4-8)
In 2-D model, it becomes:

𝜎 𝑉𝑀
= √
1
2
(( 𝜎 𝑥𝑥
− 𝜎 𝑦𝑦
)
2
+ ( 𝜎 𝑥𝑥
)
2
+ ( 𝜎 𝑦𝑦
)
2
) + 3( 𝜏 𝑥𝑦
2
) , (4-9)
Where 𝜎 𝑥𝑥
, 𝜎 𝑦𝑦
, and 𝜎 𝑧𝑧
are orthogonal normal stresses in x, y, and z directions respectively, and
𝜏 𝑥𝑦
, 𝜏 𝑥𝑧
, and 𝜏 𝑦𝑧
are orthogonal shear stresses in x, y, and z directions. The result for low resolution
and high resolution elements is presented at the following parts.
4.2.1 Effect of Resolution
Coarse meshing of the model decreases the accuracy of calculations. On the other hand, the
application of high resolution modeling is limited due to the time consuming computation. In this
section, we studied the effect of resolution on stress field modeling. In the first case, we used
coarse grids to calculate stress fields and fracture propagation in a model. The problem is solved
using 11 by 6 element meshes. It divides the plate into 66 cells.
For the next test, model of high resolution 111 by 56 elements was designed, and a stress
field is calculated. The result is shown in Figure 4.5.
As seen in Figure 4.5(a)-(d), the plate shows the linear behavior in response to the applied
stress. The stress field follows a similar trend, even though the applied load varies the degree of
magnitude. The system is designed symmetrically, but we observe an asymmetrical stress field in
upper and lower sides of the fracture. It is a problem in graphical representation of the stress field;
meaning that the stress at both sides of the crack is symmetrical. The stress concentration is
obviously at the crack tip, which is in accord with fracture mechanics model.
69

In Figure 4.5(e), fracture propagates through the plate and extends to the right side of it.
That indicates that the stress at the fracture tip has exceeded the maximum principal stress, and

(a)

(b)

(c)

(d)

(e)
Figure 4.4: Stress field distribution in a low resolution rectangular plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by factors of
10,000, 1,000, 100, 10, and 2, respectively.

the fracturing criteria has been met. Therefore, element at the crack tip becomes totally damaged
and can be removed from the analysis.
70

(a)

(b)

(c)

(d)

(e)
Figure 4.5: Stress field distribution in a high resolution rectangular plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by factors of
10,000, 1,000, 100, 10, and 2, respectively.

In contrast to Figure 4.4, the boundary of the stress field is smooth, and a bean-shaped
stress field can be observed around the crack tip. Stress fields for Figure 4.5(a)-(d) obey the
elasticity law. In Figure 4.5(d), the stress around the fracture tip looks different because of the
graphical representation of the stress in the plate. Fracture growth inside the model is shown in
Figure 4.5(e). The boundary effect changes the shape of the stress field in the fracture tip. It no
longer resembles a bean shape. The area in front of the crack resists fracture propagation, hence
71

we see the bean shaped stress field. If the area disappears as a result of fracture propagation, the
stress field shape changes, and maximum stress is transferred to the right wall.
Comparing high and low resolution models provides some insight into the effects of
resolution and grid size. In low resolution model, material failed after load exceeds 5 MPa.
However, in high resolution model, the structure was not damaged completely. We concluded that
grid size should be selected carefully to avoid an inaccurate modeling of the fracturing process.
We need to emphasize that the model is developed to simulate stress field under the given
boundary conditions. Therefore, the results do not necessary reflect the reality of a reservoir
condition. One of the main features we noticed is the stress field at the opposite side of the plate.
There is a trace of boundary effect at the right edge, which may not be observed in the reservoir.
To diminish such effects, the infinite elements need to be used. Abaqus is a powerful tool for
modeling fracture initiation and growth. However, it has limitations that one should keep in mind.
A case in point -- we performed similar analysis considering 100 by 50 elements. The result is
shown in Figure 4.6.
The stress field in Figure 4.6 is distributed randomly, which suggests an invalid calculation
of the stress in the model because the crack is located at the element boundary. Abaqus fails to
model that circumstance. Special attention should be paid to the initialization of the model to avoid
such situations.

72

Figure 4.6: The Von Mises stress distribution for 100 by 50 meshing scheme.

4.2.2 Effect of Plane Stress and Plane Strain Condition
Hydraulic fracturing simulation demands high capacity computation space, as well as well-defined
models. Having the luxury of time and space enables one to increase the accuracy of the simulation.
We are somewhat limited in resources, hence the need to make assumptions to downsize the
problem. One of the early methods for reducing the size of the problem was suggested in the KGD
approach. It is assumed that the fracture height was much greater than the fracture length.
Therefore, the effects of strain is ignored in a vertical direction. Many other researchers adapted
such assumption to solve the fracture propagation problem in 2D models (Adachi and Detournay,
2008, Garagash, 2006, Bao et al, 2015). In this part, we simulate the fracture propagation in a two
dimensional model, based on an assumption of plane strain, and compare it with the result of plane
stress that was shown in previous part.
Low resolution and high resolution models were tested and Von Mises stress was
calculated. The results are shown in Figure 4.7.
73

Stress field distribution in Figure 4.7(a)-(d) follows a similar pattern, which can be
attributed to the elastic behavior of the material before damage onset. As soon as stress exceeds
the maximum principal stress, damage occurs at the element containing the crack tip. The theory
of elasticity is then no longer valid. As seen in Figure 4.7(e), the fracture propagates along the
plate and divides the material into two parts.

(a)

(b)

(c)

(d)

(e)
Figure 4.7: Stress field distribution in a high resolution rectangular plate in a plane strain condition. Fluid
Pa. The width of fractures are
6
Pa, and (e) 5×10
6
Pa, (d) 10
5
Pa, (c) 10
4
Pa, (b) 10
3
pressure is (a) 10
magnified by factors 2.

74

The element ahead of the crack tip for the case of 1000 KPa fluid pressure was selected as
a benchmark to compare the plane stress and plane strain assumptions. The Von Mises stress for
plane stress and plane strain is 2114 KPa and 1290 KPa, respectively. Therefore, in similar
conditions, if we assume that Von Mises is a valid criteria with which to compare the cases, crack
propagation in the plane stress condition occurs in less stress than the plane strain. However, stress
in x and y directions are close in both cases. σxx and σyy are 1249 KPa and 2441 KPa in plane strain,
and are 1245 KPa and 2464 KPa in plane stress. The considerable difference between the two is
explained by looking at the σ zz. For plane stress, σzz is zero but for plane strain it is 1113 KPa. σzz
is included in the calculation of the Von Mises stress. Therefore, σzz can significantly change Von
Mises stress for plane strain and causes the difference between Von Mises stress between plane
stress and plane strain conditions.
The analysis showed that plane stress/strain assumptions could affect the model response
to the applied stress. Therefore, one should investigate to ascertain that the assumption is in accord
with the reality before making any assumptions.
4.2.3 Effect of Type of Analysis
The process of hydraulic fracturing propagation can be complex and scale dependent. Therefore,
it might be possible to have a wide range of processes in which fractures show different behavior.
In this section, we compare the implicit dynamic response of rock in a plane strain condition, to
that of a static response. Several tests were performed to study the effect of tensile load on stress
fields and fracture propagation. Similar conditions were taken into account; the only change is the
material response. In this example the load was increased to 10 MPa in 1, 10, and 100 seconds,
respectively. The material response for various load magnitudes for the case of 10MPa/100s load
rate is shown in Figure 4.8.
75

Figure 4.8(d) shows a different pattern from the Figure 4.8(a)-(c). To ensure that the
material still obeys the elasticity law, Figure 4.9 presents the same stress field, but with a different
legend. The stress distribution in the plate remains identical before damage starts to expand as it
can be seen in the figure. Therefore, there is a linear relation between stress and strain. Otherwise,
plasticity could affect the stress response of elements which leads to change in stress distribution.

(a)

(b)

(c)

(d)

(e)
Figure 4.8: Stress field distribution in a high resolution rectangular plate for 10MPa/100s load rate. Fluid
pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa, (d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are
magnified by factors of 10,000, 1,000, 100, 10, and 2, respectively.

76

Static and dynamic analysis of the model demonstrates the influence of inertia in fracture
propagation. When the stress is less than maximum horizontal stress, behavior of the fracture is
the same as depicted in Figure 4.5(a)-(d) and Figure 4.8(a)-(d). Once the crack growth occurs, the
stress field around the fracture tip will change. Figure 4.5(e) and Figure 4.8(e) depict a similar
process from a different approach. As seen from the Figures, fracture size is different at the end of
analysis as a result of the solution strategy. Incorporating any effect of inertia generates a different
result. It means that in the modeling of hydraulic fracturing, strain rate should be included in the
calculations, if applicable. The inertia might be an important factor during the early stages of
fracking, in which high pressure fluid is injected into a small perforation. The effect of inertia is
well documented for rock and metal (Doan and Gary, 2009, Rajendran et al, 1991). However, more
research is required for hydraulic fracturing.

Figure 4.9: Stress field distribution in a rectangular plate. The fracture is magnified by a factor of 2.
Legend is modified for improved representation.

Quantitative analysis of the model can give us valuable information on how inertia might
influence the system response. At the start of this analysis, we take the element in front of crack
77

tip as a reference to compare Von Mises stress between dynamic and static analysis. In static, the
crack tip intersects the element after 4 steps. It needs to be pointed out that the step in static analysis
is different than step in the dynamic. Step in static means the incremental increase in load.
However, in dynamic, step means incremental increase in time. The time-step for the dynamic load
is 21.38 seconds, which is very close to the step in the static analysis (20.38/100×5=1.069). The
Von Mises stress for the dynamic analysis is 1227.2 KPa, and for static analysis is 1216.6 KPa,
which shows they are fairly close. We performed similar analysis for an element near the right
wall. The step for static was close to the time of the dynamic model, and shows similar trend for
both cases. Small differences between the two can be attributed to the minor effect of the inertia.
A similar approach was performed in a case with a higher stress rate. In that example, the
load rate was increased from 10 MPa/100 s to 10 MPa/10 s and 10 MPa/1 s, so as to investigate
the effect of inertia on the analysis.
The dynamic approach includes the effect of inertia, which means that the time that load is
applied onto the object is integrated into the calculation. Therefore, the effect of the load time
period may have influence on the final results. To examine such a relation between time and crack
growth, the last analysis has been repeated with a change in applied time. The result is shown in
Figure 4.10.
In terms of fracture configuration, both Figure 4.10(e) and Figure 4.8(e) show analogous
behavior of fracture propagation. An important result of this analysis is that in the similar time
ratio between static, quasi-static modeling for 10 s, and quasi-static modeling for 100 s, the model
shows comparable behavior with very close stress. We performed the analysis for a 1 s load. The
result is shown in Figure 4.11.
78

Figure 4.10: Stress field distribution in a rectangular plate for load rate of 10 MPa/10 s. The width of
fracture is magnified by a factor of 2.

Figure 4.11: Stress field distribution in a rectangular plate for load rate of 10 MPa/1 s. The width of
fracture is magnified by a factor of 2.

Decreasing loading time to 1 s does not change the response of the plate significantly. That
can be attributed to the fact that the quasi-static model is a sequence of static analysis -- the strain
rate is still small and cannot exhibits the phenomena that may occur when the fracture propagation
speed is greater than seismic speed. However, one advantage of the implicit dynamic is that it does
79

not fail once the fracture reaches in the neighborhood of object boundary. Therefore, it does
provide a useful modeling tool for fractures.
As we are now familiar with the difference between quasi-static and static analysis, we will
continue our calculations based on the quasi-static assumption for convergence purposes.
4.2.4 Effect of Model Shape
Selection of accurate boundary conditions is one the most important tasks of the modeling process.
Boundary conditions should match real conditions for the experiment or field test, or must be close
enough that the effect to the final result is insignificant. In this part, we have changed the width of
the model to study the effect of model size on the stress field, and stress at the edges. In the first
scenario, a coarse model was created, and tensile load applied from both the top and bottom. Figure
4.3 shows the schematic of the model.
As seen in the Figure, a 25 cm initial fracture was placed at the center of the plate. The
purpose was to investigate the effect of distance between the boundary and fracture tip. In the
rectangular plate, the distance was 0.375 m. In the square, the distance is 0.75 m, which is twice
that of the previous example.
The load is applied on the square from bottom and top. It creates a tensile stress in the
elements. The result of analysis is shown in Error! Reference source not found..
In contrast to the rectangular plate (1 m by 0.5 m), the square plate does not obey the
elasticity law at 1 MPa. In other words, the principal stress at the fracture tip exceeds the maximum
allowable stress, and damage occurs. Two factors contribute to early fracture propagation: (1) the
80

area of applied load is increased by a factor of two, meaning that the energy concentration at the
fracture tip is going to increase, and (2) the size of the fracture increased by a factor of two as well.
In LEFM, the stress intensity factor at the fracture tip has a direct relation with the square
root of the flaw size (Roylance, 2001). This relation in another form is:
K
I
= σ
∞
√ πa, (4-10)

Where 𝐾 𝐼 is stress intensity factor, and a is the flaw size. As the size of the fracture increases, the
stress intensity factor increases. Therefore, in our case, the possibility of fracture propagation
increases with an increase in the crack size. On the other hand, increase in width of a model is not
in favor of crack growth and increases the amount of energy needed for fracture initiation. The
overall effect of these three phenomena is to reduce the stress needed for fracture propagation,
hence causes the crack to grow at lower load.
Also, the stress field distribution along the fracture path can provide useful information
about fracture behavior. Figure 4.13 shows a plot of stress along the fracture path for the load of
100 KPa.
81

(a)

(b)

(c)

(d)

(e)
Figure 4.12: Stress field distribution in a square plate. Fluid pressure is (a) 10
3
Pa, (b) 10
4
Pa, (c) 10
5
Pa,
(d) 10
6
Pa, and (e) 5×10
6
Pa. The width of fractures are magnified by factors of 10,000, 1,000, 100, 10,
and 2, respectively.

82

In Figure 4.13, we studied effect of model width size on the stress distribution in the model.
The region ahead of crack, between the crack tip and the right edge is selected for the analysis. It
can be seen that for a similar load, the smaller plate shows lower stress at the fracture tip. It
consequently delays the fracture initiation and propagation. In the Figure 4.13(b), Von Mises stress
is plotted, versus the normalized distance of the fracture tip from the boundary at the right edge.
The stress distribution patterns are similar, but not identical. In both sides, near the fracture tip and
at the right edge, the stress is higher in wider plate. However, stress is higher in the middle part of
the narrower plate.

83

Figure 4.13: Von Mises stress distribution before fracture, (a) the x-axis represents the true distance from
fracture tip, (b) the x-axis represents the normalized distance from fracture tip.

4.2.5 Effect of Initial Fracture Location
In hydraulic fracturing, fluid-driven fractures grow in a reservoir. The reservoir stresses act
upon the fracture, and limit its propagation. The fluid pressure should be sufficiently high to
dominate the stresses and break the rock. Sometimes there is a network of natural fractures in the
0.00E+00
5.00E+04
1.00E+05
1.50E+05
2.00E+05
2.50E+05
3.00E+05
3.50E+05
4.00E+05
4.50E+05
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Von Mises Stress (MPa)
Distance from tip (m)
Von Mises Stress for Rectangulare and Square
Plates
Von Mises Stress For Rectangular Plate
Von Mises Stress For Square Plate
0.00E+00
5.00E+04
1.00E+05
1.50E+05
2.00E+05
2.50E+05
3.00E+05
3.50E+05
4.00E+05
4.50E+05
0 0.2 0.4 0.6 0.8 1
Von Mises Stress (MPa)
Normalized distance from tip (m)
Von Mises Stress for Rectangulare and Square
Plates(normalized distance)
Von Mises Stress For Rectangular Plate
Von Mises Stress For Square Plate
(a)
(b)
84

rock which need to be opened in order for fracture propagation to begin. Interaction of the
pressures inside the fracture, reservoir stresses, and the presence of any natural fractures determine
the pattern of the fracturing propagation. Leaving the latter to the next section, we will now focus
on the effect of the reservoir stress field on fracture propagation. The results may shed more light
on the process.
As one sees in Figure 4.3, compressive loads are applied from the top and bottom edges of
the model. Pressure inside the fracture increases, and applies load on the fracture wall. One
important point to mention is the step design. During the hydraulic fracturing process, the reservoir
stress is already present, and fluid pressure inside the fracture increases. Similar thinking should
be applied to step design, in that the fluid pressure and in-situ stresses should not increase
simultaneously but sequentially. First, the load should increase to reach the reservoir level; then
the fluid pressure should increase until fracture propagation occurs.
Four cases are studied to investigate the effects of pressure inside the fracture, as well as
the location of the initial fracture with respect to the boundary. Some conclusions from the previous
parts are also applicable here. For example, model size could be a factor when performing the
sensitivity analysis on the position of a fracture inside the plate.
We assume that the fracture is 0.5 m, 0.34 m, 0.25 m, and 0.16 m from the top. These
numbers are 1/2, 1/3, 1/4, and 1/6 of the fracture height, respectively. We started the analysis with
a basic test in which the fracture was placed at the center of the crack. Then, we changed the
position of the fracture, and reported the results.
A square plate was created and a fracture placed at the center of the model. The fracture
was filled with fluid that was given constant pressure. Abaqus does not update the pressure load
85

inside the fracture in each step. It is equivalent to the assumption that the propagation occurs so
fast that the fluid does not move within the fracture during analysis. The result of pressure buildup
in the fracture is given in Figure 4.14.

(a)

(b)

(c)

(d)
Figure 4.14: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid pressure
is (a) 10
6
Pa, (b) 2×10
6
Pa, (c) 3×10
6
Pa, and (d) 3.308×10
6
Pa. The width of fractures are magnified by
factors of 10,000, 100, 100, and 20, respectively. Initial fracture is placed 50 cm from the top.

Error! Reference source not found. and Figure 4.14 both model the same material, but
under different load conditions. In contrast to Error! Reference source not found., in the case of
the model examined here, the stress distribution in the vicinity of a fracture tip does not form a
86

bean shaped figure. Instead, it leans toward the fracture, and generates a complicated pattern.
Figure 4.14(b) and (c) are the best examples of such deformation. In Figure 4.14(c), a secondary
bean shape is formed, which does not agree with the LEFM principles. The secondary bean is a
result of invariable location pressure, which was applied on the initial fracture and did not progress.
The material under tensile load fails once the stress is 2631 KPa (Error! Reference source
not found.). Here, the fluid pressure is 3308 KPa, which is 677 KPa higher. The salient point is
that there are two loads applied from the top of the plate, and they need to be dominate to be able
to grow a fracture in the plate. If we include these loads and subtract 1000 KPa from 3308 KPa,
we notice that the net pressure is around 2308 KPa, which is lower than the pressure for the model
of under tensile load. That can be attributed to the complex stress field around the fracture, and a
non-linear relation between stress and load.
We have repeated the last analysis, except we assume an initial fracture is located 0.34 m
from the top edge. The result is shown in Figure 4.15.
Figure 4.15(a) shows the only case in which no fracture propagation occurred. With an
increase in fluid pressure, the damage progresses and the crack starts to grow, as shown in Figure
4.15(b) and Figure 4.15(c). Once the fracture reaches the right edge, at fluid pressure of 3204 KPa,
the plate breaks. Comparing Figure 4.14(d) to Figure 4.15(d), one notices that lower pressure is
needed in the latter in order to break the plate. That might be attributed to the location of the initial
crack. Effects of the crack location can change the stress intensity factor. A study by Sha and Yang
(1998) showed the change in the stress intensity factor as a function of the relative location of a
crack in a rectangular plate. The qualitative comparison of last two cases offers a similar result.
87

(a)

(b)

(c)

(d)
Figure 4.15: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid pressure
is (a) 10
6
Pa, (b) 2×10
6
Pa, (c) 3×10
6
Pa, and (d) 3.204×10
6
Pa. The width of fractures are magnified by
factors of 10,000, 100, 100, and 20, respectively. Initial fracture is placed 34 cm from the top.

In Figure 4.16, the initial fracture is placed 0.25 m from the top edge, and pressure is
applied inside the fracture.
88

(a)

(b)

(c)
Figure 4.16: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid pressure
is (a) 10
6
Pa, (b) 2×10
6
Pa, and (c) 2.844×10
6
Pa. The width of fractures are magnified by factors of 100,
100, and 10, respectively. Initial fracture is placed 25 cm from the top.

As it can be seen in Figure 4.16, the breaking pressure is lower than in the last two cases.
It suggests that the failure is a strong function of the relative position of the fracture. The common
feature among the Figure 4.16(a)-(c) is the rotation of the plate. As discussed earlier, the plate is
allowed to rotate or move. In this example, the inequality of acting forces causes the plate rotation,
which is magnified by a factor of 100 in Figure 4.16(a) and Figure 4.16(b). Also, the bean shaped
stress field is no longer symmetrical because of the asymmetrical shape of the configuration.
89

As shown in previous sections, the stress intensity factor decreases if the fracture is placed
close to the plate boundary. In this case, we position the fracture very close to the boundary, and
then compare the results to last three cases. The result is shown in Figure 4.17.
Figure 4.17 shows the different patterns of fracture propagation relative to the previously
studied cases. In the last cases, we saw fracture path evolved in the same trend. However, fracture
in the Figure 4.17(c) deviated from fracture in Figure 4.17(b) and showed a different course of
propagation, which was attributed to the resolution dependency of the problem. The time-step
selection varies from case to case and it does influence the propagation path of a fracture.
Evidently, if the step size is small enough, the crack should indeed follow the path if the load
increases. Also, the higher mesh resolution could improve the accuracy of the modeling process.
Drawing a general conclusion requires more study on analysis of step and mesh resolution.
As argued before, the position of a fracture modifies the breaking point of the plate.
Quantitative studies of such an influence have been performed on the crack initiation pressure, as
well as rock failure pressure. Figure 4.18 and Figure 4.19 show two plots of such effects for
different fracture positions.
There is a strong correlation between the position of the crack and the minimum fluid
pressure for fracture initiation and integrity failure, as noticed in both figures. Decreasing the
distance to the top edge causes a reduction in the fluid pressure needed for cracking. In both cases,
the second degree polynomial function was used to fit the data. In the second case, the model
matched the data perfectly. Analytical analysis should be conducted to assure the presence of such
relationships.
90

(a)

(b)

(c)
Figure 4.17: Stress field distribution in a square plate under 1000 KPa compressive load. Fluid pressure
is (a) 10
6
Pa, (b) 2×10
6
Pa, and (c) 2.27×10
6
Pa. The width of fractures are magnified by factors of 100,
100, and 10, respectively. Initial fracture is placed 17 cm from the top.

91

Figure 4.18: Fracture initiation pressure in different fluid pressures

Figure 4.19: Plate failure pressure for different fluid pressures

0
0.5
1
1.5
2
2.5
0 10 20 30 40 50 60
Stress, MPa
Distance of Fracture from Edge of Plate, cm
0
0.5
1
1.5
2
2.5
3
3.5
4
0 10 20 30 40 50 60
Stress, MPa
Distance of Fracture from Edge of Plate, cm
92

In this section, various analyses were performed to study the effects of diverse parameters
on single fracture propagation. Mesh resolution, size of the model, crack position, load type,
analysis approach, and solution strategy were all discussed in detail. Next we are going to examine
the case were there are two cracks in the model.
4.3 Multiple Crack Propagation
Hydraulic fracturing in unconventional reservoirs can be considered successful when it creates a
network of interconnected fractures. In the process, the fracture propagates, and its pattern is
impacted by the presence of natural fractures as well as hydraulic fractures from other stages. The
hydraulic fractures meet the natural ones, and depending upon the type of natural fractures and
stress fields, they may cross those fractures, induce remote shear slippage, or intersect and continue
propagation through the natural fractures.
Accurate modeling of the interaction between natural fractures and hydraulic fractures and
two hydraulic fractures could be performed in presence of sufficient data. In this section, we will
model multiple fracture propagation in a simple model, as further data and a more advanced model
are needed to enable us to accurately model the entire process.
It is important to have a knowledge of reservoir properties to model the fracture interaction.
In the reservoir, the compressive forces are acting on the fractures, and limit their propagation.
The fractures are filled with the fluids and fluid pressure gives them energy to break the rock and
propagate. We will endeavor to implement the same environment, but in a laboratory scale, so that
the infinite boundary condition will not be applied. And in some cases, fluid can migrate from a
hydraulic fracture to the natural one, which makes it very realistic. To the best of author’s
knowledge, no similar analysis has been conducted in the Abaqus environment.
93

In this section, first, interaction of two parallel fracture is investigated. We start with a
rectangular plate and continue with a square plate. Then, interaction of natural and hydraulic
fractures is studied for both rectangular and square plates.
4.3.1 Effect of Fracture Spacing In Rectangular Plate
In this part, a rectangular plate will be used to model parallel fracture propagation. Two
compressive loads are applied to limit fracture propagation. The schematic of the model is
illustrated in Figure 4.20.

Figure 4.20: Plate is subjected to the two compressive loads acting from the top and bottom of the plate.
Two fractures are placed at the plate and are separated by a distance of h. The fluid pressure force the
fractures to open and propagate.

The operation is performed in two stages: (1) the compressive loads increase up to the
maximum level, and (2) the fluid pressure increases gradually to reach the maximum level. The
maximum compressive load is equal to 1000 KPa. No maximum level was set for the fluid
pressure. Figure 4.21 shows the step-by-step process of fracture propagation in the plate.
12.5 cm
100 cm
50 cm
h
P
94

(a)

(b)

(c)

(d)

(e)
Figure 4.21: Stress field distribution in a rectangular plate. Distance between two fractures is 10 cm.
Fluid pressure is (a) 2.677×10
6
Pa, (b) 3.628×10
6
Pa, (c) 3.785×10
6
Pa, (d) 3.878×10
6
Pa, and (e)
3.879×10
6
Pa. The width of fractures are magnified by factor 10.

Fluid pressure inside the fracture applies load to the fracture wall. However, the fracture
wall opening is constrained by the compressive stresses exerted by the load on the edges, as well
as effects of the other fracture. The mild diversion of two fractures is observed in Figure 4.21(a)-
(c). But, in Figure 4.21(d), there is a sudden change of pattern and we see in Figure 4.21(e), the
two fractures propagate in different directions.
95

In next, the distance between fractures is increased to 30 cm. We see in last test that
distance does influence fracture growth when it is under tensile stress. This test was conducted to
illustrate compressive load. The result is in Figure 4.22.

(a)

(b)

(c)

(d)

(e)
Figure 4.22: Stress field distribution in a rectangular plate. Distance between two fractures is 30 cm.
Fluid pressure is (a) 2.183×10
6
Pa, (b) 3.174×10
6
Pa, (c) 3.614×10
6
Pa, (d) 3.723×10
6
Pa, and (e)
3.735×10
6
Pa. The width of fractures are magnified by factor 10.

It catches one’s attention to realize that the pattern observed in under compressive load and
fluid driven fracture is similar to that under tensile load available in Figure B.4. Both show the
same material and model; the only difference is the boundary condition. We can see that boundary
96

condition does not change the pattern of fracture growth in the laboratory scale. However, in the
reservoir scale model, we may expect different results. Therefore, the realistic input data should
be used to generate a reliable results for the field purposes.
In Figure 4.22, two low stress areas are generated, possibly because Abaqus does not
update the fracture load during analysis.
For the next test, the distance between two fractures is increased to 50 cm. The result is
seen in Figure 4.23. The figure shows the step-by-step propagation of a fracture in the plate. In
contrast to last test, this fracture does not follow a straight path. It diverts slightly and continues
directly toward the boundary. In the previous case, the mechanism of fracture propagation was
completely different. The fracture grew due to an external force; thus, being close to the boundary
was a parameter. On the contrary, an internal load caused the fracture growth in this scenario,
which created a completely different pattern.
It is important to realize the limitations of software, and be able to reject results that do not
follow the physical pattern. Figure 4.23(d)-(e) show that the fracture is not symmetrical. There is
a different pattern in the upper fracture as compared to the lower one, which might be
mathematically true. It means that the asymmetrical meshing of the model has created an
asymmetrical solution. But in reality, because of symmetrical load, boundary condition, and model
configuration, such features do not exist. It could be a numerical artifact that should be avoided by
changing the meshing approach.
97

(a)

(b)

(c)

(d)

(e)
Figure 4.23: Stress field distribution in a rectangular plate. Distance between two fractures is 50 cm.
Fluid pressure is (a) 2.130×10
6
Pa, (b) 3.146×10
6
Pa, (c) 3.3127×10
6
Pa, (d) 3.3128×10
6
Pa, and (e)
3.3129×10
6
Pa. The width of fractures are magnified by factor 10.

4.3.2 Effect of Fracture Spacing In Square Plate
The boundary effect is an important parameter when modeling of hydraulic fracturing is
conducted in a small plate. The larger models have the advantage of a lower boundary effect.
Except for cases of infinite elements, the boundary effect could change the response of the system.
A large plate is used to reduce the effect of boundary in a model under compressive load. In this
98

example, the effect of the load size, and distance between two fractures were studied. Schematic
of the two fractures is in Figure 4.24.

Figure 4.24: Plate is subjected to the two compressive loads acting from the top and bottom of the plate.
Initial fractures are placed at the left part of the plate. Distance between the fractures is shown by h.

Figure 4.24 shows the relative position of two fractures; the distance between them is a
variable, and is updated for different scenarios. In the first instance, the distance between the two
is equal to 5 cm. The result of that analysis is shown in Figure 4.25.
Figure 4.25 illustrates how an increase in pressure can lead to fracture propagation. It was
seen in Figure 4.21 that when the distance between two fractures was 5 cm, they diverged, and the
material failed. The only difference here is that there was not the same sharp divergence as in
Figure 4.25. The fracture propagation pattern resembles that of tensile load which is available in
Appendix B. In that case, there is divergence at the beginning, but after passing the middle of the
plate, they once again became closer. Such similarity suggests a relationship
25 cm
100 cm
100 cm h
P
99

(a)

(b)

(c)

(d)

(e)
Figure 4.25: Stress field distribution in a square plate. Distance between two fractures is 5 cm. Fluid
pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.728×10
6
Pa, and (e) 3.744×10
6
Pa. The
width of fractures at (a) to (d) are magnified by factors 10 and 2 for figure (e).

100

between the size of the model and the type of pattern that is shown. To recap, the models behave
the same under a certain width of the plate. Fracture behavior at Figure 4.21(e) was an exceptional
case in which the divergence increased sharply, and led to total damage.
Following, the distance between cracks was increased to 10 cm. The other features were
not unlike the previous model. The result of fracture propagation is shown in Figure 4.26.
Figure 4.26(a) and (b) show that the plate remains undamaged at the beginning. However,
an increase in pressure leads to fracture propagation and total damage, as represented by Figure
4.26(e). The similarity between Figure 4.26 and Figure 4.25 is worth mentioning. The increase in
distance had an insignificant effect on the fracturing process. To qualify the effect of distance, a
new test should be designed.
The last two tests demonstrated that the distance between fractures has only a minor effect
on the fracturing pattern. This test was instigated to examine the hypothesis. Figure 4.27 shows
the fracture propagation when the h is 20 cm.
Figure 4.27(a)-(d) represent similar patterns to those shown in Figure 4.26(a)-(d), which is
that increasing the distance did not have a considerable effect on the pattern. However, the
minimum fluid pressure to damage the material declined by almost 5%, which is noteworthy. Also,
the fractures do not get closer, which is a direct effect of the top edge boundary condition.
It is important to note that they are other mechanisms that may occur which numerical
model cannot capture due to its limitation. However, we need to identify such mechanisms. For
this purpose, we will have a closer look on the fracture behavior when the fractures are 20 cm
apart. The model is shown in Figure 4.28.
101

(a)

(b)

(c)

(d)

(e)
Figure 4.26: Stress field distribution in a square plate. Distance between two fractures is 10 cm. Fluid
pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.480×10
6
Pa, and (e) 3.734×10
6
Pa. The
width of fractures are magnified by factors 10.

102

(a)

(b)

(c)

(d)

(e)
Figure 4.27: Stress field distribution in a square plate. Distance between two fractures is 20 cm. Fluid
pressure is (a) 1.0×10
6
Pa, (b) 2.0×10
6
Pa, (c) 3.0×10
6
Pa, (d) 3.551×10
6
Pa, and (e) 3.583×10
6
Pa. The
width of fractures are magnified by factors 10.

103

Figure 4.28: Fracture propagation in a plate in the direction of existing fracture

Generally, fracture propagation should occur in the direction that has maximum principal
stress. Figure 4.28 shows a case in which fracture branching is possible; the arrow indicates the
maximum principal stress. Due to limitations of XFEM in Abaqus, modeling of branching cannot
happen.
4.3.3 Hydraulic and Natural Fracture Interaction
Modeling of hydraulic and natural fracture interaction is different than other cases that have been
studied. In this process, hydraulic fracture is pressurized and propagates due to tensile load on the
fracture wall. However, there is not any fluid inside natural fracture therefore propagation in
natural fracture occurs if it becomes activated. The interaction of two fracture was modeled using
rectangular and square plate. For rectangular plate, no sensitivity analysis on the distance of the
fractures is performed since the model size is too small to draw a reliable conclusion. But
sensitivity analysis is carried out for square plate.
For the rectangular plate, the result of the test is shown in Figure 4.29.
104

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)
Figure 4.29: Stress field distribution in a rectangular plate. Distance between two fractures is 33.75 cm.
Fluid pressure is (a) 2.34×10
6
Pa, (b) 2.35×10
6
Pa, (c) 2.37×10
6
Pa, (d) 2.381×10
6
Pa, (e) 2.382×10
6
Pa,
(f) 2.383×10
6
Pa, (g) 2.384×10
6
Pa, (h) 2.385×10
6
Pa, and (i) 2.390×10
6
Pa. The width of fractures are
magnified by factors 10.
105

Figure 4.29 exhibits one of the major phenomena in hydraulic fracturing. When hydraulic
fractures grow, they meet natural fractures. The natural ones can either arrest the fracture
propagation, or let it cross over the natural fractures, or activate the natural fractures, or even
completely divert the hydraulic fractures. These mechanisms have been widely studied;
experimental and analytical models were developed to describe the behavior (Blanton, 1982,
Potluri et al, 2005). However, only a handful of studies were able to quantify the process. Figure
4.29 shows how fractures may interact in a laboratory setting. Only hydraulic fractures propagate
in Figure 4.29(a)-(d). But, when the hydraulic fracture is close enough to the natural one, it changes
the stress field around the crack tips, and thereby remotely activates the natural fracture. Figure
4.30 shows the shear stress in the model at pressure 2340 KPa.

Figure 4.30: Shear stress around the crack tips during hydraulic fracture growth.

As seen in the Figure, shear stress is concentrated near the natural fracture tip, even though
fluid has not yet reached there. The remote interaction of hydraulic fractures with natural ones
trigger the opening of the natural fractures. It could cause microseismicity in the rock even though
106

no new fracture is added to the network of natural fractures. Such mechanisms can lead to false
understanding of the inter-connected fracture networks.
In addition to knowing that fractures can grow remotely, we need to consider the effect of
boundary condition. We observe the interaction of fractures, but that does not necessarily reflect
real field behavior. Natural fracture propagation is in an opposite direction to that of hydraulic
fracture, which is unlikely in the reservoir. That behavior can be attributed to the boundary
condition effect. Natural fracture is close to the boundary and therefore it is impacted.
For the square model, effect of minimum and maximum horizontal stress are taken into
account. Maximum and minimum horizontal stresses play an important role in the hydraulic
fracturing process. The minimum pressure to begin hydraulic fracturing is commonly assumed to
be equal to the minimum horizontal stress. Therefore, in a comprehensive study, the effect of in
situ stresses should be taken into account. Here, the previous models have been modified; the size
of the plate is 1 m by 1 m, and both stresses are acting on the model. The distance between fractures
is a variable, and subject to change. The study aims to model the behavior of intersecting fractures.
Schematic is shown in Figure 4.31.
Including the effect of minimum and maximum horizontal stresses is one of the advantages
of this model. However, it is a small plate, and the effect of boundary is not negligible.
The distance between two fractures is shown by d in Figure 4.31. The model is used to
study the effect of distance between two fractures on the plate response to an applied load. Three
cases were implemented to examine the effect of distances.
107

Figure 4.31: Plate is subjected to the four compressive loads acting from the four edges. Distance
between two fractures is shown by d.

For the first test, we assign distance to be 25 cm. However, in this model, the shortest
distance between two fractures is 12.5 cm. The distance between the central fracture and right edge
is 62.5 cm, which is the farthest distance compared to the other tests. Therefore, the effect of a
right boundary is minimum compared to other tests. Figure 4.32 shows the result.
In Figure 4.32(a)-(b), the red area represents the region that has the highest stress. It is
interesting to notice that the area is not in front of the fracture, but along the crack tip. That is an
indication of a presence of natural fracture. Comparing Figure 4.14(a) to Figure 4.32(a), the only
change in the model is addition on a natural fracture. Therefore, the effect of natural fractures can
be examined by studying stress distribution.
One of the issues that needs to be taken into account is that Abaqus does not update load
on the crack, nor does it include the load in the entire fracture region. Therefore, two high stress
regions are created which may be possible in case of high viscosity fluid injection.
25 cm
100 cm
100 cm
d
P
25 cm
108

Figure 4.32: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two fractures is
25 cm. The width of fractures are magnified by factor 10. The pressure inside fracture is 2.7 MPa.

In Figure 4.32(b), the low stress region emerges, which is a direct effect of natural fractures.
It prevents the hydraulic fracture from pushing the region behind the natural fracture. And as we
notice in Figure 4.32(c), the highest stressed node is no longer in the vicinity of the hydraulic
109

fracture, but is in the natural fracture damage zone. However, it does not stay there, but returns to
the hydraulic fracture as shown in Figure 4.32(d) and (e). Abaqus is not able to model crack
intersection hence, no fracture propagation occurs when the maximum stress is at the hydraulic
fracture tip.
For the next test, the natural fracture is shifted 12.5 cm to the right and is placed 50 cm
from the right edge. The fluid pressure increases and the fracture starts to grow. Figure 4.33
illustrates the fracture propagation and interaction of the two fractures.
If one compares Figure 4.32(a) to Figure 4.33(a), similar patterns can be observed. The
high stress region is along the crack tip and not in the process zone. However, Figure 4.33(b) shows
a different pattern. The high stress zone is following the crack tip, and is stuck in a fixed position.
That indicates that the initial load is not getting updated; therefore, a high pressure zone is created
in the location of the initial fracture tip.
One might mention that the observed pattern is quite similar to the example of the plate
being under tensile load shown in Appendix B. The main factor here is the effect of remote stress.
When the fracture is in the natural fracture area, it remotely causes the growth at the fracture tip.
When the pressure reaches a natural fracture, it propagates in a created pathway which causes it to
go in a direction parallel to the hydraulic fracture. In the previous case, the small fracture deviated
from the remotely created path, and changed course. Further investigation is needed to enable us
to accurately model the process.
For the final test, the fracture is put 25 cm from the right edge, which makes the distance
between the centers of hydraulic and natural fractures 62.5 cm. The effect of the boundary is
110

minimized compared to the prior cases. The result of constant pressure fluid injection is shown in
Figure 4.34.

Figure 4.33: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two fractures is
37.5 cm. The width of fractures are magnified by factor 10. The pressure inside fracture is 2.7 MPa.

111

Figure 4.34: Stress field distribution in a 1 m by 1 m rectangular plate. Distance between two fractures is
62.5 cm. The width of fractures are magnified by factor 10. The pressure inside fracture is 2.7 MPa.

A similar pattern is repeated in Figure 4.34(a). The initial fracture propagation is not
significantly influenced by the distance between cracks. As the hydraulic fracture get closer to the
natural fracture, the high stress element location changes. It is interesting to see that in this case,
112

the central element of the natural fracture (element 6189) is under the highest stress, as shown in
Figure 4.34(d). The element is not in the tip, so the fracture cannot grow from that element.
Increasing the stress at the natural fracture tip in Figure 4.34(e) causes natural fracture propagation.
Also, the hydraulic fracture is arrested at the natural fracture, and no propagation of that fracture
is likely to occur.
4.4 Conclusion Remarks
Complexity of the created fracture network in the reservoir depends on the interaction of both
hydraulic and natural fractures, as well as the stress field distribution near the fractured area. In
this chapter, we showed how a single fracture and two interacting fractures behave in different
conditions. Sensitivity analysis was performed on various parameters to have better insight on the
hydraulic fracturing process. For single fracture, effect of resolution, plane stress and plane strain
condition, analysis type, model shape and initial fracture location is investigated. It is shown that
material response is different when resolution of model increases. For this case, higher resolution
leads to higher resistance of material. Plane stress and plane strain condition however does not
change the material behavior but change the stress field which subsequently alter the material
response to an external load. Effect of inertia also studied in this Chapter. It is shown that inertia
does not affect the load magnitude for fracture initiation. It could be attributed to the fact that in
our analysis the loading time was high thus the momentum was not significant.
Two plates with different width are designed to study the model size effect. In our study,
we show that larger plate has higher resistance to the fracture propagation and higher pressure
needs to be applied to break the material. It is also shown that the proximity of fracture to the plate
boundary reduces the pressure needed for crack initiation.
113

Multiple fracture model is designed to study the interacting fractures. For the case one, the
interaction of propagating hydraulic fractures is modeled. The sensitivity analysis is performed on
the distance between fractures. It is demonstrated that the distance between two fractures plays a
key role for fracture initiation and propagation and the size of model does not significantly affect
the pattern of fracture growth. For the second case, interaction of natural fracture and hydraulic
fracturing is investigated. The model shows that there are two reasons for natural fracture
activation. It can be either opened through fluid transfer from hydraulic fracture to the natural
fracture or it can be activated through remote interaction of natural and hydraulic fracturing.
Simulation of this phenomena is carried out by using the developed algorithm with constant fluid
pressure injection condition.
In this Chapter, we have tested the algorithm to model complicated processes that are
encountered in hydraulic fracturing. The algorithm is capable of modeling these phenomena and
provides an insight into the hydraulic fracturing process. This Chapter paved a way for
comprehensive studying of hydraulic fracturing in the field scale. In the next Chapter, the
algorithm is tested for 3-D model. The assumption of constant pressure is relaxed and lubrication
theory is used to model fluid flow in the fracture.
114

Chapter 5 : THREE DIMENSIONAL SIMULATION OF HYDRAULIC
FRACTURING IN FIELD SCALE
In the previous chapter, we explained how Abaqus can be coupled with Matlab to model fracture
initiation and propagation in a simple 2D model. In this chapter, we will describe the results of
tests on a 3D model that was built, based upon a reservoir model, to simulate fracture behavior in
the field scale. The constant pressure assumption was lessened in the developed model. The fluid
flow has been incorporated into the computation by the integration of a fluid flow simulator and
solid deformation simulator.
A 3D rock model was constructed and analogous flow model created in order to model the
deformation of rock, and fluid pressure distribution. The Lagrangian and Eulerian approaches were
used to represent solid deformation and fluid flow, respectively. The pressure was transferred to
the rock in the form of load, and stress calculated in the solid modeling tool. The model was tested
in several cases and the results are presented in this chapter.
5.1 Introduction
Numerical simulation of hydraulic fracturing could offer useful insights into the interaction of
various mechanisms, and provide an estimation of hydraulic fracturing network, which can be later
modified using the production and seismic data. In order to gain maximum benefit from the
simulation results, a mathematical model should be equipped with an accurate rock and fluid
characterization model. Too often it is difficult to obtain such accurate models, thus simpler
approaches have been applied to extract a network of fractures. One major assumptions is that of
a 2D reservoir. Many reservoirs consist of multilayer formations with a wide range of properties,
115

which nullifies that assumption. Furthermore, there are number of features in the reservoir which
are limited to certain regions and do not extend throughout the reservoir layers. Consequently,
adding a third dimension is essential to the development of a realistic model. Our model is capable
of incorporating local features such as natural fractures in a 3D modeling which give it an
exceptional advantages over current hydraulic fracturing simulators.
Many attempts have been made for 3D simulation of hydraulic fracturing. The first
research on this topic was begun by Daneshy (1978) in which the experimental and theoretical
propagation of hydraulic fracturing was investigated for a multilayer sample. He concluded that
the strength of interface between layers determines hydraulic fracturing extension from the original
layer to the adjacent one. To avoid the complexity of 3D modeling of hydraulic fracturing, a pseudo
three dimensional model (P3DH) was introduced for the simulation of hydraulic fractures in
layered reservoirs (Settari and Cleary, 1984, Settari and Cleary, 1986). The model was based on
the one dimensional lateral propagation, and two dimensional vertical modeling of fractures. The
Finite Difference Method (FDM) was utilized to solve equations numerically. The model was able
to simulate a simple case of hydraulic fractures in natural fracture-free formations. One of the
advantages of the model is multiphase simulation of fluid flow during fracturing.
Limitation of computational capacities has always been an important issue. In one study,
full 3D modeling of hydraulic fracturing was selected as the first approach, but was terminated
due to the lack of computation space (Cleary et al, 1983). However, a new formulation was
presented to model hydraulic fracturing in a layered media. It considered 2D fluid flow in the rock,
and developed an integral method to link the opening of a fracture to the pressure.
In addition to multilayer reservoirs, deviation of fractures from the conventional planar
shape was a subject of research. Vandamme et al (1989) developed a numerical method based on
116

the Displacement Discontinuity (DD), which is capable of including nonplanar fracture
propagation. It assumed Newtonian fluid flow between two parallel plates. The important
deficiency of DD is its limitation for modeling of a heterogeneous environment.
Integral method was a popular approach to model hydraulic fracturing in a pseudo or full
manner. In that method, double or triple integral of the related equations is taken and the numerical
solution is calculated from the result of integration. Meyer (1989) utilized the method to solve the
hydraulic fracturing problem in a 3D environment. According to Meyer, the method falls into the
area between full 3D and pseudo 3D models. It was tested against a few other 3D models and
proved it could generate acceptable results.
As the computation power of computers increase, the application of full 3D modeling of
hydraulic fracturing becomes viable. Abaqus provides a useful platform for numerical analysis of
hydraulic fracturing. Zhang et al (2010) studied multistage hydraulic fracturing using Abaqus. A
3D model was constructed and fully coupled fluid flow, and rock deformation equations were
solved in order to calculate hydraulic fracturing parameters. FEM was used for numerical solution,
and Cohesive Element Method was utilized to describe the element behavior under the stress. Their
model was in good agreement with the reported field data. In recent studies by Zhou et al (2015),
a 2D hydraulic fracturing simulator was developed to predict fracturing in a 2D reservoir. Hybrid
of XFEM and Finite Volume Method (FVM) was utilized to model rock deformation and fluid
flow respectively. The model was tested against KGD model and it showed a reasonable match.
Microseismic monitoring of hydraulic fracturing becomes a vital tool for hydraulic
fracturing characterization and monitoring. The microseismic cloud commonly covers an area near
a horizontal well, which could not be modeled using previously developed methods. Conventional
models consider hydraulic fractures and neglect the presence of natural fractures. Therefore, the
117

conventional models cannot match microseismic cloud. To make the model more geologically
accurate, Discrete Fracture Network (DFN) is generated to represent the network of natural and
hydraulic fractures. There have been a number of studies in this area in recent years. Meyer et al
(2011) defined a DFN for a 3D model, whereby a simple stress field is taken into account and the
fluid flow is calculated. The model was tested for Marcellus Field. Effects of simultaneous fracture
propagation was quantified for a network of fractures.
Advances in computation power contributed to the modeling of more complex structures,
although it was still not sufficient. Extended Finite Element (XFEM) was introduced to alleviate
the issue of computing speed because of remeshing in FEM. As discussed previously, XFEM
explicitly includes the fractures, hence no remeshing of the model is needed for each time step.
Chen (2013) implemented the XFEM in Abaqus. The UEL user subroutine was employed to
change the element stiffness matrix. The model was tested in a small number of cases and was able
to successfully predict the dimensionless crack opening in comparison to the analytical solutions.
Zielonka et al (2014) have studied hydraulic fracture initiation and propagation using
Abaqus. XFEM was used for numerical formulation of the model and calculation of stress field.
Cohesive Element Method was applied to describe the element behavior. The model was compared
to the classical KGD and penny-shaped hydraulic fracture models and a good match was obtained
in the case of refined meshing.
Wang (2015) used Abaqus to model fracture propagation in the reservoir; a similar
approach to Zienlonka’s model was taken to calculate important parameters. The model was tested
for brittle and ductile rocks. Also, the Mohr-Columb theory of plasticity was utilized, in
conjunction with CZM, to increase the accuracy of the model.
118

Haddad et al (2015) employed Abaqus to perform a sensitivity analysis on the major
numerical and physical parameters. Various cases were studied and the effects of boundary
condition selection and number of fractures were investigated. A closed spaced multi-stage
hydraulic fracturing provided the best performance according to their analysis.
Khodabakhshnejad and Aminzadeh (2016) develop an algorithm by coupling fluid flow
and rock deformation. XFEM used to model rock deformation and coupled with an in0gouse FDM
for fluid flow modeling. The results show that the hydraulic fracture is penny-shape at the early-
times of operation and becomes PKN like fracture at transient and late time.
Based on the given discussion, three dimensional simulation of hydraulic fracturing is
necessary for accurate modeling of fractures. Multi layering of the reservoir and non-planar
hydraulic fracturing growth are the major factors to switch from conventional 2D hydraulic
fracturing simulation to an advanced 3D modeling. The recent studies could be carried out in 3D,
owing to the advances in numerical simulation, especially in geomechanical modeling. The
collaboration of Abaqus developers with oil industry pioneers originated the new refinements in
fracture modeling in reservoirs. It provides researchers a sophisticated tool with which to model
hydraulic fracturing initiation and propagation.
Abaqus was utilized for geomechanical modeling in 3D reservoir. Abaqus is widely used
for the numerical modeling of solid; we combined that capability with a fluid flow simulator that
we developed based on the Finite Difference Method (FDM). It enables us to design numerous
scenarios of fluid flow and relaxes some assumptions that were considered for developing fluid
flow models in Abaqus.
119

In this chapter, I will first present the characteristics of the model that has been constructed
for hydraulic fracturing simulation. Then, a few scenarios will be tested for the model and the
results will be detailed. The effects of few parameters such reservoir size and damage model on
fracturing initiation and propagation will also be analyzed.
5.2 Description of Model
Reservoir generally shows the properties of infinite domain. They are typically huge structures
and their geomechanical properties are controlled by far field stresses. Simulating such gigantic
models demands an extensive computational time and space which is usually not available. Taking
this into consideration, in this study, the smaller scale model is constructed to model hydraulic
fracturing in a reservoir, which is pictured below.

Figure 5.1: Schematic of the reservoir constructed for modeling of hydraulic fracturing.

Figure 5.1 shows the reservoir model that is constructed for the analysis of hydraulic
fracturing. One sees that the main layer is enclosed between two surrounding layers. The main
21 m
7 m
40 m
2 m
2 m
120

layer is a shale formation and the property is similar to the sample that has been tested by Busetti
et al (2012). Two surrounding layers are sand and are harder than the main layer. Considering their
higher strength, it is very unlikely that fractures will grow into adjacent sand layers. Even so, no
extra boundary condition is considered in the model that would limit fracture propagation into the
adjacent layers. That does not mean that fracture grows into the rocks, as the rock is sufficiently
hard most of the times and able to withstand the high stress at the interface with the main layer.
However, the physical properties of the rock could be changed, depending upon the field data, to
allow the fracture growth in the contiguous layer. In our case, the property of the rock is given in
Table 5.1.
Maximum principal stress determines the stress at which the damage starts in the rock.
Further explanation of the parameters of traction-separation law is available in Section 2.4.
After peak stress, the softening law determines the element behavior in the fracture tip. It
describes the element model before total degradation. Depending on the material type, softening
law can be linear, exponential, or user-defined. Considering the lack of data, we selected a simple
linear model to represent the material behavior. The parameters of the model can be easily
modified to adapt to the realistic data, if available. In the next section, we detail a sensitivity
analysis reservoir geometrical parameters and rock damage properties.

121

Table 5.1: Reservoir, geometrical, and damage properties of the model
Rock Properties
Shale Young’s Modulus (GPa) 5
Shale Poisson’s Ratio (fraction) 0.3
Sand Young’s Modulus (GPa) 27
Sand Poisson’s Ratio (fraction) 0.3
Geometrical Properties
Matrix Dimension (m×m×m) 40×21×7
Adjacent Layers Dimension (m×m×m) 40×2×7
Initial Fracture Dimension (x(m)×y(m)) 3×2
Damage Type
Damage Model Damage for Traction-Separation Law
Damage Criterion Maximum Principle Stress
Maximum Principal Stress for Shale (MPa) 2
Maximum Principal Stress for Sand (MPa) 5
Damage Evolution
Type Displacement
Softening Linear
Degradation Maximum
Mode Mixed mode
Displacement at the Failure for Shale 0.01
Displacement at the Failure for Sand 0.005

5.3 Parameter Sensitivity Analysis
Construction of a real size reservoir model for hydraulic fracturing simulation requires access to
accurate reservoir characterization data such as geomechanical property of rock and location of
natural fractures, although having such a model does not necessarily guarantee a reliable reservoir
simulation. Fine scale reservoir modeling needs to be performed in order to generate a plausible
122

result. Therefore, high performance computation space is required for this purpose. As neither are
available much of the time, a coarse scale model of a small scale reservoir was created. The
sensitivity analysis was performed on rock tensile strength, reservoir stress field, and reservoir
size. The effect of fluid parameters was neglected and the pressure assumed to be constant for the
entire fracture.
5.3.1 Sensitivity Analysis on Rock Strength and Stress Field
In theoretical models like PKN and KGD, the tensile strength is neglected and it is assumed that
there is no barrier for fracture growth except minimum horizontal stress. In this part, that
assumption is evaluated by including the effect of the rock strength. Also, two scenarios for stress
field are considered, and for each scenario, strength of rock is changed. The results are in Figure
5.2.

Figure 5.2: Effect of rock tensile strength on fracture initiation pressure

0
5
10
15
20
25
30
0 1 2 3 4 5 6
Fluid Pressure. MPa
Rock Tensile Strength, MPa
Non zero S1 and S2
Zero S1 and S2
123

This Figure shows the effect of rock tensile strength on the pressure in which fracture starts
to grow. Two cases are studied. In the first, it is assumed that the only active reservoir stress is
minimum horizontal stress, which acts on the fracture reservoir wall, and limits fracture initiation
and growth. In the second case, the maximum horizontal and vertical stresses are 5 MPa and 12
MPa respectively. For both examples, the minimum horizontal stress is 2 MPa.
As seen, when only S3 is available and rock strength is low, the fracture starts to grow once
the fluid pressure exceeds the S1. However, the presence of S2 and S3 changes the starting
pressure. It is important to note that the S3 is no longer 2 MPa when S2 and S1 are present; it is
around 8.5 MPa. When fluid pressure is equal to the minimum horizontal stress, fracture initiation
occurs if the rock strength is zero. However, if the rock strength is not zero and there is not an
open natural fracture network, fluid initiation pressure increases by a factor of 3.5, relative to the
rock strength.
It is interesting to realize that minimum horizontal stress is an important factor as long as
the rock strength is less than a threshold, 1 MPa in this case. Once the rock strength exceeds the
threshold, it becomes the important factor for determination of the fracture initiation onset. We
conclude that the traditional theoretical models need to be modified to incorporate the effect of
rock tensile strength and reservoir stress fields, as both can alter the pressure of fracture initiation.
5.3.2 Sensitivity Analysis on Reservoir Size
As mentioned, simulation of hydraulic fracturing in a real size model requires fine scale reservoir
modeling, and an accurate reservoir characterization model. In this study, we have used a
downsized reservoir and simulated hydraulic fracturing in the model. The computation time
reduces if the dimension of the model decreases, but the reliability of results is reduced. In this
part, the effect of reservoir dimensions on fracture initiation pressure is investigated. Reservoir has
124

three dimensions; length, width, and height. The height in the model is selected based on realistic
values of reservoir thickness; therefore, no sensitivity analysis is needed.
Other factors are the height of overlaying and underlying layers, so a sensitivity analysis
on the thickness of those layers was performed. Width of reservoir was also subjected to analysis.
In addition to width, the distance to the nearest fracture can be an important parameter. However,
a real size reservoir needs to be constructed, which is beyond our scope, but may be the subject of
further studies.
Two tests were conducted to evaluate the effect of reservoir width and the thickness of
underlying and overlaying layers. In the first, reservoir width was changed and the fluid pressure
for fracture initiation was calculated. The results are in Figure 5.3.

Figure 5.3: Effect of reservoir width on fracture initiation pressure

As illustrated, the required fluid pressure to start the fracturing increases with the width of
the reservoir. The trend is ascending, but no specific relationship can be extracted from the figure.
9
10
11
12
13
14
15
16
17
18
0 5 10 15 20 25
FLUID PRESSURE. MPA
RESERVOIR WIDTH, M
125

In this case, however, the effect of increase in reservoir width is much higher at lower reservoir
widths compared to the higher one. In this study, the reservoir width is 7 m, but that does not mean
that the results are less sensitive at that width. The reason for selecting 7 m was simply to achieve
a balance between the model accuracy and computation speed.
In the second test, the effect of the overlaying and underlying layers was examined.

Figure 5.4: Effect of overlaying and underlying sand layers on fracture initiation pressure

As demonstrated by Figure 5.4, the sand layer thickness does not significantly change the
stress field around the fracture. Therefore, the arbitrary selection of sand layer does effect the fluid
pressure. It is notable that in the field, the sand layer thickness changes the stress field around the
fracture -- it has higher bulk density and thus alters the overburden pressure. However, in this test,
we assumed that the pressure is applied on top of the sand layer. According to this scenario, sand
layer thickness would not change the fracture initiation onset. Therefore, the selection of 2m for
sand layer thickness is an acceptable assumption.
10
11
12
13
14
15
16
0 1 2 3 4 5 6
Fluid Pressure, MPa
Sand Layer Thickness, m
126

5.4 Field Scale Numerical Analysis of Hydraulic Fracturing
In the previous sections, I explained how different parameters influence the model response in
hydraulic fracturing. Now, we test the model that was developed based on the formulations given
in Chapter 3. The parameters were selected in such a way as to optimize the simulation time and
accuracy. Description of the model is available in Part 5.2.
The fluid properties are given in Table 5.2.
Table 5.2: Hydraulic fracturing fluid properties
Viscosity (cp) 2
Compressibility (1/psi) 0.00001
Flow rate (m
3
/s) 0.004
Density (Kg/m
3
) 1000
Initial Reservoir Pressure (psi) 3000

In order to start the simulation, an initial fracture is needed to provide a space into which
the fluid is injected. Hence, an initial fracture is placed in the center of reservoir. It is a 3 m by 4
m fracture and is shown in Figure 5.5.

Figure 5.5: Initial fracture is place at the center of reservoir. The center of fracture and reservoir
coincide.
127

The fluid is injected into the two central elements to mimic hydraulic fracturing. Discharge
is divided equally between two elements; we therefore expect balanced propagation of the
fracturing. The model is symmetrical in xy, yz, and xz plane. Although the model size could be
reduced to 1/8 of its original size, the full size simulation of the model was used to avoid the
introduction of any extraneous errors.
The reservoir is fixed at the eight corners. Minimum horizontal stress, maximum horizontal
stress, and vertical stress are assumed 22.41 MPa, 25.41 MPa, and 32.41 MPa respectively.
According to Terzaghi’s effective stress relationship:

𝑃 𝑒𝑓𝑓 = 𝜎 − 𝑃 ,
(5-1)
Where 𝑃 𝑒𝑓𝑓 is effective stress, 𝜎 total stress, and P is reservoir pressure. The minimum, maximum,
and vertical stresses are translated to 2 MPa, 5 MPa, and 12 MPa respectively.
Simulation starts with the injection of fluid into central elements that are already filled with
a small amount of fluid. The fracture propagation is shown Figure 5.6.
Figure 5.6 demonstrates a radial pattern of fracture initiation and propagation in the early
stages. As one can see, the direction of fracture propagation is perpendicular to the minimum
horizontal stress. The fracture pattern is penny-shaped since the fracture has not reached the
boundary. Pattern change is expected once the boundary effect is amplified, which occurs when
the fracture is close to the reservoir boundary or near layer interfaces. Figure 5.6(d) depicts a
transient state at which the crack damage occurred in the elements in the tip. However, the
displacement in the cohesive element is not large enough to cause total damage.

128

(a) (b)

(c) (d)

Figure 5.6: Fracture propagation inside the reservoir after, (a) 9.30 s, (b) 19.48 s, (c) 26.43 s, (d)
27.17 s of fluid injection. Red zone represents the fracture.

As mentioned, the fracture is initially filled with a fluid. The initial length and width of
fracture are 4 m and 3 m respectively. Initially, the fracture propagation occurs in two directions
which is similar to radial propagation of fracture. Therefore, the length and width of fracture is
substituted by the geometric mean of length and width. The result of fracture propagation is in
Figure 5.7.
129

Figure 5.7: Fracture propagation in time. Y-axis shows the equivalent radius of the fracture.

Figure 5.7, shows that the fracture propagates in a discrete form because it opens element
by element. The fracture propagates to elements at the boundary once the damage criteria is met.
Then, the fracture jumps an entire element. Therefore, fracture propagation has a discrete form. In
a high resolution model, discrete jump disappears.
5.5 Model Verification
The approach presented in this research is fundamentally different than the existing methods of
modeling hydraulic fracturing. It relies on a new formulation, and relaxes some assumptions that
cause nonphysical results. PKN, KGD, and penny-shaped models are not appropriate for
comparison, as they have different approaches for modeling of fractures. However, that does not
mean that all of the standard fracture mechanics methods are incapable of generating similar
results. For example, Sneddon (1945) developed a method to estimate fracture width based on the
fluid pressure. It can be used here to verify the accuracy of the model in calculating the fracture
width.
0 5 10 15 20 25 30
1
2
3
4
5
Time, t (s)
Equivalent Radius, Req (m)
130

In this section, we have compared the model to analytical solutions to determine the
accuracy of our model. Fracture width, length, volume, and fluid pressure are the parameters that
were selected for that comparison.
5.5.1 Fracture Width
A reliable model should be able to accurately predict the fracture width and its distribution.
Fracture width is maximum at the injection point, then reduces toward the tip. It can be estimated
using the classic hydraulic fracturing models such as PKN and KGD, or it can be predicted from
the fracture mechanics models. In the classic models, there is a pressure distribution inside the
fracture which is different than the developed model. In such approaches, it is allowed for fluid
inside the fractures to take negative pressure which might not be physical. In our method we allow
the condition of the conservation of mass to be relaxed. Therefore, pressure will not take negative
values.
According to the Sneddon solution:

𝑊 =
8( 𝑝 −𝜎 ) 𝐶 ( 1−𝑣 2
)
𝜋𝐸
√
( 1 − (
𝑟 𝐶 )
2
) ,
(5-2)
Where P is a uniform pressure acting on the fracture walls, σ is minimum horizontal stress, C is
fracture radius, v is Poisson’s ratio, E is Young’s modulus, and r is the distance from the fracture
center.
In the beginning of hydraulic fracturing, an initial pressure of 5 MPa is assigned to the
fluid. It keeps the fracture open and prevents fluctuation in the fluid pressure after the injection
has begun. Therefore, fluid pressure inside the fracture remains uniform before fracture
propagation occurs. We conclude that the relation suggested by Sneddon is applicable for
verification of the hydraulic fracturing model.
131

The result of analysis is as follows:

Figure 5.8: Analytical and numerical calculation of fracture width for a penny-shaped fracture.

Figure 5.8 presents the fracture width change for the central element in time. Both
analytical and numerical methods are used to calculate the width. As demonstrated in the Figure,
at the early stages of the fracture initiation, the analytical and numerical models match perfectly,
but that does not continue throughout the entire simulation. After ~9s, the fracture propagation
starts and there is a sudden decline in the fracture width, because once the fracture starts to
propagate, there is a larger area for the same amount of fluid. Therefore, considering
compressibility of the fluid, fluid pressure suddenly drops, which causes the fracture to become
narrower. All these processes happen in a fraction of a second.
After the initial fracture growth, the gap between numerical and analytical solution grows
due to the boundary condition. In the beginning, the fracture is far from the boundary layers.
132

Fracture propagation causes the amplified boundary effect, which does change the fracture width
over time. The relationship between numerical and analytical solutions is outlined in Figure 5.9.

Figure 5.9: Ratio of numerical to analytical solution in a penny-shaped fracture

One can see the ratio of numerical to analytical values is fairly stable. Jumps occur during
fracture propagation but remain quite constant during the fluid pressurizing phase. Generally, there
is good agreement between numerical simulation and analytical solution. That proves that the
model generates accurate results and can be used for further analysis.
5.5.2 Fracture Volume
Both fracture width and fracture volume are indicators of a successful hydraulic fracturing
process. To verify fracture width, our focus was at the central element. However, fracture
volume is an integral part of fracture width, and is therefore more important in the sense that it
involves all of the fractures. The Griffith-Sneddon equation was used to verify the model for
fracture volume calculation. Given:

𝑊𝑉 =
16( 1−𝑣 2
) 𝐶 3
( 𝑃 −𝜎 )
3𝐸 ,
(5-3)
133

Where V is total volume of fracture.
Comparison of numerical and analytical model is shown in Figure 5.10.

Figure 5.10: Analytical and numerical representation of the fracture volume

As seen, the numerical modeling of the fracture volume is in accord with the analytical
model suggested by Griffith-Sneddon. Both methods match at the early stages of the fracture
propagation, but the numerical model deviates from the analytical solution. The reasons for that
behavior are twofold: (1) the fracture does not have a circular shape, so it does not represent a
perfect penny-shaped fracture; and (2) the boundary condition effect is amplified when the fracture
approaches the boundary of the model. More explanation is available in Part 5.3.2.
5.5.3 Fluid Pressure
As explained in Section 3.5, the algorithm consists of a fluid flow simulator and rock deformation
simulator. In this Part, the fluid flow simulator was tested using the conservation of fluid approach.
Figure 5.11 shows that fracture propagation has a discrete pattern. In other words, fluid is
pressurized and the stress on fracture walls increases until rock damage occurs. In each step,
134

conservation of mass governs the equations. In our model, we tested the conservation of mass in
every step to predict fluid pressure, and then compared the results to the numerical model. The
results are below:

Figure 5.11: Analytical and numerical calculation of fluid pressure

Figure 5.11 shows the fluid pressure over time for the entire fracture. In the model, there
is a small difference between fluid pressure at the central elements and the elements near the
fracture tip. Therefore, an assumption of constant pressure for the entire reservoir does not create
a miscalculation. One notices in the Figure that the gap between numerical and analytical solution
expands every time fracture propagation occurs. It needs to be taken into account that the pressure
at the beginning of each fracture propagation step is selected as an initial pressure for that step.
Hence, the pressure in each step is evaluated independently from all other steps of the propagation.
Two major factors contribute to the fluid pressure -- fracture width and fluid injection
flowrate. Both terms are available in the right side of Equation (3-2). A small change of any factor
generates a huge change in pressure. A quick calculation shows that a 1 percent change in fracture
width, assuming compressibility factor of 10
-5
, causes 1000 psi or 6.8 MPa pressure change.
135

5.5.4 Fracture Length
The developed model is able to simulate different patterns that might occur in hydraulic fracturing.
In the first phase, the dominate mechanism of hydraulic fracturing is based on the radial
propagation of the fracture. However, in the later stages, the fracture is contained between two
limiting zones. In this section, the first part is considered to validate the model -- meaning that we
do not strive to predict the fracture propagation when the fracture is no longer a penny-shape. This
is a valid approach since a fracture loses its shape once it encounters the reservoir boundary. The
calculation of fracture width is presented by Geertsma and de Klark (1969):

𝑤 𝑚𝑎𝑥
= 2√
𝜇𝑄𝑅 𝐺 4
,
(5-4)
Where G is shear modulus, and R is fracture radius. Now, we convert the equation to give us radius:

𝑅 =
𝐺 𝜇𝑄
(
𝑤 𝑚𝑎𝑥 2
)
4
,
(5-5)
Using that equation, we are able to estimate fracture radius (or length) analytically. The result is
shown below:

Figure 5.12: Analytical and numerical calculation of the fracture radius in the early time of the
hydraulic fracturing
136

The classic hydraulic fracturing models and the new model treat rock damage and fracture
propagation criteria quite differently; the disparity between the two is illustrated in Figure 5.12.
There is a three orders of magnitude difference between the analytical and numerical models for
the calculation of fracture length. In the classic model, the fracture length is in the order of
thousands of meters, which may not be true, because in that approach, radius is a strong function
of the width:

𝑅 𝐿 ∝ 𝑊 4
,
(5-6)
Therefore, a small change in fracture width creates a large change in the fracture length (radius).
Fracture width at the outset is effected by the rock tensile strength, which is sometimes ignored in
the classic methods, and can lead to miscalculations. In other words, there is a difference between
analytical and the proposed method, which is related to the formulation and assumptions.
5.6 Discussion
In the previous section, we saw that the model matches the analytical solutions available in the
fracture mechanics. Fracture width and volume, and fluid pressure are calculated to produce a
realistic early and transient time hydraulic fracturing operation. The results are close to the
analytical solution. However, fracture length was not within the range of the analytical solution.
Now we will discuss the discrepancy between the numerical and analytical models.
Pressure distribution inside the reservoir and the effect of leak-off during early stages of
the fracturing will be analyzed in this section. In our model, in contrast to the classic models, we
assumed that there is no filter cake on the fracture wall. It could be a valid assumption at the
beginning of hydraulic fracturing, as there has not been enough time for the layer to form. Our
model can be improved to include the filter cake formation into the calculations.
137

5.6.1 Fracture Length Analysis in Hydraulic Fracturing
Fracture length (radius) was determined using analytical and numerical approaches. There is a
substantial gap between the two methods which needs to be explained. Fracture length and width
are two major parameters for hydraulic fracturing assessment; accurate evaluation of those factors
is vital for production forecast.
The analytical method that was used in Part 5.5.1-3 was based on the static analysis of the
crack. Those techniques were developed for the static stress analysis of cracks. Traditional
hydraulic fracturing simulation methods were based on those formulations, with the addition of
the lubrication theory, which incorporates the fluid pressure distribution inside the fracture. Zero
toughness rock was assumed for the modeling, which means the fracture constantly grows. Two
factors help resist fracture propagation -- minimum horizontal stress, and negative pressure in
vicinity of the crack tip (Morgan, 2014, Salmizadeh and Khalili, 2015). The existence of the lag
between fluid front and fracture tip leads to the creation of a negative pressure area. Therefore, the
assumption of no lag between fracture tip and fluid front is translated to the negative pressure of
fluid in the vicinity of the tip to satisfy the conservation of mass condition. However, we have
lessened such conditions during fracture propagation (and not expansion phase) to avoid any
creation of negative pressure. Meshing is another factor in the existence of a sizeable negative
pressure area in our model. In each step of propagation, the fracture opens and causes a dramatic
decline in fluid pressure. In order to avoid negative pressure, smaller sized elements need to be
generated. The effect of mesh size was studied by Zielonka et al (2014).
5.6.2 Reservoir Pressure and Leak-Off
Different factors precipitate the change in reservoir pressure during hydraulic fracturing. One such
event is fluid leak-off from the fracture to the reservoir which in turn pressurizes the reservoir. On
138

the other hand, the fracture has no pressure when newly formed, which creates suction of the fluid
from reservoir to the fracture. A combination of the above determines the reservoir pressure. The
difference between reservoir pressure and the fluid pressure in said reservoir causes the leak-off.
Therefore, understanding reservoir pressure provides an acceptable indicator for leak-off. For this
purpose, the reservoir pressure change is shown in Figure 5.13.

Figure 5.13: Pressure change in the reservoir elements adjacent to the fracture wall during
fracture propagation. The scale covers -400 psi to 100 psi change of the fluid pressure compared to initial
3000 psi.

Figure 5.13 shows elements inside the reservoir that are in vicinity of fracture. These
elements have the most interaction with the fracture and pressure change inside the fracture directly
affects their pressure. One of the interesting results of such interaction is that the pressure within
139

the reservoir may drop during the fracture propagation process. Fluid migration from a fracture to
the reservoir at the central elements increases the pressure in those regions. Suction in the vicinity
of the fracture tip prompted the pressure drop in the reservoir elements. In contrast to the fracture
elements, elevated contrast between reservoir elements does not generate enough pressure gradient
for the fluid to move from the high pressure region to a low pressure region. The reason for that is
the reservoir permeability is so low that fluid migration in the reservoir is almost zero. However,
the fracture has high permeability, which allows fluid to move from the fracture to the formation
and vice versa.
5.7 Conclusion Remarks
In this chapter, a model was constructed to test the developed algorithm for simulation of hydraulic
fracturing. A sensitivity analysis was performed on the model to identify the important parameters
and their effect on the simulation. Reservoir size, and rock damage properties were analyzed. It
was demonstrated in three dimensional modeling, the model width plays an important role in the
hydraulic fracturing simulation. A wider reservoir model shows higher strength toward the fluid
injection pressure. Therefore, three dimensional model should be large enough to reduce the
sensitivity of the response to the model width. To diminish the sensitivity of reservoir on size,
utilization of infinite elements is recommended.
In addition to sensitivity analysis, this chapter aimed to demonstrate the capabilities of the
model to identify the mechanisms that may not be perceived in the other approaches. In first step,
our model was verified against the analytical solutions. It was shown that the numerical evaluation
of fracture width and volume, and fluid pressure are slightly different than the analytical solution.
However, the fracture length does not match the analytical solution. That discrepancy can be
140

attributed to the fundamental differences between the analytical and numerical models. It is
recommended that high resolution meshing be used to lessen the difference.
Fluid leak-off and fluid suction near a fracture tip were identified as two factors that alter
reservoir pressure. It was demonstrated that the effect of those parameters is greater than
previously thought. They can alter the reservoir pressure near hydraulic fracture.
In summary, it is necessary to consider many parameters in hydraulic fracturing simulation
that previously were neglected in order to simplify the analysis. This chapter introduced some of
those parameters and explained their importance in the simulation.
141

Chapter 6 : SUMMARY AND FUTURE RESEARCH DIRECTION
6.1 Summary
This research introduces a new algorithm for simulation of hydraulic fracturing. It integrates the
capabilities of Abaqus, a powerful FEM tool, and a developed in-house FDM code to predict
hydraulic fracturing in a field scale model. In contrast to previously developed models, this model
can transcend the single fracture propagation assumption, and simulates stress shadowing and
hydraulic and natural fractures interaction. Examples of such models was introduced in Chapter 4.
Traditional models were developed to simulate hydraulic fracturing propagation by
considering only a single fracture. However, as stated in Chapter 2, various factors need to be
incorporated to ensure an accurate simulation. Stress-strain relationships in rock depends upon the
various components which influence rock behavior under stress. In addition, the standard LEFM
may lead to incorrect modeling of fractures since it cannot accurately characterize fracture
response. CZM was suggested to address some of the issues in LEFM, and was discussed in detail
in Section 2.4.
To resolve the issues of theoretical models, numerical formulations were suggested by
researchers to address several hydraulic fracturing problems. FDM, FVM, FEM, and XFEM as
well as hybrid methods have been widely used for solutions. Among the numerical approaches,
XFEM is able to explicitly incorporate fractures and therefore reduce the computation time of the
simulation.
In Chapter 3, an algorithm was presented that takes advantage of XFEM for modeling
fracturing. The algorithm utilized Abaqus and incorporated a realistic rock damage model to
simulate rock response under load. The fluid flow modeling was performed using an in-house
142

FDM code. It integrated fluid flow inside the fracture and fluid flow inside the reservoir, as well
as fluid transfer at the fracture-rock interface. This model provides a reliable tool with which to
monitor pressure changes in a reservoir, which has been too often neglected in hydraulic fracturing
simulations.
An iterative approach was utilized to couple fluid flow and rock deformation. Strong
nonlinear relation between fluid flow and rock deformation dictates the development of a new
scheme which facilitates a convergence towards finding the solution. It is based on implicit
coupling of fluid flow and rock deformation.
Before testing the algorithm, Abaqus was tested to evaluate the capabilities of the software
program to model cracks. Any limitations in the FEM tool is reflected in the algorithm. Therefore,
understanding of such limitations is necessary for the identification of algorithm capabilities.
Fractures propagate in hydraulic fracturing because of the tensile stress in the vicinity of a
fracture. Tensile stress is generated by the presence of high pressure fluid inside the fracture.
Similar conditions were assumed for the analysis of crack propagation in a 2D plate. We noted
that the size of the model and the boundary condition influence fracture propagation. The stress
field applied on the boundary of a fracture introduces extra stress on the fracture wall that should
be overcome to allow propagation. Plane strain and plane stress conditions should be determined
before analysis, as they can alter the results of simulation. It was demonstrated that the distance
between two fractures plays an important role in hydraulic fracturing initiation and propagation.
Distance to the boundary is another important factor that should be taken into consideration when
a finite model is simulated. Mesh resolution was analyzed and it was shown that low resolution
meshing can generate an erroneous result. The interaction of two fractures, natural and hydraulic,
could be modeled using Abaqus. A code was developed for this purpose to enhance the capabilities
143

of the software program. It enabled us to study the behavior of interacting fractures and to perform
a sensitivity analysis on the distance between those fractures.
In Chapter 5, the algorithm was tested on a field scale reservoir model. A model was
constructed and sensitivity analysis was performed on reservoir geometric parameters, as well as
rock model. It demonstrated that change in a rock damage model disturbs the formation response.
Rock with higher strength requires higher pressure to break; therefore, the subsequent fracture
properties are influenced by such pressure. This means that the typical models are not reliable if
they are developed based on no tensile strength. The model was verified against the analytical
solutions of fracture mechanics presented by Griffith-Sneddon. Our model was able to accurately
generate fracture width and volume as well as fluid pressure. However, there was a discrepancy in
fracture length between the analytical and numerical models when the theoretical modeling of
hydraulic fracturing was used. That disparity was attributed to the fundamental differences
between the archetype model and the new algorithm in terms of handling boundary conditions.
The fluid pressure in the reservoir changes over time and is not constant, in contrast to the
theoretical models. We tested the model to study reservoir behavior. The fracture interacts with
matrix both in the central elements and tip. In the newly formed elements near the tip, the presence
of zero pressure zone reduces pressure of fluid in the matrix. In central elements, reservoir pressure
increases due to fluid leak-off to the formation. Such interactive mechanisms have not been
investigated in any of the previously developed models.
6.2 Recommendation and Future Work
Accurate simulation of hydraulic fracturing requires a comprehensive characterization of reservoir,
a valid theoretical model in fracture mechanics, and a rigorous numerical model for simulation.
Characterization of reservoir usually does not generate a detailed model of fracture location, and
144

reservoir and rock properties. Therefore, for every hydraulic fracturing simulation, such
uncertainty in the initial model needs to be taken into account. Apart from the characterization, a
valid theoretical model for fracture mechanics should be developed to model natural and hydraulic
fracturing interaction accurately. In this study, we presented several examples of interaction
between fractures, which provide a useful tool in understanding the mechanism. However, both
laboratory and field studies need to be conducted to validate numerical models.
In addition to the need for increased accuracy of input data and theoretical models, a
numerical model requires modification to increase its capabilities for hydraulic fracturing
simulation. XFEM is an innovative approach for modeling of fractures, and offers a reliable
solution for fracture problems. However, implementing XFEM in Abaqus requires modification.
Recent development in parallel computing urges developers to adapt such techniques in Abaqus
for XFEM modeling of fractures. Moreover, XFEM in Abaqus is limited to the cubical elements
in 3D and rectangular elements in 2D. Removal of this limitation is necessary as well. In our
simulation, we modified code to model constant fluid pressure inside the propagating fracture.
Such capability needs to be incorporated internally into Abaqus.
Our developed algorithm has the foundation to be used in the field scale for hydraulic
fracturing simulation. Following are a few modifications that we propose before embarking on
practical applications in a field scale:
1- To avoid negative pressure in the fractures, the conservation of mass is relaxed during
fracture propagation phase. However, it needs to be implemented in the model to offer an
accurate result. In such conditions, using high resolution meshing is necessary to avoid any
unphysical result.
145

2- Hydraulic fracturing fluid shows non Newtonian behavior. The fluid flow simulator needs
to take into account this behavior.
3- Fluid injection is very high in hydraulic fracturing and it likely exceeds the laminar flow
threshold, so that turbulence fluid flow develops. The model should be modified to adapt
to such behavior.
4- 3D simulation of hydraulic fracturing in a fractured reservoir could offer valuable
information on the behavior of hydraulic fracturing in unconventional reservoirs. The
modeling can start with a single fracture and become more complex by adding more
fractures with different orientations.
5- Construction of a field scale model and simulation of hydraulic fracturing could provide
an invaluable insight into different large-scale phenomena such as stress shadowing,
induced seismicity, and possible water contamination. However, implementation and
simulation in such models demands parallel computing, as well as high performance
computers. Large scale models can be run once the appropriate foundation is provided.

146

Appendix A Derivation of the Fluid Flow Equation
This part is dedicated to the derivation of the fluid flow that is used for hydraulic fracturing
simulation. We start with derivation of conservation of mass and then, discretization of
equations are presented.
𝑚 ̇ 𝑖𝑛
− 𝑚 ̇ 𝑜𝑢𝑡 =
∆𝑚 ∆𝑡
( 𝜌𝑢𝐴 )
𝑖𝑛
− ( ρuA)
𝑜𝑢𝑡 =
𝑑 ( ρVϕ)
𝑑𝑡
In matrix we have:
( 𝜌 𝑢 𝑥 𝐴 𝑥 )
𝑥 − ( ρu
x
A
x
)
𝑥 +∆𝑥 + ( 𝜌 𝑢 𝑦 𝐴 𝑦 )
𝑦 − ( ρu
y
A
y
)
𝑦 +∆𝑦 + ( 𝜌 𝑢 𝑧 𝐴 𝑧 )
𝑧 − ( ρu
z
A
z
)
𝑧 +∆𝑧
=
𝑑 ( ρdxdydz ϕ)
𝑑𝑡
𝐴 𝑥 = 𝑑𝑦𝑑𝑧
𝐴 𝑦 = 𝑑𝑥𝑑𝑧
𝐴 𝑧 = 𝑑𝑥𝑑𝑦
( 𝜌 𝑢 𝑥 𝑑𝑦𝑑𝑧 )
𝑥 − ( ρu
x
𝑑𝑦𝑑𝑧 )
𝑥 +∆𝑥 + ( 𝜌 𝑢 𝑦 𝑑𝑥𝑑𝑧 )
𝑦 − ( ρu
y
𝑑𝑥𝑑𝑧 )
𝑦 +∆𝑦 + ( 𝜌 𝑢 𝑧 𝑑𝑥𝑑𝑦 )
𝑧 − ( ρu
z
𝑑𝑥𝑑𝑦 )
𝑧 +∆𝑧 =
𝑑 ( ρdxdydz ϕ)
𝑑𝑡
−
𝜕 ρu
x
𝜕𝑥
−
𝜕 ρu
y
𝜕𝑦
−
𝜕 ρu
z
𝜕𝑧
=
𝑑 ( ρϕ)
𝑑𝑡
𝑢 𝑥 = −
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥

147

𝑢 𝑦 = −
𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦

𝑢 𝑧 = −
𝑘 𝑧 𝜇 𝜕 𝑃 𝑧 𝜕𝑧

−
𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
)−
𝜕 𝜕𝑦
( ρ
𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦
)−
𝜕 𝜕𝑧
( ρ
𝑘 𝑧 𝜇 𝜕𝑃
𝜕𝑧
) =
𝑑 ( ρϕ)
𝑑𝑡
Assumptions:
𝜇 = Constant, 𝑘 𝑥 = 𝑘 𝑦 = 𝑘 𝑧 = 𝐾 ,
We have:
𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
) = ρ
𝑘 𝑥 𝜇 𝜕 𝜕𝑥
(
𝜕𝑃
𝜕𝑥
) +
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥

𝜕 𝜕𝑥
( ρ)
𝑐 𝑓 =
1
𝜌 𝑑𝜌 𝑑𝑝
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥

𝜕 𝜕𝑥
( ρ) = 𝜌 𝑐 𝑓 𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
(
𝜕 P
𝜕𝑥
)
Assuming
𝜕 𝜕𝑥
(
𝜕𝑃
𝜕𝑥
) ≫ 𝑐 𝑓 (
𝜕 P
𝜕𝑥
)
2
:
𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
) = ρ
𝑘 𝑥 𝜇 𝜕 2
𝑃 𝜕 𝑥 2

We can perform similar operation to y and z directions.
So:
𝑑 ( ρϕ)
𝑑𝑡 = ρ
𝑑 ( ϕ)
𝑑𝑡 + ϕ
𝑑 ( ρ)
𝑑𝑡
148

𝑑 ( ϕ)
𝑑𝑡 = 𝑐 𝑟 𝜙 𝑑𝑝 𝑑𝑡
𝑑 ( ρ)
𝑑𝑡 = 𝑐 𝑓 𝜌 𝑑𝑝 𝑑𝑡
Substituting these to EQ. we have:
𝑑 ( ρϕ)
𝑑𝑡 = ρ𝑐 𝑟 𝜙 𝑑𝑝 𝑑𝑡 + 𝜌 𝑐 𝑓 ϕ
𝑑𝑝 𝑑𝑡
= ρ( 𝑐 𝑟 + 𝑐 𝑓 ) 𝜙 𝑑𝑝 𝑑𝑡
Total compressibility is:
𝑐 𝑡 = 𝑐 𝑟 + 𝑐 𝑓
ρ
𝐾 𝜇 𝜕 2
𝑃 𝜕 𝑥 2
+ ρ
𝐾 𝜇 𝜕 2
𝑃 𝜕 𝑦 2
+ ρ
𝐾 𝜇 𝜕 2
𝑃 𝜕 𝑧 2
= ρc
t
𝜙 𝑑𝑝 𝑑𝑡
𝜕 2
𝑃 𝜕 𝑥 2
+
𝜕 2
𝑃 𝜕 𝑦 2
+
𝜕 2
𝑃 𝜕 𝑧 2
=
𝐾 c
t
𝜙 𝜇 𝑑𝑝 𝑑𝑡
For fracture, the equations are more complicated. We need to address effect of anisotropy which
is encountered in the calculations.
In fracture, we may recall the general equation:
( 𝜌𝑢𝐴 )
𝑖𝑛
− ( ρuA)
𝑜𝑢𝑡 =
𝑑 ( ρVϕ)
𝑑𝑡
Fracture is not filled with any rock structure, therefore we have:
ϕ
f
= 1
Fracture volume can be
149

V = wdydz
Where w is fracture width. And
w = w( x, t)
If we assume w = dx,
V = dxdydz
Equation cane be written as follows:
( 𝜌 𝑢 𝑥 𝐴 𝑥 )
𝑥 − ( ρu
x
A
x
)
𝑥 +∆𝑥 + ( 𝜌 𝑢 𝑦 𝐴 𝑦 )
𝑦 − ( ρu
y
A
y
)
𝑦 +∆𝑦 + ( 𝜌 𝑢 𝑧 𝐴 𝑧 )
𝑧 − ( ρu
z
A
z
)
𝑧 +∆𝑧
=
𝑑 ( ρdxdydz )
𝑑𝑡
𝑑 ( ρdx)
𝑑𝑡 = dydz
𝑑 ( ρdx)
𝑑𝑡 = dydz
𝑑 ( ρw)
𝑑𝑡
𝑑 ( ρw)
𝑑𝑡 = ρ
𝑑 ( w)
𝑑𝑡 + w
𝑑 ( ρ)
𝑑𝑡
So we have:
( 𝜌 𝑢 𝑥 𝐴 𝑥 )
𝑥 − ( ρu
x
A
x
)
𝑥 +∆𝑥 + ( 𝜌 𝑢 𝑦 𝐴 𝑦 )
𝑦 − ( ρu
y
A
y
)
𝑦 +∆𝑦 + ( 𝜌 𝑢 𝑧 𝐴 𝑧 )
𝑧 − ( ρu
z
A
z
)
𝑧 +∆𝑧
= dydzdx
𝑑 ( ρ)
𝑑𝑡 + ρdydz
𝑑 ( w)
𝑑𝑡
Assuming fracture is small compare to reservoir, and
𝑢 𝑥 = −
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥

𝑢 𝑦 = −
𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦

150

𝑢 𝑧 = −
𝑘 𝑧 𝜇 𝜕𝑃
𝜕𝑧

1
dydzdx
( ( 𝜌 𝑢 𝑥 𝐴 𝑥 )
𝑥 − ( ρu
x
A
x
)
𝑥 +∆𝑥 + ( 𝜌 𝑢 𝑦 𝐴 𝑦 )
𝑦 − ( ρu
y
A
y
)
𝑦 +∆𝑦 + ( 𝜌 𝑢 𝑧 𝐴 𝑧 )
𝑧 − ( ρu
z
A
z
)
𝑧 +∆𝑧 )
=
𝑑 ( ρ)
𝑑𝑡 + ρ
1
𝑑𝑥 𝑑 ( w)
𝑑𝑡

𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
) +
1
𝑤 𝜕 𝜕𝑦
( ρ𝑤 𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦
) +
1
𝑤 𝜕 𝜕𝑧
( ρ𝑤 𝑘 𝑧 𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌 𝑑 ( P)
𝑑𝑡 + ρ
1
w
𝑑 ( w)
𝑑𝑡
We multiply both sides by w:
𝑤 𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
( ρ𝑤 𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
( ρ𝑤 𝑘 𝑧 𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌𝑤
𝑑 ( P)
𝑑𝑡 + ρ
𝑑 ( w)
𝑑𝑡
We have:
𝑘 𝑥 = 𝐾 ,
According to lubrication theory:
𝐾 𝑓 = 𝑘 𝑦 = 𝑘 𝑧 =
𝑤 2
12

𝑤 𝜕 𝜕𝑥
( ρ
𝑘 𝑥 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
( ρ𝑤 𝑘 𝑦 𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
( ρ𝑤 𝑘 𝑧 𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌𝑤
𝑑 ( P)
𝑑𝑡 + ρ
𝑑 ( w)
𝑑𝑡
𝑤 𝜕 𝜕𝑥
( ρ
𝐾 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
( ρ𝑤 𝑤 2
12𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
( ρ𝑤 𝑤 2
12𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌𝑤
𝑑 ( P)
𝑑𝑡 + ρ
𝑑 ( w)
𝑑𝑡
𝑤 𝜕 𝜕𝑥
( ρ
𝐾 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
( ρ
𝑤 3
12𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
( ρ
𝑤 3
12𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌𝑤
𝑑 ( P)
𝑑𝑡 + ρ
𝑑 ( w)
𝑑𝑡
151

Equation can also be written in form of:
𝑤 𝜕 𝜕𝑥
( ρ
𝐾 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
( ρ𝑤 𝐾 𝑓 𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
( ρ𝑤 𝐾 𝑓 𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝜌𝑤
𝑑 ( P)
𝑑𝑡 + ρ
𝑑 ( w)
𝑑𝑡
Assuming
𝜕 𝜕𝑥
(
𝜕 𝑃 𝑖 𝜕𝑥
) ≫ 𝑐 𝑓 (
𝜕 P
i
𝜕𝑥
)
2
and i = x,y,z:
𝑤 𝜕 𝜕 𝑥 (
𝐾 𝜇 𝜕𝑃
𝜕𝑥
) +
𝜕 𝜕𝑦
(𝑤 𝐾 𝑓 𝜇 𝜕𝑃
𝜕𝑦
) +
𝜕 𝜕𝑧
(𝑤 𝐾 𝑓 𝜇 𝜕𝑃
𝜕𝑧
) = 𝑐 𝑓 𝑤 𝑑 ( P)
𝑑𝑡 +
𝑑 ( w)
𝑑𝑡
The equation describes the fluid flow inside the fracture. The first term in the equation is very
small and can be assumed to zero. However, it does not result in negligible flow perpendicular to
the fracture and into the reservoir. Fluid flow into matrix is modeled through an integrated
approach of fluid flow inside the fracture and from fracture to the matrix.
152

Appendix B Multiple Fracture Propagation under Tensile Load
We will focus on the effect of tensile load on the fracture initiation and propagation. Two fractures
are available in the model, with different configurations. Two general kinds of tests were
conducted to model the fracture behavior. In the first, we assume that two the fractures are parallel,
and there is a tensile force perpendicular to the fractures. In the second case, we assume the
fractures are perpendicular to each other; that discussion is available in the next section.
The size of the model was included in our considerations, and a sensitivity analysis is
evaluated on the distance between the two fractures. We begin with a rectangular plate.
For the analysis of parallel fracturing, a rectangular model is designed. The model
dimension is similar to the one shown in Figure 4.1. The schematic is shown in Figure B.1.

B.1: Plate is subjected to the two tensile loads acting from the top and bottom of the plate. Two fractures
are placed at the plate, and are separated by a distance of h.

12.5 cm
100 cm
50 cm
h
153

As shown in Figure B.1, two loads are acting in opposite directions. The plate is not
supported with any external force. Therefore, every part of the plate has the freedom to move or
rotate in a different direction.
In Figure B.1, we applied the load gradually, and allowed the fractures to propagate. For
this scenario, we assumed the initial distance between two fractures was 0.1 m. The result of the
test is shown in Figure B.2.
One can see from Figure B.2 that the fractures starts to grow when the load reaches 1.091
MPa. In that event, the fractures start to divert away from each other. If, for example, we refer to
single fracture propagation under tensile load, we observe similar behavior of the fracture. Similar
patterns could be seen here, considering that the fractures are under tensile load from both sides.
Figure B.3 shows the stress in y direction for the load equal to 1.091 MPa.
We notice in Figure B.3 the tensile force acting in the process zone; it is concentered in the
fracture tip. It is the dominant stress around the fracture, and far from the tip. The stress field
around the fracture tip is such that it pushes the fractures away from each other. However, this is
only effective when the fractures are too close together.

154

(a) (b) (c)

(d) (e)
B.2: Stress field distribution in a rectangular plate. Distance between two fractures is 10 cm. Fluid
pressure is (a) 1.091×10
6
Pa, (b) 1.83×10
6
Pa, (c) 2.197×10
6
Pa, (d) 2.505×10
6
Pa, and (e) 2.743×10
6
Pa.
The width of fractures are magnified by factor 10.

155

B.3: Normal stress distribution y direction.

For the next test, the distance between fractures increases, hence, the interaction between
them decreases. The distance between two fractures is 30 cm in this test. The result is shown in
Figure B.4.
Figure B.4(a) shows how cracks influence the stress field, and change the stress distribution
around other fractures. It also represents the moment that the fracture starts to grow. It is of interest
that in Figure B.4(a)-(c), there is a strong tensile stress between two fracture tips. This stress
diminishes gradually as the tips progress in the plate. The reason for that behavior is that the
boundary effect is augmented when the tip is close to the right edge. That changes the stress
distribution, and limits the stress increase in front of the crack tip.

156

(a) (b) (c)

(d) (e)
B.4: Stress field distribution in a rectangular plate. Distance between two fractures is 30 cm. Fluid
pressure is (a) 1.105×10
6
Pa, (b) 1.954×10
6
Pa, (c) 2.268×10
6
Pa, (d) 2.671×10
6
Pa, and (e) 2. 694×10
6
Pa.
The width of fractures are magnified by factor 30.

One can compare Figure B.4(d) to e, and conclude the effects of the boundary. It is
noteworthy that the difference between the two loads is only 23 KPa, which is less than one percent
of the loads, but there is an order of magnitude difference between Von Mises stresses. It can be
explained by looking at the size of undamaged elements in front of the fracture tip. In Figure
157

B.4(e), a small element is responsible to transfer the stress across the fracture. Therefore, there is
a stress concentration on the element which eventually causes the element failure.
The distance between two fractures is increased to 50 cm. The two fracture are located 25
cm from the edges. The result of the test is illustrated in Figure B.5.

(a) (b) (c)

(d) (e)
B.5: Stress field distribution in a rectangular plate. Distance between two fractures is 50 cm. Fluid
pressure is (a) 0.934×10
6
Pa, (b) 1.56×10
6
Pa, (c) 2.268×10
6
Pa, (d) 2.697×10
6
Pa, and (e) 2. 729×10
6
Pa.
The width of fractures are magnified by factor 50.
158

Figure B.5 illustrates how the increase in distance between the two fractures can change
the stress distribution, as well as failure initiation and total damage. It can be readily seen that
Figure B.4(a) and Figure B.5(a) have different types of stress fields. In the former, the stress field
at the crack tip of both fractures is tied together. However, in the latter, there is no obvious
interaction between the two stress fields. It is also worth mentioning how the stress fields between
the fractures in Figure B.5 remain fairly similar compared to Figure B.4.
To further explain the importance of the fracture distance, Figure B.6 is plotted to show
the initiation of the crack, and total damage of the material in contrast to the distance between
fractures.

B.6: Fracture initiation and plate failure stress for different tensile loads in different fracture positions

One limitation of XFEM in Abaqus is that damaged elements should be at least 3 elements
apart. Therefore, in order to properly understand the effect of fracture distance within the
0
0.5
1
1.5
2
2.5
3
0 10 20 30 40 50 60
Load, MPa
Distance between fractures, cm
Initiation
TotalDamage
159

propagation pattern, a higher resolution meshing scheme should be utilized. Nonetheless, Figure
B.6 presents an invaluable result regarding the effect of distance. It shows how the trend of fracture
initiation and breaking point is influenced by the distance between fractures. Figure 4.18 and
Figure 4.19 show that the closer a single fracture is to the edge, a lower load is needed to damage
the element. However, it was not true in the case of two fractures. One expected a higher load for
the fracture initiation when the distance is minimum, which was not the case. The interaction of
two fractures cause a lower pressure to be needed for fracture propagation. That can happen
because the second fracture acts as the boundary, and because of the short distance between two
fractures, so the damage occurs from a lower load.
Next, I studied the effect of fracture orientation on the propagation. The test is dedicated
to investigate the effect of tensile load on the fracture behavior.
At this juncture I will discuss the effect of fracture orientation on the propagation pattern
when two orthogonal fractures are present in a model. Fractures were placed in two sides of the
model and tensile load was applied on the top and bottom edges. The schematic of the model is
shown in Figure B.7.
Application of load causes the horizontal fracture to grow. Figure B.8 illustrates the
pattern of fracture propagation with a load increase.

160

B.7: Plate is subjected to the two tensile loads acting from the top and bottom of the plate. The initial
fracture is placed at the center left part of the plate and the second one in the right side, next to the
boundary.

Figure B.8 illustrates an essential aspect of the analysis by showing how the element, under
maximum load, changes its location from around the crack tip to the natural fracture tip. It therefore
triggers fracture propagation from the natural fractures, which are not connected to the horizontal
fractures. This examples demonstrates how natural fractures can be activated, although it does not
represent the boundary condition in the reservoir. These results cannot be extended to the field
case, because pressure is in the propagating fracture, and there is no tensile load acting in the
reservoir. Nevertheless, the test can provide us with clues to how fractures may behave under
different circumstances.

12.5 cm
100 cm
50 cm
10 cm
161

(a) (b) (c)

(d) (e) (f)
B.8: Stress field distribution in a rectangular plate. Distance between two fractures is 33.75 cm. Fluid
pressure is (a) 1.067×10
6
Pa, (b) 2.389×10
6
Pa, (c) 2.550×10
6
Pa, (d) 2.624×10
6
Pa, (e) 2.680×10
6
Pa, and
(f) 3.123×10
6
Pa. The width of fractures are magnified by factor 10.

162

References
Abe, H., T. Mura, and L. M. Keer. "Growth rate of a penny‐shaped crack in hydraulic fracturing
of rocks." Journal of Geophysical Research 81.29 (1976): 5335-5340.
Adachi, J., et al. "Computer simulation of hydraulic fractures." International Journal of Rock
Mechanics and Mining Sciences 44.5 (2007): 739-757.
Adachi, José I., and Emmanuel Detournay. "Plane strain propagation of a hydraulic fracture in a
permeable rock." Engineering Fracture Mechanics 75.16 (2008): 4666-4694.
Al-Shayea, Naser A. "Crack propagation trajectories for rocks under mixed mode I–II
fracture." Engineering Geology 81.1 (2005): 84-97.
Anderson, Ted L. Fracture mechanics: fundamentals and applications. CRC press, 2005.
Atkinson, B. K., P. G. Meredith, The theory of subcritical crack growth with applications to
minerals and rocks, Fracture Mechanics of Rock B. K. Atkinson, 111–166, Academic, San
Diego, Calif., 1987a.
Azim, Reda Abdel, and Sherif S. Abdelmoneim. "Modeling hydraulic fractures in finite difference
simulators using amalgam local grid refinement (LGR)." Journal of Petroleum Exploration
and Production Technology 3.1 (2013): 21-35.
Bachman, R. C., et al. "Coupled simulation of reservoir flow, geomechanics, and formation
plugging with application to high-rate produced water reinjection." SPE Reservoir
Simulation Symposium. 2003.
Backers, Tobias. Fracture toughness determination and micromechanics of rock under mode I and
mode II loading. Geoforschungszentrum, 2005.
Barenblatt, G_I_. "The formation of equilibrium cracks during brittle fracture. General ideas and
hypotheses. Axially-symmetric cracks." Journal of Applied Mathematics and Mechanics
23.3 (1959): 622-636.
Bazant, Zdenek P., and Jaime Planas. Fracture and size effect in concrete and other quasibrittle
materials. Vol. 16. CRC press, 1997.
Benzeggagh, M. L., and M. Kenane. "Measurement of mixed-mode delamination fracture
toughness of unidirectional glass/epoxy composites with mixed-mode bending
apparatus." Composites science and technology 56.4 (1996): 439-449.
Bao, J. Q., E. Fathi, and S. Ameri. "Uniform investigation of hydraulic fracturing propagation
regimes in the plane strain model." International Journal for Numerical and Analytical
Methods in Geomechanics 39.5 (2015): 507-523.
163

Blanton, Thomas L. "An experimental study of interaction between hydraulically induced and pre-
existing fractures." SPE Unconventional Gas Recovery Symposium. Society of Petroleum
Engineers, 1982.
Boone, Thomas J., Paul A. Wawrzynek, and Anthony R. Ingraffea. "Finite element modelling of
fracture propagation in orthotropic materials." Engineering fracture mechanics 26.2
(1987): 185-201.
Boone, Thomas James, Ingraffea, A. R and Cornell University. School of Civil and Environmental
Engineering Simulation and visualization of hydraulic fracture propagation in poroelastic
rock. School of Civil and Environmental Engineering, Cornell University, Ithaca, N.Y,
1989.
Buehler, Markus J., and Theodor Ackbarow. "Fracture mechanics of protein materials." Materials
Today 10.9 (2007): 46-58.
Busetti, Seth. Fracturing of layered reservoir rocks. Vol. 71. No. 03. 2009.
Busetti, Seth, Kyran Mish, and Ze'ev Reches. "Damage and plastic deformation of reservoir rocks:
Part 1. Damage fracturing." AAPG bulletin 96.9 (2012): 1687-1709.
Chen, Zuorong. "An ABAQUS implementation of the XFEM for hydraulic fracture
problems." ISRM International Conference for Effective and Sustainable Hydraulic
Fracturing. International Society for Rock Mechanics, 2013.
Cipolla, Craig, et al. "The relationship between fracture complexity, reservoir properties, and
fracture treatment design." SPE Annual Technical Conference and Exhibition. 2008.
Cleary, M. P. "Rate and structure sensitivity in hydraulic fracturing of fluid-saturated porous
formations." 20th US Symposium on Rock Mechanics (USRMS). 1979.
Cleary, Michael P., Michael Kavvadas, and Khin Yong Lam. “Development of a fully three-
dimensional simulator for analysis and design of hydraulic fracturing.” No. CONF-
830305-. Massachusetts Inst. of Technology, 1983.
Cox, S. J. D., and P. G. Meredith. "Microcrack formation and material softening in rock measured
by monitoring acoustic emissions." International journal of rock mechanics and mining
sciences & geomechanics abstracts. Vol. 30. No. 1. Pergamon, 1993.
Crook, Anthony JL, Jian-Guo Yu, and Stephen M. Willson. "Development of an orthotropic 3D
elastoplastic material model for shale." SPE/ISRM Rock Mechanics Conference. Society
of Petroleum Engineers, 2002.
Cundall, Peter A., and Roger D. Hart. "Numerical modelling of discontinua." Engineering
computations 9.2 (1992): 101-113.
Dalguer, L. A., K. Irikura, and J. D. Riera. "Simulation of tensile crack generation by three‐
dimensional dynamic shear rupture propagation during an earthquake." Journal of
Geophysical Research: Solid Earth (1978 –2012) 108.B3 (2003).
164

Daneshy, Abbas. "On the design of vertical hydraulic fractures." Journal of Petroleum
Technology 25.1 (1973): 83-97.
Daneshy, A. A. "Hydraulic fracture propagation in layered formations." Society of Petroleum
Engineers Journal 18.01 (1978): 33-41.
De Pater, C. J., and L. J. L. Beugelsdijk. "Experiments and numerical simulation of hydraulic
fracturing in naturally fractured rock." Alaska Rocks 2005, the 40th US Symposium on Rock
Mechanics (USRMS). 2005.
Deng, S., Podgorney, R., and Huang, H.: Discrete element modeling of rock deformation, fracture
network development and permeability evolution under hydraulic stimulation. Paper
presented at the thirty-sixth workshop on geothermal reservoir engineering, Stanford.
University, http://www.geothermal-energy.org/pdf/IGAstandard/SGW/2011/deng.pdf?
(2011)
DiGiovanni, A. A., et al. "Microscale damage evolution in compacting sandstone." Geological
Society, London, Special Publications 289.1 (2007): 89-103.
Doan, Mai-Linh, and Gérard Gary. "Rock pulverization at high strain rate near the San Andreas
fault." Nature geoscience 2.10 (2009): 709-712.
Dolbow, John, Nicolas Moës, and Ted Belytschko. "An extended finite element method for
modeling crack growth with frictional contact." Computer Methods in Applied Mechanics
and Engineering 190.51 (2001): 6825-6846.
Dow, John O., Michael S. Jones, and Shawn A. Harwood. "A generalized finite difference method
for solid mechanics." Numerical Methods for Partial Differential Equations 6.2 (1990):
137-152
Dugdale, D. S. "Yielding of steel sheets containing slits." Journal of the Mechanics and Physics
of Solids 8.2 (1960): 100-104.
Stimulation, Reservoir. "MJ Economides and KG Nolte." (2000): 10-47.
Economides, M. J. "Fluid-leakoff delineation in high-permeability fracturing." SPE production &
facilities 14.02 (1999): 110-116.
Eichhubl, Peter. "Growth of ductile opening-mode fractures in geomaterials."Geological Society,
London, Special Publications 231.1 (2004): 11-24.
Eichhubl, Peter, and Atilla Aydin. "Ductile opening-mode fracture by pore growth and coalescence
during combustion alteration of siliceous mudstone." Journal of Structural Geology 25.1
(2003): 121-134.
Ertekin, T., Abou-Kassem, J. H. and King, G. R. 2001. Basic Applied Reservoir Simulation, 406
Richardson, TX: Society of Petroleum Engineers
Faruque, M.D. "Development of a generalized constitutive model and its implementation in soil-
structure interaction (plasticity)." (1983).
165

Fernandez, F. A., and Lukasz Kulas. "A simple finite difference approach using unstructured
meshes from FEM mesh generators." Microwaves, Radar and Wireless Communications,
2004. MIKON-2004. 15th International Conference on. Vol. 2. IEEE, 2004.
Foss, David R., and John F. Brady. "Structure, diffusion and rheology of Brownian suspensions
by Stokesian dynamics simulation." Journal of Fluid Mechanics 407.1 (2000): 167-200.
Gale, J. E., et al. "Flow in rocks with deformable fractures." Finite Element Methods in Flow
Problems (1974): 583-598.
Garagash, Dmitry I. "Propagation of a plane-strain hydraulic fracture with a fluid lag: Early-time
solution." International journal of solids and structures 43.18 (2006): 5811-5835.
Geertsma, J., and F. De Klerk. "A rapid method of predicting width and extent of hydraulically
induced fractures." Journal of Petroleum Technology 21.12 (1969): 1571-1581.
Golshani, Aliakbar, et al. "Numerical simulation of the excavation damaged zone around an
opening in brittle rock." International Journal of Rock Mechanics and Mining
Sciences 44.6 (2007): 835-845.
Gu, H., and X. Weng. "Criterion for fractures crossing frictional interfaces at non-orthogonal
angles." 44th US Rock Mechanics Symposium and 5th US-Canada Rock Mechanics
Symposium. American Rock Mechanics Association, 2010.
Haddad, Mahdi, and Kamy Sepehrnoori. "Simulation of hydraulic fracturing in quasi-brittle shale
formations using characterized cohesive layer: Stimulation controlling factors." Journal of
Unconventional Oil and Gas Resources 9 (2015): 65-83.
Hagoort, Jacques, Brian D. Weatherill, and Antonin Settari. "Modeling the propagation of
waterflood-induced hydraulic fractures." Society of Petroleum Engineers Journal 20.04
(1980): 293-303.
Haimson, B. and C. Fairhurst. 1967. Initiation and extension of hydraulic fractures in rocks. Soc.
Petr. Eng. J. 7:310-318.
Hanson, M. E., et al. "Some effects of stress, friction, and fluid flow on hydraulic fracturing." Old
SPE Journal 22.3 (1982): 321-332.
Hillerborg, Arne, Mats Modéer, and P-E. Petersson. "Analysis of crack formation and crack
growth in concrete by means of fracture mechanics and finite elements." Cement and
concrete research 6.6 (1976): 773-781.
Hoagland, Richard G., George T. Hahn, and Alan R. Rosenfield. "Influence of microstructure on
fracture propagation in rock." Rock Mechanics 5.2 (1973): 77-106.
Homand-Etienne, F., Dashnor Hoxha, and Jian-Fu Shao. "A continuum damage constitutive law
for brittle rocks." Computers and Geotechnics 22.2 (1998): 135-151.
Horsrud, P., R. Risnes and R.K. Bratli. 1982. Fracture initiation pressures in permeable poorly
consolidated sands. Int. J. Rock Mech. Min. Sci. Geomech. Abs. 19:255-266.
166

Hossain, Md Mofazzal, and M. K. Rahman. "Numerical simulation of complex fracture growth
during tight reservoir stimulation by hydraulic fracturing." Journal of Petroleum Science
and Engineering 60.2 (2008): 86-104.
Howard, George C., and C. R. Fast. "Optimum fluid characteristics for fracture
extension." Drilling and Production Practice (1957).
Huang, Jiefan, Zhongyan Wang, and Yonghong Zhao. "The development of rock fracture from
microfracturing to main fracture formation." International journal of rock mechanics and
mining sciences & geomechanics abstracts. Vol. 30. No. 7. Pergamon, 1993.
Irwin, G. R. "Analysis of stress and strain near the end of a crack traversing a plate". Journal of
Applied Mechanics, 24-361, 1957.
Jeffrey, R. G., and A. Settari. "A comparison of hydraulic fracture field experiments, including
mineback geometry data, with numerical fracture model simulations." SPE Annual
Technical Conference and Exhibition. 1995.
R. G., Jeffrey, et al. "An analysis of hydraulic fracture and mineback data for a treatment in the
German creek coal seam." SPE Rocky Mountain Regional Meeting. 1992.
Jeffrey, Robert, et al. "Measuring hydraulic fracture growth in naturally fractured rock." SPE
Annual Technical Conference and Exhibition. 2009.
Ji, Lujun, Antonin Settari, and Richard Sullivan. "A Novel Hydraulic Fracturing Model Fully
Coupled with Geomechanics and Reservoir Simulator." SPE Annual Technical Conference
and Exhibition. 2007.
Jin, Z-H., and C. T. Sun. "A comparison of cohesive zone modeling and classical fracture
mechanics based on near tip stress field." International journal of solids and
structures 43.5 (2006): 1047-1060.
Jing, Lanru. "A review of techniques, advances and outstanding issues in numerical modelling for
rock mechanics and rock engineering." International Journal of Rock Mechanics and
Mining Sciences 40.3 (2003): 283-353.
Kaiser, Peter K., et al. "Hydraulic Fracturing Mine Back Trials—Design Rationale and Project
Status." (2013).
Khodabakhshnejad, Arman, and Fred Aminzadeh, "An Extended Finite Element Method Based
Modeling of Hydraulic Fracturing." AAPG/SEG/SPWLA Hedberg Conference, December
7-11, 2014.
Khodabakhshnejad, A., and F. Aminzadeh. "Three Dimensional Modeling of Hydraulic Fracturing
Using a Coupling Approach." (to be submitted).
Khodaverdian, Mohamad, and Paul McElfresh. "Hydraulic fracturing stimulation in poorly
consolidated sand: Mechanisms and consequences." SPE Annual Technical Conference
and Exhibition. 2000.
167

Khristianovich, S. A., and Yu P. Zheltov. "Formation of vertical fractures by means of highly
viscous liquid." Proc. Vol. 4. 1955.
Kim, Jihoon, Hamdi A. Tchelepi, and Ruben Juanes. "Stability, accuracy and efficiency of
sequential methods for coupled flow and geomechanics." SPE Reservoir Simulation
Symposium. Society of Petroleum Engineers, 2009.
Krueger, Ronald. "Virtual crack closure technique: history, approach, and applications." Applied
Mechanics Reviews 57.2 (2004): 109-143.
Kuzkin, Vitaly A., and Anton M. Krivtsov Aleksandr M. Linkov. "Proppant transport in hydraulic
fractures: computer simulation of effective properties and movement of the
suspension." Proc. 41 Summer-School Conference" Advanced Problems in Mechanics.
2013.
Ladd, A. J. C., and R. Verberg. "Lattice-Boltzmann simulations of particle-fluid
suspensions." Journal of Statistical Physics 104.5-6 (2001): 1191-1251.
Lade, P. V., and M. K. Kim. "Single hardening constitutive model for frictional materials II. Yield
critirion and plastic work contours." Computers and geotechnics 6.1 (1988): 13-29.
Lawn, Brian. Fracture of brittle solids. Cambridge university press, 1993.
Li, Yanchao, et al. "Numerical Simulation of Hydraulically Induced Fracture Network Propagation
in Shale Formation." IPTC 2013: International Petroleum Technology Conference. 2013.
Lore, Jason S., Peter Eichhubl, and Atilla Aydin. "Alteration and fracturing of siliceous mudstone
during in situ combustion, Orcutt field, California." Journal of Petroleum Science and
Engineering 36.3 (2002): 169-182.
Louis, Claude. A study of groundwater flow in jointed rock and its influence on the stability of rock
masses. Imperial College of Science and Technology, 1969.
Lu, Yinlong, Derek Elsworth, and Lianguo Wang. "A dual-scale approach to model time-
dependent deformation, creep and fracturing of brittle rocks." Computers and
Geotechnics 60 (2014): 61-76.
Martys, Nicos S. "Study of a dissipative particle dynamics based approach for modeling
suspensions." Journal of Rheology 49 (2005): 401.
Martys, Nicos S., et al. "A smoothed particle hydrodynamics-based fluid model with a spatially
dependent viscosity: application to flow of a suspension with a non-Newtonian fluid
matrix." Rheologica acta 49.10 (2010): 1059-1069.
Mayerhofer, M. J., C. A. Ehlig-Economides, and M. J. Economides. "Pressure transient analysis
of fracture calibration tests." Journal of Petroleum Technology 47.03 (1995): 229-234.
McClure, Mark W., and Roland N. Horne. "Discrete fracture modeling of hydraulic stimulation in
enhanced geothermal systems." Proceedings of the 35th Workshop on Geothermal
Reservoir Engineering. 2010.
168

Medlin, W.L. and L. Masse. 1979. Laboratory investigation of fracture initiation pressure and
orientation. Soc. Petr. Eng. J. 129-144.
Melenk, Jens Markus, and Ivo Babuška. "The partition of unity finite element method: basic theory
and applications." Computer methods in applied mechanics and engineering 139.1 (1996):
289-314.
Meyer, B. R. "Three-dimensional hydraulic fracturing simulation on personal computers: theory
and comparison studies." SPE Eastern Regional Meeting. Society of Petroleum Engineers,
1989.
Meyer, Bruce Roman, and Lucas W. Bazan. "A discrete fracture network model for hydraulically
induced fractures-theory, parametric and case studies." SPE Hydraulic Fracturing
Technology Conference. Society of Petroleum Engineers, 2011.
Mohammadi, Soheil. Extended finite element method: for fracture analysis of structures. John
Wiley & Sons, 2008.
Morgan, William Edmund. "A fully implicit stochastic model for hydraulic fracturing based on
the discontinuous deformation analysis." (2014).
Neuman, Shlomo P. "Theoretical derivation of Darcy's law." Acta Mechanica25.3-4 (1977): 153-
170.
Nicksiar, Mohsen, and C. D. Martin. "Evaluation of Methods for Determining Crack Initiation in
Compression Tests on Low-Porosity Rocks." Rock mechanics and rock engineering 45.4
(2012): 607-617.
Nolte, Kenneth. "Fluid flow considerations in hydraulic fracturing." SPE Eastern Regional
Meeting. 1988.
Nordgren, R. P. "Propagation of a vertical hydraulic fracture." Old SPE Journal12.4 (1972): 306-
314.
Olson, J. E. "Multi-fracture propagation modeling: Applications to hydraulic fracturing in shales
and tight gas sands." The 42nd US rock mechanics symposium (USRMS). 2008.
Osher, Stanley, and James A. Sethian. "Fronts propagating with curvature-dependent speed:
algorithms based on Hamilton-Jacobi formulations." Journal of computational
physics 79.1 (1988): 12-49.
Papamichos E. Constitutive laws for geomaterials. Oil Gas Science and Technology Rev IFP 1999;
54(6):759–71
Perkins, T. K., and L. R. Kern. "Widths of hydraulic fractures." Journal of Petroleum
Technology 13.9 (1961): 937-949.
Pommier, Sylvie, et al. Extended finite element method for crack propagation. John Wiley & Sons,
2013.
169

Potluri, N. K., D. Zhu, and Alfred Daniel Hill. "The effect of natural fractures on hydraulic fracture
propagation." SPE European Formation Damage Conference. Society of Petroleum
Engineers, 2005.
Rajendran, A. M., et al. "Effects of strain rate on plastic flow and fracture in pure
tantalum." Journal of Materials Shaping Technology 9.1 (1991): 7-20.
Ranjan, Gopal, and A. S. R. Rao. Basic and applied soil mechanics. New Age International, 2007.
Remmers, J. J. C., Rene de Borst, and A. Needleman. "A cohesive segments method for the
simulation of crack growth." Computational mechanics 31.1-2 (2003): 69-77.
Raymond, S., A. Khodabakhshnejad, A. Ouenes, and F. Aminzadeh. "A comparison of the
Material Point Method and Finite Element Method in modelling cracks for hydraulic
fracturing processes." Computer Methods in Applied Mechanics and Engineering (to be
submitted).
Renshaw, C. E., and D. D. Pollard. "An experimentally verified criterion for propagation across
unbounded frictional interfaces in brittle, linear elastic materials." International journal of
rock mechanics and mining sciences & geomechanics abstracts. Vol. 32. No. 3. Pergamon,
1995.
Riahi, Azadeh, and Brank Damjanac. "Numerical study of interaction between hydraulic fracture
and discrete fracture network." ISRM International Conference for Effective and
Sustainable Hydraulic Fracturing. International Society for Rock Mechanics, 2013.
Rojek, Jerzy, and Eugenio Oñate. "Multiscale analysis using a coupled discrete/finite element
model." Interaction and Multiscale Mechanics 1.1 (2007): 1-31.
Romm, E. S. (1966), Hydraulic Properties of Fractured Rocks (in Russian), Nedra, Moscow.
Roylance, David. "Introduction to fracture mechanics." Department of Materials Science and
Engineering, Massachusetts Institute of Technology (2001).
Salehi, S., and R. Nygaard. "Numerical Study of Fracture Initiation, Propagation, Sealing to
Enhance Wellbore Fracture Gradient." 45th US Rock Mechanics/Geomechanics
Symposium. 2011.
Salimzadeh, Saeed, and Nasser Khalili. "A three-phase XFEM model for hydraulic fracturing with
cohesive crack propagation." Computers and Geotechnics 69 (2015): 82-92.
Saouma, Victor Edouard, et al. Interactive finite element analysis of reinforced concrete: A
fracture mechanics approach. No. 81/5. 1981.
Saouma, Victor E., and Mary Jane Kleinosky. "Finite element simulation of rock cutting: a fracture
mechanics approach." The 25th US Symposium on Rock Mechanics (USRMS). 1984.
Schutjens, P. M. T. M., et al. "Compaction-induced porosity/permeability reduction in sandstone
reservoirs." SPE Annual Technical Conference and Exhibition. 2001.
170

Schutjens, P. M. T. M., et al. "Compaction-induced porosity/permeability reduction in sandstone
reservoirs: Data and model for elasticity-dominated deformation." SPE Reservoir
Evaluation & Engineering 7.3 (2004): 202-216.
Segura, J. M., and I. Carol. "Zero-thickness interface elements for hydraulic fracture
simulation." XX Anales de la Mecánica de Fractura, Benicàssim, Spain (2003): 143-148.
Settari, A. "Simulation of hydraulic fracturing processes." Old SPE Journal 20.6 (1980): 487-500.
Settari, Antonin, and Michael Cleary. "Three-dimensional simulation of hydraulic
fracturing." Journal of petroleum technology 36.7 (1984): 1177-1190.
Settari, Antonin, and Michael Cleary. "Development and testing of a pseudo-three-dimensional
model of hydraulic fracture geometry." SPE Production Engineering 1.6 (1986): 449-466.
Sha, George T., and Chien-Tung Yang. "Determination of mixed mode stress intensity factors
using explicit weigh functions." Fracture mechanics. 18th symposium, ASTM STP. Vol.
945. 1988.
Shah, Surendra P. Fracture mechanics of concrete: applications of fracture mechanics to concrete,
rock and other quasi-brittle materials. John Wiley & Sons, 1995.
Sharp, John C., and Y. N. T. Maini. "Fundamental considerations on the hydraulic characteristics
of joints in rock." Proceedings Symposium on Percolation through Fissured Rock. 1972.
Sheibani, Farrokh, and Jon Olson. "Stress intensity factor determination for three-dimensional
crack using the displacement discontinuity method with applications to hydraulic fracture
height growth and non-planar propagation paths." ISRM International Conference for
Effective and Sustainable Hydraulic Fracturing. International Society for Rock
Mechanics, 2013.
Shet, C., and N. Chandra. "Analysis of energy balance when using cohesive zone models to
simulate fracture processes." Journal of engineering materials and technology 124.4
(2002): 440-450.
Simo, J. C., and J. W. Ju. "Strain-and stress-based continuum damage models—I.
Formulation." International journal of solids and structures 23.7 (1987): 821-840.
Sneddon, INa. "The distribution of stress in the neighbourhood of a crack in an elastic
solid." Proceedings of the Royal Society of London. Series A. Mathematical and Physical
Sciences 187.1009 (1946): 229-260.
Snow, David Tunison. A parallel plate model of fractured permeable media. Diss. University of
California, Berkeley, 1965.
Song, Jeong‐Hoon, Pedro Areias, and Ted Belytschko. "A method for dynamic crack and shear
band propagation with phantom nodes." International Journal for Numerical Methods in
Engineering 67.6 (2006): 868-893.
171

Stolarska, M., et al. "Modelling crack growth by level sets in the extended finite element
method." International journal for numerical methods in Engineering51.8 (2001): 943-
960.
Sulem, J., and H. Ouffroukh. "Hydromechanical behaviour of Fontainebleau sandstone." Rock
mechanics and rock engineering 39.3 (2006): 185-213.
Valencia, Karen Joy, et al. "Optimizing stimulation of coalbed methane reservoir using multi-stage
hydraulic fracturing treatment and integrated fracture modeling." SPE Asia Pacific Oil and
Gas Conference and Exhibition. 2005.
Vandamme, L., and J. H. Curran. "A three‐dimensional hydraulic fracturing
simulator." International Journal for Numerical Methods in Engineering 28.4 (1989): 909-
927.
Wang, HanYi. "Numerical modeling of non-planar hydraulic fracture propagation in brittle and
ductile rocks using XFEM with cohesive zone method." Journal of Petroleum Science and
Engineering 135 (2015): 127-140.
Wangen, Magnus. "Finite element modeling of hydraulic fracturing in 3D."Computational
Geosciences 17.4 (2013): 647-659.
Warpinski, N. R., and L. W. Teufel. "Influence of Geologic Discontinuities on Hydraulic Fracture
Propagation (includes associated papers 17011 and 17074)."Journal of Petroleum
Technology 39.2 (1987): 209-220.
Warpinski, N. R., and Michael Berry Smith. "Rock mechanics and fracture geometry." Recent
Advances in Hydraulic Fracturing 12 (1989): 57-80.
Weinberg, Roberto F., and Klaus Regenauer-Lieb. "Ductile fractures and magma migration from
source." Geology 38.4 (2010): 363-366.
Weng, Xiaowei, et al. "Modeling of hydraulic-fracture-network propagation in a naturally
fractured formation." SPE Production & Operations 26.4 (2011): 368-380.
Wilson, Charles R., and Paul A. Witherspoon. "Steady state flow in rigid networks of
fractures." Water Resources Research 10.2 (1974): 328-335.
WinterShall Company Website, at http://www.wintershall.com/technology/production-
technology/hydraulic-fracturing.html, Accessed on January 28, 2014.
Xu, Bin, and Ron Wong. "A 3D finite element model for history matching hydraulic fracturing in
unconsolidated sands formation." Journal of Canadian Petroleum Technology 49.4
(2010): 58-66.
Zhang, G. M., et al. "Three-dimensional finite element simulation and parametric study for
horizontal well hydraulic fracture." Journal of Petroleum Science and Engineering 72.3
(2010): 310-317.
172

Zhang, Xi, Marc Thiercelin, and Robert Jeffrey. "Effects of frictional geological discontinuities
on hydraulic fracture propagation." SPE Hydraulic Fracturing Technology Conference.
2007.
Zhao, Yonghong. "Crack pattern evolution and a fractal damage constitutive model for
rock." International Journal of Rock Mechanics and Mining Sciences35.3 (1998): 349-366.
Zhou, Jian, et al. "Analysis of fracture propagation behavior and fracture geometry using a tri-axial
fracturing system in naturally fractured reservoirs." International Journal of Rock
Mechanics and Mining Sciences 45.7 (2008): 1143-1152.
Zhou, Lei, et al. "Numerical modeling and investigation of fracture propagation with arbitrary
orientation through fluid injection in tight gas reservoirs with combined XFEM and
FVM." Environmental Earth Sciences 73.10 (2015): 5801-5813.
Zielonka, M.G., Searles, K.H., Ning, J. and Buechler, S.R.: “Development and Validation of Fully-
coupled Hydraulic Fracturing Simulation Capabilities.” SIMULIA Community
Conference, SCC2014, Providence, Rhode Island, May 19-21, 2014.
Zienkiewicz, Olgierd Cecil, et al. The finite element method. Vol. 3. London: McGraw-hill, 1977.
Zhu, W. C., and C. A. Tang. "Numerical simulation of Brazilian disk rock failure under static and
dynamic loading." International Journal of Rock Mechanics and Mining Sciences 43.2
(2006): 236-252.

Abstract (if available)

Abstract Maximizing hydrocarbon production from unconventional reservoirs requires the implementation of best practices in order to create a connected network of hydraulic and natural fractures. The desired result will be achieved if the operational hydraulic fracturing parameters are optimized. Determining the appropriate number of fracturing stages needed for peak performance demands a well calibrated model. Such a model allows for a more accurate prediction of fracture distribution. To this end, the industry currently employs several forward simulation tools. However, owing to the complexity of this problem, and the simplified assumptions made, existing methods often provide a poor representation of reality, leading to inaccurate predictions. To address these challenges, this work presents a novel, coupled approach to model hydraulically fractured subsurface cracks, providing a dynamic evolution of the subsurface fracture network. ❧ In this study, I introduce a new algorithm to model the initiation and propagation of hydraulic fracturing. It integrates rock deformation and fluid flow using a coupling approach. The algorithm enables us to accurately simulate fluid pressure inside the fracture and rock matrix and to calculate the in situ stress distribution within the reservoir. ❧ An in-house fluid flow simulator is developed based on a modified diffusivity equation that takes the effect of fracture width changes into account. I calculate fluid pressure using Finite Difference Method (FDM), and export it to the geo-mechanical simulator as a load. To calculate the resulting stress distribution inside the model and to predict fracture initiation and propagation, I apply the Extended Finite Element Method (XFEM). The workflow involves integrating the fluid flow and rock deformation equations

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

PDF

Modeling and simulation of complex recovery processes

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Integrated reservoir characterization for unconventional reservoirs using seismic, microseismic and well log data

PDF

A method for characterizing reservoir heterogeneity using continuous measurement of tracer data

$Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data$

PDF

Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data

PDF

Subsurface model calibration for complex facies models

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Stress and deformation analysis on fluid-exposed reservoir rocks

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers

Asset Metadata

Creator Khodabakhshnejad, Arman (author)

Core Title An extended finite element method based modeling of hydraulic fracturing

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 07/19/2016

Defense Date 12/02/2015

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag coupled approach,fluid flow,hydraulic fracturing,numerical modeling,OAI-PMH Harvest,XFEM

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Aminzadeh, Fred (committee chair), Ershaghi, Iraj (committee member), Sammis, Charles (committee member)

Creator Email khodabak@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-270650

Unique identifier UC11279630

Identifier etd-Khodabakhs-4556.pdf (filename),usctheses-c40-270650 (legacy record id)

Legacy Identifier etd-Khodabakhs-4556.pdf

Dmrecord 270650

Document Type Dissertation

Format application/pdf (imt)

Rights Khodabakhshnejad, Arman

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA