From 426c2bc56ed21ee55ed20207252778b481351918 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 29 Aug 2018 23:12:26 -0400 Subject: [PATCH] [Transport] Move calculations into WaterTransport class Calculating viscosity and thermal conductivity in the WaterProps class was just an unnecessary level of indirection. --- include/cantera/thermo/WaterProps.h | 34 ------- src/thermo/WaterProps.cpp | 149 ---------------------------- src/transport/WaterTransport.cpp | 135 ++++++++++++++++++++++++- 3 files changed, 131 insertions(+), 187 deletions(-) diff --git a/include/cantera/thermo/WaterProps.h b/include/cantera/thermo/WaterProps.h index ec74fc263..05c289be9 100644 --- a/include/cantera/thermo/WaterProps.h +++ b/include/cantera/thermo/WaterProps.h @@ -236,40 +236,6 @@ public: */ doublereal isothermalCompressibility_IAPWS(doublereal T, doublereal P); - //! Returns the viscosity of water at the current conditions - //! (kg/m/s) - /*! - * This function calculates the value of the viscosity of pure water at the - * current T and P. - * - * The formulas used are from the paper: J. V. Sengers, J. T. R. Watson, - * "Improved International Formulations for the Viscosity and Thermal - * Conductivity of Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 - * (1986). - * - * The formulation is accurate for all temperatures and pressures, for steam - * and for water, even near the critical point. Pressures above 500 MPa and - * temperature above 900 C are suspect. - */ - doublereal viscosityWater() const; - - //! Returns the thermal conductivity of water at the current conditions - //! (W/m/K) - /*! - * This function calculates the value of the thermal conductivity of - * water at the current T and P. - * - * The formulas used are from the paper: J. V. Sengers, J. T. R. Watson, - * "Improved International Formulations for the Viscosity and Thermal - * Conductivity of Water Substance", J. Phys. Chem. Ref. Data, 15, 1291 - * (1986). - * - * The formulation is accurate for all temperatures and pressures, for steam - * and for water, even near the critical point. Pressures above 500 MPa and - * temperature above 900 C are suspect. - */ - doublereal thermalConductivityWater() const; - protected: //! Pointer to the WaterPropsIAPWS object WaterPropsIAPWS* m_waterIAPWS; diff --git a/src/thermo/WaterProps.cpp b/src/thermo/WaterProps.cpp index 0244e7b8e..01ed150b0 100644 --- a/src/thermo/WaterProps.cpp +++ b/src/thermo/WaterProps.cpp @@ -263,153 +263,4 @@ doublereal WaterProps::isothermalCompressibility_IAPWS(doublereal temp, doublere return m_waterIAPWS->isothermalCompressibility(); } -static const doublereal H[4] = {1., 0.978197, 0.579829, -0.202354}; - -static const doublereal Hij[6][7] = { - { 0.5132047, 0.2151778, -0.2818107, 0.1778064, -0.04176610, 0., 0.}, - { 0.3205656, 0.7317883, -1.070786 , 0.4605040, 0., -0.01578386, 0.}, - { 0., 1.241044 , -1.263184 , 0.2340379, 0., 0., 0.}, - { 0., 1.476783 , 0., -0.4924179, 0.1600435, 0., -0.003629481}, - {-0.7782567, 0.0 , 0., 0. , 0., 0., 0.}, - { 0.1885447, 0.0 , 0., 0. , 0., 0., 0.}, -}; - -static const doublereal rhoStar = 317.763; // kg / m3 -static const doublereal presStar = 22.115E6; // Pa - -doublereal WaterProps::viscosityWater() const -{ - static const doublereal TStar = 647.27; // Kelvin - static const doublereal muStar = 55.071E-6; //Pa s - doublereal temp = m_waterIAPWS->temperature(); - doublereal dens = m_waterIAPWS->density(); - - doublereal rhobar = dens/rhoStar; - doublereal tbar = temp / TStar; - doublereal tbar2 = tbar * tbar; - doublereal tbar3 = tbar2 * tbar; - doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); - - doublereal tfac1 = 1.0 / tbar - 1.0; - doublereal tfac2 = tfac1 * tfac1; - doublereal tfac3 = tfac2 * tfac1; - doublereal tfac4 = tfac3 * tfac1; - doublereal tfac5 = tfac4 * tfac1; - - doublereal rfac1 = rhobar - 1.0; - doublereal rfac2 = rfac1 * rfac1; - doublereal rfac3 = rfac2 * rfac1; - doublereal rfac4 = rfac3 * rfac1; - doublereal rfac5 = rfac4 * rfac1; - doublereal rfac6 = rfac5 * rfac1; - - doublereal sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 + - Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 + - Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 + - Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 + - Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 + - Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6 - ); - doublereal mu1bar = std::exp(rhobar * sum); - - // Apply the near-critical point corrections if necessary - doublereal mu2bar = 1.0; - if (tbar >= 0.9970 && tbar <= 1.0082 && rhobar >= 0.755 && rhobar <= 1.290) { - doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho(); - drhodp *= presStar / rhoStar; - doublereal xsi = rhobar * drhodp; - if (xsi >= 21.93) { - mu2bar = 0.922 * std::pow(xsi, 0.0263); - } - } - - doublereal mubar = mu0bar * mu1bar * mu2bar; - return mubar * muStar; -} - -doublereal WaterProps::thermalConductivityWater() const -{ - static const doublereal Tstar = 647.27; - static const doublereal rhostar = 317.763; - static const doublereal lambdastar = 0.4945; - static const doublereal presstar = 22.115E6; - static const doublereal L[4] = { - 1.0000, - 6.978267, - 2.599096, - -0.998254 - }; - static const doublereal Lji[6][5] = { - { 1.3293046, 1.7018363, 5.2246158, 8.7127675, -1.8525999}, - {-0.40452437, -2.2156845, -10.124111, -9.5000611, 0.93404690}, - { 0.24409490, 1.6511057, 4.9874687, 4.3786606, 0.0}, - { 0.018660751, -0.76736002, -0.27297694, -0.91783782, 0.0}, - {-0.12961068, 0.37283344, -0.43083393, 0.0, 0.0}, - { 0.044809953, -0.11203160, 0.13333849, 0.0, 0.0}, - }; - - doublereal temp = m_waterIAPWS->temperature(); - doublereal dens = m_waterIAPWS->density(); - - doublereal rhobar = dens/rhostar; - doublereal tbar = temp / Tstar; - doublereal tbar2 = tbar * tbar; - doublereal tbar3 = tbar2 * tbar; - doublereal lambda0bar = sqrt(tbar) / (L[0] + L[1]/tbar + L[2]/tbar2 + L[3]/tbar3); - - doublereal tfac1 = 1.0 / tbar - 1.0; - doublereal tfac2 = tfac1 * tfac1; - doublereal tfac3 = tfac2 * tfac1; - doublereal tfac4 = tfac3 * tfac1; - - doublereal rfac1 = rhobar - 1.0; - doublereal rfac2 = rfac1 * rfac1; - doublereal rfac3 = rfac2 * rfac1; - doublereal rfac4 = rfac3 * rfac1; - doublereal rfac5 = rfac4 * rfac1; - - doublereal sum = (Lji[0][0] + Lji[0][1]*tfac1 + Lji[0][2]*tfac2 + Lji[0][3]*tfac3 + Lji[0][4]*tfac4 + - Lji[1][0]*rfac1 + Lji[1][1]*tfac1*rfac1 + Lji[1][2]*tfac2*rfac1 + Lji[1][3]*tfac3*rfac1 + Lji[1][4]*tfac4*rfac1 + - Lji[2][0]*rfac2 + Lji[2][1]*tfac1*rfac2 + Lji[2][2]*tfac2*rfac2 + Lji[2][3]*tfac3*rfac2 + - Lji[3][0]*rfac3 + Lji[3][1]*tfac1*rfac3 + Lji[3][2]*tfac2*rfac3 + Lji[3][3]*tfac3*rfac3 + - Lji[4][0]*rfac4 + Lji[4][1]*tfac1*rfac4 + Lji[4][2]*tfac2*rfac4 + - Lji[5][0]*rfac5 + Lji[5][1]*tfac1*rfac5 + Lji[5][2]*tfac2*rfac5 - ); - doublereal lambda1bar = exp(rhobar * sum); - doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); - doublereal tfac5 = tfac4 * tfac1; - doublereal rfac6 = rfac5 * rfac1; - - sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 + - Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 + - Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 + - Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 + - Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 + - Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6 - ); - doublereal mu1bar = std::exp(rhobar * sum); - doublereal t2r2 = tbar * tbar / (rhobar * rhobar); - doublereal drhodp = 1.0 / m_waterIAPWS->dpdrho(); - drhodp *= presStar / rhoStar; - doublereal xsi = rhobar * drhodp; - doublereal xsipow = std::pow(xsi, 0.4678); - doublereal rho1 = rhobar - 1.; - doublereal rho2 = rho1 * rho1; - doublereal rho4 = rho2 * rho2; - doublereal temp2 = (tbar - 1.0) * (tbar - 1.0); - - // beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho - // - // Note for ideal gases this is equal to one. - // - // beta = delta (phi0_d() + phiR_d()) - // - tau delta (phi0_dt() + phiR_dt()) - doublereal beta = m_waterIAPWS->coeffPresExp(); - doublereal dpdT_const_rho = beta * GasConstant * dens / 18.015268; - dpdT_const_rho *= Tstar / presstar; - doublereal lambda2bar = 0.0013848 / (mu0bar * mu1bar) * t2r2 * dpdT_const_rho * dpdT_const_rho * - xsipow * sqrt(rhobar) * exp(-18.66*temp2 - rho4); - return (lambda0bar * lambda1bar + lambda2bar) * lambdastar; -} - } diff --git a/src/transport/WaterTransport.cpp b/src/transport/WaterTransport.cpp index 3919373fb..6e4f283a0 100644 --- a/src/transport/WaterTransport.cpp +++ b/src/transport/WaterTransport.cpp @@ -10,6 +10,25 @@ using namespace std; +namespace { + +const double Tstar = 647.27; +const double rhoStar = 317.763; // kg / m3 +const double presStar = 22.115E6; // Pa +const double muStar = 55.071E-6; //Pa s + +const double H[4] = {1., 0.978197, 0.579829, -0.202354}; +const double Hij[6][7] = { + { 0.5132047, 0.2151778, -0.2818107, 0.1778064, -0.04176610, 0., 0.}, + { 0.3205656, 0.7317883, -1.070786 , 0.4605040, 0., -0.01578386, 0.}, + { 0., 1.241044 , -1.263184 , 0.2340379, 0., 0., 0.}, + { 0., 1.476783 , 0., -0.4924179, 0.1600435, 0., -0.003629481}, + {-0.7782567, 0.0 , 0., 0. , 0., 0., 0.}, + { 0.1885447, 0.0 , 0., 0. , 0., 0., 0.}, +}; + +} + namespace Cantera { @@ -53,14 +72,122 @@ void WaterTransport::init(thermo_t* thermo, int mode, int log_level) } } -doublereal WaterTransport::viscosity() +double WaterTransport::viscosity() { - return m_waterProps->viscosityWater(); + static const double TStar = 647.27; // Kelvin + double temp = m_thermo->temperature(); + double dens = m_thermo->density(); + + double rhobar = dens/rhoStar; + double tbar = temp / TStar; + double tbar2 = tbar * tbar; + double tbar3 = tbar2 * tbar; + double mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); + + double tfac1 = 1.0 / tbar - 1.0; + double tfac2 = tfac1 * tfac1; + double tfac3 = tfac2 * tfac1; + double tfac4 = tfac3 * tfac1; + double tfac5 = tfac4 * tfac1; + + double rfac1 = rhobar - 1.0; + double rfac2 = rfac1 * rfac1; + double rfac3 = rfac2 * rfac1; + double rfac4 = rfac3 * rfac1; + double rfac5 = rfac4 * rfac1; + double rfac6 = rfac5 * rfac1; + + double sum = Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 + + Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 + + Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 + + Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 + + Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 + + Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6; + double mu1bar = std::exp(rhobar * sum); + + // Apply the near-critical point corrections if necessary + double mu2bar = 1.0; + if (tbar >= 0.9970 && tbar <= 1.0082 && rhobar >= 0.755 && rhobar <= 1.290) { + double drhodp = m_thermo->isothermalCompressibility() * dens; + drhodp *= presStar / rhoStar; + double xsi = rhobar * drhodp; + if (xsi >= 21.93) { + mu2bar = 0.922 * std::pow(xsi, 0.0263); + } + } + + double mubar = mu0bar * mu1bar * mu2bar; + return mubar * muStar; } -doublereal WaterTransport::thermalConductivity() +double WaterTransport::thermalConductivity() { - return m_waterProps->thermalConductivityWater(); + static const double lambdastar = 0.4945; + static const double L[4] = { + 1.0000, + 6.978267, + 2.599096, + -0.998254 + }; + static const double Lji[6][5] = { + { 1.3293046, 1.7018363, 5.2246158, 8.7127675, -1.8525999}, + {-0.40452437, -2.2156845, -10.124111, -9.5000611, 0.93404690}, + { 0.24409490, 1.6511057, 4.9874687, 4.3786606, 0.0}, + { 0.018660751, -0.76736002, -0.27297694, -0.91783782, 0.0}, + {-0.12961068, 0.37283344, -0.43083393, 0.0, 0.0}, + { 0.044809953, -0.11203160, 0.13333849, 0.0, 0.0}, + }; + + double temp = m_thermo->temperature(); + double dens = m_thermo->density(); + + double rhobar = dens / rhoStar; + double tbar = temp / Tstar; + double tbar2 = tbar * tbar; + double tbar3 = tbar2 * tbar; + double lambda0bar = sqrt(tbar) / (L[0] + L[1]/tbar + L[2]/tbar2 + L[3]/tbar3); + + double tfac1 = 1.0 / tbar - 1.0; + double tfac2 = tfac1 * tfac1; + double tfac3 = tfac2 * tfac1; + double tfac4 = tfac3 * tfac1; + double tfac5 = tfac4 * tfac1; + + double rfac1 = rhobar - 1.0; + double rfac2 = rfac1 * rfac1; + double rfac3 = rfac2 * rfac1; + double rfac4 = rfac3 * rfac1; + double rfac5 = rfac4 * rfac1; + double rfac6 = rfac5 * rfac1; + + double sum = (Lji[0][0] + Lji[0][1]*tfac1 + Lji[0][2]*tfac2 + Lji[0][3]*tfac3 + Lji[0][4]*tfac4 + + Lji[1][0]*rfac1 + Lji[1][1]*tfac1*rfac1 + Lji[1][2]*tfac2*rfac1 + Lji[1][3]*tfac3*rfac1 + Lji[1][4]*tfac4*rfac1 + + Lji[2][0]*rfac2 + Lji[2][1]*tfac1*rfac2 + Lji[2][2]*tfac2*rfac2 + Lji[2][3]*tfac3*rfac2 + + Lji[3][0]*rfac3 + Lji[3][1]*tfac1*rfac3 + Lji[3][2]*tfac2*rfac3 + Lji[3][3]*tfac3*rfac3 + + Lji[4][0]*rfac4 + Lji[4][1]*tfac1*rfac4 + Lji[4][2]*tfac2*rfac4 + + Lji[5][0]*rfac5 + Lji[5][1]*tfac1*rfac5 + Lji[5][2]*tfac2*rfac5 + ); + double lambda1bar = exp(rhobar * sum); + double mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); + + sum = (Hij[0][0] + Hij[1][0]*tfac1 + Hij[4][0]*tfac4 + Hij[5][0]*tfac5 + + Hij[0][1]*rfac1 + Hij[1][1]*tfac1*rfac1 + Hij[2][1]*tfac2*rfac1 + Hij[3][1]*tfac3*rfac1 + + Hij[0][2]*rfac2 + Hij[1][2]*tfac1*rfac2 + Hij[2][2]*tfac2*rfac2 + + Hij[0][3]*rfac3 + Hij[1][3]*tfac1*rfac3 + Hij[2][3]*tfac2*rfac3 + Hij[3][3]*tfac3*rfac3 + + Hij[0][4]*rfac4 + Hij[3][4]*tfac3*rfac4 + + Hij[1][5]*tfac1*rfac5 + Hij[3][6]*tfac3*rfac6 + ); + double mu1bar = std::exp(rhobar * sum); + double t2r2 = tbar2 / (rhobar * rhobar); + double kappa = m_thermo->isothermalCompressibility(); + double xsi = rhobar * rhobar * kappa * presStar; + double xsipow = std::pow(xsi, 0.4678); + double temp2 = (tbar - 1.0) * (tbar - 1.0); + double dpdT_const_rho = m_thermo->thermalExpansionCoeff() / kappa; + dpdT_const_rho *= Tstar / presStar; + double lambda2bar = 0.0013848 / (mu0bar * mu1bar) * t2r2 * dpdT_const_rho * dpdT_const_rho * + xsipow * sqrt(rhobar) * exp(-18.66*temp2 - rfac4); + return (lambda0bar * lambda1bar + lambda2bar) * lambdastar; } }