research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Transport] Move calculations into WaterTransport class

Browse source 
Calculating viscosity and thermal conductivity in the WaterProps class
was just an unnecessary level of indirection.
This commit is contained in:
Ray Speth 2018-08-29 23:12:26 -04:00
parent ff4958a720
commit 426c2bc56e
3 changed files with 131 additions and 187 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

34
include/cantera/thermo/WaterProps.h
View file

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

149
src/thermo/WaterProps.cpp
View file

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

135
src/transport/WaterTransport.cpp
View file

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