diff --git a/Cantera/src/thermo/WaterProps.cpp b/Cantera/src/thermo/WaterProps.cpp index 76a21b495..662fed741 100644 --- a/Cantera/src/thermo/WaterProps.cpp +++ b/Cantera/src/thermo/WaterProps.cpp @@ -596,6 +596,117 @@ namespace Cantera { return mubar * muStar; } + //! 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. + */ + double WaterProps::thermalConductivityWater() const { + static const double Tstar = 647.27; + static const double rhostar = 317.763; + static const double lambdastar = 0.4945; + static const double presstar = 22.115E6; + 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_waterIAPWS->temperature(); + double dens = m_waterIAPWS->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 lambdagas = lambda0bar * lambdastar * 1.0E3; + + double tfac1 = 1.0 / tbar - 1.0; + double tfac2 = tfac1 * tfac1; + double tfac3 = tfac2 * tfac1; + double tfac4 = tfac3 * tfac1; + + double rfac1 = rhobar - 1.0; + double rfac2 = rfac1 * rfac1; + double rfac3 = rfac2 * rfac1; + double rfac4 = rfac3 * rfac1; + double rfac5 = rfac4 * 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); + + double tfac5 = tfac4 * tfac1; + double 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 + ); + double mu1bar = std::exp(rhobar * sum); + + double t2r2 = tbar * tbar / (rhobar * rhobar); + double drhodp = 1.0 / m_waterIAPWS->dpdrho(); + drhodp *= presStar / rhoStar; + double xsi = rhobar * drhodp; + double xsipow = std::pow(xsi, 0.4678); + double rho1 = rhobar - 1.; + double rho2 = rho1 * rho1; + double rho4 = rho2 * rho2; + double 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()) + */ + double beta = m_waterIAPWS->coeffPresExp(); + + double dpdT_const_rho = beta * GasConstant * dens / 18.015268; + 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 - rho4); + + double lambda = ( lambda0bar * lambda1bar + lambda2bar) * lambdastar; + return lambda; + } } diff --git a/Cantera/src/thermo/WaterProps.h b/Cantera/src/thermo/WaterProps.h index 457ae9955..f3ca51cbb 100644 --- a/Cantera/src/thermo/WaterProps.h +++ b/Cantera/src/thermo/WaterProps.h @@ -316,7 +316,27 @@ namespace Cantera { * for steam and for water, even near the critical point. * Pressures above 500 MPa and temperature above 900 C are suspect. */ - double viscosityWater() const; + double 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. + */ + double thermalConductivityWater() const; + + + protected: diff --git a/Cantera/src/transport/WaterTransport.cpp b/Cantera/src/transport/WaterTransport.cpp index 75dab2e55..83ed5c1f3 100644 --- a/Cantera/src/transport/WaterTransport.cpp +++ b/Cantera/src/transport/WaterTransport.cpp @@ -106,10 +106,44 @@ namespace Cantera { } } + // 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 WaterTransport::viscosity() { + doublereal visc = m_waterProps->viscosityWater(); + return visc; + } - double WaterTransport::viscosity() { - double visc = m_waterProps->viscosityWater(); - return visc; + // 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 WaterTransport::thermalConductivity() { + doublereal lambda = m_waterProps->thermalConductivityWater(); + return lambda; } } diff --git a/Cantera/src/transport/WaterTransport.h b/Cantera/src/transport/WaterTransport.h index 4b24b3912..6e3b10b4a 100644 --- a/Cantera/src/transport/WaterTransport.h +++ b/Cantera/src/transport/WaterTransport.h @@ -103,6 +103,37 @@ namespace Cantera { */ virtual doublereal viscosity(); + + //! The bulk viscosity in Pa-s. + /*! + * The bulk viscosity is only + * non-zero in rare cases. Most transport managers either + * overload this method to return zero, or do not implement + * it, in which case an exception is thrown if called. + */ + virtual doublereal bulkViscosity() + { + return 0.0; + } + + + //! 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. + */ + virtual doublereal thermalConductivity(); + private: diff --git a/test_problems/cathermo/wtWater/output_blessed.txt b/test_problems/cathermo/wtWater/output_blessed.txt index 955c4bd41..6cd0e74c1 100644 --- a/test_problems/cathermo/wtWater/output_blessed.txt +++ b/test_problems/cathermo/wtWater/output_blessed.txt @@ -1,12 +1,12 @@ -------------------------------------------------------------------------- - T(C) MPa Phase Visc Visc(paper) - 10-6 kg/m/s -------------------------------------------------------------------------- - 25 0.1 L 890.496 890.5 - 100 0.1 L 281.807 281.9 - 100 10 L 284.457 284.5 - 250 5 L 106.405 106.4 - 250 50 L 117.43 117.5 - 350 17.5 L 66.9916 67.0 - 400 15 SC 24.9278 24.93 -------------------------------------------------------------------------- +------------------------------------------------------------------------------------ + T(C) MPa Phase Visc Visc(paper) lambda lambda(paper) + 10-6 kg/m/s 10-3 W/m/s +------------------------------------------------------------------------------------ + 25 0.1 L 890.496 890.5 607.155 607.2 + 100 0.1 L 281.807 281.9 679.062 679.1 + 100 10 L 284.457 284.5 684.47 684.5 + 250 5 L 106.405 106.4 622.487 622.7 + 250 50 L 117.43 117.5 671.917 672.1 + 350 17.5 L 66.9916 67.0 452.197 452.3 + 400 15 SC 24.9278 24.93 80.6817 80.68 +--------------------------------------------------------------------------------- diff --git a/test_problems/cathermo/wtWater/wtWater.cpp b/test_problems/cathermo/wtWater/wtWater.cpp index 8e025184b..953f91cc8 100644 --- a/test_problems/cathermo/wtWater/wtWater.cpp +++ b/test_problems/cathermo/wtWater/wtWater.cpp @@ -29,29 +29,31 @@ double tvalue(double val, double atol = 1.0E-9) { int main () { try { - + double lambda; WaterSSTP * w = new WaterSSTP("waterTPphase.xml", ""); WaterTransport *wtTran = new WaterTransport(w, 3); - printf("-------------------------------------------------------------------------\n"); - printf(" T(C) MPa Phase Visc Visc(paper) \n"); - printf(" 10-6 kg/m/s \n"); - printf("-------------------------------------------------------------------------\n"); + printf("------------------------------------------------------------------------------------\n"); + printf(" T(C) MPa Phase Visc Visc(paper) lambda lambda(paper)\n"); + printf(" 10-6 kg/m/s 10-3 W/m/s \n"); + printf("------------------------------------------------------------------------------------\n"); double T = 273.15 + 25.0; double pres = 1.0E5; w->setState_TP(T, pres); double visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 890.5\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 890.5 %13.6g 607.2\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); T = 273.15 + 100.0; pres = 1.0E5; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 281.9\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 281.9 %13.6g 679.1\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); @@ -59,40 +61,45 @@ int main () { pres = 1.0E7; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 284.5\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 284.5 %13.6g 684.5\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); T = 273.15 + 250.0; pres = 5.0E6; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 106.4\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 106.4 %13.6g 622.7\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); T = 273.15 + 250.0; pres = 5.0E7; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 117.5\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 117.5 %13.6g 672.1\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); T = 273.15 + 350.0; pres = 1.75E7; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g L %13.6g 67.0\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g L %13.6g 67.0 %13.6g 452.3\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); T = 273.15 + 400.0; pres = 1.50E7; w->setState_TP(T, pres); visc = wtTran->viscosity(); - printf("%8g %10.3g SC %13.6g 24.93\n", - T - 273.15, pres * 1.0E-6, visc * 1.0E6); + lambda = wtTran->thermalConductivity(); + printf("%8g %10.3g SC %13.6g 24.93 %13.6g 80.68\n", + T - 273.15, pres * 1.0E-6, visc * 1.0E6, lambda * 1.0E3); - printf("-------------------------------------------------------------------------\n"); + printf("---------------------------------------------------------------------------------\n"); delete w; } catch (CanteraError) {