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 GEOSX 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 GEOSX 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 usage 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 Finite Volume Discretization.

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:

Name Type Default Description
allowLocalCompDensityChopping integer 1 Flag indicating whether local (cell-wise) chopping of negative compositions is allowed
capPressureNames string_array {} Name of the capillary pressure constitutive model to use
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 string required Name of discretization object to use for this solver.
fluidNames string_array required Names of fluid constitutive models for each region.
initialDt real64 1e+99 Initial time-step value required by the solver to the event manager.
inputFluxEstimate real64 1 Initial estimate of the input flux used only for residual scaling. This should be essentially equivalent to the input flux * dt.
logLevel integer 0 Log level
maxCompFractionChange real64 1 Maximum (absolute) change in a component fraction between two Newton iterations
meanPermCoeff real64 1 Coefficient to move between harmonic mean (1.0) and arithmetic mean (0.0) for the calculation of permeability between elements.
name string required A name is required for any non-unique nodes
relPermNames string_array required Name of the relative permeability constitutive model to use
solidNames string_array required Names of solid constitutive models for each region.
targetRegions string_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.
temperature real64 required Temperature
useMass integer 0 Use mass formulation instead of molar
LinearSolverParameters node unique Element: LinearSolverParameters
NonlinearSolverParameters node unique Element: NonlinearSolverParameters

Example

  <Solvers>
    <CompositionalMultiphaseFlow
      name="compflow"
      logLevel="1"
      discretization="fluidTPFA"
      targetRegions="{ Channel }"
      fluidNames="{ fluid1 }"
      solidNames="{ rock }"
      relPermNames="{ relperm }"
      temperature="300"
      useMass="0">
      <NonlinearSolverParameters
        newtonTol="1.0e-10"
        newtonMaxIter="15"
        maxTimeStepCuts="2"
        lineSearchMaxCuts="2"/>
      <LinearSolverParameters
        solverType="direct"
        directParallel="0"
        logLevel="0"/>
    </CompositionalMultiphaseFlow>
  </Solvers>

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