Compositional multiphase fluid model

Overview

This model represents a full composition description of a multiphase multicomponent fluid. Phase behavior is modeled by an equation of state (EOS) and partitioning of components into phases is computed based on instantaneous chemical equilibrium via a two-phase flash. Each component (species) is characterized by molar weight and critical properties that serve as input parameters for the EOS. See Petrowiki for more information.

In this model the fluid is described by $N_c$ components with $z_c$ being the total mole fraction of component $c$ . The fluid can partition into a liquid phase, denoted $\ell$ , and a vapor phase denoted by $v$ . Therefore, by taking into account the molar phase component fractions, (which is the fraction of the molar mass of phase $p$ represented by component $c$ ), the following partition matrix establishes the component distribution within the two phases:

$\begin{bmatrix} x_{1} & x_{2} & x_{3} & \cdots & x_{N_c} \\ y_{1} & y_{2} & y_{3} & \cdots & y_{N_c} \\ \end{bmatrix}$

where $x_c$ is the mole fraction of component $c$ in the liquid phase and $y_c$ is the mole fraction of component $c$ in the vapor phase.

The fluid properties are updated through the following steps:

1) The phase fractions ( $\nu_p$ ) and phase component fractions ( $x_c$ and $y_c$ ) are computed as a function of pressure ( $p$ ), temperature ( $T$ ) and total component fractions ( $z_c$ ).

2) The phase densities ( $\rho_p$ ) and phase viscosities ( $\mu_p$ ) are computed as a function of pressure, temperature and the updated phase component fractions.

After calculating the phase fractions, phase component fractions, phase densities, phase viscosities, and their derivatives with respect to pressure, temperature, and component fractions, the Compositional Multiphase Flow Solver then moves on to assembling the accumulation and flux terms.

Step 1: Computation of the phase fractions and phase component fractions (flash)

Stability test

The first step is to determine if the provided mixture with total molar fractions $z_c$ is stable as a single phase at the current pressure $p$ and temperature $T$ . However, this can only be confirmed through stability testing.

The stability of a mixture is traditionally assessed using the Tangent Plane Distance (TPD) criterion developed by Michelsen (1982a). This criterion states that a phase with composition $z$ is stable at a specified pressure $p$ and temperature $T$ if and only if

$g(y) = \sum_{i=1}^{N_c}y_i\left(\ln y_i + \ln \phi_i(y) - \ln z_i - \ln \phi_i(z) \right) \geq 0$

for any permissible trial composition $y$ , where $\phi_i$ denotes the fugacity coefficient of component $i$ .

To determine stability of the mixture this testing in initiated from a basic starting point, based on Wilson K-values, to get both a lighter and a heavier trial mixture. The two trial mixtures are calculated as $y_i = z_i/K_i$ and $y_i = z_iK_i$ where $K_i$ are defined by

$K_i = \frac{P_{ci}}{p}\exp\left( 5.37( 1 + \omega_i ) \left(1-\frac{T_{ci}}{T}\right)\right)$

where $P_{ci}$ and $T_{ci}$ are respectively, the critical pressure and temperature of component $i$ and $\omega_i$ is the accentric factor of component $i$ .

The stability problem is solved by observing that a necessary condition is that $g(y)$ must be non-negative at all its stationary points. The stationarity criterion can be expressed as

$\ln y_i + \ln \phi_i(y) - h_i = k \hspace{1cm} i=1,2,3,\ldots,N_c$

where $h_i = \ln z_i + \ln \phi_i(z)$ is a constant parameter dependent on the feed composition $z$ and $k$ is an undetermined constant. This constant can be further incorporated into the equation by defining the unnormalized trial phase moles $Y_i$ as

$Y_i = \exp(-k)y_i$

which reduces the stationarity criterion to

$\ln Y_i + \ln \phi_i(y) - h_i = 0$

with the mole fractions $y_i$ related to the trial phase moles $Y_i$ by

$y_i = Y_i/\sum_jY_j$

With the two starting mixtures, the stationarity condition is solved using successive substitution to determine the stationary points. If both initial states converge to a solution which has $g(y)\geq 0$ then the mixture is deemed to be stable, otherwise it is deemed unstable.

Phase labeling

Once it is confirmed that the fluid with composition $z$ is stable as a single phase at the current pressure and temperature, it must be labeled as either ‘liquid’ or ‘vapor’. This is necessary only to apply the correct relative permeability function for calculating the phase’s flow properties. The properties of the fluid (density, viscosity) are unchanged by the assignment of the label.

Determining the mixture’s true critical point is the most rigorous method for labeling. It is however expensive and may not always be necessary. As such, a simple correlation for pseudo-critical temperature is used and this is expected to be sufficiently accurate for correct phase labeling, except under some specific conditions.

The Li-correlation is a weighted average of the component critical temperatures and is used to determine the label applied to the mixture. The Li pseudo-critical temperature is calcaulated as

$T_{cp} = \frac{\sum_{i=1}^{N_c}T_{ci}V_{ci}z_{i}}{\sum_{i=1}^{N_c}V_{ci}z_{i}}$

where $V_{ci}$ and $T_{ci}$ are respectively the critical volume and temperature of component $i$ . This is compared to the current temperature $T$ such that if $T_{cp}<T$ then the mixture is labeled as vapor and as liquid otherwise.

Negative two-phase flash

When a cell is identified as having an unstable mixture, it is necessary to determine the amounts in the liquid and vapor phases through phase splitting. This phase split is calculated by ensuring that the two phases are in thermodynamic equilibrium. For a system to be in thermodynamic equilibrium, the fugacities of each component in both the liquid and vapor phases must be equal:

