Air.cpp

Go to the documentation of this file.

/* Properties of Air
by Ian Bell
*/

#if defined(_MSC_VER)
#define _CRTDBG_MAP_ALLOC
#define _CRT_SECURE_NO_WARNINGS
#include <stdlib.h>
#include <crtdbg.h>
#else
#include <stdlib.h>
#endif

#include "math.h"
#include "stdio.h"
#include <string.h>
#include "CoolProp.h"
#include "FluidClass.h"
#include "Air.h"

AirClass::AirClass()
{
        static const double N[]={0,
         0.118160747229,//[1]
         0.713116392079,//[2]
        -0.161824192067e1,//[3]
         0.714140178971e-1,//[4]
        -0.865421396646e-1,//[5]
         0.134211176704,//[6]
         0.112626704218e-1,//[7]
        -0.420533228842e-1,//[8]
         0.349008431982e-1,//[9]
         0.164957183186e-3,//[10]
        -0.101365037912,//[11]
        -0.173813690970,//[12]
        -0.472103183731e-1,//[13]
        -0.122523554253e-1,//[14]
        -0.146629609713,//[15]
        -0.316055879821e-1,//[16]
         0.233594806142e-3,//[17]
         0.148287891978e-1,//[18]
        -0.938782884667e-2//[19]
        };

        static const double d[]={0,
        1,//[1]
        1,//[2]
        1,//[3]
        2,//[4]
        3,//[5]
        3,//[6]
        4,//[7]
        4,//[8]
        4,//[9]
        6,//[10]
        1,//[11]
        3,//[12]
        5,//[13]
        6,//[14]
        1,//[15]
        3,//[16]
        11,//[17]
        1,//[18]
        3//[19]
        };

        static const double t[]={0.00,
        0,//[1]
        0.33,//[2]
        1.01,//[3]
        0,//[4]
        0,//[5]
        0.15,//[6]
        0,//[7]
        0.2,//[8]
        0.35,//[9]
        1.35,//[10]
        1.6,//[11]
        0.8,//[12]
        0.95,//[13]
        1.25,//[14]
        3.6,//[15]
        6,//[16]
        3.25,//[17]
        3.5,//[18]
        15//[19]
        };

        static const double l[]={
        0,0,0,0,0,0,0,0,0,0,0, // indices [0-10]
        1,//[11]
        1,//[12]
        1,//[13]
        1,//[14]
        2,//[15]
        2,//[16]
        2,//[17]
        3,//[18]
        3,//[19]
        };

        //Constants for ideal gas expression
        static const double N0[]={0.0,
         0.605719400e-7,//[1]
        -0.210274769e-4,//[2]
        -0.158860716e-3,//[3]
        -13.841928076,//[4]
         17.275266575,//[5]
        -0.195363420e-3,//[6]
         2.490888032,//[7]
         0.791309509,//[8]
         0.212236768,//[9]
        -0.197938904,//[10]
         25.36365,//[11]
         16.90741,//[12]
         87.31279//[13]
        };

        phirlist.push_back(new phir_power(N,d,t,l,1,19,20));
        
        /*
        phi0=log(delta);
    for (k=1;k<=5;k++)
    {
        phi0=phi0+N0[k]*powInt(tau,k-4);
    }
    phi0+=N0[6]*pow(tau,1.5)+N0[7]*log(tau);
    for (k=8;k<=9;k++)
    {
        phi0+=N0[k]*log(1.0-exp(-N0[k+3]*tau));
    }
    phi0+=N0[10]*log(2.0/3.0+exp(N0[13]*tau));
        */
        std::vector<double> N0_v(N0,N0+sizeof(N0)/sizeof(double));
        std::vector<double> theta0_v(N0,N0+sizeof(N0)/sizeof(double));
        // Set some constants
        theta0_v[1]=1-4;
        theta0_v[2]=2-4;
        theta0_v[3]=3-4;
        theta0_v[4]=4-4;
        theta0_v[5]=5-4;
        theta0_v[6]=1.5;
        theta0_v[8]=N0_v[11];
        theta0_v[9]=N0_v[12];
        theta0_v[10]=N0_v[13];

        phi0list.push_back(new phi0_lead(0.0,0.0));
        phi0list.push_back(new phi0_power(N0_v,theta0_v,1,5));
        phi0list.push_back(new phi0_power(N0_v[6],1.5));
        phi0list.push_back(new phi0_logtau(N0_v[7]));
        phi0list.push_back(new phi0_Planck_Einstein(N0_v,theta0_v,8,9));
        phi0list.push_back(new phi0_Planck_Einstein2(N0_v[10],theta0_v[10],2.0/3.0));
        
        // Critical parameters (max condensing temperature)
        crit.rho = 11.8308*28.96546;
        crit.p = PressureUnit(3786.0,UNIT_KPA);
        crit.T = 132.5306;
        crit.v = 1.0/crit.rho;

        maxcondT.rho = 10.4477*28.96546;
        maxcondT.p = PressureUnit(3785.02,UNIT_KPA);
        maxcondT.T = 132.6312;
        maxcondT.v = 1.0/maxcondT.rho;

        // Other fluid parameters
        params.molemass = 28.96546;
        params.Ttriple = 59.75;
        params.ptriple = 2.44582352329;
        params.accentricfactor = 0.0816749632704;
        params.R_u = 8.31451;
        isPure = false;
        preduce = &maxcondT;

        // Limits of EOS
        limits.Tmin = 60;
        limits.Tmax = 2000.0;
        limits.pmax = 2000000.0;
        limits.rhomax = 14.21*params.molemass;
        
        EOSReference.assign("Lemmon, E.W., Jacobsen, R.T, Penoncello, S.G., and Friend, D.G.,"
                                                "\"Thermodynamic Properties of Air and Mixtures of Nitrogen, Argon, and"
                                                " Oxygen from 60 to 2000 K at Pressures to 2000 MPa,"
                                                " J. Phys. Chem. Ref. Data, 29(3):331-385, 2000.");
        TransportReference.assign("E.W. Lemmon and R. T. Jacobsen, \"Viscosity and Thermal Conductivity Equations for "
                                                          "Nitrogen, Oxygen, Argon, and Air\", "
                                                          "International Journal of Thermophysics, Vol. 25, No. 1, January 2004");

        name.assign("Air");
        aliases.push_back("air");
        aliases.push_back(std::string("AIR"));

        BibTeXKeys.EOS = "Lemmon-JPCRD-2000";
        BibTeXKeys.VISCOSITY = "Lemmon-IJT-2004";
        BibTeXKeys.CONDUCTIVITY = "Lemmon-IJT-2004";
        BibTeXKeys.SURFACE_TENSION = "Mulero-JPCRD-2012";
}

