The finite volume model for the solution of the compositional pressure equation.
More...
template<class TypeTag>
class Dumux::FV2dPressure2P2CAdaptive< TypeTag >
The finite volume model for the solution of the compositional pressure equation.
Provides a Finite Volume implementation for the pressure equation of a compressible system with two components. An IMPES-like method is used for the sequential solution of the problem. Diffusion is neglected, capillarity can be regarded. Isothermal conditions and local thermodynamic equilibrium are assumed. Gravity is included.
\[ c_{total}\frac{\partial p}{\partial t} + \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} \nabla \cdot \left( \sum_{\alpha} X^{\kappa}_{\alpha} \varrho_{\alpha} \bf{v}_{\alpha}\right) = \sum_{\kappa} \frac{\partial v_{total}}{\partial C^{\kappa}} q^{\kappa}, \]
where \(\bf{v}_{\alpha} = - \lambda_{\alpha} \bf{K} \left(\nabla p_{\alpha} + \rho_{\alpha} \bf{g} \right) \). \( c_{total} \) represents the total compressibility, for constant porosity this yields \( - \frac{\partial V_{total}}{\partial p_{\alpha}} \), \(p_{\alpha} \) denotes the phase pressure, \( \bf{K} \) the absolute permeability, \( \lambda_{\alpha} \) the phase mobility, \( \rho_{\alpha} \) the phase density and \( \bf{g} \) the gravity constant and \( C^{\kappa} \) the total Component concentration. See paper SPE 99619 or "Analysis of a Compositional Model for Fluid
Flow in Porous Media" by Chen, Qin and Ewing for derivation. The partial derivatives of the actual fluid volume \( v_{total} \) are gained by using a secant method.
This adaptive implementation uses its own initialization and assembling methods for matrix and right-hand-side vector, but solution is done by the base class FVPressure. Fluxes near hanging nodes can be calculated using an mpfa method.
- Template Parameters
-
|
| void | initializeMatrix () |
| | initializes the matrix to store the system of equations More...
|
| |
| void | assemble (bool first) |
| | function which assembles the system of equations to be solved More...
|
| |
| void | getMpfaFlux (const IntersectionIterator &, const CellData &) |
| | Compute flux over an irregular interface using a mpfa method. More...
|
| |
| int | computeTransmissibilities (const IntersectionIterator &, IntersectionIterator &, TransmissivityMatrix &, TransmissivityMatrix &, GlobalPosition &, int &) |
| | Computes the transmissibility coefficients for the MPFA-l method. More...
|
| |
| void | adaptPressure () |
| | Adapt primary variables vector after adapting the grid. More...
|
| |
| | FV2dPressure2P2CAdaptive (Problem &problem) |
| | Constructs a FV2dPressure2P2CAdaptive object. More...
|
| |
| void | getSource (EntryType &sourceEntry, const Element &elementI, const CellData &cellDataI, const bool first) |
| | Assembles the source term. More...
|
| |
| void | getStorage (EntryType &storageEntry, const Element &elementI, const CellData &cellDataI, const bool first) |
| | Assembles the storage term. More...
|
| |
| void | getFlux (EntryType &entries, const Intersection &intersection, const CellData &cellDataI, const bool first) |
| | Get flux at an interface between two cells. More...
|
| |
| void | getFluxOnBoundary (EntryType &entries, const Intersection &intersection, const CellData &cellDataI, const bool first) |
| | Get flux on Boundary. More...
|
| |
| void | updateMaterialLawsInElement (const Element &elementI, bool postTimeStep) |
| | Updates secondary variables of one cell. More...
|
| |
| void | initialMaterialLaws (bool compositional) |
| | initializes the fluid distribution and hereby the variables container More...
|
| |
| void | initialize (bool solveTwice=false) |
| | Initializes the simulation run. More...
|
| |
| void | initialize () |
| | Initialize pressure model. More...
|
| |
| void | update () |
| | Compositional pressure solution routine: update estimate for secants, assemble, solve. More...
|
| |
| void | updateMaterialLaws (bool postTimeStep=false) |
| | Updates secondary variables. More...
|
| |
| void | volumeDerivatives (const GlobalPosition &, const Element &ep) |
| | Partial derivatives of the volumes w.r.t. changes in total concentration and pressure. 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 | calculateVelocity () |
| |
| void | updateVelocity () |
| |
| 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 () |
| |
|
| Problem & | problem () |
| |
| const Problem & | problem () const |
| |
| int | transmissibilityAdapter_ (const IntersectionIterator &, const IntersectionIterator &, IntersectionIterator &, GlobalPosition &, TransmissivityMatrix &) |
| | An adapter to use the traditional 2p implementation for the second interaction region. 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 | 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...
|
| |
template<class TypeTag >
| int Dumux::FV2dPressure2P2CAdaptive< TypeTag >::computeTransmissibilities |
( |
const IntersectionIterator & |
isIt, |
|
|
IntersectionIterator & |
additionalIntersectionIt, |
|
|
TransmissivityMatrix & |
T, |
|
|
TransmissivityMatrix & |
additionalT, |
|
|
GlobalPosition & |
globalPos3, |
|
|
int & |
globalIdx3 |
|
) |
| |
Computes the transmissibility coefficients for the MPFA-l method.
Indices used in a interaction volume of the MPFA-l method
___________________________________________________
| nux: cell geometry | nx: face normal |
| | =integrationOuterNormal|
| | |
| nextIt________________| |
| | nu4 3 |
| | \;.´ | |
| | .´ <--|nu3 |
| elementnumber | .´ ^n1 | |
| 1...........x_____|_____|____________|
| `. |nu5 !`. _nu7 | |
| `.`'´ ! `./` <--| |
lastIt | `._ ! `. nu2| |
\ | nu6/`. !nu1^ ` . | |
\ | `!---|------`2 |
`| | |
|\ isIt__________________|->n12 |
| \ | |
|_face14_________________|________________________|
_
/` is a normal directing upwards to the right,
while a normal directing downwards to the right looks like \;
- Parameters
-
| isIt | Iterator to the current intersection |
| additionalIntersectionIt | Iterator to the additional interface included in the second interaction region |
| T | Transmissitivity matrix of the first unique interaction region |
| additionalT | Transmissitivity matrix of the second non-unique interaction region |
| globalPos3 | Unique interaction region: Position of the 3rd cell, the other neighbor on the hanging node |
| globalIdx3 | Unique interaction region: Index of the 3rd cell, the other neighbor on the hanging node |
1) get geometrical information of interaction triangle
2) search for face13, face23
3) Calculate omega, chi for matrices
4) Calculate A, B, C, D and solve for T
template<class TypeTag >
| void Dumux::FVPressure2P2C< TypeTag >::getFlux |
( |
EntryType & |
entries, |
|
|
const Intersection & |
intersection, |
|
|
const CellData & |
cellDataI, |
|
|
const bool |
first |
|
) |
| |
|
inherited |
Get flux at an interface between two cells.
for first == true, the flux is calculated in traditional fractional-flow forn as in FVPressure2P. for first == false, the flux thorugh \( \gamma \) is calculated via a volume balance formulation
\[ - A_{\gamma} \mathbf{n}^T_{\gamma} \mathbf{K} \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \mathbf{d}_{ij} \left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}^T \mathbf{d}_{ij} \right) \sum_{\kappa} X^{\kappa}_{\alpha} \frac{\partial v_{t}}{\partial C^{\kappa}} + V_i \frac{A_{\gamma}}{U_i} \mathbf{d}^T \mathbf{K} \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \mathbf{d}_{ij} \left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}^T \mathbf{d}_{ij} \right) \sum_{\kappa} X^{\kappa}_{\alpha} \frac{\frac{\partial v_{t,j}}{\partial C^{\kappa}_j}-\frac{\partial v_{t,i}}{\partial C^{\kappa}_i}}{\Delta x} \]
This includes a boundary integral and a volume integral, because \( \frac{\partial v_{t,i}}{\partial C^{\kappa}_i} \) is not constant. Here, \( \mathbf{d}_{ij} \) is the normalized vector connecting the cell centers, and \( \mathbf{n}_{\gamma} \) represents the normal of the face \( \gamma \).
- Parameters
-
| entries | The Matrix and RHS entries |
| intersection | Intersection between cell I and J |
| cellDataI | Data of cell I |
| first | Flag if pressure field is unknown |
template<class TypeTag >
| void Dumux::FVPressure2P2C< TypeTag >::getFluxOnBoundary |
( |
EntryType & |
entries, |
|
|
const Intersection & |
intersection, |
|
|
const CellData & |
cellDataI, |
|
|
const bool |
first |
|
) |
| |
|
inherited |
Get flux on Boundary.
for first == true, the flux is calculated in traditional fractional-flow forn as in FVPressure2P. for first == false, the flux thorugh \( \gamma \) is calculated via a volume balance formulation
\[ - A_{\gamma} \mathbf{n}^T_{\gamma} \mathbf{K} \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \mathbf{d}_{ij} \left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}^T \mathbf{d}_{ij} \right) \sum_{\kappa} \frac{\partial v_{t}}{\partial C^{\kappa}} X^{\kappa}_{\alpha} \;, \]
where we skip the volume integral assuming \( \frac{\partial v_{t,i}}{\partial C^{\kappa}_i} \) to be constant at the boundary. Here, \( \mathbf{d}_{ij} \) is the normalized vector connecting the cell centers, and \( \mathbf{n}_{\gamma} \) represents the normal of the face \( \gamma \).
If a Neumann BC is set, the given (mass-)flux is directly multiplied by the volume derivative and inserted.
- Parameters
-
| entries | The Matrix and RHS entries |
| intersection | Intersection between cell I and J |
| cellDataI | Data of cell I |
| first | Flag if pressure field is unknown |
Compute flux over an irregular interface using a mpfa method.
A mpfa l-method is applied to calculate fluxes near hanging nodes, using:
\[ - \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \left( \sum_k \tau_{2k} p^t_{\alpha,k} + \varrho_{\alpha} \sum_k \tau_{2k} \mathbf{g}^T \mathbf{x}_{k} \right) \sum_{\kappa} X^{\kappa}_{\alpha} \frac{\partial v_{t}}{\partial C^{\kappa}} + \frac{ V_i}{U_i} \sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \left( \sum_k \tau_{2k} p^t_{\alpha,k} + \varrho_{\alpha} \sum_k \tau_{2k} \mathbf{g}^T \mathbf{x}_{k} \right) \sum_{\kappa} X^{\kappa}_{\alpha} \frac{\frac{\partial v_{t,j}}{\partial C^{\kappa}_j} -\frac{\partial v_{t,i}}{\partial C^{\kappa}_i}}{\Delta x} \]
We provide two options: Calculating the flux expressed by twice the flux through the one unique interaction region on the hanging node if one halfedge is stored (eg on boundaries). Or using the second interaction region covering neighboring cells. The contribution in other cells than I or J make it necessary that the matrix and rhs entries are filled up within this function.
- Parameters
-
| intersectionIterator | Iterator of the intersection between cell I and J |
| cellDataI | Data of cell I |
get geometrical Info, transmissibility matrix
compute matrix entry: advective fluxes
template<class TypeTag >
| void Dumux::FVPressure2P2C< TypeTag >::getStorage |
( |
EntryType & |
storageEntry, |
|
|
const Element & |
elementI, |
|
|
const CellData & |
cellDataI, |
|
|
const bool |
first |
|
) |
| |
|
inherited |
Assembles the storage term.
for first == true, there is no storage contribution. for first == false, the storage term comprises the compressibility (due to a change in pressure from last timestep):
\[ V_i c_{t,i} \frac{p^t_i - p^{t-\Delta t}_i}{\Delta t} \]
and the damped error introduced by the incorrect transport of the last timestep:
\[ V_i \alpha_r \frac{v_{t} - \phi}{\Delta t} \]
. The latter is damped according to Fritz 2011.
- Parameters
-
| storageEntry | The Matrix and RHS entries |
| elementI | The element I |
| cellDataI | Data of cell I |
| first | Flag if pressure field is unknown |
template<class TypeTag >
| int Dumux::FV2dPressure2P2CAdaptive< TypeTag >::transmissibilityAdapter_ |
( |
const IntersectionIterator & |
isIt, |
|
|
const IntersectionIterator & |
face23_2p2cnaming, |
|
|
IntersectionIterator & |
additionalIntersectionIt, |
|
|
GlobalPosition & |
corner1234, |
|
|
TransmissivityMatrix & |
additionalT |
|
) |
| |
|
protected |
An adapter to use the traditional 2p implementation for the second interaction region.
The second interaction region consists of 4 cells, hence the traditional non-adaptive implementation to calculate the transmissibility coefficients is used. This, however, uses its own naming conventions and local Indices used in the work of I. Aavatsmark.
Indices used in a interaction volume of the MPFA-l method
| | |
| 4-----------3-----------3 |
| | --> nu43 | nu34 <-- | |
| | |nu41 1|--> n43 ||nu32 |
| | v ^ |0 ^ v| |
|____________4__0__|n14__|__n23_|_1__2____________|
| | 1 0 | 0 | |
| | ^ |1 nu23 ^ | |
| | |nu14 0|--> n12 | | |
| | -->nu12 | nu21<-- | |
| J = 1-----------1-----------2 = I |
| | |.
- Parameters
-
| isIt | Iterator to the current intersection |
| face23_2p2cnaming | Iterator to the second interface included in the unique interaction region = face23 |
| additionalIntersectionIt | Iterator to the intersection included in the second interaction region |
| corner1234 | Center point of non-unique interaction region. |
| additionalT | Transmissibility Matrix of non unique interaction region. |
1.9.3