Variational Methods in the Simulation of Three-Dimensional Fracture Propagation.

Introduction The need to develop accurate predictive hydraulic fracturing simulators arises from the consideration of vertical containment for the fracture. Fracture height growth from a pay zone into adjacent layers has been shown (see Cleary [I]) to depend on three sets of parameters for the different strata: - elastic properties, - closure stress variations, and, -fracturing fluid leak off, and poro-elastic stresses. Qualitative analysis of the effects of these parameters has been reported in the literature. The first effects were described by Van Eekelen, using simple elasticity considerations [2]. The effect of stress variations was described numerically by Simonson and a1. [3] and experimentally by Thiercelin [4]. Effects of back stresses induced by fluid invasion on fracture propagation were considered by Cleary [5]. In order to couple these effects, a numerical code al- lowing the simulation of three dimensional growth of a hydraulically induced fracture has been developed with the ability of representing the effects of heterogeneities on fracture lateral and vertical extension. The model is based on a new variational formulation of the boundary integral representation of three-dimensional elasticity. The fluid flow and elasticity equations are solved using a moving mesh representation which offers enough flexibility to enable the simulation of any closure stress distribution. Simpler models representing fracture height growth (P3DH models) have been developed [6] [7] that approximate the fracture by a series of one dimensional segments, and the limitation of such approaches is examined in this paper. 1 Description of the model 1.1 Fracture Geometry Hydraulic fracturing simulations involve the coupling of solutions to elastic deformations, fluid flow in the fracture, and fluid leak off into the formation. Since the use of direct solution techniques (i,e,,3D finite elements) is prohibitive, the method used in this study is baaed on a boundary integral representation of the 3D elasticity equations. 1.2 3D Models The classical representation of the problem of a fracture in an infinite medium by boundary integrals is [8]: [ ] where u is the unknown displacement field, r n is the traction operator, I? is the fracture surface, and K is the Kelvin tensor of fundamental solutions to the elasticity operator. Using the notation: [ ] for the displacement discontinuity (or fracture width), one can attempt to solve the boundary integral relation using surface integrals (91 or the displacement discontinuity technique, first introduced by Crouch [lo], and extended to three-dimensional problems [Ill. However, numerical problems arise since the Kelvin tensor has a singularity in