double AirClass::rhosatL(double T)
{
        const double ti[]={0,0.65,0.85,0.95,1.1};
    const double Ni[]={0,44.3413,-240.073,285.139,-88.3366,-0.892181};
    double summer=0; int k;
    for (k=1;k<=4;k++)
    {
        summer=summer+Ni[k]*pow(1.0-T/reduce.T,ti[k]);
    }
    return reduce.rho*(1+summer+Ni[5]*log(T/reduce.T));
}

double AirClass::rhosatV(double T)
{
        const double ti[]={0,0.41,1.0,2.8,6.5};
    const double Ni[]={0,-2.0466,-4.7520,-13.259,-47.652};
    double summer=0; int k;
    for (k=1;k<=4;k++)
    {
        summer=summer+Ni[k]*pow(1.0-T/reduce.T,ti[k]);
    }
    return reduce.rho*exp(summer);
}

double AirClass::psatL(double T)
{
        const double Ni[]={0,0.2260724,-7.080499,5.700283,-12.44017,17.81926,-10.81364};
        double summer=0; int k;
    for (k=1;k<=6;k++)
    {
        summer=summer+Ni[k]*pow(1-T/reduce.T,((double)k)/2.0);
    }
        double p = reduce.p.Pa*exp(reduce.T/T*summer);
        return p;
}

double AirClass::psatV(double T)
{
        const double Ni[]={0,-0.1567266,-5.539635,0,0,0.7567212,0,0,-3.514322};
        double summer=0; int k;
    for (k=1;k<=8;k++)
    {
        summer=summer+Ni[k]*pow(1-T/reduce.T,((double)k)/2.0);
    }
        double p =  reduce.p.Pa*exp(reduce.T/T*summer);
        return p;
}

