From a3b92eb61fc3f1223f35c566c9fcf9b7cb669425 Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Fri, 14 Jan 2005 13:41:37 +0000 Subject: [PATCH] added HFC134a --- Cantera/python/Cantera/liquidvapor.py | 3 + Cantera/src/phasereport.cpp | 33 ++++- bin/mixmaster.py | 8 +- configure | 2 +- data/inputs/liquidvapor.cti | 21 ++- ext/f2c_libs/arith.h | 5 +- ext/tpx/HFC134a.cpp | 180 +++++++++++++------------- ext/tpx/HFC134a.h | 16 ++- ext/tpx/Makefile.in | 2 +- ext/tpx/subs.h | 2 +- ext/tpx/utils.cpp | 6 +- 11 files changed, 163 insertions(+), 115 deletions(-) diff --git a/Cantera/python/Cantera/liquidvapor.py b/Cantera/python/Cantera/liquidvapor.py index 5dfaffaa0..7af613d00 100644 --- a/Cantera/python/Cantera/liquidvapor.py +++ b/Cantera/python/Cantera/liquidvapor.py @@ -20,3 +20,6 @@ def Hydrogen(): def Oxygen(): return importPhase('liquidvapor.cti','oxygen') + +def HFC134a(): + return importPhase('liquidvapor.cti','hfc134a') diff --git a/Cantera/src/phasereport.cpp b/Cantera/src/phasereport.cpp index 2f23cacbe..12b857db1 100644 --- a/Cantera/src/phasereport.cpp +++ b/Cantera/src/phasereport.cpp @@ -71,15 +71,34 @@ namespace Cantera { th.getChemPotentials(mu.begin()); doublereal rt = GasConstant * th.temperature(); int k; + if (th.nSpecies() > 1) { - sprintf(p, "\n X Y Chem. Pot. / RT \n"); - s += p; - sprintf(p, " ------------- ------------ ------------\n"); - s += p; - for (k = 0; k < kk; k++) { - sprintf(p, "%18s %12.6g %12.6g %12.6g\n", - th.speciesName(k).c_str(), x[k], y[k], mu[k]/rt); + if (show_thermo) { + sprintf(p, "\n X " + " Y Chem. Pot. / RT \n"); s += p; + sprintf(p, " ------------- " + "------------ ------------\n"); + s += p; + for (k = 0; k < kk; k++) { + sprintf(p, "%18s %12.6g %12.6g %12.6g\n", + th.speciesName(k).c_str(), x[k], y[k], mu[k]/rt); + s += p; + } + } + else { + sprintf(p, "\n X" + "Y\n"); + s += p; + sprintf(p, " -------------" + " ------------\n"); + s += p; + for (k = 0; k < kk; k++) { + sprintf(p, "%18s %12.6g %12.6g\n", + th.speciesName(k).c_str(), x[k], y[k]); + s += p; + } + } } return s; } diff --git a/bin/mixmaster.py b/bin/mixmaster.py index 0d8d755dd..ae829744c 100755 --- a/bin/mixmaster.py +++ b/bin/mixmaster.py @@ -1,5 +1,3 @@ -try: - from MixMaster import MixMaster - o = MixMaster() -except: - i = input("MixMaster failed to start!") +from MixMaster import MixMaster +o = MixMaster() + diff --git a/configure b/configure index 1a538ccf4..bf65ccb09 100755 --- a/configure +++ b/configure @@ -279,7 +279,7 @@ CXX=${CXX:=g++} CC=${CC:=gcc} # C++ compiler flags -CXXFLAGS=${CXXFLAGS:="-O2 -Wall"} +CXXFLAGS=${CXXFLAGS:="-O2 -g -Wall"} # the C++ flags required for linking. Uncomment if additional flags # need to be passed to the linker. diff --git a/data/inputs/liquidvapor.cti b/data/inputs/liquidvapor.cti index d93ceb7ce..b47d6c0a8 100644 --- a/data/inputs/liquidvapor.cti +++ b/data/inputs/liquidvapor.cti @@ -5,7 +5,7 @@ # package, which in turn take most of the equations of state from the # compilation 'Thermodynamic Properties in SI', by W. C. Reynolds. - +# the substance flag corresponds to the liquid_vapor(name = "water", elements = " O H ", species = "H2O", @@ -41,6 +41,13 @@ liquid_vapor(name = "oxygen", initial_state = state(temperature = 300.0, pressure = OneAtm) ) +liquid_vapor(name = "hfc134a", + elements = " C F H ", + species = "C2F4H2", + substance_flag = 5, + initial_state = state(temperature = 300.0, + pressure = OneAtm) ) + #------------------------------------------------------------------ # Note that these species definitions are used ONLY to set the @@ -109,3 +116,15 @@ species(name = "H2", ) ) +# these thermo values result in h = 0, s = 0 for sat. liquid at 273.15 +# K. This convention is not consistent with that followed elsewhere +# in Cantera, but data on the heat of formation of c2f4h2 was not +# found. +species(name = "C2F4H2", + atoms = " C:2 F:4 H:2 ", + thermo = const_cp( + t0 = 273.15, + h0 = 23083414.8686, + s0 = 167025.466 + ) + ) diff --git a/ext/f2c_libs/arith.h b/ext/f2c_libs/arith.h index 995e5b254..508eb414f 100644 --- a/ext/f2c_libs/arith.h +++ b/ext/f2c_libs/arith.h @@ -1,3 +1,4 @@ -#define IEEE_8087 -#define Arith_Kind_ASL 1 +#define IEEE_MC68k +#define Arith_Kind_ASL 2 #define Double_Align +#define NANCHECK diff --git a/ext/tpx/HFC134a.cpp b/ext/tpx/HFC134a.cpp index 53025acbe..d1ae9916b 100755 --- a/ext/tpx/HFC134a.cpp +++ b/ext/tpx/HFC134a.cpp @@ -3,16 +3,18 @@ #include "HFC134a.h" #include -const double - M = 102.032, - Tmn = 170.0, - Tmx = 455.0, - Tc = 374.18, - Pc = 4056290.0, - Roc = 508.0, - R = 81.48885644; +namespace tpx { -const double a134[] = { + const double + M = 102.032, + Tmn = 170.0, + Tmx = 455.0, + Tc = 374.18, + Pc = 4056290.0, + Roc = 508.0, + R = 81.48885644; + + const double a134[] = { 0.5586817e-1, 0.4982230, 0.2458698e-1, @@ -34,77 +36,77 @@ const double a134[] = { 0.6995038e-2, -0.1452184e-1, -0.1285458e-3 -}; + }; -const double t134[] = { + const double t134[] = { -0.5, 0.0, 0.0, 0.0, 1.5, 1.5, 2.0, 2.0, 1.0, 3.0, 5.0, - 1.0, 5.0, 5.0, 6.0, 10.0, 10.0, 10.0, 18.0, 22.0, 50.0 -}; + 1.0, 5.0, 5.0, 6.0, 10.0, 10.0, 10.0, 18.0, 22.0, 50.0 + }; -const int d134[] = { + const int d134[] = { 2, 1, 3, 6, 6, 1, 1, 2, 5, 2, 2, 4, 1, 4, 1, 2, 4, 1, 5, 3, 10 -}; + }; -const double b134[] = { + const double b134[] = { -1.019535, - 9.047135, + 9.047135, -1.629789, -9.723916, -3.927170 -}; + }; -double HFC134a::fp() { + double HFC134a::fp() { double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, - sum4 = 0.0, sum5 = 0.0; + sum4 = 0.0, sum5 = 0.0; double tau = Tc/T; double delta = Rho/Roc; double phi0 = b134[0] + b134[1]*tau + b134[2]*log(tau) - + log(delta) + b134[3]/sqrt(tau) + b134[4]*pow(tau,-0.75); + + log(delta) + b134[3]/sqrt(tau) + b134[4]*pow(tau,-0.75); int i; for (i = 0; i<8; i++) - sum1 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); - for (i = 8; i<11; i++) - sum2 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum1 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + for (i = 8; i<11; i++) + sum2 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); for (i = 11; i<17; i++) - sum3 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum3 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); for (i = 17; i<20; i++) - sum4 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum4 += a134[i]*pow(tau,t134[i])*pow(delta,d134[i]); sum5 = a134[20]*pow(tau,t134[20])*pow(delta,d134[20]); - double phir = sum1 + exp(-delta)*sum2 + exp(-delta*delta)*sum3 - + exp(-delta*delta*delta)*sum4 - + exp(-delta*delta*delta*delta)*sum5; + double phir = sum1 + exp(-delta)*sum2 + exp(-delta*delta)*sum3 + + exp(-delta*delta*delta)*sum4 + + exp(-delta*delta*delta*delta)*sum5; return R*T*(phir + phi0); -} + } -double HFC134a::up() { + double HFC134a::up() { double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, - sum4 = 0.0, sum5 = 0.0; + sum4 = 0.0, sum5 = 0.0; double tau = Tc/T; double delta = Rho/Roc; double phi0t = b134[1]*tau + b134[2] - - 0.5*b134[3]*pow(tau,-0.5) - 0.75*b134[4]*pow(tau,-0.75); + - 0.5*b134[3]*pow(tau,-0.5) - 0.75*b134[4]*pow(tau,-0.75); int i; for (i = 0; i<8; i++) - sum1 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); - for (i = 8; i<11; i++) - sum2 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum1 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + for (i = 8; i<11; i++) + sum2 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); for (i = 11; i<17; i++) - sum3 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum3 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); for (i = 17; i<20; i++) - sum4 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); + sum4 += a134[i]*t134[i]*pow(tau,t134[i])*pow(delta,d134[i]); sum5 = a134[20]*t134[20]*pow(tau,t134[20])*pow(delta,d134[20]); - double phirt = sum1 + exp(-delta)*sum2 + exp(-delta*delta)*sum3 - + exp(-delta*delta*delta)*sum4 - + exp(-delta*delta*delta*delta)*sum5; - return R*T*(phirt + phi0t); - } + double phirt = sum1 + exp(-delta)*sum2 + exp(-delta*delta)*sum3 + + exp(-delta*delta*delta)*sum4 + + exp(-delta*delta*delta*delta)*sum5; + return R*T*(phirt + phi0t) + m_energy_offset; + } -double HFC134a::Pp(){ + double HFC134a::Pp(){ double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, - sum4 = 0.0, sum5 = 0.0; + sum4 = 0.0, sum5 = 0.0; double tau = Tc/T; double delta = Rho/Roc; @@ -112,61 +114,63 @@ double HFC134a::Pp(){ int i; for (i = 0; i<8; i++) - sum1 += a134[i]*pow(tau,t134[i])*d134[i]*pow(delta,d134[i]-1); - for (i = 8; i<11; i++) - sum2 += a134[i]*pow(tau,t134[i])*(d134[i] - delta)*pow(delta,d134[i]-1); + sum1 += a134[i]*pow(tau,t134[i])*d134[i]*pow(delta,d134[i]-1); + for (i = 8; i<11; i++) + sum2 += a134[i]*pow(tau,t134[i])*(d134[i] - delta)*pow(delta,d134[i]-1); sum2 *= exp(-delta); double dk = delta*delta; for (i = 11; i<17; i++) - sum3 += a134[i]*pow(tau,t134[i])*(d134[i] - 2.0*dk)*pow(delta,d134[i]-1); + sum3 += a134[i]*pow(tau,t134[i])*(d134[i] - 2.0*dk)*pow(delta,d134[i]-1); sum3 *= exp(-dk); dk *= delta; for (i = 17; i<20; i++) - sum4 += a134[i]*pow(tau,t134[i])*(d134[i] - 3.0*dk)*pow(delta,d134[i]-1); + sum4 += a134[i]*pow(tau,t134[i])*(d134[i] - 3.0*dk)*pow(delta,d134[i]-1); sum4 *= exp(-dk); dk *= delta; sum5 = a134[20]*pow(tau,t134[20])*(d134[20] - 4.0*dk)*pow(delta,d134[20]-1); - sum5 *= exp(-dk); - double phird = sum1 + sum2 + sum3 + sum4 + sum5; + sum5 *= exp(-dk); + double phird = sum1 + sum2 + sum3 + sum4 + sum5; return R*T*delta*delta*Roc*(phird + phi0d); - } + } -double HFC134a::Psat(){ - if ((T < Tmn) || (T > Tc)) set_Err(TempError); - double x1 = T/Tc; - double x2 = 1.0 - x1; - double f = -7.686556*x2 + 2.311791*pow(x2,1.5) - - 2.039554*x2*x2 - 3.583758*pow(x2,4); - return Pc*exp(f/x1); - } + double HFC134a::Psat(){ + if ((T < Tmn) || (T > Tc)) set_Err(TempError); + double x1 = T/Tc; + double x2 = 1.0 - x1; + double f = -7.686556*x2 + 2.311791*pow(x2,1.5) + - 2.039554*x2*x2 - 3.583758*pow(x2,4); + return Pc*exp(f/x1); + } -/* -double HFC134a::dPsatdT(){ - if ((T < Tmn) || (T > Tc)) set_Err(TempError); - double x1 = T/Tc; - double x2 = 1.0 - x1; - double f = -7.686556*x2 + 2.311791*pow(x2,1.5) - - 2.039554*x2*x2 - 3.583758*pow(x2,4); - double fp = -7.686556 + 1.5*2.311791*pow(x2,0.5) - - 2.0*2.039554*x2 - 4.0*3.583758*pow(x2,3); - return -Pc*exp(f/x1)*(fp/T + f/(x1*T)); - } -*/ + /* + double HFC134a::dPsatdT(){ + if ((T < Tmn) || (T > Tc)) set_Err(TempError); + double x1 = T/Tc; + double x2 = 1.0 - x1; + double f = -7.686556*x2 + 2.311791*pow(x2,1.5) + - 2.039554*x2*x2 - 3.583758*pow(x2,4); + double fp = -7.686556 + 1.5*2.311791*pow(x2,0.5) + - 2.0*2.039554*x2 - 4.0*3.583758*pow(x2,3); + return -Pc*exp(f/x1)*(fp/T + f/(x1*T)); + } + */ + + double HFC134a::ldens(){ + if ((T < Tmn) || (T > Tc)) set_Err(TempError); + double x1 = T/Tc; + double x2 = 1.0 - x1; + return 518.2 + 884.13*pow(x2,1.0/3.0) + 485.84*pow(x2,2.0/3.0) + + 193.29*pow(x2,10.0/3.0); + } + + double HFC134a::Tcrit() {return 374.21;} + double HFC134a::Pcrit() {return 4059280.0;} + double HFC134a::Vcrit() {return 1.0/511.95;} + double HFC134a::Tmin() {return Tmn;} + double HFC134a::Tmax() {return Tmx;} + char * HFC134a::name() {return "HFC-134a";} + char * HFC134a::formula() {return "C2F4H2";} + double HFC134a::MolWt() {return M;} -double HFC134a::ldens(){ - if ((T < Tmn) || (T > Tc)) set_Err(TempError); - double x1 = T/Tc; - double x2 = 1.0 - x1; - return 518.2 + 884.13*pow(x2,1.0/3.0) + 485.84*pow(x2,2.0/3.0) - + 193.29*pow(x2,10.0/3.0); } - - double HFC134a::Tcrit() {return 374.21;} - double HFC134a::Pcrit() {return 4059280.0;} - double HFC134a::Vcrit() {return 1.0/511.95;} - double HFC134a::Tmin() {return Tmn;} - double HFC134a::Tmax() {return Tmx;} - char * HFC134a::name() {return "HFC-134a";} - char * HFC134a::formula() {return "C2F4H2";} - double HFC134a::MolWt() {return M;} diff --git a/ext/tpx/HFC134a.h b/ext/tpx/HFC134a.h index 4a291a683..8334597b7 100755 --- a/ext/tpx/HFC134a.h +++ b/ext/tpx/HFC134a.h @@ -3,14 +3,15 @@ #include "sub.h" -class HFC134a : public Substance{ -public: +namespace tpx { + class HFC134a : public Substance{ + public: HFC134a(){} - ~HFC134a(){} + ~HFC134a(){} double MolWt(); double Tcrit(); - double Pcrit(); + double Pcrit(); double Vcrit(); double Tmin(); double Tmax(); @@ -21,11 +22,12 @@ public: double fp(); double up(); double sp() - { return (up() - fp())/T; } + { return (up() - fp())/T + m_entropy_offset; } double Psat(); -// double dPsatdT(); -private: + // double dPsatdT(); + private: double ldens(); }; +} #endif // ! HFC134_H diff --git a/ext/tpx/Makefile.in b/ext/tpx/Makefile.in index 4d7e01616..54bc0bd82 100755 --- a/ext/tpx/Makefile.in +++ b/ext/tpx/Makefile.in @@ -8,7 +8,7 @@ OBJDIR = . CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) COBJS = Methane.o Nitrogen.o Oxygen.o Water.o Hydrogen.o RedlichKwong.o \ - lk.o Sub.o utils.o + lk.o Sub.o utils.o HFC134a.o FOBJS = diff --git a/ext/tpx/subs.h b/ext/tpx/subs.h index 2e7d9b05f..a14139b48 100755 --- a/ext/tpx/subs.h +++ b/ext/tpx/subs.h @@ -1,7 +1,7 @@ #ifndef TPX_SUBS_H #define TPX_SUBS_H -// #include "HFC134a.h" +#include "HFC134a.h" #include "Hydrogen.h" #include "Methane.h" #include "Nitrogen.h" diff --git a/ext/tpx/utils.cpp b/ext/tpx/utils.cpp index 369ac1f21..93b2d4365 100755 --- a/ext/tpx/utils.cpp +++ b/ext/tpx/utils.cpp @@ -22,6 +22,8 @@ namespace tpx { return new hydrogen; else if (lcname == "oxygen") return new oxygen; + else if (lcname == "hfc134a") + return new HFC134a; else if (lcname == "rk") return new RedlichKwong; } @@ -37,8 +39,8 @@ namespace tpx { return new hydrogen; else if (isub == 4) return new oxygen; - // else if (isub == 5) - // return new HFC134a; + else if (isub == 5) + return new HFC134a; else if (isub == 6) return new RedlichKwong; else