Mixtures.h

Go to the documentation of this file.

#ifndef MIXTURES_H
#define MIXTURES_H

#include "Helmholtz.h"
#include "FluidClass.h"

#include <vector>
#include <string>

typedef std::vector<std::vector<double> > STLMatrix;

class ReducingFunction
{
protected:
        unsigned int N;
public:
        ReducingFunction(){};
        virtual ~ReducingFunction(){};
        virtual double Tr(const std::vector<double> &x) = 0;
        virtual double dTrdxi__constxj(const std::vector<double> &x, int i) = 0;
        virtual double rhorbar(const std::vector<double> &x) = 0;
        virtual double drhorbardxi__constxj(const std::vector<double> &x, int i) = 0;

        virtual void set_coeffs_from_map(int i, int j, std::map<std::string,double >) = 0;

        virtual double d2rhorbardxi2__constxj(const std::vector<double> &x, int i) = 0;
        virtual double d2rhorbardxidxj(const std::vector<double> &x, int i, int j) = 0;
        virtual double d2Trdxi2__constxj(const std::vector<double> &x, int i) = 0;
        virtual double d2Trdxidxj(const std::vector<double> &x, int i, int j) = 0;

        double d_ndTrdni_dxj__constxi(const std::vector<double> &x, int i, int j);
        double d_ndrhorbardni_dxj__constxi(const std::vector<double> &x, int i, int j);

        double ndrhorbardni__constnj(const std::vector<double> &x, int i);
        double ndTrdni__constnj(const std::vector<double> &x, int i);
};

class GERG2008ReducingFunction : public ReducingFunction
{
protected:
        STLMatrix v_c;
        STLMatrix T_c; 
        STLMatrix beta_v; 
        STLMatrix gamma_v; 
        STLMatrix beta_T; 
        STLMatrix gamma_T; 
        std::vector<Fluid *> pFluids; 
public:
        GERG2008ReducingFunction(std::vector<Fluid *> pFluids, STLMatrix beta_v, STLMatrix gamma_v, STLMatrix beta_T, STLMatrix gamma_T)
        {
                this->pFluids = pFluids;
                this->beta_v = beta_v;
                this->gamma_v = gamma_v;
                this->beta_T = beta_T;
                this->gamma_T = gamma_T;
                this->N = pFluids.size();
        };
        GERG2008ReducingFunction(std::vector<Fluid *> pFluids)
        {
                this->pFluids = pFluids;
                this->N = pFluids.size();
                beta_v.resize(N,std::vector<double>(N,0));
                gamma_v.resize(N,std::vector<double>(N,0));
                beta_T.resize(N,std::vector<double>(N,0));
                gamma_T.resize(N,std::vector<double>(N,0));
                T_c.resize(N,std::vector<double>(N,0));
                v_c.resize(N,std::vector<double>(N,0));
                for (unsigned int i = 0; i < N; i++)
                {
                        for (unsigned int j = 0; j < N; j++)
                        {
                                T_c[i][j] = sqrt(pFluids[i]->reduce.T*pFluids[j]->reduce.T);
                                v_c[i][j] = 1.0/8.0*pow(pow(pFluids[i]->reduce.rhobar, -1.0/3.0)+pow(pFluids[j]->reduce.rhobar, -1.0/3.0),(int)3);
                        }
                }
        };
        ~GERG2008ReducingFunction(){};
        double Tr(const std::vector<double> &x);
        double dTrdxi__constxj(const std::vector<double> &x, int i);
        double rhorbar(const std::vector<double> &x);
        double drhorbardxi__constxj(const std::vector<double> &x, int i);
        double dvrbardxi__constxj(const std::vector<double> &x, int i);

        double d2vrbardxi2__constxj(const std::vector<double> &x, int i);
        double d2rhorbardxi2__constxj(const std::vector<double> &x, int i);
        double d2vrbardxidxj(const std::vector<double> &x, int i, int j);
        double d2rhorbardxidxj(const std::vector<double> &x, int i, int j);
        double d2Trdxi2__constxj(const std::vector<double> &x, int i);
        double d2Trdxidxj(const std::vector<double> &x, int i, int j);

        void set_coeffs_from_map(int i, int j, std::map<std::string,double >);

        double c_Y_ij(int i, int j, std::vector< std::vector< double> > * beta, std::vector< std::vector< double> > *gamma, std::vector< std::vector< double> > *Y_c);
        double c_Y_ji(int j, int i, std::vector< std::vector< double> > * beta, std::vector< std::vector< double> > *gamma, std::vector< std::vector< double> > *Y_c);
        double f_Y_ij(const std::vector<double> &x, int i, int j, std::vector< std::vector< double> > * beta);

