Compositional Multiphase Flow Solver

Introduction

This flow solver is in charge of implementing the finite-volume discretization (mainly, accumulation and flux terms, boundary conditions) of the equations governing compositional multiphase flow in porous media. The present solver can be combined with the Compositional Multiphase Well Solver which handles the discrete multi-segment well model and provides source/sink terms for the fluid flow solver.

Below, we first review the set of Governing Equations, followed by a discussion of the choice of Primary Variables used in the global variable formulation. Then we give an overview of the Discretization and, finally, we provide a list of the solver Parameters and an input Example.

Theory

Governing Equations

Mass Conservation Equations

Mass conservation for component $c$ is expressed as:

$\phi \frac{ \partial }{\partial t} \bigg( \sum_\ell \rho_{\ell} \, y_{c \ell} \, S_{\ell} \bigg) + \nabla \cdot \bigg( \sum_\ell \rho_{\ell} \, y_{c \ell} \, \boldsymbol{u}_{\ell} \bigg) - \sum_\ell \rho_{\ell} \, y_{c \ell} \, q_{\ell} = 0,$

where $\phi$ is the porosity of the medium, $S_{\ell}$ is the saturation of phase $\ell$ , $y_{c \ell}$ is the mass fraction of component $c$ in phase $\ell$ , $\rho_{\ell}$ is the phase density, and $t$ is time. We note that the formulation currently implemented in GEOS is isothermal.

Darcy’s Law

Using the multiphase extension of Darcy’s law, the phase velocity $\boldsymbol{u}_{\ell}$ is written as a function of the phase potential gradient $\nabla \Phi_{\ell}$ :

$\boldsymbol{u}_{\ell} := -\boldsymbol{k} \lambda_{\ell} \nabla \Phi_{\ell} = - \boldsymbol{k} \lambda_{\ell} \big( \nabla (p - P_{c,\ell}) - \rho_{\ell} g \nabla z \big).$

In this equation, $\boldsymbol{k}$ is the rock permeability, $\lambda_{\ell} = k_{r \ell} / \mu_{\ell}$ is the phase mobility, defined as the phase relative permeability divided by the phase viscosity, $p$ is the reference pressure, $P_{c,\ell}$ is the the capillary pressure, $g$ is the gravitational acceleration, and $z$ is depth. The evaluation of the relative permeabilities, capillary pressures, and viscosities is reviewed in the section about Constitutive Models.

Combining the mass conservation equations with Darcy’s law yields a set of $n_c$ equations written as:

$\phi \frac{ \partial }{\partial t} \bigg( \sum_\ell \rho_{\ell} \, y_{c \ell} \, S_{\ell} \bigg) - \nabla \cdot \boldsymbol{k} \bigg( \sum_\ell \rho_{\ell} \, y_{c \ell} \, \lambda_{\ell} \nabla \Phi_{\ell} \bigg) - \sum_\ell \rho_{\ell} \, y_{c \ell} \, q_{\ell} = 0.$

Constraints and Thermodynamic Equilibrium

The volume constraint equation states that the pore space is always completely filled by the phases. The constraint can be expressed as:

$\sum_{\ell} S_{\ell} = 1.$

The system is closed by the following thermodynamic equilibrium constraints:

$f_{c \ell} - f_{c m} = 0.$

where $f_{c \ell}$ is the fugacity of component $c$ in phase $\ell$ . The flash calculations performed to enforce the thermodynamical equilibrium are reviewed in the section about Constitutive Models.

To summarize, the compositional multiphase flow solver assembles a set of $n_c+1$ equations in each element, i.e., $n_c$ mass conservation equations and one volume constraint equation. A separate module discussed in the Constitutive Models is responsible for the enforcement of the thermodynamic equilibrium at each nonlinear iteration.

Number of equations	Equation type
$n_c$	Mass conservation equations
1	Volume constraint

Primary Variables

The variable formulation implemented in GEOS is a global variable formulation based on $n_c+1$ primary variables, namely, one pressure, $p$ , and $n_c$ component densities, $\rho_c$ . By default, we use molar component densities. A flag discussed in the section Parameters can be used to select mass component densities instead of molar component densities.

Number of primary variables	Variable type
1	Pressure
$n_c$	Component densities

Assembling the residual equations and calling the Constitutive Models requires computing the molar component fractions and saturations. This is done with the relationship:

$z_c := \frac{\rho_c}{\rho_T},$

where

$\rho_T := \sum_c \rho_c.$

These secondary variables are used as input to the flash calculations. After the flash calculations, the saturations are computed as:

$S_{\ell} := \nu_{\ell} \frac{ \rho_T }{ \rho_{\ell}},$

where $\nu_{\ell}$ is the global mole fraction of phase $\ell$ and $\rho_{\ell}$ is the molar density of phase $\ell$ . These steps also involve computing the derivatives of the component fractions and saturations with respect to the pressure and component densities.

Discretization

Spatial Discretization

The governing equations are discretized using standard cell-centered finite-volume discretization.

In the approximation of the flux term at the interface between two control volumes, the calculation of the pressure stencil is general and will ultimately support a Multi-Point Flux Approximation (MPFA) approach. The current implementation of the transmissibility calculation is reviewed in the section about Numerical Methods.

The approximation of the dynamic transport coefficients multiplying the discrete potential difference (e.g., the phase mobilities) is performed with a first-order phase-per-phase single-point upwinding based on the sign of the phase potential difference at the interface.

