Singlephase Flow Solver

Introduction

Here, we describe the single-phase flow solver. The role of this solver is to implement the fully implicit finite-volume discretization (mainly, accumulation and source terms, boundary conditions) of the equations governing compressible single-phase flow in porous media. This solver can be combined with the SinglePhaseWell class which handles the discrete multi-segment well model and provides source/sink terms for the fluid flow solver.

Theory

Governing Equations

This is a cell-centered Finite Volume solver for compressible single-phase flow in porous media. Fluid pressure as the primary solution variable. Darcy’s law is used to calculate fluid velocity from pressure gradient. The solver currently only supports Dirichlet-type boundary conditions (BC) applied on cells or faces and Neumann no-flow type BC.

The following mass balance equation is solved in the domain:

\frac{\partial}{\partial t}(\phi\rho) + \boldsymbol{\nabla} \cdot (\rho\boldsymbol{u}) + q = 0,

where

\boldsymbol{u} = -\frac{1}{\mu}\boldsymbol{k}(\nabla p - \rho \boldsymbol{g})

and \phi is porosity, \rho is fluid density, \mu is fluid viscosity, \boldsymbol{k} is the permeability tensor, \boldsymbol{g} is the gravity vector, and q is the source function (currently not supported). The details on the computation of the density and the viscosity are given in Compressible single phase fluid model.

When the entire pore space is filled by a single phase, we can substitute the Darcy’s law into the mass balance equation to obtain the single phase flow equation

\frac{\partial}{\partial t}(\phi\rho) - \boldsymbol{\nabla} \cdot \frac{\rho \boldsymbol{k}}{\mu} (\nabla p - \gamma \nabla z) + q = 0,

with \gamma \nabla z= \rho \boldsymbol{g}.

Discretization

Space Discretization

Let \Omega \subset \mathbb{R}^n, \, n =1,2,3 be an open set defining the computational domain. We consider \Omega meshed by element such that \Omega = \cup_{i}V_i and integrate the single phase flow equation, described above, over each element V_i:

\int_{V_i} \frac{\partial}{\partial t}(\phi\rho) dV - \int_{V_i} \boldsymbol{\nabla} \cdot \frac{\rho \boldsymbol{k}}{\mu} (\nabla p - \gamma \nabla z) dV + \int_{V_i} q dV  = 0.

Applying the divergence theorem to the second term leads to

\int_{V_i} \frac{\partial}{\partial t}(\phi\rho)_i - \oint_{S_i} \left(\frac{\rho \boldsymbol{k}}{\mu}(\nabla p -\gamma \nabla z)\right) \cdot \boldsymbol{n} dS + \int_{V_i} q dV  = 0.

where S_i represents the surface area of the element V_i and \boldsymbol{n} is a outward unit vector normal to the surface.

For the flux term, the (static) transmissibility is currently computed with a Two-Point Flux Approximation (TPFA) as described in Finite Volume Discretization.

The pressure-dependent mobility \lambda = \frac{\rho}{\mu} at the interface is approximated using a first-order upwinding on the sign of the potential difference.

Time Discretization

Let t_0 < t_1 < \cdots < t_N=T be a grid discretization of the time interval [t_0,T], \, t_0, T \in \mathbb{R}^+. We use the backward Euler (fully implicit) method to integrate the single phase flow equation between two grid points t_n and t_{n+1}, \, n< N to obtain the residual equation:

\int_{V_i} \frac{(\phi\rho)_i^{n+1} - (\phi\rho)_i^n}{\Delta t} - \oint_{S_i} \left(\frac{\rho \boldsymbol{k}}{\mu}(\nabla p -\gamma \nabla z)\right)^{n+1} \cdot \boldsymbol{n} dS + \int_{V_i} q^{n+1} dV = 0

where \Delta t = t_{n+1}-t_n is the time-step. The expression of this residual equation and its derivative are used to form a linear system, which is solved via the solver package.

Parameters

The solver is enabled by adding a <SinglePhaseFVM> node in the Solvers section. Like any solver, time stepping is driven by events, see Event Management.

The following attributes are supported:

Name Type Default Description
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
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
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.
LinearSolverParameters node unique Element: LinearSolverParameters
NonlinearSolverParameters node unique Element: NonlinearSolverParameters

In particular:

  • discretization must point to a Finite Volume flux approximation scheme defined in the Numerical Methods section of the input file (see Finite Volume Discretization)
  • fluidName must point to a single phase fluid model defined in the Constitutive section of the input file (see Constitutive Models)
  • solidName must point to a solid mechanics model defined in the Constitutive section of the input file (see Constitutive Models)
  • targetRegions is used to specify the regions on which the solver is applied

Primary solution field label is pressure. Initial conditions must be prescribed on this field in every region, and boundary conditions must be prescribed on this field on cell or face sets of interest.

Example

  <Solvers>
    <SinglePhaseFVM
      name="SinglePhaseFlow"
      logLevel="1"
      discretization="singlePhaseTPFA"
      fluidNames="{ water }"
      solidNames="{ rock }"
      targetRegions="{ mainRegion }">
      <NonlinearSolverParameters
        newtonTol="1.0e-6"
        newtonMaxIter="8"/>
      <LinearSolverParameters
        solverType="gmres"
        krylovTol="1.0e-10"/>
    </SinglePhaseFVM>
  </Solvers>

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