当前位置：文档库 › Hydraulic_Fracturing_Simulation

Hydraulic_Fracturing_Simulation

Development and Validation of Fully-Coupled

Hydraulic Fracturing Simulation Capabilities

Matias G. Zielonka, Kevin H. Searles, Jing Ning and Scott R. Buechler

ExxonMobil Upstream Research Company

3120 Buffalo Speedway, Houston, TX 77098

Abstract: The problem of the propagation of a hydraulically driven fracture in a fully saturated, permeable, and porous medium is investigated. Fluid driven fracture propagation in porous media is a coupled problem with four unknown fields: the flow of the fracturing fluid within the fracture, the flow of the pore fluid within the pores, the porous medium deformation, and the fracture configuration. The corresponding governing equations are the mass balance of the fracturing fluid, mass balance of the pore fluid, equilibrium of the porous medium, and fracture initiation and propagation criteria. In this work, the recently co-developed Abaqus fully-coupled hydraulic fracturing modeling capabilities are evaluated by assessing their consistency, convergence, and accuracy qualities. The Abaqus “coupled pressure/deformation cohesive elements” and “coupled pressure/deformation extended finite elements (XFEM)” are used to model the propagation of the fracture and the flow of the fracturing fluid, while the porous medium deformation and pore-fluid flow are modeled with coupled “pore-pressure/deformation” continuum finite elements. The propagation of a vertical planar fluid-driven fracture with constant height and vertically uniform width within a prismatic-shaped reservoir (KGD model), and the propagation of a horizontal, circle-shaped, planar, fluid-driven fracture within a cylindrical reservoir (“Penny-Shaped” model) are simulated in both two and three dimensions. The Abaqus numerical solution obtained with each modeling technique (cohesive and XFEM) is compared with asymptotic analytical solutions for both the KGD and Penny-shaped models in the toughness/storage dominated and viscosity/storage dominated propagation regimes. Both methods are found to accurately reproduce the analytical solutions, and converge monotonically as the mesh is refined. This validation of the newly developed hydraulic fracturing capabilities within Abaqus provides confidence in its ability and readiness to simulate fluid driven fracturing applications for the oil and gas industry including injection, stimulation, and drilling operations.

Keywords: geomechanics, soil mechanics, fracture mechanics, hydraulic fracturing, fluid-driven fracturing, geostatic, soils, pore pressure, cohesive elements, extended finite elements, XFEM, reservoir, drilling, injection.

1. Introduction

Hydraulic fracturing is a fundamental problem in Petroleum Engineering and plays a critical role in many applications within the oil and natural gas industry. The process can be generally defined as the intentional (or unintentional) initiation and propagation of a fracture due to the

2014 SIMULIA Community Conference 1 https://www.wendangku.net/doc/fb7726981.html,/simulia

pressurization of fluid that flows within the fracture. Examples of applications include (a) the stimulation of rock formations with poor or damaged permeability to increase conductivity between the reservoir and the producing wells, (b) improvement of produced water re-injection (PWRI) where water is injected to replace produced fluids and maintain reservoir pressure or provide enhanced oil recovery, (c) cuttings reinjection (CRI) where a slurry of drill cuttings is injected into a formation to mitigate the cost and risk of surface disposal, (d) in-situ stress measurement by balancing the fracturing fluid pressure in a hydraulically opened fracture with the geostatic stresses, and (e) wellbore integrity analysis of drilling operations to avoid propagating near-wellbore fractures that could result in drilling fluid losses to the formation and an inability to effectively clean the wellbore.

Knowledge of the fracture dimensions (length/width/height), fracture geometry, and wellbore pressure is crucial for both the design and integrity of hydraulic fracturing field operations. For stimulation, PWRI and CRI, one of the fundamental questions is whether or not fracture containment is achieved. This means that the injection fluid and fracture are confined to a target interval or “pay” zone for PWRI and stimulation, or a dedicated disposal domain for CRI. Other important considerations include predictions of the injection rate, pressure, or injected volume required to initiate fractures, inject under matrix conditions, or minimize the potential for inducing fractures while drilling.