        double dfYkidxi__constxk(const std::vector<double> &x, int k, int i,std::vector< std::vector< double> > * beta);
        double dfYikdxi__constxk(const std::vector<double> &x, int i, int k, std::vector< std::vector< double> > * beta);
        double d2fYkidxi2__constxk(const std::vector<double> &x, int k, int i, std::vector< std::vector< double> > * beta);
        double d2fYikdxi2__constxk(const std::vector<double> &x, int i, int k, std::vector< std::vector< double> > * beta);
        double d2fYijdxidxj(const std::vector<double> &x, int i, int k, std::vector< std::vector< double> > * beta);
};

class LemmonAirHFCReducingFunction : public GERG2008ReducingFunction
{
protected:
        std::vector<Fluid *> pFluids; 
public:
        LemmonAirHFCReducingFunction(std::vector<Fluid *> pFluids) : GERG2008ReducingFunction(pFluids)
        {
                this->pFluids = pFluids;
                this->N = pFluids.size();
        };

        void set_coeffs_from_map(int i, int j, std::map<std::string,double >);
        
};








class DepartureFunction
{
public:
        DepartureFunction(){};
        virtual ~DepartureFunction(){};
        
        virtual double phir(double tau, double delta) = 0;
        virtual double dphir_dDelta(double tau, double delta) = 0;
        virtual double d2phir_dDelta2(double tau, double delta) = 0;
        virtual double d2phir_dDelta_dTau(double tau, double delta) = 0;
        virtual double dphir_dTau(double tau, double delta) = 0;
        virtual double d2phir_dTau2(double tau, double delta) = 0;
        virtual void set_coeffs_from_map(std::map<std::string,std::vector<double> >) = 0;
};

struct DepartureFunctionCacheElement
{
        double tau, delta, cached_val;
};

struct DepartureFunctionCache
{
        DepartureFunctionCacheElement phir,dphir_dDelta,dphir_dTau,d2phir_dDelta2,d2phir_dTau2,d2phir_dDelta_dTau;
};

class GERG2008DepartureFunction : public DepartureFunction
{
protected:
        bool using_gaussian;
        DepartureFunctionCache cache;
        phir_power phi1;
        phir_GERG2008_gaussian phi2;
public:
        GERG2008DepartureFunction(){};
        ~GERG2008DepartureFunction(){};
        double phir(double tau, double delta);
        double dphir_dDelta(double tau, double delta);
        double d2phir_dDelta_dTau(double tau, double delta);
        double dphir_dTau(double tau, double delta);
        double d2phir_dDelta2(double tau, double delta);
        double d2phir_dTau2(double tau, double delta);
        void set_coeffs_from_map(std::map<std::string,std::vector<double> >);
};

class LemmonHFCDepartureFunction : public DepartureFunction
{
protected:
        bool using_gaussian;
        DepartureFunctionCache cache;
        phir_power phi1;
public:
        LemmonHFCDepartureFunction(){};
        ~LemmonHFCDepartureFunction(){};
        double phir(double tau, double delta);
        double dphir_dDelta(double tau, double delta);
        double d2phir_dDelta_dTau(double tau, double delta);
        double dphir_dTau(double tau, double delta);
        double d2phir_dDelta2(double tau, double delta);
        double d2phir_dTau2(double tau, double delta);
        void set_coeffs_from_map(std::map<std::string,std::vector<double> >);
};