$\phi_{iL} = \phi_{iV} \hspace{1cm} i=1,2,3,\ldots,N_c$

where $\phi_{iL}$ is the fugacity of component $i$ in the liquid phase and $\phi_{iL}$ is the fugacity of component $i$ in the vapor phase.

Fugacities are functions of temperature, pressure, and composition:

$\phi_{iL} = \phi_{iL}(T, p, x_i) \hspace{1cm} i=1,2,3,\ldots,N_c$

and

$\phi_{iV} = \phi_{iV}(T, p, y_i) \hspace{1cm} i=1,2,3,\ldots,N_c$

and are calculated directly from an equation of state.

Equilibrium constants, also known as K-values, are defined for each component as:

$K_i = \frac{y_i}{x_i}$

where $x_i$ is the mole fraction of component $i$ in the liquid phase and $y_i$ is the mole fraction of component $i$ in the vapor phase. If we denote $V$ as the mole fraction of the vapor phase, the material balance indicates that the mole fractions of each component in the liquid and vapor phases are given by:

$x_i = \frac{z_i}{1 + (K_i - 1)V}$

and

$y_i = \frac{K_i z_i}{1 + (K_i - 1)V}$

The value of $V$ corresponding to a given set of K-values is determined by solving the so called Rachford and-Rice equation:

$F(V) = \sum_{i=1}^{N_c} \left(x_i - y_i\right) = \sum_{i=1}^{N_c} \frac{z_i(1 - K_i)}{1 + (K_i - 1)V} = 0$

The flash calculation process is as follows:

Once the mixture is confirmed to be stable, an initial set of K-values is chosen, typically using Wilson’s formula.
Given $z_i$ and $K_i$ , the Rachford-Rice equation is solved to determine the molar fraction of vapor, $V$ . This is initially solved using successive substitution, followed by Newton iterations once the residual is sufficiently reduced.
After $V$ is calculated, the corresponding liquid and vapor mole fractions, $x_i$ and $y_i$ , are computed.
These phase compositions are then used to calculate the component fugacities $\phi_{iL}$ and $\phi_{iV}$ in the liquid and vapor phases using the equation of state.
Convergence is reached when the fugacities are equal for all components. The convergence criterion is defined as:

$\sum_{i=1}^{N_c} \left( \phi_{iL} - \phi_{iV} \right)^2 < \varepsilon$

where $\varepsilon$ is the convergence tolerance.
If convergence is not achieved, successive substitution is used to update the set of K-values for the next iteration. The new K-values at iteration $t+1$ are given by:

$K_i^{(t+1)} = K_i^{(t)} \frac{\phi_{iL}}{\phi_{iV}}$

Parameters

The model represented by <CompositionalMultiphaseFluid> node in the input. Under the hood this is a wrapper around PVTPackage library, which is included as a submodule. In order to use the model, GEOS must be built with -DENABLE_PVTPACKAGE=ON (default).

The following attributes are supported:

XML Element: CompositionalMultiphaseFluid

Name	Type	Default	Description
checkPVTTablesRanges	integer	1	Enable (1) or disable (0) an error when the input pressure or temperature of the PVT tables is out of range.
componentAcentricFactor	real64_array	required	Component acentric factors
componentBinaryCoeff	real64_array2d	{{0}}	Table of binary interaction coefficients
componentCriticalPressure	real64_array	required	Component critical pressures
componentCriticalTemperature	real64_array	required	Component critical temperatures
componentMolarWeight	real64_array	required	Component molar weights
componentNames	string_array	required	List of component names
componentVolumeShift	real64_array	{0}	Component volume shifts
constantPhaseViscosity	real64_array	{0}	Viscosity for each phase
equationsOfState	string_array	required	List of equation of state types for each phase
name	groupName	required	A name is required for any non-unique nodes
phaseNames	groupNameRef_array	required	List of fluid phases

Supported phase names are:

Value	Comment
oil	Oil phase
gas	Gas phase
water	Water phase

Supported Equation of State types:

Value	Comment
PR	Peng-Robinson EOS
SRK	Soave-Redlich-Kwong EOS

Example

<Constitutive>
  <CompositionalMultiphaseFluid name="fluid1"
                                phaseNames="{ oil, gas }"
                                equationsOfState="{ PR, PR }"
                                componentNames="{ N2, C10, C20, H2O }"
                                componentCriticalPressure="{ 34e5, 25.3e5, 14.6e5, 220.5e5 }"
                                componentCriticalTemperature="{ 126.2, 622.0, 782.0, 647.0 }"
                                componentAcentricFactor="{ 0.04, 0.443, 0.816, 0.344 }"
                                componentMolarWeight="{ 28e-3, 134e-3, 275e-3, 18e-3 }"
                                componentVolumeShift="{ 0, 0, 0, 0 }"
                                componentBinaryCoeff="{ { 0, 0, 0, 0 },
                                                      { 0, 0, 0, 0 },
                                                      { 0, 0, 0, 0 },
                                                      { 0, 0, 0, 0 } }"/>
</Constitutive>

References

M. L. Michelsen, The Isothermal Flash Problem. Part I. Stability., Fluid Phase Equilibria, vol. 9.1, pp. 1-19, 1982a.
M. L. Michelsen, The Isothermal Flash Problem. Part II. Phase-Split Calculation., Fluid Phase Equilibria, vol. 9.1, pp. 21-40, 1982b.