Hydraulic Fracture
|
A Hydraulic fractures (HF) are a class of brittle fractures that propagate in pre-stressed solid media due to the injection of a viscous fluid. These fractures occur naturally when pressurized magma from deep underground chambers form vertical intrusions driven by buoyancy forces, which can in turn form geological structures such as dykes and sills. In the oil and gas industry, HF are deliberately created in reservoirs to enhance the recovery of hydrocarbons by the creation of permeable pathways . The application of HF in geotechnical engineering is growing. For example, in the mining industry they have recently been used to weaken the rock and enhance the so-called block-caving process. Similarly, the extraction of geothermal energy equires the creation of new fractures to increase the existing fracture networks. Hydraulic fractures have also been used for waste disposal and are likely to play a key role in the gravitational trapping of CO2 in deep, low permeability, ocean sediments. The propagation of hydraulic fractures into undesirable locations can have severe safety consequences in the mining industry, can cause considerable loss of hydrocarbons and environmental damage in the oil industry. Likewise, the perforation of the caprock by HF can reverse the costly capture process involved in CO2 sequestration. It is therefore of considerable importance to have accurate models in order to be able to predict the advance of HF to achieve an effective design of the engineering parameters in the injection process. Mathematical models of HF involve a degenerate system of hyper-singular integro-partial differential equations defined on a domain with a moving boundary. The integral operator may be classified as a fractional Laplacian, which when combined with a free boundary problem forms is a class of problems that has received considerable attention in the theoretical PDE community (see for example the following link to a BIRS Workshop dedicated to this topic). The strong degeneracy of the lubrication equation yields solutions with singular pressure fields close to the tips or near "pinch points," which can occur within the fracture due to local changes in the geological structure. This class of problem has been shown to exhibit a rich multi-scale structure in which the fracture can propagate in a number of different modes each determined by the dominant physical process active at the tip of the fracture. The resolution of these problems requires a unique blend of numerical and analytic techniques to capture the free boundary of the evolving fracture as well as the complex coupling between the fluid which drives the fracture and the ambient elastic medium in which the fracture is being formed. The very different time scales on which the fluid flow and elastic response occur, poses considerable interesting challenges for numerical simulation. My interest in this research is the use asymptotic analysis to identify multi-scale singular solutions and to develop numerical tools to accommodate and exploit these asymptotic solutions.
|
by E.V. Dontsov and A. P. Peirce
Engineering Fracture mechanics, 142, p 116-139, 2015.
Comparison between the results from the classical P3D model, the enhanced P3D model (EP3D), and those of the fully planar simulator ILSA. For the zero toughness case KIc=0 the values are shown in blue and those of the experiment by Jeffrey and Bunger (2007) are depicted by the blue circles. The results shown in red and orange correspond to the toughness KIc=0.94 MPa.m^(1/2). |
Abstract: The goal of this paper is to develop a more accurate pseudo-3D
model for hydraulic fracturing. The primary weaknesses of the classical
pseudo-3D model are: (1) its inability to capture viscous height growth
and (2) its failure to include lateral fracture toughness. These flaws are addressed respectively by: (1) introducing an apparent
fracture toughness |
28. Proppant transport in hydraulic fracturing: Crack tip screen-out in KGD and P3D models
by E.V. Dontsov and A. P. Peirce
Int. J. Solids & Structures, 63, p 206-218, 2015.
|
Abstract: The aim of this study is to develop a model
for proppant transport in hydraulic fractures capable of capturing both gravitational settling and tip screen-out effects, while prohibiting the particles from reaching the crack tips by imposing a width restriction based on the particle size. First, the equations that govern the propagation of hydraulic fractures and the proppant transport inside them are formulated. They are based on the solution for the steady flow of a viscous fluid, mixed with spherical particles, in a channel, which is obtained assuming an empirical constitutive model. This proppant transport model is applied to two fracture geometries -- Khristianovich-Zheltov-Geertsma-De Klerk (KGD) and pseudo-3D (P3D). Numerical simulations show that the proposed method makes it possible to capture proppant plug formation and growth, as well as the gravitational settling for both geometries. A dimensionless parameter, whose magnitude reflects the intensity of the settling, is introduced for the P3D fracture. |
27. Slurry flow, gravitational settling and a proppant transport model for hydraulic fractures,
by E.V. Dontsov and A. P. Peirce
J. Fluid Mech. V 760, pp 567-590, 2014.
|
Abstract: The goal of this study is to analyse the steady flow of a Newtonian fluid mixed with spherical particles in a channel for the purpose of modelling proppant transport with gravitational settling in hydraulic fractures. The developments are based on a continuum constitutive model for a slurry, which is approximated by an empirical formula. It is shown that the problem under consideration features a two-dimensional flow and a boundary layer, which effectively introduces slip at the boundary and allows us to describe a transition from Poiseuille flow to Darcys law for high proppant concentrations. The expressions for both the outer (i.e. outside the boundary layer) and inner (i.e. within the boundary layer) solutions are obtained in terms of the particle concentration, particle velocity and fluid velocity. Unfortunately, these solutions require the numerical solution of an integral equation, and, as a result, the development of a proppant transport model for hydraulic fracturing based on these results is not practicable. To reduce the complexity of the problem, an approximate solution is introduced. To validate the use of this approximation, the error is estimated for different regimes of flow. The approximate solution is then used to calculate the expressions for the slurry flux and the proppant flux, which are the basis for a model that can be used to account for proppant transport with gravitational settling in a fully coupled hydraulic fracturing simulator. |
26. Enrichment strategies and convergence properties of the XFEM for hydraulic fracture problems,
by E. Gordeliy and A. P. Peirce,
Comp. Meth. in Appl. Mech. and Eng., 283,
p474-502, 2015.
Wronskian of the new enrichment basis functions |
Abstract: In two recent papers (Gordeliy and Peirce, 2013; Gordeliy and Peirce, 2013) investigating the use of Extended Finite Element Method (XFEM) for modeling hydraulic fractures (HF), two classes of boundary value problem and two distinct enrichment types were identified as being essential components in constructing successful XFEM HF algorithms. In this paper we explore the accuracy and convergence properties of these boundary value formulations and enrichment strategies. In addition, we derive a novel set of crack-tip enrichment functions that enable the XFEM to model HF with the full range of power law r^{\lambda} behavior of the displacement field and the corresponding r^{\lambda-1} singularity in the stress field, for \frac{1}{2}\le\lambda<1. This novel crack-tip enrichment enables the XFEM to achieve the optimal convergence rate, which is not achieved by existing enrichment functions used for this range of power law. The two XFEM boundary value problem classes are as follows: i) a Neumann to Dirichlet map in which the pressure applied to the crack faces is the specified boundary condition and the XFEM is used to solve for the corresponding crack width (P->W); and ii) a mixed hybrid formulation of the XFEM that makes it possible to incorporate the singular behavior of the crack width in the fracture tip and uses a pressure boundary condition away from it (P&W). The two enrichment schemes considered are: i) the XFEM-t scheme with full singular crack-tip enrichment and ii) a simpler, more efficient, XFEM-s scheme in which the singular tip behavior is only imposed in a weak sense. If enrichment is applied to all the nodes of tip-enriched elements, then the resulting XFEM stiffness matrix is singular due to a linear dependence among the set of enrichment shape functions, which is a situation that also holds for the classic set of square-root enrichment functions. For the novel set of enrichment functions we show how to remove this rank deficiency by eliminating those enrichment shape functions associated with the null space of the stiffness matrix. Numerical experiments indicate that the XFEM-t scheme, with the new tip enrichment, achieves the optimal O(h^{2}) convergence rate we expect of the underlying piece-wise linear FEM discretization, which is superior to the enrichment functions currently available in the literature for \frac{1}{2}<\lambda<1. The XFEM-s scheme, with only signum enrichment to represent the crack geometry, achieves an O(h) convergence rate. It is also demonstrated that the standard P->W formulation, based on the variational principle of minimum potential energy, is not suitable for modeling hydraulic fractures in which the fluid and the fracture fronts coalesce, while the mixed hybrid P\&W formulation based on the Hellinger-Reissner variational principle does not have this disadvantage. |
by A. Peirce
Comp. Meth. in Appl. Mech and Eng., 283, p881-908, 2015.
|
Abstract:
In this paper we describe an implicit level set algorithm (ILSA)
suitable for modeling multi-scale behavior in planar hydraulic fractures
propagating in three dimensional elastic media. This multi-scale behavior
is typically encountered when multiple physical processes compete to determine
the location of the fracture free boundary. Instead of having to match the
mesh size to the finest active length scale, or having to re-mesh as the
dominant length scales change in space and time, the novel ILSA scheme is
able to represent the required multi-scale behavior on a relatively coarse rectangular mesh. This is achieved by using the local front velocity to construct, for each point of a set of control points, a mapping that adaptively identifies the dominant length scale at which the appropriate multi-scale universal asymptotic solution needs to be sampled. Finer-scale behavior is captured in a weak sense by integrating the universal asymptotic solution for the fracture width over partially filled tip elements and using these integrals to set the average values of the widths in all tip elements. The ILSA solution shows good agreement with a multi-scale reference solution comprising a radial solution that transitions from viscosity to toughness dominated ropagation regimes. The ILSA scheme is also used to model blade-like hydraulic fractures that break through stress barriers located symmetrically with respect to the injection point. For the zero toughness case, the ILSA solutionshows close agreement to experimental results. The multi-scale ILSA scheme is also used to provide results when the material toughness K_{Ic} is non-zero. In this case different parts of the fracture-free-boundary can be propagating in different regimes. It is hoped that the multi-scale ILSA solutions presented here will form a set of reference results that can be used to benchmark simulators that use a propagation criterion based on only one dissipative process (either toughness or viscosity). The multi-scale ILSA solutions at larger times (for which plane strain conditions develop in vertical cross sections) are compared with and show close agreement to plane strain exact solutions for height-growth and the fracture width in vertical cross sections. This comparison provides some measure of the accuracy of the multi-scale ILSA scheme. The multi-scale ILSA solutions are also used to identify the regimes of applicability of pseudo 3D (P3D) approximate solutions. These ILSA solutions can also be used to design improved P3D models. |
24. A New Technique for Proppant Schedule Design,
by E. V. Dontsov and A. P. Peirce
Hydraulic Fracturing Journal, V1, No. 3, 2014.
Proppant Schedule and placment
|
Abstract: This study
introduces a novel methodology for the design of the proppant pumping
schedule for a hydraulic fracture, in which the final proppant distribution
along the crack is prescribed. The method is based on the assumption that
the particles have relatively small impact on the fracture propagation,
unless they reach the tip region. This makes it possible to relate the
proppant velocity to the clear fluid velocity inside the fracture, which
is calculated assuming no proppant. Having the history of the clear fluid
velocity distribution, the prospective proppant motion can be computed.
Then, volume balance is used to relate the final concentration at some
point inside the fracture to the corresponding input concentration at
a specic time instant, which helps to avoid solving an inverse problem.
One exceptional feature of the approach lies in the fact that it is applicable
to multiple fracture geometries and can be implemented using various hydraulic
fracturing simulators. To verify the technique, two fracture geometries
are considered - Khristianovich-Zheltov-Geertsma-De Klerk (KGD) and pseudo-3D |
23. On
the moving boundary conditions for a Hydraulic Fracture
by E. Detournay and A. Peirce
International Journal of Engineering Science, 84, p 147155, 2014.
Evolving Fluid and Fracture Fronts |
Abstract: This paper re-examines the boundary conditions at the moving front of a hydraulic fracture when the fluid front has coalesced with the crack edge. This practically important particular case is treated as the zero fluid lag limit of the general case when the two fronts are distinct. The limiting process shows what becomes of the two boundary conditions on the fluid front: a pressure condition and a Stefan condition equating the front velocity to the average fluid velocity, when the lag vanishes. On the one hand, the pressure condition disappears as the net pressure (the difference between the fluid pressure and the magnitude of the far-field stress normal to the fracture) becomes singular. On the other hand, the Stefan condition transforms into a zero flux boundary condition at the front, and, as a consequence, the velocity of the coalesced front does not appear explicitly in the boundary conditions. However, the front velocity can still be extracted from the near-tip aperture field by a nonlinear asymptotic analysis. The paper concludes with a description of an algorithm to propagate the combined front, which explicitly uses the known multiscale asymptotics of the fracture aperture. |
by A.P. Bunger and A. P. Peirce, submitted to
the ASCE Shale Energy Conference, Pittsburgh, 21-23 July, 2014.
|
Abstract:
The technique of multistage hydraulic fracturing from horizontal
wells is universally credited with enabling the economical production of
hydrocarbon resources from shale formations. The method almost always entails
the injection of fluid through the wellbore with the potential to create
hydraulic fractures from multiple reservoir entry points, typically clusters
of wellbore perforations, that are spaced out along the wellbore within
a section that is colloquially referred to as a stage. Arguably
the most basic question about this situation is how many perforation clusters
within a given fracturing stage can be expected to produce growing hydraulic
fractures. This paper presents a numerical investigation of this issue that
employs a newly-developed, fully coupled parallel planar 3D hydraulic fracturing simulator that features: implicit time stepping, an implicit level set scheme to locate the propagating hydraulic fracture fronts that respond to their regimes of propagation and enables highly accurate simulations using a very coarse mesh, and the capability to dynamically partition the fluid among multiple, simultaneously growing hydraulic fractures in parallel, overlapping planes. Our results demonstrate the dependence of the energetically preferred number of growing hydraulic fractures on the length of the isolated zone, the height of the reservoir, and the relative importance of the fluid viscosity. In particular, we show that reservoirs with effective height containment and injection strategies that ensure substantial viscous dissipation will promote growth of multiple simultaneous hydraulic fractures rather than localization to just one or two dominant fractures. |
(link to Posting on UBC Circle Repository)
by A. P. Peirce and A.P. Bunger, Submitted to the SPE Journal on
the 17 October, 2013.
|
Abstract:
One of the important hurdles in horizontal well stimulation is the
generation of hydraulic fractures (HFs) from all perforation clusters within
a given stage despite the challenges posed by stress shadowing and reservoir
variability. In this paper we use a newly-developed, fully coupled, parallel-planar
3D HF model to investigate the potential to minimize the negative impact
of stress shadowing and thereby to promote more uniform fracture growth
across an array of HFs by adjusting the location of the perforation clusters.
In this model the HFs are assumed to evolve in an array of parallel planes
with full 3D stress coupling while the constant fluid influx into the wellbore
is dynamically partitioned to each fracture so that the wellbore pressure
is the same throughout the array. The model confirms the phenomenon of inner
fracture suppression due to stress shadowing when the perforation clusters
are uniformly distributed. Indeed, the localization of the fracture growth
to the outer fractures is so dominant that the total fractured area generated
by uniform arrays is largely independent of the number of perforation clusters.
However, numerical experiments indicate that certain non-uniform cluster
spacings promote a profound improvement in the even development of fracture
growth. Identifying this effect relies on this new models ability
to capture the full hydro-dynamical coupling between the simultaneously
evolving HFs in their transition from radial to PKN-like geometries. |
20. Implicit level set schemes for modeling hydraulic fractures using the XFEM
by E. Gordeliy and A. P. Peirce
Comp. Meth. in Appl. Mech. and Eng., 266, p125-143, 2013.
|
Abstract: We describe two novel XFEM schemes for modeling fluid driven fractures both of which exploit an implicit level set algorithm (ILSA) for locating the singular free boundary that occurs when the fluid and fracture fronts coalesce. Both schemes use the mixed P&W XFEM formulation developed in (Gordeliy and Peirce, 2013) to incorporate the singular asymptotic solution in the fracture tips. The proposed level set strategy also exploits the asymptotic solution to provide a robust procedure to locate the free boundary, which is not restricted to symmetric growth of the fracture geometry or to a particular mode of propagation. The versatility of the ILSA-XFEM scheme is demonstrated by sampling different asymptotic behaviors along the so-called MK edge of parameter space (Detournay, 2004) by making use of a universal asymptote (Garagash, 1998; Garagash and Detournay, 2000). The two ILSA-XFEM schemes differ in the enrichment strategies that they use to represent the fracture tips: a scheme with full tip enrichment and a simpler, more efficient, scheme in which the tip asymptotic behavior is only imposed in a weak sense. Numerical experiments indicate that the XFEM-t scheme, with full tip enrichment, achieves an O(h^{2}) asymptotic convergence rate, while the XFEM-s scheme, with only signum enrichment to represent the crack geometry, achieves an O(h) asymptotic convergence rate. |
by Lecampion B., Peirce A.P, Detournay E., Zhang X., Chen Z., Bunger
A.P., Detournay C., Napier J., Abbas S., Garagash D., Cundall, P.,
In: Effective and Sustainable Hydraulic Fracturing, AP Bunger, J McLennan and R Jeffrey (eds.), ISBN 978-953-51-1137-5, (Intech), Chapter 43, p855-873, 2013.
Evolution of the relative error in the fracture radius as a function of the mesh size for the different simulators - radial fracture in the viscosity dominated regime (i.e. zero toughness / early time, t < 1). All simulations displayed here are for t < 10E-6. |
We
benchmark a series of simulators against available reference solutions for
propagating plane-strain and radial hydraulic fractures. In particular,
we focus on the accuracy and convergence of the numerical solutions in the
important practical case of viscosity dominated propagation. The simulators are based on different propagation criteria: linear elastic fracture mechanics (LEFM), cohesive zone models/tensile strength criteria, and algorithms accounting for the multi-scale nature of hydraulic fracture propagation in the near-tip region. All the simulators tested here are able to capture the analytical solutions of the different configurations tested, but at vastly different computational costs. Algorithms based on the classical LEFM propagation condition require a fine mesh in order to capture viscosity dominated hydraulic fracture evolution. Cohesive zone models, which model the fracture process zone, require even finer meshes to obtain the same accuracy. By contrast, when the algorithms use the appropriate multi-scale hydraulic fracture asymptote in the near-tip region, the exact solution can be matched accurately with a very coarse mesh. The different analytical reference solutions used in this paper provide a crucial series of benchmark tests that any successful hydraulic fracturing simulator should pass. |
18. Coupling schemes for modeling hydraulic fracture propagation using the XFEM.
by E. Gordeliy and A. P. Peirce
Comp. Meth. in Appl. Mech. and Eng., 266, p125-143, 2013.
Top: The crack path obtained using a DD cde OribiC,
XFEM-t, and XFEM |
Abstract:
We describe coupled algorithms that use the Extended Finite Element
Method (XFEM) to solve the elastic crack component of the elasto-hydrodynamic
equations that govern the propagation of hydraulic fractures in an elastic
medium. With appropriate enrichment, the XFEM resolves the Neumann to Dirichlet
(ND) map for crack problems with O(h^{2}) accuracy and the Dirichlet to
Neumann (DN) map with O(h) accuracy. For hydraulic fracture problems with
a lag separating the fluid front from the fracture front, we demonstrate
that the finite pressure field makes it possible to use a scheme based on
the O(h^{2}) XFEM solution to the ND map. To treat problems in which there
is a coalescence of the fluid and fracture fronts, resulting in singular
tip pressures, we developed a novel mixed algorithm that combines the tip
width asymptotic solution with the O(h^{2}) XFEM solution of the ND map
away from the tips. Enrichment basis functions required for these singular
pressure fields correspond to width power law indices \lambda >=\frac{1}{2),
which are different from the index \lambda=\frac{1}{2} of linear elastic fracture mechanics. The solutions obtained from the new coupled XFEM schemes agree extremely well with those of published reference solutions. |
17. An
Integrated Extended Kalman Filter Implicit Level Set Algorithm for monitoring
Planar Hydraulic Fracture.
by A.P. Peirce and F. Rochinha
Inverse Problems, 28,(2012) 015009, (22 pp).
|
Abstract:
We describe a novel approach to the inversion of elasto-static tiltmeter
measurements to monitor a planar hydraulic fracture propagating within a
three dimensional elastic medium. The technique involves integrating the
Extended Kalman Filter (EKF) for state prediction and updates using tiltmeter
measurement time-series combined with a novel Implicit Level Set Algorithm
(ILSA) for the solution of the coupled elasto-hydrodynamic equations in order tolocate the unknown fracture free boundary. By means of a scaling argument a strategy is derived to tune the parameters in the algorithm to enable measurement information to compensate for significant un-modeled dynamics. A sequence of numerical experiments are considered in which significant changes are introduced to the fracture geometry by altering the confining geological stress field. Even though there is no confining stress field in the dynamic model used by the new EKF-ILSA scheme, it is able to use synthetic measurements from tiltmeters in order to arrive at remarkably accurate predictions of the fracture opening and the fracture footprint. These experiments also explore the robustness of the algorithm to noise and to placement of tiltmeter arrays operating in the near-field and far-field regimes. In these experiments the appropriate parameter choices and strategies to improve the robustness of the algorithm to significant measurement noise are explored. |
Hermite Basis Functions and the stresses they generate
|
Abstract: We describe a novel cubic Hermite collocation scheme for the solution of the coupled integro-partial differential equations governing the propagation of a hydraulic fracture in a state of plane strain. Special blended cubic Hermite-power law basis functions, with arbitrary index 0<a<1, are developed to treat the singular behavior of the solution that typically occurs at the tips of a hydraulic fracture. The implementation of blended infinite elements to model semi-infinite crack problems is also described. Explicit formulae for the integrated kernels associated with the cubic Hermite and blended basis functions are provided. The cubic Hermite collocation algorithm is used to solve a number of different test problems with two distinct propagation regimes and the results are shown to converge to published similarity and asymptotic solutions. The convergence rate of the cubic Hermite scheme is determined by the order of accuracy of the tip asymptotic expansion as well as the O(h^4) error due to the Hermite cubic interpolation. The errors due to these two approximations need to be matched in order to achieve optimal convergence. Backward Euler time-stepping yields a robust algorithm that, along with geometric increments in the time-step, can be used to explore the transition between propagation regimes over many orders of magnitude in time. |
|
Abstract: There is considerable interest in using remote elastostatic deformations to identify the evolving geometry of underground fractures that are forced to propagate by the injection of a high pressure viscous fluid. These so-called hydraulic fractures are used to increase the permeability in oil and gas reservoirs as well as to pre-fracture ore-bodies in preparation for a mining process known as block-caving. The undesirable intrusion of these hydraulic fractures into environmentally sensitive areas or into regions in mines which might pose safety hazards, has stimulated the search for techniques to enable the evolving hydraulic fracture geometries to be monitored. Previous approaches to this problem have involved the inversion of the elastostatic data at isolated time steps in the time series provided by tiltmeter measurements of the displacement gradient field at selected points in the elastic medium. At each time step, parameters in simple static models of the fracture (e.g., a single displacement discontinuity or low order moments of the fracture opening) are identified. The approach adopted in this paper is not to regard the sequence of sampled elastostatic data as independent, but rather to treat the data as linked by the coupled elastic-lubrication equations that govern the propagation of the evolving hydraulic fracture. We combine the Extended Kalman Filter (EKF) with features of a recently developed implicit numerical scheme to solve the coupled free boundary problem in order to form a novel algorithm to identify the evolving fracture geometry. Numerical experiments demonstrate that, despite excluding significant physical processes in the forward numerical model that substantially alter the geometry of the evolving hydraulic fracture, the EKF-numerical algorithm is able to compensate for the un-modeled processes by using the information fed back from tiltmeter data. Indeed, the proposed algorithm is able to provide reasonably faithful estimates of the fracture geometry, which are shown to converge to the actual hydraulic fracture geometry as the number of tiltmeters is increased. The location of tiltmeters can affect the resolution of the technique, which opens the possibility of using the algorithm to design the deployment of tiltmeters to optimize the resolution in regions of particular interest. |
The MATLAB code that implements the 4 th Order Collocation Scheme is given below. or |
Abstract:
This paper deals with the so-called pseudo three-dimensional
(P3D) model for a hydraulic fracture with equilibrium height growth across
two symmetric stress barriers. The key simplifying assumptions behind the
P3D model, which was originally introduced by A. Settari and M. P. Cleary
(Development and Testing of a Pseudo-Three-Dimensional Model of Hydraulic
Fracture Geometry (P3DH), SPE 10505, Proc. 6th SPE Symposium on Reservoir
Simulation, pp. 185-214, 1982) and by I. D. Palmer and H. B. Carrol (Numerical
Solution for Height and Elongated Hydraulic Fractures, SPE/DOE 11627, Proc.
SPE/DOE Symposium on Low Permeability, pp. 249-878, 1983), are that (i)
each cross section perpendicular to the main propagation direction is in a condition of plane-strain, and (ii) the local fracture height is determined by a balance between the effect of the stress jump across the barriers and that of the rock toughness. Furthermore, in the equilibrium height growth P3D models, the pressure is assumed uniform in each vertical cross-section. We revisit this particular model by first formulating the non-linear differential equations governing the evolution of the length, height, and aperture of the hydraulic fracture, in contrast to the numerical formulations adopted in many previous studies. Scaling of these equations shows that the solution depends, besides the dimensionless space and time coordinates, on only two numbers representing a scaled toughness and a scaled leak-off coefficient. Analysis of the governing equations enables us to determine explicitly the conditions under which breakthrough takes place (i.e., the onset of growth into the bounding layers), as well as the conditions of unstable height growth (i.e., the conditions of runaway height when the main assumptions of the equilibrium height model become invalid). The mathematical model is solved numerically using a novel implicit fourth order collocation scheme on a moving mesh, which makes explicit use of the fracture tip asymptotics. A complete listing of the MATLAB code for this algorithm is provided in an appendix and can be copied for experimentation. We then report the results of several numerical simulations conducted for different values of the dimensionless toughness and the dimensionless leak-off coefficients, as well as a comparison with closed-form small and large time similarity solutions that are valid under conditions where the fracture remains contained within the reservoir layer. Finally, we discuss possible extensions of the model, including layers of different elastic moduli and power law fluids. |
|
Abstract:
We present a uniform asymptotic solution (UAS) for a displacement
discontinuity (DD) that lies within the middle layer of a three-layer elastic
medium in which relative shear deformation between parallel interfaces is
controlled by linear springs. The DD is assumed to be normal to the two
interfaces between the elastic media. Using the Fourier transform (FT) method
we construct the leading term in the asymptotic expansion for the spectral
coefficient functions for a DD in a three layer-spring medium. Although
a closed form solution will require a solution in terms of an infinite series,
we demonstrate how this UAS can be used to construct highly efficient and
accurate solutions even in the case in which the DD actually touches the
interface. We compare the results using the Green's function UAS solution
for a crack crossing a soft interface with results obtained using a multi-layer
boundary element method. We also present results from an implementation
of the UAS Green's function approach in a pseudo 3D hydraulic fracturing
simulator to analyze the effect of interface shear deformation on the fracture
propagation process. These results are compared with field measurements. |
Exact and numerical Width
and Pressure solutions
|
Abstract: The coupled equations for a propagating hydraulic fracture exhibit a multi-scale structure for which resolution on all length scales becomes computationally prohibitive. Asymptotic analysis is able to identify the dominant physical process active at the computational length scale. This paper describes a novel algorithm which uses weak-form tip asymptotics on a rectangular Eulerian mesh to solve the problem of a propagating hydraulic fracture. The location of the fracture front is determined within each tip element by matching the volume associated with the known asymptotic solution to the flux of fluid into the given element. Even if the fracture front is curved, the algorithm is able to capture the solution on a relatively coarse rectangular mesh by implementing a weak-form of the tip asymptotic solution, based on averaging the volume over a tip element to determine the fracture widths at tip element centres. The fracture is divided into a ``channel'' region made up of elements that are filled with fluid and a ``tip'' region comprising partially filled elements. An iterative procedure is used to determine the fracture widths and fluid pressures in the channel region and the fracture font locations in the tip regions. The algorithm is tested by analyzing a hydraulic fracture propagating in a viscosity dominated regime within an impermeable homogeneous elastic material for which a similarity solution is available. The numerical results show close agreement to the exact solution even though a relatively coarse mesh is used. |
Movies
Comparisons with similarity solutions Viscosity dominated fracture propagation
Toughness dominated fracture propagation
|
Abstract: We describe a novel implicit level set algorithm to locate the free boundary for a propagating hydraulic fracture. A number of characteristics of the governing equations for hydraulic fractures and their coupling present considerable challenges for numerical modeling, namely: the degenerate lubrication equation; the hypersingular elastic integral equation; the indeterminate form of the velocity of the unknown fracture front, which precludes the implementation of established front evolution strategies that require an explicit velocity field; and the computationally prohibitive cost of resolving all the length scales. An implicit algorithm is also necessary for the efficient solution of the stiff evolution equations that involve fully populated matrices associated with the coupled non-local elasticity and degenerate lubrication equations. The implicit level set algorithm that we propose exploits the local tip asymptotic behavior, applicable at the computational length scale, in order to locate the free boundary. Local inversion of this tip asymptotic relation yields the boundary values for the Eikonal equation, whose solution gives the fracture front location as well as the front velocity field. The efficacy of the algorithm is tested by comparing the level set solution to analytic solutions for hydraulic fractures propagating in a number of distinct regimes. The level set algorithm is shown to resolve the free boundary problem with first order accuracy. Further it captures the field variables, such as the fracture width, with the first order accuracy that is consistent with the piecewise constant discretization that is used. |
|
Abstract: We derive a novel integral equation relating the fluid pressure in a fingerlike hydraulic fracture to the fracture width. By means of an asymptotic analysis in the small height to length ratio limit we are able to establish the action of the integral operator for receiving points that lie within three distinct regions: (1) an outer expansion region in which the dimensionless pressure is shown to be equal to the dimensionless width plus a small correction term that involves the second derivative of the width, which accounts for the nonlocal effects of the integral operator. The leading order term in this expansion is the classic local elasticity equation in the PKN model that is widely used in the oil and gas industry; (2) an inner expansion region close to the fracture tip within which the action of the elastic integral operator is shown to be the same as that of a finite Hilbert transform associated with a state of plane strain. This result will enable pressure singularities and stress intensity factors to be incorporated into analytic models of these finger-like fractures in order to model the effect of material toughness; (3) an intermediate region within which the action of the Fredholm integral operator of the first kind is reduced to a second kind operator in which the integral term appears as a small perturbation which is associated with a convergent Neumann series. These results are important for deriving analytic models of finger-like hydraulic fractures that are consistent with linear elastic fracture mechanics. |
|
Abstract: Hydraulic fracturing involves the propagation of a fracture in brittle rock by the intrusion of a high pressure viscous fluid. There is considerable interest in identifying characteristics of these evolving underground fractures via the passive monitoring of remote elastostatic deformations. In this paper, we present a far-field multipole expansion procedure to identify the harmonic moments of the fracture. The harmonic moments are related to fundamental quantities such as fracture volume and fracture asymmetries. We illustrate the efficacy of the multipole moment expansion technique by inverting synthetic displacement data from a hydraulic fracture simulator in order to identify the harmonic moments up to second order. These results are compared to those obtained by identifying the parameters of a dislocation model with a prescribed geometrya procedure which is commonly used for such problems. The multipole moment expansion technique has the following features: it provides significantly more accurate fracture volume information; it provides accurate estimates of first-order moments that can be used to identify asymmetric fractures; it is possible to adapt the truncation process to optimize the information content of a given set of measurements; it can, in some cases, provide estimates of the higher order moments which can be used to determine geometric attributes of the fracture. Given this last feature, we explore the possibility of using up to second-order harmonic moments to identify the dimensions of a simple polygonal model of the fracture footprint. This procedure is tested by attempting to identify fracture footprints from synthesized hydraulic fracture data. |
Propppant concentration (left) before and after shut-in and the corresponding fracture widths (right) |
Abstract: We provide a brief historical background of the development of hydraulic fracturing models for use in the petroleum and other industries. We discuss scaling laws and the propagation regimes that control the growth of hydraulic fractures from the laboratory to the field scale. We introduce the mathematical equations and boundary conditions that govern the hydraulic fracturing process, and discuss numerical implementation issues including: tracking of the fracture footprint, the control of the growth of the hydraulic fracture as a function of time, coupling of the equations, and time-stepping schemes. We demonstrate the complexity of hydraulic fracturing by means of an application example based on real data. Finally, we highlight some key areas of research that need to be addressed in order to improve current models. |
|
Abstract:
The dynamics of hydraulic fracture, described by a system of nonlinear
integrodifferential equations, is studied through the development and application
of a multiparameter singular perturbation analysis. We present a new single
expansion framework which describes the interaction between several physical
processes, namely viscosity, toughness, and leak-off. The problem has nonlocal
and nonlinear effects which give a complex solution structure involving
transitions on small scales near the tip of the fracture. Detailed solutions
obtained in the crack tip region vary with the dominant physical processes.
The parameters quantifying these processes can be identified from critical
scaling relationships, which are then used to construct a smooth solution
for the fracture depending on all three processes. Our work focuses on plane
strain hydraulic fractures on long time scales, and this methodology shows
promise for related models with additional time scales, fluid lag, or different
geometries, such as radial (penny-shaped) fractures and the classical PerkinsKernNordgren (PKN) model. |
|
Abstract:
This paper presents a novel asymptotic framework to obtain detailed
solutions describing the propagation of hydraulic fractures in an elastic
material. The problem consists of a system of nonlinear integro-differential
equations and a free boundary problem. This combination of local and nonlocal
effects leads to transitions on a small scale near the crack tip, which
control the behavior across the whole fracture profile. These transitions
depend upon the dominant physical process(es) and are identified by simultaneously
scaling the associated parameters with the distance from the tip. A smooth
analytic solution incorporating several physical processes in the crucial tip region can be constructed using this new framework. In order to clarify the exposition of the new methodology, this paper is confined to considering the impermeable case in which only the two physical processes of viscous dissipation and structure energy release compete. |
|
Abstract: We discuss the properties of a class of sparse localized approximations to the Jacobian operator that arises in modelling the evolution of a hydraulically driven fracture in a multi-layered elastic medium. The governing equations involve a highly non-linear coupled system of integro-partial differential equations along with the fracture front free boundary problem. We demonstrate that an incomplete LU factorization of these localized Jacobians yields an efficient preconditioner for the fully populated, stiff, non-symmetric system of algebraic equations that need to be solved multiple times for every growth increment of the fracture. The performance characteristics of this class of preconditioners is explored via spectral analysis and numerical experiment. |
|
Abstract: We present a novel multigrid (MG) procedure for the efficient solution of the large non-symmetric system of algebraic equations used to model the evolution of a hydraulically driven fracture in a multi-layered elastic medium. The governing equations involve a highly non-linear coupled system of integro-partial differential equations along with the fracture front free boundary problem. The conditioning of the algebraic equations typically degrades as O(N^3). A number of characteristics of this problem present significant new challenges for designing an effective MG strategy. Large changes in the coefficients of the PDE are dealt with by taking the appropriate harmonic averages of the discrete coefficients. Coarse level Greens functions for multiple elastic layers are constructed using a single dual mesh and superposition. Coarse grids that are sub-sets of the finest grid are used to treat mixed variable problems associated with pinch points. Localized approximations to the Jacobian at each MG level are used to devise efficient GaussSeidel smoothers and preferential line iterations are used to eliminate grid anisotropy caused by large aspect ratio elements. The performance of the MG preconditioner is demonstrated in a number of numerical experiments. |
|
Abstract: We present a planar three-dimensional (3D) fracture growth simulator, based on a displacement discontinuity (DD) method for multi-layer elasticity problems. The method uses a fixed mesh approach, with rectangular panel elements to represent the planar fracture surface. Special fracture tip logic is included that allows a tip element to be partially fractured in the tip region. The fracture perimeter is modelled in a piece-wise linear manner. The algorithm can model any number of interacting fractures that are restricted to lie on a single planar surface, located orthogonal to any number of parallel layers. The multiple layers are treated using a Fourier transform (FT) approach that provides a numerical Greens function for the DD scheme. The layers are assumed to be fully bonded together. Any fracture growth rule can be postulated for the algorithm. We demonstrate this approach on a number of test problems to verify its accuracy and efficiency, before showing some more general results. |
|
Abstract: We present a method that extends the flexibility matrix method for multilayer elasticity problems to include problems with very thin layers. This method is particularly important for solving problems in which one or a number of very thin layers are juxtaposed with very thick layers. The standard flexibility matrix method suffers from round-off errors and poor scaling of the flexibility equations which occur when one of the layers in the multilayered medium becomes much smaller than the others. The method proposed in this paper makes use of power series expansions of the various components of the flexibility matrix in order to arrive at a system of equations that is appropriately scaled. The effectiveness of the scaled flexibility matrix method is demonstrated on a number of test problems. |
|
Abstract:
We present uniform asymptotic solutions (UAS) for displacement discontinuities
(DD) that lie within the middle layer of a three layer elastic medium. The
DDs are assumed to be normal to the two parallel interfaces between the
leastic media, and solutions will be presented for both 2D and 3D elastic
media. Using the Fourier transform (FT) method we construct the leading
term in the asymptotic expansion for the spectral coefficient functions
for a DD in a three layer medium. Although a closed form solution will require
an infinite series solution, we demonstrate how this UAS can be used to
construct highly efficient and accurate solutions even in the case in which
the DD actually touches the interface. We present an explicit UAS for elements
in which the DD fields are assumed to be piecewise constant throughout a line segment in 2D and a rectangular element in 3D. We demonstrate the usefulness of this UAS by providing a number of examples in which the UAS is used to solve problems in which cracks just touch or cross an interface. The accuracy and efficiency of the algorithm is demonstrated and compared with other numerical methods such as the finite element method and the boudary integral method. |