Currently, there are no reliable techniques to measure fracture geometries during or after the hydraulic fracturing process. Furthermore, direct solutions of the underlying differential equations representing the different physical processes occurring during fracturing are difficult to construct, even in their most simplified forms. Therefore, the development of a numerical simulator with accurate predictive capability is of paramount importance.

The computational modelling of hydraulic fracturing of porous media is a challenging endeavor. The difficulty originates primarily from the strong non-linear coupling between the governing equations, as the process involves at least the interaction between four different phenomena: (i) the flow of the fracturing fluid within the fracture, (ii) the flow of the pore fluid and seepage of fracturing fluid within the pores, (iii) the deformation of a porous medium induced by both the hydraulic pressurization of the fracture and the compression/expansion and transport of pore fluid within the pores, and (iv) the fracture propagation which is an inherently an irreversible and singular process. Additionally, fracture propagation typically occurs in heterogeneous formations consisting of multiple layers of different rock types, subjected to in-situ confining stresses with non-uniform magnitudes and orientations. Furthermore, fracturing fluids typically exhibit nonlinear rheologies and the leakoff of these fluids from the fracture into the surrounding rock is often history dependent.

There are a number of commercial hydraulic fracture simulators used in the oil and natural gas industry for rapid design, analysis and prediction of fracture size, treating pressures, and flows (Clearly 1980, Meyer 1989, Warpinski 1994). These simulators rely in strong simplifying assumptions to render the problem solvable in realistic times:

?Fractures are assumed to be planar and symmetric with respect to the wellbore

?Fracture geometries are represented with few geometric parameters

2 2014 SIMULIA Community Conference

https://www.wendangku.net/doc/fb7726981.html,/simulia

?The formation is assumed to be unbounded and modeled using linear elasticity theory resulting in an integral equation relating fracture opening and fluid pressure ?The fracture propagation is modeled within the framework of linear elastic fracture mechanics without any consideration of pore fluid pressure effects

?Leakage of fracturing fluid from the fracture into the rock is modeled as one dimensional and decoupled from the porous medium deformation.

Although these simulators are useful in predicting broad trends and upper/lower bounds in operational parameters, their reliability and accuracy are restricted to unrealistic scenarios intrinsic simplistic assumptions apply, i.e., situations where some of the coupling between the many different processes involved can be neglected, and with strong symmetry in confinement stresses and geology.

The accurate modelling of the hydraulic fracturing process under realistic geologies, wellbore configurations, confining stress states, and operational conditions calls for a more advanced, multi-physics numerical simulator that incorporates the complex coupling between the injected fluid, the pore fluid, the rock deformation, and the fracture configuration, thus overcoming the limitations of currently available commercial simulation tools.

To this end, fully-coupled hydraulic fracturing simulation capabilities that leverage (i) the existing Abaqus non-linear soil consolidation analysis solver, (ii) Abaqus cohesive elements for modelling interface decohesion, and (iii) Abaqus extended finite element method (XFEM) for modelling propagating discontinuities, are being co-developed between ExxonMobil Upstream Research Company and Dassault Systemes Simulia Corporation.

Specifically, two new element classes have been integrated into the existing Abaqus/Standard coupled pore fluid diffusion and solid stress porous media analysis solver:

i. A coupled pressure/deformation cohesive element that models the progressive damage of

normal mechanical strength and normal hydraulic conductivity as well as the flow of

fracturing fluid within the opening fracture.

ii.An enriched version of the continuum coupled pore fluid diffusion/stress elements capable of activating arbitrarily oriented discontinuities in both displacements and pore

pressures while simultaneously modelling the fracturing fluid flow along the fracture. This work describes and validates these two new formulations for hydraulic fracturing modeling by assessing consistency, accuracy and convergence qualities. The propagation of a fluid-driven vertical planar fracture of uniform width and constant height within a prismatic-shaped rock formation (Khristianovich-Geertsma-de Klerk, or KGD model) and the propagation of a horizontal, circle-shaped, planar, fluid-driven fracture within a cylindrical reservoir (radial or “Penny-Shaped” model) are simulated for both two and three dimensions (Clearly 1980, Geertsma 1969, Yew 1997). The numerical solution obtained with each new modeling technique (cohesive and XFEM) are then compared with available asymptotic analytical solutions for both the KGD and Penny-shaped models in the toughness/storage dominated and viscosity/storage dominated propagation regimes. Finally, the consistency, accuracy and convergence attributes are assessed for both methods.

