flash.cpp

Go to the documentation of this file.

// C++ includes
#include <cmath>
#include <iostream>
#include <iomanip>

// FlashCCS includes
#include "flash.h"
#include "flashessential.h"

namespace Flash
{
// the nonlinear equation for the computation of the (theta)-equation
struct ThetaEquation
{
        double zNaCl, zH2O, A, a, b, c;
        double operator()(double theta) { return (A * theta/(zH2O * theta + zNaCl * (A - 1)) - theta - 1) * std::exp(b * theta + c * theta*theta) - a; }
};

// the nonlinear equation for the computation of the derivative of the (theta)-equation
struct ThetaDerivativeEquation
{
        double zNaCl, zH2O, A, a, b, c;
        double operator()(double theta)
        {
                // the denominator (a common factor that appears in almost each term)
                double den = zH2O * theta + zNaCl * (A - 1);

                return expf(b * theta + c * theta * theta) *
                                ( (A/den - A * theta * zH2O/(den * den) - 1) + (A * theta/den - theta - 1) * (b + 2 * c * theta) );
        }
};

// it computes the nonlinear (theta)-equation
double ComputeTheta(double zNaCl, double zH2O, double A, double a, double b, double c)
{
        // start value for the nonlinear solver
        double theta0 = zNaCl/zH2O;
        // the nonlinear equation
        ThetaEquation function;
        // the derivative of the nonlinear function representing the equation
        ThetaDerivativeEquation derivative;
        // setting the parameters of the equation and its derivative
        function.zNaCl = zNaCl; derivative.zNaCl = zNaCl;
        function.zH2O = zH2O;   derivative.zH2O = zH2O;
        function.A = A;         derivative.A = A;
        function.a = a;         derivative.a = a;
        function.b = b;         derivative.b = b;
        function.c = c;         derivative.c = c;
        // solving the nonlinear equation
        return NewtonRaphson(function, derivative, theta0);
}

void FlashBrineCO2H2S(BrineCO2H2SData& data, const BrineCO2Data& co2BrineData, double T, double P, double zH2O, double zCO2, double zH2S)
{
        // determine the state of the CO2 and H2S
        data.stateCO2 = co2BrineData.stateCO2;
        data.stateH2S = DeterminePhysicalStateH2S(T, P);

        // auxiliary variable to determine if the H2S is fully solubilized in the aqueous phase
        bool isH2SFullySolubilized;

        // verify the existence of an aqueous phase
        bool isThereAqueousPhase = (co2BrineData.betaA != 0.0);

        // compute the molar fractions in the aqueous phase
        if(isThereAqueousPhase)
        {
                // compute the NaCl molality in the aqueous phase
                double mNaCl = 55.508 * co2BrineData.xNaCl/co2BrineData.xH2O;

                // compute the CO2 molality in the aqueous phase
                double mCO2 = 55.508 * co2BrineData.xCO2/co2BrineData.xH2O;

                // compute the saturated H2S molality
                double saturatedMolalityH2S = (P - ComputeSaturatedPressureH2O(T))/std::exp(uRTH2S(T, P) - std::log(ComputeFugacityH2S(ComputeVolumeH2S(T, P, data.stateH2S), T, P)) + 2 * lambdaH2S(T, P) * mNaCl + zetaH2S * mNaCl * mNaCl);

                // compute the H2S molality assuming fully solubilization
                double mH2SFullySolubilized = 55.508 * zH2S/zH2O;

                // auxiliary variable for the H2S molality in the aqueous phase
                double mH2S;

                // compute the H2S molality in the aqueous phase and also determine if the H2S is fully solubilized
                if(mH2SFullySolubilized < saturatedMolalityH2S)
                {
                        mH2S = mH2SFullySolubilized;
                        isH2SFullySolubilized = true;
                }
                else
                {
                        mH2S = saturatedMolalityH2S;
                        isH2SFullySolubilized = false;
                }

                // compute the molar fractions in the aqueous phase
                data.xH2O  = 55.508/(55.508 + mNaCl + mCO2 + mH2S);
                data.xNaCl = data.xH2O * mNaCl/55.508;
                data.xCO2  = data.xH2O * mCO2/55.508;
                data.xH2S  = 1 - data.xNaCl - data.xH2O - data.xCO2;
        }
        else
        {
                // the molar fractions in the aqueous phase are all zero
                data.xH2O  = 0.0;
                data.xNaCl = 0.0;
                data.xCO2  = 0.0;
                data.xH2S  = 0.0;

                // verify if H2S is fully solubilized in the aqueous phase
                isH2SFullySolubilized = false;
        }

        // verify the existence of a CO2-rich phase
        bool isThereCO2RichPhase = (co2BrineData.betaC != 0.0);

        // verify the existence of a H2S-rich phase (CO2 or H2S has to be in the liquid state and H2S cannot be fully solubilized in the aqueous phase)
        bool isThereH2SRichPhase = (data.stateCO2 == LIQUID || data.stateH2S == LIQUID) && !isH2SFullySolubilized;

        // compute the phase molar fractions
        data.betaS = co2BrineData.betaS * (1.0 - zH2S);
        data.betaA = co2BrineData.betaA * (1.0 - zH2S)/(1 - data.xH2S);

        if(isThereH2SRichPhase)
        {
                data.betaC = co2BrineData.betaC * (1.0 - zH2S);
                data.betaH = 1.0 - data.betaS - data.betaA - data.betaC;
        }
        else
        {
                data.betaC = 1.0 - data.betaS - data.betaA;
                data.betaH = 0.0;
        }

        // compute the molar fractions in the CO2-rich phase
        if(isThereCO2RichPhase)
        {
                data.yH2S = (isThereH2SRichPhase) ? 0.0 : (zH2S - data.xH2S * data.betaA)/data.betaC;
                data.yCO2 = (zCO2 - data.xCO2 * data.betaA)/data.betaC;
                data.yH2O = 1 - data.yCO2 - data.yH2S;
        }
        else
        {
                data.yH2S = 0.0;
                data.yCO2 = 0.0;
                data.yH2O = 0.0;
        }
}

void ProgressiveFlashBrineCO2(BrineCO2Data& data, double T, double P, double zNaCl, double zH2O, double zCO2, unsigned int numCorrections)
{
        // determine the physical state of CO2
        data.stateCO2 = DeterminePhysicalStateCO2(T, P);

        // compute the volume of CO2
        double v = ComputeVolumeCO2(T, P, data.stateCO2);

        // calculate the constants A and B Spycher et al. (2003)
        double A = koH2O(T)/(ComputeFugacityH2O(v, T, P) * P) * exp((P - Po)/(R * T) * vH2O);
        double B = ComputeFugacityCO2(v, T, P) * P/(koCO2(T, data.stateCO2) * 55.508) * exp((Po - P)/(R * T) * vCO2(data.stateCO2));

        // calculate the parameters (lambda) and (zeta) from the activity model of Duan and Sun (2003)
        double lambda = lambdaCO2(T, P);
        double zeta   = zetaCO2(T, P);

        // calculate the saturated solubility of NaCl in water Sawamura et al. (2007)
        double mNaClSaturated = ComputeSaturatedMolalityNaCl(T, P);

        // calculate the salt molality if all salt would be solubilized in the aqueous phase
        double mNaClFullySolubilized = 55.508 * zNaCl/zH2O;

        // calculate the CO2 molality if all CO2 would be solubilized in the aqueous phase
        double mCO2FullySolubilized = 55.508 * zCO2/zH2O;

        // auxiliary variables to determine if there exists a solid phase, aqueous phase and CO2-rich phase
        bool isThereSolidPhase, isThereAqueousPhase, isThereCO2RichPhase;

        // the salt and CO2 molality in the aqueous phase
        double mNaCl, mCO2;

        // calculate the salt molality in the aqueous phase by assuming that the water transference to the CO2-rich phase was sufficiently small
        if(mNaClFullySolubilized < mNaClSaturated)
        {
                mNaCl = mNaClFullySolubilized;
                isThereSolidPhase = false;
        }
        else
        {
                mNaCl = mNaClSaturated;
                isThereSolidPhase = true;
        }

        // calculate a common factor that appears for each iteration in the next loop
        const double commonFactor = 55.508 * (1.0 - A)/(1.0 - B) * B;

        // for each iteration, determine the equilibrium phases and calculate the molar fractions
        for(unsigned int counter = 0; counter < numCorrections; ++counter)
        {
                // calculate the saturated CO2 solubility in brine, Spycher et al. (2005).
                double mCO2Saturated = commonFactor * std::exp(- mNaCl * (2.0 * lambda + zeta * mNaCl));

                // calculate the CO2 molality in the aqueous phase with the current salt molality (mNaCl)
                if(mCO2FullySolubilized < mCO2Saturated)
                {
                        mCO2 = mCO2FullySolubilized;
                        isThereCO2RichPhase = false;
                }
                else
                {
                        mCO2 = mCO2Saturated;
                        isThereCO2RichPhase = true;
                }

                // determine if there exists an aqueous phase only if a CO2-rich phase exists
                if(isThereCO2RichPhase)
                {
                        // assume an aqueous phase exists and compute its molar fractions
                        double xH2O = 55.508/(55.508 + mNaCl + mCO2);
                        double xCO2 = mCO2/(55.508 + mNaCl + mCO2);

                        // compute the molar fractions in the CO2-rich phase
                        double yH2O = A * xH2O;
                        double yCO2 = 1.0 - yH2O;

                        // compute betaA and check if it is greater than zero, which implies the existence of the aqueous phase
                        isThereAqueousPhase = ((yCO2 * zH2O - yH2O * zCO2)/(yCO2 * xH2O - yH2O * xCO2) > 0.0) ? true : false;
                }
                else
                {
                        // an aqueous phase exists whenever a CO2-rich phase is not present
                        isThereAqueousPhase = true;
                }

                // (AC) phase equilibrium configuration
                if(isThereAqueousPhase && isThereCO2RichPhase && !isThereSolidPhase)
                {
                        data.xNaCl = mNaCl/(55.508 + mNaCl + mCO2);
                        data.xCO2  = mCO2/(55.508 + mNaCl + mCO2);
                        data.xH2O  = 1.0 - data.xNaCl - data.xCO2;

                        data.yH2O = A * data.xH2O;
                        data.yCO2 = 1.0 - data.yH2O;

                        data.betaA = (data.yCO2 * zH2O - data.yH2O * zCO2)/(data.yCO2 * data.xH2O - data.yH2O * data.xCO2);
                        data.betaC = 1.0 - data.betaA;
                        data.betaS = 0.0;

                        // correct the concentration of NaCl in the aqueous phase due to the amount of transfered water to the CO2-rich phase if we have the AC phase equilibrium configuration.
                        double mNaClFullySolubilizedCorrected = mNaClFullySolubilized * (1.0 + A * data.betaC/data.betaA);

                        if(mNaClFullySolubilizedCorrected < mNaClSaturated)
                        {
                                mNaCl = mNaClFullySolubilizedCorrected;
                                isThereSolidPhase = false;
                        }
                        else
                        {
                                mNaCl = mNaClSaturated;
                                isThereSolidPhase = true;
                        }
                }
                // (SAC) phase equilibrium configuration
                else if(isThereAqueousPhase && isThereCO2RichPhase && isThereSolidPhase)
                {
                        data.xNaCl = mNaCl/(55.508 + mNaCl + mCO2);
                        data.xCO2  = mCO2/(55.508 + mNaCl + mCO2);
                        data.xH2O  = 1.0 - data.xNaCl - data.xCO2;

                        data.yH2O = A * data.xH2O;
                        data.yCO2 = 1.0 - data.yH2O;

                        data.betaA = (data.yCO2 * zH2O - data.yH2O * zCO2)/(data.yCO2 * data.xH2O - data.yH2O * data.xCO2);
                        data.betaC = (data.xCO2 * zH2O - data.xH2O * zCO2)/(data.xCO2 * data.yH2O - data.xH2O * data.yCO2);
                        data.betaS = 1.0 - data.betaA - data.betaC;

                        break; // abandon the loop, since this phase equilibrium configuration does not require corrections
                }
                // (SA) phase equilibrium configuration
                else if(isThereAqueousPhase && isThereSolidPhase && !isThereCO2RichPhase)
                {
                        data.xNaCl = mNaCl/(55.508 + mNaCl + mCO2);
                        data.xCO2  = mCO2/(55.508 + mNaCl + mCO2);
                        data.xH2O  = 1.0 - data.xNaCl - data.xCO2;

                        data.yH2O = 0.0;
                        data.yCO2 = 0.0;

                        data.betaA = zH2O/data.xH2O;
                        data.betaS = 1.0 - data.betaA;
                        data.betaC = 0.0;

                        break; // abandon the loop, since this phase equilibrium configuration does not require corrections
                }
                // (SC) phase equilibrium configuration
                else if(isThereCO2RichPhase &&  !isThereAqueousPhase)
                {
                        data.xNaCl = 0.0;
                        data.xCO2  = 0.0;
                        data.xH2O  = 0.0;

                        data.yH2O = zH2O/(1.0 - zNaCl);
                        data.yCO2 = 1.0 - data.yH2O;

                        data.betaS = zNaCl;
                        data.betaC = 1.0 - zNaCl;
                        data.betaA = 0.0;

                        break; // abandon the loop, since this phase equilibrium configuration does not require corrections
                }
                // (A) phase equilibrium configuration
                else
                {
                        data.xNaCl = zNaCl;
                        data.xCO2  = zCO2;
                        data.xH2O  = zH2O;

                        data.yH2O = 0.0;
                        data.yCO2 = 0.0;

                        data.betaA = 1.0;
                        data.betaS = 0.0;
                        data.betaC = 0.0;

                        break; // abandon the loop, since this phase equilibrium configuration does not require corrections
                }
        }
}

void GeometricFlashBrineCO2(BrineCO2Data& data, double T, double P, double zNaCl, double zH2O, double zCO2)
{
        // determine the physical state of CO2.
        data.stateCO2 = DeterminePhysicalStateCO2(T, P);

        // compute the volume of CO2.
        double v = ComputeVolumeCO2(T, P, data.stateCO2);

        // calculate the constants A and B Spycher et al. (2003).
        double A = koH2O(T)/(ComputeFugacityH2O(v, T, P) * P) * std::exp((P - Po)/(R * T) * vH2O);
        double B = ComputeFugacityCO2(v, T, P) * P/(koCO2(T, data.stateCO2) * 55.508) * std::exp((Po - P)/(R * T) * vCO2(data.stateCO2));

        // calculate the parameters (lambda) and (zeta) from the activity model of Duan and Sun (2003).
        double lambda = lambdaCO2(T, P);
        double zeta   = zetaCO2(T, P);

        // calculate the saturated solubility of NaCl in water Sawamura et al. (2007).
        double mNaClSaturated = ComputeSaturatedMolalityNaCl(T, P);

        // calculate the saturated solubility of CO2 in the brine assuming that it is saturated with salt.
        double mCO2SaturatedBrine = 55.508 * (1.0 - A)/(1.0 - B) * B * std::exp(- mNaClSaturated * (2.0 * lambda + zeta * mNaClSaturated));

        // calculate the molar fractions xNaCl*, xCO2* and yH2O*
        double xNaClStar = mNaClSaturated/(55.508 + mNaClSaturated + mCO2SaturatedBrine);
        double xCO2Star  = mCO2SaturatedBrine/(55.508 + mNaClSaturated + mCO2SaturatedBrine);
        double yH2OStar  = A * (1 - xNaClStar - xCO2Star);

        // calculate the vectors s*, x* and y* and define the vector z = (zNaCl, zH2O, zCO2)
        double s[3] = {1.0, 0.0, 0.0};
        double x[3] = {xNaClStar, 1 - xNaClStar - xCO2Star, xCO2Star};
        double y[3] = {0.0, yH2OStar, 1 - yH2OStar};
        double z[3] = {zNaCl, zH2O, zCO2};

        // determine the equilibrium phases
        if(((x[0] - y[0]) * (z[1] - y[1]) - (x[1] - y[1]) * (z[0] - y[0])) > 0) // aqueous phase (L) or aqueous phase and CO2-rich phase(LC)
        {
                // calculate mNaClA and mCO2A (see documentation)
                double mNaClA = 55.508 * zNaCl/zH2O;
                double mCO2A  = 55.508 * zCO2/zH2O;

                // calculate the saturated solubility of CO2 with mNaClA
                double mCO2SaturatedA = 55.508 * (1.0 - A)/(1.0 - B) * B * std::exp(- mNaClA * (2.0 * lambda + zeta * mNaClA));

                // aqueous phase configuration (A)
                if(mCO2A < mCO2SaturatedA)
                {
                        data.xNaCl = zNaCl;
                        data.xH2O  = zH2O;
                        data.xCO2  = zCO2;
                        data.yH2O  = 0;
                        data.yCO2  = 0;
                        data.betaS = 0;
                        data.betaA = 1;
                        data.betaC = 0;
                }
                // aqueous phase and CO2-rich phase configuration (AC)
                else
                {
                        // if zNaCl <= 0.0001, we assume the system does not contain salt and use the solubility model for CO2 in pure water
                        if(zNaCl <= 0.0001)
                        {
                                data.xH2O  = (1.0 - B)/(1.0 - A * B);
                                data.xCO2  = 1.0 - data.xH2O;
                                data.xNaCl = 0.0;
                        }
                        else
                        {
                                // calculate the constants a, b and c
                                double a = (1.0 - A)/(1.0 - B) * B;
                                double b = 111.016 * lambda;
                                double c = 3081.138064 * zeta;

                                // calculate (theta)
                                double theta = ComputeTheta(zNaCl, zH2O, A, a, b, c);

                                data.xH2O  = (zNaCl * (A - 1) + zH2O * theta)/(A * theta);
                                data.xNaCl = data.xH2O * theta;
                                data.xCO2  = 1 - data.xNaCl - data.xH2O;
                        }

                        // compute the remaining molar fractions
                        data.yH2O  = A * data.xH2O;
                        data.yCO2  = 1 - data.yH2O;
                        data.betaS = 0;
                        data.betaA = (data.yCO2 * zH2O - data.yH2O * zCO2)/(data.yCO2 * data.xH2O - data.yH2O * data.xCO2);
                        data.betaC = 1 - data.betaA;
                }
        }
        // solid phase and aqueous phase configuration (SA)
        else if(((s[0] - x[0]) * (z[1] - x[1]) - (s[1] - x[1]) * (z[0] - x[0])) > 0)
        {
                data.xH2O  = 1.0/(mNaClSaturated/55.508 + zCO2/zH2O + 1);
                data.xNaCl = data.xH2O * mNaClSaturated/55.508;
                data.xCO2  = 1 - data.xNaCl - data.xH2O;
                data.yH2O  = 0;
                data.yCO2  = 0;
                data.betaS = 1 - zH2O/data.xH2O;
                data.betaA = 1 - data.betaS;
                data.betaC = 0;
        }
        // solid phase and CO2-rich phase configuration (SC)
        else if(((y[0] - s[0]) * (z[1] - s[1]) - (y[1] - s[1]) * (z[0] - s[0])) > 0)
        {
                data.xH2O  = 0;
                data.xNaCl = 0;
                data.xCO2  = 0;
                data.yH2O  = zH2O/(zH2O + zCO2);
                data.yCO2  = 1 - data.yH2O;
                data.betaS = zNaCl;
                data.betaA = 0;
                data.betaC = 1 - data.betaS;
        }
        // solid phase, aqueous phase and CO2-rich phase configuration (SAC)
        else
        {
                data.xNaCl = xNaClStar;
                data.xCO2  = xCO2Star;
                data.xH2O  = 1 - xNaClStar - xCO2Star;
                data.yH2O  = yH2OStar;
                data.yCO2  = 1 - yH2OStar;
                data.betaA = (data.yCO2 * zH2O - data.yH2O * zCO2)/(data.yCO2 * data.xH2O - data.yH2O * data.xCO2);
                data.betaC = (data.xCO2 * zH2O - data.xH2O * zCO2)/(data.xCO2 * data.yH2O - data.xH2O * data.yCO2);
                data.betaS = 1 - data.betaA - data.betaC;
        }
}

void ProgressiveFlashBrineCO2H2S(BrineCO2H2SData& data, double T, double P, double zNaCl, double zH2O, double zCO2, double zH2S, unsigned int numCorrections)
{
        // compute the flash for the CO2-Brine data alone
        BrineCO2Data co2BrineData;

        ProgressiveFlashBrineCO2(co2BrineData, T, P, zNaCl, zH2O, zCO2, numCorrections);

        FlashBrineCO2H2S(data, co2BrineData, T, P, zH2O, zCO2, zH2S);
}

void GeometricFlashBrineCO2H2S(BrineCO2H2SData& data, double T, double P, double zNaCl, double zH2O, double zCO2, double zH2S)
{
        // compute the flash for the CO2-Brine system alone
        BrineCO2Data co2BrineData;

        GeometricFlashBrineCO2(co2BrineData, T, P, zNaCl, zH2O, zCO2);

        FlashBrineCO2H2S(data, co2BrineData, T, P, zH2O, zCO2, zH2S);
}

std::ostream& operator<<(std::ostream& out, BrineCO2Data& data)
{
        out << "+--------------------------+" << std::endl;
        out << "+---- Brine-CO2 Phases ----+" << std::endl;
        out << "+--------------------------+" << std::endl << std::endl;

        out << std::setprecision(6) << std::fixed << std::showpoint;

        if(data.betaS > 0.0)
        {
                out << ":: Solid Phase, " << data.betaS * 100 << "%" << std::endl << std::endl;

                out << "sNaCl = " << 100 << "%" << std::endl << std::endl;
        }

        if(data.betaA > 0.0)
        {
                out << ":: Aqueous Phase, " << data.betaA * 100 << "%" << std::endl << std::endl;

                out << "xNaCl  = " << data.xNaCl * 100 << "%" << std::endl;
                out << "xH2O   = " << data.xH2O * 100 << "%" << std::endl;
                out << "xCO2   = " << data.xCO2 * 100 << "%" << std::endl;
                out << "mNaCl  = " << 55.508 * data.xNaCl / data.xH2O << std::endl;
                out << "mCO2   = " << 55.508 * data.xCO2 / data.xH2O << std::endl << std::endl;
        }

        if(data.betaC > 0.0)
        {
                out << ((data.stateCO2 == LIQUID) ? ":: Liquid CO2-Rich Phase, " : (data.stateCO2 == GAS) ? ":: Gas CO2-Rich Phase, " : ":: Supercritical CO2-Rich Phase, ") << data.betaC * 100 << "%" << std::endl << std::endl;

                out << "yH2O  = " << data.yH2O * 100 << "%" << std::endl;
                out << "yCO2  = " << data.yCO2 * 100 << "%" << std::endl << std::endl;
        }

        return out;
}

std::ostream& operator<<(std::ostream& out, BrineCO2H2SData& data)
{
        out << "+------------------------------+" << std::endl;
        out << "+---- Brine-CO2-H2S Phases ----+" << std::endl;
        out << "+------------------------------+" << std::endl << std::endl;

        out << std::setprecision(6) << std::fixed << std::showpoint;

        if(data.betaS > 0.0)
        {
                out << ":: Solid Phase, " << data.betaS * 100 << "%" << std::endl << std::endl;

                out << "sNaCl = " << 100 << "%" << std::endl << std::endl;
        }

        if(data.betaA > 0.0)
        {
                out << ":: Aqueous Phase, " << data.betaA * 100 << "%" << std::endl << std::endl;

                out << "xNaCl  = " << data.xNaCl * 100 << "%" << std::endl;
                out << "xH2O   = " << data.xH2O * 100 << "%" << std::endl;
                out << "xCO2   = " << data.xCO2 * 100 << "%" << std::endl;
                out << "xH2S   = " << data.xH2S * 100 << "%" << std::endl << std::endl;
                out << "mNaCl  = " << 55.508 * data.xNaCl / data.xH2O << std::endl;
                out << "mCO2   = " << 55.508 * data.xCO2 / data.xH2O << std::endl;
                out << "mH2S   = " << 55.508 * data.xH2S / data.xH2O << std::endl << std::endl;
        }

        if(data.betaC > 0.0)
        {
                out << ((data.stateCO2 == LIQUID) ? ":: Liquid CO2-Rich Phase, " : (data.stateCO2 == GAS) ? ":: Gas CO2-Rich Phase, " : ":: Supercritical CO2-Rich Phase, ") << data.betaC * 100 << "%" << std::endl << std::endl;

                out << "yH2O  = " << data.yH2O * 100 << "%" << std::endl;
                out << "yCO2  = " << data.yCO2 * 100 << "%" << std::endl;
                out << "yH2S  = " << data.yH2S * 100 << "%" << std::endl << std::endl;
        }

        if(data.betaH > 0.0)
        {
                out << ((data.stateH2S == LIQUID) ? ":: Liquid H2S-Rich Phase, " : (data.stateH2S == GAS) ? ":: Gas H2S-Rich Phase, " : ":: Supercritical H2S-Rich Phase, ") << data.betaH * 100 << "%" << std::endl << std::endl;

                out << "hH2S  = " << 100 << "%" << std::endl;
        }

        return out;
}
} // namespace Flash

---