GeneralizedCubic.h

Go to the documentation of this file.

 
#ifndef CUBIC_H
#define CUBIC_H
 
#include <vector>
#include <cmath>
#include "crossplatform_shared_ptr.h"
#include "Exceptions.h"
 
class AbstractCubicAlphaFunction
{
   protected:
    double a0;              
    double Tr_over_Tci,     
      sqrt_Tr_Tci;          
    std::vector<double> c;  
   public:
    virtual ~AbstractCubicAlphaFunction(){};
    virtual double term(double tau, std::size_t itau) = 0;
    void set_Tr_over_Tci(double Tr_over_Tci) {
        this->Tr_over_Tci = Tr_over_Tci;
        this->sqrt_Tr_Tci = sqrt(Tr_over_Tci);
    };
    AbstractCubicAlphaFunction(double a0, double Tr_over_Tci) : a0(a0), Tr_over_Tci(Tr_over_Tci), sqrt_Tr_Tci(sqrt(Tr_over_Tci)){};
};
 
class BasicMathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
{
    double m;  
   public:
    BasicMathiasCopemanAlphaFunction(double a0, double m_ii, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
        this->m = m_ii;
    };
    double term(double tau, std::size_t itau);
};
 
class TwuAlphaFunction : public AbstractCubicAlphaFunction
{
   public:
    TwuAlphaFunction(double a0, double L, double M, double N, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
        c.resize(3);
        c[0] = L;
        c[1] = M;
        c[2] = N;
    };
    double term(double tau, std::size_t itau);
};
 
class MathiasCopemanAlphaFunction : public AbstractCubicAlphaFunction
{
   public:
    MathiasCopemanAlphaFunction(double a0, double c1, double c2, double c3, double Tr_over_Tci) : AbstractCubicAlphaFunction(a0, Tr_over_Tci) {
        c.resize(3);
        c[0] = c1;
        c[1] = c2;
        c[2] = c3;
    };
    double term(double tau, std::size_t itau);
};
 
class AbstractCubic
{
 
   protected:
    double rho_r,                                               
      T_r;                                                      
    std::vector<double> Tc,                                     
      pc,                                                       
      acentric;                                                 
    double R_u;                                                 
    double Delta_1,                                             
      Delta_2;                                                  
    int N;                                                      
    std::vector<std::vector<double>> k;                         
    double cm;                                                  
    std::vector<shared_ptr<AbstractCubicAlphaFunction>> alpha;  
   public:
    AbstractCubic(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, double Delta_1, double Delta_2,
                  std::vector<double> C1 = std::vector<double>(), std::vector<double> C2 = std::vector<double>(),
                  std::vector<double> C3 = std::vector<double>());
    virtual ~AbstractCubic(){};
    void set_alpha(const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3);
    void set_alpha_function(std::size_t i, shared_ptr<AbstractCubicAlphaFunction>& acaf) {
        alpha[i] = acaf;
    };
    shared_ptr<AbstractCubicAlphaFunction> get_alpha_function(std::size_t i) {
        return alpha[i];
    };
    const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& get_all_alpha_functions() {
        return this->alpha;
    };
    void set_all_alpha_functions(const std::vector<shared_ptr<AbstractCubicAlphaFunction>>& alpha) {
        this->alpha = alpha;
    };
 
    const std::vector<std::vector<double>>& get_kmat() {
        return k;
    };
    void set_kmat(const std::vector<std::vector<double>>& k) {
        this->k = k;
    };
    void set_kij(std::size_t i, std::size_t j, double val) {
        k[i][j] = val;
        k[j][i] = val;
    }
    double get_kij(std::size_t i, std::size_t j) {
        return k[i][j];
    }
    std::vector<double>& get_Tc() {
        return Tc;
    }
    std::vector<double>& get_pc() {
        return pc;
    }
    std::vector<double>& get_acentric() {
        return acentric;
    }
    double get_Delta_1() {
        return Delta_1;
    }
    double get_Delta_2() {
        return Delta_2;
    }
    double get_R_u() {
        return R_u;
    }
    void set_Tr(double Tr) {
        T_r = Tr;
        for (std::size_t i = 0; i < alpha.size(); ++i) {
            alpha[i]->set_Tr_over_Tci(T_r / Tc[i]);
        }
    }
    void set_rhor(double rhor) {
        rho_r = rhor;
    }
    double get_Tr() {
        return T_r;
    }
    double get_rhor() {
        return rho_r;
    }
 
