Immiscible Multiphase Flow Solver

Introduction

This flow solver is used to implement the finite-volume discretization for the problem of modeling multiphase flow in porous media under the influence of viscous, gravity, and capillary forces while neglecting miscibility and taking into account rock and fluid compressibility.

In here, we go over the governing equations Governing Equations that covers two different formulation options, followed by the Discretization, and we conclude by providing a list of the solver Parameters and an input Example.

Theory

Governing Equations

Mass Conservation Equations

We consider a two-component system, say gas and water, flow in a compressible porous medium, in which both components can exist only in their corresponding phases of vapor and liquid. The gas and water components are denoted by the subscripts $g$ and $w$ , respectively. Moreover, the liquid, which is the wetting phase, and the vapor, the non-wetting phase, are denoted by the subscripts $\ell$ and $v$ , respectively. The mass conservation laws are expressed as:

$\frac{\partial}{\partial t} (\phi\rho_v S_v) + \nabla \cdot (\rho_v \boldsymbol{u}_v) = \rho_v q_v,$

and

$\frac{\partial}{\partial t} (\phi\rho_\ell S_\ell) + \nabla \cdot (\rho_\ell \boldsymbol{u}_\ell) = \rho_\ell q_\ell,$

where $\phi(\mathbf{x})$ is the porosity of the medium which is a function of pressure, $S_\ell(\mathbf{x},t)$ is the saturation of the phase $\ell$ and similarly for the phase $v$ , and $t$ is the time. The source/sink terms $q_{\ell}$ and $q_{v}$ are positive for injection and negative for production. The phase velocity, $\boldsymbol{u}_\ell$ and $\boldsymbol{u}_v$ , are defined using the multiphase extension of Darcy’s law (conservation of momentum) as

$\boldsymbol{u}_\ell := -k\lambda_\ell(\nabla p_\ell - \rho_\ell g \nabla z),$

and

$\boldsymbol{u}_v := -k\lambda_v(\nabla p_v - \rho_v g \nabla z).$

Here, $k(\mathbf{x})$ is the scalar absolute permeability of the medium, $\lambda_\ell$ is the phase mobility of the liquid phase defined as $k_{r\ell}/\mu_\ell$ , where $k_{r\ell}(\mathbf{x},S_\ell)$ is the phase relative permeability, $\mu_\ell$ is the phase viscosity, and $\rho_{\ell}$ is the phase density. These are also defined similarly for the vapor phase. In both cases we assume that the relative permeabilities are strictly increasing functions of their own saturation. The gravitational acceleration is denoted by $g$ , and the depth by $z$ (positive going downward). The conservation of mass equations are constrained by the volume contraint equation:

$S_{\ell} + S_v = 1,$

Moreover, the capillary pressure constraint relates the two phase pressures with

$P_{c}(S_{\ell}) = p_{v} - p_{\ell}.$

We assume that capillary pressure is a strictly decreasing function of the wetting-phase saturation.

The evaluation of the relative permeabilities, capillary pressures, and viscosities is reviewed in the section about Constitutive Models.

We note that the formulation currently implemented in GEOS is isothermal.

To summarize, the Immiscible multiphase flow solver assembles a set of $n_p+1$ equations in each element, i.e., $n_p$ mass conservation equations and one volume constraint equation.

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

Primary Variables

There are two formulations implemented in GEOS for the Immiscible multiphsae solver and both formulations are based on $n_p+1$ primary variables, namely, one pressure, $p$ , and $n_p$ phase volume fractions, $S_{p}$ .

Number of primary variables	Variable type
1	Pressure
$n_p - 1$	Phase volume fractions

The main formulation is the standard formulation which solves the individual components mass conservation equations. Also, another formulation based on the total mass flux is implemented which is useful for multiple purposes such as hybrid upwinding techniques and sequential finite volume methods. This latter formulation is explained next.

Flow and Transport Equations

To develop this formulation we use a flux approximation as required by the finite-volume numerical solution scheme. Thus, we choose to construct this approximation in fractional flow form, and with this we will be able to show the coupling between the different physical processes. This formulation is obtained by decomposing the governing equations into a flow problem for both phases and a transport problem for one of the two phases. To obtain this decomposition, we use a total-mass balance formulation by summing both components mass conversation equations and then using the mass constraint to result in the following elliptic PDE governing the temporal evolution of the pressure field:

