VLERoutines.h

Go to the documentation of this file.

#ifndef VLEROUTINES_H
#define VLEROUTINES_H
 
#include "HelmholtzEOSMixtureBackend.h"
 
namespace CoolProp {
 
namespace SaturationSolvers {
struct saturation_T_pure_Akasaka_options
{
    bool use_guesses;  
    CoolPropDbl omega, rhoL, rhoV, pL, pV;
    saturation_T_pure_Akasaka_options(bool use_guesses = false)
      : use_guesses(use_guesses), omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE) {}
};
struct saturation_T_pure_options
{
    CoolPropDbl omega, rhoL, rhoV, pL, pV, p, T;
    saturation_T_pure_options() : omega(_HUGE), rhoL(_HUGE), rhoV(_HUGE), pL(_HUGE), pV(_HUGE), p(_HUGE), T(_HUGE) {}
};
 
struct saturation_D_pure_options
{
    enum imposed_rho_options
    {
        IMPOSED_RHOL,
        IMPOSED_RHOV
    };
    bool use_guesses,  
      use_logdelta;    
    CoolPropDbl omega, rhoL, rhoV, pL, pV;
    int imposed_rho;
    int max_iterations;
    saturation_D_pure_options()
      : use_guesses(false),
        use_logdelta(true),
        omega(1.0),
        rhoL(_HUGE),
        rhoV(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        imposed_rho(0),
        max_iterations(200) {}  // Defaults
};
 
enum sstype_enum
{
    imposed_T,
    imposed_p
};
struct mixture_VLE_IO
{
    sstype_enum sstype;
    int Nstep_max;
    CoolPropDbl rhomolar_liq, rhomolar_vap, p, T, beta;
    std::vector<CoolPropDbl> x, y, K;
};
 
static CoolPropDbl Wilson_lnK_factor(const HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl p, std::size_t i) {
    const double pci = HEOS.get_fluid_constant(i, iP_critical);
    const double Tci = HEOS.get_fluid_constant(i, iT_critical);
    const double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
    return log(pci / p) + 5.373 * (1 + omegai) * (1 - Tci / T);
};
 
void saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options);
void saturation_T_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);
void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);
void saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_Akasaka_options& options);
 
struct saturation_PHSU_pure_options
{
    enum specified_variable_options
    {
        IMPOSED_HL,
        IMPOSED_HV,
        IMPOSED_PL,
        IMPOSED_PV,
        IMPOSED_SL,
        IMPOSED_SV,
        IMPOSED_UL,
        IMPOSED_UV,
        IMPOSED_INVALID_INPUT
    };
    bool use_guesses,  
      use_logdelta;    
    specified_variable_options specified_variable;
    CoolPropDbl omega, rhoL, rhoV, pL, pV, T, p;
    saturation_PHSU_pure_options()
      : use_guesses(true),
        use_logdelta(true),
        specified_variable(IMPOSED_INVALID_INPUT),
        omega(1.0),
        rhoL(_HUGE),
        rhoV(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        T(_HUGE),
        p(_HUGE) {}
};
 
#if !defined(NO_FMTLIB) && FMT_VERSION >= 90000
inline int format_as(saturation_PHSU_pure_options::specified_variable_options option) {
    return fmt::underlying(option);
}
#endif
 
void saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl specified_value, saturation_PHSU_pure_options& options);
 
/* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution
     *
     * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param p Imposed pressure in kPa
     * @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided)
     */
void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, saturation_PHSU_pure_options& options);
 
/* \brief This is a backup saturation_T solver for the case where the Newton solver cannot approach closely enough the solution
     *
     * This is especially a problem at low pressures where catastrophic truncation error occurs, especially in the saturated vapor side
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param T Imposed temperature in K
     * @param options Options to be passed to the function (at least p, rhoL and rhoV must be provided)
     */
void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, saturation_T_pure_options& options);
 