    void set_C_MC(std::size_t i, double c1, double c2, double c3) {
        alpha[i].reset(new MathiasCopemanAlphaFunction(a0_ii(i), c1, c2, c3, T_r / Tc[i]));
    }
    void set_C_Twu(std::size_t i, double L, double M, double N) {
        alpha[i].reset(new TwuAlphaFunction(a0_ii(i), L, M, N, T_r / Tc[i]));
    }
    virtual double a0_ii(std::size_t i) = 0;
    virtual double b0_ii(std::size_t i) = 0;
    virtual double m_ii(std::size_t i) = 0;
 
    virtual double alphar(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
    virtual double d_alphar_dxi(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
                                bool xN_independent);
    virtual double d2_alphar_dxidxj(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
                                    std::size_t j, bool xN_independent);
    virtual double d3_alphar_dxidxjdxk(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
                                       std::size_t j, std::size_t k, bool xN_independent);
 
    virtual double am_term(double tau, const std::vector<double>& x, std::size_t itau);
    virtual double d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
    virtual double d2_am_term_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, bool xN_independent);
    virtual double d3_am_term_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
                                        bool xN_independent);
 
    virtual double bm_term(const std::vector<double>& x);
    virtual double d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent);
    virtual double d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent);
    virtual double d3_bm_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent);
    virtual double cm_term();
    void set_cm(double val) {
        cm = val;
    }
    double get_cm() {
        return cm;
    }
 
    virtual void set_Q_k(const size_t sgi, const double value) {
        throw CoolProp::ValueError("set_Q_k not defined for AbstractCubic");
    };
    virtual double get_Q_k(const size_t sgi) const {
        throw CoolProp::ValueError("get_Q_k not defined for AbstractCubic");
    };
 
    double aij_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
    double u_term(double tau, std::size_t i, std::size_t j, std::size_t itau);
    double aii_term(double tau, std::size_t i, std::size_t itau);
 
    double psi_minus(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta);
    double d_psi_minus_dxi(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, bool xN_independent);
    double d2_psi_minus_dxidxj(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j,
                               bool xN_independent);
    double d3_psi_minus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i, std::size_t j,
                                  std::size_t k, bool xN_independent);
 
    double PI_12(double delta, const std::vector<double>& x, std::size_t idelta);
    double d_PI_12_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
    double d2_PI_12_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent);
    double d3_PI_12_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
                              bool xN_independent);
 
    double tau_times_a(double tau, const std::vector<double>& x, std::size_t itau);
    double d_tau_times_a_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent);
    double d2_tau_times_a_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, bool xN_independent);
    double d3_tau_times_a_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
                                    bool xN_independent);
 
    double psi_plus(double delta, const std::vector<double>& x, std::size_t idelta);
    double d_psi_plus_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent);
    double d2_psi_plus_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent);
    double d3_psi_plus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
                                 bool xN_independent);
 
    double c_term(const std::vector<double>& x) {
        return 1 / bm_term(x);
    };
    double d_c_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
        return -d_bm_term_dxi(x, i, xN_independent) / pow(bm_term(x), 2);
    };
    double d2_c_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
        double bm = bm_term(x);
        return (2 * d_bm_term_dxi(x, i, xN_independent) * d_bm_term_dxi(x, j, xN_independent) - bm * d2_bm_term_dxidxj(x, i, j, xN_independent))
               / pow(bm, 3);
    };
    double d3_c_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
        double bm = bm_term(x);
        return 1 / pow(bm, 4)
               * (2 * bm
                    * (d_bm_term_dxi(x, i, xN_independent) * d2_bm_term_dxidxj(x, j, k, xN_independent)
                       + d_bm_term_dxi(x, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
                       + d_bm_term_dxi(x, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent))
                  - pow(bm, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
                  - 6 * d_bm_term_dxi(x, i, xN_independent) * d_bm_term_dxi(x, j, xN_independent) * d_bm_term_dxi(x, k, xN_independent));
    };
 
    double A_term(double delta, const std::vector<double>& x) {
        double bm = bm_term(x);
        double cm = cm_term();
        return log((delta * rho_r * (Delta_1 * bm + cm) + 1) / (delta * rho_r * (Delta_2 * bm + cm) + 1));
    };
    double d_A_term_dxi(double delta, const std::vector<double>& x, std::size_t i, bool xN_independent) {
        std::size_t idelta = 0;
        return delta * rho_r * d_bm_term_dxi(x, i, xN_independent) * (Delta_1 - Delta_2) / PI_12(delta, x, idelta);
    };
    double d2_A_term_dxidxj(double delta, const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
        std::size_t idelta = 0;
        double PI12 = PI_12(delta, x, idelta);
        return delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 2)
               * (PI12 * d2_bm_term_dxidxj(x, i, j, xN_independent)
                  - d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_bm_term_dxi(x, i, xN_independent));
    };
    double d3_A_term_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
        std::size_t idelta = 0;
        double PI12 = PI_12(delta, x, idelta);
        // The leading factor
        double lead = delta * rho_r * (Delta_1 - Delta_2) / pow(PI12, 3);
        return lead
               * (-PI12
                    * (d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d2_bm_term_dxidxj(x, i, k, xN_independent)
                       + d_PI_12_dxi(delta, x, idelta, k, xN_independent) * d2_bm_term_dxidxj(x, i, j, xN_independent)
                       + d_bm_term_dxi(x, i, xN_independent) * d2_PI_12_dxidxj(delta, x, idelta, j, k, xN_independent))
                  + pow(PI12, 2) * d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent)
                  + 2 * d_PI_12_dxi(delta, x, idelta, j, xN_independent) * d_PI_12_dxi(delta, x, idelta, k, xN_independent)
                      * d_bm_term_dxi(x, i, xN_independent));
    };
    // Allows to modify the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
    virtual void set_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter, const double value) {
        throw CoolProp::NotImplementedError("set_interaction_parameter is not implemented for this backend");
    }
    // Allows to get the unifac interaction parameters aij, bij and cij. Only for use with VTPR backend.
    virtual double get_interaction_parameter(const std::size_t mgi1, const std::size_t mgi2, const std::string& parameter) {
        throw CoolProp::NotImplementedError("get_interaction_parameter is not implemented for this backend");
    }
};
 