2014 SIMULIA Community Conference 3 https://www.wendangku.net/doc/fb7726981.html,/simulia

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

Section 2 describes the governing equations for each of the coupled processes as well as the

constitutive and kinetic relations assumed for the porous medium, pore fluid and fracturing fluid, including:

Equilibrium equation for the porous medium ii.

Constitutive equation for the porous medium (Biot’s theory of poroelasticity) iii.

Continuity equation for the pore fluid iv.

Continuity equation for the fracturing fluid v.

Momentum equation for the pore fluid (Darcy’s Law) vi. Momentum equation for the fracturing (Lubrication Equation)

Section 3 details the procedures employed by both formulations (cohesive and extended finite element methods) and the fracture initiation and propagation criteria. Section 4 defines the test models (KGD plane-strain and the Penny-Shaped models) and the model set-up and assumptions used within Abaqus, while Section 5 presents numerical results and assesses accuracy and convergence by comparing the main solution variables obtained with meshes of different

resolutions with available asymptotic analytical solutions. Finally, some concluding remarks are summarized in Section 6.

2. Governing Equations

As stated in the introduction, hydraulic fracturing involves the interaction between four different phenomena:

Porous medium deformation ii.

Pore fluid flow iii.

Fracturing fluid flow iv. Fracture propagation

The equations and constitutive relation governing these coupled processes, i.e., Biot’s theory of poroelasticity for porous media, Darcy’s Law for pore fluid flow, Reynold’s lubrication theory for fracturing fluid flow and the cohesive zone model to characterize fracturing (Abaqus 2013, Charlez 1997)) are summarized in what follows.

2.1 Porous Media Deformation

Porous media can be modelled as an isotropic, poroelastic material undergoing quasistatic

deformation. The equilibrium equation enforced by Abaqus, when body forces are neglected is,

, 0

(1) while the poroelastic constitutive relation, assuming small strains, is given by, 2 (2) 2 1

in which is Biot’s coefficient, and are the dry elastic shear and bulk moduli, is the dry Young’s modulus, and is the dry Poisson’s ratio. Abaqus is formulated in terms of Terzaghi effective stresses ′, defined for fully saturated media as (Abaqus 2013, Charlez 1997)

In terms of the latter, the constitutive relation takes the form

2 1

Defining effective strains as

the constitutive relation simplifies to

This identity is identical to the constitutive relation for linear elastic materials, but expressed in terms of Terzaghi effective stresses ′ and effective strains ′. Abaqus internally translates total stresses and strains into Terzaghi effective stresses and strains to leverage this equivalence (Abaqus 2013).

2.2 Pore Fluid Flow

The continuity equation for the pore fluid is, assuming small volumetric strains, given by

, 0

where is the pore fluid seepage velocity, and and are Biot’s modulus and Biot’s coefficient, respectively. These two poroelastic constants are defined by the identities

where is the pore fluid bulk modulus, is the porous medium solid grain bulk modulus, and is the initial porosity. In Abaqus, the two compressibilities and are specified using the *POROUS BULK MODULI keyword. Pore fluid is assumed to flow through an interconnected pore network according to Darcy’s law

in which is the permeability, is the pore fluid viscosity, is the hydraulic conductivity and is the pore fluid specific weight. Combining with the continuity equation, the pore fluid diffusion equation is obtained

2014 SIMULIA Community Conference 5 https://www.wendangku.net/doc/fb7726981.html,/simulia

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

1 , (3) Within Abaqus, the hydraulic conductivity and specific weight are specified through the *PERMEABILITY keyword.

2.3 Fracturing Fluid Flow

Longitudinal fluid flow within the fracture is governed by Reynold’s lubrication theory defined by the continuity equation

0 and the momentum equation for incompressible flow and Newtonian fluids through narrow parallel plates (i.e., Poiseuille flow)

where is the fracture gap (Figure 1), ? is the fracturing fluid flow (per unit width) across the fracture, and are the normal flow velocities of fracturing fluid leaking through the top and bottom faces of the fracture into the porous medium, is the fracturing fluid viscosity, and is the fracturing fluid pressure along the fracture surface parameterized with the curvilinear coordinate, .