class ExcessTerm
{
public:
        unsigned int N;
        std::vector<std::vector<DepartureFunction*> > DepartureFunctionMatrix;
        std::vector<std::vector<double> > F;
        ExcessTerm(int N);
        ~ExcessTerm();
        double phir(double tau, double delta, const std::vector<double> &x);
        double dphir_dDelta(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta2(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta_dTau(double tau, double delta, const std::vector<double> &x);
        double dphir_dTau(double tau, double delta, const std::vector<double> &x);
        double d2phir_dTau2(double tau, double delta, const std::vector<double> &x);
        double dphir_dxi(double tau, double delta, const std::vector<double> &x, unsigned int i);
        double d2phirdxidxj(double tau, double delta, const std::vector<double> &x, unsigned int i, unsigned int j);
        double d2phir_dxi_dTau(double tau, double delta, const std::vector<double> &x, unsigned int i);
        double d2phir_dxi_dDelta(double tau, double delta, const std::vector<double> &x, unsigned int i);
        void set_coeffs_from_map(int i, int j, std::map<std::string,std::vector<double> >);
};


class ResidualIdealMixture
{
protected:
        unsigned int N;
        std::vector<Fluid*> pFluids;
public:
        ResidualIdealMixture(std::vector<Fluid*> pFluids);
        double phir(double tau, double delta, const std::vector<double> &x);
        double dphir_dDelta(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta2(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta_dTau(double tau, double delta, const std::vector<double> &x);
        double d2phir_dTau2(double tau, double delta, const std::vector<double> &x);
        double dphir_dTau(double tau, double delta, const std::vector<double> &x);
};

class Mixture;  // Forward declaration since some classes that are members of Mixture take pointers to Mixture

struct SuccessiveSubstitutionStep
{
        double T, 
                   p; 
        std::vector<double> K, 
                                                lnK, 
                                                ln_phi_liq, 
                                                ln_phi_vap; 
};

class SuccessiveSubstitutionVLE
{
public:
        bool useNR; 
        bool logging; 
        int Nsteps; 
        int Nstep_max; 
        Mixture *Mix; 
        double rhobar_liq, 
                   rhobar_vap; 
        std::vector<double> K, 
                                            ln_phi_liq, 
                                                ln_phi_vap, 
                                                x, 
                                                y;
        std::vector<SuccessiveSubstitutionStep> step_logger;

        SuccessiveSubstitutionVLE(){useNR = false; logging = false; Nstep_max = 10;};

        double call(double beta, double T, double p, const std::vector<double> &z, std::vector<double> &K);
};

class NewtonRaphsonVLE
{
public:
        double error_rms, rhobar_liq, rhobar_vap, T, p;
        unsigned int N;
        bool logging;
        int Nsteps;
        int Nsteps_max;
        Mixture *Mix;
        STLMatrix J;
        std::vector<double> K, x, y, phi_ij_liq, phi_ij_vap, r, dXdS, neg_dFdS;
        std::vector<SuccessiveSubstitutionStep> step_logger;

        NewtonRaphsonVLE(){Nsteps_max = 20;};

        void resize(unsigned int N);
        
        // Reset the state of all the internal variables
        void pre_call()
        {
                K.clear();
                x.clear();
                y.clear();
                phi_ij_liq.clear();
                phi_ij_vap.clear();
                Nsteps = 0;
                step_logger.clear();
                error_rms = 1e99;
                rhobar_liq = _HUGE;
                rhobar_vap = _HUGE;
                T = _HUGE;
                p = _HUGE;
        };

        double call(double beta, double T, double p, double rhobar_liq, double rhobar_vap, const std::vector<double> &z, std::vector<double> &K, int spec_index, double spec_value);

        void build_arrays(double beta, double T, double p, const std::vector<double> &z, std::vector<double> & K, int spec_index, double spec_value);
};

struct PhaseEnvelopeLog
{
        std::vector< std::vector<double> > K, lnK, x, y;
        std::vector<double> T, p, lnT, lnp, rhobar_liq, rhobar_vap, lnrhobar_liq, lnrhobar_vap;
        std::vector<int> iS;

        void store_variables(const double T, 
                                 const double p, 
                                                 const double rhobar_liq, 
                                                 const double rhobar_vap, 
                                                 const std::vector<double> & K, 
                                                 const int iS, 
                                                 const std::vector<double> & x, 
                                                 const std::vector<double> & y, 
                                                 const unsigned int N)
        {
                this->p.push_back(p);
                this->T.push_back(T);
                this->lnT.push_back(log(T));
                this->lnp.push_back(log(p));
                this->rhobar_liq.push_back(rhobar_liq);
                this->rhobar_vap.push_back(rhobar_vap);
                this->lnrhobar_liq.push_back(log(rhobar_liq));
                this->lnrhobar_vap.push_back(log(rhobar_vap));
                this->iS.push_back(iS);
                for (unsigned int i = 0; i < N; i++)
                {
                        this->K[i].push_back(K[i]);
                        this->x[i].push_back(x[i]);
                        this->y[i].push_back(y[i]);
                        this->lnK[i].push_back(log(K[i]));
                }
        };
};
class PhaseEnvelope
{
public:
        PhaseEnvelopeLog bubble, dew;
        double rhobar_liq, rhobar_vap;
        std::vector<double> K;
        Mixture *Mix;
        void build(double p0, const std::vector<double> &z, double beta_envelope);
        void to_python_files(std::string fname);
};

class Mixture
{

public:

        unsigned int N;
        Mixture(std::vector<Fluid *> pFluids);
        Mixture(std::string FluidsString);
        ~Mixture();
        void initialize();

        double Rbar(const std::vector<double> &x);

        std::vector<Fluid *> pFluids;
        ReducingFunction * pReducing;
        ExcessTerm * pExcess;
        ResidualIdealMixture * pResidualIdealMix;

        PhaseEnvelope Envelope;

        SuccessiveSubstitutionVLE SS;
        NewtonRaphsonVLE NRVLE;

        double Wilson_lnK_factor(double T, double p, int i);

        double phir(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta_dTau(double tau, double delta, const std::vector<double> &x);
        double d2phir_dTau2(double tau, double delta, const std::vector<double> &x);
        double dphir_dDelta(double tau, double delta, const std::vector<double> &x);
        double d2phir_dDelta2(double tau, double delta, const std::vector<double> &x);
        double dphir_dTau(double tau, double delta, const std::vector<double> &x);
        double dphir_dxi(double tau, double delta, const std::vector<double> &x, int i);
        double d2phir_dxi_dTau(double tau, double delta, const std::vector<double> &x, int i);
        double d2phir_dxi_dDelta(double tau, double delta, const std::vector<double> &x, int i);
        double d2phirdxidxj(double tau, double delta, const std::vector<double> &x, int i, int j);

        double fugacity(double tau, double delta, const std::vector<double> &x, int i);
        
        double rhobar_Tpz(double T, double p, const std::vector<double> &z, double rhobar0);

        void TpzFlash(double T, double p, const std::vector<double> &z, double &rhobar, std::vector<double> &x, std::vector<double> &y);

        double g_RachfordRice(const std::vector<double> &z, const std::vector<double> &lnK, double beta);

        double dgdbeta_RachfordRice(const std::vector<double> &z, const std::vector<double> &lnK, double beta);

        double ndphir_dni__constT_V_nj(double tau, double delta, const std::vector<double> &x, int i);

        double partial_molar_volume(double tau, double delta, const std::vector<double> &x, int i);

        double dpdT__constV_n(double tau, double delta, const std::vector<double> &x, int i);

        double ndpdV__constT_n(double tau, double delta, const std::vector<double> &x, int i);

        double ndpdni__constT_V_nj(double tau, double delta, const std::vector<double> &x, int i);

        double ln_fugacity_coefficient(double tau, double delta, const std::vector<double> &x, int i);

        double dln_fugacity_coefficient_dT__constrho(double tau, double delta, const std::vector<double> &x, int i);

        double dln_fugacity_coefficient_dT__constp_n(double tau, double delta, const std::vector<double> &x, int i);

        // GERG Equation 7.42
        double dnphir_dni__constT_V_nj(double tau, double delta, const std::vector<double> &x, int i);

        double dln_fugacity_coefficient_dp__constT_n(double tau, double delta, const std::vector<double> &x, int i);

        double ndln_fugacity_coefficient_dnj__constT_p(double tau, double delta, const std::vector<double> &x, int i, int j);


        // Catch test provided
        double dln_fugacity_coefficient_dxj__constT_p_xi(double tau, double delta, const std::vector<double> &x, int i, int j);
        
        double dpdxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int j);

        double d2nphir_dni_dxj__constT_V(double tau, double delta, const std::vector<double> &x, int i, int j);

        double dphir_dxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int j);

        double d_ndphirdni_dxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int i, int j);

        double d_dphirddelta_dxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int j);

        double ddelta_dxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int j);

