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" (value 1 stands for a sum of all fluxes, 2 adds per-flux details and 3 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:

pic1 pic2

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