Viscosity-Storage-Dominated PKN Hydraulic Fracture

Context

In this example, we simulate the propagation of a Perkins–Kern–Nordgren (PKN) fracture in a viscosity-storage-dominated regime, a classic benchmark in hydraulic fracturing. The developed planar fracture displays an elliptical vertical cross-section. Unlike KGD and penny-shaped fractures, the growth height of a PKN fracture is constrained by mechanical barriers (such as bedding layers, sedimentary laminations, or weak interfaces), thus promoting lateral propagation. This problem is solved using the hydrofracture solver in GEOS to obtain the temporal evolutions of the fracture characteristics (length, aperture, and pressure). We validate these simulated values against existing analytical solutions (Kovalyshen and Detournay, 2010; Economides and Nolte, 2000).

Input file

This example uses no external input files. Everything we need is contained within two GEOS input files:

inputFiles/hydraulicFracturing/pknViscosityDominated_base.xml

inputFiles/hydraulicFracturing/pknViscosityDominated_benchmark.xml

Python scripts for post-processing and visualizing the simulation results are also prepared:

inputFiles/hydraulicFracturing/scripts/hydrofractureQueries.py

inputFiles/hydraulicFracturing/scripts/hydrofractureFigure.py

Description of the case

In this example, a hydraulic fracture initiates and propagates from the center of a 20m-thick layer. This layer is homogeneous and bounded by neighboring upper and lower layers. For viscosity-dominated fractures, more energy is necessary to move the fracturing fluid than to split the intact rock. If fluid leak-off is neglected, storage-dominated propagation occurs with most of the injected fluid confined within the open surfaces. To meet the requirements of the viscosity-storage-dominated assumptions, impermeable domain (no fluid leak-off), incompressible fluid with constant viscosity ( $1.0 cp$ ) and ultra-low rock toughness ( $0.1 MPa { \sqrt{m} }$ ) are chosen in the GEOS simulation. With these parameters, the fracture stays within the target layer; it extends horizontally and meets the conditions of the PKN fracture in a viscosity-storage-dominated regime.

We assume that the fluid injected in the fracture follows the lubrication equation resulting from mass conservation and Poiseuille’s law. The fracture propagates by creating new surfaces if the stress intensity factor exceeds the local rock toughness $K_{IC}$ . As the geometry of the PKN fracture exhibits symmetry, the simulation is reduced to a quarter-scale. For verification purposes, a plane strain deformation is considered in the numerical model.

We set up and solve a hydraulic fracture model to obtain the temporal solutions of the fracture half length $l$ , the net pressure $p_0$ and the fracture aperture $w_0$ at the fluid inlet for the PKN fracture propagating in viscosity-storage-dominated regime. Kovalyshen and Detournay (2010) and Economides and Nolte (2000) derived the analytical solutions for this classic hydraulic fracture problem, used here to verify the results of the GEOS simulations:

$l(t) = 0.3817 (\frac{ E_p Q_0^3 t^4 }{ \mu h^4 })^{ 1/5 }$

$w_0(t) = 3 (\frac{ \mu Q_0 l }{ E_p })^{ 1/4 }$

$p_0(t) = (\frac{ 16 \mu Q_0 E_p^3 l }{ \pi h^4 })^{ 1/4 }$

where the plane modulus $E_p$ is related to Young’s modulus $E$ and Poisson’s ratio $\nu$ :

$E_p = \frac{ E }{ 1-\nu^2 }$

For this example, we focus on the Mesh, the Constitutive, and the FieldSpecifications tags.

Mesh

The following figure shows the mesh used in this problem.

We use the internal mesh generator to create a computational domain ( $400\, m \, \times 400 \, m \, \times 800 \, m$ ), as parametrized in the InternalMesh XML tag. The structured mesh contains 105 x 105 x 60 eight-node brick elements in the x, y, and z directions respectively. Such eight-node hexahedral elements are defined as C3D8 elementTypes, and their collection forms a mesh with one group of cell blocks named here cb1. Local refinement is performed for the elements in the vicinity of the fracture plane.

  <Mesh>
    <InternalMesh
      name="mesh1"
      elementTypes="{ C3D8 }"
      xCoords="{ 0, 150, 200, 400 }"
      yCoords="{ 0, 150, 200, 400 }"
      zCoords="{ -400, -100, -20, 20, 100, 400 }"
      nx="{ 75, 10, 20 }"
      ny="{ 75, 10, 20 }"
      nz="{ 10, 10, 20, 10, 10 }"
      cellBlockNames="{ cb1 }"/> 
  </Mesh>

The fracture plane is defined by a nodeset occupying a small region within the computational domain, where the fracture tends to open and propagate upon fluid injection:

    <Box
      name="core"
      xMin="{ -500.1, -500.1, -0.1 }"
      xMax="{ 500.1, 10.1, 0.1 }"/>

