Hydraulic Fracturing
Context
In this example, we use a fully coupled hydrofracture solver from GEOS to solve for the propagation of a single fracture within a reservoir with heterogeneous in-situ properties. Advanced xml features will be used throughout the example.
Objectives
At the end of this example you will know:
how to use multiple solvers for hydraulic fracturing problems,
how to specify pre-existing fractures and where new fractures can develop,
how to construct a mesh with bias,
how to specify heterogeneous in-situ properties and initial conditions,
how to use parameters, symbolic math, and units in xml files.
Input files
This example uses a set of input files and table files located at:
inputFiles/hydraulicFracturing
Because the input files use the advanced xml features, they must be preprocessed using the geosx_xml_tools package. If you have not already done so, setup these features by following the instructions here: Advanced XML Features .
Description of the case
Here, our goal is to demonstrate how hydraulic fractures are modeled in a typical environment. The in-situ properties and initial conditions are based upon a randomly generated, fractal, 1D layer-cake model.
The inputs for this case are contained inside a case-specific (heterogeneousInSitu_benchmark.xml) and base (heterogeneousInSitu_base.xml) XML files.
The tables directory contains the pre-constructed geologic model.
This example will first focus on the case-specific input file, which contains the key parameter definitions, then consider the base xml file.
Included: including external xml files
At the head of the case-specific xml file is a block that will instruct GEOS to include an external file. In our case, this points to the base hydraulic fracturing input file.
<Included>
<File
name="./heterogeneousInSitu_base.xml"/>
</Included>
Parameters: defining variables to be used throughout the file
The Parameters block defines a series of variables that can be used throughout the input file.
These variables allow a given input file to be easily understood and/or modified for a specific environment, even by non-expert users. Parameters are specified in pairs of names and values.
The names should only contain alphanumeric characters and underlines.
The values can contain any type (strings, doubles, etc.).
Parameters can be used throughout the input file (or an included input file) by placing them in-between dollar signs.
Barring any circular-definition errors, parameters can be used within other parameters.
For example, see the parameter mu_upscaled.
The value of this parameter is a symbolic expression, which is denoted by the surrounding back-ticks, and is dependent upon two other parameters.
During pre-processing, geosx_xml_tools will substitute the parameter definitions, and evaluate the symbolic expression using a python-derived syntax.
A number of the input parameters include optional unit definitions, which are denoted by the square brackets following a value.
For example, the parameter t_max is used to set the maximum time for the simulation to 20 minutes.
<Parameters>
<!-- Use the swarm upscaling law -->
<Parameter
name="Nperf"
value="5"/>
<Parameter
name="Nswarm"
value="5"/>
<Parameter
name="mu_init"
value="0.001"/>
<Parameter
name="K_init"
value="1e6"/>
<Parameter
name="mu_upscaled"
value="`$mu_init$*($Nswarm$**2)`"/>
<Parameter
name="K_upscaled"
value="`$K_init$*($Nswarm$**0.5)`"/>
<Parameter
name="ContactStiffness"
value="1e10"/>
<!-- Event timing -->
<Parameter
name="pump_start"
value="1 [min]"/>
<Parameter
name="pump_ramp"
value="5 [s]"/>
<Parameter
name="pump_ramp_dt_limit"
value="0.2 [s]"/>
<Parameter
name="dt_max"
value="30 [s]"/>
<Parameter
name="t_max"
value="20 [min]"/>
<!-- Etc. -->
<Parameter
name="table_root"
value="./tables"/>
<Parameter
name="t_allocation"
value="28 [min]"/>
</Parameters>
Mesh with biased boundaries
The mesh block for this example uses a biased mesh along the simulation boundaries to reduce the size of the problem, while maintaining the desired spatial extents. For a given element region with bias, the left-hand element will be x% smaller and the right-hand element will be x% larger than the average element size. Along the x-axis of the mesh, we select a value of zero for the first region to indicate that we want a uniform-sized mesh, and we select a bias of -0.6 for the second region to indicate that we want the element size to smoothly increase in size as it moves in the +x direction. The other dimensions of the mesh follow a similar pattern.
<Mesh>
<InternalMesh
name="mesh1"
xCoords="{ 0, 200, 250 }"
yCoords="{ -100, 0, 100 }"
zCoords="{ -150, -100, 0, 100, 150 }"
nx="{ 50, 5 }"
ny="{ 10, 10 }"
nz="{ 5, 25, 25, 5 }"
xBias="{ 0, -0.6 }"
yBias="{ 0.6, -0.6 }"
zBias="{ 0.6, 0, 0, -0.6 }"
cellBlockNames="{ cb1 }"
elementTypes="{ C3D8 }"/>
</Mesh>
Defining a fracture nodeset
For this example, we want to propagate a single hydraulic fracture along the plane defined by y = 0. To achieve this, we need to define three nodesets:
source_a: The location where we want to inject fluid. Typically, we want this to be a single face in the x-z plane.
perf_a: This is the initial fracture for the simulation. This nodeset needs to be at least two-faces wide in the x-z plane (to represent the fracture at least one internal node needs to be open).
fracturable_a: This is the set of faces where we will allow the fracture to grow. For a problem where we expect the fracture to curve out of the plane defined by y = 0 , this could be replaced.
<Geometry>
<Box
name="source_a"
xMin="{ -0.1, -0.1, -0.1 }"
xMax="{ 4.1, 0.1, 4.1 }"/>
<Box
name="perf_a"
xMin="{ -4.1, -0.1, -4.1 }"
xMax="{ 4.1, 0.1, 4.1 }"/>
<ThickPlane
name="fracturable_a"
normal="{ 0, 1, 0 }"
origin="{ 0, 0, 0 }"
thickness="0.1"/>
</Geometry>
Boundary conditions
The boundary conditions for this problem are defined in the case-specific and the base xml files. The case specific block includes four instructions:
frac: this marks the initial perforation.
separableFace: this marks the set of faces that are allowed to break during the simulation.
waterDensity: this initializes the fluid in the perforation.
sourceTerm: this instructs the code to inject fluid into the source_a nodeset. Note the usage of the symbolic expression and parameters in the scale. This boundary condition is also driven by a function, which we will define later.
<FieldSpecifications>
<!-- Fracture-related nodesets -->
<FieldSpecification
name="frac"
fieldName="ruptureState"
initialCondition="1"
objectPath="faceManager"
scale="1"
setNames="{ perf_a }"/>
<FieldSpecification
name="separableFace"
fieldName="isFaceSeparable"
initialCondition="1"
objectPath="faceManager"
scale="1"
setNames="{ fracturable_a }"/>
<!-- Fluid Initial Conditions -->
<FieldSpecification
name="waterDensity"
fieldName="water_density"
initialCondition="1"
objectPath="ElementRegions"
scale="1000"
setNames="{ perf_a }"/>
<!-- Fluid BC's -->
<!-- Note: the units of the source flux BC are in kg/s not m3/s -->
<SourceFlux
name="sourceTerm"
objectPath="ElementRegions/Fracture"
scale="`-1000.0/$Nperf$`"
functionName="flow_rate"
setNames="{ source_a }"/>
</FieldSpecifications>
The base block includes instructions to set the initial in-situ properties and stresses. It is also used to specify the external mechanical boundaries on the system. In this example, we are using roller-boundary conditions (zero normal-displacement). Depending upon how close they are to the fracture, they can significantly affect its growth. Therefore, it is important to test whether the size of the model is large enough to avoid this.
<FieldSpecifications>
<FieldSpecification
name="bulk_modulus"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions"
fieldName="rock_bulkModulus"
functionName="bulk_modulus"
scale="1.0"/>
<FieldSpecification
name="shear_modulus"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions"
fieldName="rock_shearModulus"
functionName="shear_modulus"
scale="1.0"/>
<FieldSpecification
name="sigma_xx"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions"
fieldName="rock_stress"
component="0"
functionName="sigma_xx"
scale="1.0"/>
<FieldSpecification
name="sigma_yy"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions"
fieldName="rock_stress"
component="1"
functionName="sigma_yy"
scale="1.0"/>
<FieldSpecification
name="sigma_zz"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions"
fieldName="rock_stress"
component="2"
functionName="sigma_zz"
scale="1.0"/>
<!-- Mechanical BC's -->
<FieldSpecification
name="x_constraint"
component="0"
fieldName="totalDisplacement"
objectPath="nodeManager"
scale="0.0"
setNames="{ xneg, xpos }"/>
<FieldSpecification
name="y_constraint"
component="1"
fieldName="totalDisplacement"
objectPath="nodeManager"
scale="0.0"
setNames="{ yneg, ypos }"/>
<FieldSpecification
name="z_constraint"
component="2"
fieldName="totalDisplacement"
objectPath="nodeManager"
scale="0.0"
setNames="{ zneg, zpos }"/>
</FieldSpecifications>
Coupled hydraulic fracturing solver
The Solvers block is located in the base xml file.
Note that the gravityVector attribute indicates that we are applying gravity in the z-direction in this problem.
Similar to other coupled physics solvers, the Hydrofracture solver is specified in three parts:
Hydrofracture: this is the primary solver, which will be called by the event manager. Two of its key attributes are the names of the dependent solid and fluid solvers.
SolidMechanicsLagrangianSSLE: this is the solid mechanics solver.
SinglePhaseFVM: this is the fluid solver.
The final solver present in this example is the SurfaceGenerator, which manages how faces in the model break.
<Solvers
gravityVector="{ 0.0, 0.0, -9.81 }">
<Hydrofracture
name="hydrofracture"
solidSolverName="lagsolve"
flowSolverName="SinglePhaseFlow"
surfaceGeneratorName="SurfaceGen"
logLevel="1"
targetRegions="{ Fracture }"
contactRelationName="fractureContact"
maxNumResolves="2">
<NonlinearSolverParameters
newtonTol="1.0e-5"
newtonMaxIter="40"
lineSearchMaxCuts="3"/>
<LinearSolverParameters
logLevel="1"
solverType="gmres"
preconditionerType="mgr"/>
</Hydrofracture>
<SolidMechanicsLagrangianSSLE
name="lagsolve"
logLevel="1"
timeIntegrationOption="QuasiStatic"
discretization="FE1"
targetRegions="{ Domain, Fracture }"
contactRelationName="fractureContact">
<NonlinearSolverParameters
newtonTol="1.0e-6"/>
<LinearSolverParameters
solverType="gmres"
krylovTol="1.0e-10"/>
</SolidMechanicsLagrangianSSLE>
<SinglePhaseFVM
name="SinglePhaseFlow"
logLevel="1"
discretization="singlePhaseTPFA"
targetRegions="{ Fracture }">
<NonlinearSolverParameters
newtonTol="1.0e-5"
newtonMaxIter="10"/>
<LinearSolverParameters
solverType="gmres"
krylovTol="1.0e-12"/>
</SinglePhaseFVM>
<SurfaceGenerator
name="SurfaceGen"
targetRegions="{ Domain }"
rockToughness="$K_upscaled$"
nodeBasedSIF="1"
mpiCommOrder="1"/>
</Solvers>
Events
Rather than explicitly specify the desired timestep behavior, this example uses a flexible approach for timestepping.
The hydrofracture solver is applied in the solverApplications event, which request a maxEventDt = 30 s.
To maintain stability during the critical early phase of the model, we delay turning on the pump by pump_start.
We then use the pumpStart event to limit the further limit the timestep to pump_ramp_dt_limit as the fracture experiences rapid development (pump_start to pump_start + pump_ramp).
Note that while this event does not have a target, it can still influence the time step behavior.
After this period, the hydraulic fracture solver will attempt to increase / decrease the requested timestep to maintain stability.
Other key events in this problem include:
preFracture: this calls the surface generator at the beginning of the problem and helps to initialize the fracture.
outputs_vtk and outputs_silo: these produces output vtk and silo files.
restarts (inactive): this is a HaltEvent, which tracks the external clock. When the runtime exceeds the specified value (here $t_allocation$=28 minutes), the code will call the target (which writes a restart file) and instruct the code to exit.
<Events
maxTime="$t_max$"
logLevel="1">
<!-- Generate the initial fractures -->
<SoloEvent
name="preFracture"
target="/Solvers/SurfaceGen"/>
<!-- Primary outputs -->
<PeriodicEvent
name="outputs_vtk"
timeFrequency="1 [min]"
targetExactTimestep="0"
target="/Outputs/vtkOutput"/>
<PeriodicEvent
name="outputs_silo"
timeFrequency="1 [min]"
targetExactTimestep="0"
target="/Outputs/siloOutput"/>
<!-- Apply the hydrofracture solver -->
<PeriodicEvent
name="solverApplications"
maxEventDt="$dt_max$"
target="/Solvers/hydrofracture"/>
<!-- Limit dt during the pump ramp-up -->
<PeriodicEvent
name="pumpStart"
beginTime="$pump_start$"
endTime="`$pump_start$+$pump_ramp$`"
maxEventDt="$pump_ramp_dt_limit$"/>
<!-- Watch the wall-clock, write a restart, and exit gracefully if necessary -->
<!-- <HaltEvent
name="restarts"
maxRuntime="$t_allocation$"
target="/Outputs/restartOutput"/> -->
</Events>
Functions to set in-situ properties
The function definitions are in the base xml file, and rely upon the files in the tables directory. The functions in this example include the flow rate over time, the in-situ principal stress, and the bulk/shear moduli of the rock. Note the use of the table_root parameter, which contains the root path to the table files.
The flow_rate TableFunction is an example of a 1D function. It has a single input, which is time. The table is defined using a single coordinateFile:
0.00000e+00
6.00000e+01
1.20000e+02
3.72000e+03
3.78000e+03
1.00000e+09
And a single voxelFile:
0.00000e+00
0.00000e+00
5.00000e-02
5.00000e-02
0.00000e+00
0.00000e+00
Given the specified linear interpolation method, these values define a simple trapezoidal function.
Note: since this is a 1D table, these values could alternately be given within the xml file using the coordinates and values attributes.
The sigma_xx TableFunction is an example of a 3D function.
It uses elementCenter as its input, which is a vector.
It is specified using a set of three coordinate files (one for each axis), and a single voxel file.
The geologic model in this example is a layer-cake, which was randomly generated, so the size of the x and y axes are 1.
The interpolation method used here is upper, so the values in the table indicate those at the top of each layer.
<Functions>
<!-- Pumping Schedule -->
<TableFunction
name="flow_rate"
inputVarNames="{ time }"
coordinateFiles="{ $table_root$/flowRate_time.csv }"
voxelFile="$table_root$/flowRate.csv"/>
<!-- Geologic Model -->
<TableFunction
name="sigma_xx"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/sigma_xx.csv"
interpolation="upper"/>
<TableFunction
name="sigma_yy"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/sigma_yy.csv"
interpolation="upper"/>
<TableFunction
name="sigma_zz"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/sigma_zz.csv"
interpolation="upper"/>
<TableFunction
name="init_pressure"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/porePressure.csv"
interpolation="upper"/>
<TableFunction
name="bulk_modulus"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/bulkModulus.csv"
interpolation="upper"/>
<TableFunction
name="shear_modulus"
inputVarNames="{ elementCenter }"
coordinateFiles="{ $table_root$/x.csv, $table_root$/y.csv, $table_root$/z.csv }"
voxelFile="$table_root$/shearModulus.csv"
interpolation="upper"/>
<TableFunction
name="apertureTable"
coordinates="{ -1.0e-3, 0.0 }"
values="{ 1.0e-6, 1.0e-4 }"/>
</Functions>
Running GEOS
Assuming that the preprocessing tools have been correctly installed (see Advanced XML Features ), there will be a script in the GEOS build/bin directory called geosx_preprocessed. Replacing geosx with geosx_preprocessed in an input command will automatically apply the preprocessor and send the results to GEOS.
Before beginning, we reccomend that you make a local copy of the example and its tables. Because we are using advanced xml features in this example, the input file must be pre-processed before running. For example, this will run the code on a debug partition using a total of 36 cores.
cp -r examples/hydraulicFracturing ./hf_example
cd hf_example
srun -n 36 -ppdebug geosx_preprocessed -i heterogeneousInSitu_benchmark.xml -x 6 -y 2 -z 3 -o hf_results
Note that as part of the automatic preprocessing step a compiled xml file is written to the disk (by default ‘[input_name].preprocessed’). When developing an xml with advanced features, we reccomend that you check this file to ensure its accuracy.
Inspecting results
In the above example, we requested vtk- and silo-format output files every minute. We can therefore import these into VisIt, Paraview, or python and visualize the outcome. The following figure shows the extents of the generated fracture over time:
Notes for visualization tools:
In Visit, we currently recommend that you look at the silo-format files (due to a compatibility issue with vtk)
In Paraview, you may need to use the Multi-block Inspector (on the right-hand side of the screen by default) to limit the visualization to the fracture. In addition, the Properties inspector (on the left-hand side of the sceen by default) may not include some of the parameters present on the fracture. Instead, we recommend that you use the property dropdown box at the top of the screen.
Because we did not explicitly specify any fracture barriers in this example, the fracture dimensions are controlled by the in-situ stresses. During the first couple of minutes of growth, the fracture quickly reaches its maximum/minimum height, which corresponds to a region of low in-situ minimum stress.
The following figures show the aperture and pressure of the hydraulic fracture (near the source) over time:
Modifying Parameters Via the Command-Line
The advanced xml feature preprocessor allows parameters to be set or overriden by specifying any number of -p name value arguments on the command-line. Note that if the parameter value has spaces, it needs to be enclosed by quotation marks.
To illustrate this feature, we can re-run the previous analysis with viscosity increased from 1 cP to 5 cP:
srun -n 36 -ppdebug geosx_preprocessed -i heterogeneousInSitu_benchmark.xml -p mu 0.005 -x 6 -y 2 -z 3 -o hf_results_lower_mu
To go further
Feedback on this example
This concludes the hydraulic fracturing example. For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.
For more details
More on advanced xml features, please see Advanced XML Features.
More on functions, please see Functions.
More on biased meshes, please see Mesh Bias.