Bug fixes for calculating the standard state volume of

HKFT standard states.
This commit is contained in:
Harry Moffat 2008-08-30 00:56:56 +00:00
parent 0fb7da6c05
commit 23cb714cc1
3 changed files with 40 additions and 17 deletions

View file

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

View file

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

View file

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