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.
|
18. Coupling schemes for modeling hydraulic fracture propagation using the XFEM.
by E. Gordeliy and A. P. Peirce
accepted by Comp. Meth. in Appl. Mech. and Eng.
on 18 Aug.2012.
Top: The crack path obtained using a DD cde OribiC,
XFEM-t, and XFEM |
Abstract:
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. |
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 geometry—a 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 Perkins–Kern–Nordgren (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 Green’s 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 Gauss–Seidel 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 Green’s 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. |
