FlashRoutines.h

Go to the documentation of this file.

 
// ***************************************************************
// *******************  FLASH ROUTINES  **************************
// ***************************************************************
 
#ifndef FLASHROUTINES_H
#define FLASHROUTINES_H
 
#include "HelmholtzEOSMixtureBackend.h"
#include "CoolProp/numerics/Solvers.h"
#include <optional>
 
namespace CoolProp {
 
class FlashRoutines
{
   public:
    template <class T>
    T static g_RachfordRice(const std::vector<T>& z, const std::vector<T>& lnK, T beta) {
        // g function from Rachford-Rice
        T summer = 0;
        for (std::size_t i = 0; i < z.size(); i++) {
            T Ki = exp(lnK[i]);
            summer += z[i] * (Ki - 1) / (1 - beta + beta * Ki);
        }
        return summer;
    }
    template <class T>
    T static dgdbeta_RachfordRice(const std::vector<T>& z, const std::vector<T>& lnK, T beta) {
        // derivative of g function from Rachford-Rice with respect to beta
        T summer = 0;
        for (std::size_t i = 0; i < z.size(); i++) {
            T Ki = exp(lnK[i]);
            summer += -z[i] * pow((Ki - 1) / (1 - beta + beta * Ki), 2);
        }
        return summer;
    }
 
    static void PQ_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void PQ_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static void QT_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static void PT_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static void QT_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void QS_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void DQ_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void DQ_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static void HQ_flash(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl Tguess = -1);
 
    static void HQ_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static void QS_flash_with_guesses(HelmholtzEOSMixtureBackend& HEOS, const GuessesStructure& guess);
 
    static double resolve_T_via_superancillary(HelmholtzEOSMixtureBackend& HEOS, parameters key, double target_value, std::optional<double> guess_T,
                                               const char* fn_name);
 
    static bool sat_superanc_path_applies(HelmholtzEOSMixtureBackend& HEOS);
 
    static std::pair<double, double> alpha0_offset_total(HelmholtzEOSMixtureBackend& HEOS);
 
    static void PT_Q_flash_mixtures(HelmholtzEOSMixtureBackend& HEOS, parameters other, CoolPropDbl value);
 
    static void PT_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void PT_flash_mixtures(HelmholtzEOSMixtureBackend& HEOS);
 
    static double T_DP_PengRobinson(HelmholtzEOSMixtureBackend& HEOS, double rhomolar, double p);
 
    static void DP_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void solver_for_rho_given_T_oneof_HSU(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl value, parameters other);
 
    static void DHSU_T_flash(HelmholtzEOSMixtureBackend& HEOS, parameters other);
 
    static void HSU_P_flash(HelmholtzEOSMixtureBackend& HEOS, parameters other);
 
    static void HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend& HEOS, parameters other, CoolPropDbl T0, CoolPropDbl rhomolar0);
 
    static void HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend& HEOS, parameters other, CoolPropDbl value, CoolPropDbl Tmin,
                                              CoolPropDbl Tmax, phases phase);
 
    static void HSU_D_flash_twophase(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar_spec, parameters other, CoolPropDbl value);
 
    static void HSU_D_flash(HelmholtzEOSMixtureBackend& HEOS, parameters other);
 
    static void HS_flash(HelmholtzEOSMixtureBackend& HEOS);
 
    static void HS_flash_generate_TP_singlephase_guess(HelmholtzEOSMixtureBackend& HEOS, double& T, double& p);
 
    struct HS_flash_singlephaseOptions
    {
        double omega;
        HS_flash_singlephaseOptions() : omega(1.0) {}
    };
    static void HS_flash_singlephase(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl hmolar_spec, CoolPropDbl smolar_spec,
                                     HS_flash_singlephaseOptions& options);
 
    struct HS_flash_twophaseOptions
    {
        double omega;
        HS_flash_twophaseOptions() : omega(1.0) {}
    };
    static void HS_flash_twophase(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl hmolar_spec, CoolPropDbl smolar_spec,
                                  HS_flash_twophaseOptions& options);
 
    static bool hs_corrector_probe(HelmholtzEOSMixtureBackend& H, double T0, double rho0, double h_t, double s_t, double& T_out, double& rho_out,
                                   double Tlo_override);
};
 
class solver_TP_resid : public FuncWrapper1DWithDeriv
{
   public:
    HelmholtzEOSMixtureBackend* HEOS;
    CoolPropDbl T, p, rhor, tau, R_u, delta;
 
    solver_TP_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl T, CoolPropDbl p)
      : HEOS(&HEOS),
        T(T),
        p(p),
        rhor(HEOS.get_reducing_state().rhomolar),
        tau(HEOS.get_reducing_state().T / T),
        R_u(HEOS.gas_constant()),
        delta(-_HUGE) {}
    double call(double rhomolar) override {
        delta = rhomolar / rhor;  // needed for derivative
        HEOS->update_DmolarT_direct(rhomolar, T);
        CoolPropDbl peos = HEOS->p();
        return (peos - p) / p;
    };
    double deriv(double rhomolar) override {
        // dp/drho|T / pspecified
        return R_u * T * (1 + 2 * delta * HEOS->dalphar_dDelta() + pow(delta, 2) * HEOS->d2alphar_dDelta2()) / p;
    };
};
 
class PY_singlephase_flash_resid : public FuncWrapper1D
{
   public:
    HelmholtzEOSMixtureBackend* HEOS;
    CoolPropDbl p;
    parameters other;
    CoolPropDbl value;
    PY_singlephase_flash_resid(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl p, parameters other, CoolPropDbl value)
      : HEOS(&HEOS), p(p), other(other), value(value) {
        // Specify the state to avoid saturation calls, but only if phase is subcritical
        if (HEOS.phase() == iphase_liquid || HEOS.phase() == iphase_gas) {
            HEOS.specify_phase(HEOS.phase());
        }
    };
    double call(double T) override {
 
        // Run the solver with T,P as inputs;
        HEOS->update(PT_INPUTS, p, T);
 
        CoolPropDbl rhomolar = HEOS->rhomolar();
        HEOS->update(DmolarT_INPUTS, rhomolar, T);
        // Get the value of the desired variable
        CoolPropDbl eos = HEOS->keyed_output(other);
 
        // Difference between the two is to be driven to zero
        return eos - value;
    };
};
 
} /* namespace CoolProp */
#endif /* FLASHROUTINES_H */