Added in a model for the thermal conductivity of pure water

at all T and P.
Added in a test for the new model.
This commit is contained in:
Harry Moffat 2009-09-16 02:13:43 +00:00
parent 5009a3817f
commit 2bd7816b65
6 changed files with 239 additions and 36 deletions

View file

@ -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;
}
}

View file

@ -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:

View file

@ -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;
}
}

View file

@ -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:

View file

@ -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
---------------------------------------------------------------------------------

View file

@ -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) {