/* \brief A robust but slow solver in the very-near-critical region
     *
     * This solver operates in the following fashion:
     * 1. Using a bounded interval for rho'':[rhoc, rhoc-??], guess a value for rho''
     * 2. For guessed value of rho'' and given value of T, calculate p
     * 3. Using a Brent solver on the other co-existing phase (rho'), calculate the (bounded) value of rho' that yields the same pressure
     * 4. Use another outer Brent solver on rho'' to enforce the same Gibbs function between liquid and vapor
     * 5. Fin.
     *
     * @param HEOS The Helmholtz EOS backend instance to be used
     * @param ykey The CoolProp::parameters key to be imposed - one of iT or iP
     * @param y The value for the imposed variable
     */
void saturation_critical(HelmholtzEOSMixtureBackend& HEOS, CoolProp::parameters ykey, CoolPropDbl y);
 
void successive_substitution(HelmholtzEOSMixtureBackend& HEOS, const CoolPropDbl beta, CoolPropDbl T, CoolPropDbl p,
                             const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K, mixture_VLE_IO& options);
void x_and_y_from_K(CoolPropDbl beta, const std::vector<CoolPropDbl>& K, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& x,
                    std::vector<CoolPropDbl>& y);
 
void successive_substitution_guessrho(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y,
                                      CoolPropDbl& rhomolar_liq, CoolPropDbl& rhomolar_vap, const std::vector<CoolPropDbl>& z, int num_steps,
                                      double tol = 1e-6);
 
bool guess_split_from_wilson(HelmholtzEOSMixtureBackend& HEOS, std::vector<CoolPropDbl>& x, std::vector<CoolPropDbl>& y, CoolPropDbl& rhomolar_liq,
                             CoolPropDbl& rhomolar_vap, const std::vector<CoolPropDbl>& z, CoolPropDbl T, CoolPropDbl p, int num_steps,
                             bool require_bracket = true);
 
class WilsonK_resid : public FuncWrapper1D
{
   public:
    sstype_enum input_type;
    double T, p, beta;
    const std::vector<CoolPropDbl>& z;
    std::vector<CoolPropDbl>& K;
    const HelmholtzEOSMixtureBackend& HEOS;
 
