Heptane.cpp

/* FILE: Heptane.cpp
 * DESCRIPTION:
 *  representation of substance Heptane
 *  values and functions are from
 *  "Thermodynamic Properties in SI" bu W.C. Reynolds
 * AUTHOR: jrh@stanford.edu: GCEP, Stanford University
 *
 */

#include "Heptane.h"
#include <math.h>
#include <string.h>

namespace tpx
{

/*
 * Heptane constants
 */
static const double Tmn = 182.56;   // [K] minimum temperature for which calculations are valid
static const double Tmx = 1000.0;   // [K] maximum temperature for which calculations are valid
static const double Tc=537.68;      // [K] critical temperature
static const double Roc=197.60;        // [kg/m^3] critical density
static const double To=300;            // [K] reference Temperature
static const double R=82.99504;        // [J/(kg*K)] gas constant (for this substance)
static const double Gamma=9.611604E-6;    // [??]
static const double u0=3.4058439E5;    // [] internal energy at To
static const double s0=1.1080254E3;    // [] entropy at To
static const double Tp=400;            // [K] ??
static const double Pc=2.6199E6;    // [Pa] critical pressure
static const double M=100.20;        // [kg/kmol] molar density

/*
 * array Ahept is used by the function Pp
 */
static const double Ahept[]= {
    2.246032E-3,
    2.082990E2,
    5.085746E7,
    3.566396E9,
    1.622168E9,
    1.065237E-5,
    5.987922E-1,
    7.736602,
    1.929386E5,
    5.291379E-9
};


/*
 * array F is used by  Psat
 */
static const double F[]= {
    -7.2298764,
    3.8607475E-1,
    -3.4216472,
    4.6274432E-1,
    -9.7926124,
    -4.2058094E1,
    7.5468678E1,
    3.1758992E2
};


/*
 * array D is used by the function ldens
 */
static const double D[]= {
    1.9760405E2,
    8.9451237E2,
    -1.1462908E3,
    1.7996947E3,
    -1.7250843E3,
    9.7088329E2
};


/*
 * array G is used by the function sp
 */
static const double G[]= {
    1.1925213E5,
    -7.7231363E2,
    7.4463527,
    -3.0888167E-3,
    0.0,
    0.0
};


/*
 * C returns a multiplier in each term of the sum
 * in P-2, used in conjunction with C in the function Pp
 * j is used to represent which of the values in the summation to calculate
 * j=0 is the second additive in the formula in reynolds
 * j=1 is the third...
 */
double Heptane::C(int j,double Tinverse, double T2inverse, double T3inverse, double T4inverse)
{
    switch (j) {
    case 0 :
        return    Ahept[0] * R * T -
                  Ahept[1] -
                  Ahept[2] * T2inverse +
                  Ahept[3] * T3inverse -
                  Ahept[4] * T4inverse;
    case 1 :
        return    Ahept[5] * R * T -
                  Ahept[6] -
                  Ahept[7] * Tinverse;
    case 2 :
        return    Ahept[9] * (Ahept[6] + Ahept[7] * Tinverse);
    case 3 :
        return    Ahept[8] * T2inverse;
    default :
        return    0.0;
    }
}


/* cprime
 * derivative of C(i)
 */
inline double Heptane::Cprime(int j, double T2inverse, double T3inverse, double T4inverse)
{
    switch (j) {
    case 0 :
        return    Ahept[0] * R -
                  -2 * Ahept[2] * T3inverse +
                  -3 * Ahept[3] * T4inverse -
                  -4 * Ahept[4] * pow(T, -5.0);
    case 1 :
        return    Ahept[5] * R  -
                  -1 * Ahept[7] * T2inverse;
    case 2 :
        return    Ahept[9] * (-1 * Ahept[7] * T2inverse);
    case 3 :
        return    -2 * Ahept[8] * T3inverse;
    default :
        return 0.0;
    }
}


/*
 * I = integral from o-rho { 1/(rho^2) * H(i, rho) d rho }
 * ( see section 2 of Reynolds TPSI )
 */
inline double Heptane::I(int j, double ergho, double Gamma)
{
    switch (j) {
    case 0:
        return Rho;
    case 1:
        return Rho * Rho / 2;
    case 2:
        return pow(Rho, 5.0)/ 5;
    case 3:
        return 1 / Gamma - (Gamma * Rho * Rho + 2) * ergho / (2 * Gamma);
    default:
        return 0.0;
    }
}


/* H returns a multiplier in each term of the sum
 * in P-2
 * this is used in conjunction with C in the function Pp
 * this represents the product rho^n
 * i=0 is the second additive in the formula in reynolds
 * i=1 is the third ...
 */
double Heptane::H(int i, double egrho)
{
    if (i < 2) {
        return pow(Rho,i+2);
    } else if (i == 2) {
        return pow(Rho,6.0);
    } else if (i == 3) {
        return pow(Rho,3) * (1 + Gamma * Rho * Rho) * egrho;
    } else {
        return 0;
    }
}


/*
 * internal energy
 *  see Reynolds eqn (15) section 2
 *  u = (the integral from T to To of co(T)dT) +
 *         sum from i to N ([C(i) - T*Cprime(i)] + uo
 */
double Heptane::up()
{
    double Tinverse = 1.0/T;
    double T2inverse = pow(T, -2);
    double T3inverse = pow(T, -3);
    double T4inverse = pow(T, -4);
    double egrho = exp(-Gamma*Rho*Rho);

    double sum = 0.0;
    int i;
    for (i=1; i<=5; i++) {
        sum += G[i]*(pow(T,i) - pow(To,i))/double(i);
    }

    sum += G[0]*log(T/To);

    for (i=0; i<=6; i++) {
        sum += (C(i, Tinverse, T2inverse, T3inverse, T4inverse) - T*Cprime(i,T2inverse, T3inverse, T4inverse))*I(i,egrho, Gamma);
    }

    sum += u0;

    return sum + m_energy_offset;
}


/*
 *  entropy
 *  see Reynolds eqn (16) section 2
 */
double Heptane::sp()
{
    double T2inverse = pow(T, -2);
    double T3inverse = pow(T, -3);
    double T4inverse = pow(T, -4);
    double egrho = exp(-Gamma*Rho*Rho);

    double sum = 0.0;

    for (int i=2; i<=5; i++) {
        sum += G[i]*(pow(T,i-1) - pow(To,i-1))/double(i-1);
    }

    sum += G[1]*log(T/To);
    sum -= G[0]*(1.0/T - 1.0/To);

    for (int i=0; i<=6; i++) {
        sum -= Cprime(i,T2inverse, T3inverse, T4inverse)*I(i,egrho, Gamma);
    }

    sum += s0 - R*log(Rho);

    return sum + m_entropy_offset;
}


/*
 * Equation P-2 in Reynolds
 * P - rho - T
 * returns P (pressure)
 */
double Heptane::Pp()
{
    double Tinverse = pow(T,-1);
    double T2inverse = pow(T, -2);
    double T3inverse = pow(T, -3);
    double T4inverse = pow(T, -4);
    double egrho = exp(-Gamma*Rho*Rho);

    double P = Rho*R*T;

    for (int i=0; i<=3; i++) {
        P += C(i,Tinverse, T2inverse, T3inverse, T4inverse)*H(i,egrho);
    }

    return P;
}


/*
 * Equation S-2 in Reynolds
 * Pressure at Saturation
 */
double Heptane::Psat()
{
    double log, sum=0,P;
    if ((T < Tmn) || (T > Tc)) {

        set_Err(TempError); // Error("Heptane::Psat",TempError,T);
    }
    for (int i=1; i<=8; i++) {
        sum += F[i-1] * pow((T/Tp -1),double(i-1));
    }

    log = ((Tc/T)-1)*sum;
    P=exp(log)*Pc;

    return P;
}


/*
 * Equation D2 in Reynolds
 * liquid density, of rho_f
 */
double Heptane::ldens()
{
    double xx=1-(T/Tc), sum=0;
    if ((T < Tmn) || (T > Tc)) {
        set_Err(TempError);
    }
    for (int i=1; i<=6; i++) {
        sum+=D[i-1]*pow(xx,double(i-1)/3.0);
    }

    return sum;
}


/*
 * the following functions allow users
 * to get the properties of Heptane
 * that are not dependent on the state
 */
double Heptane::Tcrit()
{
    return Tc;
}
double Heptane::Pcrit()
{
    return Pc;
}
double Heptane::Vcrit()
{
    return 1.0/Roc;
}
double Heptane::Tmin()
{
    return Tmn;
}
double Heptane::Tmax()
{
    return Tmx;
}
char* Heptane::name()
{
    return (char*) m_name.c_str();
}
char* Heptane::formula()
{
    return (char*) m_formula.c_str();
}
double Heptane::MolWt()
{
    return M;
}
}