Temporal Discretization

The compositional multiphase solver uses a fully implicit (backward Euler) temporal discretization.

Solution Strategy

The nonlinear solution strategy is based on Newton’s method. At each Newton iteration, the solver assembles a residual vector, $R$ , collecting the $n_c$ discrete mass conservation equations and the volume constraint for all the control volumes.

Parameters

The following attributes are supported:

XML Element: CompositionalMultiphaseFVM

Name	Type	Default	Description
allowLocalCompDensityChopping	integer	1	Flag indicating whether local (cell-wise) chopping of negative compositions is allowed
allowNegativePressure	integer	1	Flag indicating if negative pressure is allowed
cflFactor	real64	0.5	Factor to apply to the CFL condition when calculating the maximum allowable time step. Values should be in the interval (0,1]
contMultiplierDBC	real64	0.5	Factor by which continuation parameter is changed every newton when DBC is used
continuationDBC	integer	1	Flag for enabling continuation parameter
discretization	groupNameRef	required	Name of discretization object (defined in the Numerical Methods) to use for this solver. For instance, if this is a Finite Element Solver, the name of a Finite Element Discretization should be specified. If this is a Finite Volume Method, the name of a Finite Volume Discretization discretization should be specified.
initialDt	real64	1e+99	Initial time-step value required by the solver to the event manager.
isThermal	integer	0	Flag indicating whether the problem is thermal or not.
kappaminDBC	real64	1e-20	Factor that controls how much dissipation is kept in the system when continuation is used
logLevel	integer	0	Sets the level of information to write in the standard output (the console typically). Level 0 outputs no specific information for this solver. Higher levels require more outputs. 1 - Line search information - Solution information (scaling, maximum changes, quality check) - Convergence information - Time step information - Linear solver information - Nonlinear solver information - Solver timers information 2 - The summary of declared fields and coupling
maxAbsolutePressureChange	real64	-1	Maximum (absolute) pressure change in a Newton iteration
maxCompFractionChange	real64	0.5	Maximum (absolute) change in a component fraction in a Newton iteration
maxRelativeCompDensChange	real64	1.79769e+208	Maximum (relative) change in a component density in a Newton iteration
maxRelativePressureChange	real64	0.5	Maximum (relative) change in pressure in a Newton iteration
maxRelativeTemperatureChange	real64	0.5	Maximum (relative) change in temperature in a Newton iteration
maxSequentialCompDensChange	real64	1	Maximum (absolute) component density change in a sequential iteration, used for outer loop convergence check
maxSequentialPressureChange	real64	100000	Maximum (absolute) pressure change in a sequential iteration, used for outer loop convergence check
maxSequentialTemperatureChange	real64	0.1	Maximum (absolute) temperature change in a sequential iteration, used for outer loop convergence check
minCompDens	real64	1e-10	Minimum allowed global component density
minScalingFactor	real64	0.01	Minimum value for solution scaling factor
miscibleDBC	integer	0	Flag for enabling DBC formulation with/without miscibility
name	groupName	required	A name is required for any non-unique nodes
omegaDBC	real64	1	Factor by which DBC flux is multiplied
scalingType	geos_CompositionalMultiphaseFVM_ScalingType	Global	Solution scaling type.Valid options: * Global * Local
solutionChangeScalingFactor	real64	0.5	Damping factor for solution change targets
targetFlowCFL	real64	-1	Target CFL condition CFL condition when computing the next timestep.
targetPhaseVolFractionChangeInTimeStep	real64	0.2	Target (absolute) change in phase volume fraction in a time step
targetRegions	groupNameRef_array	required	Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.
targetRelativeCompDensChangeInTimeStep	real64	1.79769e+308	Target (relative) change in component density in a time step
targetRelativePressureChangeInTimeStep	real64	0.2	Target (relative) change in pressure in a time step (expected value between 0 and 1)
targetRelativeTemperatureChangeInTimeStep	real64	0.2	Target (relative) change in temperature in a time step (expected value between 0 and 1)
temperature	real64	required	Temperature
useDBC	integer	0	Enable Dissipation-based continuation flux
useMass	integer	0	Use mass formulation instead of molar. Warning : Affects SourceFlux rates units.
useNewGravity	integer	0	Flag indicating whether new gravity treatment is used
useSimpleAccumulation	integer	1	Flag indicating whether simple accumulation form is used
useTotalMassEquation	integer	1	Flag indicating whether total mass equation is used
writeLinearSystem	integer	0	Write matrix, rhs, solution to screen ( = 1) or file ( = 2).
LinearSolverParameters	node	unique	XML Element: LinearSolverParameters
NonlinearSolverParameters	node	unique	XML Element: NonlinearSolverParameters

Example

  <Solvers>
    <CompositionalMultiphaseFVM
      name="compflow"
      logLevel="1"
      discretization="fluidTPFA"
      targetRelativePressureChangeInTimeStep="1"
      targetPhaseVolFractionChangeInTimeStep="1"      
      targetRegions="{ Channel }"
      temperature="300">
      <NonlinearSolverParameters
        newtonTol="1.0e-10"
        newtonMaxIter="10"/>
      <LinearSolverParameters
        directParallel="0"/>
    </CompositionalMultiphaseFVM>
  </Solvers>

We refer the reader to Multiphase Flow for a complete tutorial illustrating the use of this solver.