    WilsonK_resid(const HelmholtzEOSMixtureBackend& HEOS, double beta, double imposed_value, sstype_enum input_type,
                  const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& K)
      : input_type(input_type),
        T(imposed_value),
        p(imposed_value),
        beta(beta),
        z(z),
        K(K),
        HEOS(HEOS) {}  // if input_type == imposed_T -> use T, else use p; init both
    double call(double input_value) override {
        double summer = 0;
        if (input_type == imposed_T) {
            p = input_value;  // Iterate on pressure
        } else {
            T = input_value;  // Iterate on temperature, pressure imposed
        }
        for (unsigned int i = 0; i < z.size(); i++) {
            K[i] = exp(Wilson_lnK_factor(HEOS, T, p, i));
            summer += z[i] * (K[i] - 1) / (1 - beta + beta * K[i]);
        }
        return summer;
    }
};
inline double saturation_preconditioner(HelmholtzEOSMixtureBackend& HEOS, double input_value, sstype_enum input_type,
                                        const std::vector<CoolPropDbl>& z) {
    double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;
 
    if (HEOS.get_components().empty()) {
        return -1;
    }
 
    for (unsigned int i = 0; i < z.size(); i++) {
        Tcrit += HEOS.get_fluid_constant(i, iT_critical) * z[i];
        pcrit += HEOS.get_fluid_constant(i, iP_critical) * z[i];
        Ttriple += HEOS.get_fluid_constant(i, iT_triple) * z[i];
        ptriple += HEOS.get_fluid_constant(i, iP_triple) * z[i];
    }
    // Return an invalid number if either triple point temperature or pressure are not available
    if (!ValidNumber(Ttriple) || !ValidNumber(ptriple)) {
        return _HUGE;
    }
 
    if (input_type == imposed_T) {
        return exp(log(pcrit / ptriple) / (Tcrit - Ttriple) * (input_value - Ttriple) + log(ptriple));
    } else if (input_type == imposed_p) {
        return 1 / (1 / Tcrit - (1 / Ttriple - 1 / Tcrit) / log(pcrit / ptriple) * log(input_value / pcrit));
    } else {
        throw ValueError();
    }
}
inline double saturation_Wilson(HelmholtzEOSMixtureBackend& HEOS, double beta, double input_value, sstype_enum input_type,
                                const std::vector<CoolPropDbl>& z, double guess) {
    double out = 0;
    std::string errstr;
 
    // If T is input and beta = 0 or beta = 1, explicit solution for p is possible
    if (input_type == imposed_T && (std::abs(beta) < 1e-12 || std::abs(beta - 1) < 1e-12)) {
        const std::vector<double> z = HEOS.get_mole_fractions_ref();
        bool beta0 = std::abs(beta) < 1e-12;  // True is beta is approx. zero
        for (int i = 0; i < static_cast<int>(z.size()); ++i) {
            double pci = HEOS.get_fluid_constant(i, iP_critical);
            double Tci = HEOS.get_fluid_constant(i, iT_critical);
            double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
            if (beta0) {
                out += z[i] * pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
            } else {
                out += z[i] / (pci * exp(5.373 * (1 + omegai) * (1 - Tci / input_value)));
            }
        }
        if (!beta0) {       // beta = 1
            out = 1 / out;  // summation is for 1/p, take reciprocal to get p
        }
        std::vector<CoolPropDbl>& K = HEOS.get_K();
        for (int i = 0; i < static_cast<int>(z.size()); ++i) {
            double pci = HEOS.get_fluid_constant(i, iP_critical);
            double Tci = HEOS.get_fluid_constant(i, iT_critical);
            double omegai = HEOS.get_fluid_constant(i, iacentric_factor);
            K[i] = pci / out * exp(5.373 * (1 + omegai) * (1 - Tci / input_value));
        }
    } else {
        // Find first guess for output variable using Wilson K-factors.
        // Always try Brent (bounded) first — Secant (unbounded) can diverge
        // for multi-component mixtures at certain Q values.
        // Fix from jakobreichert PR CoolProp/CoolProp#2720.
        WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K());
        try {
            if (input_type == imposed_T) {
                out = Brent(Resid, 1.0, 1e9, 1e-10, 1e-10, 100);
            } else {
                out = Brent(Resid, 50, 10000, 1e-10, 1e-10, 100);
            }
        } catch (const std::exception&) {
            // Brent failed to bracket; fall back to Secant from the provided guess
            if (!ValidNumber(guess) || guess < 0) {
                throw;
            }
            out = Secant(Resid, guess, 0.001, 1e-10, 100);
        }
        if (!ValidNumber(out)) {
            throw ValueError("saturation_p_Wilson failed to get good output value");
        }
    }
    return out;
}
struct SuccessiveSubstitutionStep
{
    CoolPropDbl T, p;
};
 
struct newton_raphson_twophase_options
{
    enum imposed_variable_options
    {
        NO_VARIABLE_IMPOSED = 0,
        P_IMPOSED,
        T_IMPOSED
    };
    int Nstep_max;
    std::size_t Nsteps;
    CoolPropDbl beta, omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
    imposed_variable_options imposed_variable;
    std::vector<CoolPropDbl> x, y, z;
    newton_raphson_twophase_options()
      : Nstep_max(30),
        Nsteps(0),
        beta(-1),
        omega(1),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        p(_HUGE),
        T(_HUGE),
        hmolar_liq(_HUGE),
        hmolar_vap(_HUGE),
        smolar_liq(_HUGE),
        smolar_vap(_HUGE),
        imposed_variable(NO_VARIABLE_IMPOSED) {}  // Defaults
};
 