Figure 1: Fracture aperture, width and fracturing fluid flow

Abaqus computes the normal fracturing fluid velocities as

where and are the pore fluid pressures on the top and bottom surface of the fracture and and are the so-called “leakoff coefficients”. This simple leakoff model simulates a layer of filtrate that might accumulate and reduce the effective normal permeability of the fracture surfaces.

Inserting the Poiseuille flow equation and the simplified leakoff model into the continuity equation for the fracturing fluid yields the final form: ?

2014 SIMULIA Community Conference 7 https://www.wendangku.net/doc/fb7726981.html,/simulia

(5) Abaqus specifies the fracturing flow viscosity and leakoff coefficients with the *GAPFLOW and the *LEAKOFF keywords, respectively.

2.4 Fracture Initiation and Propagation

Fracturing can be conceptualized as the transition between two limiting states: the undamaged state with continuous displacements and non-zero tractions in all directions and the fully damaged state characterized by the presence of a displacement discontinuity along a material interface with zero tractions in the direction normal to the interface. In Abaqus, this transition process is modeled as a progressive degradation of cohesive strength along a zero-thickness interface whose

orientation and extent is either predefined (cohesive element method) or calculated during the simulation (extended finite element method). The gradual loss of strength in the interface with increasing separation is defined with an interface traction/interface separation relation or cohesive law (Abaqus 2013, Ortiz 1999).

For the purpose of this study, a traction-separation cohesive law with linear softening (Figure 2) is assumed, defined by the cohesive energy (area under the softening part of the traction separation curve) and the cohesive strength . For the cohesive element procedure, it is also required to define the traction-separation behavior prior to damage initiation, which is assumed to be linear with initial stiffness . The cohesive traction of the interface thus evolves from a maximum tensile strength at damage initiation, down to zero when the interface is fully damaged and free to open beyond the total separation . If the interface is unloaded prior to complete damage, the traction is assumed to ramp down linearly with a damaged stiffness . The interface effective tractions are therefore given by

Figure 2: Cohesive law for Cohesive and XFEM procedures

Upon damage initiation, the fracture is pressurized by instantaneously applying the fracturing fluid pressure, , calculated from the fracturing fluid equations in Equation 5. The total tractions resisted by (and acting on) the interface elements are therefore given by T r a c t i o n

Separation

COHESIVE

XFEM

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

As stated in Section 2.1, the Abaqus porous media analysis solver is formulated in terms of Terzaghi effective stresses. Therefore, the cohesive strength defining the onset of interface decohesion must be understood in terms of an effective strength (and not total strength).

3. Cohesive Element and Extended Finite Element Methods

The evolution of a fracture is modelled in Abaqus through zero-thickness interface elements with separation resisted by gradually decreasing tensile tractions. For the cohesive element procedure, these interface elements are defined a priori and placed between continuum element faces, whereas in the extended finite element method (XFEM), they are inserted and oriented

automatically during the course of the simulation within existing continuum elements.

3.1 Cohesive Element Method

The coupled pressure/deformation cohesive elements implemented in Abaqus (COHPE4P , COHAX4P , COH3D8P ) are standard linear isoparametric elements with displacement and pore pressure degrees of freedom associated with their corner nodes, as depicted in Figure 3 (nodes 1,2,3,4). These

elements must be inserted a priori between the faces of adjacent pressure diffusion/stress elements (CPE4P , CAX4P , C3D8P ) in order to model the yet to open fracture. To accommodate the coupling of the fracturing fluid flow equations, the elements are equipped with additional pressure degrees of freedom (attached to the center of the element edges perpendicular to the fracture) to interpolate the fracturing fluid pressure after damage initiation (nodes 5 and 6, Figure 3).

Figure 3: Coupled pressure/deformation cohesive elements for hydraulic fracturing The cohesive elements can have arbitrary undeformed geometric thickness as the instantaneous gap coupled in the fracturing fluid flow equation (Equation 5) is defined in Abaqus as the

