Finite Volume discretization of a two-phase flow pressure equation of the sequential IMPES model. More...
#include <dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressure.hh>
Finite Volume discretization of a two-phase flow pressure equation of the sequential IMPES model.
This model implements two-phase flow of two immiscible fluids \(\alpha \in \{ w, n \}\) using a standard multi-phase Darcy approach as the equation for the conservation of momentum, i.e.
\[ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} {\textbf g} \right). \]
By inserting this into the equation for the conservation of the phase mass, one gets
\[ \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t} - \text{div} \left\{ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right) \right\} = q_\alpha \;. \]
In the incompressible case the phase densities \( \varrho_\alpha \) can be eliminated from the equations. The two equations are then added up, and because \(S_w + S_n = 1\), the first term cancels out. This leads to the so-called pressure equation. For the wetting ( \(w\)) phase pressure as primary variable this yields
\[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_w + f_n \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q, \]
for the non-wetting ( \( n \)) phase pressure as primary variable it yields
\[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_n - f_w \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q, \]
and for the global pressure as primary variable it leads to
\[ - \text{div}\, \left[\lambda \boldsymbol K \left(\textbf{grad}\, p_{global} - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q. \]
Here, \( p_\alpha \) is a phase pressure, \( p_ {global} \) the global pressure of a classical fractional flow formulation (see e.g. Binning and Celia (1999) [11]), \( p_c = p_n - p_w \) is the capillary pressure, \( \boldsymbol K \) the absolute permeability tensor, \( \lambda = \lambda_w + \lambda_n \) the total mobility depending on the saturation ( \( \lambda_\alpha = k_{r_\alpha} / \mu_\alpha \)), \( f_\alpha = \lambda_\alpha / \lambda \) the fractional flow function of a phase, \( \varrho_\alpha \) a phase density, \( {\textbf g} \) the gravitational acceleration vector and \( q = q_w + q_n \) the total source term. Depending on the primary variable chosen, one of the pressure equations above is solved sequentially together with the conservation of the phase mass for one phase.
For all cases, \( p = p_D \) on \( \Gamma_{Dirichlet} \), and \( \boldsymbol v_{total} \cdot \boldsymbol n = q_N \) on \( \Gamma_{Neumann} \).
The slightly compressible case is only implemented for phase pressures! In this case for a wetting \((w) \) phase pressure as primary variable the pressure equation is formulated as
\[ \phi \left( \varrho_w \frac{\partial S_w}{\partial t} + \varrho_n \frac{\partial S_n}{\partial t}\right) - \text{div}\, \left[\lambda \boldsymbol{K} \left(\textbf{grad}\, p_w + f_n \, \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q, \]
and for a non-wetting ( \( n \)) phase pressure as primary variable as
\[ \phi \left( \varrho_w \frac{\partial S_w}{\partial t} + \varrho_n \frac{\partial S_n}{\partial t}\right) - \text{div}\, \left[\lambda \boldsymbol{K} \left(\textbf{grad}\, p_n - f_w \textbf{grad}\, p_c - \sum f_\alpha \varrho_\alpha {\textbf g}\right)\right] = q. \]
In this slightly compressible case the following definitions are valid: \( \lambda = \varrho_w \lambda_w + \varrho_n \lambda_n \), \( f_\alpha = (\varrho_\alpha \lambda_\alpha) / \lambda \) This model assumes that temporal changes in density are very small and thus terms of temporal derivatives are negligible in the pressure equation. Depending on the formulation the terms including time derivatives of saturations are simplified by inserting \( S_w + S_n = 1 \).
In the IMPES models the default setting is:
TypeTag | The Type Tag |
Public Member Functions | |
void | getSource (EntryType &entry, const Element &element, const CellData &cellData, const bool first) |
Function which calculates the source entry. More... | |
void | getStorage (EntryType &entry, const Element &element, const CellData &cellData, const bool first) |
Function which calculates the storage entry. More... | |
void | getFlux (EntryType &entry, const Intersection &intersection, const CellData &cellData, const bool first) |
Function which calculates the flux entry. More... | |
void | getFluxOnBoundary (EntryType &entry, const Intersection &intersection, const CellData &cellData, const bool first) |
Function which calculates the boundary flux entry. More... | |
void | updateMaterialLaws () |
Updates and stores constitutive relations. More... | |
void | initialize (bool solveTwice=true) |
Initializes the pressure model. More... | |
void | update () |
Pressure update. More... | |
void | updateVelocity () |
Velocity update. More... | |
void | storePressureSolution () |
Globally stores the pressure solution. More... | |
void | storePressureSolution (const Element &element) |
Stores the pressure solution of a cell. More... | |
template<class MultiWriter > | |
void | addOutputVtkFields (MultiWriter &writer) |
Adds pressure output to the output file. More... | |
FVPressure2P (Problem &problem) | |
Constructs a FVPressure2P object. More... | |
const Scalar | pressure (const int eIdxGlobal) const |
Public access function for the primary pressure variable. More... | |
const Matrix & | globalMatrix () |
Returns the global matrix of the last pressure solution step. More... | |
const RHSVector & | rightHandSide () |
Returns the right hand side vector of the last pressure solution step. More... | |
void | initialize () |
Initialize pressure model. More... | |
void | calculateVelocity () |
void | serializeEntity (std::ostream &outstream, const Element &element) |
Function for serialization of the pressure field. More... | |
void | deserializeEntity (std::istream &instream, const Element &element) |
Function for deserialization of the pressure field. More... | |
void | setFixPressureAtIndex (Scalar pressure, int eIdxGlobal) |
Set a pressure to be fixed at a certain cell. More... | |
void | unsetFixPressureAtIndex (int eIdxGlobal) |
Reset the fixed pressure state. More... | |
void | resetFixPressureAtIndex () |
Protected Types | |
enum | { rhs = 1 , matrix = 0 } |
Indices of matrix and rhs entries. More... | |
enum | { pressEqIdx = Indices::pressureEqIdx } |
using | EntryType = Dune::FieldVector< Scalar, 2 > |
Protected Member Functions | |
void | initializeMatrix () |
Initialize the global matrix of the system of equations to solve. More... | |
void | initializeMatrixRowSize () |
Initialize the global matrix of the system of equations to solve. More... | |
void | initializeMatrixIndices () |
Initialize the global matrix of the system of equations to solve. More... | |
void | assemble (bool first) |
Function which assembles the system of equations to be solved. More... | |
void | solve () |
Solves the global system of equations to get the spatial distribution of the pressure. More... | |
PressureSolution & | pressure () |
Returns the vector containing the pressure solution. More... | |
const PressureSolution & | pressure () const |
Returns the vector containing the pressure solution. More... | |
void | initializePressure () |
Initialization of the pressure solution vector: Initialization with meaningful values may. More... | |
Protected Attributes | |
Matrix | A_ |
Global stiffnes matrix (sparse matrix which is build by the initializeMatrix() function) More... | |
RHSVector | f_ |
Right hand side vector. More... | |
|
protectedinherited |
Type of the vector of entries
Contains the return values of the get*()-functions (matrix or right-hand side entry).
|
protectedinherited |
Indices of matrix and rhs entries.
During the assembling of the global system of equations get-functions are called (getSource(), getFlux(), etc.), which return global matrix or right hand side entries in a vector. These can be accessed using following indices:
Enumerator | |
---|---|
rhs | index for the right hand side entry |
matrix | index for the global matrix entry |
|
inline |
Constructs a FVPressure2P object.
problem | A problem class object |
|
inline |
Adds pressure output to the output file.
Adds the phase pressures or a global pressure (depending on the formulation) as well as the capillary pressure to the output. In the compressible case additionally density and viscosity are added.
MultiWriter | Class defining the output writer |
writer | The output writer (usually a VTKMultiWriter object) |
|
protectedinherited |
Function which assembles the system of equations to be solved.
This function assembles the Matrix and the right hand side (RHS) vector to solve for a pressure field with a Finite-Volume (FV) discretization. Implementations of this base class have to provide the methods getSource()
, getStorage()
, getFlux()
and getFluxOnBoundary()
if the assemble() method is called!
first | Indicates if function is called at the initialization step or during the simulation (If first is true , no pressure field of previous iterations is required) |
|
inlineinherited |
|
inlineinherited |
Function for deserialization of the pressure field.
Function needed for restart option. Reads the pressure of a grid element from a restart file.
instream | Stream from the restart file. |
element | Grid element |
void Dumux::FVPressure2P< TypeTag >::getFlux | ( | EntryType & | entry, |
const Intersection & | intersection, | ||
const CellData & | cellData, | ||
const bool | first | ||
) |
Function which calculates the flux entry.
Function computes the inter-cell flux term and writes it to the corresponding entry of the entry vector
entry | Vector containing return values of the function |
intersection | Intersection of two grid elements |
cellData | Object containing all model relevant cell data |
first | Indicates if function is called in the initialization step or during the simulation |
void Dumux::FVPressure2P< TypeTag >::getFluxOnBoundary | ( | EntryType & | entry, |
const Intersection & | intersection, | ||
const CellData & | cellData, | ||
const bool | first | ||
) |
Function which calculates the boundary flux entry.
Function which calculates the flux entry at a boundary.
Function computes the boundary-flux term and writes it to the corresponding entry of the entry vector
entry | Vector containing return values of the function |
intersection | Intersection of two grid elements |
cellData | Object containing all model relevant cell data |
first | Indicates if function is called in the initialization step or during the simulation |
Dirichlet boundary condition for pressure depends on the formulation ( \(p_w\) (default), \(p_n\), \(p_{global}\)), Neumann boundary condition are the phase mass fluxes ( \(q_w\) and \(q_n\), [ \(\text{kg}/(\text{m}^2 \text{s}\)])
void Dumux::FVPressure2P< TypeTag >::getSource | ( | EntryType & | entry, |
const Element & | element, | ||
const CellData & | cellData, | ||
const bool | first | ||
) |
Function which calculates the source entry.
Function computes the source term and writes it to the corresponding entry of the entry vector
entry | Vector containing return values of the function |
element | Grid element |
cellData | Object containing all model relevant cell data |
first | Indicates if function is called in the initialization step or during the simulation |
Source of each fluid phase has to be added as mass flux ( \(\text{kg}/(\text{m}^3 \text{s}\)).
void Dumux::FVPressure2P< TypeTag >::getStorage | ( | EntryType & | entry, |
const Element & | element, | ||
const CellData & | cellData, | ||
const bool | first | ||
) |
Function which calculates the storage entry.
Function computes the storage term and writes it to the corresponding entry of the entry vector
entry | Vector containing return values of the function |
element | Grid element |
cellData | Object containing all model relevant cell data |
first | Indicates if function is called in the initialization step or during the simulation |
If compressibility is enabled this functions calculates the term
\[ \phi \sum_\alpha \varrho_\alpha \frac{\partial S_\alpha}{\partial t} V \]
In the incompressible case an volume correction term is calculated which corrects for non-physical saturation overshoots/undershoots. These can occur if the estimated time step for the explicit transport was too large. Correction by an artificial source term allows to correct this errors due to wrong time-stepping without losing mass conservation. The error term looks as follows:
\[ q_{error} = \begin{cases} S < 0 & a_{error} \frac{S}{\Delta t} V \\ S > 1 & a_{error} \frac{(S - 1)}{\Delta t} V \\ 0 \le S \le 1 & 0 \end{cases} \]
where \(a_{error}\) is a weighting factor (default: \(a_{error} = 0.5\))
|
inlineinherited |
Returns the global matrix of the last pressure solution step.
|
inlineinherited |
Initialize pressure model.
Function initializes the sparse matrix to solve the global system of equations and sets/calculates the initial pressure
|
inline |
Initializes the pressure model.
Function initializes the sparse matrix to solve the global system of equations and sets/calculates the initial pressure
solveTwice | indicates if more than one iteration is allowed to get an initial pressure solution |
|
protectedinherited |
Initialize the global matrix of the system of equations to solve.
|
protectedinherited |
Initialize the global matrix of the system of equations to solve.
|
protectedinherited |
Initialize the global matrix of the system of equations to solve.
|
inlineprotectedinherited |
Initialization of the pressure solution vector: Initialization with meaningful values may.
|
inlineprotectedinherited |
Returns the vector containing the pressure solution.
|
inlineprotectedinherited |
Returns the vector containing the pressure solution.
|
inlineinherited |
Public access function for the primary pressure variable.
Function returns the cell pressure value at index eIdxGlobal
eIdxGlobal | Global index of a grid cell |
|
inlineinherited |
|
inlineinherited |
Returns the right hand side vector of the last pressure solution step.
|
inlineinherited |
Function for serialization of the pressure field.
Function needed for restart option. Writes the pressure of a grid element to a restart file.
outstream | Stream into the restart file. |
element | Grid element |
|
inlineinherited |
Set a pressure to be fixed at a certain cell.
Allows to fix a pressure somewhere (at one certain cell) in the domain. This can be necessary e.g. if only Neumann boundary conditions are defined. The pressure is fixed until the unsetFixPressureAtIndex()
function is called
pressure | Pressure value at eIdxGlobal |
eIdxGlobal | Global index of a grid cell |
|
protectedinherited |
Solves the global system of equations to get the spatial distribution of the pressure.
|
inline |
Globally stores the pressure solution.
|
inline |
Stores the pressure solution of a cell.
Calculates secondary pressure variables and stores pressures.
element | Grid element |
|
inlineinherited |
Reset the fixed pressure state.
No pressure is fixed inside the domain until setFixPressureAtIndex()
function is called again.
eIdxGlobal | Global index of a grid cell |
|
inline |
Pressure update.
Function reassembles the system of equations and solves for a new pressure solution.
void Dumux::FVPressure2P< TypeTag >::updateMaterialLaws |
Updates and stores constitutive relations.
Updates constitutive relations and stores them in the variable class.
Stores mobility, fractional flow function and capillary pressure for all grid cells. In the compressible case additionally the densities and viscosities are stored.
TypeTag | The problem TypeTag |
|
inline |
Velocity update.
Resets the velocities in the cellData.
|
protectedinherited |
Global stiffnes matrix (sparse matrix which is build by the initializeMatrix()
function)
|
protectedinherited |
Right hand side vector.