class newton_raphson_twophase
{
   public:
    HelmholtzEOSMixtureBackend* HEOS;
    newton_raphson_twophase_options::imposed_variable_options imposed_variable;
    CoolPropDbl error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta;
    std::size_t N;
    bool logging;
    int Nsteps;
    Eigen::MatrixXd J;
    // r and err_rel hold 2N-1 residuals/relative-errors (see class docstring), so they MUST
    // be dynamically sized.  They were Eigen::Vector2d (fixed size 2), which overflowed for
    // every mixture (N>=2 -> 2N-1>=3) and hard-crashed via out-of-bounds writes (GH #3192).
    Eigen::VectorXd r, err_rel;
    std::vector<CoolPropDbl> K, x, y, z;
    std::vector<SuccessiveSubstitutionStep> step_logger;
 
    newton_raphson_twophase()
      : HEOS(nullptr),
        imposed_variable(newton_raphson_twophase_options::NO_VARIABLE_IMPOSED),
        error_rms(_HUGE),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        T(_HUGE),
        p(_HUGE),
        min_rel_change(_HUGE),
        beta(_HUGE),
        N(0),
        logging(false),
        Nsteps(0) {};
 
    void resize(unsigned int N);
 
    // Reset the state of all the internal variables
    void pre_call() {
        K.clear();
        x.clear();
        y.clear();
        step_logger.clear();
        error_rms = 1e99;
        Nsteps = 0;
        rhomolar_liq = _HUGE;
        rhomolar_vap = _HUGE;
        T = _HUGE;
        p = _HUGE;
    };
 
    void call(HelmholtzEOSMixtureBackend& HEOS, newton_raphson_twophase_options& IO);
 
    /* \brief Build the arrays for the Newton-Raphson solve
         *
         */
    void build_arrays();
};
 
struct newton_raphson_saturation_options
{
    enum imposed_variable_options
    {
        NO_VARIABLE_IMPOSED = 0,
        P_IMPOSED,
        RHOV_IMPOSED,
        T_IMPOSED
    };
    int Nstep_max;
    bool bubble_point;
    std::size_t Nsteps;
    CoolPropDbl omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, hmolar_liq, hmolar_vap, smolar_liq, smolar_vap;
    imposed_variable_options imposed_variable;
    std::vector<CoolPropDbl> x, y;
    newton_raphson_saturation_options()
      : Nstep_max(30),
        bubble_point(false),
        Nsteps(0),
        omega(_HUGE),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        p(_HUGE),
        T(_HUGE),
        hmolar_liq(_HUGE),
        hmolar_vap(_HUGE),
        smolar_liq(_HUGE),
        smolar_vap(_HUGE),
        imposed_variable(NO_VARIABLE_IMPOSED) {}  // Defaults
};
 
class newton_raphson_saturation
{
   public:
    newton_raphson_saturation_options::imposed_variable_options imposed_variable = newton_raphson_saturation_options::NO_VARIABLE_IMPOSED;
    CoolPropDbl error_rms = _HUGE, rhomolar_liq = _HUGE, rhomolar_vap = _HUGE, T = _HUGE, p = _HUGE, min_rel_change = _HUGE;
    std::size_t N = 0;
    bool logging = false;
    bool bubble_point = false;
    int Nsteps = 0;
    Eigen::MatrixXd J;
    HelmholtzEOSMixtureBackend* HEOS = nullptr;
    CoolPropDbl dTsat_dPsat = _HUGE, dPsat_dTsat = _HUGE;
    std::vector<CoolPropDbl> K, x, y;
    Eigen::VectorXd r, err_rel;
    std::vector<SuccessiveSubstitutionStep> step_logger;
 
    newton_raphson_saturation() = default;
 
    void resize(std::size_t N);
 
    // Reset the state of all the internal variables
    void pre_call() {
        step_logger.clear();
        error_rms = 1e99;
        Nsteps = 0;
        rhomolar_liq = _HUGE;
        rhomolar_vap = _HUGE;
        T = _HUGE;
        p = _HUGE;
    };
 
    void call(HelmholtzEOSMixtureBackend& HEOS, const std::vector<CoolPropDbl>& z, std::vector<CoolPropDbl>& z_incipient,
              newton_raphson_saturation_options& IO);
 