difference between the deformed and underfed thickness, i.e., . Prior to damage, the top and bottom faces of the unopened fracture are subjected to the pore fluid pressure acting towards increasing separation and the cohesive effective tractions resisting separation,

where is the stiffness of the cohesive element prior to failure (Figure 2). After damage

initiation, the pore fluid is displaced by the fracturing fluid pressurizing the interface. The total tractions acting on the top and bottom faces of the opening fracture are then substituted by

5 64Undeformed

Configuration Deformed Configuration 243 1

2 563

2014 SIMULIA Community Conference 9 https://www.wendangku.net/doc/fb7726981.html,/simulia

where is the damaging stiffness (Figure 2). A coupled cohesive element method for hydraulic fracturing similar to the formulation just outlined is described by Boone, 1990 and Carrier, 2012.

3.2 Extended Finite element method (XFEM)

The Extended Finite Element Method (XFEM) is implemented within Abaqus using the so called “phantom node” approach (Abaqus 2013, Remmers 2008, Song 2006, Sukumar 2003). In this implementation, each enriched pressure diffusion/stress element (CPE4P , CAX4P , C3D8P ) is

internally duplicated with the addition of corner phantom nodes, as depicted in Figure 4, in which original nodes are represented with full circles and corner phantom nodes with hollow circles. Prior to damage initiation only one copy of the element is active. Upon damage initiation the displacement and pore pressure degrees of freedom associated with the corner phantom nodes are activated and both copies of the element are allowed to deform independently, pore pressures are allowed to diffuse independently, and the created interface behavior is enforced with a traction-separation cohesive law.

Figure 4: Implementation of the XFEM with “corner” and “edge” phantom nodes In order to enable the solution of the fracturing fluid flow equations, the enriched elements also incorporate new “edge-phantom nodes” (depicted as red triangles in Figure 4) that interpolate the fracturing fluid pressure within the fracture. The pore fluid pressure and at the top and bottom faces of the fracture are interpolated from the pore pressure degrees of freedom at the corner real nodes and phantom nodes. The difference with the fracturing fluid pressure

(interpolated at the edge-phantom node) is the driving force that controls the leakage of fracturing fluid into the porous medium (Equation 4).

The fracture is extended to a new element ahead of the fracture tip when the maximum effective principal stress at this element (interpolated to the tip) in a given iteration is equal to the cohesive strength . The orientation of the fracture segment to be extended into the tip element is set to the Undeformed

configuration

after damage Deformed configuration after damage 2,2 4,4 3,3 1,1

2 1

3 4 1 2

2 4

3 1 Undeformed configuration before damage

direction perpendicular to the minimum principal stress of the current iteration. This fracture initiation/orientation criterion is defined in Abaqus through the keyword

*DAGAMAGE INITATION, CRITERION=MAXPS, POSITION=CRACKTIP

As in the cohesive element formulation, the fracturing fluid pressure is applied to the top and bottom faces of the fracture and superposed to the cohesive tractions.

4. Benchmark Models

In this section, the two formulations previously outlined (coupled pressure/displacement cohesive and extended finite elements) are applied to model the propagation of a hydraulically driven fracture in two different configurations:

i.Horizontal, circle-shaped, planar, fracture within a cylindrical domain, (radial or “Penny-

Shaped” model (Clearly 1980, Charlez 1997, Yew 1997))

ii.Vertical, rectangle-shaped, planar fracture within a prismatic-shaped domain (Khristianovich-Geertsma-de Klerk, or KGD model (Charlez 1997, Geertsma 1969, Yew 1997)).

These models serve as benchmark examples to assess the consistency, convergence and accuracy of the numerical solution obtained with Abaqus.

4.1 Fracture Propagation Regimes