$\frac{\partial}{\partial t}(\phi \rho_t) + \nabla \cdot (\rho_{\ell} \boldsymbol{u}_{\ell} + \rho_{v} \boldsymbol{u}_{v}) = \rho_{\ell} q_{\ell} + \rho_{v} q_{v},$

where

$\rho_t = \rho_{\ell}S_{\ell}+\rho_{v}S_{v}$

and we defined a total mass flux as

$\boldsymbol{U}_T := \rho_{\ell} \boldsymbol{u}_{\ell} + \rho_{v} \boldsymbol{u}_{v}= -k (\rho_{\ell} \lambda_{\ell} + \rho_{v} \lambda_{v}) \nabla p + k ( \lambda_{\ell} \rho^2_{\ell} + \lambda_{v} \rho^2_{v}) g \nabla z + k \rho_{\ell}\lambda_{\ell} \nabla P_{c}.$

Next, the highly nonlinear parabolic transport equation is obtained by using this total mass flux to formally eliminate the pressure variable from the individual components mass conservation equations, yielding

$\frac{\partial}{\partial t}(\phi\rho_v S_v) + \nabla \cdot F_v = \rho_v q_v,$

and

$\frac{\partial}{\partial t}(\phi\rho_\ell S_\ell) + \nabla \cdot F_\ell = \rho_\ell q_\ell,$

where the flow flux for each phase is defined as

${ F_{\ell} := \frac{\rho_\ell \lambda_\ell}{\rho_\ell \lambda_\ell+\rho_v \lambda_v}\boldsymbol{U}_T} + k \frac{\rho_\ell \lambda_\ell\rho_v \lambda_v}{\rho_\ell \lambda_\ell+\rho_v \lambda_v}(\rho_\ell - \rho_v) g\nabla z + k \frac{\rho_\ell \lambda_\ell\rho_v \lambda_v}{\rho_\ell \lambda_\ell+\rho_v \lambda_v} ( \nabla P_{c})$

and

${ F_{v} := \frac{\rho_v \lambda_v}{\rho_\ell \lambda_\ell+\rho_v \lambda_v}\boldsymbol{U}_T} + k \frac{\rho_\ell \lambda_\ell\rho_v \lambda_v}{\rho_\ell \lambda_\ell+\rho_v \lambda_v}(\rho_v - \rho_\ell) g\nabla z - k \frac{\rho_\ell \lambda_\ell\rho_v \lambda_v}{\rho_\ell \lambda_\ell+\rho_v \lambda_v} ( \nabla P_{c})$

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 immiscible 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_p$ discrete mass conservation equations and the volume constraint for all the control volumes.

Parameters

The following attributes are supported:

XML Element: ImmiscibleMultiphaseFlow

Name	Type	Default	Description
allowNegativePressure	integer	0	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]
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.
gravityDensityScheme	geos_GravityDensityScheme	ArithmeticAverage	Scheme for density treatment in gravity
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.
logLevel	integer	0	Sets the level of information to write in the standard output (the console typically). Information output from lower logLevels is added with the desired log level 1 - Linear solver information - Solution information (scaling, maximum changes, quality check) - Convergence information - Time step information - Solver timers information 2 - The summary of declared fields and coupling
maxAbsolutePressureChange	real64	-1	Maximum (absolute) pressure change in a Newton iteration
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
name	groupName	required	A name is required for any non-unique nodes
solutionChangeScalingFactor	real64	0.5	Damping factor for solution change targets
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.
targetRelativePressureChangeInTimeStep	real64	0.2	Target (relative) change in pressure in a time step (expected value between 0 and 1)
temperature	real64	required	Temperature
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
    gravityVector="{ 0.0, 0.0, -9.81 }">
    <ImmiscibleMultiphaseFlow
      name="FlowSolver"
      discretization="TPFA"
      targetRegions="{ Domain }"
      logLevel="3"
      writeLinearSystem="2"
      temperature="300">
      <NonlinearSolverParameters
        newtonTol="1.0e-6"
        newtonMaxIter="8"/>
      <LinearSolverParameters
        directParallel="0"/>
    </ImmiscibleMultiphaseFlow>
  </Solvers>