double AirClass::viscosity_Trho(double T, double rho)
{
        /*
        E.W. Lemmon and R.T Jacobsen, Viscosity and Thermal Conductivity Equations for Nitrogen, Oxygen, Argon and Air
        International Journal of Thermophysics, Vol. 25, No. 1, January 2004, p.28 
        */

        double e_k=103.3, //[K]
                   sigma=0.360; //[nm]
        double eta0,etar,OMEGA,delta,tau,Tstar;
        double b[]={0.431,-0.4623,0.08406,0.005341,-0.00331};
        
        // Ideal-gas part
        Tstar = T/(e_k);
        OMEGA=exp(b[0]*powInt(log(Tstar),0)
                         +b[1]*powInt(log(Tstar),1)
                     +b[2]*powInt(log(Tstar),2)
                         +b[3]*powInt(log(Tstar),3)
                     +b[4]*powInt(log(Tstar),4));

        eta0=0.0266958*sqrt(28.9586*T)/(sigma*sigma*OMEGA);

        // Residual part
        double N[]={0,10.72,1.122,0.002019,-8.876,-0.02916};
        double t[]={0,0.2,0.05,2.4,0.6,3.6};
        double d[]={0,1,4,9,1,8};
        double l[]={0,0,0,0,1,1};
        double g[]={0,0,0,0,1,1};
        delta = rho/(10.4477*28.9586);
        tau = 132.6312/T;
        etar=N[1]*pow(tau,t[1])*pow(delta,d[1])*exp(-g[1]*pow(delta,l[1]))
                +N[2]*pow(tau,t[2])*pow(delta,d[2])*exp(-g[2]*pow(delta,l[2]))
                +N[3]*pow(tau,t[3])*pow(delta,d[3])*exp(-g[3]*pow(delta,l[3]))
                +N[4]*pow(tau,t[4])*pow(delta,d[4])*exp(-g[4]*pow(delta,l[4]))
                +N[5]*pow(tau,t[5])*pow(delta,d[5])*exp(-g[5]*pow(delta,l[5]));

        return (eta0+etar)/1e6; // uPa-s to Pa-s
}

double AirClass::X_tilde(double T,double tau,double delta)
{
        // X_tilde is dimensionless
        // Equation 11 slightly rewritten
        double drho_dp,R_Air;
    R_Air = params.R_u/params.molemass;
        drho_dp=1.0/(R_Air*T*(1+2*delta*dphir_dDelta(tau,delta)+delta*delta*d2phir_dDelta2(tau,delta)));
        return crit.p.Pa*delta/crit.rho*drho_dp;
}

double AirClass::conductivity_Trho(double T, double rho)
{
        /*
        E.W. Lemmon and R.T Jacobsen, Viscosity and Thermal Conductivity Equations for Nitrogen, Oxygen, Argon and Air
        International Journal of Thermophysics, Vol. 25, No. 1, January 2004, p.28 
        */
        double e_k=103.3, //[K]
                   sigma=0.360, //[nm]
                   Tref=265.262, //[K]
                   zeta0=0.11, //[nm]
                   LAMBDA=0.055,
                   q_D=0.31; //[nm]
        double eta0,OMEGA,delta,tau,Tstar,lambda0,lambdar,num,
                cp,cv,OMEGA_tilde,OMEGA_tilde0,zeta,nu,gamma,R0,lambdac,k,
                pi=3.141592654,mu;
        double b[]={0.431,-0.4623,0.08406,0.005341,-0.00331};

        double N[]={0,1.308,1.405,-1.036,8.743,14.76,-16.62,3.793,-6.142,-0.3778};
        double t[]={0,0,-1.1,-0.3,0.1,0.0,0.5,2.7,0.3,1.3};
        double d[]={0,0,0,0,1,2,3,7,7,11};
        double l[]={0,0,0,0,0,0,2,2,2,2};
        double g[]={0,0,0,0,0,0,1,1,1,1};
        
        delta=rho/(10.4477*28.9586);
        tau=132.6312/T;
        Tstar=T/(e_k);

        OMEGA=exp(b[0]*powInt(log(Tstar),0)
                         +b[1]*powInt(log(Tstar),1)
                     +b[2]*powInt(log(Tstar),2)
                         +b[3]*powInt(log(Tstar),3)
                     +b[4]*powInt(log(Tstar),4));

        eta0=0.0266958*sqrt(28.9586*T)/(sigma*sigma*OMEGA);
        lambda0=N[1]*eta0+N[2]*pow(tau,t[2])+N[3]*pow(tau,t[3]);

        lambdar=N[4]*pow(tau,t[4])*pow(delta,d[4])*exp(-g[4]*pow(delta,l[4]))
                   +N[5]*pow(tau,t[5])*pow(delta,d[5])*exp(-g[5]*pow(delta,l[5]))
                   +N[6]*pow(tau,t[6])*pow(delta,d[6])*exp(-g[6]*pow(delta,l[6]))
                   +N[7]*pow(tau,t[7])*pow(delta,d[7])*exp(-g[7]*pow(delta,l[7]))
                   +N[8]*pow(tau,t[8])*pow(delta,d[8])*exp(-g[8]*pow(delta,l[8]))
                   +N[9]*pow(tau,t[9])*pow(delta,d[9])*exp(-g[9]*pow(delta,l[9]));

        R0=1.01;
        nu=0.63;
        gamma=1.2415;
        k=1.380658e-23; //[J/K]

        num=X_tilde(T,crit.T/T,delta)-X_tilde(Tref,crit.T/Tref,delta)*Tref/T;

        // no critical enhancement if numerator of Eq. 10 is negative
        if (num<0)
                return (lambda0+lambdar)/1e3;
    
        cp = specific_heat_p_Trho(T,rho); //[J/kg/K]
        cv = specific_heat_v_Trho(T,rho); //[J/kg/K]
        mu = viscosity_Trho(T,rho)*1e6; //[uPa-s]

        zeta=zeta0*pow(num/LAMBDA,nu/gamma); //[nm]
        OMEGA_tilde=2.0/pi*((cp-cv)/cp*atan(zeta/q_D)+cv/cp*(zeta/q_D));
        OMEGA_tilde0=2.0/pi*(1.-exp(-1./(q_D/zeta+1.0/3.0*(zeta/q_D)*(zeta/q_D)/delta/delta)));
        lambdac=rho*cp*k*R0*T/(6*pi*zeta*mu)*(OMEGA_tilde-OMEGA_tilde0)*1e18; // 1e18 is conversion to mW/m-K (not described in paper)

        return (lambda0+lambdar+lambdac)/1e3;
}

