diff --git a/ext/tpx/CarbonDioxide.cpp b/ext/tpx/CarbonDioxide.cpp new file mode 100755 index 000000000..51a089085 --- /dev/null +++ b/ext/tpx/CarbonDioxide.cpp @@ -0,0 +1,360 @@ +/* FILE: CarbonDioxide.cpp + * DESCRIPTION: + * representation of substance Carbon Dioxide + * values and functions are from + * "Thermodynamic Properties in SI" bu W.C. Reynolds + * AUTHOR: me@rebeccahhunt.com: GCEP, Stanford University + * + */ + +#include "CarbonDioxide.h" +#include +#include + +namespace tpx { + +/* + * Carbon Dioxide constants + */ +static const double Tmn = 216.54; // [K] minimum temperature for which calculations are valid +static const double Tmx = 1500.0; // [K] maximum temperature for which calculations are valid +static const double Tc=304.21; // [K] critical temperature +static const double Roc=464.00; // [kg/m^3] critical density +static const double To=216.54; // [K] reference Temperature +static const double R=188.918; // [] gas constant for CO2 J/kg/K +static const double Gamma=5.0E-6; // [??] +static const double u0=3.217405E5; // [] internal energy at To +static const double s0=2.1396056E3; // [] entropy at To +static const double Tp=250; // [K] ?? +static const double Pc=7.38350E6; // [Pa] critical pressure +static const double M=44.01; // [kg/kmol] molar density + +/* + * array Acarbdi is used by the function named Pp + */ +static const double Acarbdi[]={ + 2.2488558E-1, +-1.3717965E2, +-1.4430214E4, +-2.9630491E6, +-2.0606039E8, + 4.5554393E-5, + 7.7042840E-2, + 4.0602371E1, + 4.0029509E-7, +-3.9436077E-4, + 1.2115286E-10, + 1.0783386E-7, + 4.3962336E-11, +-3.6505545E4, + 1.9490511E7, +-2.9186718E9, + 2.4358627E-2, + -3.7546530E1, + 1.1898141E4 +}; + + +/* + * array F is used by the function named Psat + */ +static const double F[]={ + -6.5412610, + -2.7914636E-1, + -3.4716202, + -3.4989637, + -1.9770948E1, + 1.3922839E2, + -2.7670389E2, + -7.0510251E3 +}; + + +/* + * array D is used by the function ldens + */ +static const double D[]={ + 4.6400009E2, + 6.7938129E2, + 1.4776836E3, + -3.1267676E3, + 3.6397656E3, + -1.3437098E3 +}; + + +/* + * array G is used by the function sp + */ +static const double G[]={ + 8.726361E3, + 1.840040E2, + 1.914025, + -1.667825E-3, + 7.305950E-7, + -1.255290E-10, + 3.2174105E5, + 2.1396056E3 + +}; + +/* + * C returns a multiplier in each term of the sum + * in P-3, 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... + * (this part does not include the multiplier rho^n) + */ + double CarbonDioxide::C(int j,double Tinverse, double T2inverse, double T3inverse, double T4inverse) { + switch(j) { + case 0 : + return Acarbdi[0]*T + + Acarbdi[1] + + Acarbdi[2] * Tinverse + + Acarbdi[3] * T2inverse + + Acarbdi[4] * T3inverse ; + case 1 : + return Acarbdi[5] *T + + Acarbdi[6] + + Acarbdi[7] * Tinverse ; + case 2 : + return Acarbdi[8]*T + Acarbdi[9]; + case 3 : + return Acarbdi[10]*T + Acarbdi[11]; + case 4 : + return Acarbdi[12]; + case 5 : + return Acarbdi[13] *T2inverse + + Acarbdi[14] *T3inverse + + Acarbdi[15] *T4inverse; + case 6 : + return Acarbdi[16] *T2inverse + + Acarbdi[17] *T3inverse + + Acarbdi[18] *T4inverse; + default : + return 0.0; + } +} + + /* cprime + * derivative of C(i) + */ +inline double CarbonDioxide::Cprime(int j, double T2inverse, double T3inverse, double T4inverse) { + + switch(j) { + case 0 : + return Acarbdi[0] + + - Acarbdi[2] * T2inverse + + -2 * Acarbdi[3] * T3inverse + + -3 * Acarbdi[4] * T4inverse ; + case 1 : + return Acarbdi[5] - + Acarbdi[7] * T2inverse; + case 2 : + return Acarbdi[8] ; + case 3 : + return Acarbdi[10] ; + case 4 : + return 0; + case 5 : + return + -2 *Acarbdi[13] *T3inverse + + -3 *Acarbdi[14] *T4inverse + + -4 *Acarbdi[15]* pow(T,-5); + case 6 : + return + -2 *Acarbdi[16] *T3inverse + + -3 *Acarbdi[17] *T4inverse + + -4 *Acarbdi[18] *pow(T,-5); + 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 CarbonDioxide::I(int j, double ergho, double Gamma) { + switch (j) { + + case 0: + return Rho; + case 1: + return pow(Rho, 2)/2; + case 2: + return pow(Rho, 3)/ 3; + case 3: + return pow(Rho, 4)/ 4; + case 4: + return pow(Rho, 5)/ 5; + case 5: + return (1 - ergho ) / double(2 * Gamma); + case 6: + return ( 1 - ergho * double( Gamma * pow(Rho,2) + double(1) ) )/ double(2 * Gamma * Gamma); + default: + return 0.0; + } +} + + +/* H returns a multiplier in each term of the sum + * in P-3 + * 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 CarbonDioxide::H(int i, double egrho) { + if (i < 5) + return pow(Rho,i+2); + else if (i == 5) + return pow(Rho,3)*egrho; + else if (i == 6) + return pow(Rho,5)*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 CarbonDioxide::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; + + // Equation C-6 integrated + sum += G[0]*log(T/To); + for (int i=1; i<=5; i++) + sum += G[i]*(pow(T,i) - pow(To,i))/double(i); + + + for (i=0; i<=6; i++) { + sum += I(i,egrho, Gamma) * + ( C(i, Tinverse, T2inverse, T3inverse, T4inverse) - T*Cprime(i,T2inverse, T3inverse, T4inverse) ); + } + + sum += u0; + return sum + m_energy_offset; + + } + +/* +* entropy + * see Reynolds eqn (16) section 2 +*/ + +double CarbonDioxide::sp() { + 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; + + 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/To - 1.0/T); + + + 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-3 in Reynolds + * P - rho - T + * returns P (pressure) + */ +double CarbonDioxide::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; + + // when i=0 we are on second sum of equation (where rho^2) + for(int i=0; i<=6; i++) { + P += C(i,Tinverse, T2inverse, T3inverse, T4inverse)*H(i,egrho); + } + return P; +} + + +/* + * Equation S-2 in Reynolds + * Pressure at Saturation + */ +double CarbonDioxide::Psat(){ + + double log, sum=0,P; + if ((T < Tmn) || (T > Tc)) { + cout << " error in Psat " << TempError << endl; + set_Err(TempError); // Error("CarbonDioxide::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; + + //cout << "Psat is returning " << P << " at T " << T << " and Pc " << Pc << " and Tp " << Tp << endl; + return P; + +} + +/* + * Equation D2 in Reynolds + * liquid density, of rho_f + */ +double CarbonDioxide::ldens() { + double xx=1-(T/Tc), sum=0; + if ((T < Tmn) || (T > Tc)) { + cout << " error in ldens " << TempError << endl; + 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 CarbonDioxide + * that are not dependent on the state + */ +double CarbonDioxide::Tcrit() {return Tc;} +double CarbonDioxide::Pcrit() {return Pc;} +double CarbonDioxide::Vcrit() {return 1.0/Roc;} +double CarbonDioxide::Tmin() {return Tmn;} +double CarbonDioxide::Tmax() {return Tmx;} +char * CarbonDioxide::name() {return "CarbonDioxide";} +char * CarbonDioxide::formula() {return "CO2";} +double CarbonDioxide::MolWt() {return M;} + +} + + + diff --git a/ext/tpx/CarbonDioxide.h b/ext/tpx/CarbonDioxide.h new file mode 100755 index 000000000..018c55649 --- /dev/null +++ b/ext/tpx/CarbonDioxide.h @@ -0,0 +1,49 @@ +#ifndef TPX_CARBONDIOXIDE_H +#define TPX_CARBONDIOXIDE_H + +#include "Sub.h" + + + +/* FILE: CarbonDioxide.h + * DESCRIPTION: + * representation of substance Carbon Dioxide + * values and functions are from + * "Thermodynamic Properties in SI" bu W.C. Reynolds + * AUTHOR: me@rebeccahhunt.com: GCEP, Stanford University + * + */ +namespace tpx { + +class CarbonDioxide : public Substance{ +public: + CarbonDioxide(){} + virtual ~CarbonDioxide() {} + + double MolWt(); + double Tcrit(); + double Pcrit(); + double Vcrit(); + double Tmin(); + double Tmax(); + char * name(); + char * formula(); + + double Pp(); + double up(); + double sp(); + double Psat(); + +private: + double ldens(); + double C(int jm, double, double, double, double); + double Cprime(int i, double, double, double); + double I(int i, double, double); + double H(int i, double egrho); + }; + +} + +#endif // ! TPX_CARBONDIOXIDE_H + + diff --git a/ext/tpx/Heptane.cpp b/ext/tpx/Heptane.cpp new file mode 100755 index 000000000..8824f2d71 --- /dev/null +++ b/ext/tpx/Heptane.cpp @@ -0,0 +1,311 @@ +/* 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 +#include + +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; + + for (int 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 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; + + 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 "Heptane";} +char * Heptane::formula() {return "C7H16";} +double Heptane::MolWt() {return M;} +} diff --git a/ext/tpx/Heptane.h b/ext/tpx/Heptane.h new file mode 100755 index 000000000..121732552 --- /dev/null +++ b/ext/tpx/Heptane.h @@ -0,0 +1,50 @@ +#ifndef TPX_HEPTANE_H +#define TPX_HEPTANE_H + +#include "Sub.h" + + + +/* FILE: Heptane.h + * DESCRIPTION: + * representation of substance Heptane + * values and functions are from + * "Thermodynamic Properties in SI" bu W.C. Reynolds + * AUTHOR: me@rebeccahhunt.com: GCEP, Stanford University + * AUTHOR: jrh@stanford.edu: GCEP, Stanford University + * + */ +namespace tpx { + +class Heptane : public Substance{ +public: + Heptane(){} + virtual ~Heptane() {} + + double MolWt(); + double Tcrit(); + double Pcrit(); + double Vcrit(); + double Tmin(); + double Tmax(); + char * name(); + char * formula(); + + double Pp(); + double up(); + double sp(); + double Psat(); + +private: + double ldens(); + double C(int jm, double, double, double, double); + double Cprime(int i, double, double, double); + double I(int i, double, double); + double H(int i, double egrho); + }; + +} + +#endif // ! TPX_HEPTANE_H + + diff --git a/ext/tpx/Makefile.in b/ext/tpx/Makefile.in index 13bb044cd..5ad2af5bd 100755 --- a/ext/tpx/Makefile.in +++ b/ext/tpx/Makefile.in @@ -7,7 +7,7 @@ do_ranlib = @DO_RANLIB@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) COBJS = Methane.o Nitrogen.o Oxygen.o Water.o Hydrogen.o RedlichKwong.o \ - lk.o Sub.o utils.o HFC134a.o + CarbonDioxide.o Heptane.o lk.o Sub.o utils.o HFC134a.o FOBJS = diff --git a/ext/tpx/subs.h b/ext/tpx/subs.h index a14139b48..1f6f06713 100755 --- a/ext/tpx/subs.h +++ b/ext/tpx/subs.h @@ -1,6 +1,8 @@ #ifndef TPX_SUBS_H #define TPX_SUBS_H +#include "CarbonDioxide.h" +#include "Heptane.h" #include "HFC134a.h" #include "Hydrogen.h" #include "Methane.h" @@ -11,3 +13,4 @@ // #include "lk.h" #endif + diff --git a/ext/tpx/utils.cpp b/ext/tpx/utils.cpp index 3c272b3f8..be81b09ac 100755 --- a/ext/tpx/utils.cpp +++ b/ext/tpx/utils.cpp @@ -26,6 +26,10 @@ namespace tpx { return new HFC134a; else if (lcname == "rk") return new RedlichKwong; + else if (lcname == "carbondioxide") + return new CarbonDioxide; + else if (lcname == "heptane") + return new Heptane; else return 0; } @@ -33,20 +37,25 @@ namespace tpx { Substance * GetSub(int isub) { if (isub == 0) return new water; - else if (isub == 1) + else if (isub == 1) return new nitrogen; - else if (isub == 2) + else if (isub == 2) return new methane; - else if (isub == 3) + else if (isub == 3) return new hydrogen; - else if (isub == 4) + else if (isub == 4) return new oxygen; - else if (isub == 5) + else if (isub == 5) return new HFC134a; - else if (isub == 6) + else if (isub == 6) return new RedlichKwong; - else + else if (isub == 7) + return new CarbonDioxide; + else if (isub == 8) + return new Heptane; + else return 0; - } + } } +