Multiphase Flow
Context
In this example, we set up a multiphase, multicomponent test case (see Compositional Multiphase Flow Solver). The permeability field corresponds to the two bottom layers (layers 84 and 85) of the SPE10 test case. The thermodynamic behavior of the fluid mixture is specified using a simple immiscible two-phase (Dead-Oil) model. Injection and production are simulated using boundary conditions.
Objective
The objective of this example is to review the main elements of a simple two-phase simulation in GEOS, including:
the compositional multiphase flow solver,
the multiphase constitutive models,
the specifications of multiphase boundary conditions.
Input file
This example is based on the XML file located at
inputFiles/compositionalMultiphaseFlow/benchmarks/SPE10/deadOilSpe10Layers84_85_base_iterative.xml
The XML file considered here follows the typical structure of the GEOS input files:
Multiphase flow solver
In GEOS, the setup of a multiphase simulation starts in the Solvers XML block of the input file. This example relies on a solver of type CompositionalMultiphaseFVM that implements a fully implicit finite-volume scheme based on the standard two-point approximation of the flux (TPFA). More information on this solver can be found at Compositional Multiphase Flow Solver.
Let us have a closer look at the Solvers XML block displayed below.
The solver has a name (here, compflow
) that can be chosen by the user and is not imposed by GEOS.
Note that this name is used in the Events XML block to trigger the application of the solver.
Using the targetRegions
attribute, the solver defines target regions on which it is applied.
In this example, there is only one region, named reservoir
.
The CompositionalMultiphaseFVM block contains two important sub-blocks, NonlinearSolverParameters and LinearSolverParameters.
In NonlinearSolverParameters, we can tune the nonlinear tolerance, the application of the linear search algorithm, and the heuristics that control time step sizes.
In LinearSolverParameters, we can specify the linear tolerance, the type of (direct or iterative) linear solver, and the
type of preconditioner, if any.
For large multiphase flow problems, we recommend using an iterative linear solver (solverType="gmres"
or solverType="fgmres"
) combined
with the multigrid reduction (MGR) preconditioner (preconditionerType="mgr"
). More information about the MGR preconditioner can be found in Linear Solvers.
Note
For non-trivial simulations, we recommend setting the initialDt
attribute to a small value (relative to the time scale of the problem) in seconds. If the simulation appears to be slow, use logLevel="1"
in CompositionalMultiphaseFVM to detect potential Newton convergence problems. If the Newton solver struggles, please set lineSearchAction="Attempt"
in NonlinearSolverParameters. If the Newton convergence is good, please add logLevel="1"
in the LinearSolverParameters block to detect linear solver problems, especially if an iterative linear solver is used.
Note
To use the linear solver options of this example, you need to ensure that GEOS is configured to use the Hypre linear solver package.
<Solvers>
<CompositionalMultiphaseFVM
name="compflow"
logLevel="1"
discretization="fluidTPFA"
targetRegions="{ reservoir }"
temperature="300"
useMass="1"
initialDt="1e3"
maxCompFractionChange="0.1">
<NonlinearSolverParameters
newtonTol="1.0e-4"
newtonMaxIter="40"
maxTimeStepCuts="10"
lineSearchAction="None"/>
<LinearSolverParameters
solverType="fgmres"
preconditionerType="mgr"
krylovTol="1.0e-5"/>
</CompositionalMultiphaseFVM>
</Solvers>
Mesh
In this simulation, we define a mesh generated internally using the InternalMesh generator, illustrated in the first tutorial (Tutorial 1: First Steps). The mesh dimensions and cell sizes are chosen to be those specified in the the two bottom layers of the SPE10 test case. The mesh dimensions are specified in meters.
<Mesh>
<InternalMesh
name="mesh"
elementTypes="{ C3D8 }"
xCoords="{ 0, 365.76 }"
yCoords="{ 0, 670.56 }"
zCoords="{ 0, 1.22 }"
nx="{ 60 }"
ny="{ 220 }"
nz="{ 2 }"
cellBlockNames="{ block }"/>
</Mesh>
Geometry
The Geometry XML block is used to select cells where boundary conditions are applied. To mimic the setup of the original SPE10 test case, we place a box-shaped source term in the middle of the domain and a sink term in each corner. The specification of the boundary conditions applied to the selected mesh cells is done in the FieldSpecifications block of the XML file using the names of the boxes defined here.
<Geometry>
<Box
name="source"
xMin="{ 182.85, 335.25, -0.01 }"
xMax="{ 189.00, 338.35, 2.00 }"/>
<Box
name="sink1"
xMin="{ -0.01, -0.01, -0.01 }"
xMax="{ 6.126, 3.078, 2.00 }"/>
<Box
name="sink2"
xMin="{ -0.01, 667.482, -0.01 }"
xMax="{ 6.126, 670.60, 2.00 }"/>
<Box
name="sink3"
xMin="{ 359.634, -0.01, -0.01 }"
xMax="{ 365.8, 3.048, 2.00 }"/>
<Box
name="sink4"
xMin="{ 359.634, 667.482, -0.01 }"
xMax="{ 365.8, 670.60, 2.00 }"/>
</Geometry>
Events
In the Events XML block of this example, we specify two types of PeriodicEvents: one for solver application, and one for result output.
The periodic event named solverApplications
triggers the application of the solver on its target region.
This event must point to the solver by name.
In this example, the name of the solver is compflow
and was defined in the Solvers block.
The time step is initialized using the initialDt
attribute of the flow solver.
Then, if the solver converges in less than a certain number of nonlinear iterations (by default, 40% of the
maximum number of nonlinear iterations), the time step will increase until it reaches the maximum
time step size specified with maxEventDt
.
If convergence fails, the time step will be cut. Parameters defining the time stepping strategy
can be adjusted by the user in the flow solver block.
All times are in seconds.
The output event instructs GEOS to write out results in a given format at the frequency specified by the attribute
timeFrequency
.
Here, we output the results using the VTK format (see Tutorial 2: External Meshes
for a tutorial that uses the Silo output file format).
Here, by setting targetExactTimestep=1
, we make sure that the time stepping strategy is constrained
to generate an output exactly at the chosen frequency, even if solvers could take larger time steps.
In the target
attribute, we use the name defined in the VTK XML tag
inside the Output XML section, as documented at the end of this example (here, vtkOutput
).
More information about Events can be found at Event Management.
<Events
maxTime="2e6">
<PeriodicEvent
name="outputs"
timeFrequency="5e5"
targetExactTimestep="1"
target="/Outputs/vtkOutput"/>
<PeriodicEvent
name="solverApplications"
maxEventDt="5e5"
target="/Solvers/compflow"/>
<PeriodicEvent name="sourceFluxStatsEvent"
timeFrequency="5e5"
target="/Tasks/fluxesStats" />
<PeriodicEvent
name="restarts"
timeFrequency="1e6"
targetExactTimestep="0"
target="/Outputs/restartOutput"/>
</Events>
Numerical methods
In the NumericalMethods XML block, we select a two-point flux approximation (TPFA) finite-volume scheme to discretize the governing equations on the reservoir mesh.
<NumericalMethods>
<FiniteVolume>
<TwoPointFluxApproximation
name="fluidTPFA"/>
</FiniteVolume>
</NumericalMethods>
Reservoir region
In the ElementRegions XML block, we define a CellElementRegion named reservoir
corresponding to the
reservoir mesh.
cellBlocks
is set to { * }
to automatically include all cells of the mesh.
The CellElementRegion must point to the constitutive models that update
the dynamic rock and fluid properties in the reservoir cells.
The names fluid
, rock
, and relperm
used for this in the materialList
correspond to the Constitutive blocks.
<ElementRegions>
<CellElementRegion
name="reservoir"
cellBlocks="{ * }"
materialList="{ fluid, rock, relperm }"/>
</ElementRegions>
Constitutive models
For a simulation performed with the CompositionalMultiphaseFVM physics solver, at least four types of constitutive models must be specified in the Constitutive XML block:
a fluid model describing the thermodynamics behavior of the fluid mixture,
a relative permeability model,
a rock permeability model,
a rock porosity model.
All these models use SI units exclusively. A capillary pressure model can also be specified in this block but is omitted here for simplicity.
Here, we use a simplified fluid model corresponding to an immiscible two-phase model (Dead Oil) by using the XML tag DeadOilFluid. Other fluid models can be used with the CompositionalMultiphaseFVM solver, as explained in Fluid Models.
We define a relative permeability model with the tag BrooksCoreyRelativePermeability. A list of available relative permeability models can be found at Relative Permeability Models.
The properties are chosen to match those of the original SPE10 test case.
Note
The names and order of the phases listed for the attribute phaseNames
must be identical in the fluid model
(here, DeadOilFluid) and the relative permeability model (here, BrooksCoreyRelativePermeability).
Otherwise, GEOS will throw an error and terminate.
We introduce models to define rock compressibility and permeability. This step is similar to what is described in the previous examples (see for instance Tutorial 1: First Steps).
The attribute name
of the constitutive models defined here
must be used in the ElementRegions and Solvers XML blocks to point the element
regions and the physics solvers to their respective constitutive models.
<Constitutive>
<DeadOilFluid
name="fluid"
phaseNames="{ oil, water }"
surfaceDensities="{ 800.0, 1022.0 }"
componentMolarWeight="{ 114e-3, 18e-3 }"
hydrocarbonFormationVolFactorTableNames="{ B_o_table }"
hydrocarbonViscosityTableNames="{ visc_o_table }"
waterReferencePressure="30600000.1"
waterFormationVolumeFactor="1.03"
waterCompressibility="0.00000000041"
waterViscosity="0.0003"/>
<CompressibleSolidConstantPermeability
name="rock"
solidModelName="nullSolid"
porosityModelName="rockPorosity"
permeabilityModelName="rockPerm"/>
<NullModel
name="nullSolid"/>
<PressurePorosity
name="rockPorosity"
defaultReferencePorosity="0.1"
referencePressure="1.0e7"
compressibility="1e-10"/>
<BrooksCoreyRelativePermeability
name="relperm"
phaseNames="{ oil, water }"
phaseMinVolumeFraction="{ 0.0, 0.0 }"
phaseRelPermExponent="{ 2.0, 2.0 }"
phaseRelPermMaxValue="{ 1.0, 1.0 }"/>
<ConstantPermeability
name="rockPerm"
permeabilityComponents="{ 1.0e-14, 1.0e-14, 1.0e-18 }"/>
</Constitutive>
Initial and boundary conditions
In the FieldSpecifications section, we define initial and boundary conditions, as well as geological properties (porosity, permeability). Everything is specified using SI units. We focus on the specification of the initial and boundary conditions for a simulation performed with the CompositionalMultiphaseFVM solver. We refer to Tutorial 1: First Steps for a more general discussion on the FieldSpecification XML blocks.
For a simulation performed with the CompositionalMultiphaseFVM solver,
we set the initial pressure and the initial global component
fractions (here, the oil and water component fractions).
The component
attribute of the FieldSpecification XML block must use the
order in which the phaseNames
have been defined in the DeadOilFluid
XML block. In other words, component=0
is used to initialize the oil
global component fraction and component=1
is used to initialize the water global
component fraction, because we previously set phaseNames="{oil, water}"
in the DeadOilFluid XML block.
To specify sink terms, we use the FieldSpecification mechanism in a similar fashion to impose the sink pressure and composition. This mimics a pressure-controlled well (before breakthrough). To specify the source term, we use a SourceFlux block to impose a fixed mass injection rate of component 1 (water) to mimic a rate-controlled well.
<FieldSpecifications>
<FieldSpecification
name="permx"
component="0"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="rockPerm_permeability"
functionName="permxFunc"
scale="9.869233e-16"/>
<FieldSpecification
name="permy"
component="1"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="rockPerm_permeability"
functionName="permyFunc"
scale="9.869233e-16"/>
<FieldSpecification
name="permz"
component="2"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="rockPerm_permeability"
functionName="permzFunc"
scale="9.869233e-16"/>
<FieldSpecification
name="referencePorosity"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="rockPorosity_referencePorosity"
functionName="poroFunc"
scale="1.0"/>
<FieldSpecification
name="initialPressure"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="pressure"
scale="4.1369e+7"/>
<FieldSpecification
name="initialComposition_oil"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="globalCompFraction"
component="0"
scale="0.9995"/>
<FieldSpecification
name="initialComposition_water"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/reservoir/block"
fieldName="globalCompFraction"
component="1"
scale="0.0005"/>
<SourceFlux
name="sourceTerm"
objectPath="ElementRegions/reservoir/block"
scale="-0.07279"
component="1"
setNames="{ source }"/>
<FieldSpecification
name="sinkPressure"
setNames="{ sink1, sink2, sink3, sink4 }"
objectPath="ElementRegions/reservoir/block"
fieldName="pressure"
scale="2.7579e+7"/>
<FieldSpecification
name="sinkComposition_oil"
setNames="{ sink1, sink2, sink3, sink4 }"
objectPath="ElementRegions/reservoir/block"
fieldName="globalCompFraction"
component="0"
scale="0.9995"/>
<FieldSpecification
name="sinkComposition_water"
setNames="{ sink1, sink2, sink3, sink4 }"
objectPath="ElementRegions/reservoir/block"
fieldName="globalCompFraction"
component="1"
scale="0.0005"/>
</FieldSpecifications>
Source flux statistics
GEOS reports source flux statistics to help monitor the overall mass balance and flow rates of a group of wells defined using SourceFlux
.
This functionality is implemented through a SourceFluxStatistics
Task that reports statistics at regular intervals via Events to:
a CSV file (comma-separated values) that reports all fluxes and regions details and sum, when specifying
writeCSV="1"
.console log output, when specifying
logLevel="1"
(value1
stands for a sum of all fluxes,2
adds per-flux details and3
adds per-region details).
The group of target SourceFlux is defined in fluxNames
by specifying their names, or by using *
to report all source fluxes.
The task requires a flowSolverName
referring to the solver used for flow calculations.
<Tasks>
<SourceFluxStatistics name="fluxesStats"
flowSolverName="compflow"
fluxNames="{*}"
logLevel="2"
writeCSV="1" />
</Tasks>
To output these statistics during the simulation, an event must be defined (previously shown sourceFluxStatsEvent
).
In addition to the complete information reported in the fluxesStats.csv
file, the source flux values are reported in table format in the console log output:
-----------------------------------------------------------------------------------------------------
| fluxesStats, flux statistics for: sourceTerm |
|---------------------------------------------------------------------------------------------------|
| Flux(es) | Region | Element Count | Prod. mass [kg] | Prod. rate [kg/s] |
|--------------|---------------|-----------------|----------------------------|---------------------|
| sourceTerm | all_regions | 2 | [0, -57048.521447934894] | [0, -0.07279] |
| flux_set | all_regions | 2 | [0, -57048.521447934894] | [0, -0.07279] |
-----------------------------------------------------------------------------------------------------
Output
In this section, we request an output of the results in VTK format. Note that the name defined here must match the names used in the Events XML block to define the output frequency.
<Outputs>
<VTK
name="vtkOutput"/>
<Restart
name="restartOutput"/>
</Outputs>
All elements are now in place to run GEOS.
Running GEOS
The simulation can be launched with:
geosx -i deadOilSpe10Layers84_85_benchmark.xml
The first few lines appearing to the console are indicating that the XML elements are read and registered correctly:
Adding Solver of type CompositionalMultiphaseFVM, named compflow
Adding Mesh: InternalMesh, mesh
Adding Geometric Object: Box, source
Adding Geometric Object: Box, sink1
Adding Geometric Object: Box, sink2
Adding Geometric Object: Box, sink3
Adding Geometric Object: Box, sink4
Adding Event: PeriodicEvent, outputs
Adding Event: PeriodicEvent, solverApplications
TableFunction: permxFunc
TableFunction: permyFunc
TableFunction: permzFunc
TableFunction: poroFunc
TableFunction: B_o_table
TableFunction: visc_o_table
Adding Output: VTK, vtkOutput
Adding Object CellElementRegion named region from ObjectManager::Catalog.
region/block/fluid is allocated with 1 quadrature points.
region/block/rock is allocated with 1 quadrature points.
aaregion/block/relperm is allocated with 1 quadrature points.
At this point, we are done with the case set-up and the code steps into the execution of the simulation itself:
Time: 0s, dt:1000s, Cycle: 0
Attempt: 0, NewtonIter: 0
( Rfluid ) = (2.28e+00) ; ( R ) = ( 2.28e+00 ) ;
Attempt: 0, NewtonIter: 1
( Rfluid ) = (8.83e-03) ; ( R ) = ( 8.83e-03 ) ;
Last LinSolve(iter,res) = ( 2, 2.74e-03 ) ;
Attempt: 0, NewtonIter: 2
( Rfluid ) = (8.86e-05) ; ( R ) = ( 8.86e-05 ) ;
Last LinSolve(iter,res) = ( 2, 8.92e-03 ) ;
compflow: Max phase CFL number: 0.00399585
compflow: Max component CFL number: 0.152466
compflow: Newton solver converged in less than 16 iterations, time-step required will be doubled.
Visualization
A file compatible with Paraview is produced in this example. It is found in the output folder, and usually has the extension .pvd. More details about this file format can be found here. We can load this file into Paraview directly and visualize results:
To go further
Feedback on this example
This concludes the example on setting up an immiscible two-phase flow simulation in a channelized permeability field. For any feedback on this example, please submit a GitHub issue on the project’s GitHub page.
For more details
A complete description of the reservoir flow solver is found here: Compositional Multiphase Flow Solver.
The available constitutive models are listed at Constitutive Models.