Despite the simplicity of the fracture geometry and strong symmetry in the chosen benchmark problems, there are no available closed-form analytical solutions for these problems when all coupled processes are considered in the analysis, i.e., when the formation is assumed to be porous and permeable with pore fluid flow and fracturing fluid is leaking into the pore space displacing the pore fluid. However, using the more restrictive theoretical framework resulting from assuming (i) an infinite domain, (ii) material fully impermeable, (iii), material linear elastic, (iv) linear-elastic fracture mechanics, and (v) Carter’s leakoff model (Howard 1957, Charlez 1997), approximate analytical solutions exist in the form of regular asymptotic expansions (Bunger 2005, Detournay 2006, Garagash 2006, Hu 2010, Garagash 2011, Peirce 2008, Savitski 2002). The governing equations then simplify to (i) the equilibrium equation for the linear elastic material, that for an infinite domain can be represented as an singular integral equation relating fracture opening and fluid pressure, (ii) the local and global mass balance equations for the fracturing fluid, and (iii) the fracture propagation criterion, also expressible as a singular integral equation relating fracturing pressure and fracture toughness. A non-dimensional analysis of this reduced system of equations uncovers the presence of two pairs of competing physical processes. The first pair consists of competing dissipative mechanisms: (a) energy dissipated due to fluid viscosity and (b) energy dissipated due to fracture propagation; the second pair consists of competing components of fluid balance: (a) fluid storage within the fracture and (b) fluid leakage from the fracture into the surrounding material. Depending on which of the two dissipative mechanisms and which of the two storage mechanisms dominate, four primary limiting regimes of propagation emerge: ?Viscosity dominated and storage dominated propagation regime ( ).

?Toughness dominated and storage dominated propagation regime ( ).

10 2014 SIMULIA Community Conference

https://www.wendangku.net/doc/fb7726981.html,/simulia

2014 SIMULIA Community Conference 11 https://www.wendangku.net/doc/fb7726981.html,/simulia

Viscosity dominated and leak-off dominated ( ). ? Toughness dominated and leak-off dominated regime (

). These four fracture propagation regimes can be conceptualized in a rectangular parametric space where each limiting regime corresponds to each of the vertices of the rectangle with one dissipation mechanism dominating and the other being neglected, and one component of fluid global balance dominating with the other also neglected (Figure 5).

Figure 5: Parametric diagram representing the four limiting propagation regimes of

hydraulically driven fractures

This work analyzes each benchmark problem in both the toughness/storage-dominated (near- ) and the viscosity/storage-dominated (near- ) propagation regimes. The near- and the near- asymptotic solutions (small time solutions in the toughness and viscosity regimes) are used to compare to Abaqus numerical solutions for each formulation (cohesive element method and XFEM) with material parameters, loads, and boundary conditions that reproduce each of these propagation regimes.

In order to render the Abaqus solution comparable with the asymptotic solutions, the model

dimensions and material properties are selected such that the more restrictive conditions for which these solutions apply are adequately recreated. Specifically, the dimensions of the domain of analysis are much bigger than the fracture aperture and length, the permeability is defined to

minimize the influence of poroelastic effects ahead of the fracture tip, and cohesive properties are selected to ensure a small cohesive zone relative to the size of the fracture.

4.2 Radial (Penny-Shaped) Model

The first benchmark problem consists of an axisymmetric, penny-shaped, hydraulically-driven fracture propagating in a cylindrically shaped poroelastic formation as illustrated in Figure 6.

∞

0 ∞

∞

12 2014 SIMULIA Community Conference

https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 6: Cylindrical domain with a horizontal, circular-shaped, hydraulically driven

fracture

The domain of analysis is characterized by the inner radius , outer radius , and height . The porous medium is characterized by Young’s modulus , Poisson ratio , fracture toughness , porosity , Biot’s coefficient , Biot’s modulus , and hydraulic conductivity

. An incompressible Newtonian fluid with viscosity is injected at a constant rate at the center of the fracture from a vertical wellbore. The fracture aperture , , the net pressure , (defined as the difference between the fracturing fluid pressure , and the confining stress ), and the fracture radius are the sought quantities.

4.3 Plane Strain (KGD) Model

The second benchmark problem considers a hydraulically-driven vertical fracture propagating in a poroelastic prismatic-shaped formation of length L, depth R and height H as illustrated in Figure 7. 0

2014 SIMULIA Community Conference 13 https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 7: Prismatic domain with a planar vertical, hydraulically driven fracture Again, the porous medium is characterized by Young’s modulus , Poisson ratio and fracture toughness , porosity , Biot’s coefficient , Biot’s modulus , and hydraulic conductivity with an incompressible Newtonian fluid with viscosity injected along a vertical wellbore at a constant rate per unit of vertical length . The unknowns of the problem are the fracture aperture , , the net pressure , , , and the fracture radius .

