From 847b1b82378ab185d69a8a56e57f1d35b30ee0e2 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Mon, 23 Jun 2025 22:40:50 -0500 Subject: [PATCH 01/13] Compositional fluid model doc --- .../docs/CompositionalMultiphaseFluid.rst | 255 ------------------ .../constitutive/docs/FluidModels.rst | 2 +- .../CompositionalMultiphaseFluid.rst | 12 + .../constitutive/docs/compositional/Flash.rst | 8 + .../compositional/NegativeTwoPhaseFlash.rst | 72 +++++ .../docs/compositional/Overview.rst | 40 +++ .../docs/compositional/References.rst | 14 + .../docs/compositional/StabilityTest.rst | 75 ++++++ 8 files changed, 222 insertions(+), 256 deletions(-) delete mode 100644 src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/Flash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/Overview.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/References.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/StabilityTest.rst diff --git a/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst deleted file mode 100644 index 66723ee555c..00000000000 --- a/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst +++ /dev/null @@ -1,255 +0,0 @@ -.. _CompositionalMultiphaseFluid: - -############################################ -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 :math:`N_c` components with :math:`z_c` being the total -mole fraction of component :math:`c`. The fluid can partition into a liquid phase, denoted :math:`\ell`, -and a vapor phase denoted by :math:`v`. Therefore, by taking into account the molar phase component -fractions, (which is the fraction of the molar mass of phase :math:`p` represented by component -:math:`c`), the following partition matrix establishes the component distribution within the two -phases: - -.. math:: - \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 :math:`x_c` is the mole fraction of component :math:`c` in the liquid phase and :math:`y_c` -is the mole fraction of component :math:`c` in the vapor phase. - -The fluid properties are updated through the following steps: - -1) The phase fractions (:math:`\nu_p`) and phase component fractions (:math:`x_c` and :math:`y_c`) are -computed as a function of pressure (:math:`p`), temperature (:math:`T`) and total component fractions -(:math:`z_c`). - -2) The phase densities (:math:`\rho_p`) and phase viscosities (:math:`\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 -:ref:`CompositionalMultiphaseFlow` 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 :math:`z_c` is stable -as a single phase at the current pressure :math:`p` and temperature :math:`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 :math:`z` is stable -at a specified pressure :math:`p` and temperature :math:`T` if and only if - -.. math:: - 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 :math:`y`, where :math:`\phi_i` denotes the fugacity -coefficient of component :math:`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 :math:`y_i = z_i/K_i` and :math:`y_i = z_iK_i` where :math:`K_i` are defined by - -.. math:: - K_i = \frac{P_{ci}}{p}\exp\left( 5.37( 1 + \omega_i ) \left(1-\frac{T_{ci}}{T}\right)\right) - -where :math:`P_{ci}` and :math:`T_{ci}` are respectively, the critical pressure and temperature of -component :math:`i` and :math:`\omega_i` is the accentric factor of component :math:`i`. - -The stability problem is solved by observing that a necessary condition is that :math:`g(y)` must -be non-negative at all its stationary points. The stationarity criterion can be expressed as - -.. math:: - \ln y_i + \ln \phi_i(y) - h_i = k \hspace{1cm} i=1,2,3,\ldots,N_c - -where :math:`h_i = \ln z_i + \ln \phi_i(z)` is a constant parameter dependent on the feed composition -:math:`z` and :math:`k` is an undetermined constant. This constant can be further incorporated into -the equation by defining the unnormalized trial phase moles :math:`Y_i` as - -.. math:: - Y_i = \exp(-k)y_i - -which reduces the stationarity criterion to - -.. math:: - \ln Y_i + \ln \phi_i(y) - h_i = 0 - -with the mole fractions :math:`y_i` related to the trial phase moles :math:`Y_i` by - -.. math:: - 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 :math:`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 :math:`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 - -.. math:: - T_{cp} = \frac{\sum_{i=1}^{N_c}T_{ci}V_{ci}z_{i}}{\sum_{i=1}^{N_c}V_{ci}z_{i}} - -where :math:`V_{ci}` and :math:`T_{ci}` are respectively the critical volume and temperature of component -:math:`i`. This is compared to the current temperature :math:`T` such that if :math:`T_{cp}`` 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: - -.. include:: /docs/sphinx/datastructure/CompositionalMultiphaseFluid.rst - -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 -========================= - -.. code-block:: xml - - - - - -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. - -.. _Petrowiki: https://petrowiki.spe.org/Phase_behavior_in_reservoir_simulation#Equation-of-state_models diff --git a/src/coreComponents/constitutive/docs/FluidModels.rst b/src/coreComponents/constitutive/docs/FluidModels.rst index 1aed1ce387e..4c5ad023db3 100644 --- a/src/coreComponents/constitutive/docs/FluidModels.rst +++ b/src/coreComponents/constitutive/docs/FluidModels.rst @@ -15,6 +15,6 @@ single fluids and fluid mixtures. BlackOilFluid - CompositionalMultiphaseFluid + compositional/CompositionalMultiphaseFluid CO2BrineFluid diff --git a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst new file mode 100644 index 00000000000..466ac8f8e5d --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst @@ -0,0 +1,12 @@ +.. _CompositionalMultiphaseFluid: + +############################################ +Compositional multiphase fluid model +############################################ + +.. toctree:: + :maxdepth: 1 + + Overview + Flash + References diff --git a/src/coreComponents/constitutive/docs/compositional/Flash.rst b/src/coreComponents/constitutive/docs/compositional/Flash.rst new file mode 100644 index 00000000000..b0606ffa325 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Flash.rst @@ -0,0 +1,8 @@ +.. _Flash: + +Step 1: Computation of the phase fractions and phase component fractions (flash) +================================================================================ + +.. toctree:: + StabilityTest + NegativeTwoPhaseFlash diff --git a/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst b/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst new file mode 100644 index 00000000000..ba10ad454bb --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst @@ -0,0 +1,72 @@ +.. _NegativeTwoPhaseFlash: + +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: + +.. math:: + \phi_{iL} = \phi_{iV} \hspace{1cm} i=1,2,3,\ldots,N_c + +where :math:`\phi_{iL}` is the fugacity of component :math:`i` in the liquid phase and :math:`\phi_{iL}` is the +fugacity of component :math:`i` in the vapor phase. + +Fugacities are functions of temperature, pressure, and composition: + +.. math:: + \phi_{iL} = \phi_{iL}(T, p, x_i) \hspace{1cm} i=1,2,3,\ldots,N_c + +and + +.. math:: + \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: + +.. math:: + K_i = \frac{y_i}{x_i} + +where :math:`x_i` is the mole fraction of component :math:`i` in the liquid phase and :math:`y_i` is the +mole fraction of component :math:`i` in the vapor phase. If we denote :math:`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: + +.. math:: + x_i = \frac{z_i}{1 + (K_i - 1)V} + +and + +.. math:: + y_i = \frac{K_i z_i}{1 + (K_i - 1)V} + +The value of :math:`V` corresponding to a given set of K-values is determined by solving the +so called Rachford and-Rice equation: + +.. math:: + 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 :math:`z_i` and :math:`K_i`, the Rachford-Rice equation is solved to determine the molar fraction of vapor, :math:`V`. This is initially solved using successive substitution, followed by Newton iterations once the residual is sufficiently reduced. + +#. After :math:`V` is calculated, the corresponding liquid and vapor mole fractions, :math:`x_i` and :math:`y_i`, are computed. + +#. These phase compositions are then used to calculate the component fugacities :math:`\phi_{iL}` and :math:`\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: + + .. math:: + \sum_{i=1}^{N_c} \left( \phi_{iL} - \phi_{iV} \right)^2 < \varepsilon + + where :math:`\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 :math:`t+1` are given by: + + .. math:: + K_i^{(t+1)} = K_i^{(t)} \frac{\phi_{iL}}{\phi_{iV}} diff --git a/src/coreComponents/constitutive/docs/compositional/Overview.rst b/src/coreComponents/constitutive/docs/compositional/Overview.rst new file mode 100644 index 00000000000..4602ed79ab9 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Overview.rst @@ -0,0 +1,40 @@ +.. _Overview: + +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 :math:`N_c` components with :math:`z_c` being the total +mole fraction of component :math:`c`. The fluid can partition into a liquid phase, denoted :math:`\ell`, +and a vapor phase denoted by :math:`v`. Therefore, by taking into account the molar phase component +fractions, (which is the fraction of the molar mass of phase :math:`p` represented by component +:math:`c`), the following partition matrix establishes the component distribution within the two +phases: + +.. math:: + \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 :math:`x_c` is the mole fraction of component :math:`c` in the liquid phase and :math:`y_c` +is the mole fraction of component :math:`c` in the vapor phase. + +The fluid properties are updated through the following steps: + +1) The phase fractions (:math:`\nu_p`) and phase component fractions (:math:`x_c` and :math:`y_c`) are +computed as a function of pressure (:math:`p`), temperature (:math:`T`) and total component fractions +(:math:`z_c`). + +2) The phase densities (:math:`\rho_p`) and phase viscosities (:math:`\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 +:ref:`CompositionalMultiphaseFlow` then moves on to assembling the accumulation and flux terms. diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst new file mode 100644 index 00000000000..32bdd33fd18 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -0,0 +1,14 @@ +.. _References + +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. + +.. _Petrowiki: https://petrowiki.spe.org/Phase_behavior_in_reservoir_simulation#Equation-of-state_models diff --git a/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst b/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst new file mode 100644 index 00000000000..03f4f0d1947 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst @@ -0,0 +1,75 @@ +.. _StabilityTest: + +Stability test +------------------------------- +The first step is to determine if the provided mixture with total molar fractions :math:`z_c` is stable +as a single phase at the current pressure :math:`p` and temperature :math:`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 :math:`z` is stable +at a specified pressure :math:`p` and temperature :math:`T` if and only if + +.. math:: + 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 :math:`y`, where :math:`\phi_i` denotes the fugacity +coefficient of component :math:`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 :math:`y_i = z_i/K_i` and :math:`y_i = z_iK_i` where :math:`K_i` are defined by + +.. math:: + K_i = \frac{P_{ci}}{p}\exp\left( 5.37( 1 + \omega_i ) \left(1-\frac{T_{ci}}{T}\right)\right) + +where :math:`P_{ci}` and :math:`T_{ci}` are respectively, the critical pressure and temperature of +component :math:`i` and :math:`\omega_i` is the accentric factor of component :math:`i`. + +The stability problem is solved by observing that a necessary condition is that :math:`g(y)` must +be non-negative at all its stationary points. The stationarity criterion can be expressed as + +.. math:: + \ln y_i + \ln \phi_i(y) - h_i = k \hspace{1cm} i=1,2,3,\ldots,N_c + +where :math:`h_i = \ln z_i + \ln \phi_i(z)` is a constant parameter dependent on the feed composition +:math:`z` and :math:`k` is an undetermined constant. This constant can be further incorporated into +the equation by defining the unnormalized trial phase moles :math:`Y_i` as + +.. math:: + Y_i = \exp(-k)y_i + +which reduces the stationarity criterion to + +.. math:: + \ln Y_i + \ln \phi_i(y) - h_i = 0 + +with the mole fractions :math:`y_i` related to the trial phase moles :math:`Y_i` by + +.. math:: + 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 :math:`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 :math:`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 + +.. math:: + T_{cp} = \frac{\sum_{i=1}^{N_c}T_{ci}V_{ci}z_{i}}{\sum_{i=1}^{N_c}V_{ci}z_{i}} + +where :math:`V_{ci}` and :math:`T_{ci}` are respectively the critical volume and temperature of component +:math:`i`. This is compared to the current temperature :math:`T` such that if :math:`T_{cp} Date: Tue, 24 Jun 2025 09:02:59 -0500 Subject: [PATCH 02/13] Remove link --- src/coreComponents/constitutive/docs/compositional/Overview.rst | 1 - .../constitutive/docs/compositional/References.rst | 2 -- 2 files changed, 3 deletions(-) diff --git a/src/coreComponents/constitutive/docs/compositional/Overview.rst b/src/coreComponents/constitutive/docs/compositional/Overview.rst index 4602ed79ab9..b2798445630 100644 --- a/src/coreComponents/constitutive/docs/compositional/Overview.rst +++ b/src/coreComponents/constitutive/docs/compositional/Overview.rst @@ -8,7 +8,6 @@ Phase behavior is modeled by an equation of state (EOS) and partitioning of comp 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 :math:`N_c` components with :math:`z_c` being the total mole fraction of component :math:`c`. The fluid can partition into a liquid phase, denoted :math:`\ell`, diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst index 32bdd33fd18..3d12a38867d 100644 --- a/src/coreComponents/constitutive/docs/compositional/References.rst +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -10,5 +10,3 @@ References - M. L. Michelsen, `The Isothermal Flash Problem. Part II. Phase-Split Calculation. `__, Fluid Phase Equilibria, vol. 9.1, pp. 21-40, 1982b. - -.. _Petrowiki: https://petrowiki.spe.org/Phase_behavior_in_reservoir_simulation#Equation-of-state_models From ea6b32a08c676a3366a4a6689648dc5f725cafd3 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Wed, 3 Jun 2026 16:51:52 -0500 Subject: [PATCH 03/13] Phase split calculations --- .../CompositionalMultiphaseFluid.rst | 19 +- .../docs/compositional/Density.rst | 46 ++++ .../docs/compositional/Enthalpy.rst | 34 +++ .../docs/compositional/EquationsOfState.rst | 224 ++++++++++++++++++ .../constitutive/docs/compositional/Flash.rst | 8 - .../compositional/ImmiscibleWaterFlash.rst | 86 +++++++ .../docs/compositional/KValueFlash.rst | 74 ++++++ .../docs/compositional/Miscellaneous.rst | 17 ++ .../docs/compositional/NegativeFlash.rst | 75 ++++++ .../compositional/NegativeTwoPhaseFlash.rst | 72 ------ .../docs/compositional/Overview.rst | 72 +++--- .../docs/compositional/PhaseSplit.rst | 47 ++++ .../docs/compositional/References.rst | 18 +- .../docs/compositional/StabilityTest.rst | 75 ------ .../docs/compositional/Viscosity.rst | 37 +++ 15 files changed, 710 insertions(+), 194 deletions(-) create mode 100644 src/coreComponents/constitutive/docs/compositional/Density.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/Enthalpy.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst delete mode 100644 src/coreComponents/constitutive/docs/compositional/Flash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/KValueFlash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst delete mode 100644 src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst delete mode 100644 src/coreComponents/constitutive/docs/compositional/StabilityTest.rst create mode 100644 src/coreComponents/constitutive/docs/compositional/Viscosity.rst diff --git a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst index 466ac8f8e5d..75adc853688 100644 --- a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst +++ b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst @@ -1,12 +1,23 @@ -.. _CompositionalMultiphaseFluid: +.. _CompositionalFluidModel: ############################################ Compositional multiphase fluid model ############################################ +Some text + .. toctree:: - :maxdepth: 1 + :maxdepth: 2 + + PhaseSplit + + Density + + Viscosity + + Enthalpy + + Miscellaneous - Overview - Flash References + diff --git a/src/coreComponents/constitutive/docs/compositional/Density.rst b/src/coreComponents/constitutive/docs/compositional/Density.rst new file mode 100644 index 00000000000..47e8aa30b33 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Density.rst @@ -0,0 +1,46 @@ +======= +Density +======= + +The accurate determination of phase density is critical for modelling fluid flow, as it directly impacts volumetric calculations, gravity drainage, and phase mobilities. + +Compositional Density Model +--------------------------- + +The compositional density model derives the molar volume of a phase using the compressibility factor, :math:`Z`, computed directly from the cubic Equation of State (e.g., Peng-Robinson). The molar density, :math:`\rho_{molar}`, is given by: + +.. math:: + \rho_{molar} = \frac{P}{Z R T} + +Because standard cubic equations of state notoriously under-predict liquid phase densities, a Peneloux volume shift correction is mathematically applied. The corrected molar volume, :math:`v_{corrected}`, is adjusted by a sum of component-specific volume shift parameters, :math:`c_i`, weighted by the phase composition: + +.. math:: + v_{corrected} = v_{EOS} - \sum_{i=1}^{N_c} x_i c_i + +The molar density is then inverted to mass density using the molecular weight of the phase mixture. This model is recommended for all hydrocarbon phases. + +* **Catalog Name:** ``CompositionalDensity`` +* **Parameters:** * ``componentVolumeShift``: Component-specific volume shift parameters. + +Phillips Brine Density Model +---------------------------- + +For aqueous phases containing dissolved salts, cubic equations of state often struggle to capture the complex ionic interactions. The Phillips Brine Density model employs a robust empirical correlation developed by Phillips et al. (1981). + +This model calculates the mass density of the brine as a direct polynomial and exponential function of pressure, temperature, and salinity (molality of dissolved salts). For low salinities, the model interpolates smoothly towards pure water properties to maintain physical continuity. This is highly recommended for geothermal reservoirs or saline aquifers. + +* **Catalog Name:** ``PhillipsBrineDensity`` +* **Parameters:** + * ``salinity``: Brine salinity. + * ``saltMolarWeight``: The molar weight of the salt component. + +Immiscible Water Density +------------------------ + +A simplified exponential density model is available for pure, immiscible water phases. It evaluates density based on isothermal compressibility and thermal expansion from a reference state: + +.. math:: + \rho = \rho_{ref} \exp(c_w (P - P_{ref})) \exp(-\alpha_w (T - T_{ref})) + +* **Catalog Name:** ``ImmiscibleWaterDensity`` +* **Parameters:** ``waterReferencePressure``, ``waterReferenceTemperature``, ``waterDensity``, ``waterCompressibility``, ``waterExpansionCoefficient``. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst new file mode 100644 index 00000000000..7be1e969ef7 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst @@ -0,0 +1,34 @@ +======== +Enthalpy +======== + +For non-isothermal simulations, such as thermal recovery processes or geothermal operations, the accurate calculation of phase enthalpies and internal energies is required to resolve the energy balance equation. + +Compositional Enthalpy Model +---------------------------- + +The enthalpy of a compositional fluid phase is computed using a rigorous thermodynamic departure function approach. The true enthalpy at a given pressure and temperature, :math:`H(P, T)`, is expressed as the sum of an ideal gas enthalpy at the same temperature, :math:`H_{ideal}(T)`, and an EoS-derived departure enthalpy, :math:`H_{departure}(P, T, x)`: + +.. math:: + H(P, T) = H_{ideal}(T) + H_{departure}(P, T, x) + +**Ideal Gas Contribution** + +The ideal gas enthalpy is derived by integrating a temperature-dependent specific heat capacity polynomial, :math:`C_p(T)`, from a defined reference temperature, :math:`T_{ref}`, to the current system temperature: + +.. math:: + C_p(\Delta T) = a_0 + a_1 \Delta T + a_2 \Delta T^2 + a_3 \Delta T^3 + a_4 \Delta T^4 + +where :math:`\Delta T = T - T_{ref}`. The integrated polynomial is added to a user-supplied reference enthalpy for each component. + +**Departure Function** + +The departure enthalpy accounts for real-gas behaviour and molecular interactions. It is calculated exactly from the Equation of State by evaluating the partial derivative of the EoS attractive mixture parameter, :math:`a`, with respect to temperature. + +This model is recommended whenever non-isothermal compositional physics are required. + +* **Catalog Name:** ``CompositionalEnthalpy`` +* **Parameters:** + * ``enthalpyReferenceTemperature``: The reference temperature :math:`T_{ref}`. + * ``referenceEnthalpy``: The enthalpy of each component at the reference temperature. + * ``componentHeatCapacityCoefficients``: The polynomial coefficients (:math:`a_0` to :math:`a_4`) for the heat capacity of each component. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst new file mode 100644 index 00000000000..cc701567f07 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst @@ -0,0 +1,224 @@ +Equations of state +================== + +An Equation of State (EoS) mathematically relates the pressure, volume, temperature, and composition of a fluid. In compositional modelling, the EoS serves a dual purpose: +1. It calculates the fugacity coefficients, :math:`\phi_i`, required to establish thermodynamic phase equilibrium. +2. It calculates the compressibility factor, :math:`Z`, required to determine the phase density. + +The fugacity coefficient of a component is derived from the exact thermodynamic relationship involving the integration of the EoS volume departure with respect to pressure. + +Cubic equation of state +----------------------- + +The Peng-Robinson (PR) and Soave-Redlich-Kwong (SRK) equations of state are widely used for hydrocarbon systems encompassing natural gas, gas condensates, and volatile oils. The fundamental pressure-volume-temperature (PVT) relationships for these equations are defined as: + +* Peng-Robinson (PR): + + .. math:: + P = \frac{RT}{v - b} - \frac{a(T)}{v(v + b) + b(v - b)} + +* Soave-Redlich-Kwong (SRK): + + .. math:: + P = \frac{RT}{v - b} - \frac{a(T)}{v(v + b)} + +In these PVT relationships, several intermediate parameters hold specific physical meanings: + +* :math:`v` represents the molar volume of the fluid phase. +* :math:`R` is the universal gas constant, serving as the fundamental scaling factor relating macroscopic thermodynamic state variables to the molar energy. +* :math:`a(T)` is the temperature-dependent attractive parameter, representing the intermolecular attractive forces between the fluid molecules. +* :math:`b` is the repulsive co-volume parameter, representing the effective physical volume occupied by the fluid molecules themselves. + +To efficiently solve these equations in a multicomponent mixture, they are rewritten into a generalized cubic equation for the compressibility factor, :math:`Z`: + +.. math:: + Z^3 + E_2 Z^2 + E_1 Z + E_0 = 0 + +where the coefficients are defined by the dimensionless mixture parameters :math:`A` and :math:`B`: + +.. math:: + E_2 = (\delta_1 + \delta_2 - 1)B - 1 + +.. math:: + E_1 = A + \delta_1 \delta_2 B^2 - (\delta_1 + \delta_2)B(B + 1) + +.. math:: + E_0 = -[AB + \delta_1 \delta_2 B^2(B + 1)] + +The constants :math:`\delta_1` and :math:`\delta_2` determine the specific EoS used: + +* Peng-Robinson (PR): :math:`\delta_1 = 1 + \sqrt{2}`, :math:`\delta_2 = 1 - \sqrt{2}` +* Soave-Redlich-Kwong (SRK): :math:`\delta_1 = 0`, :math:`\delta_2 = 1` + +The mixture parameters :math:`A` and :math:`B` are calculated using standard Van der Waals mixing rules over the phase mole fractions, :math:`x_i`: + +.. math:: + A = \sum_i \sum_j x_i x_j (1 - k_{ij}) \sqrt{A_i A_j} + +.. math:: + B = \sum_i x_i B_i + +where :math:`k_{ij}` are the binary interaction coefficients (BICs). The pure component parameters :math:`A_i` and :math:`B_i` are defined as: + +.. math:: + A_i = \Omega_a \frac{P_{r,i}}{T_{r,i}^2} \alpha_i + +.. math:: + B_i = \Omega_b \frac{P_{r,i}}{T_{r,i}} + +where :math:`P_{r,i} = P/P_{c,i}` and :math:`T_{r,i} = T/T_{c,i}` are the reduced pressure and temperature of component :math:`i`. + +The dimensionless constants :math:`\Omega_a` and :math:`\Omega_b` are derived mathematically by applying the critical point constraints (where the first and second pressure-volume derivatives vanish) to dictate the weight of the attractive and repulsive forces. The exact values implemented in the code are: + +* Peng-Robinson: :math:`\Omega_a = 0.457235529`, :math:`\Omega_b = 0.077796074` +* Soave-Redlich-Kwong: :math:`\Omega_a = 0.42748`, :math:`\Omega_b = 0.08664` + +The alpha function, :math:`\alpha_i`, accounts for the temperature dependence of the attractive parameter and relies on the acentric factor, :math:`\omega_i`, of each component. For both equations, it takes the form: + +.. math:: + \alpha_i = \left( 1 + m_i \left( 1 - \sqrt{T_{r,i}} \right) \right)^2 + +The parameter :math:`m_i` is computed differently depending on the chosen equation of state: + +* Peng-Robinson: + If :math:`\omega_i < 0.49`: + + .. math:: + m_i = 0.37464 + 1.54226 \omega_i - 0.26992 \omega_i^2 + + If :math:`\omega_i \ge 0.49`: + + .. math:: + m_i = 0.3796 + 1.485 \omega_i - 0.164423 \omega_i^2 + 0.016666 \omega_i^3 + +* Soave-Redlich-Kwong: + + .. math:: + m_i = 0.480 + 1.574 \omega_i - 0.176 \omega_i^2 + +Finally, the fugacity coefficient, :math:`\phi_i`, for each component is calculated using the derived roots for :math:`Z`: + +.. math:: + \ln \phi_i = (Z - 1) \frac{B_i}{B} - F - G \left( 2 \sum_j x_j (1 - k_{ij}) \sqrt{A_i A_j} - A \frac{B_i}{B} \right) E + +where the intermediate terms are defined as: + +.. math:: + E = \ln \left( \frac{Z + \delta_1 B}{Z + \delta_2 B} \right), \quad F = \ln(Z - B), \quad G = \frac{1}{(\delta_1 - \delta_2)B} + +Soreide-whitson equation of state +--------------------------------- + +The Soreide-Whitson model represents a specialized modification of the Peng-Robinson EoS tailored specifically for aqueous-hydrocarbon mixtures, such as those encountered in carbon sequestration (CO2-brine systems) or water-flooding in oil reservoirs. It adapts the standard cubic EoS by capturing the "salting-out" effect, where increased brine salinity reduces the solubility of hydrocarbons and CO2 in the aqueous phase. It is highly recommended whenever accurate modelling of gas dissolution in saline water is critical. + +This model introduces two primary modifications to the standard Peng-Robinson formulation, applying special treatment exclusively to the water component and its interactions within the aqueous phase: + +First, the standard alpha function, :math:`\alpha_i`, is replaced specifically for the water component by a salinity-dependent formulation. If :math:`c_{sw}` represents the brine salinity, the water alpha function is evaluated as: + +.. math:: + \alpha_w = \left[ 1 + 0.4530 \left( 1 - T_{r,w} \left( 1 - 0.0103 c_{sw}^{1.1} \right) \right) + 0.0034 \left( T_{r,w}^{-3} - 1 \right) \right]^2 + +where :math:`T_{r,w}` is the reduced temperature of water. + +Second, unlike standard equations of state that use constant binary interaction coefficients (BICs), this model calculates dynamic aqueous phase BICs, :math:`k_{ij}^{AQ}`, exclusively for pairs where one component is water and the other is a hydrocarbon or specific non-hydrocarbon. + +For water-hydrocarbon pairs, the interaction coefficient is modeled as a polynomial function of the hydrocarbon's reduced temperature, :math:`T_{r,i}`, its acentric factor, :math:`\omega_i`, and the brine salinity: + +.. math:: + k_{ij}^{AQ} = A_0 (1 + \alpha_0 c_{sw}) + A_1 T_{r,i} (1 + \alpha_1 c_{sw}) + A_2 T_{r,i}^2 (1 + \alpha_2 c_{sw}) + +The constants for this polynomial depend on the acentric factor. A special treatment is applied for methane, identified by an acentric factor :math:`\omega_i \le 0.02` (methane's acentric factor is approximately 0.0108). To achieve higher accuracy for methane-brine systems, a specifically tuned set of empirical constants is used instead of the generalized hydrocarbon constants. + +For standard hydrocarbons (:math:`\omega_i > 0.02`): + +* :math:`A_0 = 1.1120 - 1.7369 \omega_i^{-0.1}` +* :math:`A_1 = 1.1001 + 0.8360 \omega_i` +* :math:`A_2 = -0.15742 - 1.0988 \omega_i` +* :math:`\alpha_0 = 0.017407` +* :math:`\alpha_1 = 0.033516` +* :math:`\alpha_2 = 0.011478` + +For methane (:math:`\omega_i \le 0.02`): + +* :math:`A_0 = 1.616705 - 1.884534 \omega_i^{-0.1}` +* :math:`A_1 = 0.8157014 + 0.8723632 \omega_i` +* :math:`A_2 = -0.0887821 - 0.8767864 \omega_i` +* :math:`\alpha_0 = 0.09988448` +* :math:`\alpha_1 = 0.1485516` +* :math:`\alpha_2 = 0.2111324` + +For specific non-hydrocarbon components interacting with water, the dynamic BICs are calculated using explicitly tuned correlations based on the component type: + +For carbon dioxide (CO2): + +.. math:: + k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i} + A_2 a_2 + +where :math:`A_0 = -0.31092`, :math:`A_1 = 0.23580`, :math:`A_2 = -21.2566`, :math:`a_0 = 1.0 + 0.15587 c_{sw}^{0.7505}`, :math:`a_1 = 1.0 + 0.17837 c_{sw}^{0.979}`, and :math:`a_2 = \exp(-6.7222 T_{r,i} - c_{sw})`. + +For nitrogen (N2): + +.. math:: + k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i} + +where :math:`A_0 = -1.70235`, :math:`A_1 = 0.44338`, :math:`a_0 = 1.0 + 0.025587 c_{sw}^{0.75}`, and :math:`a_1 = 1.0 + 0.08126 c_{sw}^{0.75}`. + +For hydrogen sulfide (H2S): + +.. math:: + k_{ij}^{AQ} = A_0 + A_1 T_{r,i} + +where :math:`A_0 = -0.20441` and :math:`A_1 = 0.23426`. + +For hydrogen (H2): + +.. math:: + k_{ij}^{AQ} = D_0 (1 + a_0 c_{sw}) + D_1 T_{r,i} (1 + a_1 c_{sw}) + D_2 \exp(D_3 T_{r,i}) + +where :math:`D_0 = -2.11917`, :math:`D_1 = 0.14888`, :math:`D_2 = -13.01835`, :math:`D_3 = -0.43946`, :math:`a_0 = -0.0226322`, and :math:`a_1 = -0.0044736`. + +Component naming requirements +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Because the Soreide-Whitson EoS applies highly specific correlations based on the exact type of component interacting with water, it relies on the component's assigned name to identify its type. This name matching is case-insensitive. + +To trigger the correct dynamic binary interaction coefficients, you must use the following exact names for the special components: + +* Water: ``H2O`` +* Carbon dioxide: ``CO2`` +* Nitrogen: ``N2`` +* Hydrogen sulfide: ``H2S`` +* Hydrogen: ``H2`` + +If a component's name does not match any of these predefined strings (for example, if you name it "CarbonDioxide", "C1", or "Methane"), the model will automatically default to categorising it as a generic hydrocarbon (``hc``). While this is the intended and correct behavior for actual hydrocarbons (where methane is safely identified via its acentric factor), using an unrecognized name for a special non-hydrocarbon (like CO2) will cause the simulator to mistakenly apply the generic hydrocarbon-water correlation instead of its highly tuned specific correlation. + +Model parameters +---------------- + +Equations of state are assigned per-phase within the specific fluid model XML block (e.g., ````). + +The ``equationsOfState`` attribute takes a list specifying the EoS type for each phase corresponding to the order defined in the ``phaseNames`` list. Standard component properties necessary for the EoS calculations must also be provided in matching component order. + +Here is an example demonstrating how to specify the parameters for a two-phase system using the Soreide-Whitson EoS for the aqueous phase and Peng-Robinson for the gas phase: + +.. code-block:: xml + + + +* ``equationsOfState``: List specifying the EoS type (e.g., ``PengRobinson``, ``SoaveRedlichKwong``, ``SoreideWhitson``) corresponding to each phase. +* ``componentCriticalPressure``: Critical pressure of each component [Pa]. +* ``componentCriticalTemperature``: Critical temperature of each component [K]. +* ``componentAcentricFactor``: Acentric factor of each component. +* ``componentVolumeShift``: Volume shift parameter for each component (optional). +* ``componentBinaryCoeff``: Flattened :math:`N \times N` matrix of binary interaction coefficients (optional). +* ``salinity``: Brine salinity [mol/kg] utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Flash.rst b/src/coreComponents/constitutive/docs/compositional/Flash.rst deleted file mode 100644 index b0606ffa325..00000000000 --- a/src/coreComponents/constitutive/docs/compositional/Flash.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _Flash: - -Step 1: Computation of the phase fractions and phase component fractions (flash) -================================================================================ - -.. toctree:: - StabilityTest - NegativeTwoPhaseFlash diff --git a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst new file mode 100644 index 00000000000..1e1ebc0e4e6 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst @@ -0,0 +1,86 @@ +Immiscible water flash +====================== + +The immiscible water flash is an efficient three-phase equilibrium model designed for systems containing an aqueous phase alongside hydrocarbon liquid and vapour phases. Instead of solving a fully coupled three-phase equation of state problem, this model vastly simplifies the phase split by assuming water forms a strictly pure aqueous phase and is completely immiscible in the hydrocarbon phases, while hydrocarbons are entirely insoluble in the aqueous phase. + +Methodology +----------- + +In this formulation, the mixture composition is explicitly separated into a water component and a hydrocarbon mixture before any thermodynamic equilibrium calculations occur. The water component is identified using its assigned component name (e.g., ``H2O``). + +Given a total feed mole fraction :math:`z_i` for each component, the mole fraction of the aqueous phase, :math:`V_{aq}`, is determined to be exactly equal to the total feed mole fraction of the water component: + +.. math:: + V_{aq} = z_{w} + +The composition of the aqueous phase is fixed at 100% water. The remaining feed defines the total hydrocarbon mole fraction, :math:`z_{hc}`: + +.. math:: + z_{hc} = 1.0 - z_{w} + +The feed composition of the remaining hydrocarbon mixture is then normalised: + +.. math:: + z_{i,hc} = \frac{z_i}{z_{hc}} + +for all non-water components, while the normalised hydrocarbon feed fraction for water is set to zero. + +Equilibrium calculation +----------------------- + +During the simulation step, if the total hydrocarbon fraction :math:`z_{hc}` is negligibly small, the model bypasses the flash and assigns the system as a single-phase aqueous fluid. + +If hydrocarbons are present, the model performs a rigorous two-phase negative flash exclusively on the normalised hydrocarbon mixture (:math:`z_{i,hc}`) using the chosen equations of state for the liquid and vapour phases. + +This embedded flash calculation determines the intermediate hydrocarbon vapour fraction, :math:`V_{V,hc}`, and the corresponding compositions of the hydrocarbon liquid (:math:`x_i`) and hydrocarbon vapour (:math:`y_i`) phases using the standard Rachford-Rice equation and iterative fugacity updates. + +Once the two-phase hydrocarbon flash converges, the overall phase fractions for the entire three-phase system are scaled back using the total hydrocarbon mole fraction: + +.. math:: + V_{V} = z_{hc} V_{V,hc} + +.. math:: + V_{L} = z_{hc} (1.0 - V_{V,hc}) + +The derivatives of the phase fractions and compositions with respect to pressure, temperature, and total composition are analytically transformed using the chain rule to account for this scaling and normalisation. + +Recommendations and pitfalls +---------------------------- + +The immiscible water flash is highly computationally efficient and is recommended for multi-phase reservoir simulations where water solubility in the hydrocarbon phases (and vice versa) is negligible. By avoiding a full three-phase fugacity solver, it drastically reduces simulation time and improves non-linear solver stability. + +Pitfalls: Because this model strictly enforces total immiscibility, it cannot accurately model scenarios where mutual solubility is physically significant, such as high-pressure CO2 injection into saline aquifers (where CO2 dissolves heavily into brine, and water vaporises into the gas stream). For such complex mutual solubility problems, using a two-phase flash paired with the Soreide-Whitson equation of state might be required. + +Model parameters +---------------- + +The immiscible water flash fluid models are assigned using specific XML blocks that pair the three-phase flash model with distinct viscosity and density models for the hydrocarbon and aqueous phases. + +Here is an example demonstrating how to specify the parameters for a three-phase system using the immiscible water flash paired with the Lohrenz-Bray-Clark viscosity model for hydrocarbons and a distinct constant-compressibility formulation for the pure water phase: + +.. code-block:: xml + + + +* catalog options: ``CompositionalThreePhaseLohrenzBrayClarkViscosity`` +* parameters: + * ``equationsOfState``: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment. + * ``waterReferencePressure``: The reference pressure for water density and viscosity [Pa]. + * ``waterReferenceTemperature``: The reference temperature for water density and viscosity [K] (optional, default: 293.15). + * ``waterDensity``: The water density at the reference pressure and temperature [kg/m^3]. + * ``waterViscosity``: The water viscosity at the reference pressure and temperature [Pa.s]. + * ``waterCompressibility``: The constant isothermal compressibility of water [1/Pa]. + * ``waterExpansionCoefficient``: The volumetric coefficient of thermal expansion of water [1/K] (optional, default: 0.0). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst new file mode 100644 index 00000000000..291b084beb6 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst @@ -0,0 +1,74 @@ +K-value flash +============= + +The K-value flash model offers a computationally simplified alternative to the rigorous equation of state flash calculation. Rather than iteratively evaluating complex equations of state to find fugacity coefficients, this model determines phase splits using pre-tabulated K-values. + +The model supports both two-phase and three-phase variants. The two-phase flash calculates equilibrium partitioning between a liquid and a vapour phase. The three-phase flash introduces an additional aqueous phase, allowing for the predefined partitioning of components into water alongside the hydrocarbon phases. + +Methodology +----------- + +In this formulation, the equilibrium ratio for each component in the secondary phases relative to the primary liquid phase, :math:`K_i = y_i / x_i`, is treated strictly as a predefined function of pressure and temperature. + +Users input these functions as either tables or symbolic functions. To optimize performance during the simulation, the model discretizes these functions during initialization to construct a multi-dimensional hypercube of K-values spanning the phase, component, pressure, and temperature dimensions. + +The pressure and temperature coordinates for this hypercube grid can be explicitly chosen and specified by the user. If they are omitted, the model automatically inspects the provided table functions, extracts all unique coordinate points from the underlying data, and constructs the grid dynamically. If only a single point is found, an artificial offset is added to permit interpolation. + +Equilibrium calculation +----------------------- + +During the simulation step, the model performs bilinear interpolation on the pre-built hypercube to evaluate the K-values and their partial derivatives at the current pressure and temperature. + +Because the K-values are explicitly known, the iterative fugacity updates are bypassed. Instead, the vapour fraction, :math:`V`, is found by directly solving the Rachford-Rice equation using successive substitution and Newton iterations: + +.. math:: + \sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V(K_i - 1)} = 0 + +Once the vapour fraction is determined, the phase compositions are directly evaluated: + +.. math:: + x_i = \frac{z_i}{1 + V(K_i - 1)}, \quad y_i = K_i x_i + +This direct solution approach intrinsically supports a "negative" flash. The mathematical solver allows the vapour fraction, :math:`V`, to temporarily converge to values outside the physically meaningful bounds of :math:`[0, 1]`. If the resulting :math:`V \le 0` or :math:`V \ge 1`, the system cleanly truncates the compositions to represent a single-phase liquid or single-phase vapour solution, respectively. + +Recommendations and pitfalls +---------------------------- + +The K-value flash is highly computationally efficient and is recommended for thermal recovery processes (like steam injection) or relatively stable black-oil systems where fluid compositions remain within a narrow, predictable envelope. + +Pitfalls: In reality, thermodynamic K-values are heavily dependent on the overall mixture composition, not just pressure and temperature. Relying on tabulated K-values implies assuming a fixed compositional path. Care should be taken when using this model for highly variable compositional processes such as miscible gas injection (e.g., CO2 or hydrocarbon gas flooding), or strong depletion of gas condensates, where composition dependency is critical to accurate phase behaviour. + +Model parameters +---------------- + +The K-value flash fluid models are assigned using specific XML blocks that pair the flash model with viscosity and density models. + +Here is an example demonstrating how to specify the parameters for a two-phase system using the K-value flash paired with the Phillips brine model: + +.. code-block:: xml + + + +* catalog options: ``CompositionalKValueLohrenzBrayClarkViscosity``, ``CompositionalKValuePhillipsBrine`` +* parameters: + * ``kValueTables``: List of function names linking to the pre-tabulated K-values. For an :math:`N`-phase system, this requires :math:`(N-1) \times N_c` function names. + * ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). + * ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst new file mode 100644 index 00000000000..b82c7442a8f --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst @@ -0,0 +1,17 @@ +============= +Miscellaneous +============= + +Unit and Variable Conversions +----------------------------- + +Thermodynamic models, including Equations of State and Flash calculations, are fundamentally derived from molar balances and molar properties. Consequently, the internal computations within the compositional fluid framework are strictly executed using molar variables (e.g., molar fractions, molar densities). + +**Mass Formulation Support** + +While the core thermodynamics rely on molar formulations, many reservoir simulators frame their governing conservation equations in terms of mass to ensure strict mass conservation across numerical schemes. + +To bridge this gap seamlessly, the fluid model supports an automated conversion layer: +1. **Entry Conversion:** Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. +2. **Internal Computation:** All stability, flash, density, and viscosity calculations proceed purely in molar space. +3. **Exit Transformation:** Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst new file mode 100644 index 00000000000..2068968eb75 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst @@ -0,0 +1,75 @@ +Negative flash +============== + +The negative flash is a robust, rigorous approach to solving the phase equilibrium problem using an Equation of State (EoS). + +Stability analysis +------------------ + +Before executing a flash calculation, the model must ascertain if the homogeneous single-phase mixture is thermodynamically unstable at the specified pressure :math:`P` and temperature :math:`T`. This is performed via a tangent plane distance (TPD) analysis developed by Michelsen. The TPD assesses whether the Gibbs free energy of the mixture can be minimised by splitting into multiple phases. A phase with feed composition :math:`z` is stable if and only if: + +.. math:: + g(\mathbf{y}) = \sum_{i=1}^{N_c} y_i \left( \ln y_i + \ln \phi_i(\mathbf{y}) - \ln z_i - \ln \phi_i(\mathbf{z}) \right) \ge 0 + +for any permissible trial composition :math:`\mathbf{y}`, where :math:`\phi_i` denotes the fugacity coefficient of component :math:`i`. + +To determine stability, testing is initiated from basic starting points using Wilson's K-values to generate a lighter trial mixture (:math:`y_i = z_i / K_i`) and a heavier trial mixture (:math:`y_i = z_i K_i`). Wilson's K-values are defined as: + +.. math:: + K_i = \frac{P_{ci}}{P} \exp \left( 5.37 (1 + \omega_i) \left( 1 - \frac{T_{ci}}{T} \right) \right) + +where :math:`P_{ci}`, :math:`T_{ci}`, and :math:`\omega_i` are the critical pressure, critical temperature, and acentric factor of component :math:`i`, respectively. + +The stationarity condition, :math:`\ln Y_i + \ln \phi_i(\mathbf{y}) - h_i = 0` (where :math:`h_i = \ln z_i + \ln \phi_i(\mathbf{z})` and :math:`Y_i` are unnormalized trial phase moles), is solved using successive substitution. If both initial states converge to a solution where :math:`g(\mathbf{y}) \ge 0`, the mixture is stable; otherwise, it is unstable. + +Phase labeling +-------------- + +If the stability test confirms a single stable phase, the two-phase flash calculation is bypassed. To appropriately label the single phase as liquid or vapour for relative permeability evaluations, the model utilises a "Li-temperature" correlation. This pseudo-critical temperature is weighted by the critical volumes (:math:`V_{ci}`): + +.. math:: + T_{cp} = \frac{\sum_{i=1}^{N_c} T_{ci} V_{ci} z_i}{\sum_{i=1}^{N_c} V_{ci} z_i} + +If :math:`T_{cp} < T`, the mixture is logically categorised as a vapour; otherwise, it is labelled as a liquid. + +Flash calculation +----------------- + +If the mixture is deemed unstable, the system must determine the phase split by ensuring thermodynamic equilibrium, meaning the fugacities in the liquid and vapour phases must be equal (:math:`\phi_{iL} = \phi_{iV}`). + +Based on the material balance and the equilibrium constants (:math:`K_i = y_i/x_i`), the mole fractions in the liquid (:math:`x_i`) and vapour (:math:`y_i`) phases are given by: + +.. math:: + x_i = \frac{z_i}{1 + (K_i - 1)V}, \quad y_i = \frac{K_i z_i}{1 + (K_i - 1)V} + +The vapour molar fraction, :math:`V`, is determined by solving the Rachford-Rice equation alongside the EoS: + +.. math:: + F(V) = \sum_{i=1}^{N_c} (x_i - y_i) = \sum_{i=1}^{N_c} \frac{z_i(1 - K_i)}{1 + (K_i - 1)V} = 0 + +The flash algorithm: + +1. Initialize: An initial set of K-values is chosen using Wilson's formula. +2. Solve for V: The Rachford-Rice equation is solved for the vapour fraction :math:`V` (using successive substitution, followed by Newton iterations). +3. Compute compositions: The corresponding liquid (:math:`x_i`) and vapour (:math:`y_i`) mole fractions are computed. +4. Evaluate EoS: These compositions are used to calculate the component fugacities :math:`\phi_{iL}` and :math:`\phi_{iV}` via the Equation of State. +5. Check convergence: Convergence is reached when :math:`\sum_{i=1}^{N_c} (\phi_{iL} - \phi_{iV})^2 < \epsilon`. +6. Update K-values: If not converged, the algorithm employs successive substitution iterations, constantly updating the K-values using fugacity coefficients derived from the EoS: + +.. math:: + K_i^{(new)} = K_i^{(old)} \frac{\phi_{iL}}{\phi_{iV}} = K_i^{(old)} \exp \left( \ln \phi_i^L - \ln \phi_i^V \right) + +The term "negative" flash arises because the algorithm temporarily permits the vapour fraction, :math:`V`, to converge to values slightly outside the physically meaningful bounds of :math:`[0, 1]`. This mathematical relaxation prevents the solver from getting artificially trapped at phase boundaries, vastly improving convergence robustness near the critical point or saturation envelopes. Upon convergence, if :math:`V` is negative or greater than unity, the system truncates to a single-phase solution. + +Recommendations +--------------- + +The negative two-phase flash is the recommended standard for rigorous compositional modelling of miscible gas injection, primary depletion of volatile oils, and scenarios where fluid compositions change drastically over time. An extended three-phase variant is available for immiscible water systems. + +* Catalog name: ``TwoPhase``, ``ThreePhase`` +* Parameters: + * ``stabilityThreshold``: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). + * ``stabilityTolerance``: Tolerance for stationarity in the stability test. + * ``stabilityMaxIterations``: Maximum successive substitution steps for stability analysis. + * ``flashTolerance``: Convergence tolerance for the fugacity ratio error. + * ``flashMaxIterations``: Maximum successive substitution steps for the flash solve. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst b/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst deleted file mode 100644 index ba10ad454bb..00000000000 --- a/src/coreComponents/constitutive/docs/compositional/NegativeTwoPhaseFlash.rst +++ /dev/null @@ -1,72 +0,0 @@ -.. _NegativeTwoPhaseFlash: - -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: - -.. math:: - \phi_{iL} = \phi_{iV} \hspace{1cm} i=1,2,3,\ldots,N_c - -where :math:`\phi_{iL}` is the fugacity of component :math:`i` in the liquid phase and :math:`\phi_{iL}` is the -fugacity of component :math:`i` in the vapor phase. - -Fugacities are functions of temperature, pressure, and composition: - -.. math:: - \phi_{iL} = \phi_{iL}(T, p, x_i) \hspace{1cm} i=1,2,3,\ldots,N_c - -and - -.. math:: - \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: - -.. math:: - K_i = \frac{y_i}{x_i} - -where :math:`x_i` is the mole fraction of component :math:`i` in the liquid phase and :math:`y_i` is the -mole fraction of component :math:`i` in the vapor phase. If we denote :math:`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: - -.. math:: - x_i = \frac{z_i}{1 + (K_i - 1)V} - -and - -.. math:: - y_i = \frac{K_i z_i}{1 + (K_i - 1)V} - -The value of :math:`V` corresponding to a given set of K-values is determined by solving the -so called Rachford and-Rice equation: - -.. math:: - 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 :math:`z_i` and :math:`K_i`, the Rachford-Rice equation is solved to determine the molar fraction of vapor, :math:`V`. This is initially solved using successive substitution, followed by Newton iterations once the residual is sufficiently reduced. - -#. After :math:`V` is calculated, the corresponding liquid and vapor mole fractions, :math:`x_i` and :math:`y_i`, are computed. - -#. These phase compositions are then used to calculate the component fugacities :math:`\phi_{iL}` and :math:`\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: - - .. math:: - \sum_{i=1}^{N_c} \left( \phi_{iL} - \phi_{iV} \right)^2 < \varepsilon - - where :math:`\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 :math:`t+1` are given by: - - .. math:: - K_i^{(t+1)} = K_i^{(t)} \frac{\phi_{iL}}{\phi_{iV}} diff --git a/src/coreComponents/constitutive/docs/compositional/Overview.rst b/src/coreComponents/constitutive/docs/compositional/Overview.rst index b2798445630..7e99cc8c6a0 100644 --- a/src/coreComponents/constitutive/docs/compositional/Overview.rst +++ b/src/coreComponents/constitutive/docs/compositional/Overview.rst @@ -1,39 +1,57 @@ -.. _Overview: - Overview -========================= +======== + +The compositional fluid models implemented in this framework provide a robust and flexible approach for modelling the thermodynamic behaviour of multicomponent, multiphase fluid mixtures. These models are intended to support complex subsurface flow simulations, including primary hydrocarbon recovery, miscible gas injection, carbon sequestration, and geothermal energy extraction. -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. +Fluid composition and phase partitioning +---------------------------------------- -In this model the fluid is described by :math:`N_c` components with :math:`z_c` being the total -mole fraction of component :math:`c`. The fluid can partition into a liquid phase, denoted :math:`\ell`, -and a vapor phase denoted by :math:`v`. Therefore, by taking into account the molar phase component -fractions, (which is the fraction of the molar mass of phase :math:`p` represented by component -:math:`c`), the following partition matrix establishes the component distribution within the two -phases: +In this model, the fluid is described by :math:`N_c` components, with :math:`z_c` being the total mole fraction of component :math:`c`. The fluid can partition into up to three phases: a liquid phase (denoted :math:`\ell`), a vapour phase (denoted :math:`v`), and potentially an aqueous phase (denoted :math:`a`). + +Therefore, by taking into account the molar phase component fractions (which is the fraction of the molar mass of phase :math:`p` represented by component :math:`c`), the following partition matrix establishes the component distribution within a two-phase (liquid-vapour) system: .. math:: + \begin{bmatrix} - x_{1} & x_{2} & x_{3} & \cdots & x_{N_c} \\ - y_{1} & y_{2} & y_{3} & \cdots & y_{N_c} \\ + x_1 & x_2 & x_3 & \cdots & x_{N_c} \\ + y_1 & y_2 & y_3 & \cdots & y_{N_c} \end{bmatrix} -where :math:`x_c` is the mole fraction of component :math:`c` in the liquid phase and :math:`y_c` -is the mole fraction of component :math:`c` in the vapor phase. +where :math:`x_c` is the mole fraction of component :math:`c` in the liquid phase and :math:`y_c` is the mole fraction of component :math:`c` in the vapour phase. + +If a third aqueous phase is also present, with component mole fractions denoted by :math:`w_c`, this partition matrix naturally extends to three rows to fully describe the component distribution across all three phases: + +.. math:: + + \begin{bmatrix} + x_1 & x_2 & x_3 & \cdots & x_{N_c} \\ + y_1 & y_2 & y_3 & \cdots & y_{N_c} \\ + w_1 & w_2 & w_3 & \cdots & w_{N_c} + \end{bmatrix} + +Structure of the fluid property calculations +-------------------------------------------- + +The evaluation of fluid properties follows a sequential thermodynamic structure designed to ensure consistency across all phase properties. At a given pressure, temperature, and overall composition, the framework executes the following overarching evaluations: + +1. State Conversion: If the system is formulated using mass variables, the overall mass fractions are converted into overall mole fractions. +2. Phase Split (Flash Calculation): The thermodynamic stability of the mixture is evaluated. If unstable, a flash calculation splits the mixture into separate phases (e.g., liquid, vapour, and potentially an aqueous phase), determining the phase fractions and the composition of each individual phase. +3. Phase Properties: Using the computed phase compositions, the model evaluates properties for each specific phase: + + * Density: Evaluated typically via a cubic Equation of State (EoS) combined with volume shift corrections. + * Viscosity: Evaluated using established compositional correlations or empirical models. + * Thermal Properties: If the simulation is non-isothermal, the phase enthalpies and internal energies are computed using ideal gas heat capacities and EoS departure functions. -The fluid properties are updated through the following steps: +4. Total Fluid Properties: The individual phase properties are homogenised to yield total fluid properties (e.g., total density). Derivatives with respect to pressure, temperature, and composition are rigorously propagated using the chain rule. -1) The phase fractions (:math:`\nu_p`) and phase component fractions (:math:`x_c` and :math:`y_c`) are -computed as a function of pressure (:math:`p`), temperature (:math:`T`) and total component fractions -(:math:`z_c`). +Structure of this Documentation +------------------------------- -2) The phase densities (:math:`\rho_p`) and phase viscosities (:math:`\mu_p`) are computed as a function -of pressure, temperature and the updated phase component fractions. +This documentation is divided into the following sections: -After calculating the phase fractions, phase component fractions, phase densities, phase viscosities, -and their derivatives with respect to pressure, temperature, and component fractions, the -:ref:`CompositionalMultiphaseFlow` then moves on to assembling the accumulation and flux terms. +* `Phase Split `_: Details the thermodynamic stability and flash calculations (Negative Flash and K-Value Flash), including the underlying Equations of State. +* `Density `_: Describes the models used to compute molar and mass phase densities. +* `Viscosity `_: Details the models used to evaluate the resistance to flow for each phase. +* `Enthalpy `_: Explains the calculation of thermal properties. +* `Miscellaneous `_: Outlines internal unit conversions and mass-to-molar variable handling. +* `References `_: Provides the bibliography for the scientific correlations implemented in the code. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst b/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst new file mode 100644 index 00000000000..54d970ccdb8 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst @@ -0,0 +1,47 @@ +Phase split +=========== + +The primary purpose of the phase split (or flash) calculation is to determine the number of fluid phases present at equilibrium, the fraction of the total mixture occupying each phase, and the detailed molar composition of these phases. This is a fundamental requirement for compositional modelling, as the physical properties (density, viscosity) of a fluid depend heavily on its phase state. + +Phase split equations +--------------------- + +The phase split relies on two fundamental physical principles: the conservation of mass (material balance) and thermodynamic equilibrium. + +**Material Balance (The Rachford-Rice Equation)** + +For a two-phase system (liquid and vapour), the overall mole fraction of a component, :math:`z_i`, must equal the sum of its molar contributions in the liquid phase, :math:`x_i`, and the vapour phase, :math:`y_i`. Introducing the vapour phase fraction, :math:`V`, and the equilibrium ratio (K-value), :math:`K_i = y_i / x_i`, the material balance is expressed by the Rachford-Rice equation: + +.. math:: + \sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V (K_i - 1)} = 0 + +The root of this equation provides the vapour fraction :math:`V`. The phase compositions are subsequently derived via: + +.. math:: + x_i = \frac{z_i}{1 + V(K_i - 1)} + +.. math:: + y_i = K_i x_i + +**Thermodynamic Equilibrium** + +At thermodynamic equilibrium, the chemical potential, or equivalently the fugacity, of each component must be identical across all existing phases: + +.. math:: + f_i^L(P, T, x) = f_i^V(P, T, y) + +Using the fugacity coefficient, :math:`\phi_i`, where :math:`f_i = x_i \phi_i P`, the equilibrium condition dictates the K-values: + +.. math:: + K_i = \frac{\phi_i^L}{\phi_i^V} + +The phase split models in this framework either rigorously solve for these fugacity coefficients using an Equation of State or use simplified, pre-tabulated K-values. + +.. include:: EquationsOfState.rst + +.. include:: NegativeFlash.rst + +.. include:: KValueFlash.rst + +.. include:: ImmiscibleWaterFlash.rst + diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst index 3d12a38867d..0dc072d143e 100644 --- a/src/coreComponents/constitutive/docs/compositional/References.rst +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -1,12 +1,14 @@ -.. _References - +========== References ========== -- M. L. Michelsen, `The Isothermal Flash Problem. Part I. Stability. - `__, Fluid Phase Equilibria, - vol. 9.1, pp. 1-19, 1982a. +The compositional fluid models implemented in this framework are based heavily on established thermodynamic and empirical literature. The following references document the mathematical correlations and mixing rules utilised: -- M. L. Michelsen, `The Isothermal Flash Problem. Part II. Phase-Split Calculation. - `__, Fluid Phase Equilibria, - vol. 9.1, pp. 21-40, 1982b. +* **Stiel, L. I., & Thodos, G. (1961).** The viscosity of nonpolar gases at normal pressures. *AIChE Journal*, 7(4), 611-615. +* **Herning, F., & Zipperer, L. (1936).** Calculation of the viscosity of technical gas mixtures from the viscosity of individual gases. *Gas u. Wasserfach*, 79(49), 69. +* **Wilke, C. R. (1950).** A viscosity equation for gas mixtures. *The Journal of Chemical Physics*, 18(4), 517-519. +* **Brokaw, R. S. (1968).** Viscosity of gas mixtures. *National Aeronautics and Space Administration*. +* **Lohrenz, J., Bray, B. G., & Clark, C. R. (1964).** Calculating viscosities of reservoir fluids from their compositions. *Journal of Petroleum Technology*, 16(10), 1171-1176. +* **Phillips, S. L., Igbene, A., Fair, J. A., Ozbek, H., & Tavana, M. (1981).** A technical databook for geothermal energy utilization. *Lawrence Berkeley National Laboratory*. +* **Soreide, I., & Whitson, C. H. (1992).** Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2S with tie water. *Fluid Phase Equilibria*, 77, 217-240. +* **Peneloux, A., Rauzy, E., & Freze, R. (1982).** A consistent correction for Redlich-Kwong-Soave volumes. *Fluid Phase Equilibria*, 8(1), 7-23. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst b/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst deleted file mode 100644 index 03f4f0d1947..00000000000 --- a/src/coreComponents/constitutive/docs/compositional/StabilityTest.rst +++ /dev/null @@ -1,75 +0,0 @@ -.. _StabilityTest: - -Stability test -------------------------------- -The first step is to determine if the provided mixture with total molar fractions :math:`z_c` is stable -as a single phase at the current pressure :math:`p` and temperature :math:`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 :math:`z` is stable -at a specified pressure :math:`p` and temperature :math:`T` if and only if - -.. math:: - 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 :math:`y`, where :math:`\phi_i` denotes the fugacity -coefficient of component :math:`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 :math:`y_i = z_i/K_i` and :math:`y_i = z_iK_i` where :math:`K_i` are defined by - -.. math:: - K_i = \frac{P_{ci}}{p}\exp\left( 5.37( 1 + \omega_i ) \left(1-\frac{T_{ci}}{T}\right)\right) - -where :math:`P_{ci}` and :math:`T_{ci}` are respectively, the critical pressure and temperature of -component :math:`i` and :math:`\omega_i` is the accentric factor of component :math:`i`. - -The stability problem is solved by observing that a necessary condition is that :math:`g(y)` must -be non-negative at all its stationary points. The stationarity criterion can be expressed as - -.. math:: - \ln y_i + \ln \phi_i(y) - h_i = k \hspace{1cm} i=1,2,3,\ldots,N_c - -where :math:`h_i = \ln z_i + \ln \phi_i(z)` is a constant parameter dependent on the feed composition -:math:`z` and :math:`k` is an undetermined constant. This constant can be further incorporated into -the equation by defining the unnormalized trial phase moles :math:`Y_i` as - -.. math:: - Y_i = \exp(-k)y_i - -which reduces the stationarity criterion to - -.. math:: - \ln Y_i + \ln \phi_i(y) - h_i = 0 - -with the mole fractions :math:`y_i` related to the trial phase moles :math:`Y_i` by - -.. math:: - 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 :math:`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 :math:`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 - -.. math:: - T_{cp} = \frac{\sum_{i=1}^{N_c}T_{ci}V_{ci}z_{i}}{\sum_{i=1}^{N_c}V_{ci}z_{i}} - -where :math:`V_{ci}` and :math:`T_{ci}` are respectively the critical volume and temperature of component -:math:`i`. This is compared to the current temperature :math:`T` such that if :math:`T_{cp} Date: Wed, 3 Jun 2026 17:12:15 -0500 Subject: [PATCH 04/13] Fix index --- .../docs/CompositionalTwoPhaseFluid.rst | 231 ------------------ .../CompositionalMultiphaseFluid.rst | 2 + 2 files changed, 2 insertions(+), 231 deletions(-) delete mode 100644 src/coreComponents/constitutive/docs/CompositionalTwoPhaseFluid.rst diff --git a/src/coreComponents/constitutive/docs/CompositionalTwoPhaseFluid.rst b/src/coreComponents/constitutive/docs/CompositionalTwoPhaseFluid.rst deleted file mode 100644 index 4bdde45e0c3..00000000000 --- a/src/coreComponents/constitutive/docs/CompositionalTwoPhaseFluid.rst +++ /dev/null @@ -1,231 +0,0 @@ -.. _CompositionalMultiphaseFluid: - -############################################ -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 :math:`N_c` components with :math:`z_c` being the total -mole fraction of component :math:`c`. The fluid can partition into a liquid phase, denoted :math:`\ell`, -and a vapor phase denoted by :math:`v`. Therefore, by taking into account the molar phase component -fractions, (which is the fraction of the molar mass of phase :math:`p` represented by component -:math:`c`), the following partition matrix establishes the component distribution within the two -phases: - -.. math:: - \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 :math:`x_c` is the mole fraction of component :math:`c` in the liquid phase and :math:`y_c` -is the mole fraction of component :math:`c` in the vapor phase. - -The fluid properties are updated through the following steps: - -1) The phase fractions (:math:`\nu_p`) and phase component fractions (:math:`x_c` and :math:`y_c`) are -computed as a function of pressure (:math:`p`), temperature (:math:`T`) and total component fractions -(:math:`z_c`). - -2) The phase densities (:math:`\rho_p`) and phase viscosities (:math:`\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 -:ref:`CompositionalMultiphaseFlow` 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 :math:`z_c` is stable -as a single phase at the current pressure :math:`p` and temperature :math:`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 :math:`z` is stable -at a specified pressure :math:`p` and temperature :math:`T` if and only if - -.. math:: - 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 :math:`y`, where :math:`\phi_i` denotes the fugacity -coefficient of component :math:`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 :math:`y_i = z_i/K_i` and :math:`y_i = z_iK_i` where :math:`K_i` are defined by - -.. math:: - K_i = \frac{P_{ci}}{p}\exp\left( 5.37( 1 + \omega_i ) \left(1-\frac{T_{ci}}{T}\right)\right) - -where :math:`P_{ci}` and :math:`T_{ci}` are respectively, the critical pressure and temperature of -component :math:`i` and :math:`\omega_i` is the accentric factor of component :math:`i`. - -The stability problem is solved by observing that a necessary condition is that :math:`g(y)` must -be non-negative at all its stationary points. The stationarity criterion can be expressed as - -.. math:: - \ln y_i + \ln \phi_i(y) - h_i = k \hspace{1cm} i=1,2,3,\ldots,N_c - -where :math:`h_i = \ln z_i + \ln \phi_i(z)` is a constant parameter dependent on the feed composition -:math:`z` and :math:`k` is an undetermined constant. This constant can be further incorporated into -the equation by defining the unnormalized trial phase moles :math:`Y_i` as - -.. math:: - Y_i = \exp(-k)y_i - -which reduces the stationarity criterion to - -.. math:: - \ln Y_i + \ln \phi_i(y) - h_i = 0 - -with the mole fractions :math:`y_i` related to the trial phase moles :math:`Y_i` by - -.. math:: - 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 :math:`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 :math:`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 - -.. math:: - T_{cp} = \frac{\sum_{i=1}^{N_c}T_{ci}V_{ci}z_{i}}{\sum_{i=1}^{N_c}V_{ci}z_{i}} - -where :math:`V_{ci}` and :math:`T_{ci}` are respectively the critical volume and temperature of component -:math:`i`. This is compared to the current temperature :math:`T` such that if :math:`T_{cp}`` node in the input. - -The following attributes are supported: - -.. include:: /docs/sphinx/datastructure/CompositionalTwoPhaseFluid.rst - -Example -========================= - -.. code-block:: xml - - - - - -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. - -.. _Petrowiki: https://petrowiki.spe.org/Phase_behavior_in_reservoir_simulation#Equation-of-state_models diff --git a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst index 75adc853688..742cbfe3224 100644 --- a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst +++ b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst @@ -9,6 +9,8 @@ Some text .. toctree:: :maxdepth: 2 + Overview + PhaseSplit Density From f82b1163b8082639b04762831877a45608a50901 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Wed, 3 Jun 2026 17:37:38 -0500 Subject: [PATCH 05/13] Density models --- .../CompositionalMultiphaseFluid.rst | 2 +- .../docs/compositional/Density.rst | 110 ++++++++++++++---- 2 files changed, 91 insertions(+), 21 deletions(-) diff --git a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst index 742cbfe3224..df91acdac8a 100644 --- a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst +++ b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst @@ -7,7 +7,7 @@ Compositional multiphase fluid model Some text .. toctree:: - :maxdepth: 2 + :maxdepth: 1 Overview diff --git a/src/coreComponents/constitutive/docs/compositional/Density.rst b/src/coreComponents/constitutive/docs/compositional/Density.rst index 47e8aa30b33..3185c190835 100644 --- a/src/coreComponents/constitutive/docs/compositional/Density.rst +++ b/src/coreComponents/constitutive/docs/compositional/Density.rst @@ -1,46 +1,116 @@ -======= Density ======= The accurate determination of phase density is critical for modelling fluid flow, as it directly impacts volumetric calculations, gravity drainage, and phase mobilities. -Compositional Density Model +Compositional density model --------------------------- -The compositional density model derives the molar volume of a phase using the compressibility factor, :math:`Z`, computed directly from the cubic Equation of State (e.g., Peng-Robinson). The molar density, :math:`\rho_{molar}`, is given by: +The compositional density model first calculates the uncorrected molar volume, :math:`v_{EOS}`, of a phase using the compressibility factor, :math:`Z`, obtained from the cubic Equation of State: + +.. math:: + v_{EOS} = \frac{Z R T}{P} + +To improve the accuracy of liquid phase densities, a Peneloux volume shift s described by Péneloux et al. (1982). correction is applied. The corrected molar volume, :math:`v`, is calculated by subtracting a composition-weighted sum of component-specific dimensional volume shifts, :math:`c_i`, from the EOS molar volume: .. math:: - \rho_{molar} = \frac{P}{Z R T} + v = v_{EOS} - \sum_{i=1}^{N_c} x_i c_i + +where :math:`x_i` is the mole fraction of component :math:`i` in the phase. -Because standard cubic equations of state notoriously under-predict liquid phase densities, a Peneloux volume shift correction is mathematically applied. The corrected molar volume, :math:`v_{corrected}`, is adjusted by a sum of component-specific volume shift parameters, :math:`c_i`, weighted by the phase composition: +The dimensional volume shifts, :math:`c_i`, are computed from the user-provided non-dimensional volume shift parameters, :math:`s_i` (specified as ``componentVolumeShift``), along with the critical temperature, :math:`T_{c,i}`, and critical pressure, :math:`P_{c,i}`, of each component: .. math:: - v_{corrected} = v_{EOS} - \sum_{i=1}^{N_c} x_i c_i + c_i = s_i \Omega_B R \frac{T_{c,i}}{P_{c,i}} + +The EOS-specific constant :math:`\Omega_B` takes the value :math:`0.077796074` for the Peng-Robinson equation of state and :math:`0.08664` for the Soave-Redlich-Kwong equation of state. -The molar density is then inverted to mass density using the molecular weight of the phase mixture. This model is recommended for all hydrocarbon phases. +The molar density of the phase, :math:`\rho_{molar}`, is the inverse of the corrected molar volume: + +.. math:: + \rho_{molar} = \frac{1}{v} + +Finally, the mass density of the phase, :math:`\rho_{mass}`, is obtained by multiplying the molar density by the apparent molecular weight of the mixture: + +.. math:: + \rho_{mass} = \rho_{molar} \sum_{i=1}^{N_c} x_i MW_i -* **Catalog Name:** ``CompositionalDensity`` -* **Parameters:** * ``componentVolumeShift``: Component-specific volume shift parameters. +where :math:`MW_i` is the molecular weight of component :math:`i`. This model is recommended for all hydrocarbon phases. -Phillips Brine Density Model +Parameters: + +* ``componentVolumeShift``: Component-specific volume shift parameters (dimensionless). + +Phillips brine density model ---------------------------- -For aqueous phases containing dissolved salts, cubic equations of state often struggle to capture the complex ionic interactions. The Phillips Brine Density model employs a robust empirical correlation developed by Phillips et al. (1981). +For aqueous phases containing dissolved salts, cubic equations of state often struggle to capture complex ionic interactions. To address this, the Phillips brine density model utilizes the empirical correlation developed by Phillips et al. (1981) to pre-compute a 2D lookup table of volume shifts specifically for the water component as a function of pressure and temperature. + +First, the apparent molar weight of the brine, :math:`MW_{brine}`, is calculated from the salinity, :math:`S` (molality), the salt molar weight, :math:`MW_{salt}`, and the water molar weight, :math:`MW_{H2O}`: + +.. math:: + MW_{brine} = \frac{1}{S + \frac{1 - S \cdot MW_{salt}}{MW_{H2O}}} + +The target mass density of the brine, :math:`\rho_{brine}`, is then evaluated. The pressure must be in Pascal and must be less than :math:`5 \times 10^7` Pascal. The temperature must be in Kelvin and must be between :math:`283.15` and :math:`623.15` Kelvin. Using these parameters, a preprocessing step constructs a two-dimensional table storing the brine density for the specified salinity as a function of pressure and temperature. For salinities above a threshold of :math:`0.25` mol/kg, the density is calculated directly using the Phillips correlation, :math:`\rho_{Phillips}(P, T, S)`, using the expression: -This model calculates the mass density of the brine as a direct polynomial and exponential function of pressure, temperature, and salinity (molality of dissolved salts). For low salinities, the model interpolates smoothly towards pure water properties to maintain physical continuity. This is highly recommended for geothermal reservoirs or saline aquifers. +.. math:: + \rho_{Phillips}(P, T, S) = A + Bx + Cx^2 + Dx^3 -* **Catalog Name:** ``PhillipsBrineDensity`` -* **Parameters:** - * ``salinity``: Brine salinity. - * ``saltMolarWeight``: The molar weight of the salt component. +.. math:: + x = c_1 \exp(a_1 S) + c_2 \exp(a_2 T) + c_3 \exp(a_3 P) -Immiscible Water Density +We refer the reader to Phillips et al. (1981), equations (4) and (5), pages 14 and 15 for the definition of the coefficients involved in the previous equation. For low salinities (:math:`S < 0.25`), the model linearly interpolates between the pure water density, :math:`\rho_{pure}(P, T)`, and the Phillips correlation to maintain physical continuity towards the pure water limit: + +.. math:: + \rho_{brine}(P, T) = \begin{cases} + \rho_{Phillips}(P, T, S) & \text{if } S \ge 0.25 \\ + \frac{S}{0.25} \rho_{Phillips}(P, T, S) + \left(1 - \frac{S}{0.25}\right) \rho_{pure}(P, T) & \text{if } S < 0.25 + \end{cases} + +The 2D table of dimensional volume shifts for water, :math:`c_{H2O}(P, T)`, is generated by taking the difference between the target molar volume of the brine and the molar volume of pure water predicted by the chosen Equation of State, :math:`v_{EOS,H2O}`: + +.. math:: + c_{H2O}(P, T) = \frac{MW_{brine}}{\rho_{brine}(P, T)} - v_{EOS,H2O}(P, T) + +During the simulation, this pre-computed volume shift couples the empirical correlation to the overall phase EOS density calculation. The corrected molar volume of the aqueous phase, :math:`v`, is obtained by adding the water component's volume shift weighted by its mole fraction, :math:`x_{H2O}`, to the uncorrected phase EOS molar volume, :math:`v_{EOS}`: + +.. math:: + v = v_{EOS} + x_{H2O} c_{H2O}(P, T) + +The phase molar density, :math:`\rho_{molar}`, is calculated as the inverse of the corrected molar volume: + +.. math:: + \rho_{molar} = \frac{1}{v} + +The phase mass density, :math:`\rho_{mass}`, is then computed by multiplying the phase molar density by the apparent molar weight of the phase mixture: + +.. math:: + \rho_{mass} = \rho_{molar} \left( x_{H2O} MW_{brine} + \sum_{i \neq H2O} x_i MW_i \right) + +*(Note: While the volume shift tightly couples the Phillips correlation to the EOS density, the brine viscosity is not derived from this density. Instead, the phase viscosity is modelled independently using a distinct temperature- and salinity-dependent multiplier applied to the pure water viscosity).* + +Parameters: + +* ``salinity``: Brine salinity (mol/kg). +* ``saltMolarWeight``: The molar weight of the salt component (kg/mol). + +Immiscible water density ------------------------ -A simplified exponential density model is available for pure, immiscible water phases. It evaluates density based on isothermal compressibility and thermal expansion from a reference state: +A simplified exponential density model is available for pure, immiscible water phases. It evaluates the mass density, :math:`\rho`, based on isothermal compressibility and thermal expansion from a reference state: .. math:: \rho = \rho_{ref} \exp(c_w (P - P_{ref})) \exp(-\alpha_w (T - T_{ref})) -* **Catalog Name:** ``ImmiscibleWaterDensity`` -* **Parameters:** ``waterReferencePressure``, ``waterReferenceTemperature``, ``waterDensity``, ``waterCompressibility``, ``waterExpansionCoefficient``. \ No newline at end of file +The phase molar density, :math:`\rho_{molar}`, is then calculated by dividing the mass density by the molecular weight of water, :math:`MW_{H2O}`: + +.. math:: + \rho_{molar} = \frac{\rho}{MW_{H2O}} + +Parameters: + +* ``waterReferencePressure``: Reference pressure for the water density calculation (Pa). +* ``waterReferenceTemperature``: Reference temperature for the water density calculation (K). +* ``waterDensity``: Water mass density at the reference pressure and temperature (kg/m^3). +* ``waterCompressibility``: Isothermal compressibility of water (1/Pa). +* ``waterExpansionCoefficient``: Volumetric coefficient of thermal expansion of water (1/K). \ No newline at end of file From 781c44ce305510e9f26e1b69a90c46fc72a2a4cd Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 11:18:39 -0500 Subject: [PATCH 06/13] Update density text --- scripts/SchemaToRSTDocumentation.py | 2 +- .../constitutive/docs/FluidModels.rst | 2 +- .../CompositionalMultiphaseFluid.rst | 25 --------------- .../docs/compositional/Density.rst | 11 ++++--- .../docs/compositional/Enthalpy.rst | 1 - .../docs/compositional/EquationsOfState.rst | 18 +++++------ .../compositional/ImmiscibleWaterFlash.rst | 29 ++++++++--------- .../docs/compositional/KValueFlash.rst | 23 +++++++------- .../docs/compositional/Miscellaneous.rst | 1 - .../docs/compositional/NegativeFlash.rst | 31 ++++++++++--------- .../docs/compositional/Overview.rst | 12 ------- .../docs/compositional/PhaseSplit.rst | 12 ++----- .../docs/compositional/References.rst | 1 - .../docs/compositional/Viscosity.rst | 1 - 14 files changed, 63 insertions(+), 106 deletions(-) delete mode 100644 src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst diff --git a/scripts/SchemaToRSTDocumentation.py b/scripts/SchemaToRSTDocumentation.py index 8a7b1a8cda7..ee4330b1eb6 100644 --- a/scripts/SchemaToRSTDocumentation.py +++ b/scripts/SchemaToRSTDocumentation.py @@ -56,7 +56,7 @@ def writeTableRST(type_name, file_name, title_prefix, values): formatted_lines.append(multiple_row_format_b.format(d)) # Build table - with open(file_name, 'w') as f: + with open(file_name, 'w', encoding="utf-8") as f: f.write('%s\n' % (element_header)) f.write('=' * len(element_header) + '\n') f.write('\n') diff --git a/src/coreComponents/constitutive/docs/FluidModels.rst b/src/coreComponents/constitutive/docs/FluidModels.rst index 4c5ad023db3..1aed1ce387e 100644 --- a/src/coreComponents/constitutive/docs/FluidModels.rst +++ b/src/coreComponents/constitutive/docs/FluidModels.rst @@ -15,6 +15,6 @@ single fluids and fluid mixtures. BlackOilFluid - compositional/CompositionalMultiphaseFluid + CompositionalMultiphaseFluid CO2BrineFluid diff --git a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst deleted file mode 100644 index df91acdac8a..00000000000 --- a/src/coreComponents/constitutive/docs/compositional/CompositionalMultiphaseFluid.rst +++ /dev/null @@ -1,25 +0,0 @@ -.. _CompositionalFluidModel: - -############################################ -Compositional multiphase fluid model -############################################ - -Some text - -.. toctree:: - :maxdepth: 1 - - Overview - - PhaseSplit - - Density - - Viscosity - - Enthalpy - - Miscellaneous - - References - diff --git a/src/coreComponents/constitutive/docs/compositional/Density.rst b/src/coreComponents/constitutive/docs/compositional/Density.rst index 3185c190835..d63bee00b86 100644 --- a/src/coreComponents/constitutive/docs/compositional/Density.rst +++ b/src/coreComponents/constitutive/docs/compositional/Density.rst @@ -37,7 +37,8 @@ Finally, the mass density of the phase, :math:`\rho_{mass}`, is obtained by mult where :math:`MW_i` is the molecular weight of component :math:`i`. This model is recommended for all hydrocarbon phases. -Parameters: +Parameters +~~~~~~~~~~~~~~~~ * ``componentVolumeShift``: Component-specific volume shift parameters (dimensionless). @@ -49,7 +50,7 @@ For aqueous phases containing dissolved salts, cubic equations of state often st First, the apparent molar weight of the brine, :math:`MW_{brine}`, is calculated from the salinity, :math:`S` (molality), the salt molar weight, :math:`MW_{salt}`, and the water molar weight, :math:`MW_{H2O}`: .. math:: - MW_{brine} = \frac{1}{S + \frac{1 - S \cdot MW_{salt}}{MW_{H2O}}} + MW_{brine} = \frac{MW_{H2O}}{1 + S\left(MW_{H2O} - MW_{salt}\right)} The target mass density of the brine, :math:`\rho_{brine}`, is then evaluated. The pressure must be in Pascal and must be less than :math:`5 \times 10^7` Pascal. The temperature must be in Kelvin and must be between :math:`283.15` and :math:`623.15` Kelvin. Using these parameters, a preprocessing step constructs a two-dimensional table storing the brine density for the specified salinity as a function of pressure and temperature. For salinities above a threshold of :math:`0.25` mol/kg, the density is calculated directly using the Phillips correlation, :math:`\rho_{Phillips}(P, T, S)`, using the expression: @@ -89,7 +90,8 @@ The phase mass density, :math:`\rho_{mass}`, is then computed by multiplying the *(Note: While the volume shift tightly couples the Phillips correlation to the EOS density, the brine viscosity is not derived from this density. Instead, the phase viscosity is modelled independently using a distinct temperature- and salinity-dependent multiplier applied to the pure water viscosity).* -Parameters: +Parameters +~~~~~~~~~~~~~~~~ * ``salinity``: Brine salinity (mol/kg). * ``saltMolarWeight``: The molar weight of the salt component (kg/mol). @@ -107,7 +109,8 @@ The phase molar density, :math:`\rho_{molar}`, is then calculated by dividing th .. math:: \rho_{molar} = \frac{\rho}{MW_{H2O}} -Parameters: +Parameters +~~~~~~~~~~~~~~~~ * ``waterReferencePressure``: Reference pressure for the water density calculation (Pa). * ``waterReferenceTemperature``: Reference temperature for the water density calculation (K). diff --git a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst index 7be1e969ef7..94a00d7d0a7 100644 --- a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst +++ b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst @@ -1,4 +1,3 @@ -======== Enthalpy ======== diff --git a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst index cc701567f07..9f213d2fe94 100644 --- a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst +++ b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst @@ -1,5 +1,5 @@ Equations of state -================== +------------------ An Equation of State (EoS) mathematically relates the pressure, volume, temperature, and composition of a fluid. In compositional modelling, the EoS serves a dual purpose: 1. It calculates the fugacity coefficients, :math:`\phi_i`, required to establish thermodynamic phase equilibrium. @@ -8,9 +8,9 @@ An Equation of State (EoS) mathematically relates the pressure, volume, temperat The fugacity coefficient of a component is derived from the exact thermodynamic relationship involving the integration of the EoS volume departure with respect to pressure. Cubic equation of state ------------------------ +~~~~~~~~~~~~~~~~~~~~~~~ -The Peng-Robinson (PR) and Soave-Redlich-Kwong (SRK) equations of state are widely used for hydrocarbon systems encompassing natural gas, gas condensates, and volatile oils. The fundamental pressure-volume-temperature (PVT) relationships for these equations are defined as: +The Peng-Robinson (PR) (Peng and Robinson, 1976; Robinson and Peng, 1978) and Soave-Redlich-Kwong (SRK) (Soave, 1972) equations of state are widely used for hydrocarbon systems encompassing natural gas, gas condensates, and volatile oils. The fundamental pressure-volume-temperature (PVT) relationships for these equations are defined as: * Peng-Robinson (PR): @@ -106,12 +106,12 @@ where the intermediate terms are defined as: .. math:: E = \ln \left( \frac{Z + \delta_1 B}{Z + \delta_2 B} \right), \quad F = \ln(Z - B), \quad G = \frac{1}{(\delta_1 - \delta_2)B} -Soreide-whitson equation of state ---------------------------------- +Soreide-Whitson equation of state +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The Soreide-Whitson model represents a specialized modification of the Peng-Robinson EoS tailored specifically for aqueous-hydrocarbon mixtures, such as those encountered in carbon sequestration (CO2-brine systems) or water-flooding in oil reservoirs. It adapts the standard cubic EoS by capturing the "salting-out" effect, where increased brine salinity reduces the solubility of hydrocarbons and CO2 in the aqueous phase. It is highly recommended whenever accurate modelling of gas dissolution in saline water is critical. -This model introduces two primary modifications to the standard Peng-Robinson formulation, applying special treatment exclusively to the water component and its interactions within the aqueous phase: +This model introduces two primary modifications to the standard Peng-Robinson formulation, applying special treatment exclusively to the water component and its interactions within the aqueous phase and follows closely the description given by Soreide and Whitson (1992): First, the standard alpha function, :math:`\alpha_i`, is replaced specifically for the water component by a salinity-dependent formulation. If :math:`c_{sw}` represents the brine salinity, the water alpha function is evaluated as: @@ -170,7 +170,7 @@ For hydrogen sulfide (H2S): where :math:`A_0 = -0.20441` and :math:`A_1 = 0.23426`. -For hydrogen (H2): +For hydrogen (H2) the correlation due to Chabab et al. (2024) is used: .. math:: k_{ij}^{AQ} = D_0 (1 + a_0 c_{sw}) + D_1 T_{r,i} (1 + a_1 c_{sw}) + D_2 \exp(D_3 T_{r,i}) @@ -178,7 +178,7 @@ For hydrogen (H2): where :math:`D_0 = -2.11917`, :math:`D_1 = 0.14888`, :math:`D_2 = -13.01835`, :math:`D_3 = -0.43946`, :math:`a_0 = -0.0226322`, and :math:`a_1 = -0.0044736`. Component naming requirements -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Because the Soreide-Whitson EoS applies highly specific correlations based on the exact type of component interacting with water, it relies on the component's assigned name to identify its type. This name matching is case-insensitive. @@ -193,7 +193,7 @@ To trigger the correct dynamic binary interaction coefficients, you must use the If a component's name does not match any of these predefined strings (for example, if you name it "CarbonDioxide", "C1", or "Methane"), the model will automatically default to categorising it as a generic hydrocarbon (``hc``). While this is the intended and correct behavior for actual hydrocarbons (where methane is safely identified via its acentric factor), using an unrecognized name for a special non-hydrocarbon (like CO2) will cause the simulator to mistakenly apply the generic hydrocarbon-water correlation instead of its highly tuned specific correlation. Model parameters ----------------- +~~~~~~~~~~~~~~~~ Equations of state are assigned per-phase within the specific fluid model XML block (e.g., ````). diff --git a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst index 1e1ebc0e4e6..2134908ade0 100644 --- a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst @@ -1,10 +1,10 @@ Immiscible water flash -====================== +---------------------- The immiscible water flash is an efficient three-phase equilibrium model designed for systems containing an aqueous phase alongside hydrocarbon liquid and vapour phases. Instead of solving a fully coupled three-phase equation of state problem, this model vastly simplifies the phase split by assuming water forms a strictly pure aqueous phase and is completely immiscible in the hydrocarbon phases, while hydrocarbons are entirely insoluble in the aqueous phase. Methodology ------------ +~~~~~~~~~~~ In this formulation, the mixture composition is explicitly separated into a water component and a hydrocarbon mixture before any thermodynamic equilibrium calculations occur. The water component is identified using its assigned component name (e.g., ``H2O``). @@ -26,7 +26,7 @@ The feed composition of the remaining hydrocarbon mixture is then normalised: for all non-water components, while the normalised hydrocarbon feed fraction for water is set to zero. Equilibrium calculation ------------------------ +~~~~~~~~~~~~~~~~~~~~~~~ During the simulation step, if the total hydrocarbon fraction :math:`z_{hc}` is negligibly small, the model bypasses the flash and assigns the system as a single-phase aqueous fluid. @@ -45,14 +45,14 @@ Once the two-phase hydrocarbon flash converges, the overall phase fractions for The derivatives of the phase fractions and compositions with respect to pressure, temperature, and total composition are analytically transformed using the chain rule to account for this scaling and normalisation. Recommendations and pitfalls ----------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The immiscible water flash is highly computationally efficient and is recommended for multi-phase reservoir simulations where water solubility in the hydrocarbon phases (and vice versa) is negligible. By avoiding a full three-phase fugacity solver, it drastically reduces simulation time and improves non-linear solver stability. Pitfalls: Because this model strictly enforces total immiscibility, it cannot accurately model scenarios where mutual solubility is physically significant, such as high-pressure CO2 injection into saline aquifers (where CO2 dissolves heavily into brine, and water vaporises into the gas stream). For such complex mutual solubility problems, using a two-phase flash paired with the Soreide-Whitson equation of state might be required. Model parameters ----------------- +~~~~~~~~~~~~~~~~ The immiscible water flash fluid models are assigned using specific XML blocks that pair the three-phase flash model with distinct viscosity and density models for the hydrocarbon and aqueous phases. @@ -75,12 +75,13 @@ Here is an example demonstrating how to specify the parameters for a three-phase waterViscosity="1.002e-3" waterCompressibility="4.5e-10" /> -* catalog options: ``CompositionalThreePhaseLohrenzBrayClarkViscosity`` -* parameters: - * ``equationsOfState``: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment. - * ``waterReferencePressure``: The reference pressure for water density and viscosity [Pa]. - * ``waterReferenceTemperature``: The reference temperature for water density and viscosity [K] (optional, default: 293.15). - * ``waterDensity``: The water density at the reference pressure and temperature [kg/m^3]. - * ``waterViscosity``: The water viscosity at the reference pressure and temperature [Pa.s]. - * ``waterCompressibility``: The constant isothermal compressibility of water [1/Pa]. - * ``waterExpansionCoefficient``: The volumetric coefficient of thermal expansion of water [1/K] (optional, default: 0.0). \ No newline at end of file +Parameters +~~~~~~~~~~~~~~~~ + +* ``equationsOfState``: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment. +* ``waterReferencePressure``: The reference pressure for water density and viscosity [Pa]. +* ``waterReferenceTemperature``: The reference temperature for water density and viscosity [K] (optional, default: 293.15). +* ``waterDensity``: The water density at the reference pressure and temperature [kg/m^3]. +* ``waterViscosity``: The water viscosity at the reference pressure and temperature [Pa.s]. +* ``waterCompressibility``: The constant isothermal compressibility of water [1/Pa]. +* ``waterExpansionCoefficient``: The volumetric coefficient of thermal expansion of water [1/K] (optional, default: 0.0). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst index 291b084beb6..fa2763d2126 100644 --- a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst @@ -1,12 +1,12 @@ K-value flash -============= +------------- The K-value flash model offers a computationally simplified alternative to the rigorous equation of state flash calculation. Rather than iteratively evaluating complex equations of state to find fugacity coefficients, this model determines phase splits using pre-tabulated K-values. The model supports both two-phase and three-phase variants. The two-phase flash calculates equilibrium partitioning between a liquid and a vapour phase. The three-phase flash introduces an additional aqueous phase, allowing for the predefined partitioning of components into water alongside the hydrocarbon phases. Methodology ------------ +~~~~~~~~~~~ In this formulation, the equilibrium ratio for each component in the secondary phases relative to the primary liquid phase, :math:`K_i = y_i / x_i`, is treated strictly as a predefined function of pressure and temperature. @@ -15,11 +15,11 @@ Users input these functions as either tables or symbolic functions. To optimize The pressure and temperature coordinates for this hypercube grid can be explicitly chosen and specified by the user. If they are omitted, the model automatically inspects the provided table functions, extracts all unique coordinate points from the underlying data, and constructs the grid dynamically. If only a single point is found, an artificial offset is added to permit interpolation. Equilibrium calculation ------------------------ +~~~~~~~~~~~~~~~~~~~~~~~ During the simulation step, the model performs bilinear interpolation on the pre-built hypercube to evaluate the K-values and their partial derivatives at the current pressure and temperature. -Because the K-values are explicitly known, the iterative fugacity updates are bypassed. Instead, the vapour fraction, :math:`V`, is found by directly solving the Rachford-Rice equation using successive substitution and Newton iterations: +Because the K-values are explicitly known, the iterative fugacity updates are bypassed. Instead, the vapour fraction, :math:`V`, is found by directly solving the Rachford-Rice equation (Rachford and Rice, 1952) using successive substitution and Newton iterations: .. math:: \sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V(K_i - 1)} = 0 @@ -32,14 +32,14 @@ Once the vapour fraction is determined, the phase compositions are directly eval This direct solution approach intrinsically supports a "negative" flash. The mathematical solver allows the vapour fraction, :math:`V`, to temporarily converge to values outside the physically meaningful bounds of :math:`[0, 1]`. If the resulting :math:`V \le 0` or :math:`V \ge 1`, the system cleanly truncates the compositions to represent a single-phase liquid or single-phase vapour solution, respectively. Recommendations and pitfalls ----------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The K-value flash is highly computationally efficient and is recommended for thermal recovery processes (like steam injection) or relatively stable black-oil systems where fluid compositions remain within a narrow, predictable envelope. Pitfalls: In reality, thermodynamic K-values are heavily dependent on the overall mixture composition, not just pressure and temperature. Relying on tabulated K-values implies assuming a fixed compositional path. Care should be taken when using this model for highly variable compositional processes such as miscible gas injection (e.g., CO2 or hydrocarbon gas flooding), or strong depletion of gas condensates, where composition dependency is critical to accurate phase behaviour. Model parameters ----------------- +~~~~~~~~~~~~~~~~ The K-value flash fluid models are assigned using specific XML blocks that pair the flash model with viscosity and density models. @@ -67,8 +67,9 @@ Here is an example demonstrating how to specify the parameters for a two-phase s waterCompressibility="4.1483E-10" salinity="2.3" /> -* catalog options: ``CompositionalKValueLohrenzBrayClarkViscosity``, ``CompositionalKValuePhillipsBrine`` -* parameters: - * ``kValueTables``: List of function names linking to the pre-tabulated K-values. For an :math:`N`-phase system, this requires :math:`(N-1) \times N_c` function names. - * ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). - * ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). \ No newline at end of file +Parameters +~~~~~~~~~~~~~~~~ + +* ``kValueTables``: List of function names linking to the pre-tabulated K-values. For an :math:`N`-phase system, this requires :math:`(N-1) \times N_c` function names. +* ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). +* ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst index b82c7442a8f..dad45161a65 100644 --- a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst +++ b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst @@ -1,4 +1,3 @@ -============= Miscellaneous ============= diff --git a/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst index 2068968eb75..ef84be59441 100644 --- a/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst @@ -1,19 +1,19 @@ -Negative flash -============== +Negative two-phase flash +------------------------ The negative flash is a robust, rigorous approach to solving the phase equilibrium problem using an Equation of State (EoS). Stability analysis ------------------- +~~~~~~~~~~~~~~~~~~ -Before executing a flash calculation, the model must ascertain if the homogeneous single-phase mixture is thermodynamically unstable at the specified pressure :math:`P` and temperature :math:`T`. This is performed via a tangent plane distance (TPD) analysis developed by Michelsen. The TPD assesses whether the Gibbs free energy of the mixture can be minimised by splitting into multiple phases. A phase with feed composition :math:`z` is stable if and only if: +Before executing a flash calculation, the model must ascertain if the homogeneous single-phase mixture is thermodynamically unstable at the specified pressure :math:`P` and temperature :math:`T`. This is performed via a tangent plane distance (TPD) analysis developed by Michelsen (Michelsen, 1982a; Michelsen, 1982b). The TPD assesses whether the Gibbs free energy of the mixture can be minimised by splitting into multiple phases. A phase with feed composition :math:`z` is stable if and only if: .. math:: g(\mathbf{y}) = \sum_{i=1}^{N_c} y_i \left( \ln y_i + \ln \phi_i(\mathbf{y}) - \ln z_i - \ln \phi_i(\mathbf{z}) \right) \ge 0 for any permissible trial composition :math:`\mathbf{y}`, where :math:`\phi_i` denotes the fugacity coefficient of component :math:`i`. -To determine stability, testing is initiated from basic starting points using Wilson's K-values to generate a lighter trial mixture (:math:`y_i = z_i / K_i`) and a heavier trial mixture (:math:`y_i = z_i K_i`). Wilson's K-values are defined as: +To determine stability, testing is initiated from basic starting points using Wilson's K-values (Wilson, 1969; Whitson and Torp, 1981) to generate a lighter trial mixture (:math:`y_i = z_i / K_i`) and a heavier trial mixture (:math:`y_i = z_i K_i`). Wilson's K-values are defined as: .. math:: K_i = \frac{P_{ci}}{P} \exp \left( 5.37 (1 + \omega_i) \left( 1 - \frac{T_{ci}}{T} \right) \right) @@ -23,7 +23,7 @@ where :math:`P_{ci}`, :math:`T_{ci}`, and :math:`\omega_i` are the critical pres The stationarity condition, :math:`\ln Y_i + \ln \phi_i(\mathbf{y}) - h_i = 0` (where :math:`h_i = \ln z_i + \ln \phi_i(\mathbf{z})` and :math:`Y_i` are unnormalized trial phase moles), is solved using successive substitution. If both initial states converge to a solution where :math:`g(\mathbf{y}) \ge 0`, the mixture is stable; otherwise, it is unstable. Phase labeling --------------- +~~~~~~~~~~~~~~ If the stability test confirms a single stable phase, the two-phase flash calculation is bypassed. To appropriately label the single phase as liquid or vapour for relative permeability evaluations, the model utilises a "Li-temperature" correlation. This pseudo-critical temperature is weighted by the critical volumes (:math:`V_{ci}`): @@ -33,7 +33,7 @@ If the stability test confirms a single stable phase, the two-phase flash calcul If :math:`T_{cp} < T`, the mixture is logically categorised as a vapour; otherwise, it is labelled as a liquid. Flash calculation ------------------ +~~~~~~~~~~~~~~~~~ If the mixture is deemed unstable, the system must determine the phase split by ensuring thermodynamic equilibrium, meaning the fugacities in the liquid and vapour phases must be equal (:math:`\phi_{iL} = \phi_{iV}`). @@ -62,14 +62,15 @@ The flash algorithm: The term "negative" flash arises because the algorithm temporarily permits the vapour fraction, :math:`V`, to converge to values slightly outside the physically meaningful bounds of :math:`[0, 1]`. This mathematical relaxation prevents the solver from getting artificially trapped at phase boundaries, vastly improving convergence robustness near the critical point or saturation envelopes. Upon convergence, if :math:`V` is negative or greater than unity, the system truncates to a single-phase solution. Recommendations ---------------- +~~~~~~~~~~~~~~~ The negative two-phase flash is the recommended standard for rigorous compositional modelling of miscible gas injection, primary depletion of volatile oils, and scenarios where fluid compositions change drastically over time. An extended three-phase variant is available for immiscible water systems. -* Catalog name: ``TwoPhase``, ``ThreePhase`` -* Parameters: - * ``stabilityThreshold``: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). - * ``stabilityTolerance``: Tolerance for stationarity in the stability test. - * ``stabilityMaxIterations``: Maximum successive substitution steps for stability analysis. - * ``flashTolerance``: Convergence tolerance for the fugacity ratio error. - * ``flashMaxIterations``: Maximum successive substitution steps for the flash solve. \ No newline at end of file +Parameters +~~~~~~~~~~~~~~~~ + +* ``stabilityThreshold``: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). +* ``stabilityTolerance``: Tolerance for stationarity in the stability test. +* ``stabilityMaxIterations``: Maximum successive substitution steps for stability analysis. +* ``flashTolerance``: Convergence tolerance for the fugacity ratio error. +* ``flashMaxIterations``: Maximum successive substitution steps for the flash solve. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Overview.rst b/src/coreComponents/constitutive/docs/compositional/Overview.rst index 7e99cc8c6a0..7ef4327ec21 100644 --- a/src/coreComponents/constitutive/docs/compositional/Overview.rst +++ b/src/coreComponents/constitutive/docs/compositional/Overview.rst @@ -43,15 +43,3 @@ The evaluation of fluid properties follows a sequential thermodynamic structure * Thermal Properties: If the simulation is non-isothermal, the phase enthalpies and internal energies are computed using ideal gas heat capacities and EoS departure functions. 4. Total Fluid Properties: The individual phase properties are homogenised to yield total fluid properties (e.g., total density). Derivatives with respect to pressure, temperature, and composition are rigorously propagated using the chain rule. - -Structure of this Documentation -------------------------------- - -This documentation is divided into the following sections: - -* `Phase Split `_: Details the thermodynamic stability and flash calculations (Negative Flash and K-Value Flash), including the underlying Equations of State. -* `Density `_: Describes the models used to compute molar and mass phase densities. -* `Viscosity `_: Details the models used to evaluate the resistance to flow for each phase. -* `Enthalpy `_: Explains the calculation of thermal properties. -* `Miscellaneous `_: Outlines internal unit conversions and mass-to-molar variable handling. -* `References `_: Provides the bibliography for the scientific correlations implemented in the code. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst b/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst index 54d970ccdb8..0020a6f642d 100644 --- a/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst +++ b/src/coreComponents/constitutive/docs/compositional/PhaseSplit.rst @@ -8,9 +8,9 @@ Phase split equations The phase split relies on two fundamental physical principles: the conservation of mass (material balance) and thermodynamic equilibrium. -**Material Balance (The Rachford-Rice Equation)** +**Material Balance** -For a two-phase system (liquid and vapour), the overall mole fraction of a component, :math:`z_i`, must equal the sum of its molar contributions in the liquid phase, :math:`x_i`, and the vapour phase, :math:`y_i`. Introducing the vapour phase fraction, :math:`V`, and the equilibrium ratio (K-value), :math:`K_i = y_i / x_i`, the material balance is expressed by the Rachford-Rice equation: +For a two-phase system (liquid and vapour), the overall mole fraction of a component, :math:`z_i`, must equal the sum of its molar contributions in the liquid phase, :math:`x_i`, and the vapour phase, :math:`y_i`. Introducing the vapour phase fraction, :math:`V`, and the equilibrium ratio (K-value), :math:`K_i = y_i / x_i`, the material balance is expressed by the Rachford-Rice equation (Rachford and Rice, 1952): .. math:: \sum_{i=1}^{N_c} \frac{z_i (K_i - 1)}{1 + V (K_i - 1)} = 0 @@ -37,11 +37,3 @@ Using the fugacity coefficient, :math:`\phi_i`, where :math:`f_i = x_i \phi_i P` The phase split models in this framework either rigorously solve for these fugacity coefficients using an Equation of State or use simplified, pre-tabulated K-values. -.. include:: EquationsOfState.rst - -.. include:: NegativeFlash.rst - -.. include:: KValueFlash.rst - -.. include:: ImmiscibleWaterFlash.rst - diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst index 0dc072d143e..5fbcfd263f6 100644 --- a/src/coreComponents/constitutive/docs/compositional/References.rst +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -1,4 +1,3 @@ -========== References ========== diff --git a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst index b2f32851d7d..786b26b32f8 100644 --- a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst +++ b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst @@ -1,4 +1,3 @@ -========= Viscosity ========= From f6c793ba305248260cc91dfb95c296a2eaf22781 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 11:21:29 -0500 Subject: [PATCH 07/13] Add main file --- .../docs/CompositionalMultiphaseFluid.rst | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst diff --git a/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst new file mode 100644 index 00000000000..8f4a8a9f2f3 --- /dev/null +++ b/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst @@ -0,0 +1,39 @@ +.. _CompositionalMultiphaseFluid: + +############################################ +Compositional multiphase fluid model +############################################ + +.. include:: compositional/Overview.rst + +Structure of this Documentation +------------------------------- + +This documentation is divided into the following sections: + +* `Phase Split`_: Details the thermodynamic stability and flash calculations (Negative Flash and K-Value Flash), including the underlying Equations of State. +* `Density`_: Describes the models used to compute molar and mass phase densities. +* `Viscosity`_: Details the models used to evaluate the resistance to flow for each phase. +* `Enthalpy`_: Explains the calculation of thermal properties. +* `Miscellaneous`_: Outlines internal unit conversions and mass-to-molar variable handling. +* `References`_: Provides the bibliography for the scientific correlations implemented in the code. + +.. include:: compositional/PhaseSplit.rst + +.. include:: compositional/EquationsOfState.rst + +.. include:: compositional/NegativeFlash.rst + +.. include:: compositional/KValueFlash.rst + +.. include:: compositional/ImmiscibleWaterFlash.rst + +.. include:: compositional/Density.rst + +.. include:: compositional/Viscosity.rst + +.. include:: compositional/Enthalpy.rst + +.. include:: compositional/Miscellaneous.rst + +.. include:: compositional/References.rst From b0840d629f6ec61ef68650f402d5f35e1050c6b4 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 12:40:19 -0500 Subject: [PATCH 08/13] Viscosity and enthalpy --- .../docs/compositional/Enthalpy.rst | 72 +++++++-- .../docs/compositional/Miscellaneous.rst | 10 +- .../docs/compositional/Viscosity.rst | 146 ++++++++++++++++-- 3 files changed, 192 insertions(+), 36 deletions(-) diff --git a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst index 94a00d7d0a7..7eccd18927f 100644 --- a/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst +++ b/src/coreComponents/constitutive/docs/compositional/Enthalpy.rst @@ -3,31 +3,73 @@ Enthalpy For non-isothermal simulations, such as thermal recovery processes or geothermal operations, the accurate calculation of phase enthalpies and internal energies is required to resolve the energy balance equation. -Compositional Enthalpy Model +Compositional enthalpy model ---------------------------- -The enthalpy of a compositional fluid phase is computed using a rigorous thermodynamic departure function approach. The true enthalpy at a given pressure and temperature, :math:`H(P, T)`, is expressed as the sum of an ideal gas enthalpy at the same temperature, :math:`H_{ideal}(T)`, and an EoS-derived departure enthalpy, :math:`H_{departure}(P, T, x)`: +The compositional enthalpy model for a multiphase fluid system computes the total enthalpy of the phase as the sum of an ideal (polynomial) enthalpy and a real fluid correction (or departure enthalpy), which is derived from the equation of state. + +The total enthalpy of a fluid phase is defined as: + +.. math:: + H_{total} = H_{ideal} + H_{dep} + +Ideal gas contribution +~~~~~~~~~~~~~~~~~~~~~~ + +The ideal enthalpy represents the enthalpy of the mixture assuming it behaves as an ideal gas, relying on the heat capacity of the individual components. The specific heat capacity :math:`C_{p,i}` for each component :math:`i` is modeled as a 4th-order Poling polynomial (Poling et al., 2001) based on the temperature difference from a reference state, :math:`\Delta T = T - T_{ref}`: + +.. math:: + C_{p,i}(T) = a_{i,0} + a_{i,1}\Delta T + a_{i,2}\Delta T^2 + a_{i,3}\Delta T^3 + a_{i,4}\Delta T^4 + +To find the ideal enthalpy for a component, this heat capacity is integrated with respect to temperature from the reference temperature to the current temperature, and then added to a baseline reference enthalpy :math:`H_{ref,i}`: + +.. math:: + H_{ideal,i}(T) = H_{ref,i} + \Delta T \left( a_{i,0} + \frac{a_{i,1}}{2}\Delta T + \frac{a_{i,2}}{3}\Delta T^2 + \frac{a_{i,3}}{4}\Delta T^3 + \frac{a_{i,4}}{5}\Delta T^4 \right) + +The total ideal enthalpy for the phase is simply the mole-fraction weighted sum of the individual component ideal enthalpies: .. math:: - H(P, T) = H_{ideal}(T) + H_{departure}(P, T, x) + H_{ideal} = \sum_i x_i H_{ideal,i}(T) -**Ideal Gas Contribution** +Departure function +~~~~~~~~~~~~~~~~~~ -The ideal gas enthalpy is derived by integrating a temperature-dependent specific heat capacity polynomial, :math:`C_p(T)`, from a defined reference temperature, :math:`T_{ref}`, to the current system temperature: +The dimensionless enthalpy departure is given by: .. math:: - C_p(\Delta T) = a_0 + a_1 \Delta T + a_2 \Delta T^2 + a_3 \Delta T^3 + a_4 \Delta T^4 + \frac{H_{dep}}{RT} = Z - 1 + \frac{A + T\frac{\partial A}{\partial T}}{B(\delta_1 - \delta_2)} \ln \left( \frac{Z + \delta_1 B}{Z + \delta_2 B} \right) -where :math:`\Delta T = T - T_{ref}`. The integrated polynomial is added to a user-supplied reference enthalpy for each component. +where :math:`Z` is the compressibility factor of the mixture. -**Departure Function** +This model is recommended whenever non-isothermal compositional physics are required. + +Mass formulation conversion +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +When the simulator is configured to use the mass formulation, the underlying inputs and thermodynamic calculations remain entirely in the molar formulation. The final phase molar enthalpy :math:`H_{total}` [J/mol] is converted to a mass-based phase enthalpy :math:`H_{mass}` [J/kg] by dividing by the phase molar weight :math:`M_{phase}`: + +.. math:: + H_{mass} = \frac{H_{total}}{M_{phase}} + +where the phase molar weight is the mole-fraction weighted sum of the component molar weights :math:`M_i`: + +.. math:: + M_{phase} = \sum_i x_i M_i + +Because all inputs must be provided in molar units, if a user's specific heat capacity coefficients or reference enthalpies are originally defined in a mass formulation (e.g., :math:`\hat{a}_{i,k}` in [J/(kg K^n)]), they must be scaled by the respective component's molar weight :math:`M_i` [kg/mol] prior to input to yield the molar coefficients :math:`a_{i,k}`: + +.. math:: + a_{i,k} = \hat{a}_{i,k} M_i -The departure enthalpy accounts for real-gas behaviour and molecular interactions. It is calculated exactly from the Equation of State by evaluating the partial derivative of the EoS attractive mixture parameter, :math:`a`, with respect to temperature. +Parameters +~~~~~~~~~~ -This model is recommended whenever non-isothermal compositional physics are required. +* ``enthalpyReferenceTemperature``: (Optional) A scalar defining the reference temperature :math:`T_{ref}` with a default of 298.15 K (25 C). Unit: [K]. +* ``referenceEnthalpy``: (Optional) A 2D array mapping the reference enthalpy for each component within each phase at the reference temperature. The dimensions must match ``[number of phases, number of components]``. Defaults to zero. Unit: [J/mol]. +* ``componentHeatCapacityCoefficients``: (Required) A 2D array containing the 5 polynomial coefficients (:math:`a_0` to :math:`a_4`) for the specific heat capacity of each component. The dimensions must strictly be ``[number of components, 5]``. To maintain dimensional consistency with a heat capacity in [J/(mol K)] and temperature difference in [K], the units for the specific coefficients are: -* **Catalog Name:** ``CompositionalEnthalpy`` -* **Parameters:** - * ``enthalpyReferenceTemperature``: The reference temperature :math:`T_{ref}`. - * ``referenceEnthalpy``: The enthalpy of each component at the reference temperature. - * ``componentHeatCapacityCoefficients``: The polynomial coefficients (:math:`a_0` to :math:`a_4`) for the heat capacity of each component. \ No newline at end of file + * :math:`a_{i,0}`: [J/(mol K)] + * :math:`a_{i,1}`: [J/(mol K^2)] + * :math:`a_{i,2}`: [J/(mol K^3)] + * :math:`a_{i,3}`: [J/(mol K^4)] + * :math:`a_{i,4}`: [J/(mol K^5)] diff --git a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst index dad45161a65..904cc9fd0a8 100644 --- a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst +++ b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst @@ -6,11 +6,13 @@ Unit and Variable Conversions Thermodynamic models, including Equations of State and Flash calculations, are fundamentally derived from molar balances and molar properties. Consequently, the internal computations within the compositional fluid framework are strictly executed using molar variables (e.g., molar fractions, molar densities). -**Mass Formulation Support** +Mass Formulation Support +----------------------------- While the core thermodynamics rely on molar formulations, many reservoir simulators frame their governing conservation equations in terms of mass to ensure strict mass conservation across numerical schemes. To bridge this gap seamlessly, the fluid model supports an automated conversion layer: -1. **Entry Conversion:** Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. -2. **Internal Computation:** All stability, flash, density, and viscosity calculations proceed purely in molar space. -3. **Exit Transformation:** Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. \ No newline at end of file + + 1. Entry Conversion: Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. + 2. Internal Computation: All stability, flash, density, and viscosity calculations proceed purely in molar space. + 3. Exit Transformation: Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. diff --git a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst index 786b26b32f8..e353d97d91f 100644 --- a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst +++ b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst @@ -3,34 +3,146 @@ Viscosity Viscosity dictates the fluid's resistance to shear flow and fundamentally controls phase mobilities in porous media. The framework provides multiple approaches to evaluate viscosity. -Constant Viscosity Model +Constant viscosity model ------------------------ -The simplest approach applies a static, user-defined viscosity for a specified phase, entirely independent of pressure, temperature, or composition. This is recommended solely for idealised simulations, synthetic testing, or when fluid flow is dominated entirely by pressure gradients rather than mobility contrasts. +The simplest approach applies a static, user-defined viscosity. While the viscosity does not change with state variables such as pressure, temperature, or composition, it is possible to assign different constant viscosities to the different fluid phases. This is recommended solely for idealised simulations, synthetic testing, or when fluid flow is dominated entirely by pressure gradients rather than mobility contrasts. -* **Catalog Name:** Selected by passing ``ConstantViscosity`` into the phase model formulation. -* **Parameters:** ``constantPhaseViscosity``. +Parameters +~~~~~~~~~~~~ -Lohrenz-Bray-Clark (LBC) Viscosity Model +* ``constantPhaseViscosity``: The constant viscosity value applied to each phase. Unit: [Pa.s]. + +Lohrenz-Bray-Clark (LBC) viscosity model ---------------------------------------- -The Lohrenz, Bray, and Clark (1964) correlation is the industry-standard methodology for calculating the viscosity of compositional hydrocarbon mixtures. +The Lohrenz, Bray, and Clark correlation (Lohrenz et al., 1964) is the industry-standard methodology for calculating the viscosity of compositional hydrocarbon mixtures. + +The model evaluates viscosity by combining a low-pressure dilute gas contribution with a dense fluid residual term. + +Dilute gas viscosity +~~~~~~~~~~~~~~~~~~~~~ + +The low-pressure, dilute gas viscosity :math:`\mu_{i}^*` of each pure component :math:`i` is estimated as a function of temperature using the correlations of Stiel and Thodos (1961). This relies on the viscosity-reducing parameter :math:`\xi_i`, defined as: + +.. math:: + + \xi_i = \frac{T_{c,i}^{1/6}}{M_i^{1/2} P_{c,i}^{2/3}} + +where :math:`T_{c,i}` is the critical temperature [K], :math:`P_{c,i}` is the critical pressure [atm], and :math:`M_i` is the molar weight [g/mol]. + +For standard non-polar components, the dilute viscosity (in centipoise) depends on the reduced temperature :math:`T_{r,i} = T / T_{c,i}`: + +.. math:: + + \mu_{i}^* = \frac{34.0 \times 10^{-5} T_{r,i}^{0.94}}{\xi_i} \quad \text{for } T_{r,i} \le 1.5 + +.. math:: + + \mu_{i}^* = \frac{17.78 \times 10^{-5} (4.58 T_{r,i} - 1.67)^{0.625}}{\xi_i} \quad \text{for } T_{r,i} > 1.5 + +For hydrogen gas (:math:`M_i < 2.1 \times 10^{-3}` kg/mol), the specific correlation is: + +.. math:: + + \mu_{H_2}^* = 90.71 \times 10^{-5} (0.1375 T - 1.67)^{0.625} + +These pure viscosities are then blended to find the dilute mixture viscosity, :math:`\mu^*`. The user may select between the Herning-Zipperer, Wilke, or Brokaw mixing rules to aggregate the individual components with mole fractions :math:`y_i`. + +The Herning-Zipperer mixing rule (Herning and Zipperer, 1936) evaluates the mixture viscosity as: + +.. math:: + + \mu^* = \frac{\sum_i y_i \mu_i^* \sqrt{M_i}}{\sum_i y_i \sqrt{M_i}} + +The Wilke mixing rule (Wilke, 1950) evaluates the mixture viscosity as: + +.. math:: + + \mu^* = \sum_i \frac{y_i \mu_i^*}{\sum_j y_j \phi_{ij}} + +where the interaction parameter :math:`\phi_{ij}` is defined as: + +.. math:: + + \phi_{ij} = \frac{\left[1 + (\mu_i^*/\mu_j^*)^{1/2} (M_i/M_j)^{-1/4}\right]^2}{\sqrt{8(1 + M_i/M_j)}} + +The Brokaw mixing rule (Brokaw, R. S., 1968, Viscosity of Gas Mixtures, National Aeronautics and Space Administration) evaluates the mixture viscosity as: + +.. math:: + + \mu^* = \sum_i \frac{y_i \mu_i^*}{\sum_j y_j \phi_{ij} \sqrt{\mu_i^*/\mu_j^*}} -The model evaluates viscosity in two distinct steps: -1. **Dilute Gas Viscosity:** The viscosity of each pure component at low-pressure (dilute) conditions is estimated as a function of temperature using the Stiel and Thodos (1961) correlations. These pure viscosities are then blended to find the dilute mixture viscosity. The user may select between the Herning-Zipperer, Wilke, or Brokaw mixing rules. -2. **Dense Fluid Contribution:** A residual viscosity, representing the effects of elevated pressure and density, is calculated using a fourth-order polynomial dependent on the reduced density of the phase. This residual is added to the dilute gas viscosity to obtain the final phase viscosity. +where the interaction parameter :math:`\phi_{ij}` depends on the molar weight ratios :math:`A_{ij} = M_i / M_j` and :math:`B_{ij} = [4 M_i M_j / (M_i + M_j)^2]^{0.25}`: + +.. math:: + + \phi_{ij} = \frac{B_{ij}}{\sqrt{A_{ij}}} \left[ 1 + \frac{A_{ij} - A_{ij}^{0.45}}{2 + 2A_{ij} + B_{ij}(1+A_{ij}^{0.45})/(1+B_{ij})} \right] + +Dense fluid contribution +~~~~~~~~~~~~~~~~~~~~~~~~~ + +A residual viscosity, representing the effects of elevated pressure and density, is calculated using the Lohrenz et al. (1964) correlation. This evaluates a fourth-order polynomial dependent on the reduced density of the phase. + +First, the mixture critical parameters are calculated using Kay's mixing rule: + +.. math:: + + T_c = \sum_i y_i T_{c,i}, \quad P_c = \sum_i y_i P_{c,i}, \quad V_c = \sum_i y_i V_{c,i}, \quad M = \sum_i y_i M_i + +The mixture viscosity-reducing parameter :math:`\xi` and reduced density :math:`\rho_r` are defined using the phase molar density :math:`\rho_m`: + +.. math:: + + \xi = \frac{T_c^{1/6}}{M^{1/2} P_c^{2/3}}, \quad \rho_r = \rho_m V_c + +The dense fluid polynomial :math:`f(\rho_r)` is then evaluated: + +.. math:: + + f(\rho_r) = 0.1023 + 0.023364 \rho_r + 0.058533 \rho_r^2 - 0.040758 \rho_r^3 + 0.0093324 \rho_r^4 + +This residual is added to the dilute gas viscosity to obtain the final phase viscosity :math:`\mu`: + +.. math:: + + \mu = \mu^* + \frac{f(\rho_r)^4 - 10^{-4}}{\xi} This model is strongly recommended for all miscible gas, volatile oil, and general compositional hydrocarbon simulations. -* **Catalog Name:** ``LohrenzBrayClark`` -* **Parameters:** - * ``viscosityMixingRule``: Defines the dilute gas mixing rule (``HerningZipperer``, ``Wilke``, or ``Brokaw``). - * ``componentCriticalVolume``: Critical volumes required for the reduced density calculation. +Parameters +~~~~~~~~~~~~ -Phillips Brine Viscosity Model +* ``viscosityMixingRule``: Defines the dilute gas mixing rule. Valid options are ``HerningZipperer``, ``Wilke``, or ``Brokaw``. Unit: Dimensionless (String). +* ``componentCriticalVolume``: The critical volumes of each component, required for the reduced density calculation. Unit: [m^3/mol]. + +Phillips brine viscosity model ------------------------------ -For aqueous phases, the Phillips Brine Viscosity model uses empirical correlations from Phillips et al. (1981) to calculate viscosity as a function of temperature, pressure, and salinity. The model relies on tabulated pure water saturation viscosities and scales them using a salinity-dependent multiplier. This model is recommended for saline aquifers and CO2 sequestration. +For aqueous phases, the Phillips brine viscosity model uses empirical correlations from Phillips et al. (1981) to calculate viscosity as a function of temperature and salinity. In this specific implementation, the brine viscosity is evaluated independent of pressure. This model is recommended for saline aquifers and CO2 sequestration. + +The model determines the final brine viscosity :math:`\mu` by applying a salinity- and temperature-dependent multiplier to the pure water viscosity :math:`\mu_w`: + +.. math:: + + \mu = \mu_w(T) \left[ c_0(S) + c_1(S) T \right] + +where :math:`T` is the temperature in degrees Celsius and :math:`S` is the salinity in mol/kg. The pure water viscosity :math:`\mu_w(T)` is not evaluated from a closed-form analytical equation; rather, it is interpolated from a pre-computed tabulated function of pure water saturation viscosities corresponding to the given temperature. + +The coefficients :math:`c_0(S)` and :math:`c_1(S)` dictate the scaling effect of the dissolved salts and are evaluated as follows: + +.. math:: + + c_0(S) = 1 + a S + b S^2 + c S^3 + +.. math:: + + c_1(S) = d (1 - \exp(k S)) + +The correlation uses the specific empirical constants :math:`a = 0.0816`, :math:`b = 0.0122`, :math:`c = 0.000128`, :math:`d = 0.000629`, and :math:`k = -0.7`. + +Parameters +~~~~~~~~~~~~ -* **Catalog Name:** ``PhillipsBrine`` -* **Parameters:** Tabulated via function managers; relies on the global ``salinity`` parameter. \ No newline at end of file +* ``salinity``: The salinity of the brine. Unit: [mol/kg]. +* ``temperatureCoordinates``: List of temperature values for the interpolation of tabulated pure water properties. Unit: [K]. From ae59072e93c808fdbe3c71560f3e95df3d5b2482 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 12:56:45 -0500 Subject: [PATCH 09/13] Update units --- .../docs/compositional/Density.rst | 16 ++--- .../docs/compositional/EquationsOfState.rst | 16 ++--- .../compositional/ImmiscibleWaterFlash.rst | 10 +--- .../docs/compositional/KValueFlash.rst | 4 +- .../docs/compositional/Miscellaneous.rst | 6 +- .../docs/compositional/NegativeFlash.rst | 10 ++-- .../docs/compositional/References.rst | 59 +++++++++++++++---- .../docs/compositional/Viscosity.rst | 2 +- 8 files changed, 78 insertions(+), 45 deletions(-) diff --git a/src/coreComponents/constitutive/docs/compositional/Density.rst b/src/coreComponents/constitutive/docs/compositional/Density.rst index d63bee00b86..6598096cf59 100644 --- a/src/coreComponents/constitutive/docs/compositional/Density.rst +++ b/src/coreComponents/constitutive/docs/compositional/Density.rst @@ -40,7 +40,7 @@ where :math:`MW_i` is the molecular weight of component :math:`i`. This model is Parameters ~~~~~~~~~~~~~~~~ -* ``componentVolumeShift``: Component-specific volume shift parameters (dimensionless). +* ``componentVolumeShift``: Component-specific volume shift parameters. Unit: [dimensionless]. Phillips brine density model ---------------------------- @@ -93,8 +93,8 @@ The phase mass density, :math:`\rho_{mass}`, is then computed by multiplying the Parameters ~~~~~~~~~~~~~~~~ -* ``salinity``: Brine salinity (mol/kg). -* ``saltMolarWeight``: The molar weight of the salt component (kg/mol). +* ``salinity``: Brine salinity. Unit: [mol/kg]. +* ``saltMolarWeight``: The molar weight of the salt component. Unit: [kg/mol]. Immiscible water density ------------------------ @@ -112,8 +112,8 @@ The phase molar density, :math:`\rho_{molar}`, is then calculated by dividing th Parameters ~~~~~~~~~~~~~~~~ -* ``waterReferencePressure``: Reference pressure for the water density calculation (Pa). -* ``waterReferenceTemperature``: Reference temperature for the water density calculation (K). -* ``waterDensity``: Water mass density at the reference pressure and temperature (kg/m^3). -* ``waterCompressibility``: Isothermal compressibility of water (1/Pa). -* ``waterExpansionCoefficient``: Volumetric coefficient of thermal expansion of water (1/K). \ No newline at end of file +* ``waterReferencePressure``: Reference pressure for the water density calculation. Unit: [Pa]. +* ``waterReferenceTemperature``: Reference temperature for the water density calculation. Unit: [K]. +* ``waterDensity``: Water mass density at the reference pressure and temperature. Unit: [kg/m^3]. +* ``waterCompressibility``: Isothermal compressibility of water. Unit: [1/Pa]. +* ``waterExpansionCoefficient``: Volumetric coefficient of thermal expansion of water. Unit: [1/K]. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst index 9f213d2fe94..e2aeaa4d35d 100644 --- a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst +++ b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst @@ -154,14 +154,14 @@ For carbon dioxide (CO2): .. math:: k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i} + A_2 a_2 -where :math:`A_0 = -0.31092`, :math:`A_1 = 0.23580`, :math:`A_2 = -21.2566`, :math:`a_0 = 1.0 + 0.15587 c_{sw}^{0.7505}`, :math:`a_1 = 1.0 + 0.17837 c_{sw}^{0.979}`, and :math:`a_2 = \exp(-6.7222 T_{r,i} - c_{sw})`. +where :math:`A_0 = -0.31092`, :math:`A_1 = 0.23580`, :math:`A_2 = -21.2566`, :math:`a_0 = 1 + 0.15587 c_{sw}^{0.7505}`, :math:`a_1 = 1 + 0.17837 c_{sw}^{0.979}`, and :math:`a_2 = \exp(-6.7222 T_{r,i} - c_{sw})`. For nitrogen (N2): .. math:: k_{ij}^{AQ} = A_0 a_0 + A_1 a_1 T_{r,i} -where :math:`A_0 = -1.70235`, :math:`A_1 = 0.44338`, :math:`a_0 = 1.0 + 0.025587 c_{sw}^{0.75}`, and :math:`a_1 = 1.0 + 0.08126 c_{sw}^{0.75}`. +where :math:`A_0 = -1.70235`, :math:`A_1 = 0.44338`, :math:`a_0 = 1 + 0.025587 c_{sw}^{0.75}`, and :math:`a_1 = 1 + 0.08126 c_{sw}^{0.75}`. For hydrogen sulfide (H2S): @@ -216,9 +216,9 @@ Here is an example demonstrating how to specify the parameters for a two-phase s salinity="0.6" /> * ``equationsOfState``: List specifying the EoS type (e.g., ``PengRobinson``, ``SoaveRedlichKwong``, ``SoreideWhitson``) corresponding to each phase. -* ``componentCriticalPressure``: Critical pressure of each component [Pa]. -* ``componentCriticalTemperature``: Critical temperature of each component [K]. -* ``componentAcentricFactor``: Acentric factor of each component. -* ``componentVolumeShift``: Volume shift parameter for each component (optional). -* ``componentBinaryCoeff``: Flattened :math:`N \times N` matrix of binary interaction coefficients (optional). -* ``salinity``: Brine salinity [mol/kg] utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). \ No newline at end of file +* ``componentCriticalPressure``: Critical pressure of each component. Unit: [Pa]. +* ``componentCriticalTemperature``: Critical temperature of each component. Unit: [K]. +* ``componentAcentricFactor``: Acentric factor of each component. Unit: [dimensionless]. +* ``componentVolumeShift``: Volume shift parameter for each component (optional). Unit: [dimensionless]. +* ``componentBinaryCoeff``: Flattened :math:`N \times N` matrix of binary interaction coefficients (optional). Unit: [dimensionless]. +* ``salinity``: Brine salinity utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). Unit: [mol/kg]. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst index 2134908ade0..b9c71abbe0b 100644 --- a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst @@ -16,7 +16,7 @@ Given a total feed mole fraction :math:`z_i` for each component, the mole fracti The composition of the aqueous phase is fixed at 100% water. The remaining feed defines the total hydrocarbon mole fraction, :math:`z_{hc}`: .. math:: - z_{hc} = 1.0 - z_{w} + z_{hc} = 1 - z_{w} The feed composition of the remaining hydrocarbon mixture is then normalised: @@ -40,7 +40,7 @@ Once the two-phase hydrocarbon flash converges, the overall phase fractions for V_{V} = z_{hc} V_{V,hc} .. math:: - V_{L} = z_{hc} (1.0 - V_{V,hc}) + V_{L} = z_{hc} (1 - V_{V,hc}) The derivatives of the phase fractions and compositions with respect to pressure, temperature, and total composition are analytically transformed using the chain rule to account for this scaling and normalisation. @@ -79,9 +79,3 @@ Parameters ~~~~~~~~~~~~~~~~ * ``equationsOfState``: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment. -* ``waterReferencePressure``: The reference pressure for water density and viscosity [Pa]. -* ``waterReferenceTemperature``: The reference temperature for water density and viscosity [K] (optional, default: 293.15). -* ``waterDensity``: The water density at the reference pressure and temperature [kg/m^3]. -* ``waterViscosity``: The water viscosity at the reference pressure and temperature [Pa.s]. -* ``waterCompressibility``: The constant isothermal compressibility of water [1/Pa]. -* ``waterExpansionCoefficient``: The volumetric coefficient of thermal expansion of water [1/K] (optional, default: 0.0). \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst index fa2763d2126..cc8c676ed33 100644 --- a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst @@ -71,5 +71,5 @@ Parameters ~~~~~~~~~~~~~~~~ * ``kValueTables``: List of function names linking to the pre-tabulated K-values. For an :math:`N`-phase system, this requires :math:`(N-1) \times N_c` function names. -* ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). -* ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). \ No newline at end of file +* ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). Unit: [Pa]. +* ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). Unit: [K]. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst index 904cc9fd0a8..51fecb8babd 100644 --- a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst +++ b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst @@ -13,6 +13,6 @@ While the core thermodynamics rely on molar formulations, many reservoir simulat To bridge this gap seamlessly, the fluid model supports an automated conversion layer: - 1. Entry Conversion: Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. - 2. Internal Computation: All stability, flash, density, and viscosity calculations proceed purely in molar space. - 3. Exit Transformation: Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. +1. Entry Conversion: Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. +2. Internal Computation: All stability, flash, density, and viscosity calculations proceed purely in molar space. +3. Exit Transformation: Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. diff --git a/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst index ef84be59441..5f3f676ed65 100644 --- a/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/NegativeFlash.rst @@ -69,8 +69,8 @@ The negative two-phase flash is the recommended standard for rigorous compositio Parameters ~~~~~~~~~~~~~~~~ -* ``stabilityThreshold``: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). -* ``stabilityTolerance``: Tolerance for stationarity in the stability test. -* ``stabilityMaxIterations``: Maximum successive substitution steps for stability analysis. -* ``flashTolerance``: Convergence tolerance for the fugacity ratio error. -* ``flashMaxIterations``: Maximum successive substitution steps for the flash solve. \ No newline at end of file +* ``stabilityThreshold``: Tangent plane distance below which a mixture is unstable (default: -1.0e-8). Unit: [dimensionless]. +* ``stabilityTolerance``: Tolerance for stationarity in the stability test. Unit: [dimensionless]. +* ``stabilityMaxIterations``: Maximum successive substitution steps for stability analysis. Unit: [dimensionless]. +* ``flashTolerance``: Convergence tolerance for the fugacity ratio error. Unit: [dimensionless]. +* ``flashMaxIterations``: Maximum successive substitution steps for the flash solve. Unit: [dimensionless]. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst index 5fbcfd263f6..dee85471d57 100644 --- a/src/coreComponents/constitutive/docs/compositional/References.rst +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -1,13 +1,52 @@ References ========== -The compositional fluid models implemented in this framework are based heavily on established thermodynamic and empirical literature. The following references document the mathematical correlations and mixing rules utilised: - -* **Stiel, L. I., & Thodos, G. (1961).** The viscosity of nonpolar gases at normal pressures. *AIChE Journal*, 7(4), 611-615. -* **Herning, F., & Zipperer, L. (1936).** Calculation of the viscosity of technical gas mixtures from the viscosity of individual gases. *Gas u. Wasserfach*, 79(49), 69. -* **Wilke, C. R. (1950).** A viscosity equation for gas mixtures. *The Journal of Chemical Physics*, 18(4), 517-519. -* **Brokaw, R. S. (1968).** Viscosity of gas mixtures. *National Aeronautics and Space Administration*. -* **Lohrenz, J., Bray, B. G., & Clark, C. R. (1964).** Calculating viscosities of reservoir fluids from their compositions. *Journal of Petroleum Technology*, 16(10), 1171-1176. -* **Phillips, S. L., Igbene, A., Fair, J. A., Ozbek, H., & Tavana, M. (1981).** A technical databook for geothermal energy utilization. *Lawrence Berkeley National Laboratory*. -* **Soreide, I., & Whitson, C. H. (1992).** Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2S with tie water. *Fluid Phase Equilibria*, 77, 217-240. -* **Peneloux, A., Rauzy, E., & Freze, R. (1982).** A consistent correction for Redlich-Kwong-Soave volumes. *Fluid Phase Equilibria*, 8(1), 7-23. \ No newline at end of file +- Brokaw, R. S. (1968). + Viscosity of gas mixtures. *Washington, D.C.: National Aeronautics and Space Administration*. +- Chabab, S., Kerkache, H., Bouchkira, I., Poulain, M., Baudouin, O., Moine, É., Ducousso, M., Hoang, H., Galliéro, G., and Cézac, P. (2024). + Solubility of H2 in water and NaCl brine under subsurface storage conditions: Measurements and thermodynamic modeling. *International Journal of Hydrogen Energy*, **50**, 648–658. + ``__ +- Herning, F., and Zipperer, L. (1936). + Beitrag zur Berechnung der Zähigkeit Technischer Gasgemische aus den Zähigkeitswerten der Einzelbestandteile; Calculation of the Viscosity of Technical Gas Mixtures from Viscosity of the individual Gases. *Das Gas- und Wasserfach*, **79**, 49–54 and 69–73. +- Lohrenz, J., Bray, B. G., and Clark, C. R. (1964). + Calculating Viscosities of Reservoir Fluids From Their Compositions. *Journal of Petroleum Technology*, **16** (10), 1171–1176. + ``__ +- Michelsen, M. L. (1982a). + The isothermal flash problem. Part I. Stability. *Fluid Phase Equilibria*, **9** (1), 1–19. + ``__ +- Michelsen, M. L. (1982b). + The isothermal flash problem. Part II. Phase-split calculation. *Fluid Phase Equilibria*, **9** (1), 21–40. + ``__ +- Péneloux, A., Rauzy, E., and Fréze, R. (1982). + A consistent correction for Redlich-Kwong-Soave volumes. *Fluid Phase Equilibria*, **8** (1), 7–23. + ``__ +- Peng, D.-Y., and Robinson, D. B. (1976). + A New Two-Constant Equation of State. *Industrial & Engineering Chemistry Fundamentals*, **15** (1), 59–64. + ``__ +- Phillips, S. L., Igbene, A., Fair, J. A., Ozbek, H., and Tavana, M. (1981). + A technical databook for geothermal energy utilization. *Lawrence Berkeley National Laboratory*. +- Poling, B. E., Prausnitz, J. M., and O'Connell, J. P. (2001). + The Properties of Gases and Liquids. *McGraw-Hill*. +- Rachford, H. H., Jr., and Rice, J. D. (1952). + Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. *Journal of Petroleum Technology*, **4** (10), 19–3. + ``__ +- Robinson, D. B., and Peng, D. (1978). + The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs. *United States: Gas Processors Association*. +- Soave, G. (1972). + Equilibrium constants from a modified Redlich-Kwong equation of state. *Chemical Engineering Science*, **27** (6), 1197–1203. + ``__ +- Søreide, I., and Whitson, C. H. (1992). + Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2 S with pure water and NaCI brine. *Fluid Phase Equilibria*, **77**, 217–240. + ``__ +- Stiel, L. I., and Thodos, G. (1961). + The viscosity of nonpolar gases at normal pressures. *AIChE Journal*, **7** (4), 611–615. + ``__ +- Whitson, C. H., and Torp, S. B. (1981). + Evaluating Constant Volume Depletion Data. *SPE Annual Technical Conference and Exhibition*. + ``__ +- Wilke, C. R. (1950). + A Viscosity Equation for Gas Mixtures. *The Journal of Chemical Physics*, **18** (4), 517–519. + ``__ +- Wilson, G. M. (1969). + A modified Redlich-Kwong equation of state, application to general physical data calculations. *65th National AIChE Meeting, Cleveland, OH*, **15**. + \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst index e353d97d91f..cc509a92317 100644 --- a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst +++ b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst @@ -113,7 +113,7 @@ This model is strongly recommended for all miscible gas, volatile oil, and gener Parameters ~~~~~~~~~~~~ -* ``viscosityMixingRule``: Defines the dilute gas mixing rule. Valid options are ``HerningZipperer``, ``Wilke``, or ``Brokaw``. Unit: Dimensionless (String). +* ``viscosityMixingRule``: Defines the dilute gas mixing rule. Valid options are ``HerningZipperer``, ``Wilke``, or ``Brokaw``. * ``componentCriticalVolume``: The critical volumes of each component, required for the reduced density calculation. Unit: [m^3/mol]. Phillips brine viscosity model From 086d0571636fdfd2b0d7261dfc69355b87b0549e Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 13:01:44 -0500 Subject: [PATCH 10/13] Add immisciblel viscosity model --- .../docs/compositional/Viscosity.rst | 20 +++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst index cc509a92317..abc8020bf80 100644 --- a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst +++ b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst @@ -146,3 +146,23 @@ Parameters * ``salinity``: The salinity of the brine. Unit: [mol/kg]. * ``temperatureCoordinates``: List of temperature values for the interpolation of tabulated pure water properties. Unit: [K]. + +Immiscible water viscosity +-------------------------- + +A simplified exponential viscosity model is available for pure, immiscible water phases. It evaluates the viscosity, :math:`\mu`, based on pressure and temperature variations from a reference state: + +.. math:: + + \mu = \mu_{ref} \exp(c_{\mu} (P - P_{ref})) \exp(-\alpha_{\mu} (T - T_{ref})) + +where :math:`\mu_{ref}` is the reference viscosity, :math:`c_{\mu}` is the viscosity compressibility, and :math:`\alpha_{\mu}` is the viscosity expansion coefficient. + +Parameters +~~~~~~~~~~~~ + +* ``waterReferencePressure``: Reference pressure for the water viscosity calculation. Unit: [Pa]. +* ``waterReferenceTemperature``: Reference temperature for the water viscosity calculation. Unit: [K]. +* ``waterViscosity``: Water viscosity at the reference pressure and temperature. Unit: [Pa.s]. +* ``waterViscosityCompressibility``: The compressibility (normalized derivative with respect to pressure) of the water viscosity. Unit: [1/Pa]. +* ``waterViscosityExpansionCoefficient``: The coefficient of thermal expansion (normalized derivative with respect to temperature) of water viscosity. Unit: [1/K]. From 8161e8f1f6d4df8b5564b874f8da9bfc34d9b1e1 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 13:13:05 -0500 Subject: [PATCH 11/13] Ihmels law --- .../docs/compositional/References.rst | 28 ++++++++++--------- .../docs/compositional/Viscosity.rst | 8 +++++- 2 files changed, 22 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/constitutive/docs/compositional/References.rst b/src/coreComponents/constitutive/docs/compositional/References.rst index dee85471d57..57c7a2ff9dc 100644 --- a/src/coreComponents/constitutive/docs/compositional/References.rst +++ b/src/coreComponents/constitutive/docs/compositional/References.rst @@ -4,49 +4,51 @@ References - Brokaw, R. S. (1968). Viscosity of gas mixtures. *Washington, D.C.: National Aeronautics and Space Administration*. - Chabab, S., Kerkache, H., Bouchkira, I., Poulain, M., Baudouin, O., Moine, É., Ducousso, M., Hoang, H., Galliéro, G., and Cézac, P. (2024). - Solubility of H2 in water and NaCl brine under subsurface storage conditions: Measurements and thermodynamic modeling. *International Journal of Hydrogen Energy*, **50**, 648–658. + Solubility of H2 in water and NaCl brine under subsurface storage conditions: Measurements and thermodynamic modeling. *International Journal of Hydrogen Energy*, **50**, 648-658. ``__ - Herning, F., and Zipperer, L. (1936). - Beitrag zur Berechnung der Zähigkeit Technischer Gasgemische aus den Zähigkeitswerten der Einzelbestandteile; Calculation of the Viscosity of Technical Gas Mixtures from Viscosity of the individual Gases. *Das Gas- und Wasserfach*, **79**, 49–54 and 69–73. + Beitrag zur Berechnung der Zähigkeit Technischer Gasgemische aus den Zähigkeitswerten der Einzelbestandteile; Calculation of the Viscosity of Technical Gas Mixtures from Viscosity of the individual Gases. *Das Gas- und Wasserfach*, **79**, 49-54 and 69-73. +- Ihmels, E. C. (2010). + The Critical Surface. *Journal of Chemical & Engineering Data*, **55** (9), 3474-3480. + ``__ - Lohrenz, J., Bray, B. G., and Clark, C. R. (1964). - Calculating Viscosities of Reservoir Fluids From Their Compositions. *Journal of Petroleum Technology*, **16** (10), 1171–1176. + Calculating Viscosities of Reservoir Fluids From Their Compositions. *Journal of Petroleum Technology*, **16** (10), 1171-1176. ``__ - Michelsen, M. L. (1982a). - The isothermal flash problem. Part I. Stability. *Fluid Phase Equilibria*, **9** (1), 1–19. + The isothermal flash problem. Part I. Stability. *Fluid Phase Equilibria*, **9** (1), 1-19. ``__ - Michelsen, M. L. (1982b). - The isothermal flash problem. Part II. Phase-split calculation. *Fluid Phase Equilibria*, **9** (1), 21–40. + The isothermal flash problem. Part II. Phase-split calculation. *Fluid Phase Equilibria*, **9** (1), 21-40. ``__ - Péneloux, A., Rauzy, E., and Fréze, R. (1982). - A consistent correction for Redlich-Kwong-Soave volumes. *Fluid Phase Equilibria*, **8** (1), 7–23. + A consistent correction for Redlich-Kwong-Soave volumes. *Fluid Phase Equilibria*, **8** (1), 7-23. ``__ - Peng, D.-Y., and Robinson, D. B. (1976). - A New Two-Constant Equation of State. *Industrial & Engineering Chemistry Fundamentals*, **15** (1), 59–64. + A New Two-Constant Equation of State. *Industrial & Engineering Chemistry Fundamentals*, **15** (1), 59-64. ``__ - Phillips, S. L., Igbene, A., Fair, J. A., Ozbek, H., and Tavana, M. (1981). A technical databook for geothermal energy utilization. *Lawrence Berkeley National Laboratory*. - Poling, B. E., Prausnitz, J. M., and O'Connell, J. P. (2001). The Properties of Gases and Liquids. *McGraw-Hill*. - Rachford, H. H., Jr., and Rice, J. D. (1952). - Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. *Journal of Petroleum Technology*, **4** (10), 19–3. + Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. *Journal of Petroleum Technology*, **4** (10), 19-3. ``__ - Robinson, D. B., and Peng, D. (1978). The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs. *United States: Gas Processors Association*. - Soave, G. (1972). - Equilibrium constants from a modified Redlich-Kwong equation of state. *Chemical Engineering Science*, **27** (6), 1197–1203. + Equilibrium constants from a modified Redlich-Kwong equation of state. *Chemical Engineering Science*, **27** (6), 1197-1203. ``__ - Søreide, I., and Whitson, C. H. (1992). - Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2 S with pure water and NaCI brine. *Fluid Phase Equilibria*, **77**, 217–240. + Peng-Robinson predictions for hydrocarbons, CO2, N2, and H2 S with pure water and NaCI brine. *Fluid Phase Equilibria*, **77**, 217-240. ``__ - Stiel, L. I., and Thodos, G. (1961). - The viscosity of nonpolar gases at normal pressures. *AIChE Journal*, **7** (4), 611–615. + The viscosity of nonpolar gases at normal pressures. *AIChE Journal*, **7** (4), 611-615. ``__ - Whitson, C. H., and Torp, S. B. (1981). Evaluating Constant Volume Depletion Data. *SPE Annual Technical Conference and Exhibition*. ``__ - Wilke, C. R. (1950). - A Viscosity Equation for Gas Mixtures. *The Journal of Chemical Physics*, **18** (4), 517–519. + A Viscosity Equation for Gas Mixtures. *The Journal of Chemical Physics*, **18** (4), 517-519. ``__ - Wilson, G. M. (1969). A modified Redlich-Kwong equation of state, application to general physical data calculations. *65th National AIChE Meeting, Cleveland, OH*, **15**. - \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst index abc8020bf80..f45418608d9 100644 --- a/src/coreComponents/constitutive/docs/compositional/Viscosity.rst +++ b/src/coreComponents/constitutive/docs/compositional/Viscosity.rst @@ -67,7 +67,7 @@ where the interaction parameter :math:`\phi_{ij}` is defined as: \phi_{ij} = \frac{\left[1 + (\mu_i^*/\mu_j^*)^{1/2} (M_i/M_j)^{-1/4}\right]^2}{\sqrt{8(1 + M_i/M_j)}} -The Brokaw mixing rule (Brokaw, R. S., 1968, Viscosity of Gas Mixtures, National Aeronautics and Space Administration) evaluates the mixture viscosity as: +The Brokaw mixing rule (Brokaw, 1968) evaluates the mixture viscosity as: .. math:: @@ -110,6 +110,12 @@ This residual is added to the dilute gas viscosity to obtain the final phase vis This model is strongly recommended for all miscible gas, volatile oil, and general compositional hydrocarbon simulations. +The model relies on critical component volumes; when these are not provided, they are estimated internally using a correlation due to Ihmels (2010) + +.. math:: + + V_c = \frac{2.215 \times 10^{-6} T_c}{0.025 + 10^{-6} P_c} + Parameters ~~~~~~~~~~~~ From ddc0f197a0fbf4a903557f2b3c12795ccfa24e94 Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 15:34:04 -0500 Subject: [PATCH 12/13] Tidy links --- .../docs/CompositionalMultiphaseFluid.rst | 2 ++ .../docs/compositional/Density.rst | 3 ++- .../docs/compositional/EquationsOfState.rst | 22 ++------------- .../compositional/ImmiscibleWaterFlash.rst | 27 +++---------------- .../docs/compositional/KValueFlash.rst | 27 ------------------- .../docs/compositional/Miscellaneous.rst | 6 +++++ .../constitutiveDrivers/docs/PVTDriver.rst | 2 ++ 7 files changed, 18 insertions(+), 71 deletions(-) diff --git a/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst b/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst index 8f4a8a9f2f3..ce1c2c3a628 100644 --- a/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst +++ b/src/coreComponents/constitutive/docs/CompositionalMultiphaseFluid.rst @@ -36,4 +36,6 @@ This documentation is divided into the following sections: .. include:: compositional/Miscellaneous.rst +.. include:: compositional/Models.rst + .. include:: compositional/References.rst diff --git a/src/coreComponents/constitutive/docs/compositional/Density.rst b/src/coreComponents/constitutive/docs/compositional/Density.rst index 6598096cf59..cb1e2e47d2c 100644 --- a/src/coreComponents/constitutive/docs/compositional/Density.rst +++ b/src/coreComponents/constitutive/docs/compositional/Density.rst @@ -88,7 +88,8 @@ The phase mass density, :math:`\rho_{mass}`, is then computed by multiplying the .. math:: \rho_{mass} = \rho_{molar} \left( x_{H2O} MW_{brine} + \sum_{i \neq H2O} x_i MW_i \right) -*(Note: While the volume shift tightly couples the Phillips correlation to the EOS density, the brine viscosity is not derived from this density. Instead, the phase viscosity is modelled independently using a distinct temperature- and salinity-dependent multiplier applied to the pure water viscosity).* +.. note:: + While the volume shift tightly couples the Phillips correlation to the EOS density, the brine viscosity is not derived from this density. Instead, the phase viscosity is modelled independently using a distinct temperature- and salinity-dependent multiplier applied to the pure water viscosity. Parameters ~~~~~~~~~~~~~~~~ diff --git a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst index e2aeaa4d35d..84588820740 100644 --- a/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst +++ b/src/coreComponents/constitutive/docs/compositional/EquationsOfState.rst @@ -195,25 +195,7 @@ If a component's name does not match any of these predefined strings (for exampl Model parameters ~~~~~~~~~~~~~~~~ -Equations of state are assigned per-phase within the specific fluid model XML block (e.g., ````). - -The ``equationsOfState`` attribute takes a list specifying the EoS type for each phase corresponding to the order defined in the ``phaseNames`` list. Standard component properties necessary for the EoS calculations must also be provided in matching component order. - -Here is an example demonstrating how to specify the parameters for a two-phase system using the Soreide-Whitson EoS for the aqueous phase and Peng-Robinson for the gas phase: - -.. code-block:: xml - - +Equations of state are assigned per-phase within the specific fluid model XML block. The ``equationsOfState`` attribute takes a list specifying the EoS type for each phase corresponding to the order defined in the ``phaseNames`` list. Standard component properties necessary for the EoS calculations must also be provided in matching component order. * ``equationsOfState``: List specifying the EoS type (e.g., ``PengRobinson``, ``SoaveRedlichKwong``, ``SoreideWhitson``) corresponding to each phase. * ``componentCriticalPressure``: Critical pressure of each component. Unit: [Pa]. @@ -221,4 +203,4 @@ Here is an example demonstrating how to specify the parameters for a two-phase s * ``componentAcentricFactor``: Acentric factor of each component. Unit: [dimensionless]. * ``componentVolumeShift``: Volume shift parameter for each component (optional). Unit: [dimensionless]. * ``componentBinaryCoeff``: Flattened :math:`N \times N` matrix of binary interaction coefficients (optional). Unit: [dimensionless]. -* ``salinity``: Brine salinity utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). Unit: [mol/kg]. \ No newline at end of file +* ``salinity``: Brine salinity utilized by the Soreide-Whitson aqueous model (optional, defaults to 0.0). Unit: [mol/kg]. diff --git a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst index b9c71abbe0b..4532bab29b0 100644 --- a/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/ImmiscibleWaterFlash.rst @@ -3,10 +3,13 @@ Immiscible water flash The immiscible water flash is an efficient three-phase equilibrium model designed for systems containing an aqueous phase alongside hydrocarbon liquid and vapour phases. Instead of solving a fully coupled three-phase equation of state problem, this model vastly simplifies the phase split by assuming water forms a strictly pure aqueous phase and is completely immiscible in the hydrocarbon phases, while hydrocarbons are entirely insoluble in the aqueous phase. +.. note:: + When using the immiscibile water flash, the water component still needs to be explicitly specified as one of the components with the name ``H2O``. + Methodology ~~~~~~~~~~~ -In this formulation, the mixture composition is explicitly separated into a water component and a hydrocarbon mixture before any thermodynamic equilibrium calculations occur. The water component is identified using its assigned component name (e.g., ``H2O``). +In this formulation, the mixture composition is explicitly separated into a water component and a hydrocarbon mixture before any thermodynamic equilibrium calculations occur. The water component is identified using its assigned component name (``H2O``). Given a total feed mole fraction :math:`z_i` for each component, the mole fraction of the aqueous phase, :math:`V_{aq}`, is determined to be exactly equal to the total feed mole fraction of the water component: @@ -56,26 +59,4 @@ Model parameters The immiscible water flash fluid models are assigned using specific XML blocks that pair the three-phase flash model with distinct viscosity and density models for the hydrocarbon and aqueous phases. -Here is an example demonstrating how to specify the parameters for a three-phase system using the immiscible water flash paired with the Lohrenz-Bray-Clark viscosity model for hydrocarbons and a distinct constant-compressibility formulation for the pure water phase: - -.. code-block:: xml - - - -Parameters -~~~~~~~~~~~~~~~~ - * ``equationsOfState``: List specifying the EoS type for each phase (liquid, vapour, aqueous). Note that the aqueous EoS entry is largely ignored by the pure-water density model but must be provided for list alignment. diff --git a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst index cc8c676ed33..a621cee2cb0 100644 --- a/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst +++ b/src/coreComponents/constitutive/docs/compositional/KValueFlash.rst @@ -43,33 +43,6 @@ Model parameters The K-value flash fluid models are assigned using specific XML blocks that pair the flash model with viscosity and density models. -Here is an example demonstrating how to specify the parameters for a two-phase system using the K-value flash paired with the Phillips brine model: - -.. code-block:: xml - - - -Parameters -~~~~~~~~~~~~~~~~ - * ``kValueTables``: List of function names linking to the pre-tabulated K-values. For an :math:`N`-phase system, this requires :math:`(N-1) \times N_c` function names. * ``pressureCoordinates``: List of pressure values for interpolation grid generation (optional). Unit: [Pa]. * ``temperatureCoordinates``: List of temperature values for interpolation grid generation (optional). Unit: [K]. \ No newline at end of file diff --git a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst index 51fecb8babd..e65a5587d79 100644 --- a/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst +++ b/src/coreComponents/constitutive/docs/compositional/Miscellaneous.rst @@ -16,3 +16,9 @@ To bridge this gap seamlessly, the fluid model supports an automated conversion 1. Entry Conversion: Upon entering the fluid property update routine, if the overall fluid state is provided in mass fractions, the framework dynamically converts these inputs into molar fractions using the designated component molecular weights. 2. Internal Computation: All stability, flash, density, and viscosity calculations proceed purely in molar space. 3. Exit Transformation: Before the update routine concludes, the resulting phase compositions, phase fractions, and molar densities are transformed back into mass variables. Additionally, rigorous chain-rule transformations are applied to all derivatives (with respect to pressure, temperature, and composition) to ensure the Jacobian matrices populated by the simulator remain mathematically consistent with the mass-based primary variables. + +Note on Component Limits +-------------------------- + +When compositional fluid models are evaluated in standalone mode (for example, using the ``PVTDriver`` task (:ref:`PVTDriver`)), the maximum number of allowed components is 9. However, when these fluid models are actively coupled with a flow solver during a full simulation, the maximum number of components is restricted to 5. + diff --git a/src/coreComponents/constitutiveDrivers/docs/PVTDriver.rst b/src/coreComponents/constitutiveDrivers/docs/PVTDriver.rst index 322f2e32fae..7546fffa1d4 100644 --- a/src/coreComponents/constitutiveDrivers/docs/PVTDriver.rst +++ b/src/coreComponents/constitutiveDrivers/docs/PVTDriver.rst @@ -1,3 +1,5 @@ +.. _PVTDriver: + PVT Driver =============== From 749c229449b671c1804131f56f8a9203141d61db Mon Sep 17 00:00:00 2001 From: Dickson Kachuma Date: Thu, 4 Jun 2026 15:57:53 -0500 Subject: [PATCH 13/13] Upload missing file --- .../docs/compositional/Models.rst | 57 +++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 src/coreComponents/constitutive/docs/compositional/Models.rst diff --git a/src/coreComponents/constitutive/docs/compositional/Models.rst b/src/coreComponents/constitutive/docs/compositional/Models.rst new file mode 100644 index 00000000000..a7e01ac1c29 --- /dev/null +++ b/src/coreComponents/constitutive/docs/compositional/Models.rst @@ -0,0 +1,57 @@ +List of available compositional fluid models +============================================= + +The following table outlines the available compositional fluid models (specified via their XML elements), detailing the number of phases they support, their underlying thermodynamic phase split (flash) methods, the phase density and viscosity models used, and their recommended applications: + +.. list-table:: + :widths: 100 + :header-rows: 1 + + * - Compositional fluid models + * - ``CompositionalTwoPhaseFluid`` :ref:`XML_CompositionalTwoPhaseFluid` + + * Number of phases: 2. + * Possible applications: Idealised simulations, synthetic testing, or when fluid flow is dominated entirely by pressure gradients. + * Phase split method: Negative Two-Phase Flash. + * Phase density models: Compositional density. + * Phase viscosity models: Constant viscosity. + + * - ``CompositionalTwoPhaseFluidLohrenzBrayClark`` :ref:`XML_CompositionalTwoPhaseFluidLohrenzBrayClark` + + * Number of phases: 2. + * Possible applications: Miscible gas injection, primary depletion of volatile oils, and scenarios where fluid compositions change drastically. + * Phase split method: Negative Two-Phase Flash. + * Phase density models: Compositional density. + * Phase viscosity models: Lohrenz-Bray-Clark viscosity. + + * - ``CompositionalTwoPhaseFluidPhillipsBrine`` :ref:`XML_CompositionalTwoPhaseFluidPhillipsBrine` + + * Number of phases: 2. + * Possible applications: Carbon sequestration (CO2-brine systems) or water-flooding in oil reservoirs. + * Phase split method: Negative Two-Phase Flash. + * Phase density models: Phillips brine density model for the aqueous phase; Compositional density for the gas phase. + * Phase viscosity models: Phillips brine viscosity model for the aqueous phase; Lohrenz-Bray-Clark viscosity for the gas phase. + + * - ``CompositionalThreePhaseFluidLohrenzBrayClark`` :ref:`XML_CompositionalThreePhaseFluidLohrenzBrayClark` + + * Number of phases: 3. + * Possible applications: Multi-phase reservoir simulations where water solubility in the hydrocarbon phases is negligible. + * Phase split method: Immiscible Water Flash. + * Phase density models: Immiscible water density for the aqueous phase; Compositional density for the hydrocarbon phases. + * Phase viscosity models: Immiscible water viscosity for the aqueous phase; Lohrenz-Bray-Clark viscosity for the hydrocarbon phases. + + * - ``CompositionalTwoPhaseKValueFluidLohrenzBrayClark`` :ref:`XML_CompositionalTwoPhaseKValueFluidLohrenzBrayClark` + + * Number of phases: 2. + * Possible applications: Thermal recovery processes or relatively stable black-oil systems. + * Phase split method: K-Value Flash. + * Phase density models: Compositional density. + * Phase viscosity models: Lohrenz-Bray-Clark viscosity. + + * - ``CompositionalTwoPhaseKValueFluidPhillipsBrine`` :ref:`XML_CompositionalTwoPhaseKValueFluidPhillipsBrine` + + * Number of phases: 2. + * Possible applications: Saline aquifers and CO2 sequestration with stable compositions where pre-tabulated K-values are sufficient. + * Phase split method: K-Value Flash. + * Phase density models: Phillips brine density model; Compositional density. + * Phase viscosity models: Phillips brine viscosity model; Lohrenz-Bray-Clark viscosity. \ No newline at end of file