class PengRobinson : public AbstractCubic
{
   public:
    PengRobinson(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u,
                 std::vector<double> C1 = std::vector<double>(), std::vector<double> C2 = std::vector<double>(),
                 std::vector<double> C3 = std::vector<double>())
      : AbstractCubic(Tc, pc, acentric, R_u, 1 + sqrt(2.0), 1 - sqrt(2.0), C1, C2, C3) {
        set_alpha(C1, C2, C3);
    };
 
    PengRobinson(double Tc, double pc, double acentric, double R_u)
      : AbstractCubic(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u, 1 + sqrt(2.0), 1 - sqrt(2.0)) {
        set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
    };
 
    double a0_ii(std::size_t i);
    double b0_ii(std::size_t i);
    double m_ii(std::size_t i);
};
 
class SRK : public AbstractCubic
{
   public:
    SRK(std::vector<double> Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, std::vector<double> C1 = std::vector<double>(),
        std::vector<double> C2 = std::vector<double>(), std::vector<double> C3 = std::vector<double>())
      : AbstractCubic(Tc, pc, acentric, R_u, 1, 0, C1, C2, C3) {
        set_alpha(C1, C2, C3);
    };
    SRK(double Tc, double pc, double acentric, double R_u)
      : AbstractCubic(std::vector<double>(1, Tc), std::vector<double>(1, pc), std::vector<double>(1, acentric), R_u, 1, 0) {
        set_alpha(std::vector<double>(), std::vector<double>(), std::vector<double>());
    };
 
    double a0_ii(std::size_t i);
    double b0_ii(std::size_t i);
    double m_ii(std::size_t i);
};
 
#endif