        double dtau_dxj__constT_V_xi(double tau, double delta, const std::vector<double> &x, int j);

        double d2nphir_dni_dT(double tau, double delta, const std::vector<double> &x, int i);

        double d_ndphirdni_dTau(double tau, double delta, const std::vector<double> &x, int i);

        double d_ndphirdni_dDelta(double tau, double delta, const std::vector<double> &x, int i);

        double nd2nphirdnidnj__constT_V(double tau, double delta, const std::vector<double> &x, int i, int j);

        double nddeltadni__constT_V_nj(double tau, double delta, const std::vector<double> &x, int i);
        double ndtaudni__constT_V_nj(double tau, double delta, const std::vector<double> &x, int i);

        double d_ndphirdni_dxj__constdelta_tau_xi(double tau, double delta, const std::vector<double> &x, int i, int j);

        double saturation_p_preconditioner(double p, const std::vector<double> &z);

        double saturation_p_Wilson(double beta, double p, const std::vector<double> &z, double T_guess, std::vector<double> &K);

        double saturation_p(double beta, double p, const std::vector<double> &z, std::vector<double> &x, std::vector<double> &y);

        double rhobar_pengrobinson(double T, double p, const std::vector<double> &x, int solution);
        
        void x_and_y_from_K(double beta, const std::vector<double> &K, const std::vector<double> &z, std::vector<double> &x, std::vector<double> &y);


        void load_excess_values(int i, int j);

        std::map<std::string, double> load_reducing_values(int i, int j);

        void check_MethaneEthane();
        void check_WaterEthanol();

        void test();
};



#endif