4.4 Model Setup

4.4.1 Geometry and Mesh for the Radial (Penny-Shaped) Fracture Model

For each propagation regime (toughness/storage and viscosity/storage dominated) and modeling procedure (cohesive element method and XFEM) two models are constructed: (a) 2D

axisymmetric model (Figure 8) and (b) 3D model (Figure 9). Both models are expected to render similar solutions. Due to the axial symmetry of the problem only one quarter of the 3D cylindrical domain is considered.

, ,

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 8: Axisymmetric and plane strain model meshes, cohesive element and

XFEM procedures

Figure 9: 3D model meshes, cohesive element and XFEM procedures

The formation is discretized with linear coupled pore fluid diffusion/stress elements (CAX4P for the axisymmetric model, and C3D8P for the 3D model). For the cohesive element method, coupled pressure/deformation cohesive elements (COHAX4P for the axisymmetric model, and COH3D8P for the 3D model) are inserted along the horizontal mid-plane of the domain in which the fracture is expected to propagate. Continuum elements and cohesive elements are connected with shared nodes.

Coupled pressure/deformation

cohesive elements (COHAX4P ) Coupled pressure/deformation continuum elements (CAX4P ) enriched

Coupled pressure/deformation

continuum elements (C3D8P ) enriched

Coupled pressure/deformation

cohesive elements (COH3D8P )

For the XFEM model, an equivalent mesh is used where the layer of cohesive elements is replaced by a set of enriched coupled pore fluid diffusion/stress elements (CAX4P and C3D8P), as depicted in Figures 8 and 9. To declare a set of continuum elements as enriched in which a fluid-driven fracture may propagate the following keyword is required:

*ENRICHMENT NAME=enrichmentname, TYPE=PROPAGATION CRACK, ELSET=elsetname

The dimensions of the analysis domain are chosen large enough to minimize boundary effects (Table 1).

Table 1: Penny-Shaped fracture model dimensions

Dimension Value

Inner radius 0.01

Outer radius 45

Height 30

Angle 90°

In order to analyze the convergence of the two modeling procedures as the mesh is refined, three meshes of decreasing element size are recreated for the 2D axisymmetric and 3D continuum cases (Figures 10 and 11). To minimize model size, only the elements in the vicinity of the fracture plane are refined. These elements are generated by subdividing elements in the parent (coarse) mesh into four (axisymmetric model) or eight (3D model) child elements. The compatibility of displacements and pore fluid pressures on the nodes lying at the intersection between coarse and fine regions of the resultant meshes are enforced with multi point constraints (*MPC, LINEAR, and *MPC, BILINEAR). This refinement technique is chosen to minimize element distortion. Table 3 lists the number of elements and nodes for the respective meshes.

2014 SIMULIA Community Conference 15 https://www.wendangku.net/doc/fb7726981.html,/simulia

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 10: Meshes of different resolution used for the 2D penny-shaped and planar

fracture models (axisymmetric and plane strain)

Coarse

Intermediate

Fine COH E XFEM

2014 SIMULIA Community Conference 17 https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 11: Meshes of different resolution used for the Penny-shaped fracture

model

4.4.2 Geometry and Mesh for the Planar (KGK) Fracture Model

For the second benchmark problem (planar fracture geometry) two pairs of models are also

constructed: (a) 2D plane-strain model with nodal coordinates and element connectivity identical to the axisymmetric model (Figures 8 and 10) and (b) 3D model as depicted in Figure 12. As in the case of the penny-shaped mesh, the formation is discretized with linear coupled pore fluid diffusion/stress elements (CPE4P for the plane strain model and C3D8P for the 3D model), while the fracture plane is modeled with either plane strain (COH2D4P ) or 3D (COH3D8P ) coupled

pressure/deformation cohesive elements (i.e. cohesive element method), or the enriched version of the coupled pore fluid diffusion/stress continuum elements (i.e. XFEM).

Table 2 summarizes the planar fracture model dimensions chosen to minimize boundary effects.

Table 2: Planar fracture model dimensions Dimension Value