//
//double phi0_Air(double tau, double delta)
//{
//    double phi0=0;
//    int k;
//    
//    phi0=log(delta);
//    for (k=1;k<=5;k++)
//    {
//        phi0=phi0+N0[k]*powInt(tau,k-4);
//    }
//    phi0+=N0[6]*pow(tau,1.5)+N0[7]*log(tau);
//    for (k=8;k<=9;k++)
//    {
//        phi0+=N0[k]*log(1.0-exp(-N0[k+3]*tau));
//    }
//    phi0+=N0[10]*log(2.0/3.0+exp(N0[13]*tau));
//    return phi0;
//}
//
//double dphi0_dDelta_Air(double tau, double delta)
//{
//    return 1/delta;
//}
//
//double dphi02_dDelta2_Air(double tau, double delta)
//{
//    return -1.0/powInt(delta,2);
//}
//
//double dphi0_dTau_Air(double tau, double delta)
//{
//    double dphi0_dTau=0; int k;
//    for (k=1;k<=5;k++)
//    {
//        dphi0_dTau+=N0[k]*(k-4)*powInt(tau,k-5);
//    }
//    dphi0_dTau+=1.5*N0[6]*sqrt(tau)+N0[7]/tau;
//    dphi0_dTau+=N0[8]*N0[11]/(exp(N0[11]*tau)-1)+N0[9]*N0[12]/(exp(N0[12]*tau)-1);
//    dphi0_dTau+=N0[10]*N0[13]/(2.0/3.0*exp(-N0[13]*tau)+1);
//    
//    return dphi0_dTau;
//}
//
//double dphi02_dTau2_Air(double tau, double delta)
//{
//    double dphi02_dTau2=0;
//    int k;
//    
//    for (k=1;k<=5;k++)
//    {
//        dphi02_dTau2+=N0[k]*(k-4)*(k-5)*powInt(tau,k-6);
//    }
//    dphi02_dTau2+=0.75*N0[6]*pow(tau,-0.5)-N0[7]/(tau*tau);
//    dphi02_dTau2+=-N0[8]*N0[11]*N0[11]*exp(N0[11]*tau)/powInt(exp(N0[11]*tau)-1,2);
//    dphi02_dTau2+=-N0[9]*N0[12]*N0[12]*exp(N0[12]*tau)/powInt(exp(N0[12]*tau)-1,2);
//    dphi02_dTau2+=-(2.0/3.0)*N0[10]*N0[13]*N0[13]*exp(-N0[13]*tau)/powInt((2.0/3.0)*exp(-N0[13]*tau)+1,2);
//    return dphi02_dTau2;
//}