US7509245B2 - Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator - Google Patents

Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator Download PDF

Info

Publication number
US7509245B2
US7509245B2 US11/095,104 US9510405A US7509245B2 US 7509245 B2 US7509245 B2 US 7509245B2 US 9510405 A US9510405 A US 9510405A US 7509245 B2 US7509245 B2 US 7509245B2
Authority
US
United States
Prior art keywords
layer
fracture
defines
stress
low
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related, expires
Application number
US11/095,104
Other versions
US20050171751A1 (en
Inventor
Eduard Siebrits
Anthony Peirce
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Schlumberger Technology Corp
Original Assignee
Schlumberger Technology Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from US09/301,961 external-priority patent/US6876959B1/en
Priority claimed from US10/831,799 external-priority patent/US7063147B2/en
Application filed by Schlumberger Technology Corp filed Critical Schlumberger Technology Corp
Priority to US11/095,104 priority Critical patent/US7509245B2/en
Assigned to SCHLUMBERGER TECHNOLOGY CORPORATION reassignment SCHLUMBERGER TECHNOLOGY CORPORATION ASSIGNMENT OF ASSIGNORS INTEREST (SEE DOCUMENT FOR DETAILS). Assignors: PEIRCE, ANTHONY, SIEBRITS, EDUARD
Publication of US20050171751A1 publication Critical patent/US20050171751A1/en
Priority to CA2541103A priority patent/CA2541103C/en
Priority to EA200600516A priority patent/EA010456B1/en
Priority to US12/368,549 priority patent/US8428923B2/en
Publication of US7509245B2 publication Critical patent/US7509245B2/en
Application granted granted Critical
Adjusted expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • EFIXED CONSTRUCTIONS
    • E21EARTH DRILLING; MINING
    • E21BEARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B43/00Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
    • E21B43/25Methods for stimulating production
    • E21B43/26Methods for stimulating production by forming crevices or fractures

Definitions

  • the subject matter set forth in this specification relates to a Hydraulic Fracturing Simulator Software adapted for designing and monitoring and evaluating reservoir, for example petroleum reservoir, fracturing, and, in particular, to a Hydraulic Fracturing Simulator Software adapted for simulating a multilayer reservoir and partially active elements of a mesh overlaying a fracture footprint when designing and monitoring and evaluating reservoir fracturing.
  • Proppant or propping agent is carried into the fracture by a viscosified fluid, and deposited into the fracture.
  • Proppant provides a permeable flow channel for formation fluids, such as oil and gas, to travel up the wellbore and to the Earth's surface.
  • Fracturing involves many variables, including: viscosity of the fracturing fluid, rate of leak-off of fracturing fluid into the reservoir, proppant carrying capacity of the fluid, viscosity of the fluid as a function of temperature, time history of fluid volumes (i.e., the amount of fluid pumped over a given period of time), time history of proppant volumes, fluid physical constants, proppant properties, and the geological properties of various zones in the reservoir.
  • a Hydraulic Fracturing Simulator software is capable of modeling the ‘fracturing event’.
  • the Hydraulic Fracturing Simulator software will design and monitor and evaluate the ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’ from a time extending before, during, and after the ‘fracturing event’.
  • the Hydraulic Fracturing Simulator software designs and monitors and evaluates the ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’
  • the Hydraulic Fracturing Simulator software should function to model a ‘multilayered reservoir’ when modeling the ‘fracturing event’.
  • the Hydraulic Fracturing Simulator software should also function to model ‘partially active elements’ (partially enclosed by the fracture footprint, as opposed to ‘fully active elements’ fully enclosed by the fracture footprint), among a plurality of elements of a mesh overlaying a fracture footprint, when the Hydraulic Fracturing Simulator software models the ‘fracturing event’.
  • a Hydraulic Fracturing Simulator Software is needed that is capable of modeling or simulating a ‘multilayered reservoir’ during a time when the Hydraulic Fracturing Simulator Software is designing and monitoring and evaluating ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’.
  • a Hydraulic Fracturing Simulator Software is needed that is capable of modeling or simulating ‘partially active elements’, among a plurality of elements of a mesh overlaying a fracture footprint, during a time when the Hydraulic Fracturing Simulator Software is designing and monitoring and evaluating ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’.
  • One aspect of the present invention involves a program storage device readable by a machine tangibly embodying a program of instructions executable by the machine to perform method steps of simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, the method steps comprising: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the multilayered reservoir.
  • Another aspect of the present invention involves a method of simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, comprising the step of: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the multilayered reservoir.
  • Another aspect of the present invention involves a program storage device readable by a machine tangibly embodying a program of instructions executable by the machine to perform method steps of simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, the method steps comprising: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
  • Another aspect of the present invention involves a method of simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, comprising the step of: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
  • a further aspect of the present invention involves a system adapted for simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, comprising: apparatus adapted for calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the system adapted for simulating the hydraulic fracture in the formation will take into account an existence of the multilayered reservoir.
  • a further aspect of the present invention involves a system adapted for simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, comprising: apparatus adapted for calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the system adapted for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
  • FIGS. 1 through 3 illustrate a typical Hydraulic Fracturing (HF) job in a wellbore
  • FIGS. 4 through 6 illustrate the fracture footprint created in the formation penetrated by the wellbore when the HF job is pumped
  • FIGS. 7 through 9 illustrate how a mesh comprised of a plurality of grid cells will overlay on top of the fracture footprint of FIGS. 4 through 6 , each grid cell of the mesh having a width and pressure, some of the grids cells called ‘tip elements’ being intersected by a perimeter of the fracture footprint, the tip elements having a width and a pressure (w, p), a portion of each tip element having fracturing fluid disposed therein;
  • FIG. 10 illustrates a mesh overlayed on top of a fracture footprint, one or more ‘fully active elements’ being enclosed by the fracture footprint, one or more ‘partially active elements’ being partially enclosed by the fracture footprint, the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘partially active elements’;
  • FIGS. 11 and 12 illustrate, in greater detail, the ‘partially active elements’ of FIG. 10 , the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘partially active elements’;
  • FIGS. 13 and 14 illustrate a ‘multi-layered reservoir’, the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘multi-layered reservoir’;
  • FIG. 15 illustrates an apparatus used in connection with a Hydraulic Fracturing (HF) job adapted for fracturing a formation penetrated by a wellbore, the apparatus including a computer system for storing a software called a ‘Hydraulic Fracturing Simulator software’;
  • HF Hydraulic Fracturing
  • FIG. 16 illustrates the computer system of FIG. 15 which stores the software called a ‘Hydraulic Fracturing Simulator software’ adapted for modeling a ‘multilayered reservoir’ and for modeling ‘partially active elements’ of a mesh overlaying a fracture footprint during a ‘petroleum reservoir fracturing’ event;
  • a ‘Hydraulic Fracturing Simulator software’ adapted for modeling a ‘multilayered reservoir’ and for modeling ‘partially active elements’ of a mesh overlaying a fracture footprint during a ‘petroleum reservoir fracturing’ event
  • FIG. 17 illustrates in greater detail the ‘Other Data’ of FIG. 16 ;
  • FIG. 18 illustrates the ‘Other Software Instructions’ of FIG. 16 ;
  • FIG. 19 illustrates a construction of the ‘Hydraulic Fracturing Simulator software’ of FIG. 16 ;
  • FIG. 20 illustrates the function associated with the ‘Set Up Influence Coefficient Matrix [C]’ step associated with the construction of the ‘Hydraulic Fracturing Simulator software’ which is illustrated in FIG. 19 ;
  • FIGS. 21 and 22 illustrate a more detailed construction of the ‘Hydraulic Fracturing Simulator software’ which is illustrated in FIG. 19 ;
  • FIGS. 23 and 24 illustrate a construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIGS. 19 and 21 ;
  • FIGS. 25-28 are used during a discussion of a detailed construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 in FIGS. 19 and 21 .
  • This specification discloses a Hydraulic Fracturing Simulator software adapted to be stored in a memory of a program storage device of a computer system for modeling and simulating a multilayered reservoir and for modeling and simulating partially active elements of a mesh overlaying a fracture footprint when the Hydraulic Fracturing Simulator software is designing and monitoring and evaluating petroleum reservoir fracturing.
  • the Hydraulic Fracturing Simulator software includes a first step further including the step of Setting Up an Influence Coefficient Matrix and a second step further including two iteration loops whereby, in a first iteration loop at a first time step, a second iteration loop will continuously calculate fracture width for each element of the mesh given a previously determined fluid pressure and will continuously calculate fluid pressure for each element of the mesh given a previously determined fracture width until convergence of the solution of width and pressure is reached at which time the fracture footprint is updated to a first value and ‘output data’ is generated, then, in the first iteration loop at a second time step, the second iteration loop calculates fracture width for each element of the mesh given a previously determined fluid pressure, and fluid pressure for each element of the mesh is calculated given a previously determined fracture width at which time the fracture footprint is updated to a second value and ‘output data’ is generated, and the process repeats, at which point, ‘additional output data’ is generated.
  • the ‘output data’ includes the previously determined first value of the fracture footprint and the previously determined second value of the fracture footprint.
  • the ‘output data’ representative of the updated fracture footprint at each of the incremented time steps, is recorded or displayed on a recorder or display device.
  • the Influence Coefficient Matrix is calculated and determined in a special way such that, when the Influence Coefficient Matrix is determined, and during any ‘fracturing event’, the Hydraulic Fracturing Simulator software will model and simulate the following: (1) a multilayered reservoir of the type illustrated in FIGS. 13 and 14 , and (2) partially active elements of a mesh which overlays a fracture footprint of the type illustrated in FIGS. 10 , 11 , and 12 .
  • a perforating gun 10 is disposed in a wellbore 12 and a packer 14 isolates a plurality of shaped charges 16 of the perforating gun 10 downhole in relation to the environment uphole.
  • the shaped charges 16 detonate and a corresponding plurality of perforations 18 are produced in a formation 20 penetrated by the wellbore 12 .
  • a fracturing fluid 22 is pumped downhole into the perforations 18 in accordance with a particular pumping schedule 24 .
  • the formation 20 surrounding the perforations 18 is fractured.
  • oil or other hydrocarbon deposits 26 begin to flow from the fractures, into the perforations 18 , into the wellbore 12 , and uphole to the surface.
  • the oil or other hydrocarbon deposits flow at a certain ‘production rate’ 28 (e.g., in m 3 /day).
  • FIG. 4 when the wellbore 12 of FIG. 2 is fractured, pump trucks 30 situated at the surface of the wellbore will pump fracturing fluid down a tubing and into the perforations 18 in the formation 20 penetrated by the wellbore, as shown in FIG. 2 .
  • the formation 20 includes different layers, such as the different layers 42 , one such layer being a perforated interval 40 .
  • a fracture footprint 32 will be formed in the perforated interval 40 (and possibly in additional adjacent intervals 42 ) of a formation 20 penetrated by the wellbore 12 .
  • the fracture footprint 32 will change to a new footprint 34 in the perforated interval 40 (and possibly in additional intervals 42 ) of a formation 20 penetrated by the wellbore 12 .
  • the fracture footprint 34 will change to a new footprint 36 in the perforated interval 40 (and possibly in additional intervals 42 ) of a formation 20 penetrated by the wellbore 12 .
  • the fracture footprint 36 will change to a new footprint 38 in the perforated interval 40 (and possibly in additional intervals 42 ) of a formation 20 penetrated by the wellbore 12 .
  • FIG. 5 a schematic three dimensional view of a fracture footprint, such as the fracture footprints 32 - 38 of FIG. 4 , is illustrated.
  • each fracture footprint 32 - 38 has a length ‘L’ and a height ‘H’ and a width ‘w’.
  • the wellbore 12 is illustrated again, and a plurality of perforations 18 are shown in the formation 20 penetrated by the wellbore 12 , as illustrated in FIGS. 1-3 .
  • the formation 20 includes a plurality of different layers 42 .
  • FIG. 6 when the pump trucks 30 of FIG.
  • a ‘fracture footprint’ 46 is created in the formation 20 which is identical to the fracture footprints 32 , 34 , 36 , and 38 shown in FIG. 4 that are created, respectively, over the different periods of time t 1 , t 2 , t 3 , and t 4 .
  • the ‘fracture footprint’ 46 in FIG. 6 has a cross section 44 , the cross section 44 having a fracture width ‘w’ similar to the width ‘w’ of the fracture footprint 32 - 38 shown in FIG. 5 .
  • FIG. 7 recalling the fracture footprint 46 shown in FIG. 6 , in FIG. 7 , a mesh 48 comprised of a plurality of grid elements 48 a or grid cells 48 a is illustrated.
  • the fracture footprint 46 is assumed (by the model of this specification) to lie in the 2D plane, although, in principle and in reality, the fracture footprint 46 does lie in the 3D plane.
  • the mesh 48 is overlayed over the top of the fracture footprint 46 of FIG. 6 .
  • the mesh 48 includes a plurality of active elements or active grid cells 48 a 1 and a plurality of inactive elements or inactive grid cells 48 a 2 .
  • each of the active grid cells 48 a 1 will overlay the fracture footprint 46 , whereas the inactive grid cells 48 a 2 will not overlay the fracture footprint 46 .
  • Each of the active grid elements or grid cells 48 a 1 of the mesh 48 has a width ‘w’ and a pressure ‘p’ assigned thereto, denoted by the symbol: (w, p). Therefore, each active grid cell 48 a 1 of the mesh 48 has a width/pressure value (w, p) assigned to that active grid cell.
  • FIG. 6 since the fracturing fluid propagating down the wellbore 12 enters the perforations 18 and creates the fracture footprint 46 , in FIG. 7 , each of the active grid cells 48 a 1 in the mesh 48 has a fracturing fluid disposed therein.
  • FIG. 6 since the fracturing fluid propagating down the wellbore 12 enters the perforations 18 and creates the fracture footprint 46 , in FIG. 7 , each of the active grid cells 48 a 1 in the mesh 48 has a fracturing
  • active grid cells 48 a 1 there are two types of active grid cells 48 a 1 : (1) an active grid cell 48 a 1 which is intersected by a perimeter 46 a of the fracture footprint 46 , and (2) an active grid cell 48 a 1 which is not intersected by the perimeter 46 a of the fracture footprint 46 .
  • An active grid cell 48 a 1 that is intersected by the perimeter 46 a of the fracture footprint 46 is known as a ‘tip element’.
  • ‘tip element’ 50 is an active grid cell 48 a 1 that is intersected by the perimeter 46 a of the fracture footprint 46 .
  • An active grid cell 48 a 1 which is not intersected by the perimeter 46 a of the fracture footprint 46 has a volume which is wholly occupied by the fracturing fluid (i.e., 100% of the volume of the active grid cell is occupied by the fracturing fluid), where the fracturing ‘fluid’ may or may not include proppant.
  • active grid cell 52 is not intersected by the perimeter 46 a of the fracture footprint 46 and its volume is wholly occupied by the fracturing fluid (100% of the volume of the active grid cell 52 is occupied by fracturing fluid).
  • an active grid cell 48 a 1 that is intersected by the perimeter 46 a has a volume that is occupied by ‘less than 100%’ of the fracturing fluid.
  • the active grid cell or ‘tip element’ 54 is intersected by the perimeter 46 a of the fracture footprint 46 , however, only 45% of the volume of the active grid cell 54 is occupied by the fracturing fluid.
  • an inactive grid cell 48 a 2 such as inactive grid cell 55 , has a volume which is wholly devoid of any fracturing fluid (0% of the volume of the inactive grid cell 55 is occupied by fracturing fluid).
  • the ‘VOF software’ of that prior pending application calculates, over a series of time steps, the ‘amount of fracturing fluid that is contained within each of the active grid cells 48 a 1 that are intersected by the perimeter 46 a of the fracture footprint 46 ’. That is, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 now U.S. Pat. No. 7,063,147 calculates, over the series of time steps, the ‘amount of fracturing fluid that is contained within each of the tip elements 50 ’.
  • the ‘amount of fracturing fluid that is contained within each of the tip elements 50 ’ is calculated from the ‘fill fraction’, the ‘fill fraction’ being denoted by the letter ‘F’.
  • the ‘fill fraction’ F for the ‘tip element’ 54 is 45%. Therefore, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 calculates, over a series of time steps, the ‘fill fraction’ (F) for each of the ‘active grid cells 48 a 1 of the mesh 48 that are intersected by the perimeter 46 a of the fracture footprint 46 ’. That is, the ‘VOF software’ of prior pending application Ser. No.
  • 10/831,799 calculates, over the series of time steps, the ‘fill fraction’ (F) for each of the ‘tip elements’ 50 of the mesh 48 of FIG. 7 .
  • the amount or degree by which the perimeter 46 a of the footprint 46 of the fracture is expanding (or contracting), as a result of the pumping of the fracturing fluid into the perforations 18 in the formation 20 by the pump trucks 30 of FIG. 4 can be determined.
  • the ‘VOF software’ of prior pending application Ser. No. 10/831,799 now U.S. Pat. No. 7,063,147 is embodied in steps 138 and 96 of FIG. 22 .
  • the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21 (and the ‘Set Up Influence Coefficient Matrix [C]’ 102 of FIG. 19 ) is discussed in detail.
  • FIGS. 8 and 9 two more examples of a mesh 48 similar to the mesh 48 of FIG. 7 are illustrated.
  • a mesh 48 is illustrated as overlaying the fracture footprint 46 having a perimeter 46 a . Fracturing fluid is disposed inside the perimeter 46 a , but the fracturing fluid is not disposed outside the perimeter 46 a .
  • the active grid cell 56 let us analyze the active grid cell 56 .
  • the ‘VOF software’ of prior pending application Ser. No. 10/831,799 also calculates the volume of an active grid cell occupied by a first type of fluid and the volume of that same active grid cell occupied by a second type of fluid.
  • active grid cell 58 includes a first volume of 80% occupied by a first type of fluid, and a second volume of 20% occupied by a second type of fluid.
  • 10/831,799 calculates, over a series of time steps, the ‘fill fraction’ (F) for each of the ‘active grid cells 48 a 1 that are intersected by the perimeter 46 a of the fracture footprint 46 ’ in the mesh 48 ; that is, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 will calculate, over the series of time steps, the ‘fill fraction’ (F) for each of the ‘tip elements’ in the mesh 48 shown in FIGS. 7 , 8 , and 9 .
  • the amount or degree by which the perimeter 46 a of the footprint 46 is expanding (or contracting), in response to the pumping of fracturing fluid into the perforations 18 can be determined.
  • FIG. 10 ‘fully active elements’ and ‘inactive elements’ and ‘partially active elements’ are illustrated.
  • a wellbore 12 is illustrated, and a fracture 17 grows in a direction that is directed away from the wellbore 12 .
  • a mesh 19 overlays the fracture 17 .
  • an ‘inactive element’ 21 is illustrated, and a ‘fully active element’ 23 is also illustrated.
  • FIG. 10 relates to ‘fracture growth’ of the fracture 17 where the ‘growth’ of the fracture is in a direction that is either directed away from the wellbore 12 (i.e., expansion) or is directed toward the wellbore 12 (contraction).
  • the fracture 17 represents a moving boundary that is growing in a direction which is either directed away from the wellbore 12 (expansion) or is directed toward the wellbore 12 (contraction).
  • ‘partially active elements’ 25 and 27 and 29 and 43 are illustrated.
  • the ‘partially active element’ 25 and 27 only a portion 25 a and 27 a of the ‘partially active element’ 25 and 27 is disposed inside of the boundary 17 of the fracture 17 . Since only portions 25 a and 27 a of ‘partially active elements’ 25 and 27 are disposed inside the boundary 17 of fracture 17 , the elements 25 and 27 of mesh 19 of FIG. 10 are each known as a ‘partially active element’.
  • FIG. 11 illustrates a close up view of the ‘partially active element’ 29 of FIG. 10 .
  • the ‘partially active element’ 29 has a fracture leading edge 31 with crossing points 33 and 35 , respectively.
  • Straight line 37 is erected between the crossing points, and it forms the boundary for the active portion 39 of the element 29 and the inactive portion 41 of the element 29 .
  • FIG. 12 illustrates the same features for the ‘partially active element’ 43 of FIG. 10 as did FIG. 11 with respect to the ‘partially active element’ 29 of FIG. 10 .
  • a ‘partially active element’ 43 has a fracture leading edge 31 with crossing points 33 and 35 , respectively.
  • Straight line 37 is erected between the crossing points 33 , 35 , and it forms the boundary for the active portion 39 of the element 43 and the inactive portion 41 of the element 43 .
  • FIG. 13 a reservoir or Earth formation 20 is shown.
  • pumping truck 47 provides fluid at high pressures and flow rates to wellhead 49 , which is operably connected to the wellbore 12 at or near the ground surface 53 .
  • FIG. 13 illustrates the fracture boundary 55 at a particular time. Two other fracture fluid boundaries 57 and 59 also are indicated in FIG. 13 .
  • the reservoir 20 represents a ‘multi-layered reservoir’ 20 because the reservoir 20 in FIG.
  • the fluid boundaries in FIG. 13 reveal separate types or compositions of pumped fluid.
  • the fracture preferably is stopped prior to the water bearing zone 71 seen at the lower portion of FIG. 13 .
  • a ‘multilayered’, hydraulically fractured reservoir 20 or Earth formation 20 is penetrated by a wellbore 12 .
  • the wellbore 12 penetrates the Earth formation 20 and the hydraulic fracture 59 of FIG. 13 is shown between layers 20 f and 20 h .
  • the fracture 59 has a fracture width 59 a which was created in response to the pressurization of the formation 20 by a fracturing fluid.
  • the ‘multilayered’ reservoir/Earth formation 20 includes a plurality of Earth formation layers 20 a , 20 b , 20 c , 20 d , 20 e , 20 f , 20 g , 20 h , 20 i , and 20 j .
  • Each of the layers 20 a - 20 j can be characterized by a Young's Modulus (E) and a Poisson's Ratio ( ⁇ ).
  • the Young's Modulus (E) and the Poisson's Ratio ( ⁇ ) describe the elastic properties of each of the layers 20 a - 20 j in FIG. 14 .
  • the elastic properties of layer 20 a can be characterized by Young's Modulus and Poisson's Ratio (E a , ⁇ a ), the elastic properties of layer 20 b being characterized by Young's Modulus and Poisson's Ratio (E b , ⁇ b ), the elastic properties of layer 20 c being characterized by Young's Modulus and Poisson's Ratio (E c , ⁇ c ), the elastic properties of layer 20 d being characterized by Young's Modulus and Poisson's Ratio (E d , ⁇ d ), the elastic properties of layer 20 e being characterized by Young's Modulus and Poisson's Ratio (E e , ⁇ e ), the elastic properties of layer 20 f being characterized by Young's Modulus and Poisson's Ratio (E f , ⁇ f ), the elastic properties of layer 20 g being characterized by Young's Modulus and Poisson's Rati
  • the pump trucks 30 of FIG. 4 will pump a fracturing fluid 62 (usually including proppant) down the wellbore 12 of FIG. 4 in accordance with a pumping schedule 60 (an example used in connection with this discussion).
  • the fracturing fluid 62 will enter the perforations 18 , and, responsive thereto, create a ‘fracture footprint’ 46 , similar to the fracture footprint 46 shown in FIG. 6 .
  • Micro-seismic data sensor(s) 64 and tiltmeter data or other sensor(s) 66 will monitor the approximate geometry of the fracture footprint 46 and, responsive thereto, the sensor(s) 64 and 66 will generate output signals, the micro-seismic data sensor(s) 64 generating a micro-seismic data output signal 64 a , the tiltmeter data sensor(s) 66 generating a tiltmeter data output signal 66 a , and the pumping schedule 60 generating a pumping schedule output signal 60 a representative of the pumping schedule 60 .
  • the pumping schedule output signal 60 a , the tiltmeter data output signal 66 a , and the micro-seismic data output signal 64 a are time line merged via a ‘time line merging’ step 68 .
  • the pumping schedule output signal 60 a , the tiltmeter data output signal 66 a , and the micro-seismic data output signal 64 a are ‘time synchronized’ in a particular manner such that the tiltmeter data output signal(s) 66 a and the micro-seismic data output signal(s) 64 a are synchronized with the times present in the pumping schedule 60 .
  • a ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 is generated which is provided as ‘input data’ to a ‘computer system’ 72 optionally disposed within a truck 74 situated at the site of the wellbore 12 , such as a monitoring truck 74 or a ‘FracCAT vehicle’ 74 (which is a vehicle with software that monitors and measures the fracture and controls the fracture treatment).
  • a ‘computer system’ 72 optionally disposed within a truck 74 situated at the site of the wellbore 12 , such as a monitoring truck 74 or a ‘FracCAT vehicle’ 74 (which is a vehicle with software that monitors and measures the fracture and controls the fracture treatment).
  • the ‘computer system’ 72 which is optionally disposed within the truck 74 of FIG. 15 , such as the ‘FracCAT vehicle’ 74 , is illustrated.
  • the ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 of FIG. 15 is provided as ‘input data’ to the computer system 72 disposed within the truck 74 , the output signal 70 being comprised of a time line merged pumping schedule and tiltmeter data and micro-seismic data plus other data including downhole temperature and bottom hole pressure.
  • the 16 includes a processor 72 a operatively connected to a system bus, a memory or other program storage device 72 b operatively connected to the system bus, and a recorder or display device 72 c operatively connected to the system bus.
  • the memory or other program storage device 72 b stores the following software ( 76 , 78 , and 80 ): a Planar 3D User Interface software 76 , a Planar 3D ‘engine’ or software 78 , and a Hydraulic Fracturing Simulator Software 80 .
  • the software 76 , 78 , and 80 which is stored in the memory 72 b of FIG. 16 , can be initially stored on a CD-ROM, where that CD-ROM is also a ‘program storage device’.
  • the computer system 72 of FIG. 16 receives Input Data 82 , including: (1) the time line merged pumping schedule, tiltmeter data, and micro-seismic data 84 (which corresponds to the ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 of FIG. 15 ), and (2) Other Data 86 .
  • the computer system 72 also receives Other Software Instructions 88 .
  • the processor 72 a will execute the Hydraulic Fracturing Simulator Software 80 (including the Planar 3D User Interface software 76 and the Planar 3D ‘engine’ 78 ), while simultaneously using the Input Data 82 and Other Software Instructions 88 ; and, responsive thereto, the recorder or display device 72 c will generate a set of ‘Output Data’ 72 c 1 which is adapted to be recorded by or displayed on the recorder or display device 72 c .
  • the computer system 72 may be a personal computer (PC), a workstation, or a mainframe.
  • the memory or program storage device 72 b is a computer readable medium or a program storage device which is readable by a machine, such as the processor 72 a .
  • the processor 72 a may be, for example, a microprocessor, microcontroller, or a mainframe or workstation processor.
  • the memory or program storage device 72 b which stores the Hydraulic Fracturing Simulator Software 80 , may be, for example, a hard disk, ROM, CD-ROM, DRAM, or other RAM, flash memory, magnetic storage, optical storage, registers, or other volatile and/or non-volatile memory.
  • the Other Data 86 and the Other Software Instructions 88 of FIG. 16 are illustrated.
  • the Other Data 86 will include: layer confining stresses and properties, perforating interval and depth, wellbore data, fluid and proppant properties, time history of fluid volumes to be pumped, time history of proppant volumes to be pumped, and logs identifying the identity, properties, and location of geological zones.
  • the Other Software Instructions 88 include instructions to calculate values representing physical dimensions of the fracture and pressures inside the fracture.
  • the Hydraulic Fracturing Simulator Software 80 includes: an initial ‘time stepping’ step 90 , a second step 92 which calculates fracture width (w) given the fluid pressure (p); a third step 94 which calculates fluid pressure (p) given the fracture width (w); and a fourth step 96 adapted to update the fracture footprint.
  • the fracture ‘footprint’ is shown in FIGS. 4 and 13 at times t 1 , t 2 , t 3 , . . . , t n (for example at times t 1 , t 2 , t 3 , and t 4 in FIG. 4 and at times corresponding to the three boundaries 59 , 57 and 55 in FIG. 13 ).
  • the time stepping loop 98 in FIG. 19 wherein the ‘time stepping’ step 90 will increment from a first time step (t 1 ), to a second time step (t 2 ), to a third time step (t 3 ), . . . , and to an nth time step (t n ).
  • the second step 92 and the third step 94 will iterate on each other (as indicated by arrow 100 ) until the fracture width (w) and the fluid pressure (p) are solved at every point (i.e., within each ‘element’ or ‘grid cell’) in the fracture geometry of FIGS. 4 and 13 .
  • the fracture footprint is updated in step 96 .
  • the proppant concentration is calculated for each grid cell in the updated fracture footprint.
  • the ‘output data’ 72 c 1 associated with the last iteration loop 100 , is generated and stored.
  • the time step 90 is then incremented to the second time step t 2 , and, responsive thereto, the second step 92 and the third step 94 will then iterate again on each other (as indicated by arrow 100 ) until the fracture width (w) and the fluid pressure (p) are solved at every grid cell in the fracture geometry of FIGS. 4 and 13 .
  • the fracture footprint is updated in step 96 .
  • the proppant concentration is calculated for each grid cell in the updated fracture footprint.
  • the ‘output data’ 72 c 1 associated with the last iteration loop 100 , is generated and stored.
  • the time step 90 is then incremented once again to the third time step t 3 , the second step 92 and the third step 94 will then iterate again on each other (as indicated by arrow 100 ) until the fracture width (w) and the fluid pressure (p) are solved at every point (i.e., within each ‘element’ or ‘grid cell’) in the fracture geometry of FIGS. 4 and 13 .
  • the fracture footprint is updated in step 96 .
  • step 97 the proppant concentration is calculated for each grid cell in the updated fracture footprint.
  • the ‘output data’ 72 c 1 associated with the last iteration loop 100 , is generated and stored.
  • the time step 90 is then incremented once again, and the above process repeats until it has reached the end of the pumping schedule.
  • ‘Convergence’ takes place when the ‘solution does not change from one iteration to the next’.
  • the inner iteration loop 100 is solving for two things: (1) fracture width (w) using the ‘elasticity equation’, and (2) fluid pressure (p) using the ‘fluid flow equation’.
  • Steps 92 and 94 can be solved in various ways, such as by iteration of two equations as shown here, or by direct substitution of the one equation into the other equation, or vice-versa.
  • the Output Data 72 c 1 is generated at the end of each time step. However, ‘additional output data’ 99 is generated when the time stepping in the outer iteration loop 98 is complete, the ‘additional output data’ 99 being used in subsequent calculations, such as in the generation of other graphical plots.
  • a first step 102 is practiced, the first step 102 being called ‘Set Up Influence Coefficient Matrix [C]’ 102 .
  • the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 will be discussed below with reference to FIG. 20 of the drawings.
  • step 92 will calculate the fracture width (w) and step 94 will calculate the fluid pressure (p).
  • the fracture width (w) of step 92 is actually calculated by using the ‘Elasticity Equation’ 104 .
  • FIG. 20 a generic form of the ‘Influence Coefficient Matrix [C]’ 108 is shown, where the ‘Influence Coefficient Matrix [C]’ 108 has ‘M’ rows and ‘M’ columns.
  • the ‘Influence Coefficient Matrix [C]’ 108 is fully populated with numbers 114 and is used to calculate fracture width (w) for step 92 in FIG. 19 .
  • the ‘elasticity equation’ 104 did not previously take into account the existence of a ‘multilayered reservoir’ of the type illustrated in FIGS. 13 and 14 , and the ‘elasticity equation’ 104 did not previously take into account the existence of ‘partially active elements’, such as the ‘partially active elements’ 25 , 27 , 29 , and 43 of FIG. 10 .
  • ‘partially active elements’ such as the ‘partially active elements’ 25 , 27 , 29 , and 43 of FIG. 10 .
  • a ‘multilayered reservoir’ does exist.
  • ‘partially active elements’ also do exist.
  • the Hydraulic Fracturing Simulator Software 80 of FIG. 16 disclosed in this specification will include the effect of, and take into account the existence of, ‘multilayered reservoirs’ and ‘partially active elements’. Consequently, in order to take into account the existence of ‘multilayered reservoirs’ and ‘partially active elements’, it is necessary to execute the ‘set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 , associated with the Hydraulic Fracturing Simulator Software 80 of FIGS. 16 and 19 .
  • the ‘Set Up Influence Coefficient Matrix [C]’ of Step 102 of FIG. 19 will calculate the ‘Influence Coefficient Matrix [C]’ 108 in the ‘Elasticity Equation’ 104 of FIG.
  • the numbers 114 in the ‘Influence Coefficient Matrix [C]’ 108 will be changed accordingly; and, when the numbers 114 of the ‘Influence Coefficient Matrix [C]’ 108 are changed accordingly, the Hydraulic Fracturing Simulator Software 80 will then include the effect of, or take into account the existence of, the ‘multilayered reservoirs’ and the ‘partially active elements’.
  • the Hydraulic Fracturing Simulator Software 80 will then include the effect of, or take into account the existence of, the ‘multilayered reservoirs’ and the ‘partially active elements’.
  • FIGS. 21 and 22 a more detailed construction of the Hydraulic Fracturing Simulator Software 80 of FIG. 19 , which is disposed within the Planar 3D software 78 and the Planar 3D User Interface 76 , is illustrated.
  • the input data 82 is provided from the human interface, such as the injection rate and the pumping schedule, which includes the injection rate as a function of time, proppant concentration, fluid viscosity, the geology or the properties of the elastic layers of the reservoir including the elastic constants comprising the Young's Modulus and the Poisson's Ratio (E, ⁇ ) as previously described, and the leakoff behavior.
  • the depths of each of the ‘interfaces’ of the ‘layers’ of FIG. 13 are calculated, such as the end of the sandstone layer, the end of the gas layer, etc.
  • the Young's Modulus, Poisson's Ratio, Toughness, Leakoff Coefficients, Stresses are assigned to each of the ‘layers’ (of FIG. 13 ) the depths of which were calculated in step 116 . Therefore, a series of numbers (comprising the Young's Modulus, Poisson's Ratio, Toughness, Leakoff Coefficients, and Stresses) is assigned to each Earth formation ‘layer’ that is shown in FIG. 13 .
  • the maximum possible ‘length’ to which the ‘fracture’ will propagate and the maximum possible ‘height’ to which the ‘fracture’ will propagate is assumed or introduced.
  • a parent mesh is assigned to the ‘fracture surface’, where the parent mesh is divided into rectangular ‘elements’. See FIGS. 7 and 10 for examples of the parent mesh.
  • the parent mesh is broken down into rows and columns comprising ‘grids’ or ‘elements’ which are generally rectangular in shape. It is assumed that the ‘fracture’ will grow into the ‘elements’ of the parent mesh, but no further.
  • the numerical parent mesh is generated; that is, in step 122 , the ‘coordinates’ of each ‘grid cell’ or ‘element’ of the parent mesh are generated.
  • the ‘coordinates’ will define where each ‘grid cell’ or ‘element’ exists within an axis system, such as the (x, y) axis system. Therefore, in step 122 of FIG. 21 , the ‘coordinates’ of each of the ‘grid cells’ (such as ‘grid cells’ 23 , 25 , 29 , and 43 shown in FIG. 10 or each of the ‘grid cells’ 48 a 1 and 52 shown in FIG. 7 ) are allocated.
  • the next step is the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 which is the ‘Set up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 .
  • the Influence Coefficient Matrix [C] ( 108 of FIG. 20 ) is generated.
  • the Influence Coefficient Matrix [C] (as shown in FIG.
  • each of the ‘numbers’ 114 of the Influence Coefficient Matrix [C] will indicate how one ‘element’ of the mesh relates to another ‘element’ of the mesh in an ‘elastic manner’. For example, if a fracture exists within one ‘element’ of the mesh, what kind of stress or strain exists within another ‘element’ of the mesh in response thereto (a phenomenon that is called an ‘elasticity behavior’). Each of the ‘numbers’ 114 in the Influence Coefficient Matrix [C] describes this ‘elasticity behavior’.
  • the Influence Coefficient Matrix [C] ( 108 of FIG. 20 ) is a matrix that indicates how one ‘element’ of the parent mesh (such as the mesh shown in FIGS. 7 and 10 ) talks to another ‘element’ of the parent mesh in an ‘elastic manner’. Therefore, step 102 of FIGS. 19 and 21 (i.e., the ‘Set Up Influence Coefficient Matrix’ step 102 of FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21 ), which determines the Influence Coefficient Matrix 108 of FIG. 20 , is responsible for allowing the Hydraulic Fracturing Simulator Software 80 of FIG. 16 to simulate the ‘multilayered reservoir’ of FIG.
  • Step 124 of FIG. 21 entitled the ‘Restart Option’ will allow the user to restart the simulator without having to recalculate the Influence Coefficient Matrix (it is stored on the hard drive). Such situations will occur if the reservoir properties are unchanged, but the user wants to change the injection schedule, for example.
  • the loop 98 which represents a time stepping loop. The loop 98 controls how the fracture grows as a function of time during the duration of the ‘pumping schedule’, where the term ‘pumping schedule’ includes the ‘shut-in time’.
  • step 126 at any particular time step, the fracture will have propagated to some extent in the parent mesh, which is called the ‘footprint’ of the fracture; step 126 will sweep through all the ‘grid cells’ in the entire parent mesh, and then step 126 will determine if a particular ‘element’ or ‘grid cell’, within the perimeter of the ‘fracture footprint’, is an ‘inactive element’, or an ‘active element’, or a ‘partially active element’.
  • step 128 entitled ‘Extract Elastic Influence Coefficient Matrix for Current Fracture Geometry’ the Influence Coefficient Matrix [C] 108 , as previously described with reference to step 102 , contains ‘all possible combinations within the parent mesh’; and, at any particular time step, we need a subset of the aforementioned ‘all possible combinations within the parent mesh’; thus, step 128 extracts, from the Influence Coefficient Matrix [C] 108 , a ‘submatrix’ which contains essential information regarding the current size of the ‘fracture footprint’.
  • step 130 entitled ‘Assign Special Characteristics to Coefficients for Tip Elements’ a ‘particular subset’ of the ‘grid cells’ or ‘elements’ of the parent mesh will cross a ‘perimeter of the fracture footprint’; in connection with only the ‘particular subset’ of the ‘grid cells’ which cross the ‘perimeter of the fracture footprint’, it is necessary to change some of the properties of the ‘submatrix’ (and recall that the ‘submatrix’ contains essential information regarding the current size of the ‘fracture footprint’); step 130 will make this change to the ‘submatrix’; consequently, step 130 in FIG.
  • step 130 in FIG. 22 will change the ‘partially active element submatrix’ (i.e., a submatrix which pertains specifically to the ‘partially active elements’) by ‘assigning special characteristics to the coefficients’ in the ‘partially active element submatrix’.
  • step 132 entitled ‘Assign Loading (Fluid Pressure) to Each Element in Current Mesh’ we now have a current time step and a current fracture footprint; in addition, each ‘active element’ has a ‘fluid pressure’; thus, step 132 will assign a ‘fluid pressure’ to each ‘element’ of the parent mesh.
  • step 92 entitled ‘Solve Elastic Equation System for Fracture Widths’ we use the ‘elasticity equation’ 104 of FIG. 20 (including their Fourier Transforms) to solve for the ‘fracture width’ in each ‘element’ or ‘grid cell’ of the mesh, at the current time step, given the ‘fluid pressure’ that was assigned to each ‘element’ of the mesh in step 132 .
  • step 134 entitled ‘Calculate Influence Matrix for Fluid Flow in Current Mesh Using Widths from Above’
  • step 94 entitled ‘Solve Fluid Flow Equations for Fluid Pressure in Current Mesh’
  • step 94 when the ‘fluid pressures’ are determined in step 94 , those ‘fluid pressures’ determined from step 94 are used to solve for ‘fracture widths’ in step 92 , and that loop 100 will continue to process as indicated until ‘global mass balance’ is achieved in step 136 .
  • step 136 to step 132 an inner iteration will take place; in that inner iteration, we continue to iterate on fluid pressure and fracture width until convergence of the solution takes place at that time step; convergence takes place when the solution does not change from one iteration to the next.
  • the inner iteration loop 100 is solving for two things: (1) fracture width using the ‘elasticity equation’, and (2) the fluid pressure using the ‘fluid flow equation’; when the change in each of these ‘equations’ is below a tolerance, we know that the inner iteration 100 has converged.
  • Steps 138 and 96 involve updating for the next time step.
  • Step 138 (‘Calculate Local Fracture Tip Velocity of Propagation’)
  • step 96 (‘Grow Fracture: Update New Layout’) are described in prior pending application Ser. No. 10/831,799, filed Apr. 26, 2004, now U.S. Pat. No.
  • step 97 the proppant concentration is calculated for each grid cell in the updated fracture footprint.
  • the ‘output data’ 72 c 1 associated with the last iteration loop 100 , is generated and stored.
  • step 140 if the time is not greater than the maximum allowed, time step 90 is incremented once again, and the above process repeats until it has reached the end of the pumping schedule.
  • step 140 if the ‘time’ is greater than the maximum allowed, ‘additional output data’ 99 is generated.
  • the ‘output data’ 72 c 1 and the ‘additional output data’ 99 is stored in a storage device, such as the memory or program storage device 72 b of FIG. 16 .
  • the Output Data 72 c 1 includes, at each time step, the fluid pressure (p) and the fracture width (w) at each of the ‘elements’ or ‘grid cells’ 23 , 25 , 29 , and 43 shown in FIG. 10 or each of the ‘elements’ or ‘grid cells’ 48 a 1 and 52 shown in FIG. 7 (again, at each time step). Since the Output Data 72 c 1 includes fluid pressure (p) and fracture width (w) for each ‘grid cell’ of FIGS.
  • the Output Data 72 c 1 can also include a number of 2D or 3D plots representative of the fluid pressure (p) and the fracture width (w) at each of the ‘elements’ or ‘grid cells’ 23 , 25 , 29 , and 43 shown in FIG. 10 or each of the ‘elements’ or ‘grid cells’ 48 a 1 and 52 shown in FIG. 7 at each time step.
  • Steps 92 , 134 , and 94 can be solved in various ways, such as by iteration of the two equations shown here (i.e., the ‘elasticity equation’ and the ‘fluid flow equation’), or by direct substitution of one equation into the other equation, or vice-versa.
  • FIGS. 23 and 24 a construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21 , which would allow the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’ (of the type illustrated in FIG. 13 ), is illustrated.
  • FIGS. 23 and 24 recall that the ‘Set Up Influence Coefficient Matrix [C]’ step 102 in FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG.
  • FIGS. 23 and 24 a ‘construction of the Set Up Influence Coefficient Matrix step 102 ’ of FIGS. 19 and 21 , for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’ (of the type illustrated in FIG. 13 ), is illustrated.
  • the ‘Hydraulic Fracturing Simulator software’ 80 for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 to model or take into account the existence of a ‘multi-layer reservoir, includes a ‘plurality of calculation steps’ 102 , the ‘plurality of calculation steps’ 102 being practiced by the ‘Set Up Influence Coefficient Matrix’ step 102 in FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 when the Influence Coefficient Matrix 108 of FIG. 20 is calculated.
  • the ‘construction of the Set Up Influence Coefficient Matrix step 102 ’ of FIGS. 19 and 21 for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, includes ‘three basic steps’ 102 a , 102 b , and 102 c , as follows: (1) Find Spectral Coefficients in Each Layer, step 102 a , (2) Perform Exponential Approximation, step 102 b , and (3) Assemble an ‘Influence Coefficient Matrix’ (i.e., matrix 108 of FIG. 20 ) Using Exponential Expansion Coefficients, step 102 c .
  • the Influence Coefficient Matrix is sent to a storage device.
  • FIG. 24 a detailed construction of each of the ‘three basic steps’ 102 a , 102 b , and 102 c of FIG. 23 will be discussed.
  • the ‘plurality of calculation steps’ 102 which are practiced by the ‘Set Up Influence Coefficient Matrix’ step 102 in FIG. 19 and by the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, include the following steps:
  • the ‘Find Spectral Coefficients in Each Layer’ step 102 a of FIG. 23 includes the following step (1), as follows:
  • the ‘Perform Exponential Approximation’ step 102 b of FIG. 23 includes the following steps (2), (3), and (4), as follows:
  • the ‘Assemble Influence Coefficient Matrix Using Exponential Expansion Coefficients’ step 102 c of FIG. 23 includes the following steps (5), (6), and (7), as follows:
  • step 102 of FIG. 19 is set forth in the following paragraphs with reference to FIGS. 25 through 28 .
  • step 102 of FIG. 19 is the same step as step 102 in FIG. 21 .
  • the numerical algorithm employed in this invention comprises an efficient technique to determine the local width of a hydraulic fracture due to local pressure applied to the fracture faces by the injection of hydraulic fluid and proppant into the fracture. Further, a method to track the dimensions and width of said fracture as it grows as a function of time is shown.
  • the hydraulic fracture(s) may span any number of layers in a laminated reservoir, with the restriction that all layers must be parallel to each other, as depicted for example in FIG. 23 . Layers may be inclined at any angle to the horizontal.
  • FIG. 25 shows a section through multiple hydraulic fractures in a layered reservoir.
  • the calculation of the fracture width due to the pressure from the injected fluids and proppant mixture is determined by taking into account, accurately and efficiently, the physical properties of each layer comprising the laminated reservoir.
  • the technique used to calculate the relationship between the layered reservoir and the growing hydraulic fracture is based on a well-established numerical technique called the Displacement Discontinuity Boundary Element (hereafter “DD”) Method.
  • DD Displacement Discontinuity Boundary Element
  • the method is extended to enable efficient and accurate calculation of the physical effects of layering in the reservoir by the use of a Fourier Transform Method, whereby the relations between stress and strain in the layered reservoir are determined.
  • the numerical method assumes that each hydraulic fracture is divided into a regular mesh of rectangular elements, as depicted in FIG. 26 , wherein each numerical element contains its own unique properties.
  • Such properties include applied fluid and proppant pressure, fluid and proppant propagation direction and velocity, local reservoir properties,
  • FIG. 26 shows a numerical mesh of elements subdividing the fracture surface for purposes of calculation.
  • certain additional information is uniquely defined for such elements.
  • such information may include the local velocity of propagation of the fracture tip, the special relationship between the fluid in the fracture and the surrounding layered reservoir, and how the fluid and reservoir physically interact with each other. This interaction is accounted for by means of special properties assigned to the tip elements, comprising the interaction between a fluid-filled fracture and the host material it is fracturing.
  • FIG. 27 reveals a fracture outline on a numerical mesh.
  • Each numerical element depicted in FIG. 26 or 27 relates to every other element in the numerical mesh by means of special mathematical relationships.
  • elements (1) sending or source elements, and (2) receiver elements.
  • a source element sends a signal representing a mathematical relationship to a receiver element.
  • the signal is the net pressure, which is equal to the fluid pressure (p) 110 minus the confining stress ( ⁇ c ) 112 as indicated by the multi-layer elasticity equation 104 in FIG. 20 , in that portion of the fracture.
  • the receiver signal comprises the stress and strain experienced at the receiver location due to the applied pressure at the source element location.
  • Many of these signals between source and receiver element are duplicated in the numerical mesh, and in these cases, special algorithms are designed to dramatically minimize the volume of storage required, so that only unique signals between different elements need to be stored.
  • the signals between each unique pair of receiver and source elements are stored in the computer memory or on a physical storage device in a matrix.
  • the hydraulic fracture propagation numerical method is designed so that the fracture propagates in a finite number of time steps.
  • the signal matrix is invoked, extracting those signals which are active over the part of the numerical mesh that is covered by the current configuration of the hydraulic fracture. This matrix is then used to build a system of numerical equations that are solved for the fracture width at the current time—at each active element location.
  • the combined system of equations must be coupled together in an efficient manner, so that they feed off each other until a balanced solution of fluid pressure and fracture width is obtained at each time step.
  • This coupling between the two equation systems is accomplished by means of a special numerical algorithm that efficiently and accurately ensures that the correct solution is obtained.
  • the entire system is designed to ensure that no fluid or proppant is unaccounted for in any time step.
  • ⁇ kr l 1 2 ⁇ ( u k , r l + u r , k l ) ( 3 )
  • FIG. 28 representing a schematic showing multiple parallel layers in a three dimensional case.
  • u ⁇ j l ⁇ ( m , n , y ) ⁇ - ⁇ ⁇ ⁇ ⁇ - ⁇ ⁇ ⁇ e i ⁇ ( mx + nz ) ⁇ u j l ⁇ ( x , y , z ) ⁇ d x ⁇ d z ( 4 )
  • d jr l and f jr l are constants that depend on the material constants of the layer
  • the ⁇ j l are the roots of the characteristic equation for the system of ordinary differential equations
  • the A j l (k) are free parameters of the solution that are determined by the forcing terms b i l in (1) and the interface conditions prescribed at the boundary between each of the layers (e.g., bonded, frictionless, etc.).
  • ⁇ ⁇ pq l ⁇ j ⁇ ( s jpq l ) ⁇ e ⁇ j l ⁇ ky ⁇ A j l ⁇ ( k ) ( 8 )
  • ⁇ ⁇ pq l ⁇ j ⁇ ( s jpq l + t jpq l ⁇ ky ) ⁇ e ⁇ j l ⁇ ky ⁇ A j l ⁇ ( k ) ( 9 )
  • u j l 1 ( 2 ⁇ ⁇ ) 2 ⁇ ⁇ - ⁇ ⁇ ⁇ ⁇ - ⁇ ⁇ ⁇ e - i ⁇ ( mx + nz ) ⁇ u ⁇ j l ⁇ d m ⁇ d n ( 10 ) of equations (6)-(9).
  • a j l ( ⁇ ) represents the “high frequency” components of the spectrum of the solution, which represents the singular part of the solution in real space.
  • the inversion process can be achieved by evaluating integrals of the form:
  • I j pl 1 ( 2 ⁇ ⁇ ) 2 ⁇ ⁇ - ⁇ ⁇ ⁇ ⁇ - ⁇ ⁇ ⁇ [ A j l ⁇ ( k ) - A j l ⁇ ( ⁇ ) ] ⁇ k p ⁇ e ⁇ j l ⁇ ky ⁇ e - i ⁇ ( mx + nz ) ⁇ d m ⁇ d n ⁇ ⁇ or ⁇ ⁇ I j pl ⁇ 1 ( 2 ⁇ ⁇ ) 2 ⁇ ⁇ j ⁇ a jr l ⁇ - ⁇ ⁇ ⁇ ⁇ - ⁇ ⁇ ⁇ k p ⁇ e k ⁇ ( ⁇ j l ⁇ y + b jr l ) ⁇ e - i ⁇ ( mx + nz ) ⁇ d m ⁇ d n ( 12 )
  • the shifted components ⁇ j l y+b jr l in (12) represent a finite number of complex images that approximate what would be an L-fold infinite Fourier Series (for L layers) that would be required to represent the Green's function in a closed form using the method of images. Typically three or four complex images suffice to give a high order of accuracy.
  • FIG. 26 representing a Schematic showing a duplicity of influence coefficients for a multiple layer problem.
  • the DD elements in layer N denoted by the solid circle 202 , the shaded circle 204 , the unshaded circle 206 , and the arrow 208 , in FIG. 26 each have the same set of exponential expansion coefficients a jr l and b jr l in layer N-1, where the arrow 208 indicates where the expansion coefficients a jr l and b jr l are identical.
  • a DD influence at a specified point within a given layer is constructed by constructing a pseudo interface parallel to the layering across which there can be a jump in the displacement field.
  • a jump in u y is specified, whereas to construct a shear or ride DD a jump in u x or u z is specified. This technique limits the orientation of DD components to be aligned with the pseudo interface that is parallel to the layering.
  • DD components that specify jumps in the displacement field which are across interfaces that are not parallel to the layering in the material.
  • DD field of a desired orientation by utilizing the duality relationship between the stresses due to a force discontinuity (or point force) and the displacements due to a displacement discontinuity.
  • the key idea here is to construct a planar Green's function or influence matrix, which represents the influences of all the DDs that lie in a vertical fracture plane.
  • the influence matrix will only represent the mutual influences on one another of DDs that lie in the fracture plane. However, it will implicitly contain the information about the variations in material properties due to the layering.
  • a reduced influence matrix can be constructed by any numerical method, including that proposed above, which can represent rigorously the changes in material properties between the layers, for example the finite element method or a boundary integral method in which elements are placed on the interfaces between material layers.
  • the in-plane influence coefficients would be calculated by means of the following procedure.
  • the value of the ij th DD would be assigned a value of unity and all the other fracture plane DDs would be set equal to zero.
  • the boundary integral or finite element solution on the interfaces between the material layers would then be determined so as to ensure that there is compatibility in the displacements between the material layers as well as a balance in the forces between the layers at the interface.
  • the numerical procedure outlined to construct the desired in-plane influence matrix would take a considerable amount of time to compute. Indeed, such a process would likely exclude the possibility of real time processing with current personal computers or workstations, but could conceivably be performed in a batch mode prior to the desired simulation being performed.
  • the semi-analytic method outlined above would be much more efficient, as the fully three-dimensional (or two-dimensional, in the case of plane strain or plane stress), problem that needs to be solved to calculate the influence of each DD element has been effectively reduced to a one-dimensional one.
  • Numerical models for multi-layered materials require that the interface between each material type is numerically “stitched” together by means of elements.
  • a boundary element method implementation would require that each interface between different material types is discretized into elements.
  • a finite element or finite difference method implementation would require that the entire domain is discretized into elements.
  • the material interfaces are indirectly taken into account without the requirement of explicitly including elements away from the surface of the crack or fracture. The implication thereof, is a dramatic reduction in the number of equations to be solved, with a commensurate dramatic decrease in computer processing time. In addition, accuracy of the solution is maintained.
  • One aspect of this invention which distinguishes it from previous work is that it is capable of solving problems in multi-layered elastic materials with arbitrarily inclined multiple cracks or fractures in two or three-dimensional space.
  • ‘elements’ can intersect layers. This is accomplished by taking special care of the mathematical stress/strain relationships across interfaces in such a way as to obtain the correct physical response for the element which is located across the interface(s).
  • references 1-3 below are classic papers that establish the Fourier scheme to solve elastic multi-layer problems, but do not utilize the inversion scheme proposed here.
  • references 1 and 2 a propagator matrix approach is suggested to solve the system of algebraic equations necessary for the Fourier scheme, but this particular scheme will become unstable for problems with many layers.
  • References 4 and 5 use exponential approximation for inversion. The methods in those references do not give rise to the complex images generated by the algorithm presented in this invention that so efficiently represent the effect of many layers.
  • Reference 5 extends the propagator approach used in references 1 and 2 to solve the algebraic equations of the Fourier method.
  • Reference 5 discloses an inversion scheme that is an integral part of the propagator method. This method involves an exponential approximation, similar to that proposed in this patent, but it is applied to only one part of the propagator equations. As a result, a least squares fit of many terms (more than 50) is required to yield reasonable results using the teachings of this reference. Apart from stability issues involved with exponential fitting, the large number of terms would probably be less efficient than using direct numerical integration for inversion. The exponential fit of the spectral coefficients that we propose involves less than five terms.
  • References 6 through 10 extend the Fourier method to transversely isotropic media.
  • References 7-10 use the propagator matrix for solving the algebraic equations, while reference 6 proposes direct solution. All these methods of solution would be numerically unstable for problems with many thick layers. While reference 10 proposes numerical inversion using continued fractions, little mention is made of the inversion process.
  • References 11 and 12 describe methodologies for multi-layer dielectric materials containing point electrical charges, or line charge distributions aligned parallel to interfaces (i.e., with different Green's functions to those used in elasticity).
  • Reference 13 describes a so-called “sweeping” algorithm to solve layered systems.
  • the method disclosed in reference 13 is essentially the classic block LU decomposition for a block tri-diagonal system.
  • this algorithm we use this algorithm to obtain a stable solution of the algebraic equations that determine the Fourier spectral coefficients in each of the layers. This method is particularly effective for problems in which the layers are thick or the wave-numbers are large.
  • references 1 through 15 are incorporated herein by reference:

Abstract

A method, practiced by a Hydraulic Fracturing Simulator software, is adapted for simulating a hydraulic fracture in an Earth formation, the formation including a multilayered reservoir, a mesh adapted to overlay the fracture in the formation thereby defining a plurality of elements, the mesh and the fracture defining one or more partially active elements, the method comprising the step of: setting up an Influence Coefficient Matrix, the matrix having a plurality of numbers, the plurality of numbers of the Influence Coefficient Matrix being calculated and determined such that the method, practiced by the Hydraulic Fracturing Simulator software for simulating the hydraulic fracture, models, or take into account an existence of, the multilayered reservoir and the partially active elements.

Description

CROSS REFERENCE TO RELATED APPLICATIONS
This application is a Continuation in Part of prior application Ser. No. 09/301,961, filed Apr. 29, 1999, now U.S. Pat. No. 6,876,959 entitled “Method and Apparatus for Hydraulic Fracturing Analysis and Design”; this application is also a Continuation in Part of prior application Ser. No. 10/831,799, filed Apr. 26, 2004, now U.S. Pat. No. 7,063,147 entitled “Method and Apparatus and Program Storage Device for Front Tracking in Hydraulic Fracturing Simulators”. This application is related to a application Ser. No. 11/096,047, filed Mar. 31, 2005, now U.S. Pat. No. 7,386,431 entitled “Method System and Program Storage Device for Simulating Interfacial Slip In a Hydraulic Fracturing Simulator Software,” having the same inventors and the same assignee as the present invention.
BACKGROUND
The subject matter set forth in this specification relates to a Hydraulic Fracturing Simulator Software adapted for designing and monitoring and evaluating reservoir, for example petroleum reservoir, fracturing, and, in particular, to a Hydraulic Fracturing Simulator Software adapted for simulating a multilayer reservoir and partially active elements of a mesh overlaying a fracture footprint when designing and monitoring and evaluating reservoir fracturing.
In hydraulic fracturing, thousands of gallons of fluid are forced under high pressure underground to split open the rock in a subterranean formation, a process which, for petroleum reservoirs, is known as ‘petroleum reservoir fracturing’ associated with ‘a fracturing event’. Proppant or propping agent is carried into the fracture by a viscosified fluid, and deposited into the fracture. Proppant provides a permeable flow channel for formation fluids, such as oil and gas, to travel up the wellbore and to the Earth's surface. Fracturing involves many variables, including: viscosity of the fracturing fluid, rate of leak-off of fracturing fluid into the reservoir, proppant carrying capacity of the fluid, viscosity of the fluid as a function of temperature, time history of fluid volumes (i.e., the amount of fluid pumped over a given period of time), time history of proppant volumes, fluid physical constants, proppant properties, and the geological properties of various zones in the reservoir.
A Hydraulic Fracturing Simulator software is capable of modeling the ‘fracturing event’. In fact, the Hydraulic Fracturing Simulator software will design and monitor and evaluate the ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’ from a time extending before, during, and after the ‘fracturing event’. However, when the Hydraulic Fracturing Simulator software designs and monitors and evaluates the ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’, the Hydraulic Fracturing Simulator software should function to model a ‘multilayered reservoir’ when modeling the ‘fracturing event’. In addition, the Hydraulic Fracturing Simulator software should also function to model ‘partially active elements’ (partially enclosed by the fracture footprint, as opposed to ‘fully active elements’ fully enclosed by the fracture footprint), among a plurality of elements of a mesh overlaying a fracture footprint, when the Hydraulic Fracturing Simulator software models the ‘fracturing event’.
Therefore, a Hydraulic Fracturing Simulator Software is needed that is capable of modeling or simulating a ‘multilayered reservoir’ during a time when the Hydraulic Fracturing Simulator Software is designing and monitoring and evaluating ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’.
In addition, a Hydraulic Fracturing Simulator Software is needed that is capable of modeling or simulating ‘partially active elements’, among a plurality of elements of a mesh overlaying a fracture footprint, during a time when the Hydraulic Fracturing Simulator Software is designing and monitoring and evaluating ‘petroleum reservoir fracturing’ associated with the ‘fracturing event’.
SUMMARY
One aspect of the present invention involves a program storage device readable by a machine tangibly embodying a program of instructions executable by the machine to perform method steps of simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, the method steps comprising: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the multilayered reservoir.
Another aspect of the present invention involves a method of simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, comprising the step of: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the multilayered reservoir.
Another aspect of the present invention involves a program storage device readable by a machine tangibly embodying a program of instructions executable by the machine to perform method steps of simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, the method steps comprising: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
Another aspect of the present invention involves a method of simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, comprising the step of: calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the method steps for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
A further aspect of the present invention involves a system adapted for simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir, a mesh overlaying the fracture thereby defining a plurality of elements, comprising: apparatus adapted for calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the system adapted for simulating the hydraulic fracture in the formation will take into account an existence of the multilayered reservoir.
A further aspect of the present invention involves a system adapted for simulating a hydraulic fracture in an Earth formation where a mesh overlays the fracture, the mesh and the fracture collectively defining one or more partially active elements, comprising: apparatus adapted for calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, the plurality of numbers of the influence coefficient matrix being calculated and determined such that the system adapted for simulating the hydraulic fracture will take into account an existence of the one or more partially active elements.
Further scope of applicability will become apparent from the detailed description presented hereinafter. It should be understood, however, that the detailed description and the specific examples set forth below are given by way of illustration only, since various changes and modifications within the spirit and scope of the ‘Hydraulic Fracturing Simulator Software’, as described and claimed in this specification, will become obvious to one skilled in the art from a reading of the following detailed description.
BRIEF DESCRIPTION OF THE DRAWINGS
A full understanding will be obtained from the detailed description presented hereinbelow, and the accompanying drawings, which are given by way of illustration only and are not intended to be limitative to any extent, and wherein:
FIGS. 1 through 3 illustrate a typical Hydraulic Fracturing (HF) job in a wellbore;
FIGS. 4 through 6 illustrate the fracture footprint created in the formation penetrated by the wellbore when the HF job is pumped;
FIGS. 7 through 9 illustrate how a mesh comprised of a plurality of grid cells will overlay on top of the fracture footprint of FIGS. 4 through 6, each grid cell of the mesh having a width and pressure, some of the grids cells called ‘tip elements’ being intersected by a perimeter of the fracture footprint, the tip elements having a width and a pressure (w, p), a portion of each tip element having fracturing fluid disposed therein;
FIG. 10 illustrates a mesh overlayed on top of a fracture footprint, one or more ‘fully active elements’ being enclosed by the fracture footprint, one or more ‘partially active elements’ being partially enclosed by the fracture footprint, the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘partially active elements’;
FIGS. 11 and 12 illustrate, in greater detail, the ‘partially active elements’ of FIG. 10, the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘partially active elements’;
FIGS. 13 and 14 illustrate a ‘multi-layered reservoir’, the Hydraulic Fracturing Simulator software of FIG. 16 modeling the ‘multi-layered reservoir’;
FIG. 15 illustrates an apparatus used in connection with a Hydraulic Fracturing (HF) job adapted for fracturing a formation penetrated by a wellbore, the apparatus including a computer system for storing a software called a ‘Hydraulic Fracturing Simulator software’;
FIG. 16 illustrates the computer system of FIG. 15 which stores the software called a ‘Hydraulic Fracturing Simulator software’ adapted for modeling a ‘multilayered reservoir’ and for modeling ‘partially active elements’ of a mesh overlaying a fracture footprint during a ‘petroleum reservoir fracturing’ event;
FIG. 17 illustrates in greater detail the ‘Other Data’ of FIG. 16;
FIG. 18 illustrates the ‘Other Software Instructions’ of FIG. 16;
FIG. 19 illustrates a construction of the ‘Hydraulic Fracturing Simulator software’ of FIG. 16;
FIG. 20 illustrates the function associated with the ‘Set Up Influence Coefficient Matrix [C]’ step associated with the construction of the ‘Hydraulic Fracturing Simulator software’ which is illustrated in FIG. 19;
FIGS. 21 and 22 illustrate a more detailed construction of the ‘Hydraulic Fracturing Simulator software’ which is illustrated in FIG. 19;
FIGS. 23 and 24 illustrate a construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIGS. 19 and 21; and
FIGS. 25-28 are used during a discussion of a detailed construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 in FIGS. 19 and 21.
DESCRIPTION
This specification discloses a Hydraulic Fracturing Simulator software adapted to be stored in a memory of a program storage device of a computer system for modeling and simulating a multilayered reservoir and for modeling and simulating partially active elements of a mesh overlaying a fracture footprint when the Hydraulic Fracturing Simulator software is designing and monitoring and evaluating petroleum reservoir fracturing. The Hydraulic Fracturing Simulator software includes a first step further including the step of Setting Up an Influence Coefficient Matrix and a second step further including two iteration loops whereby, in a first iteration loop at a first time step, a second iteration loop will continuously calculate fracture width for each element of the mesh given a previously determined fluid pressure and will continuously calculate fluid pressure for each element of the mesh given a previously determined fracture width until convergence of the solution of width and pressure is reached at which time the fracture footprint is updated to a first value and ‘output data’ is generated, then, in the first iteration loop at a second time step, the second iteration loop calculates fracture width for each element of the mesh given a previously determined fluid pressure, and fluid pressure for each element of the mesh is calculated given a previously determined fracture width at which time the fracture footprint is updated to a second value and ‘output data’ is generated, and the process repeats, at which point, ‘additional output data’ is generated. The ‘output data’ includes the previously determined first value of the fracture footprint and the previously determined second value of the fracture footprint. The ‘output data’, representative of the updated fracture footprint at each of the incremented time steps, is recorded or displayed on a recorder or display device. In the first step including the step of ‘setting Up an Influence Coefficient Matrix’, the Influence Coefficient Matrix is calculated and determined in a special way such that, when the Influence Coefficient Matrix is determined, and during any ‘fracturing event’, the Hydraulic Fracturing Simulator software will model and simulate the following: (1) a multilayered reservoir of the type illustrated in FIGS. 13 and 14, and (2) partially active elements of a mesh which overlays a fracture footprint of the type illustrated in FIGS. 10, 11, and 12.
Referring to FIG. 1, a perforating gun 10 is disposed in a wellbore 12 and a packer 14 isolates a plurality of shaped charges 16 of the perforating gun 10 downhole in relation to the environment uphole. The shaped charges 16 detonate and a corresponding plurality of perforations 18 are produced in a formation 20 penetrated by the wellbore 12.
Referring to FIG. 2, when the formation 20 is perforated, a fracturing fluid 22 is pumped downhole into the perforations 18 in accordance with a particular pumping schedule 24. In response thereto, the formation 20 surrounding the perforations 18 is fractured.
Referring to FIG. 3, when the formation 20 surrounding the perforations 18 is fractured, oil or other hydrocarbon deposits 26 begin to flow from the fractures, into the perforations 18, into the wellbore 12, and uphole to the surface. The oil or other hydrocarbon deposits flow at a certain ‘production rate’ 28 (e.g., in m3/day).
Referring to FIG. 4, when the wellbore 12 of FIG. 2 is fractured, pump trucks 30 situated at the surface of the wellbore will pump fracturing fluid down a tubing and into the perforations 18 in the formation 20 penetrated by the wellbore, as shown in FIG. 2. The formation 20 includes different layers, such as the different layers 42, one such layer being a perforated interval 40. In response thereto, at time t1, a fracture footprint 32 will be formed in the perforated interval 40 (and possibly in additional adjacent intervals 42) of a formation 20 penetrated by the wellbore 12. At time t2, the fracture footprint 32 will change to a new footprint 34 in the perforated interval 40 (and possibly in additional intervals 42) of a formation 20 penetrated by the wellbore 12. At time t3, the fracture footprint 34 will change to a new footprint 36 in the perforated interval 40 (and possibly in additional intervals 42) of a formation 20 penetrated by the wellbore 12. At time t4, the fracture footprint 36 will change to a new footprint 38 in the perforated interval 40 (and possibly in additional intervals 42) of a formation 20 penetrated by the wellbore 12.
Referring to FIGS. 5 and 6, referring initially to FIG. 5, a schematic three dimensional view of a fracture footprint, such as the fracture footprints 32-38 of FIG. 4, is illustrated. In FIG. 5, each fracture footprint 32-38 has a length ‘L’ and a height ‘H’ and a width ‘w’. In FIG. 6, the wellbore 12 is illustrated again, and a plurality of perforations 18 are shown in the formation 20 penetrated by the wellbore 12, as illustrated in FIGS. 1-3. As noted in FIG. 4, the formation 20 includes a plurality of different layers 42. In FIG. 6, when the pump trucks 30 of FIG. 4 pump the fracturing fluid into the perforations 18, a ‘fracture footprint’ 46 is created in the formation 20 which is identical to the fracture footprints 32, 34, 36, and 38 shown in FIG. 4 that are created, respectively, over the different periods of time t1, t2, t3, and t4. Note that the ‘fracture footprint’ 46 in FIG. 6 has a cross section 44, the cross section 44 having a fracture width ‘w’ similar to the width ‘w’ of the fracture footprint 32-38 shown in FIG. 5.
Referring to FIG. 7, recalling the fracture footprint 46 shown in FIG. 6, in FIG. 7, a mesh 48 comprised of a plurality of grid elements 48 a or grid cells 48 a is illustrated. The fracture footprint 46 is assumed (by the model of this specification) to lie in the 2D plane, although, in principle and in reality, the fracture footprint 46 does lie in the 3D plane. In FIG. 7, the mesh 48 is overlayed over the top of the fracture footprint 46 of FIG. 6. The mesh 48 includes a plurality of active elements or active grid cells 48 a 1 and a plurality of inactive elements or inactive grid cells 48 a 2. The active grid cells 48 a 1 will overlay the fracture footprint 46, whereas the inactive grid cells 48 a 2 will not overlay the fracture footprint 46. Each of the active grid elements or grid cells 48 a 1 of the mesh 48 has a width ‘w’ and a pressure ‘p’ assigned thereto, denoted by the symbol: (w, p). Therefore, each active grid cell 48 a 1 of the mesh 48 has a width/pressure value (w, p) assigned to that active grid cell. In FIG. 6, since the fracturing fluid propagating down the wellbore 12 enters the perforations 18 and creates the fracture footprint 46, in FIG. 7, each of the active grid cells 48 a 1 in the mesh 48 has a fracturing fluid disposed therein. In FIG. 7, there are two types of active grid cells 48 a 1: (1) an active grid cell 48 a 1 which is intersected by a perimeter 46 a of the fracture footprint 46, and (2) an active grid cell 48 a 1 which is not intersected by the perimeter 46 a of the fracture footprint 46. An active grid cell 48 a 1 that is intersected by the perimeter 46 a of the fracture footprint 46 is known as a ‘tip element’. For example, in FIG. 7, ‘tip element’ 50 is an active grid cell 48 a 1 that is intersected by the perimeter 46 a of the fracture footprint 46. An active grid cell 48 a 1 which is not intersected by the perimeter 46 a of the fracture footprint 46 has a volume which is wholly occupied by the fracturing fluid (i.e., 100% of the volume of the active grid cell is occupied by the fracturing fluid), where the fracturing ‘fluid’ may or may not include proppant. For example, in FIG. 7, active grid cell 52 is not intersected by the perimeter 46 a of the fracture footprint 46 and its volume is wholly occupied by the fracturing fluid (100% of the volume of the active grid cell 52 is occupied by fracturing fluid). However, an active grid cell 48 a 1 that is intersected by the perimeter 46 a (i.e., a ‘tip element’) has a volume that is occupied by ‘less than 100%’ of the fracturing fluid. For example, the active grid cell or ‘tip element’ 54 is intersected by the perimeter 46 a of the fracture footprint 46, however, only 45% of the volume of the active grid cell 54 is occupied by the fracturing fluid. In comparison, an inactive grid cell 48 a 2, such as inactive grid cell 55, has a volume which is wholly devoid of any fracturing fluid (0% of the volume of the inactive grid cell 55 is occupied by fracturing fluid). In prior pending application Ser. No. 10/831,799, filed Apr. 26, 2004, now U.S. Pat. No. 7,063,147 entitled “Method and Apparatus and Program Storage Device for Front Tracking in Hydraulic Fracturing Simulators” which is directed to the ‘VOF software’, the ‘VOF software’ of that prior pending application calculates, over a series of time steps, the ‘amount of fracturing fluid that is contained within each of the active grid cells 48 a 1 that are intersected by the perimeter 46 a of the fracture footprint 46’. That is, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 now U.S. Pat. No. 7,063,147 calculates, over the series of time steps, the ‘amount of fracturing fluid that is contained within each of the tip elements 50’. The ‘amount of fracturing fluid that is contained within each of the tip elements 50’ is calculated from the ‘fill fraction’, the ‘fill fraction’ being denoted by the letter ‘F’. For example, in FIG. 7, the ‘fill fraction’ F for the ‘tip element’ 54 is 45%. Therefore, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 calculates, over a series of time steps, the ‘fill fraction’ (F) for each of the ‘active grid cells 48 a 1 of the mesh 48 that are intersected by the perimeter 46 a of the fracture footprint 46’. That is, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 calculates, over the series of time steps, the ‘fill fraction’ (F) for each of the ‘tip elements’ 50 of the mesh 48 of FIG. 7. As a result, by calculating the ‘fill fraction’ (F) for each of the ‘tip elements’ 50 over a series of time steps, the amount or degree by which the perimeter 46 a of the footprint 46 of the fracture is expanding (or contracting), as a result of the pumping of the fracturing fluid into the perforations 18 in the formation 20 by the pump trucks 30 of FIG. 4, can be determined. The ‘VOF software’ of prior pending application Ser. No. 10/831,799 now U.S. Pat. No. 7,063,147 is embodied in steps 138 and 96 of FIG. 22. However, in this specification, the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21 (and the ‘Set Up Influence Coefficient Matrix [C]’ 102 of FIG. 19) is discussed in detail.
Referring to FIGS. 8 and 9, two more examples of a mesh 48 similar to the mesh 48 of FIG. 7 are illustrated. In FIG. 8, a mesh 48 is illustrated as overlaying the fracture footprint 46 having a perimeter 46 a. Fracturing fluid is disposed inside the perimeter 46 a, but the fracturing fluid is not disposed outside the perimeter 46 a. In FIG. 8, since the inactive grid cell 48 a 2 is disposed outside the perimeter 46 a, there is no fracturing fluid disposed inside the inactive grid cell 48 a 2 and, therefore, the ‘fill fraction’ F associated with the inactive grid cell 48 a 2 of FIG. 8 is ‘zero’ (F=0). In FIG. 8, the active grid cells 48 a 1 are disposed wholly within the perimeter 46 a (i.e., the perimeter 46 a does not intersect the active grid cells 48 a 1); therefore, the entire volume (i.e., 100%) of the active grid cells 48 a 1 is occupied by fracturing fluid and, as a result, the ‘fill fraction’ F associated with the active grid cells 48 a 1 in FIG. 8 is ‘1’ (F=1). However, in FIG. 8, let us analyze the active grid cell 56. The active grid cell 56 is intersected by the perimeter 46 a and, as a result, 80% of the volume of the active grid cell 56 is occupied by the fracturing fluid; therefore, the ‘fill fraction’ F for the active grid cell 56 is 0.8 (F=0.8). In FIG. 9, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 also calculates the volume of an active grid cell occupied by a first type of fluid and the volume of that same active grid cell occupied by a second type of fluid. For example, active grid cell 58 includes a first volume of 80% occupied by a first type of fluid, and a second volume of 20% occupied by a second type of fluid. The ‘VOF software’ of prior pending application Ser. No. 10/831,799 calculates, over a series of time steps, the ‘fill fraction’ (F) for each of the ‘active grid cells 48 a 1 that are intersected by the perimeter 46 a of the fracture footprint 46’ in the mesh 48; that is, the ‘VOF software’ of prior pending application Ser. No. 10/831,799 will calculate, over the series of time steps, the ‘fill fraction’ (F) for each of the ‘tip elements’ in the mesh 48 shown in FIGS. 7, 8, and 9. As a result, the amount or degree by which the perimeter 46 a of the footprint 46 is expanding (or contracting), in response to the pumping of fracturing fluid into the perforations 18, can be determined.
Referring to FIG. 10, ‘fully active elements’ and ‘inactive elements’ and ‘partially active elements’ are illustrated. In FIG. 10, a wellbore 12 is illustrated, and a fracture 17 grows in a direction that is directed away from the wellbore 12. A mesh 19 overlays the fracture 17. In the mesh 19, an ‘inactive element’ 21 is illustrated, and a ‘fully active element’ 23 is also illustrated. Note that FIG. 10 relates to ‘fracture growth’ of the fracture 17 where the ‘growth’ of the fracture is in a direction that is either directed away from the wellbore 12 (i.e., expansion) or is directed toward the wellbore 12 (contraction). That is, the fracture 17 represents a moving boundary that is growing in a direction which is either directed away from the wellbore 12 (expansion) or is directed toward the wellbore 12 (contraction). In FIG. 10, ‘partially active elements’ 25 and 27 and 29 and 43 are illustrated. In the ‘partially active element’ 25 and 27, only a portion 25 a and 27 a of the ‘partially active element’ 25 and 27 is disposed inside of the boundary 17 of the fracture 17. Since only portions 25 a and 27 a of ‘partially active elements’ 25 and 27 are disposed inside the boundary 17 of fracture 17, the elements 25 and 27 of mesh 19 of FIG. 10 are each known as a ‘partially active element’.
FIG. 11 illustrates a close up view of the ‘partially active element’ 29 of FIG. 10. The ‘partially active element’ 29 has a fracture leading edge 31 with crossing points 33 and 35, respectively. Straight line 37 is erected between the crossing points, and it forms the boundary for the active portion 39 of the element 29 and the inactive portion 41 of the element 29.
FIG. 12 illustrates the same features for the ‘partially active element’ 43 of FIG. 10 as did FIG. 11 with respect to the ‘partially active element’ 29 of FIG. 10. In FIG. 12, a ‘partially active element’ 43 has a fracture leading edge 31 with crossing points 33 and 35, respectively. Straight line 37 is erected between the crossing points 33, 35, and it forms the boundary for the active portion 39 of the element 43 and the inactive portion 41 of the element 43.
The concept of ‘partially active elements’, such as the ‘partially active elements’ 25, 27, 29, and 43 shown in FIGS. 10, 11, and 12, is set forth in the following two publications, each of which is incorporated herein by reference: (1) Ryder, J. A. and Napier, J. A. L., 1985, Error Analysis and Design of a Large Scale Tabular Mining Stress Analyzer, Proceedings of the Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Japan, [Balkema] 1549-1555, the disclosure of which is incorporated by reference into the specification of this application; and (2) Ryder, J. A., Eds.: E. G., Beer, J. R. Booker, and J. P. Carter, Optimal Iteration Schemes Suitable for General Non-linear Boundary Element Modeling Applications: Proceedings of the 7th International Conference on Computer Methods and Advances in Geomechanics, Cairns, Australia, [Balkema], 1991, the disclosure of which is incorporated by reference into the specification of this application.
Referring to FIG. 13, a reservoir or Earth formation 20 is shown. In FIG. 13, pumping truck 47 provides fluid at high pressures and flow rates to wellhead 49, which is operably connected to the wellbore 12 at or near the ground surface 53. FIG. 13 illustrates the fracture boundary 55 at a particular time. Two other fracture fluid boundaries 57 and 59 also are indicated in FIG. 13. In FIG. 13, the reservoir 20 represents a ‘multi-layered reservoir’ 20 because the reservoir 20 in FIG. 13 includes the following plurality of layers 61-71 of Earth formation where the layers 61-71 represent various zones or laminations of underground geological formation: (1) a first layer 61 of ‘shale’, (2) a second layer 63 of ‘sandstone’, (3) a third layer 65 representing an ‘oil/gas pay zone’, (4) a fourth layer 67 of ‘shale’, (5) a fifth layer 69 representing an ‘oil/gas pay zone, and (6) a sixth layer 71 being a ‘water-bearing zone’. The fluid boundaries in FIG. 13 reveal separate types or compositions of pumped fluid. In FIG. 13, the fracture preferably is stopped prior to the water bearing zone 71 seen at the lower portion of FIG. 13.
Referring to FIG. 14, a ‘multilayered’, hydraulically fractured reservoir 20 or Earth formation 20 is penetrated by a wellbore 12. The wellbore 12 penetrates the Earth formation 20 and the hydraulic fracture 59 of FIG. 13 is shown between layers 20 f and 20 h. The fracture 59 has a fracture width 59 a which was created in response to the pressurization of the formation 20 by a fracturing fluid. In FIG. 14, the ‘multilayered’ reservoir/Earth formation 20 includes a plurality of Earth formation layers 20 a, 20 b, 20 c, 20 d, 20 e, 20 f, 20 g, 20 h, 20 i, and 20 j. Each of the layers 20 a-20 j can be characterized by a Young's Modulus (E) and a Poisson's Ratio (ν). The Young's Modulus (E) and the Poisson's Ratio (ν) describe the elastic properties of each of the layers 20 a-20 j in FIG. 14. For example, the elastic properties of layer 20 a can be characterized by Young's Modulus and Poisson's Ratio (Ea, νa), the elastic properties of layer 20 b being characterized by Young's Modulus and Poisson's Ratio (Eb, νb), the elastic properties of layer 20 c being characterized by Young's Modulus and Poisson's Ratio (Ec, νc), the elastic properties of layer 20 d being characterized by Young's Modulus and Poisson's Ratio (Ed, νd), the elastic properties of layer 20 e being characterized by Young's Modulus and Poisson's Ratio (Ee, νe), the elastic properties of layer 20 f being characterized by Young's Modulus and Poisson's Ratio (Ef, νf), the elastic properties of layer 20 g being characterized by Young's Modulus and Poisson's Ratio (Eg, νg), the elastic properties of layer 20 h being characterized by Young's Modulus and Poisson's Ratio (Eh, νh), the elastic properties of layer 20 i being characterized by Young's Modulus and Poisson's Ratio (Ei, νi), and the elastic properties of layer 20 j being characterized by Young's Modulus and Poisson's Ratio (Ej, νj).
Referring to FIG. 15, the pump trucks 30 of FIG. 4 will pump a fracturing fluid 62 (usually including proppant) down the wellbore 12 of FIG. 4 in accordance with a pumping schedule 60 (an example used in connection with this discussion). The fracturing fluid 62 will enter the perforations 18, and, responsive thereto, create a ‘fracture footprint’ 46, similar to the fracture footprint 46 shown in FIG. 6. Micro-seismic data sensor(s) 64 and tiltmeter data or other sensor(s) 66 will monitor the approximate geometry of the fracture footprint 46 and, responsive thereto, the sensor(s) 64 and 66 will generate output signals, the micro-seismic data sensor(s) 64 generating a micro-seismic data output signal 64 a, the tiltmeter data sensor(s) 66 generating a tiltmeter data output signal 66 a, and the pumping schedule 60 generating a pumping schedule output signal 60 a representative of the pumping schedule 60. The pumping schedule output signal 60 a, the tiltmeter data output signal 66 a, and the micro-seismic data output signal 64 a are time line merged via a ‘time line merging’ step 68. In this ‘time line merging’ step 68, the pumping schedule output signal 60 a, the tiltmeter data output signal 66 a, and the micro-seismic data output signal 64 a are ‘time synchronized’ in a particular manner such that the tiltmeter data output signal(s) 66 a and the micro-seismic data output signal(s) 64 a are synchronized with the times present in the pumping schedule 60. When the ‘time line merging’ step 68 is complete, a ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 is generated which is provided as ‘input data’ to a ‘computer system’ 72 optionally disposed within a truck 74 situated at the site of the wellbore 12, such as a monitoring truck 74 or a ‘FracCAT vehicle’ 74 (which is a vehicle with software that monitors and measures the fracture and controls the fracture treatment).
Referring to FIG. 16, the ‘computer system’ 72 which is optionally disposed within the truck 74 of FIG. 15, such as the ‘FracCAT vehicle’ 74, is illustrated. In FIG. 16, recall that the ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 of FIG. 15 is provided as ‘input data’ to the computer system 72 disposed within the truck 74, the output signal 70 being comprised of a time line merged pumping schedule and tiltmeter data and micro-seismic data plus other data including downhole temperature and bottom hole pressure. The computer system 72 of FIG. 16 includes a processor 72 a operatively connected to a system bus, a memory or other program storage device 72 b operatively connected to the system bus, and a recorder or display device 72 c operatively connected to the system bus. The memory or other program storage device 72 b stores the following software (76, 78, and 80): a Planar 3D User Interface software 76, a Planar 3D ‘engine’ or software 78, and a Hydraulic Fracturing Simulator Software 80. The software 76, 78, and 80, which is stored in the memory 72 b of FIG. 16, can be initially stored on a CD-ROM, where that CD-ROM is also a ‘program storage device’. That CD-ROM can be inserted into the computer system 72, and, then, the software 76, 78, and 80, which includes the Hydraulic Fracturing Simulator Software 80, can be loaded from that CD-ROM and into the memory/program storage device 72 b of the computer system 72 of FIG. 16. The Hydraulic Fracture Simulator Software 80 will be described in detail in the following paragraphs. The computer system 72 of FIG. 16 receives Input Data 82, including: (1) the time line merged pumping schedule, tiltmeter data, and micro-seismic data 84 (which corresponds to the ‘time line merged pumping schedule and tiltmeter data and micro-seismic data’ output signal 70 of FIG. 15), and (2) Other Data 86. The computer system 72 also receives Other Software Instructions 88. The processor 72 a will execute the Hydraulic Fracturing Simulator Software 80 (including the Planar 3D User Interface software 76 and the Planar 3D ‘engine’ 78), while simultaneously using the Input Data 82 and Other Software Instructions 88; and, responsive thereto, the recorder or display device 72 c will generate a set of ‘Output Data’ 72 c 1 which is adapted to be recorded by or displayed on the recorder or display device 72 c. The computer system 72 may be a personal computer (PC), a workstation, or a mainframe. Examples of possible workstations include a Silicon Graphics Indigo 2 workstation or a Sun SPARC workstation or a Sun ULTRA workstation or a Sun BLADE workstation. The memory or program storage device 72 b is a computer readable medium or a program storage device which is readable by a machine, such as the processor 72 a. The processor 72 a may be, for example, a microprocessor, microcontroller, or a mainframe or workstation processor. The memory or program storage device 72 b, which stores the Hydraulic Fracturing Simulator Software 80, may be, for example, a hard disk, ROM, CD-ROM, DRAM, or other RAM, flash memory, magnetic storage, optical storage, registers, or other volatile and/or non-volatile memory.
Referring to FIGS. 17 and 18, the Other Data 86 and the Other Software Instructions 88 of FIG. 16 are illustrated. In FIG. 17, the Other Data 86 will include: layer confining stresses and properties, perforating interval and depth, wellbore data, fluid and proppant properties, time history of fluid volumes to be pumped, time history of proppant volumes to be pumped, and logs identifying the identity, properties, and location of geological zones. In FIG. 18, the Other Software Instructions 88 include instructions to calculate values representing physical dimensions of the fracture and pressures inside the fracture.
Referring to FIG. 19, a construction of the Hydraulic Fracturing Simulator Software 80 disposed within the Planar 3D software 78 and the Planar 3D User Interface 76 is illustrated. In FIG. 19, the Input Data 82 is provided to the Planar 3D User Interface 76, to the Planar 3D software 78, and to the Hydraulic Fracturing Simulator Software 80. The Hydraulic Fracturing Simulator Software 80 includes: an initial ‘time stepping’ step 90, a second step 92 which calculates fracture width (w) given the fluid pressure (p); a third step 94 which calculates fluid pressure (p) given the fracture width (w); and a fourth step 96 adapted to update the fracture footprint. The fracture ‘footprint’ is shown in FIGS. 4 and 13 at times t1, t2, t3, . . . , tn (for example at times t1, t2, t3, and t4 in FIG. 4 and at times corresponding to the three boundaries 59, 57 and 55 in FIG. 13). In operation, note the time stepping loop 98 in FIG. 19 wherein the ‘time stepping’ step 90 will increment from a first time step (t1), to a second time step (t2), to a third time step (t3), . . . , and to an nth time step (tn). Operating within the first time step t1, the second step 92 and the third step 94 will iterate on each other (as indicated by arrow 100) until the fracture width (w) and the fluid pressure (p) are solved at every point (i.e., within each ‘element’ or ‘grid cell’) in the fracture geometry of FIGS. 4 and 13. When the iteration 100 of steps 92 and 94 has converged, the fracture footprint is updated in step 96. In step 97, the proppant concentration is calculated for each grid cell in the updated fracture footprint. The ‘output data’ 72 c 1, associated with the last iteration loop 100, is generated and stored.
The time step 90 is then incremented to the second time step t2, and, responsive thereto, the second step 92 and the third step 94 will then iterate again on each other (as indicated by arrow 100) until the fracture width (w) and the fluid pressure (p) are solved at every grid cell in the fracture geometry of FIGS. 4 and 13. When the iteration 100 of steps 92 and 94 has converged, the fracture footprint is updated in step 96. In step 97, the proppant concentration is calculated for each grid cell in the updated fracture footprint. The ‘output data’ 72 c 1, associated with the last iteration loop 100, is generated and stored. The time step 90 is then incremented once again to the third time step t3, the second step 92 and the third step 94 will then iterate again on each other (as indicated by arrow 100) until the fracture width (w) and the fluid pressure (p) are solved at every point (i.e., within each ‘element’ or ‘grid cell’) in the fracture geometry of FIGS. 4 and 13. When the iteration 100 of steps 92 and 94 is complete, the fracture footprint is updated in step 96. In step 97, the proppant concentration is calculated for each grid cell in the updated fracture footprint. The ‘output data’ 72 c 1, associated with the last iteration loop 100, is generated and stored. The time step 90 is then incremented once again, and the above process repeats until it has reached the end of the pumping schedule. ‘Convergence’ takes place when the ‘solution does not change from one iteration to the next’. The inner iteration loop 100 is solving for two things: (1) fracture width (w) using the ‘elasticity equation’, and (2) fluid pressure (p) using the ‘fluid flow equation’. Thus, the ‘solution does not change from one iteration to the next’ when the change in the ‘elasticity equation’ solution and when the change in the ‘fluid flow equation’ solution is below a ‘tolerance’. When the change in the ‘elasticity equation’ solution and the change in the ‘fluid flow equation’ solution is below the ‘tolerance’, we know that the inner iteration loop 100 has converged. Steps 92 and 94 can be solved in various ways, such as by iteration of two equations as shown here, or by direct substitution of the one equation into the other equation, or vice-versa.
The Output Data 72 c 1 is generated at the end of each time step. However, ‘additional output data’ 99 is generated when the time stepping in the outer iteration loop 98 is complete, the ‘additional output data’ 99 being used in subsequent calculations, such as in the generation of other graphical plots.
In FIG. 19, however, before the second step 92, a first step 102 is practiced, the first step 102 being called ‘Set Up Influence Coefficient Matrix [C]’ 102. The ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 will be discussed below with reference to FIG. 20 of the drawings.
Referring to FIG. 20, a ‘Multi-Layer Elasticity Equation’ 104 is illustrated. Recall from FIG. 19 that step 92 will calculate the fracture width (w) and step 94 will calculate the fluid pressure (p). In FIG. 20, the fracture width (w) of step 92 is actually calculated by using the ‘Elasticity Equation’ 104. In the ‘Elasticity Equation’ 104, the fracture width (w) 106 is calculated given the inverse of a ‘Matrix of Influence Coefficients [C]’ 108 (also known as an ‘Influence Coefficient Matrix’) multiplied by (fluid pressure p 110 minus confining stress σc 112), as follows:
{w}=[C] −1 {p−σ c}, where   Elasticity Equation 104
    • w is the fracture width 106,
    • [C] is the ‘matrix of influence coefficients’ 108,
    • p is the fluid pressure 110, and
    • σc is the confining stress 112.
In FIG. 20, a generic form of the ‘Influence Coefficient Matrix [C]’ 108 is shown, where the ‘Influence Coefficient Matrix [C]’ 108 has ‘M’ rows and ‘M’ columns. In FIG. 20, the ‘Influence Coefficient Matrix [C]’ 108 is fully populated with numbers 114 and is used to calculate fracture width (w) for step 92 in FIG. 19.
The ‘elasticity equation’ 104 did not previously take into account the existence of a ‘multilayered reservoir’ of the type illustrated in FIGS. 13 and 14, and the ‘elasticity equation’ 104 did not previously take into account the existence of ‘partially active elements’, such as the ‘partially active elements’ 25, 27, 29, and 43 of FIG. 10. However, in a real reservoir, a ‘multilayered reservoir’ does exist. Furthermore, when overlaying a mesh over a fracture footprint in the manner illustrated in FIG. 10, ‘partially active elements’ also do exist.
Therefore, the Hydraulic Fracturing Simulator Software 80 of FIG. 16 disclosed in this specification will include the effect of, and take into account the existence of, ‘multilayered reservoirs’ and ‘partially active elements’. Consequently, in order to take into account the existence of ‘multilayered reservoirs’ and ‘partially active elements’, it is necessary to execute the ‘set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19, associated with the Hydraulic Fracturing Simulator Software 80 of FIGS. 16 and 19. The ‘Set Up Influence Coefficient Matrix [C]’ of Step 102 of FIG. 19 will calculate the ‘Influence Coefficient Matrix [C]’ 108 in the ‘Elasticity Equation’ 104 of FIG. 20 in a special way in order to allow the Hydraulic Fracturing Simulator software 80 to include the effect of, or take into account the existence of, the ‘multilayered reservoirs’ and the ‘partially active elements’. When the ‘Influence Coefficient Matrix [C]’ 108 of the ‘elasticity equation’ 104 of FIG. 20 is calculated (by the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 associated with the Hydraulic Fracturing Simulator Software 80 of FIGS. 16 and 19) in this special way, the numbers 114 in the ‘Influence Coefficient Matrix [C]’ 108 will be changed accordingly; and, when the numbers 114 of the ‘Influence Coefficient Matrix [C]’ 108 are changed accordingly, the Hydraulic Fracturing Simulator Software 80 will then include the effect of, or take into account the existence of, the ‘multilayered reservoirs’ and the ‘partially active elements’. As disclosed in this specification, there is a systematic way for determining and changing the numbers 114 in the ‘Influence Coefficient Matrix’ 108 of FIG. 20 in order to include or take into account the existence of the ‘multilayered reservoirs’ and the ‘partially active elements’. Consequently, the following pages of this specification will present a ‘method’ (along with an accompanying ‘system’ and ‘program storage device’) for determining how the numbers 114 in the ‘Influence Coefficient Matrix [C]’ 108 will be changed (by the ‘Set Up Influence Coefficient Matrix’ step 102 of FIG. 19) for the ultimate purpose of allowing the Hydraulic Fracturing Simulator Software 80 to include the effect of, or take into account the existence of, the ‘multilayered reservoirs’ of FIGS. 13-14 (such as the ‘multilayered reservoir’ 20 having multiple layers 61-71 of FIG. 13) and the ‘partially active elements’ of FIG. 10 (such as the ‘partially active elements 25, 27, 29, and 43 of FIG. 10).
Referring to FIGS. 21 and 22, a more detailed construction of the Hydraulic Fracturing Simulator Software 80 of FIG. 19, which is disposed within the Planar 3D software 78 and the Planar 3D User Interface 76, is illustrated.
In FIG. 21, the input data 82 is provided from the human interface, such as the injection rate and the pumping schedule, which includes the injection rate as a function of time, proppant concentration, fluid viscosity, the geology or the properties of the elastic layers of the reservoir including the elastic constants comprising the Young's Modulus and the Poisson's Ratio (E, ν) as previously described, and the leakoff behavior. In FIG. 21, in the ‘Generate Layer Interface Locations’ step 116, the depths of each of the ‘interfaces’ of the ‘layers’ of FIG. 13 are calculated, such as the end of the sandstone layer, the end of the gas layer, etc. In the ‘Assign Layer Properties (Young's Modulus, Poisson's Ratio, Toughness, Leakoff Coefficients, Stresses)’ step 118, the Young's Modulus, Poisson's Ratio, Toughness, Leakoff Coefficients, and Stresses are assigned to each of the ‘layers’ (of FIG. 13) the depths of which were calculated in step 116. Therefore, a series of numbers (comprising the Young's Modulus, Poisson's Ratio, Toughness, Leakoff Coefficients, and Stresses) is assigned to each Earth formation ‘layer’ that is shown in FIG. 13. In the ‘Assign Maximum Expected Fracture Height and Extent of Fracture’ step 120 of FIG. 21, before any simulation is performed, the maximum possible ‘length’ to which the ‘fracture’ will propagate and the maximum possible ‘height’ to which the ‘fracture’ will propagate is assumed or introduced. Then a parent mesh is assigned to the ‘fracture surface’, where the parent mesh is divided into rectangular ‘elements’. See FIGS. 7 and 10 for examples of the parent mesh. The parent mesh is broken down into rows and columns comprising ‘grids’ or ‘elements’ which are generally rectangular in shape. It is assumed that the ‘fracture’ will grow into the ‘elements’ of the parent mesh, but no further. In the ‘Generate Numerical Parent Mesh’ step 122, the numerical parent mesh is generated; that is, in step 122, the ‘coordinates’ of each ‘grid cell’ or ‘element’ of the parent mesh are generated. The ‘coordinates’ will define where each ‘grid cell’ or ‘element’ exists within an axis system, such as the (x, y) axis system. Therefore, in step 122 of FIG. 21, the ‘coordinates’ of each of the ‘grid cells’ (such as ‘grid cells’ 23, 25, 29, and 43 shown in FIG. 10 or each of the ‘grid cells’ 48 a 1 and 52 shown in FIG. 7) are allocated. At this point, we know where our ‘layers’ exist, we know the numerical mesh which includes the rectangular mesh of elements, and we know where the fracture exists within the parent mesh including all their coordinates. In FIG. 21, the next step is the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 which is the ‘Set up Influence Coefficient Matrix [C]’ step 102 of FIG. 19. In step 102 of FIGS. 19 and 21, the Influence Coefficient Matrix [C] (108 of FIG. 20) is generated. As previously mentioned, the Influence Coefficient Matrix [C] (as shown in FIG. 20) comprises rows and columns of ‘numbers’ 114, the Influence Coefficient Matrix [C] being fully populated with such ‘numbers’ 114. Each of the ‘numbers’ 114 of the Influence Coefficient Matrix [C] will indicate how one ‘element’ of the mesh relates to another ‘element’ of the mesh in an ‘elastic manner’. For example, if a fracture exists within one ‘element’ of the mesh, what kind of stress or strain exists within another ‘element’ of the mesh in response thereto (a phenomenon that is called an ‘elasticity behavior’). Each of the ‘numbers’ 114 in the Influence Coefficient Matrix [C] describes this ‘elasticity behavior’. Therefore, the Influence Coefficient Matrix [C] (108 of FIG. 20) is a matrix that indicates how one ‘element’ of the parent mesh (such as the mesh shown in FIGS. 7 and 10) talks to another ‘element’ of the parent mesh in an ‘elastic manner’. Therefore, step 102 of FIGS. 19 and 21 (i.e., the ‘Set Up Influence Coefficient Matrix’ step 102 of FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21), which determines the Influence Coefficient Matrix 108 of FIG. 20, is responsible for allowing the Hydraulic Fracturing Simulator Software 80 of FIG. 16 to simulate the ‘multilayered reservoir’ of FIG. 13 and the ‘partially active elements’ of FIGS. 10, 11, and 12. Step 124 of FIG. 21 entitled the ‘Restart Option’ will allow the user to restart the simulator without having to recalculate the Influence Coefficient Matrix (it is stored on the hard drive). Such situations will occur if the reservoir properties are unchanged, but the user wants to change the injection schedule, for example.
In FIG. 22, the ‘Time=Time+Δ t’ step 90 is a ‘time stepping’ step. Note the loop 98 which represents a time stepping loop. The loop 98 controls how the fracture grows as a function of time during the duration of the ‘pumping schedule’, where the term ‘pumping schedule’ includes the ‘shut-in time’. In the ‘Assign Latest Fracture Extent (Active Elements)’ step 126, at any particular time step, the fracture will have propagated to some extent in the parent mesh, which is called the ‘footprint’ of the fracture; step 126 will sweep through all the ‘grid cells’ in the entire parent mesh, and then step 126 will determine if a particular ‘element’ or ‘grid cell’, within the perimeter of the ‘fracture footprint’, is an ‘inactive element’, or an ‘active element’, or a ‘partially active element’. In step 128 entitled ‘Extract Elastic Influence Coefficient Matrix for Current Fracture Geometry’, the Influence Coefficient Matrix [C] 108, as previously described with reference to step 102, contains ‘all possible combinations within the parent mesh’; and, at any particular time step, we need a subset of the aforementioned ‘all possible combinations within the parent mesh’; thus, step 128 extracts, from the Influence Coefficient Matrix [C] 108, a ‘submatrix’ which contains essential information regarding the current size of the ‘fracture footprint’. In step 130 entitled ‘Assign Special Characteristics to Coefficients for Tip Elements’, a ‘particular subset’ of the ‘grid cells’ or ‘elements’ of the parent mesh will cross a ‘perimeter of the fracture footprint’; in connection with only the ‘particular subset’ of the ‘grid cells’ which cross the ‘perimeter of the fracture footprint’, it is necessary to change some of the properties of the ‘submatrix’ (and recall that the ‘submatrix’ contains essential information regarding the current size of the ‘fracture footprint’); step 130 will make this change to the ‘submatrix’; consequently, step 130 in FIG. 22 entitled ‘Assign Special Characteristics to Coefficients for Tip Elements’ will deal specifically with the ‘partially active elements’ of FIGS. 10, 11, and 12; that is, step 130 in FIG. 22 will change the ‘partially active element submatrix’ (i.e., a submatrix which pertains specifically to the ‘partially active elements’) by ‘assigning special characteristics to the coefficients’ in the ‘partially active element submatrix’. In step 132 entitled ‘Assign Loading (Fluid Pressure) to Each Element in Current Mesh’, we now have a current time step and a current fracture footprint; in addition, each ‘active element’ has a ‘fluid pressure’; thus, step 132 will assign a ‘fluid pressure’ to each ‘element’ of the parent mesh. In step 92 entitled ‘Solve Elastic Equation System for Fracture Widths’, we use the ‘elasticity equation’ 104 of FIG. 20 (including their Fourier Transforms) to solve for the ‘fracture width’ in each ‘element’ or ‘grid cell’ of the mesh, at the current time step, given the ‘fluid pressure’ that was assigned to each ‘element’ of the mesh in step 132. In step 134 entitled ‘Calculate Influence Matrix for Fluid Flow in Current Mesh Using Widths from Above’, now that we have calculated the ‘fracture width’ in each ‘element’ at the current time step from step 92, we can do the same thing for the fluid flow equations in step 94. Therefore, in step 94 entitled ‘Solve Fluid Flow Equations for Fluid Pressure in Current Mesh’, we can solve for the ‘fluid pressures’. Note the loop 100 in FIG. 22 whereby, when the ‘fluid pressures’ are determined in step 94, those ‘fluid pressures’ determined from step 94 are used to solve for ‘fracture widths’ in step 92, and that loop 100 will continue to process as indicated until ‘global mass balance’ is achieved in step 136. Thus, from step 136 to step 132, an inner iteration will take place; in that inner iteration, we continue to iterate on fluid pressure and fracture width until convergence of the solution takes place at that time step; convergence takes place when the solution does not change from one iteration to the next. The inner iteration loop 100 is solving for two things: (1) fracture width using the ‘elasticity equation’, and (2) the fluid pressure using the ‘fluid flow equation’; when the change in each of these ‘equations’ is below a tolerance, we know that the inner iteration 100 has converged. Steps 138 and 96 involve updating for the next time step. Step 138 (‘Calculate Local Fracture Tip Velocity of Propagation’) and step 96 (‘Grow Fracture: Update New Layout’) are described in prior pending application Ser. No. 10/831,799, filed Apr. 26, 2004, now U.S. Pat. No. 7,063,147 directed to the ‘VOF Algorithm’, entitled “Method and Apparatus and Program Storage Device for Front Tracking in Hydraulic Fracturing Simulators”, the disclosure of which is incorporated by reference into the specification of this application. In step 97, the proppant concentration is calculated for each grid cell in the updated fracture footprint. The ‘output data’ 72 c 1, associated with the last iteration loop 100, is generated and stored. In step 140, if the time is not greater than the maximum allowed, time step 90 is incremented once again, and the above process repeats until it has reached the end of the pumping schedule. In step 140, if the ‘time’ is greater than the maximum allowed, ‘additional output data’ 99 is generated. The ‘output data’ 72 c 1 and the ‘additional output data’ 99 is stored in a storage device, such as the memory or program storage device 72 b of FIG. 16. In FIG. 22, the Output Data 72 c 1 includes, at each time step, the fluid pressure (p) and the fracture width (w) at each of the ‘elements’ or ‘grid cells’ 23, 25, 29, and 43 shown in FIG. 10 or each of the ‘elements’ or ‘grid cells’ 48 a 1 and 52 shown in FIG. 7 (again, at each time step). Since the Output Data 72 c 1 includes fluid pressure (p) and fracture width (w) for each ‘grid cell’ of FIGS. 7 and 10 at each time step, the Output Data 72 c 1 can also include a number of 2D or 3D plots representative of the fluid pressure (p) and the fracture width (w) at each of the ‘elements’ or ‘grid cells’ 23, 25, 29, and 43 shown in FIG. 10 or each of the ‘elements’ or ‘grid cells’ 48 a 1 and 52 shown in FIG. 7 at each time step. Steps 92, 134, and 94 can be solved in various ways, such as by iteration of the two equations shown here (i.e., the ‘elasticity equation’ and the ‘fluid flow equation’), or by direct substitution of one equation into the other equation, or vice-versa.
Recall that the ‘Set Up Influence Coefficient Matrix [C]’ step 102 in FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 will recalculate the Influence Coefficient Matrix 108 in FIG. 20 in a ‘special way’ in order to allow the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model, or take into account the existence of, ‘partially active elements’ (of the type illustrated in FIGS. 10, 11, and 12). The concept of ‘partially active elements’, such as the ‘partially active elements’ 25, 27, 29, and 43 shown in FIGS. 10, 11, and 12, is set forth in the following two publications, each of which is incorporated herein by reference: (1) Ryder, J. A. and Napier, J. A. L. 1985, Error Analysis and Design of a Large Scale Tabular Mining Stress Analyzer, Proceedings of the Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Japan, [Balkema] 1549-1555, the disclosure of which is incorporated by reference into the specification of this application; and (2) J. A. Ryder, Eds.: E. G. Beer, J. R. Booker, and J. P. Carter, Optimal Iteration Schemes Suitable for General Non-linear Boundary Element Modeling Applications: Proceedings of the 7th International Conference on Computer Methods and Advances in Geomechanics, Cairns, Australia,[Balkema], 1991, the disclosure of which is incorporated by reference into the specification of this application.
Referring to FIGS. 23 and 24, a construction of the ‘Set Up Influence Coefficient Matrix [C]’ step 102 of FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21, which would allow the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’ (of the type illustrated in FIG. 13), is illustrated. In FIGS. 23 and 24, recall that the ‘Set Up Influence Coefficient Matrix [C]’ step 102 in FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 will also recalculate the Influence Coefficient Matrix 108 in FIG. 20 in a ‘special way’ in order to allow the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’ (of the type illustrated in FIG. 13).
In FIGS. 23 and 24, a ‘construction of the Set Up Influence Coefficient Matrix step 102’ of FIGS. 19 and 21, for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’ (of the type illustrated in FIG. 13), is illustrated. The ‘construction of the Set Up Influence Coefficient Matrix step 102’ shown in FIGS. 23 and 24, for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 to model or take into account the existence of a ‘multi-layer reservoir, includes a ‘plurality of calculation steps’ 102, the ‘plurality of calculation steps’ 102 being practiced by the ‘Set Up Influence Coefficient Matrix’ step 102 in FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 when the Influence Coefficient Matrix 108 of FIG. 20 is calculated.
In FIG. 23, the ‘construction of the Set Up Influence Coefficient Matrix step 102’ of FIGS. 19 and 21, for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, includes ‘three basic steps’ 102 a, 102 b, and 102 c, as follows: (1) Find Spectral Coefficients in Each Layer, step 102 a, (2) Perform Exponential Approximation, step 102 b, and (3) Assemble an ‘Influence Coefficient Matrix’ (i.e., matrix 108 of FIG. 20) Using Exponential Expansion Coefficients, step 102 c. In step 200, the Influence Coefficient Matrix is sent to a storage device.
The above referenced ‘three basic steps’ 102 a, 102 b, and 102 c in FIG. 23, adapted for constructing an Influence Coefficient Matrix 108 for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, will be discussed in greater detail below with reference to FIG. 24.
In FIG. 24, a detailed construction of each of the ‘three basic steps’ 102 a, 102 b, and 102 c of FIG. 23 will be discussed.
In FIG. 24, the ‘plurality of calculation steps’ 102, which are practiced by the ‘Set Up Influence Coefficient Matrix’ step 102 in FIG. 19 and by the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 in FIG. 21 for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, include the following steps:
Find Spectral Coefficients in Each Layer, Step 102 a of FIG. 23:
The ‘Find Spectral Coefficients in Each Layer’ step 102 a of FIG. 23 includes the following step (1), as follows:
    • (1) In order to obtain the stresses and displacements due to a prescribed source [e.g., Displacement Discontinuity (DD) element], determine the spectral coefficients Aj l(k) by solving a system of algebraic equations (7) and (9) that express the continuity of tractions and displacements across layer interfaces; Equations (7) and (9) will be discussed in detail in the following paragraphs; however, equations (7) and (9) are duplicated below, as follows:
u ^ r l = j ( d jr l + f jr l y ) α j l ky A j l ( k ) ( 7 ) σ ^ pq l = j ( s jpq l + t jpq l ky ) α j l ky A j l ( k ) ( 9 )
Perform Exponential Approximation, Step 102 b of FIG. 23:
The ‘Perform Exponential Approximation’ step 102 b of FIG. 23 includes the following steps (2), (3), and (4), as follows:
    • (2) Determine the (free space) spectral coefficients Aj l(∞) by solving the algebraic equations (7) and (9), set forth above and discussed below, for an infinite homogeneous medium corresponding to the high frequency components associated with the prescribed Displacement Discontinuity (DD); the explicit expressions are identified in FIG. 24 by numeral 102 b 1; where λ and μ are Lame constants:
λ = vE ( 1 + v ) ( 1 - 2 v ) , and μ = E 2 ( 1 + v ) .
    • (3) Evaluate the low frequency spectral components which are defined by a further expression that is identified in FIG. 24 by numeral 102 b 2;
    • (4) Approximate these low frequency components (Aj l(k))LOW by a series of exponential functions by solving for the unknown constants ajr l and bjr l in the expansion that is identified in FIG. 24 by numeral 102 b 3.
      Assemble Influence Coefficient Matrix Using Exponential Expansion Coefficients, Step 102 c in FIG. 23:
The ‘Assemble Influence Coefficient Matrix Using Exponential Expansion Coefficients’ step 102 c of FIG. 23 includes the following steps (5), (6), and (7), as follows:
    • (5) Substitute these expansions for the low frequency components into equations (7) and (9) to obtain the expressions for the ‘displacement and stress components’ that are identified in FIG. 24 by numeral 102 c 1;
    • (6) Invert the displacements and stresses associated with the low frequency Fourier Transforms by evaluating integrals of the form that is identified in FIG. 24 by numeral 102 c 2, where ‘i’ is the imaginary number, and ‘e’ is the exponential operator; and
    • (7) Combine the low frequency displacement and stress components with the infinite space displacements and stresses in the manner which is identified in FIG. 24 by numeral 102 c 3.
Referring to FIGS. 25 through 28, a detailed construction of the ‘Set Up Influence Coefficient Matrix’ step 102 of FIG. 19 and the ‘Generate Elastic Influence Coefficient Matrix for Parent Mesh’ step 102 of FIG. 21, for the purpose of allowing the ‘Hydraulic Fracturing Simulator software’ 80 of FIGS. 16 and 19 to model or take into account the existence of a ‘multi-layer reservoir’, is set forth in the following paragraphs with reference to FIGS. 25 through 28. Note that step 102 of FIG. 19 is the same step as step 102 in FIG. 21.
The numerical algorithm employed in this invention comprises an efficient technique to determine the local width of a hydraulic fracture due to local pressure applied to the fracture faces by the injection of hydraulic fluid and proppant into the fracture. Further, a method to track the dimensions and width of said fracture as it grows as a function of time is shown. The hydraulic fracture(s) may span any number of layers in a laminated reservoir, with the restriction that all layers must be parallel to each other, as depicted for example in FIG. 23. Layers may be inclined at any angle to the horizontal.
FIG. 25 shows a section through multiple hydraulic fractures in a layered reservoir. The calculation of the fracture width due to the pressure from the injected fluids and proppant mixture is determined by taking into account, accurately and efficiently, the physical properties of each layer comprising the laminated reservoir. The technique used to calculate the relationship between the layered reservoir and the growing hydraulic fracture is based on a well-established numerical technique called the Displacement Discontinuity Boundary Element (hereafter “DD”) Method. The method is extended to enable efficient and accurate calculation of the physical effects of layering in the reservoir by the use of a Fourier Transform Method, whereby the relations between stress and strain in the layered reservoir are determined. The numerical method assumes that each hydraulic fracture is divided into a regular mesh of rectangular elements, as depicted in FIG. 26, wherein each numerical element contains its own unique properties. Such properties include applied fluid and proppant pressure, fluid and proppant propagation direction and velocity, local reservoir properties, stress-strain relations, and fracture width.
FIG. 26 shows a numerical mesh of elements subdividing the fracture surface for purposes of calculation. In cases where the numerical element coincides with the fracture edge or tip (see FIG. 27), certain additional information is uniquely defined for such elements. For example, such information may include the local velocity of propagation of the fracture tip, the special relationship between the fluid in the fracture and the surrounding layered reservoir, and how the fluid and reservoir physically interact with each other. This interaction is accounted for by means of special properties assigned to the tip elements, comprising the interaction between a fluid-filled fracture and the host material it is fracturing.
FIG. 27 reveals a fracture outline on a numerical mesh. Each numerical element depicted in FIG. 26 or 27 relates to every other element in the numerical mesh by means of special mathematical relationships. We refer to elements as: (1) sending or source elements, and (2) receiver elements. A source element sends a signal representing a mathematical relationship to a receiver element. The signal is the net pressure, which is equal to the fluid pressure (p) 110 minus the confining stress (σc) 112 as indicated by the multi-layer elasticity equation 104 in FIG. 20, in that portion of the fracture. The receiver signal comprises the stress and strain experienced at the receiver location due to the applied pressure at the source element location. Many of these signals between source and receiver element are duplicated in the numerical mesh, and in these cases, special algorithms are designed to dramatically minimize the volume of storage required, so that only unique signals between different elements need to be stored.
The signals between each unique pair of receiver and source elements are stored in the computer memory or on a physical storage device in a matrix. The hydraulic fracture propagation numerical method is designed so that the fracture propagates in a finite number of time steps. At each time step, the signal matrix is invoked, extracting those signals which are active over the part of the numerical mesh that is covered by the current configuration of the hydraulic fracture. This matrix is then used to build a system of numerical equations that are solved for the fracture width at the current time—at each active element location.
During each fracture propagation step, another matrix of signals is constructed, the matrix comprising the physical behavior of the fluid in the hydraulic fracture, which relates the local fluid pressure to the local fracture width. This system of equations is also solved iteratively for local fluid pressures at each time step.
The combined system of equations must be coupled together in an efficient manner, so that they feed off each other until a balanced solution of fluid pressure and fracture width is obtained at each time step. This coupling between the two equation systems is accomplished by means of a special numerical algorithm that efficiently and accurately ensures that the correct solution is obtained. The entire system is designed to ensure that no fluid or proppant is unaccounted for in any time step.
The above process is repeated at each time step, thereby allowing the calculation of the way in which the fracture grows as a function of time. At each time step, the algorithm predicts which elements are active (i.e., filled with fluid and proppant), and the fracture width and fracture pressure on each active element. A complete description of the process of the propagation of a hydraulic fracture is thus obtained.
Solutions of the multi-layer equilibrium equations are provided.
In the following, we assume a three-dimensional body, so subscripts range from 1 to 3. The theory also applies to the two-dimensional cases (plane strain, plane stress, antiplane strain). The method provides an efficient way of determining the solution to the equilibrium equations:
σij,j l +b i l=0   (1)
for an, in general, transversely isotropic elastic layered medium, where superscript l indicates the layer number, with a stress-strain relationship given by:
σij l =C ijkr lεkr l  (2)
In the case of a transversely isotropic three-dimensional elastic medium, there are five independent material constants. The strain components in (2) are given by
ɛ kr l = 1 2 ( u k , r l + u r , k l ) ( 3 )
For a medium comprising multiple parallel layers each of which is homogeneous (see FIG. 28), it is possible to obtain a solution to the governing equations (1)-(3) by means of the Fourier Transform.
Refer now to FIG. 28 representing a schematic showing multiple parallel layers in a three dimensional case.
By substituting (3) and (2) into (1) and by taking the two-dimensional Fourier Transform with respect to x and z (where subscripts 1=x, 2=y, 3=z):
u ^ j l ( m , n , y ) = - - ( mx + nz ) u j l ( x , y , z ) x z ( 4 )
of the resulting equilibrium equations in terms of the displacements, we obtain a system of ordinary differential equations in the independent variable y. This system of ordinary differential equations determines the Fourier Transforms of displacement components ux l, uy l, and uz l:
L ( C ijkr l ) [ u ^ x u ^ y u ^ z ] l = [ b ^ x b ^ y b ^ z ] l ( 5 )
For a layered material, there is a system of differential equations of the form (5) for each layer, each of whose coefficients is determined by the material properties of the layer. It is possible to solve the system of differential equations for a typical layer l to obtain the general solution to the r th displacement components in the form:
u ^ r l = j d jr l α j l ky A j l ( k ) where k = m 2 + n 2 ( 6 )
In the case of repeated roots of the characteristic equation associated with (5), which occurs for the important case of isotropic layers, the system (5) has the general solution:
u ^ r l = j ( d jr l + f jr l y ) α j l ky A j l ( k ) ( 7 )
Here djr l and fjr l are constants that depend on the material constants of the layer, the αj l are the roots of the characteristic equation for the system of ordinary differential equations, and the Aj l(k) are free parameters of the solution that are determined by the forcing terms bi l in (1) and the interface conditions prescribed at the boundary between each of the layers (e.g., bonded, frictionless, etc.).
Substituting these displacement components into the stress strain law (2), we can obtain the corresponding stress components: {circumflex over (σ)}xx l, {circumflex over (σ)}yy l, {circumflex over (σ)}zz l, {circumflex over (σ)}xy l, {circumflex over (σ)}xz l, and {circumflex over (σ)}yz l, which can be expressed in the form:
σ ^ pq l = j ( s jpq l ) α j l ky A j l ( k ) ( 8 )
In the case of repeated roots the stress components assume the form:
σ ^ pq l = j ( s jpq l + t jpq l ky ) α j l ky A j l ( k ) ( 9 )
For each layer and for each sending DD element located at a particular y coordinate, there are a set of six parameters Aj l(k) that need to be determined from a system of algebraic equations that express the required body forces and boundary conditions in the model. Once the Aj l(k) have been determined, it is possible to calculate the influences of any DD having the same y component on any receiving point in any layer by taking the inverse Fourier Transform:
u j l = 1 ( 2 π ) 2 - - - ( mx + nz ) u ^ j l m n ( 10 )
of equations (6)-(9).
One of the major computational bottlenecks in the procedure outlined above is the inversion process represented by (10), which involves the numerical inversion of a two-dimensional Fourier Transform for each sending-receiving pair of DD elements. The scheme we propose uses an exponential approximation of the spectral solution coefficients Aj l(k) (see for example reference [11], which was constructed only for point sources in layered dielectric materials, or reference [12]—only for horizontal conducting elements in layered dielectric materials) of the form:
A j l ( k ) - A j l ( ) r a jr l b jr l k ( 11 )
Here Aj l(∞) represents the “high frequency” components of the spectrum of the solution, which represents the singular part of the solution in real space.
The inversion process can be achieved by evaluating integrals of the form:
I j pl = 1 ( 2 π ) 2 - - [ A j l ( k ) - A j l ( ) ] k p α j l ky - ( mx + nz ) m n or I j pl 1 ( 2 π ) 2 j a jr l - - k p k ( α j l y + b jr l ) - ( mx + nz ) m n ( 12 )
which can be evaluated in a closed form. The shifted components αj ly+bjr l in (12) represent a finite number of complex images that approximate what would be an L-fold infinite Fourier Series (for L layers) that would be required to represent the Green's function in a closed form using the method of images. Typically three or four complex images suffice to give a high order of accuracy.
The expressions of the form (12) are not much more complicated than the pair-wise DD influences that apply for a homogeneous elastic medium. The only difference in this case is that for each sending DD element, the expansion coefficients ajr l and bjr l for each layer need to be determined by solving the appropriate set of algebraic equations and performing the exponential fit (11) of these coefficients. Once these coefficients have been determined we have a very efficient procedure to determine the influences between DD elements.
Refer now to FIG. 26 representing a Schematic showing a duplicity of influence coefficients for a multiple layer problem.
In FIG. 26, for a regular array of DD elements there is an additional saving that can be exploited. In this case only the influence of a single sending DD element at each horizon (i.e., y level) needs to be determined in order to determine the whole influence matrix. For example, the DD elements in layer N denoted by the solid circle 202, the shaded circle 204, the unshaded circle 206, and the arrow 208, in FIG. 26, each have the same set of exponential expansion coefficients ajr l and bjr l in layer N-1, where the arrow 208 indicates where the expansion coefficients ajr l and bjr l are identical.
A DD influence at a specified point within a given layer is constructed by constructing a pseudo interface parallel to the layering across which there can be a jump in the displacement field. To construct a normal DD a jump in uy is specified, whereas to construct a shear or ride DD a jump in ux or uz is specified. This technique limits the orientation of DD components to be aligned with the pseudo interface that is parallel to the layering.
However, it is desirable to have DD components that specify jumps in the displacement field which are across interfaces that are not parallel to the layering in the material. In particular, for hydraulic fracturing problems in the petroleum industry, it is important to be able to model vertical fracture planes that are perpendicular to the horizontally layered material. In this case, and for arbitrarily oriented DDs, it is possible to construct a DD field of a desired orientation, by utilizing the duality relationship between the stresses due to a force discontinuity (or point force) and the displacements due to a displacement discontinuity. The solution to a force discontinuity in the r th direction can be constructed by taking br=δ(x, y, z)Fr, where δ(x, y, z) represents the Dirac delta function.
Having obtained the stresses due to a force discontinuity:
{circumflex over (σ)}ij={circumflex over (γ)}ijr F r   (13)
it is possible to determine the displacements due to a DD according to the following duality relation:
ûr={circumflex over (γ)}ijr D ij   (14)
The key idea here is to construct a planar Green's function or influence matrix, which represents the influences of all the DDs that lie in a vertical fracture plane. The influence matrix will only represent the mutual influences on one another of DDs that lie in the fracture plane. However, it will implicitly contain the information about the variations in material properties due to the layering.
A reduced influence matrix can be constructed by any numerical method, including that proposed above, which can represent rigorously the changes in material properties between the layers, for example the finite element method or a boundary integral method in which elements are placed on the interfaces between material layers. The in-plane influence coefficients would be calculated by means of the following procedure.
To calculate the influence of the ij th in-plane DD on the kl th DD anywhere else on the fracture plane, the value of the ij th DD would be assigned a value of unity and all the other fracture plane DDs would be set equal to zero. The boundary integral or finite element solution on the interfaces between the material layers would then be determined so as to ensure that there is compatibility in the displacements between the material layers as well as a balance in the forces between the layers at the interface. Once numerical solution has been calculated for the whole problem, the corresponding stresses on each of the in-plane DD elements can be evaluated to determine the in-plane stress influence of that unit DD on all the other DDs in the fracture plane. By repeating this process for each of the DDs that lie in the fracture plane, it is possible to determine the in-plane influence matrix which represents the effect that each DD in the plane has on any of the other in-plane DDs. By allowing the interface solution values to react to the sending DD element, the effect of the layers has been incorporated implicitly into this abbreviated set of influence coefficients.
The numerical procedure outlined to construct the desired in-plane influence matrix would take a considerable amount of time to compute. Indeed, such a process would likely exclude the possibility of real time processing with current personal computers or workstations, but could conceivably be performed in a batch mode prior to the desired simulation being performed. The semi-analytic method outlined above would be much more efficient, as the fully three-dimensional (or two-dimensional, in the case of plane strain or plane stress), problem that needs to be solved to calculate the influence of each DD element has been effectively reduced to a one-dimensional one.
Numerical models for multi-layered materials require that the interface between each material type is numerically “stitched” together by means of elements. For example, a boundary element method implementation would require that each interface between different material types is discretized into elements. A finite element or finite difference method implementation would require that the entire domain is discretized into elements. In the inventive method, the material interfaces are indirectly taken into account without the requirement of explicitly including elements away from the surface of the crack or fracture. The implication thereof, is a dramatic reduction in the number of equations to be solved, with a commensurate dramatic decrease in computer processing time. In addition, accuracy of the solution is maintained. One aspect of this invention which distinguishes it from previous work is that it is capable of solving problems in multi-layered elastic materials with arbitrarily inclined multiple cracks or fractures in two or three-dimensional space.
In this specification, note that ‘elements’ can intersect layers. This is accomplished by taking special care of the mathematical stress/strain relationships across interfaces in such a way as to obtain the correct physical response for the element which is located across the interface(s).
References 1-3 below are classic papers that establish the Fourier scheme to solve elastic multi-layer problems, but do not utilize the inversion scheme proposed here. In references 1 and 2, a propagator matrix approach is suggested to solve the system of algebraic equations necessary for the Fourier scheme, but this particular scheme will become unstable for problems with many layers.
References 4 and 5 use exponential approximation for inversion. The methods in those references do not give rise to the complex images generated by the algorithm presented in this invention that so efficiently represent the effect of many layers. Reference 5 extends the propagator approach used in references 1 and 2 to solve the algebraic equations of the Fourier method. Reference 5 discloses an inversion scheme that is an integral part of the propagator method. This method involves an exponential approximation, similar to that proposed in this patent, but it is applied to only one part of the propagator equations. As a result, a least squares fit of many terms (more than 50) is required to yield reasonable results using the teachings of this reference. Apart from stability issues involved with exponential fitting, the large number of terms would probably be less efficient than using direct numerical integration for inversion. The exponential fit of the spectral coefficients that we propose involves less than five terms.
References 6 through 10 extend the Fourier method to transversely isotropic media. References 7-10 use the propagator matrix for solving the algebraic equations, while reference 6 proposes direct solution. All these methods of solution would be numerically unstable for problems with many thick layers. While reference 10 proposes numerical inversion using continued fractions, little mention is made of the inversion process.
References 11 and 12 describe methodologies for multi-layer dielectric materials containing point electrical charges, or line charge distributions aligned parallel to interfaces (i.e., with different Green's functions to those used in elasticity).
Reference 13 describes a so-called “sweeping” algorithm to solve layered systems. The method disclosed in reference 13 is essentially the classic block LU decomposition for a block tri-diagonal system. In this specification, we use this algorithm to obtain a stable solution of the algebraic equations that determine the Fourier spectral coefficients in each of the layers. This method is particularly effective for problems in which the layers are thick or the wave-numbers are large.
It is recognized that other mathematical relationships may be used in the invention to achieve the same commercial or physical purpose. While not employing exactly the same equations, such methods are within the scope of the invention set forth in this specification.
The following references (i.e., references 1 through 15) are incorporated herein by reference:
  • 1. Ben-Menahem, A. and Singh, S. J. 1968. Multipolar elastic fields in a layered half space. Bull. Seism. Soc. Am. 58(5). 1,519-72.
  • 2. Singh, S. J. 1970. Static deformation of a multi-layered half-space by internal sources. J. Geophys. Res. 75(17). 3,257-63.
  • 3. Chen, W. T. 1971. Computation of stresses and displacements in a layered elastic medium. Int. J. Engng. Sci. vol. 9. 775-800.
  • 4. Sato, R. and Matsu'ura, M. 1973. Static deformations due to the fault spreading over several layers in a multi-layered medium Part I: Displacement. J. Phys. Earth. 21. 227-249.
  • 5. Jovanovich, D. B., Husseini, M. I. and Chinnery, M. A. 1974. Elastic dislocations in a layered half-space—I. Basic theory and numerical methods. Geophys. Jour. Roy. Astro. Soc. 39. 205-217.
  • 6. Wardle, L. J. 1980. Stress analysis of multilayered anisotropic elastic systems subject to rectangular loads. CSIRO Aust. Div. Appl. Geomech. Tech. paper no. 33. 1-37.
  • 7. Singh, S. J. 1986. Static deformation of a transversely isotropic multilayered half-space by surface loads. Physics of the Earth and Planetary Interiors. 42. 263-273.
  • 8. Pan, E. 1989. Static response of a transversely isotropic and layered half-space to general surface loads. Physics of the Earth and Planetary Interiors. 54. 353-363.
  • 9. Pan, E. 1989. Static response of a transversely isotropic and layered half-space to general dislocation sources. Physics of the Earth and Planetary Interiors. 58. 103-117.
  • 10. Pan, E. 1997. Static Green 's functions in multilayered half spaces. Appl. Math. Modeling. 21. 509-521.
  • 11. Chow, Y. L., Yang, J. J., and Howard, G. E. 1991. Complex images for electrostatic field computation in multilayered media. IEEE Trans. on Microwave Theory and Techniques. vol. 39. no. 7. 1120-25.
  • 12. Crampagne, R., Ahmadpanah, M. and Guiraud, J.-L. 1978. A simple method for determining the Green's function for a class of MIC lines having multilayered dielectric structures. IEEE Trans. on Microwave Theory and Techniques. vol. MTT-26. No. 2. 82-87.
  • 13. Linkov, A. M., Linkova, A. A., and Savitski, A. A. 1994. An effective method for multi-layered media with cracks and cavities. Int. J. of Damage Mech. 3. 338-35.
  • 14. Ryder, J. A., and Napier, J. A. L. 1985. Error Analysis and Design of a Large Scale Tabular Mining Stress Analyzer. Proceedings of the Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Japan, [Balkema] 1549-1555.
  • 15. J. A. Ryder, Eds.: E. G. Beer, J. R. Booker, and J. P. Carter, Optimal Iteration Schemes Suitable for General Non-linear Boundary Element Modeling Applications: Proceedings of the 7th International Conference on Computer Methods and Advances in Geomechanics, Cairns, Australia [Balkema], 1991.
The above description of the ‘Hydraulic Fracturing Simulator Software’ being thus described, it will be obvious that the same may be varied in many ways. Such variations are not to be regarded as a departure from the spirit and scope of the claimed method or apparatus or program storage device, and all such modifications as would be obvious to one skilled in the art are intended to be included within the scope of the following claims.

Claims (10)

1. A program storage device readable by a machine tangibly embodying a program of instructions executable by the machine to perform method steps comprising:
simulating a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir;
calculating a mesh overlaying the fracture thereby defining a plurality of elements;
calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, wherein each of the plurality of elements corresponds to at least one of the plurality of numbers, and wherein the plurality of numbers corresponds to a plurality of layers of the multilayered reservoir, wherein the calculating and determining an influence coefficient matrix comprises finding spectral coefficients in each layer of the multilayered reservoir, performing an exponential approximation, and assembling an influence coefficient matrix using exponential expansion coefficients, and wherein the finding spectral coefficients comprises determining the spectral coefficients by solving a system of algebraic equations that express the continuity of tractions and displacements across layer interfaces;
providing output data in response to the simulated hydraulic fracture to at least one of a display device and a recorder; and
wherein the system of algebraic equations comprise:
u ^ r l = j ( d jr l + f jr l y ) α j l ky A j l ( k ) , and σ ^ pq l = j ( s jpq l + t jpq l ky ) α j l ky A j l ( k ) ,
wherein ûr l defines the rth displacement component for the layer l, wherein p and q define stress tensor indicia, wherein {circumflex over (σ)}pq l defines a stress component for layer l in the direction pq, wherein Aj l(k) defines the spectral coefficient for column vector value j of the layer l, wherein y defines a spatial coordinate, wherein k defines a wave number in the Fourier-transform (FT) domain, wherein s and t define layer-dependent stress-related material property constants, wherein d and f define layer-dependent material property constants, and wherein αj lky defines a root of the characteristic equation for a system of ordinary differential equations (ODEs) relating directional stress to directional force for a layer.
2. The program storage device of claim 1, wherein the performing an exponential approximation comprises determining a set of free space spectral coefficients by solving the system of algebraic equations for an infinite homogeneous medium corresponding to a set of high frequency components associated with a prescribed displacement discontinuity, and evaluating a set of low frequency spectral components defined by

(A j l(k))LOW =A j l(k)−A j l(∞), and
approximating said low frequency spectral components by a series of exponential functions by solving for unknown constants ajr l and bjr l in the following expansion:
( A j l ( k ) ) LOW a jr l b jr l k ,
wherein (Aj l(k))LOW defines the low frequency spectral coefficient for column vector value j of layer l, wherein Aj l(∞) defines the high frequency spectral coefficient for column vector value j of layer l, and wherein ajr l and bjr l comprise exponential expansion coefficients for layer l.
3. The program storage device of claim 1, wherein the assembling an influence coefficient matrix using exponential expansion coefficients comprises:
substituting a set of expansions representing said low frequency spectral components into the system of algebraic equations to obtain a set of expressions for displacement and stress components, the set of expressions including,
u ^ p l ( k , y ) = j s jp r a jr l ( α j l + b jr l ) k , and σ ^ pq l ( k , y ) = j s jpq r a jr l ( α j l y + b jr l ) k ,
inverting the displacement and stress components associated with the low frequency spectral components in connection with a set of Fourier Transforms by evaluating integrals of the form,
I j p w l 1 ( 2 π ) 2 r a jr l - - k p w k ( α j l y + b jr l ) - ( mx + nz ) m n , and
combining the low frequency displacement and stress components with the infinite space displacements and stresses as follows

u i l=(u i l)LOW+(u i l) and σij l=(σij l)LOW+(σij l),
wherein Ij pl defines an integral equation for a power term pw of layer l, wherein ui l defines the ith displacement component for the layer l, wherein (ui l)LOW defines the ith low frequency displacement component for the layer l, wherein (ui l) defines the ith high frequency, displacement component for the layer l, wherein σij l defines the ith stress component for the layer l, wherein (σij l)LOW defines the ith low frequency stress component for the layer l, wherein (σij l) defines the ith high frequency stress component for the layer l, and wherein m and n define wave number pairs in the FT domain, and wherein x and z comprise spatial coordinates.
4. The program storage device of claim 1, wherein the method steps further comprise:
setting a time step to be a first time step;
in the first time step, iteratively calculating a fracture width in each element given a fluid pressure in said each element and iteratively calculating a fluid pressure in each element given a fracture width in each element until convergence occurs;
updating a fracture footprint in response to the iteratively calculating the fracture width;
updating the time step to a next time step;
repeating the setting a time step, the updating a fracture footprint, and the updating the time step until the time step reaches a maximum allowed time; and
generating the output data representing a fracture width and a fluid pressure in each element of the plurality of elements of the mesh.
5. A method comprising:
simulating using a processor a hydraulic fracture in an Earth formation where the formation includes a multilayered reservoir;
calculating a mesh overlaying the fracture thereby defining a plurality of elements;
calculating and determining an influence coefficient matrix, the matrix having a plurality of numbers, wherein each of the plurality of elements corresponds to at least one of the plurality of numbers, and wherein the plurality of numbers corresponds to a plurality of layers of the multilayered reservoir, wherein the calculating and determining an influence coefficient matrix comprises finding spectral coefficients in each layer of the multilayered reservoir, performing an exponential approximation, and assembling an influence coefficient matrix using exponential expansion coefficients, and wherein the finding spectral coefficients comprises determining the spectral coefficients by solving a system of algebraic equations that express the continuity of tractions and displacements across layer interfaces;
providing output data in response to the simulated hydraulic fracture to at least one of a display device and a recorder; and
wherein the system of algebraic equations comprise:
u ^ r l = j ( d jr l + f jr l y ) α j l ky A j l ( k ) , and σ ^ pq l = j ( s jpq l + t jpq l ky ) α j l ky A j l ( k ) ,
wherein ûr l defines the rth displacement component for the layer l, wherein p and q define stress tensor indicia, wherein {circumflex over (σ)}pq l defines a stress component for layer l in the direction pq, wherein Aj l(k) defines the spectral coefficient for column vector value j of the layer l, wherein y defines a spatial coordinate, wherein k defines a wave number in the Fourier-transform (FT) domain, wherein s and t define layer-dependent stress-related material property constants, wherein d and f define layer-dependent material property constants, and wherein αj lky defines a root of the characteristic equation for a system of ordinary differential equations (ODEs) relating directional stress to directional force for a layer.
6. The method of claim 5, wherein the performing an exponential approximation comprises determining a set of free space spectral coefficients by solving the system of algebraic equations for an infinite homogeneous medium corresponding to a set of high frequency components associated with a prescribed displacement discontinuity, evaluating a set of low frequency spectral components defined by

(A j l(k))LOW =A j l(k)−A j l(∞), and
approximating the low frequency spectral components by a series of exponential functions by solving for unknown constants ajr l and bjr l in the following expansion:
( A j l ( k ) ) LOW r a jr l b jr l k ,
wherein (Aj l(k))LOW defines the low frequency spectral coefficient for column vector value j of layer l, wherein Aj l(∞) defines the high frequency spectral coefficient for column vector value j of layer l, and wherein ajr l and bjr l comprise exponential expansion coefficients for layer l.
7. The method of claim 5, wherein the assembling an influence coefficient matrix using exponential expansion coefficients comprises substituting a set of expansions representing the low frequency spectral components into the system of algebraic equations to obtain a set of expressions for displacement and stress components, the set of expressions including
u ^ p l ( k , y ) = j s jp r a jr l ( α j l + b jr l ) k , and σ ^ pq l ( k , y ) = j s jpq r a jr l ( α j l y + b jr l ) k ,
inverting the displacement and stress components associated with the low frequency spectral components in connection with a set of Fourier Transforms by evaluating integrals of the form,
I j p w l 1 ( 2 π ) 2 r a jr l - - k p w k ( α j l y + b jr l ) - ( mx + nz ) m n , and
combining the low frequency displacement and stress components with the infinite space displacements and stresses as follows

u i l=(u i l)LOW+(u i l) and σij l=(σij l)LOW+(σij l),
wherein Ij pl defines an integral equation for a power term Pw of layer l, wherein ui l defines the ith displacement component for the layer l, wherein (ui l)LOW defines the ith low frequency displacement component for the layer l, wherein (ui l) defines the ith high frequency, displacement component for the layer l, wherein σij l defines the ith stress component for the layer l, wherein (σij l)LOW defines the ith low frequency stress component for the layer l, wherein (σij l) defines the ith high frequency stress component for the layer l, and wherein m and n define wave number pairs in the FT domain, and wherein x and z comprise spatial coordinates.
8. The method of claim 5, wherein the method steps further comprise:
setting a time step to be a first time step;
in the first time step, iteratively calculating a fracture width in each element given a fluid pressure in said each element and iteratively calculating a fluid pressure in each element given a fracture width in each element until convergence occurs;
updating a fracture footprint in response to the iteratively calculating the fracture width;
updating the time step to a next time step;
repeating the setting a time step, the updating a fracture footprint, and the updating the time step until the time step reaches a maximum allowed time; and
generating the output data representing a fracture width and a fluid pressure in each element of the plurality of elements of the mesh.
9. The method of claim 1, further comprising changing an injection schedule in response to the output data.
10. The method of claim 5, further comprising changing an injection schedule in response to the output data.
US11/095,104 1999-04-29 2005-03-31 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator Expired - Fee Related US7509245B2 (en)

Priority Applications (4)

Application Number Priority Date Filing Date Title
US11/095,104 US7509245B2 (en) 1999-04-29 2005-03-31 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
CA2541103A CA2541103C (en) 2005-03-31 2006-03-28 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
EA200600516A EA010456B1 (en) 2005-03-31 2006-03-30 Method, system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
US12/368,549 US8428923B2 (en) 1999-04-29 2009-02-10 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
US09/301,961 US6876959B1 (en) 1999-04-29 1999-04-29 Method and apparatus for hydraulic fractioning analysis and design
US10/831,799 US7063147B2 (en) 2004-04-26 2004-04-26 Method and apparatus and program storage device for front tracking in hydraulic fracturing simulators
US11/095,104 US7509245B2 (en) 1999-04-29 2005-03-31 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator

Related Parent Applications (2)

Application Number Title Priority Date Filing Date
US09/301,961 Continuation-In-Part US6876959B1 (en) 1999-04-29 1999-04-29 Method and apparatus for hydraulic fractioning analysis and design
US10/831,799 Continuation-In-Part US7063147B2 (en) 1999-04-29 2004-04-26 Method and apparatus and program storage device for front tracking in hydraulic fracturing simulators

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US09/301,961 Continuation US6876959B1 (en) 1999-04-29 1999-04-29 Method and apparatus for hydraulic fractioning analysis and design

Publications (2)

Publication Number Publication Date
US20050171751A1 US20050171751A1 (en) 2005-08-04
US7509245B2 true US7509245B2 (en) 2009-03-24

Family

ID=37055107

Family Applications (1)

Application Number Title Priority Date Filing Date
US11/095,104 Expired - Fee Related US7509245B2 (en) 1999-04-29 2005-03-31 Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator

Country Status (3)

Country Link
US (1) US7509245B2 (en)
CA (1) CA2541103C (en)
EA (1) EA010456B1 (en)

Cited By (29)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20080091396A1 (en) * 2006-10-13 2008-04-17 Kennon Stephen R Method and system for modeling and predicting hydraulic fracture performance in hydrocarbon reservoirs
US20090145599A1 (en) * 1999-04-29 2009-06-11 Eduard Siebrits Method System and Program Storage Device for Simulating A Multilayer Reservoir and Partially Active Elements In A Hydraulic Fracturing Simulator
US7689949B1 (en) * 2007-03-16 2010-03-30 Cadence Design Systems, Inc. Evaluating Green's functions for modeling multilayer integrated circuits
US20100224365A1 (en) * 2009-03-06 2010-09-09 Carlos Abad Method of treating a subterranean formation and forming treatment fluids using chemo-mathematical models and process control
US20110077918A1 (en) * 2009-09-25 2011-03-31 Ovunc Mutlu Method of Predicting Natural Fractures And Damage In A Subsurface Region
US8393393B2 (en) 2010-12-17 2013-03-12 Halliburton Energy Services, Inc. Coupler compliance tuning for mitigating shock produced by well perforating
US8397814B2 (en) 2010-12-17 2013-03-19 Halliburton Energy Serivces, Inc. Perforating string with bending shock de-coupler
US8397800B2 (en) 2010-12-17 2013-03-19 Halliburton Energy Services, Inc. Perforating string with longitudinal shock de-coupler
US20140067353A1 (en) * 2012-09-05 2014-03-06 Stratagen Wellbore completion and hydraulic fracturing optimization methods and associated systems
US8714251B2 (en) 2011-04-29 2014-05-06 Halliburton Energy Services, Inc. Shock load mitigation in a downhole perforation tool assembly
US8875796B2 (en) 2011-03-22 2014-11-04 Halliburton Energy Services, Inc. Well tool assemblies with quick connectors and shock mitigating capabilities
US8978817B2 (en) 2012-12-01 2015-03-17 Halliburton Energy Services, Inc. Protection of electronic devices used with perforating guns
US8978749B2 (en) 2012-09-19 2015-03-17 Halliburton Energy Services, Inc. Perforation gun string energy propagation management with tuned mass damper
US8985200B2 (en) 2010-12-17 2015-03-24 Halliburton Energy Services, Inc. Sensing shock during well perforating
US9091152B2 (en) 2011-08-31 2015-07-28 Halliburton Energy Services, Inc. Perforating gun with internal shock mitigation
US9297228B2 (en) 2012-04-03 2016-03-29 Halliburton Energy Services, Inc. Shock attenuator for gun system
US9598940B2 (en) 2012-09-19 2017-03-21 Halliburton Energy Services, Inc. Perforation gun string energy propagation management system and methods
US9618652B2 (en) 2011-11-04 2017-04-11 Schlumberger Technology Corporation Method of calibrating fracture geometry to microseismic events
US10151856B1 (en) * 2017-01-24 2018-12-11 Baker Hughes, A Ge Company, Llc Simulating a subterranean region using a finite element mesh and a boundary element mesh
US10352145B2 (en) 2011-03-11 2019-07-16 Schlumberger Technology Corporation Method of calibrating fracture geometry to microseismic events
WO2019143320A1 (en) * 2018-01-16 2019-07-25 Halliburton Energy Services, Inc. Modeling fracture closure processes in hydraulic fracturing simulators
US10422208B2 (en) 2011-11-04 2019-09-24 Schlumberger Technology Corporation Stacked height growth fracture modeling
US10544667B2 (en) 2011-11-04 2020-01-28 Schlumberger Technology Corporation Modeling of interaction of hydraulic fractures in complex fracture networks
US10572611B2 (en) * 2016-04-29 2020-02-25 Exxonmobil Upstream Research Company Method and system for characterizing fractures in a subsurface region
US10787887B2 (en) 2015-08-07 2020-09-29 Schlumberger Technology Corporation Method of performing integrated fracture and reservoir operations for multiple wellbores at a wellsite
US10794154B2 (en) 2015-08-07 2020-10-06 Schlumberger Technology Corporation Method of performing complex fracture operations at a wellsite having ledged fractures
US10920552B2 (en) 2015-09-03 2021-02-16 Schlumberger Technology Corporation Method of integrating fracture, production, and reservoir operations into geomechanical operations of a wellsite
US10920538B2 (en) 2015-08-07 2021-02-16 Schlumberger Technology Corporation Method integrating fracture and reservoir operations into geomechanical operations of a wellsite
US11578568B2 (en) 2015-08-07 2023-02-14 Schlumberger Technology Corporation Well management on cloud computing system

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7516056B2 (en) * 2005-04-26 2009-04-07 Schlumberger Technology Corporation Apparatus, method and system for improved reservoir simulation using a multiplicative overlapping Schwarz preconditioning for adaptive implicit linear systems
RU2486336C2 (en) * 2007-11-01 2013-06-27 Лоджинд Б.В. Method of formation breakdown simulation and its estimation, and computer-read carrier
US8175807B2 (en) * 2008-02-20 2012-05-08 Schlumberger Technology Corporation Estimating horizontal stress from three-dimensional anisotropy
US10060241B2 (en) * 2009-06-05 2018-08-28 Schlumberger Technology Corporation Method for performing wellbore fracture operations using fluid temperature predictions
US8498852B2 (en) * 2009-06-05 2013-07-30 Schlumberger Tehcnology Corporation Method and apparatus for efficient real-time characterization of hydraulic fractures and fracturing optimization based thereon
MX2012007687A (en) * 2009-12-30 2012-07-20 Schlumberger Technology Bv Method for controlling the trajectory of a hydraulic fracture in strata containing natural fractures.
US9366121B2 (en) * 2012-02-06 2016-06-14 Halliburton Energy Services, Inc. Modeling fracturing fluid leak-off
US9217318B2 (en) * 2013-03-14 2015-12-22 Halliburton Energy Services, Inc. Determining a target net treating pressure for a subterranean region
US9297250B2 (en) 2013-03-14 2016-03-29 Halliburton Energy Services, Inc. Controlling net treating pressure in a subterranean region
US9958572B2 (en) * 2015-03-31 2018-05-01 Halliburton Energy Services, Inc. Synthetic test beds for fracturing optimization and methods of manufacture and use thereof
CN105003254B (en) * 2015-08-05 2017-06-06 中国海洋石油总公司 HTHP pit shaft annular space confined space temperature and pressure change modeling experimental provision
CA2993293A1 (en) * 2015-09-08 2017-03-16 Halliburton Energy Services, Inc. Simulators and simulation methods using adaptive domains
CN108763841B (en) * 2018-07-24 2022-04-05 北京航空航天大学青岛研究院 Elastic fracture simulation method based on dual boundary element and strain energy optimization analysis
CN111322050B (en) * 2020-04-24 2022-02-11 西南石油大学 Shale horizontal well section internal osculating temporary plugging fracturing construction optimization method

Citations (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5310003A (en) * 1993-04-19 1994-05-10 Mobil Oil Corporation Method for acid fracturing a subterranean formation
US5442173A (en) * 1994-03-04 1995-08-15 Schlumberger Technology Corporation Method and system for real-time monitoring of earth formation fracture movement
US5771170A (en) * 1994-02-14 1998-06-23 Atlantic Richfield Company System and program for locating seismic events during earth fracture propagation
US6052520A (en) * 1998-02-10 2000-04-18 Exxon Production Research Company Process for predicting behavior of a subterranean formation
US6094619A (en) * 1997-07-04 2000-07-25 Institut Francais Du Petrole Method for determining large-scale representative hydraulic parameters of a fractured medium
US6106561A (en) * 1997-06-23 2000-08-22 Schlumberger Technology Corporation Simulation gridding method and apparatus including a structured areal gridder adapted for use by a reservoir simulator
US6192316B1 (en) * 1999-05-26 2001-02-20 Atlantic Richfield Company Fracture discrimination using borehole frequency response of stoneley waves
US6216786B1 (en) * 1998-06-08 2001-04-17 Atlantic Richfield Company Method for forming a fracture in a viscous oil, subterranean formation
US6230101B1 (en) * 1999-06-03 2001-05-08 Schlumberger Technology Corporation Simulation method and apparatus
US20020032494A1 (en) * 2000-06-29 2002-03-14 Kennon Stephen R. Feature modeling in a finite element model
US20020091502A1 (en) * 2000-04-04 2002-07-11 Anders Malthe-Sorenssen Method of modeling of faulting and fracturing in the earth
US6662146B1 (en) * 1998-11-25 2003-12-09 Landmark Graphics Corporation Methods for performing reservoir simulation
US20040016541A1 (en) * 2002-02-01 2004-01-29 Emmanuel Detournay Interpretation and design of hydraulic fracturing treatments
US6810370B1 (en) * 1999-03-31 2004-10-26 Exxonmobil Upstream Research Company Method for simulation characteristic of a physical system
US6826483B1 (en) * 1999-10-13 2004-11-30 The Trustees Of Columbia University In The City Of New York Petroleum reservoir simulation and characterization system and method
US6826520B1 (en) * 1999-06-24 2004-11-30 Exxonmobil Upstream Research Company Method of upscaling permeability for unstructured grids
US20040254730A1 (en) * 2002-01-31 2004-12-16 Najmuddin Ilyas Juzer Method and apparatus for detecting fractures using frequency data derived from seismic data
US6834233B2 (en) * 2002-02-08 2004-12-21 University Of Houston System and method for stress and stability related measurements in boreholes
US6842725B1 (en) * 1998-12-11 2005-01-11 Institut Francais Du Petrole Method for modelling fluid flows in a fractured multilayer porous medium and correlative interactions in a production well
US6876959B1 (en) * 1999-04-29 2005-04-05 Schlumberger Technology Corporation Method and apparatus for hydraulic fractioning analysis and design
US6883153B2 (en) * 2003-01-10 2005-04-19 Intel Corporation Minimization of microelectronic interconnect thickness variations
US20050115711A1 (en) * 2003-11-11 2005-06-02 Schlumberger Technology Corporation Method and system for determining an optimum pumping schedule corresponding to an optimum return on investment when fracturing a formation penetrated by a wellbore
US6938693B2 (en) * 2001-10-31 2005-09-06 Schlumberger Technology Corporation Methods for controlling screenouts
US20050222852A1 (en) * 2004-03-30 2005-10-06 Craig David P Method and an apparatus for detecting fracture with significant residual width from previous treatments
US7006959B1 (en) * 1999-10-12 2006-02-28 Exxonmobil Upstream Research Company Method and system for simulating a hydrocarbon-bearing formation
US7024342B1 (en) * 2000-07-01 2006-04-04 Mercury Marine Thermal flow simulation for casting/molding processes
US7096122B2 (en) * 2003-07-22 2006-08-22 Dianli Han Method for producing full field radial grid for hydrocarbon reservoir simulation
US7386431B2 (en) * 2005-03-31 2008-06-10 Schlumberger Technology Corporation Method system and program storage device for simulating interfacial slip in a hydraulic fracturing simulator software

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5870691A (en) * 1996-12-06 1999-02-09 Amoco Corporation Spectral decomposition for seismic interpretation
US6131071A (en) * 1996-12-06 2000-10-10 Bp Amoco Corporation Spectral decomposition for seismic interpretation
RU2119679C1 (en) * 1997-12-05 1998-09-27 Станислав Васильевич Васильев Method for searching and prospecting of oil and gas deposits

Patent Citations (31)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5310003A (en) * 1993-04-19 1994-05-10 Mobil Oil Corporation Method for acid fracturing a subterranean formation
US5771170A (en) * 1994-02-14 1998-06-23 Atlantic Richfield Company System and program for locating seismic events during earth fracture propagation
US5442173A (en) * 1994-03-04 1995-08-15 Schlumberger Technology Corporation Method and system for real-time monitoring of earth formation fracture movement
US6106561A (en) * 1997-06-23 2000-08-22 Schlumberger Technology Corporation Simulation gridding method and apparatus including a structured areal gridder adapted for use by a reservoir simulator
US6094619A (en) * 1997-07-04 2000-07-25 Institut Francais Du Petrole Method for determining large-scale representative hydraulic parameters of a fractured medium
US6052520A (en) * 1998-02-10 2000-04-18 Exxon Production Research Company Process for predicting behavior of a subterranean formation
US6216786B1 (en) * 1998-06-08 2001-04-17 Atlantic Richfield Company Method for forming a fracture in a viscous oil, subterranean formation
US6662146B1 (en) * 1998-11-25 2003-12-09 Landmark Graphics Corporation Methods for performing reservoir simulation
US6842725B1 (en) * 1998-12-11 2005-01-11 Institut Francais Du Petrole Method for modelling fluid flows in a fractured multilayer porous medium and correlative interactions in a production well
US6810370B1 (en) * 1999-03-31 2004-10-26 Exxonmobil Upstream Research Company Method for simulation characteristic of a physical system
US6876959B1 (en) * 1999-04-29 2005-04-05 Schlumberger Technology Corporation Method and apparatus for hydraulic fractioning analysis and design
US6192316B1 (en) * 1999-05-26 2001-02-20 Atlantic Richfield Company Fracture discrimination using borehole frequency response of stoneley waves
US6230101B1 (en) * 1999-06-03 2001-05-08 Schlumberger Technology Corporation Simulation method and apparatus
US6826520B1 (en) * 1999-06-24 2004-11-30 Exxonmobil Upstream Research Company Method of upscaling permeability for unstructured grids
US7006959B1 (en) * 1999-10-12 2006-02-28 Exxonmobil Upstream Research Company Method and system for simulating a hydrocarbon-bearing formation
US6826483B1 (en) * 1999-10-13 2004-11-30 The Trustees Of Columbia University In The City Of New York Petroleum reservoir simulation and characterization system and method
US20020091502A1 (en) * 2000-04-04 2002-07-11 Anders Malthe-Sorenssen Method of modeling of faulting and fracturing in the earth
US20020032494A1 (en) * 2000-06-29 2002-03-14 Kennon Stephen R. Feature modeling in a finite element model
US7260508B2 (en) * 2000-06-29 2007-08-21 Object Reservoir, Inc. Method and system for high-resolution modeling of a well bore in a hydrocarbon reservoir
US6674432B2 (en) * 2000-06-29 2004-01-06 Object Reservoir, Inc. Method and system for modeling geological structures using an unstructured four-dimensional mesh
US7024342B1 (en) * 2000-07-01 2006-04-04 Mercury Marine Thermal flow simulation for casting/molding processes
US6938693B2 (en) * 2001-10-31 2005-09-06 Schlumberger Technology Corporation Methods for controlling screenouts
US20040254730A1 (en) * 2002-01-31 2004-12-16 Najmuddin Ilyas Juzer Method and apparatus for detecting fractures using frequency data derived from seismic data
US20040016541A1 (en) * 2002-02-01 2004-01-29 Emmanuel Detournay Interpretation and design of hydraulic fracturing treatments
US6834233B2 (en) * 2002-02-08 2004-12-21 University Of Houston System and method for stress and stability related measurements in boreholes
US7006918B2 (en) * 2002-02-08 2006-02-28 University Of Houston Method for stress and stability related measurements in boreholes
US6883153B2 (en) * 2003-01-10 2005-04-19 Intel Corporation Minimization of microelectronic interconnect thickness variations
US7096122B2 (en) * 2003-07-22 2006-08-22 Dianli Han Method for producing full field radial grid for hydrocarbon reservoir simulation
US20050115711A1 (en) * 2003-11-11 2005-06-02 Schlumberger Technology Corporation Method and system for determining an optimum pumping schedule corresponding to an optimum return on investment when fracturing a formation penetrated by a wellbore
US20050222852A1 (en) * 2004-03-30 2005-10-06 Craig David P Method and an apparatus for detecting fracture with significant residual width from previous treatments
US7386431B2 (en) * 2005-03-31 2008-06-10 Schlumberger Technology Corporation Method system and program storage device for simulating interfacial slip in a hydraulic fracturing simulator software

Non-Patent Citations (39)

* Cited by examiner, † Cited by third party
Title
A. A. Linkova and A. A. Savitski, An Effective Method for Multi-Layered Media with Cracks and Cavities, p. 338-355, (1994) International Journal of Damage Mechanics, vol. 3.
A. Peirce, and E. Siebrits, "Stability Analysis and Design of Time-Stepping Schemes for General Elastodynamic Boundary Element Models", International Journal for Numerical Methods in Engineering, vol. 40, pp. 319-342, 1997. *
A. Peirce, E. Siebrits, "Stability Analysis and Design of Time-Stepping Schemes for General Elastodynamic Boundary elements Models,", International Journal for Numerical Method in Engineering, vol. 40, 319-342, 1997. *
Ben-Menahem, A. and Singh, S.J. 1968. Multipolar Elastic Fields In A Layered Half Space, Bull. Seism. Soc. Am. 58(5). 1, 519-72.
Cade Office, Solution Revolution Break Open Your Reserves with FracCADE, p. 1-6, (Jan. 1998), Schlumberger Oilfield Services, Dowell, GeoQuest.
Chen. W.T. 1971. Computation Of Stresses And Displacements In A Layered Elastic Medium. Int. J. Engng Sci. vol. 9. 775-800.
Chow, Y.L., Yang, J.J. and Howard, G.E. 1991. Complex Images For Electrostatic Field Computation In Multilayered Media. IEEE Trans. On Microwave Theory and Techniques. vol. 39.No. 7. 1120-25.
Crampagne, R. Ahmadpanah, M and Guiraud, J.L. 1978. A Simple Method For Determining The Green's Function For A Class of MIC Lines Having Multilayered Dielectric Structures. IEEE Trans. On Microwave Theory and Techniques. vol. MTT-26.No. 2. 82-87.
D. A. Mendelsohn, A Review of Hydraulic Fracture Modeling-Part I: General Concepts, 2D Models, Motivation for 3D Modeling, p. 369-376, (Sep. 1984), Journal of Energy Resources Technology.
Dr. Quan Guo, A Fully Three-Dimensional Hydraulic Fracturing Simulator, p. 1-14, TerraFrac (Jan. 1, 1995).
Dushan B. Jovanovic, Moujahed I. Husseini and Michael A. Chinnery, Elastic Dislocations in a Layered Half-Space-I. Basics Theory and Numerical Methods, p. 205-217, (1974), Geophys. J. R. astr. Soc. 39.
Ernian Pan, Static Green's Functions in MultiLayered HalfSpaces, p. 509-521, (1997), Appl. Math. Modeling, vol. 21.
Ernian Pan, Static Response of a Transversley Isotropic and Layered Half-Space to General Dislocation Sources, p. 103-115, (1989), Physics of the Earth and Planetary Interiors, 58.
FracCade Software, p. 1-4, Schlumberger Dowell, (Jan. 1, 1995).
Gohfer Grid Oriented Hydraulic Fracture Extension Replicator, p. 1-52, Appendix, p. 1-11, (1996) by Stim-Lab, Inc. and Marathon Oil Company.
J. A. Ryder and J. A. L. Napier, Error Analysis and Design of a Large-Scale Tabular Mining Stress Analyser, (1985) Chamber of Mines of South Africa Research Organization.
J.A. Ryder, Cairns, E.G. Beer, J.R. Booker, and J.P. Carter, Optimal Iteration Schemes Suitable For General Non-Linear Boundary Element Modeling Applications: Proceedings for the 7th International Conference on Computer Methods and Advances in Geomechanics, Balkema Rotterdam, 1991.
J.A. Ryder, Optimal iteration schemes suitable for general non-linear boundary element modeling applications, pp. 1079-1084, (1991), Rock Engineering Division, COMRO, South Africa.
Jovanovich, D.B. Husseini, M.I. and Chinnery, M.A. 1974. Elastic Dislocations In A Layered Half-Space-I. Basic Theory And Numerical Methods. Geophys. J.R., Astro. Soc. 39.205-217.
L. J. Wardle, Stress Analysis of Multilayered Anisotropic Elastic Systems Subject to Rectangular Loads, p. 1-37, (1980), CSIRO.
Linkov, A.M., Linkova, A.A. and Savitski, A.A. 1994. An Effective Method For Multi-layered Media With Cracks and Cavities. Int. J. of Damage Mech. 3. 338-35.
Pan, E. 1989, Static Response Of A Transversely Isotropic And Layered Half-Space To General Surface Lodas. Physic of the Earth and Planetary Interiors. 54. 353-363.
Pan, E. 1989. Static Response Of A Transverseley Istropic And Layered Half-Space To General Dislocation Sources. Physics of the Earth and Planetary Interiors,58. 103-117.
Pan, E. 1997. Static Green's Functions In Multilayered Half Spaces. Appl. Math. Modeling. 21, 509-521.
R.J. Clifton, Three-Dimensional Fracture-Propagation Models, p. 95-108, Brown U. and Terra Tek Inc. (Jan. 1, 1989).
Raymond Crampagne, Majid Ahmadpanah, and Jean-Louis Guiraud, A Simple Method for Determining the Green's Function for a Large Class of MIC Lines Having Multilayered Dielectric Structures, p. 82-87 (1978) IEEE.
Ryder, J.A. and Napier, J.A. L. 1985. Error Analysis and Design Of A Large Scale Tabular Mining Stress Analyzer. Proc: Fifith Int. Conf. On Num. Meth. In Geomech. Nagoya. 1549-1555.
Ryosuke Sato and Mitsuhiro Matsu'Ural, Static Deformations Due to the Fault Spreading over Several Layers in a Multi-Layered Medium Part I: Displacement, p. 227-249, (Sep. 20, 1973), J. Phys. Earth, 21.
Salamon M., "Some Application of the Frictionless Laminated Model". Rock Mechanics as a Guide for Efficient Utilization of Natural Resources, Khair (ed.), copyright 1989, Balkema, Rotterdam, ISBN 90 6191 871 5.
Salamon M., "Subsidence Prediction Using a Laminated Linear Model", Rock Mechanics as a Guide for Efficient Utilization of Natural Resources, Khair (ed.), copyright 1989, Balkema Rotterdam, ISBN 90 6191 8715.
Sarva Jit Singh, Static Deformation of a Multilayered Half-Space by Internal Sources, p. 3257-3263. (Jun. 10, 1970), Journal of Geophysical Research.
Sarva Jit Singh, Static Deformation of Transversely Isotropic Multilayered Half-Space by Surface Loads, p. 263-273, (1986), Physics of the Earth and Planetary Interiors, 42.
Sato, R. and Matsu'ura, M. 1973. Static Deformations Due To The Fault Spreading Over Several Layers In A Multi-Layered Medium Part I: Displacement. J.Phys. Earth. 21.
Singh, S.J. 1970. Static Deformation Of A Multi-Layered Half-Space By Internal Sources. J. Geophys. Res. 75(17). 3, 257-63.
Singh, S.J. 1986. Static Deformations Of A Transversely Isotropic Multilayered Half-Space By Surface Loads. Physics of the Earth and Planetary Interiors. 42. 263-273.
W. T. Chen, Computation of Stresses and Displacements in a Layered Elastic Medium, p. 775-800 (1971) Int. J. Engng Sci. vol. 9.
Wardle, L.J. 1980. Stress Analysis Of Multilayered Anisotropic Elastic Systems Subject To Rectangular Loads, CSIRO Aust. Div. Appl. Geomech. Tech. paper No. 33. 1-37.
Y. Leonard Chow, Jian Jun Yang, and Gregory E. Howard, Complex Images for Electrostatic Field Computation in Multilayered Media, p. 1120-1125, (1991) IEEE.
Yue, Z. and Yin J., "Backward Transfer-Matrix Method for Elastic Analysis of Layered Solids with Imperfect Bonding", Journal of Elasticity, 50: 109-128, 1998 copyright 1998 Kluwer Academic Publishers.

Cited By (43)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8428923B2 (en) * 1999-04-29 2013-04-23 Schlumberger Technology Corporation Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
US20090145599A1 (en) * 1999-04-29 2009-06-11 Eduard Siebrits Method System and Program Storage Device for Simulating A Multilayer Reservoir and Partially Active Elements In A Hydraulic Fracturing Simulator
US7925482B2 (en) * 2006-10-13 2011-04-12 Object Reservoir, Inc. Method and system for modeling and predicting hydraulic fracture performance in hydrocarbon reservoirs
US20080091396A1 (en) * 2006-10-13 2008-04-17 Kennon Stephen R Method and system for modeling and predicting hydraulic fracture performance in hydrocarbon reservoirs
US7689949B1 (en) * 2007-03-16 2010-03-30 Cadence Design Systems, Inc. Evaluating Green's functions for modeling multilayer integrated circuits
US20100224365A1 (en) * 2009-03-06 2010-09-09 Carlos Abad Method of treating a subterranean formation and forming treatment fluids using chemo-mathematical models and process control
US9085975B2 (en) 2009-03-06 2015-07-21 Schlumberger Technology Corporation Method of treating a subterranean formation and forming treatment fluids using chemo-mathematical models and process control
US20110077918A1 (en) * 2009-09-25 2011-03-31 Ovunc Mutlu Method of Predicting Natural Fractures And Damage In A Subsurface Region
US8494827B2 (en) 2009-09-25 2013-07-23 Exxonmobil Upstream Research Company Method of predicting natural fractures and damage in a subsurface region
US8985200B2 (en) 2010-12-17 2015-03-24 Halliburton Energy Services, Inc. Sensing shock during well perforating
US8408286B2 (en) 2010-12-17 2013-04-02 Halliburton Energy Services, Inc. Perforating string with longitudinal shock de-coupler
US8490686B2 (en) 2010-12-17 2013-07-23 Halliburton Energy Services, Inc. Coupler compliance tuning for mitigating shock produced by well perforating
US8397800B2 (en) 2010-12-17 2013-03-19 Halliburton Energy Services, Inc. Perforating string with longitudinal shock de-coupler
US8393393B2 (en) 2010-12-17 2013-03-12 Halliburton Energy Services, Inc. Coupler compliance tuning for mitigating shock produced by well perforating
US8397814B2 (en) 2010-12-17 2013-03-19 Halliburton Energy Serivces, Inc. Perforating string with bending shock de-coupler
US10352145B2 (en) 2011-03-11 2019-07-16 Schlumberger Technology Corporation Method of calibrating fracture geometry to microseismic events
US9206675B2 (en) 2011-03-22 2015-12-08 Halliburton Energy Services, Inc Well tool assemblies with quick connectors and shock mitigating capabilities
US8875796B2 (en) 2011-03-22 2014-11-04 Halliburton Energy Services, Inc. Well tool assemblies with quick connectors and shock mitigating capabilities
US8714252B2 (en) 2011-04-29 2014-05-06 Halliburton Energy Services, Inc. Shock load mitigation in a downhole perforation tool assembly
US8881816B2 (en) 2011-04-29 2014-11-11 Halliburton Energy Services, Inc. Shock load mitigation in a downhole perforation tool assembly
US8714251B2 (en) 2011-04-29 2014-05-06 Halliburton Energy Services, Inc. Shock load mitigation in a downhole perforation tool assembly
US9091152B2 (en) 2011-08-31 2015-07-28 Halliburton Energy Services, Inc. Perforating gun with internal shock mitigation
US10544667B2 (en) 2011-11-04 2020-01-28 Schlumberger Technology Corporation Modeling of interaction of hydraulic fractures in complex fracture networks
US9618652B2 (en) 2011-11-04 2017-04-11 Schlumberger Technology Corporation Method of calibrating fracture geometry to microseismic events
US10422208B2 (en) 2011-11-04 2019-09-24 Schlumberger Technology Corporation Stacked height growth fracture modeling
US9297228B2 (en) 2012-04-03 2016-03-29 Halliburton Energy Services, Inc. Shock attenuator for gun system
US20140067353A1 (en) * 2012-09-05 2014-03-06 Stratagen Wellbore completion and hydraulic fracturing optimization methods and associated systems
US9262713B2 (en) * 2012-09-05 2016-02-16 Carbo Ceramics Inc. Wellbore completion and hydraulic fracturing optimization methods and associated systems
US8978749B2 (en) 2012-09-19 2015-03-17 Halliburton Energy Services, Inc. Perforation gun string energy propagation management with tuned mass damper
US9598940B2 (en) 2012-09-19 2017-03-21 Halliburton Energy Services, Inc. Perforation gun string energy propagation management system and methods
US8978817B2 (en) 2012-12-01 2015-03-17 Halliburton Energy Services, Inc. Protection of electronic devices used with perforating guns
US9926777B2 (en) 2012-12-01 2018-03-27 Halliburton Energy Services, Inc. Protection of electronic devices used with perforating guns
US9909408B2 (en) 2012-12-01 2018-03-06 Halliburton Energy Service, Inc. Protection of electronic devices used with perforating guns
US9447678B2 (en) 2012-12-01 2016-09-20 Halliburton Energy Services, Inc. Protection of electronic devices used with perforating guns
US10787887B2 (en) 2015-08-07 2020-09-29 Schlumberger Technology Corporation Method of performing integrated fracture and reservoir operations for multiple wellbores at a wellsite
US10794154B2 (en) 2015-08-07 2020-10-06 Schlumberger Technology Corporation Method of performing complex fracture operations at a wellsite having ledged fractures
US10920538B2 (en) 2015-08-07 2021-02-16 Schlumberger Technology Corporation Method integrating fracture and reservoir operations into geomechanical operations of a wellsite
US11578568B2 (en) 2015-08-07 2023-02-14 Schlumberger Technology Corporation Well management on cloud computing system
US10920552B2 (en) 2015-09-03 2021-02-16 Schlumberger Technology Corporation Method of integrating fracture, production, and reservoir operations into geomechanical operations of a wellsite
US10572611B2 (en) * 2016-04-29 2020-02-25 Exxonmobil Upstream Research Company Method and system for characterizing fractures in a subsurface region
US10151856B1 (en) * 2017-01-24 2018-12-11 Baker Hughes, A Ge Company, Llc Simulating a subterranean region using a finite element mesh and a boundary element mesh
WO2019143320A1 (en) * 2018-01-16 2019-07-25 Halliburton Energy Services, Inc. Modeling fracture closure processes in hydraulic fracturing simulators
US20200401739A1 (en) * 2018-01-16 2020-12-24 Halliburton Energy Services, Inc. Modeling Fracture Closure Processes In Hydraulic Fracturing Simulators

Also Published As

Publication number Publication date
CA2541103A1 (en) 2006-09-30
EA200600516A1 (en) 2006-10-27
US20050171751A1 (en) 2005-08-04
EA010456B1 (en) 2008-08-29
CA2541103C (en) 2016-07-12

Similar Documents

Publication Publication Date Title
US7509245B2 (en) Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
US8428923B2 (en) Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
US6876959B1 (en) Method and apparatus for hydraulic fractioning analysis and design
EP3271547B1 (en) Stacked height growth fracture modeling
US7386431B2 (en) Method system and program storage device for simulating interfacial slip in a hydraulic fracturing simulator software
EP2774066B1 (en) Modeling of interaction of hydraulic fractures in complex fracture networks
Nagel et al. Coupled numerical evaluations of the geomechanical interactions between a hydraulic fracture stimulation and a natural fracture system in shale formations
Kresse et al. Numerical modeling of hydraulic fractures interaction in complex naturally fractured formations
CN105074125B (en) The method calibrated for microseismic event fracture geometry
US20160265331A1 (en) Modeling of interaction of hydraulic fractures in complex fracture networks
AU2019200654A1 (en) Modeling of interaction of hydraulic fractures in complex fracture networks
US10083254B2 (en) Method and system for stabilizing formulation methods
CN113821953B (en) Temporary plugging repeated fracturing design method based on dynamic ground stress and residual oil characteristics
Zhang et al. Hydraulic fracturing induced fault slip and casing shear in Sichuan Basin: A multi-scale numerical investigation
US7063147B2 (en) Method and apparatus and program storage device for front tracking in hydraulic fracturing simulators
Mahrer Microseismic logging: a new hydraulic fracture diagnostic method
McDaniel A Review of Design Considerations for Fracture Stimulation of Highly Deviated Wellbores
Ji et al. Modeling hydraulic fracturing fully coupled with reservoir and geomechanical simulation
Azim et al. Modeling hydraulic fractures in finite difference simulators using amalgam local grid refinement (LGR)
Yu et al. Advances in geomechanical modeling
Huang et al. Application of dilation-recompaction model in hydraulic fracturing simulation
MXPA06003481A (en) Method system and program storage device for simulating a multilayer reservoir and partially active elements in a hydraulic fracturing simulator
Nandlal High-Resolution Modeling and Visualization of Drained Rock Volume around Hydraulically Fractured Wells in Heterogeneous Reservoirs

Legal Events

Date Code Title Description
AS Assignment

Owner name: SCHLUMBERGER TECHNOLOGY CORPORATION, TEXAS

Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:SIEBRITS, EDUARD;PEIRCE, ANTHONY;REEL/FRAME:016247/0809;SIGNING DATES FROM 20050330 TO 20050411

STCF Information on status: patent grant

Free format text: PATENTED CASE

FPAY Fee payment

Year of fee payment: 4

FPAY Fee payment

Year of fee payment: 8

FEPP Fee payment procedure

Free format text: MAINTENANCE FEE REMINDER MAILED (ORIGINAL EVENT CODE: REM.); ENTITY STATUS OF PATENT OWNER: LARGE ENTITY

LAPS Lapse for failure to pay maintenance fees

Free format text: PATENT EXPIRED FOR FAILURE TO PAY MAINTENANCE FEES (ORIGINAL EVENT CODE: EXP.); ENTITY STATUS OF PATENT OWNER: LARGE ENTITY

STCH Information on status: patent discontinuation

Free format text: PATENT EXPIRED DUE TO NONPAYMENT OF MAINTENANCE FEES UNDER 37 CFR 1.362

FP Lapsed due to failure to pay maintenance fee

Effective date: 20210324