This model implements a $M$-\/phase flow of a fluid mixture composed of $N$ chemical species. The phases are denoted by lower index $\alpha\in\{1, \dots, M \}$. All fluid phases are mixtures of $N \geq M -1$ chemical species which are denoted by the upper index $\kappa\in\{1, \dots, N \}$.
By default, the standard multi-\/phase Darcy approach is used to determine the velocity, i.\-e. \[\mathbf{v}_\alpha=-\frac{k_{r\alpha}}{\mu_\alpha}\mathbf{K}\left(\text{grad}\, p_\alpha-\varrho_{\alpha}\mathbf{g}\right)\;, \] although the actual approach which is used can be specified via the {\ttfamily Velocity\-Module} property. For example, the velocity model can by changed to the Forchheimer approach by
The core of the model is the conservation mass of each component by means of the equation \[\sum_\alpha\frac{\partial\;\phi c_\alpha^\kappa S_\alpha}{\partial t}-\sum_\alpha\text{div}\left\{ c_\alpha^\kappa\mathbf{v}_\alpha\right\}- q^\kappa=0\;. \]
For the missing $M$ model assumptions, the model uses non-\/linear complementarity functions. These are based on the observation that if a fluid phase is not present, the sum of the mole fractions of this fluid phase is smaller than $1$, i.\-e. \[\forall\alpha: S_\alpha=0\implies\sum_\kappa x_\alpha^\kappa\leq1\]
Also, if a fluid phase may be present at a given spatial location its saturation must be non-\/negative\-: \[\forall\alpha: \sum_\kappa x_\alpha^\kappa=1\implies S_\alpha\geq0\]
Since at any given spatial location, a phase is always either present or not present, one of the strict equalities on the right hand side is always true, i.\-e. \[\forall\alpha: S_\alpha\left(\sum_\kappa x_\alpha^\kappa-1\right)=0\] always holds.
These three equations constitute a non-\/linear complementarity problem, which can be solved using so-\/called non-\/linear complementarity functions $\Phi(a, b)$. Such functions have the property \[\Phi(a,b)=0\iff a \geq0\land b \geq0\land a \cdot b =0\]
Several non-\/linear complementarity functions have been suggested, e.\-g. the Fischer-\/\-Burmeister function \[\Phi(a,b)= a + b -\sqrt{a^2+ b^2}\;. \] This model uses \[\Phi(a,b)=\min\{a, b \}\;, \] because of its piecewise linearity.
These equations are then discretized using a fully-\/implicit vertex centered finite volume scheme (often known as 'box'-\/scheme) for spatial discretization and the implicit Euler method as temporal discretization.
The model assumes local thermodynamic equilibrium and uses the following primary variables\-:
\begin{itemize}
\item The pressure of the first phase $p_1$
\item The component fugacities $f^1, \dots, f^{N}$
\item The saturations of the first $M-1$ phases $S_1, \dots, S_{M-1}$
\item Temperature $T$ if the energy equation is enabled