Width 45

Height 45

Length 60

Coarse

Intermediate

Fine

COH

E XFEM

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

Figure 13 details the sequence of meshes of different resolution constructed to analyze

convergence for this second benchmark model.

Figure 12: 3D mesh for the KGD model, cohesive procedures

Figure 13: Meshes of different resolution used for the KGD fracture model Coupled pressure/deformation

cohesive elements (COH3D8P ) Coupled pressure/deformation

continuum elements (C3D8P ) enriched

Table 3: Number of elements and nodes for each mesh Model Mesh Elements Nodes

2D Coarse 2745 2898 Intermediate 3330 3577

Fine 4500 4931

3D Coarse 54900 60858 Intermediate 82800 95137

FIne 194400

230291

4.4.3 Material Properties

The poroelastic material parameters used for all simulations are summarized in Table 4.

Table 4: Poroelastic material parameters

Parameter Value

Young modulus 17

Poisson’s ratio 0.2

Biot coefficient 0.75

Porosity 0.2

Pore fluid specific weight 0.0098

Hydraulic conductivity 9.8 10

Biot modulus 68.7

To define the cohesive behavior of the material, a traction-separation law with linear softening is used (*DAMAGE EVOLUTION, TYPE=ENERY, SOFTENING=LINEAR) defined by a cohesive energy

120 ? . Fracture toughness can be estimated as

1.46 √ . The quadratic nominal stress fracture initiation criterion (*DAMAGE INITIATION, CRITERION=QUADS) with a cohesive strength of 1.25 is used in all simulations. For the cohesive element method, the elastic stiffness prior to failure (*ELASTIC, TYPE=TRACTION) is defined as 100 times the Young’s modulus of the material.

The fracturing fluid is modeled as a Newtonian fluid (*GAPFLOW, TYPE=NEWTONIAN) with viscosity equal to 100 for the viscosity/storage dominated propagation regime and

0.1 for the toughness/storage dominated regime.

4.4.4 Initial Conditions, Boundary Conditions and Loads

The following initial conditions are specified for all solid sections with poroelastic constitutive behavior (Table 5):

Table 5: Initial conditions

Pore fluid Initial saturation 1.0

Pore fluid initial pressure 0.0

Rock initial void ratio 0.25

Initial principal in-situ stresses , , 0.0,0.0,0.0

The media is assumed to be initially fully saturated and with no confinement stresses. Abaqus requires an explicit specification of the set of cohesive elements or enriched elements where the fracturing fluid flow equations will be initially solved. These initial conditions are specified

2014 SIMULIA Community Conference 19 https://www.wendangku.net/doc/fb7726981.html,/simulia

2014 SIMULIA Community Conference https://www.wendangku.net/doc/fb7726981.html,/simulia

through the keyword *INITIAL CONDITIONS, TYPE=INITIALGAP , and *INITIAL CONDITIONS, TYPE=ENRICHMENT , respectively. Only the set of elements containing the nodes where the fluid is injected are needed to meet this requirement (Figure 14).

Figure 14: Initial set of cohesive elements and initial enrichment required to define

the initial solution domain for the fracturing fluid flow equations

The displacements in the normal direction to all boundary surfaces and symmetry planes are constrained. Additionally, pore fluid pressures are fixed to the uniform value 0 on all model boundaries. For the XFEM model, the corner phantom nodes on the symmetry surfaces and

boundary surfaces are constrained to move within these surfaces (Figure 15). This special type of boundary condition is enforced with the keyboard *BOUNDARY, PHANTOM=NODE .

Figure 15: Displacement boundary conditions on corner phantom nodes

In all cases, fracturing fluid is injected at a constant rate 0.001 and injection is simulated for 40 . As explained in the previous section, the fracturing fluid pressure degrees of freedom are associated with the mid-edge nodes of the cohesive elements and the edge phantom nodes of the enriched elements. Therefore, concentrated fracturing fluid flow must be applied directly to these mid-edge (*CFLOW ) and phantom edge nodes (*CFLOW, PHANTOM=EDGE ), as depicted in Figure 16. Cohesive Element Model XFEM model

Corner phantom nodes constrained Corner phantom nodes unconstrained