flashessential.h

Go to the documentation of this file.

#ifndef FLASHESSENTIAL_H
#define FLASHESSENTIAL_H

// C++ includes
#include <cmath>

// FlashCCS includes
#include "newtonraphson.h"
#include <assert.h>
namespace Flash
{
//==================================================+
//--------------------------------------------------+
// Possible Physical States of CO2 and H2S
//--------------------------------------------------+
//==================================================+
enum PhysicalState { LIQUID, GAS, SUPERCRITICAL };

//==================================================+
//--------------------------------------------------+
// Critical Temperatures and Pressures
//--------------------------------------------------+
//==================================================+
// the critical temperature and pressure of H2O
const double TcH2O = 647.30;
const double PcH2O = 221.20;

// the critical temperature and pressure of CO2
const double TcCO2 = 304.20;
const double PcCO2 =  73.83;

// the critical temperature and pressure of H2S
const double TcH2S = 373.50;
const double PcH2S =  89.63;

//==================================================+
//--------------------------------------------------+
//  Universal Gas Constant (R)
//--------------------------------------------------+
//==================================================+
const double R = 83.14472; // [R] = bar.cm^3/(mol.K)

//==================================================+
//--------------------------------------------------+
// Parameters and Constants for Solubility Models
//--------------------------------------------------+
//==================================================+
// the parameters for the Redlich-Kwong EOS Spycher et al. (2003)
inline double aCO2(double T) { return 7.54E7 - 4.13E4 * T; };
 const double bCO2    = 27.80;
 const double bH2O    = 18.18;
 const double aH2OCO2 = 7.89E7;

// the parameters for the EOS of Duan and Sun (2007) for H2S
const double a1H2S  =  5.2386075E-2;
const double a2H2S  = -2.7463906E-1;
const double a3H2S  = -9.6760173E-2;
const double a4H2S  =  1.3618104E-2;
const double a5H2S  = -8.8681753E-2;
const double a6H2S  =  4.1176908E-2;
const double a7H2S  =  3.6354018E-4;
const double a8H2S  =  2.2719194E-3;
const double a9H2S  = -7.6962514E-4;
const double a10H2S = -2.1948579E-5;
const double a11H2S = -1.1707631E-4;
const double a12H2S =  4.0756926E-5;
const double a13H2S =  5.7582260E-2;
const double a14H2S =  1.00;
const double a15H2S =  0.06;

// the equilibrium parameters Spycher et al. (2003)
inline double koH2O(double T)
{ 
        T = T - 273.15;
        return std::pow(10.0, -2.209 + 3.097E-2 * T - 1.098E-4 * T*T + 2.048E-7 * T*T*T); 
}

inline double koCO2(double T, PhysicalState stateCO2) 
{ 
        T = T - 273.15; 
        return (stateCO2 == GAS || stateCO2 == SUPERCRITICAL) ? 
                std::pow(10.0, 1.189 + 1.304E-2 * T - 5.446E-5 * T*T) : 
                std::pow(10.0, 1.169 + 1.368E-2 * T - 5.380E-5 * T*T); 
}

// the reference pressure Spycher et al. (2003)
const double Po = 1.0;

// the partial average molar volume Spycher et al. (2003)
const  double vH2O = 18.1;
inline double vCO2(PhysicalState stateCO2) { return (stateCO2 == GAS || stateCO2 == SUPERCRITICAL) ? 32.6 : 32.0; }

// the parameters (lambda)CO2 and (zeta)CO2 of the activity model of Duan and Sun (2003)
inline double lambdaCO2(double T, double P) 
{ 
        return -0.411370585 + 6.07632013E-4 * T + 97.5347708/T - 0.0237622469 * P/T + 0.0170656236 * P/(630.0 - T) + 1.41335834E-5 * T * std::log(P); 
}

inline double zetaCO2(double T, double P)
{ 
        return 3.36389723E-4 - 1.98298980E-5 * T + 2.12220830E-3 * P/T - 5.24873303E-3 * P/(630.0 - T);
}

// the parameters (mu)oH2S/RT, (lambda)H2S and (zeta)H2S of the activity model of Duan and Sun (2007)
inline double uRTH2S(double T, double P) 
{ 
        return 42.564957 - 8.6260377E-2 * T - 6084.3775/T + 6.8714437E-5 * (T*T) - 102.76849/(680 - T) + 8.4482895E-4 * P - 1.0590768 * P/(680 - T) + 3.5665902E-3 * (P*P)/T; 
}

inline double lambdaH2S(double T, double P) 
{ 
        return 8.5004999E-2 + 3.5330378E-5 * T - 1.5882605/T + 1.1894926E-5 * P;
}

const  double zetaH2S = -1.0832589E-2;

//==================================================+
//--------------------------------------------------+
// Saturated NaCl Molality Method
//--------------------------------------------------+
//==================================================+
// it computes the saturated solubility of salt in water Sawamura et al. (2007)
inline double ComputeSaturatedMolalityNaCl(double T, double P) { return std::exp( 2.86623 - 571.361/T + 77128/(T*T) + 4.0479E-4 * P - 7.445E-7 * (P*P) + 6.209E-10 * (P*P*P) ); }

//==================================================+
//--------------------------------------------------+
// Saturated Pressure Methods
//--------------------------------------------------+
//==================================================+
// it computes the vapor pressure of H2O (The properties of gases and liquids, 5th ed., McGraw-Hill, 2001)
inline double ComputeSaturatedPressureH2O(double T) 
{
        double x = 1.0 - T/TcH2O;
        return PcH2O * std::exp(1.0/(1 - x) * (-7.76451 * x + 1.45838 * std::pow(x, 1.5) - 2.77580 * std::pow(x, 3) - 1.23303 * std::pow(x, 6))); 
}

// it computes the vapor pressure of CO2 (The properties of gases and liquids, 5th ed., McGraw-Hill, 2001)
inline double ComputeSaturatedPressureCO2(double T) 
{
        double x = 1.0 - T/TcCO2;
        return PcCO2 * std::exp(1.0/(1 - x) * (-6.95626 * x + 1.19695 * std::pow(x, 1.5) - 3.12614 * std::pow(x, 3) + 2.99448 * std::pow(x, 6))); 
}

// the nonlinear equation for the computation of the vapor pressure of H2S
struct SaturatedPressureH2SEquation
{
        double T;
        double operator()(double Psat) { return 36.067 - 3132.31/T - 3.985 * std::log(T) + 653.0 * Psat/(T*T) - std::log(Psat); }       
};

// it computes the vapor pressure of H2S (The properties of gases and liquids, 5th ed., McGraw-Hill, 2001)
inline double ComputeSaturatedPressureH2S(double T)
{
        // start value for the nonlinear solver
        double Psat0 = 20.0; // vapor pressure near 300K
        // the nonlinear equation
        static SaturatedPressureH2SEquation equation;
        // setting the parameters of the equation
        equation.T = T;
        // solving the nonlinear equation
        return NewtonRaphson(equation, Psat0);
}

//==================================================+
//--------------------------------------------------+
// Physical State Determination Methods
//--------------------------------------------------+
//==================================================+
// it determines the physical state of CO2
inline PhysicalState DeterminePhysicalStateCO2(double T, double P)
{
        return (T < TcCO2) ? ( (P > ComputeSaturatedPressureCO2(T)) ? LIQUID : GAS ) : ( (P > PcCO2) ? SUPERCRITICAL : GAS );
}

// it determines the physical state of H2S
inline PhysicalState DeterminePhysicalStateH2S(double T, double P)
{
        return (T < TcH2S) ? ( (P > ComputeSaturatedPressureH2S(T)) ? LIQUID : GAS ) : ( (P > PcH2S) ? SUPERCRITICAL : GAS );
}

//==================================================+
//--------------------------------------------------+
// Volume Methods
//--------------------------------------------------+
//==================================================+
// the nonlinear equation for the computation of the volume of CO2, Spycher et al. (2003)
struct VolumeCO2Equation
{
        double T, P;
        double operator()(double V) { return V*V*V - (R*T/P)*(V*V) - (R*T*bCO2/P - aCO2(T)/(P*std::sqrt(T)) + bCO2*bCO2)*(V) - aCO2(T)*bCO2/(P*std::sqrt(T)); }
};

// the derivative of the nonlinear equation for the computation of the volume of CO2
struct VolumeCO2DerivativeEquation
{
        double T, P;
        double operator()(double V) { return 3*V*V - 2*(R*T/P)*V - (R*T*bCO2/P - aCO2(T)/(P*sqrtf(T)) + bCO2*bCO2); }
};

// it computes the volume of CO2, Spycher et al. (2003)
inline double ComputeVolumeCO2(double T, double P, PhysicalState stateCO2)
{
        // start value for the nonlinear solver
        double V0 = (stateCO2 == GAS || stateCO2 == SUPERCRITICAL) ? (R * T/P) : 59.22; // vCO2 = 59.22 mol/cm^3 at 25 °C and 80.0 bar

        // the nonlinear equation
        static VolumeCO2Equation equation;
        // the derivative of the nonlinear equation
        static VolumeCO2DerivativeEquation derivative;
        // setting the parameters of the equation
        equation.T = T; derivative.T = T;
        equation.P = P; derivative.P = P;

        // solving the nonlinear equation
        return NewtonRaphson(equation, derivative, V0);
}