Solid mechanics solver

GEOS is a multi-physics platform. Different combinations of physics solvers available in the code can be applied in different regions of the domain and be functional at different stages of the simulation. The Solvers tag in the XML file is used to list and parameterize these solvers.

Three elementary solvers are combined in the solver hydrofracture to model the coupling between fluid flow within the fracture, rock deformation, fracture deformation and propagation:

    <Hydrofracture
      name="hydrofracture"
      solidSolverName="lagsolve"
      flowSolverName="SinglePhaseFlow"
      surfaceGeneratorName="SurfaceGen"
      logLevel="1"
      targetRegions="{ Fracture }"
      maxNumResolves="5"
      initialDt="0.1">
      <NonlinearSolverParameters
        newtonTol="1.0e-4"
        newtonMaxIter="10"
        maxTimeStepCuts="5"
        maxAllowedResidualNorm="1e+15"/>
      <LinearSolverParameters
        solverType="gmres"
        preconditionerType="mgr"
        logLevel="1"
        krylovAdaptiveTol="1"/>
    </Hydrofracture>

Rock and fracture deformations are modeled by the solid mechanics solver SolidMechanicsLagrangianSSLE. In this solver, we define targetRegions that includes both the continuum region and the fracture region. The name of the contact constitutive behavior is specified in this solver by the contactRelationName.

    <SolidMechanicsLagrangianSSLE
      name="lagsolve"
      timeIntegrationOption="QuasiStatic"
      discretization="FE1"
      targetRegions="{ Domain, Fracture }"
      contactRelationName="fractureContact"
      contactPenaltyStiffness="1.0e0">
      <NonlinearSolverParameters
        newtonTol="1.0e-6"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-10"/>    
    </SolidMechanicsLagrangianSSLE>

The single-phase fluid flow inside the fracture is solved by the finite volume method in the solver SinglePhaseFVM.

    <SinglePhaseFVM
      name="SinglePhaseFlow"
      discretization="singlePhaseTPFA"
      targetRegions="{ Fracture }">
      <NonlinearSolverParameters
        newtonTol="1.0e-5"
        newtonMaxIter="10"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-12"/>     
    </SinglePhaseFVM>

The solver SurfaceGenerator defines the fracture region and rock toughness rockToughness="0.1e6". With nodeBasedSIF="1", a node-based Stress Intensity Factor (SIF) calculation is chosen for the fracture propagation criterion.

    <SurfaceGenerator
      name="SurfaceGen"
      targetRegions="{ Domain }"
      nodeBasedSIF="1"      
      rockToughness="0.1e6"
      mpiCommOrder="1"/>

Constitutive laws

For this problem, a homogeneous and isotropic domain with one solid material is assumed. Its mechanical properties and associated fluid rheology are specified in the Constitutive section. ElasticIsotropic model is used to describe the mechanical behavior of rock when subjected to fluid injection. The single-phase fluid model CompressibleSinglePhaseFluid is selected to simulate the response of water upon fracture propagation.

  <Constitutive>
    <CompressibleSinglePhaseFluid
      name="water"
      defaultDensity="1000"
      defaultViscosity="0.001"
      referencePressure="0.0"
      compressibility="5e-12"
      referenceViscosity="1.0e-3"
      viscosibility="0.0"/>

    <ElasticIsotropic
      name="rock"
      defaultDensity="2700"
      defaultBulkModulus="20.0e9"
      defaultShearModulus="12.0e9"/>

    <CompressibleSolidParallelPlatesPermeability
      name="fractureFilling"
      solidModelName="nullSolid"
      porosityModelName="fracturePorosity"
      permeabilityModelName="fracturePerm"/>

    <NullModel
      name="nullSolid"/>

    <PressurePorosity
      name="fracturePorosity"
      defaultReferencePorosity="1.00"
      referencePressure="0.0"
      compressibility="0.0"/>

    <ParallelPlatesPermeability
      name="fracturePerm"/>

    <FrictionlessContact
      name="fractureContact"/>

    <HydraulicApertureTable
      name="hApertureModel"
      apertureTableName="apertureTable"/>    

  </Constitutive>

All constitutive parameters such as density, viscosity, bulk modulus, and shear modulus are specified in the International System of Units.

Time history function

In the Tasks section, PackCollection tasks are defined to collect time history information from fields. Either the entire field or specified named sets of indices in the field can be collected. In this example, pressureCollection, apertureCollection, hydraulicApertureCollection and areaCollection are specified to output the time history of fracture characterisctics (pressure, width and area). objectPath="ElementRegions/Fracture/FractureSubRegion" indicates that these PackCollection tasks are applied to the fracure element subregion.

  <Tasks>
    <PackCollection
      name="pressureCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="pressure"/>

    <PackCollection
      name="apertureCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="elementAperture"/> 

    <PackCollection
      name="hydraulicApertureCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="hydraulicAperture"/>  

    <PackCollection
      name="areaCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="elementArea"/>

    <!-- Collect aperture, pressure at the source for curve checks -->
    <PackCollection
      name="sourcePressureCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="pressure"
      setNames="{ source }"/>

    <PackCollection
      name="sourceHydraulicApertureCollection"
      objectPath="ElementRegions/Fracture/FractureSubRegion"
      fieldName="hydraulicAperture"
      setNames="{ source }"/>
  </Tasks>