    void build_arrays();
 
    void check_Jacobian();
};
 
struct PTflash_twophase_options
{
    int Nstep_max;       
    std::size_t Nsteps;  
    CoolPropDbl omega, rhomolar_liq, rhomolar_vap, pL, pV, p, T, beta;
    std::vector<CoolPropDbl> x,  
      y,                         
      z;                         
    bool nonconvergence;
    PTflash_twophase_options()
      : Nstep_max(30),
        Nsteps(0),
        omega(_HUGE),
        rhomolar_liq(_HUGE),
        rhomolar_vap(_HUGE),
        pL(_HUGE),
        pV(_HUGE),
        p(_HUGE),
        T(_HUGE),
        beta(0.5),
        nonconvergence(false) {}
};
 
class PTflash_twophase
{
   public:
    double error_rms = _HUGE;
    bool logging = false;
    int Nsteps = 0;
    Eigen::MatrixXd J;
    Eigen::VectorXd r, err_rel;
    HelmholtzEOSMixtureBackend& HEOS;
    PTflash_twophase_options& IO;
    std::vector<SuccessiveSubstitutionStep> step_logger;
 
    PTflash_twophase(HelmholtzEOSMixtureBackend& HEOS, PTflash_twophase_options& IO) : HEOS(HEOS), IO(IO) {};
 
    void solve();
    void solve_legacy();
    void solve_michelsen();
 
    void build_arrays();
};
};  // namespace SaturationSolvers
 
namespace StabilityRoutines {
 
class StabilityEvaluationClass
{
   protected:
    HelmholtzEOSMixtureBackend& HEOS;
    std::vector<double> lnK, K, K0, x, y, xL, xH;
    const std::vector<double>& z;
    double rhomolar_liq, rhomolar_vap, beta, tpd_liq, tpd_vap, DELTAG_nRT;
    double m_T,  
      m_p;       
 
   private:
    bool _stable;
    bool _uncertain;  
    bool debug;
    bool use_michelsen;
 
   public:
    StabilityEvaluationClass(HelmholtzEOSMixtureBackend& HEOS)
      : HEOS(HEOS),
        z(HEOS.get_mole_fractions_doubleref()),
        rhomolar_liq(-1),
        rhomolar_vap(-1),
        beta(-1),
        tpd_liq(10000),
        tpd_vap(100000),
        DELTAG_nRT(10000),
        m_T(-1),
        m_p(-1),
        _stable(false),
        _uncertain(false),
        debug(false),
        use_michelsen(get_config_int(MIXTURE_STABILITY_ALGORITHM) != 0) {};
    void set_TP(double T, double p) {
        m_T = T;
        m_p = p;
    };
    void rho_TP_w_guesses();
    void rho_TP_global();
    void rho_TP_SRK_translated();
 
    void trial_compositions();
    void successive_substitution(int num_steps);
 
    void check_stability();
    void check_stability_legacy();
    void check_stability_michelsen();
 
    bool minimize_tpd(std::vector<CoolPropDbl>& Y, const std::vector<CoolPropDbl>& ln_f_z, CoolPropDbl the_T, CoolPropDbl the_p, bool& is_unstable);
 
    void set_use_michelsen(bool value) {
        use_michelsen = value;
    }
 
    bool is_stable() {
        if (!use_michelsen) {
            // Legacy path needs trial compositions and a few SS steps
            // before check_stability_legacy() can evaluate TPD
            trial_compositions();
            successive_substitution(3);
        }
        check_stability();
        return _stable;
    }
    bool is_uncertain() const {
        return _uncertain;
    }
    void get_liq(std::vector<double>& x, double& rhomolar) {
        x = this->x;
        rhomolar = rhomolar_liq;
    }
    void get_vap(std::vector<double>& y, double& rhomolar) {
        y = this->y;
        rhomolar = rhomolar_vap;
    }
};
 
}; /* namespace StabilityRoutines*/
 
}; /* namespace CoolProp*/
 
#endif