 inline double ComputeVolumeCO2(double T, double P)
{
  PhysicalState stateCO2=DeterminePhysicalStateCO2(T,P);

        // start value for the nonlinear solver
        double V0 = (stateCO2 == GAS || stateCO2 == SUPERCRITICAL) ? (R * T/P) : 59.22; // vCO2 = 59.22 mol/cm^3 at 25 °C and 80.0 bar

        // the nonlinear equation
        static VolumeCO2Equation equation;
        // the derivative of the nonlinear equation
        static VolumeCO2DerivativeEquation derivative;
        // setting the parameters of the equation
        equation.T = T; derivative.T = T;
        equation.P = P; derivative.P = P;

        double result=NewtonRaphson(equation, derivative, V0);
        assert(result>=0.0);
        return result;
}

// the nonlinear equation for the computation of the volume of CO2, Spycher et al. (2003)
struct VolumeH2SEquation
{
        double T, P;
        double operator()(double V)
        {
                static double Tr1; // 1/Tr = Tc/T, where Tc is critical temperature
                static double Tr2; // 1/Tr^2
                static double Tr3; // 1/Tr^3
                static double Vr1; // 1/Vr = R·Tc/V·Pc, where Pc is critical pressure
                static double Vr2; // 1/Vr^2
                static double Vr4; // 1/Vr^4
                static double Vr5; // 1/Vr^5

                Tr1 = TcH2S/T;
                Tr2 = Tr1 * Tr1;
                Tr3 = Tr2 * Tr1;
                Vr1 = (R * TcH2S)/(V * PcH2S);
                Vr2 = Vr1 * Vr1;
                Vr4 = Vr2 * Vr2;
                Vr5 = Vr4 * Vr1;

                return  1.0 +
                        Vr1 * (a1H2S + a2H2S * Tr2 + a3H2S * Tr3) +
                        Vr2 * (a4H2S + a5H2S * Tr2 + a6H2S * Tr3) +
                        Vr4 * (a7H2S + a8H2S * Tr2 + a9H2S * Tr3) +
                        Vr5 * (a10H2S + a11H2S * Tr2 + a12H2S * Tr3) +
                        a13H2S * Tr3 * Vr2 * (a14H2S + a15H2S * Vr2) * 
                        std::exp(-a15H2S * Vr2) - (P*V)/(R*T);
        }
};

// the derivative of the nonlinear equation for the computation of the volume of H2S, Duan et al. (2007)
struct VolumeH2SDerivativeEquation
{
        double T, P;
        double operator()(double V)
        {
                double Tr1 = TcH2S/T;
                double Tr2 = Tr1 * Tr1;
                double Tr3 = Tr2 * Tr1;
                double Vr1 = (R * TcH2S)/(V * PcH2S);
                double Vr2 = Vr1 * Vr1;
                double Vr3 = Vr2 * Vr1;
                double Vr4 = Vr2 * Vr2;
                double Vr5 = Vr4 * Vr1;
                double Vr6 = Vr5 * Vr1;

                return - (Vr2 * (a1H2S + a2H2S * Tr2 + a3H2S * Tr3) +
                        2.0 * Vr3 * (a4H2S + a5H2S * Tr2 + a6H2S * Tr3) +
                        4.0 * Vr5 * (a7H2S + a8H2S * Tr2 + a9H2S * Tr3) +
                        5.0 * Vr6 * (a10H2S + a11H2S * Tr2 + a12H2S * Tr3) +
                        4.0 * a13H2S * a15H2S * Tr3 * Vr6 * (a14H2S + 2.0 * a15H2S * Vr2) * 
                        expf(-a15H2S * Vr2)) * PcH2S/(R*TcH2S) -  P/(R*T);
        }
};

// the tolerance used for the computation of the H2S volume equation (this should not be a very small value, otherwise the flash calculation is be very slow)
const double toleranceVolumeH2SEquation = 1.0E-2; // recommended value is 0.01

// it computes the volume of H2S, Duan and Sun (2007)
inline double ComputeVolumeH2S(double T, double P, PhysicalState stateH2S)
{
        // start value for the nonlinear solver
        double V0 = (stateH2S == GAS || stateH2S == SUPERCRITICAL) ? (R * T/P) :  44.60; // vH2S = 44.60 mol/cm^3 at 60.0°C and 300.0 bar
        // the nonlinear equation
        static VolumeH2SEquation equation;
        // the derivative of the nonlinear equation
        static VolumeH2SDerivativeEquation derivative;
        // setting the parameters of the equation
        equation.T = T; derivative.T = T;
        equation.P = P; derivative.P = P;
        // solving the nonlinear equation
        return NewtonRaphson(equation, derivative, V0, toleranceVolumeH2SEquation);
}       

//==================================================+
//--------------------------------------------------+
// Fugacity Coefficient Methods
//--------------------------------------------------+
//==================================================+
// it computes the fugacity coefficient of H2O, Spycher et al. (2003)
inline double ComputeFugacityH2O(double volumeCO2, double T, double P)
{
        // auxiliary variables for performance reasons
        double a1 = std::log(volumeCO2/(volumeCO2 - bCO2));
        double a2 = std::log(1 + bCO2/volumeCO2);
        double a3 = bCO2/(volumeCO2 + bCO2);
        double a4 = std::log((P * volumeCO2)/(R * T));
        double a5 = R * std::pow(T, 1.5);
        // compute the fugacity coefficient of H2O
        return std::exp(a1 + bH2O/(volumeCO2 - bCO2) - 2 * aH2OCO2/bCO2 * a2/a5 + (aCO2(T) * bH2O)/(bCO2 * bCO2) * (a2 - a3)/a5 - a4);
}

// it computes the fugacity coefficient of CO2, Spycher et al. (2003)
inline double ComputeFugacityCO2(double volumeCO2, double T, double P)
{
        // auxiliary variables for performance reasons
        double a1 = std::log(volumeCO2/(volumeCO2 - bCO2));
        double a2 = std::log(1 + bCO2/volumeCO2);
        double a3 = bCO2/(volumeCO2 + bCO2);
        double a4 = std::log((P * volumeCO2)/(R * T));
        double a5 = R * std::pow(T, 1.5);
        // compute the fugacity coefficient of CO2
        return std::exp(a1 + bCO2/(volumeCO2 - bCO2) - 2 * aCO2(T)/bCO2 * a2/a5 + aCO2(T)/bCO2 * (a2 - a3)/a5 - a4);
}

// it computes the fugacity coefficient of H2S, Duan and Sun (2007)
inline double ComputeFugacityH2S(double volumeH2S, double T, double P)
{
        // static auxiliary variables
        static double Z;   // PV/RT
        static double Tr1; // 1/Tr = Tc/T, where Tc is critical temperature
        static double Tr2; // 1/Tr^2 
        static double Tr3; // 1/Tr^3
        static double Vr1; // 1/Vr = R·Tc/V·Pc, where Pc is critical pressure
        static double Vr2; // 1/Vr^2
        static double Vr4; // 1/Vr^4
        static double Vr5; // 1/Vr^5
        
        Z   = (P * volumeH2S)/(R * T);
        Tr1 = TcH2S/T;
        Tr2 = Tr1 * Tr1;
        Tr3 = Tr2 * Tr1;
        Vr1 = (R * TcH2S)/(volumeH2S * PcH2S);
        Vr2 = Vr1 * Vr1;
        Vr4 = Vr2 * Vr2;
        Vr5 = Vr4 * Vr1;        
        
        return std::exp(Z - 1.0 - std::log(Z) +
                Vr1 * (a1H2S + a2H2S * Tr2 + a3H2S * Tr3) + 
                Vr2 * (a4H2S + a5H2S * Tr2 + a6H2S * Tr3) * 0.50 +
                Vr4 * (a7H2S + a8H2S * Tr2 + a9H2S * Tr3) * 0.25 +
                Vr5 * (a10H2S + a11H2S * Tr2 + a12H2S * Tr3) * 0.20 + 
                a13H2S/a15H2S * Tr3 * 0.50 * 
                (a14H2S + 1.0 - (a14H2S + 1.0 + a15H2S * Vr2) * 
                std::exp(-a15H2S * Vr2)));
}
} // namespace Flash

#endif /* FLASHESSENTIAL_H */

---