These tasks are triggered using the Event manager with a PeriodicEvent defined for the recurring tasks. GEOS writes one file named after the string defined in the filename keyword and formatted as a HDF5 file (pknViscosityDominated_output.hdf5). This TimeHistory file contains the collected time history information from specified time history collector. This file includes datasets for the simulation time, fluid pressure, element aperture, hydraulic aperture and element area for the propagating hydraulic fracture. A Python script is prepared to read and query any specified subset of the time history data for verification and visualization.

Initial and boundary conditions

The next step is to specify:

The initial values: the waterDensity, separableFace and the ruptureState of the propagating fracture have to be initialized,

The boundary conditions: fluid injection rates and the constraints of the outer boundaries have to be set.

In this example, a mass injection rate SourceFlux (scale="-6.625") is applied at the surfaces of the initial fracture. Only one fourth of the total injection rate is used because only a quarter of the fracture is modeled (the problem is symmetric). The value given for scale is $Q_0 \rho_f/4$ (not $Q_0 /4$ ). All the outer boundaries are subject to roller constraints. These boundary conditions are set through the FieldSpecifications section.

  <FieldSpecifications>   
    <FieldSpecification
      name="waterDensity"
      initialCondition="1"
      setNames="{ fracture }"
      objectPath="ElementRegions"
      fieldName="water_density"
      scale="1000"/>

    <FieldSpecification
      name="separableFace"
      initialCondition="1"
      setNames="{ core }"
      objectPath="faceManager"
      fieldName="isFaceSeparable"
      scale="1"/>

    <FieldSpecification
      name="frac"
      initialCondition="1"
      setNames="{ fracture }"
      objectPath="faceManager"
      fieldName="ruptureState"
      scale="1"/>
  
    <FieldSpecification
      name="yconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="1"
      scale="0.0"
      setNames="{ yneg, ypos }"/>

    <FieldSpecification
      name="zconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="2"
      scale="0.0"
      setNames="{ zneg, zpos }"/>

    <FieldSpecification
      name="xconstraint"
      objectPath="nodeManager"
      fieldName="totalDisplacement"
      component="0"
      scale="0.0"
      setNames="{ xneg, xpos }"/>
   
    <SourceFlux
      name="sourceTerm"
      objectPath="ElementRegions/Fracture"
      scale="-6.625"
      setNames="{ source }"/>
  </FieldSpecifications>

The parameters used in the simulation are summarized in the following table.

Symbol	Parameter	Unit	Value
$K$	Bulk Modulus	[GPa]	20.0
$G$	Shear Modulus	[GPa]	12.0
$K_{IC}$	Rock Toughness	[MPa.m^1/2]	0.1
$\mu$	Fluid Viscosity	[Pa.s]	1.0x10^-3
$Q_0$	Injection Rate	[m³/s]	0.0265
$t_{inj}$	Injection Time	[s]	200
$h_f$	Fracture Height	[m]	20

Inspecting results

The following figure shows the distribution of $\sigma_{zz}$ at $t=200 s$ within the computational domain..

../../../../../../_images/szz2.png — Fig. 35 Simulation result of $\sigma_{zz}$ at $t=200 s$ 

First, by running the query script

python ./hydrofractureQueries.py pknViscosityDominated

the HDF5 output is postprocessed and temporal evolution of fracture characterisctics (fluid pressure and fracture width at fluid inlet and fracure half length) are saved into a txt file model-results.txt, which can be used for verification and visualization:

[['      time', '  pressure', '  aperture', '    length']]
         2 1.413e+06 0.0006093       5.6
         4 1.174e+06 0.0007132       8.4
         6 1.077e+06 0.0007849      10.8
         8 1.044e+06 0.0008482      12.8
        10 1.047e+06 0.0009098      14.8

Note: GEOS python tools geosx_xml_tools should be installed to run the query script (See Python Tools Setup for details).

Next, figure below shows the comparisons between the results from GEOS simulations (markers) and the corresponding analytical solutions (curves) for the example with viscosity-storage dominated assumptions, which is generated using the visualization script:

python ./pknViscosityDominatedFigure.py

The evolution in time of the fracture half-length, the near-wellbore fracture aperture, and the fluid pressure all correlate well with the analytical solutions.

(Source code)

../../../../../../_images/pknViscosityDominatedFigure.png

To go further

Feedback on this example

For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.