From 23cb714cc1d70a66b37d0090d7cae9ba5bcd379f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 30 Aug 2008 00:56:56 +0000 Subject: [PATCH] Bug fixes for calculating the standard state volume of HKFT standard states. --- Cantera/src/thermo/PDSS_HKFT.cpp | 26 ++++++++++++++------------ Cantera/src/thermo/WaterPropsIAPWS.cpp | 22 +++++++++++++++++----- Cantera/src/thermo/WaterPropsIAPWS.h | 9 +++++++++ 3 files changed, 40 insertions(+), 17 deletions(-) diff --git a/Cantera/src/thermo/PDSS_HKFT.cpp b/Cantera/src/thermo/PDSS_HKFT.cpp index c1a9f96da..f49489182 100644 --- a/Cantera/src/thermo/PDSS_HKFT.cpp +++ b/Cantera/src/thermo/PDSS_HKFT.cpp @@ -200,7 +200,6 @@ namespace Cantera { doublereal PDSS_HKFT::cp_mole() const { - double pbar = m_pres * 1.0E-5; double m_presR_bar = OneAtm * 1.0E-5; @@ -274,18 +273,21 @@ namespace Cantera { doublereal PDSS_HKFT::molarVolume() const { - double a1term = m_a1; + double pbar = m_pres * 1.0E-5; + double m_presR_bar = OneAtm * 1.0E-5; + + double a1term = m_a1 * 1.0E-5; double a2term = m_a2 / (2600.E5 + m_pres); - double a3term = m_a3 / (m_temp - 228.); + double a3term = m_a3 * 1.0E-5/ (m_temp - 228.); double a4term = m_a4 / (m_temp - 228.) / (2600.E5 + m_pres); - double nu = 166027; + double nu = 166027.; double r_e_j_pr_tr = m_charge_j * m_charge_j / (m_omega_pr_tr/nu + m_charge_j/3.082); - double gval = gstar(m_temp, m_pres, 0); + double gval = gstar(m_temp, m_pres, 0); double dgvaldP = gstar(m_temp, m_pres, 3); double r_e_j = r_e_j_pr_tr + fabs(m_charge_j) * gval; @@ -309,10 +311,10 @@ namespace Cantera { double qterm = - omega_j * Q; - double molVol_calgmolbar = a1term + a2term + a3term + a4term + wterm + qterm; + double molVol_calgmolPascal = a1term + a2term + a3term + a4term + wterm + qterm; // Convert to m**3 / kmol - double molVol = molVol_calgmolbar * 4.184 / 100.; + double molVol = molVol_calgmolPascal * 4.184 * 1.0E3; return molVol; } @@ -785,10 +787,10 @@ namespace Cantera { double T1 = (TC-155.0)/300.; double fac1; - double p2 = presBar * presBar; - double p3 = presBar * p2; + double p2 = (1000. - presBar) * (1000. - presBar); + double p3 = (1000. - presBar) * p2; double p4 = p2 * p2; - double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4; + double fac2 = af_coeff[1] * p3 + af_coeff[2] * p4; if (ifunc == 0) { fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0); return fac1 * fac2; @@ -800,7 +802,7 @@ namespace Cantera { return fac1 * fac2; } else if (ifunc == 3) { fac1 = pow(T1,4.8) + af_coeff[0] * pow(T1, 16.0); - fac2 = (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5; + fac2 = - (3.0 * af_coeff[1] * p2 + 4.0 * af_coeff[2] * p3 )/ 1.0E5; return fac1 * fac2; } else { throw CanteraError("HKFT_PDSS::gg", "unimplemented"); @@ -865,7 +867,7 @@ namespace Cantera { return dgdp; } else { - throw CanteraError("HKFT_PDSS::gg", "unimplemented"); + throw CanteraError("HKFT_PDSS::g", "unimplemented"); } return 0.0; } diff --git a/Cantera/src/thermo/WaterPropsIAPWS.cpp b/Cantera/src/thermo/WaterPropsIAPWS.cpp index b1d6e5aa8..cbaa191d6 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.cpp +++ b/Cantera/src/thermo/WaterPropsIAPWS.cpp @@ -295,28 +295,36 @@ isothermalCompressibility(double temperature, double pressure) { * Difference amount is large, because we are solving for * density underneath */ - double deltaP = -0.001 * pressure; + double deltaP; double psat_at=0.0; double rhoguess = -1; int phase = -1; if (temperature > T_c) { rhoguess = pressure * M_water / (Rgas * temperature); - deltaP = +0.0001 * pressure; phase = WATER_SUPERCRIT; } else { psat_at = psat(temperature); if (pressure >= psat_at) { phase = WATER_LIQUID; - deltaP = +0.0001 * pressure; - } else + } else { phase = WATER_GAS; - deltaP = -0.0001 * pressure; + } } double dens_base = density(temperature, pressure, phase, rhoguess); if (dens_base == -1.0) { printf("problems\n"); exit(-1); } + + if (iState == WATER_GAS) { + deltaP = -0.0001 * pressure; + } else if (iState == WATER_LIQUID) { + deltaP = +0.0001 * pressure; + } else { + deltaP = +0.0001 * pressure; + } + + double pres_del = pressure + deltaP; double dens_del = density(temperature, pres_del, phase, dens_base); double Vavg = 0.5 * (1./dens_del + 1./dens_base); @@ -421,6 +429,10 @@ double WaterPropsIAPWS::psat(double temperature) { return p; } +int WaterPropsIAPWS::phaseState() const { + return iState; +} + /** * Sets the internal state of the object to the * specified temperature and density. diff --git a/Cantera/src/thermo/WaterPropsIAPWS.h b/Cantera/src/thermo/WaterPropsIAPWS.h index 4eef0d3eb..92496b9d5 100644 --- a/Cantera/src/thermo/WaterPropsIAPWS.h +++ b/Cantera/src/thermo/WaterPropsIAPWS.h @@ -349,6 +349,15 @@ public: */ double psat(double temperature); + //! Returns the Phase State flag for the current state of the object + /*! + * There are three values: + * WATER_GAS below the critical temperature but below the critical density + * WATER_LIQUID below the critical temperature but above the critical density + * WATER_CRIT above the critical temperature + */ + int phaseState() const ; + //! Returns the critical temperature of water (Kelvin) /*! * This is hard coded to the value 647.096 Kelvin