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:

XML Element: SinglePhaseFVM

Name	Type	Default	Description
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]
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. SourceFluxes application if isThermal is enabled : - negative value (injection): the mass balance equation is modified to considered the additional source term, - positive value (production): both the mass balance and the energy balance equations are modified to considered the additional source term. For the energy balance equation, the mass flux is multipied by the enthalpy in the cell from which the fluid is being produced.
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
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
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.
temperature	real64	0	Temperature
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

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"
      targetRegions="{ mainRegion }">
      <NonlinearSolverParameters
        newtonTol="1.0e-6"
        newtonMaxIter="8"/>
      <LinearSolverParameters
        solverType="gmres"
	preconditionerType="amg" 
        krylovTol="1.0e-10"/>
    </SinglePhaseFVM>
  </Solvers>

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