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.
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 :ref:`CompressibleSinglePhaseFluid`.
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}.
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 :ref:`FiniteVolume`.
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.
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.
The solver is enabled by adding a <SinglePhaseFVM> node in the Solvers section.
Like any solver, time stepping is driven by events, see :ref:`EventManager`.
The following attributes are supported:
In particular:
discretizationmust point to a Finite Volume flux approximation scheme defined in the Numerical Methods section of the input file (see :ref:`FiniteVolume`)fluidNamemust point to a single phase fluid model defined in the Constitutive section of the input file (see :ref:`Constitutive`)solidNamemust point to a solid mechanics model defined in the Constitutive section of the input file (see :ref:`Constitutive`)targetRegionsis 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.
.. literalinclude:: ../../../../../inputFiles/singlePhaseFlow/3D_10x10x10_compressible_base.xml :language: xml :start-after: <!-- SPHINX_TUT_INT_HEX_SOLVERS --> :end-before: <!-- SPHINX_TUT_INT_HEX_SOLVERS_END -->
We refer the reader to :ref:`this page <TutorialSinglePhaseFlowWithInternalMesh>` for a complete tutorial illustrating the use of this solver.