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

Internal Upgrades:

Browse source 
Added routines to calculate the spinodal curves for the water object.
This commit is contained in:
Harry Moffat 2008-09-12 21:51:04 +00:00
parent 71bcae74dd
commit 55cd65fc93
13 changed files with 565 additions and 132 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/src/thermo/PDSS.cpp
View file

@ -431,6 +431,10 @@ namespace Cantera {
err("setState_TP()");
}
void PDSS::setState_TR(doublereal temp, doublereal rho) {
err("setState_TR()");
}
/// saturation pressure
doublereal PDSS::satPressure(doublereal t){
err("satPressure()");

7
Cantera/src/thermo/PDSS.h
View file

@ -373,6 +373,13 @@ namespace Cantera {
*/
virtual void setState_TP(doublereal temp, doublereal pres);
//! Set the internal temperature and density
/*!
* @param temp Temperature (Kelvin)
* @param rho Density (kg m-3)
*/
virtual void setState_TR(doublereal temp, doublereal rho);
/**
* @}
* @name Miscellaneous properties of the standard state

14
Cantera/src/thermo/PDSS_ConstVol.cpp
View file

@ -326,7 +326,7 @@ namespace Cantera {
m_gss_RT_ptr[m_spindex] = m_hss_RT_ptr[m_spindex] - m_sss_R_ptr[m_spindex];
}
void PDSS_ConstVol::setTemperature(double temp) {
void PDSS_ConstVol::setTemperature(doublereal temp) {
m_temp = temp;
m_spthermo->update_one(m_spindex, temp,
m_cp0_R_ptr, m_h0_RT_ptr, m_s0_R_ptr);
@ -342,11 +342,21 @@ namespace Cantera {
}
void PDSS_ConstVol::setState_TP(double temp, double pres) {
void PDSS_ConstVol::setState_TP(doublereal temp, doublereal pres) {
setTemperature(temp);
setPressure(pres);
}
void PDSS_ConstVol::setState_TR(doublereal temp, doublereal rho) {
doublereal rhoStored = m_mw / m_constMolarVolume;
if (fabs(rhoStored - rho) / (rhoStored + rho) > 1.0E-4) {
throw CanteraError("PDSS_ConstVol::setState_TR",
"Inconsistent supplied rho");
}
setTemperature(temp);
}
/// saturation pressure
doublereal PDSS_ConstVol::satPressure(doublereal t){
return (1.0E-200);

11
Cantera/src/thermo/PDSS_ConstVol.h
View file

@ -298,6 +298,14 @@ namespace Cantera {
*/
virtual void setState_TP(double temp, double pres);
//! Set the internal temperature and density
/*!
* @param temp Temperature (Kelvin)
* @param rho Density (kg m-3)
*/
virtual void setState_TR(double temp, double rho);
/**
* @}
* @name Miscellaneous properties of the standard state
@ -406,6 +414,9 @@ namespace Cantera {
private:
//! Value of the constant molar volume for the species
/*!
* m3 / kmol
*/
doublereal m_constMolarVolume;
};

4
Cantera/src/thermo/PDSS_HKFT.cpp
View file

@ -273,8 +273,8 @@ namespace Cantera {
doublereal
PDSS_HKFT::molarVolume() const {
double pbar = m_pres * 1.0E-5;
double m_presR_bar = OneAtm * 1.0E-5;
// double pbar = m_pres * 1.0E-5;
//double m_presR_bar = OneAtm * 1.0E-5;
double a1term = m_a1 * 1.0E-5;

5
Cantera/src/thermo/PDSS_IdealGas.cpp
View file

@ -340,6 +340,11 @@ namespace Cantera {
setTemperature(temp);
}
void PDSS_IdealGas::setState_TR(double temp, double rho) {
m_pres = GasConstant * temp * rho / m_mw;
setTemperature(temp);
}
/// saturation pressure
doublereal PDSS_IdealGas::satPressure(doublereal t){
throw CanteraError("PDSS_IdealGas::satPressure()", "unimplemented");

13
Cantera/src/thermo/PDSS_IdealGas.h
View file

@ -300,17 +300,24 @@ namespace Cantera {
/*!
* @param temp Temperature (Kelvin)
*/
virtual void setTemperature(double temp);
virtual void setTemperature(doublereal temp);
//! Return the current storred temperature
double temperature() const;
doublereal temperature() const;
//! Set the internal temperature and pressure
/*!
* @param temp Temperature (Kelvin)
* @param pres pressure (Pascals)
*/
virtual void setState_TP(double temp, double pres);
virtual void setState_TP(doublereal temp, doublereal pres);
//! Set the internal temperature and density
/*!
* @param temp Temperature (Kelvin)
* @param rho Density (Pascals)
*/
virtual void setState_TR(doublereal temp, doublereal rho);
/**
* @}

27
Cantera/src/thermo/PDSS_Water.cpp
View file

@ -366,15 +366,19 @@ namespace Cantera {
}
void PDSS_Water::
setPressure(doublereal p) {
// In this routine we must be sure to only find the water branch of the
// curve and not the gas branch
void PDSS_Water::setPressure(doublereal p) {
doublereal T = m_temp;
doublereal dens = m_dens;
int waterState = WATER_GAS;
doublereal rc = m_sub->Rhocrit();
if (dens > rc) {
waterState = WATER_LIQUID;
int waterState = WATER_LIQUID;
if (T > m_sub->Tcrit()) {
waterState = WATER_SUPERCRIT;
}
if (p < 1.0) {
waterState = WATER_GAS;
}
#ifdef DEBUG_HKM
//printf("waterPDSS: set pres = %g t = %g, waterState = %d\n",
// p, T, waterState);
@ -388,6 +392,9 @@ namespace Cantera {
}
m_dens = dd;
m_pres = p;
m_iState = m_sub->phaseState(true);
}
// Return the volumetric thermal expansion coefficient. Units: 1/K.
@ -405,7 +412,7 @@ namespace Cantera {
doublereal PDSS_Water::dthermalExpansionCoeffdT() const {
doublereal pres = pressure();
doublereal dens_save = m_dens;
double tt = m_temp - 0.04;
doublereal tt = m_temp - 0.04;
doublereal dd = m_sub->density(tt, pres, m_iState, m_dens);
if (dd < 0.0) {
throw CanteraError("PDSS_Water::dthermalExpansionCoeffdT",
@ -452,6 +459,12 @@ namespace Cantera {
setPressure(pres);
}
void PDSS_Water::setState_TR(doublereal temp, doublereal dens) {
m_temp = temp;
m_dens = dens;
m_sub->setState_TR(m_temp, m_dens);
}
// saturation pressure
doublereal PDSS_Water::satPressure(doublereal t){
doublereal pp = m_sub->psat(t);

16
Cantera/src/thermo/PDSS_Water.h
View file

@ -311,6 +311,14 @@ namespace Cantera {
*/
virtual void setState_TP(doublereal temp, doublereal pres);
//! Set the temperature and density in the object
/*!
* @param temp Temperature (Kelvin)
* @param rho Density (kg/m3)
*/
virtual void setState_TR(doublereal temp, doublereal rho);
//! Set the density of the water phase
/*!
* This is a non-virtual function because it specific
@ -493,9 +501,11 @@ namespace Cantera {
//! state of the fluid
/*!
* 0 gas
* 1 liquid
* 2 supercrit
* 0 WATER_GAS 0
* 1 WATER_LIQUID 1
* 2 WATER_SUPERCRIT 2
* 3 WATER_UNSTABLELIQUID 3
* 4 WATER_UNSTABLEGAS
*/
int m_iState;

480
Cantera/src/thermo/WaterPropsIAPWS.cpp
View file

@ -23,13 +23,13 @@
* Critical Point values of water in mks units
*/
//! Critical Temperature value (kelvin)
const double T_c = 647.096;
const doublereal T_c = 647.096;
//! Critical Pressure (Pascals)
static const double P_c = 22.064E6;
static const doublereal P_c = 22.064E6;
//! Value of the Density at the critical point (kg m-3)
const double Rho_c = 322.;
const doublereal Rho_c = 322.;
//! Molecular Weight of water that is consistent with the paper (kg kmol-1)
static const double M_water = 18.015268;
static const doublereal M_water = 18.015268;
/*
* Note, this is the Rgas value quoted in the paper. For consistency
@ -37,8 +37,16 @@ static const double M_water = 18.015268;
*
* The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper.
*/
//static const double Rgas = 8.314472E3; // Joules kmol-1 K-1
static const double Rgas = 8.314371E3; // Joules kmol-1 K-1
//static const doublereal Rgas = 8.314472E3; // Joules kmol-1 K-1
static const doublereal Rgas = 8.314371E3; // Joules kmol-1 K-1
#ifndef MAX
# define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#endif
#ifndef MIN
# define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
#endif
WaterPropsIAPWS:: WaterPropsIAPWS() :
@ -75,7 +83,7 @@ WaterPropsIAPWS::~WaterPropsIAPWS() {
}
void WaterPropsIAPWS::calcDim(double temperature, double rho) {
void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) {
tau = T_c / temperature;
delta = rho / Rho_c;
/*
@ -92,10 +100,10 @@ void WaterPropsIAPWS::calcDim(double temperature, double rho) {
}
}
double WaterPropsIAPWS::helmholtzFE() const {
double retn = m_phi->phi(tau, delta);
double temperature = T_c/tau;
double RT = Rgas * temperature;
doublereal WaterPropsIAPWS::helmholtzFE() const {
doublereal retn = m_phi->phi(tau, delta);
doublereal temperature = T_c/tau;
doublereal RT = Rgas * temperature;
return (retn * RT);
}
@ -105,10 +113,10 @@ double WaterPropsIAPWS::helmholtzFE() const {
* Temperature: kelvin
* rho: density in kg m-3
*/
double WaterPropsIAPWS::pressure() const {
double retn = m_phi->pressureM_rhoRT(tau, delta);
double rho = delta * Rho_c;
double temperature = T_c / tau;
doublereal WaterPropsIAPWS::pressure() const {
doublereal retn = m_phi->pressureM_rhoRT(tau, delta);
doublereal rho = delta * Rho_c;
doublereal temperature = T_c / tau;
return (retn * rho * Rgas * temperature/M_water);
}
@ -127,10 +135,10 @@ double WaterPropsIAPWS::pressure() const {
*
* If a problem is encountered, a negative 1 is returned.
*/
double WaterPropsIAPWS::
density(double temperature, double pressure, int phase, double rhoguess) {
doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
int phase, doublereal rhoguess) {
double deltaGuess = 0.0;
doublereal deltaGuess = 0.0;
if (rhoguess == -1.0) {
if (phase != -1) {
if (temperature > T_c) {
@ -161,11 +169,11 @@ density(double temperature, double pressure, int phase, double rhoguess) {
}
}
double p_red = pressure * M_water / (Rgas * temperature * Rho_c);
doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
deltaGuess = rhoguess / Rho_c;
setState_TR(temperature, rhoguess);
double delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
double density_retn;
doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
doublereal density_retn;
if (delta_retn >0.0) {
delta = delta_retn;
@ -186,7 +194,71 @@ density(double temperature, double pressure, int phase, double rhoguess) {
return density_retn;
}
double WaterPropsIAPWS::density() const {
doublereal WaterPropsIAPWS::density_const(doublereal pressure,
int phase, doublereal rhoguess) const {
doublereal temperature = T_c / tau;
doublereal deltaGuess = 0.0;
doublereal deltaSave = delta;
if (rhoguess == -1.0) {
if (phase != -1) {
if (temperature > T_c) {
rhoguess = pressure * M_water / (Rgas * temperature);
} else {
if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
rhoguess = pressure * M_water / (Rgas * temperature);
} else if (phase == WATER_LIQUID) {
/*
* Provide a guess about the liquid density that is
* relatively high -> convergnce from above seems robust.
*/
rhoguess = 1000.;
} else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
throw Cantera::CanteraError("WaterPropsIAPWS::density",
"Unstable Branch finder is untested");
} else {
throw Cantera::CanteraError("WaterPropsIAPWS::density",
"unknown state: " + Cantera::int2str(phase));
}
}
} else {
/*
* Assume the Gas phase initial guess, if nothing is
* specified to the routine
*/
rhoguess = pressure * M_water / (Rgas * temperature);
}
}
doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
deltaGuess = rhoguess / Rho_c;
delta = deltaGuess;
m_phi->tdpolycalc(tau, delta);
// setState_TR(temperature, rhoguess);
doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
doublereal density_retn;
if (delta_retn > 0.0) {
delta = delta_retn;
/*
* Dimensionalize the density before returning
*/
density_retn = delta_retn * Rho_c;
} else {
density_retn = -1.0;
}
delta = deltaSave;
m_phi->tdpolycalc(tau, delta);
return density_retn;
}
doublereal WaterPropsIAPWS::density() const {
return (delta * Rho_c);
}
@ -201,9 +273,9 @@ double WaterPropsIAPWS::density() const {
* return:
* psat (Pascals)
*/
double WaterPropsIAPWS::psat_est(double temperature) {
doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const {
static const double A[8] = {
static const doublereal A[8] = {
-7.8889166E0,
2.5514255E0,
-6.716169E0,
@ -213,20 +285,20 @@ double WaterPropsIAPWS::psat_est(double temperature) {
-148.39348E0,
48.631602E0
};
double ps;
doublereal ps;
if (temperature < 314.) {
double pl = 6.3573118E0 - 8858.843E0 / temperature
doublereal pl = 6.3573118E0 - 8858.843E0 / temperature
+ 607.56335E0 * pow(temperature, -0.6);
ps = 0.1 * exp(pl);
} else {
double v = temperature / 647.25;
double w = fabs(1.0-v);
double b = 0.0;
doublereal v = temperature / 647.25;
doublereal w = fabs(1.0-v);
doublereal b = 0.0;
for (int i = 0; i < 8; i++) {
double z = i + 1;
doublereal z = i + 1;
b += A[i] * pow(w, ((z+1.0)/2.0));
}
double q = b / v;
doublereal q = b / v;
ps = 22.093*exp(q);
}
/*
@ -240,31 +312,35 @@ double WaterPropsIAPWS::psat_est(double temperature) {
* Returns the coefficient of isothermal compressibility
* of temperature and pressure.
* kappa = - d (ln V) / dP at constant T.
*
*/
double WaterPropsIAPWS::isothermalCompressibility() const {
double retn = m_phi->dimdpdrho(tau, delta);
double temperature = T_c/tau;
double dpdrho = retn * Rgas * temperature / M_water;
double dens = delta * Rho_c;
return (1.0 / (dens * dpdrho));
doublereal WaterPropsIAPWS::isothermalCompressibility() const {
doublereal dpdrho_val = dpdrho();
doublereal dens = delta * Rho_c;
return (1.0 / (dens * dpdrho_val));
}
double WaterPropsIAPWS:: coeffPresExp() const {
double retn = m_phi->dimdpdT(tau, delta);
doublereal WaterPropsIAPWS::dpdrho() const {
doublereal retn = m_phi->dimdpdrho(tau, delta);
doublereal temperature = T_c/tau;
doublereal val = retn * Rgas * temperature / M_water;
return val;
}
doublereal WaterPropsIAPWS:: coeffPresExp() const {
doublereal retn = m_phi->dimdpdT(tau, delta);
return (retn);
}
double WaterPropsIAPWS:: coeffThermExp() const {
double kappa = isothermalCompressibility();
double beta = coeffPresExp();
double dens = delta * Rho_c;
doublereal WaterPropsIAPWS:: coeffThermExp() const {
doublereal kappa = isothermalCompressibility();
doublereal beta = coeffPresExp();
doublereal dens = delta * Rho_c;
return (kappa * dens * Rgas * beta / M_water);
}
double WaterPropsIAPWS::Gibbs() const {
double gRT = m_phi->gibbs_RT();
double temperature = T_c/tau;
doublereal WaterPropsIAPWS::Gibbs() const {
doublereal gRT = m_phi->gibbs_RT();
doublereal temperature = T_c/tau;
return (gRT * Rgas * temperature);
}
@ -273,41 +349,53 @@ double WaterPropsIAPWS::Gibbs() const {
* J kmol-1 K-1.
*/
void WaterPropsIAPWS::
corr(double temperature, double pressure, double &densLiq,
double &densGas, double &delGRT) {
corr(doublereal temperature, doublereal pressure, doublereal &densLiq,
doublereal &densGas, doublereal &delGRT) {
densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
if (densLiq <= 0.0) {
printf("error liq\n");
exit(-1);
throw Cantera::CanteraError("WaterPropsIAPWS::corr",
"Error occurred trying to find liquid density at (T,P) = "
+ Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
}
setState_TR(temperature, densLiq);
double gibbsLiqRT = m_phi->gibbs_RT();
doublereal gibbsLiqRT = m_phi->gibbs_RT();
densGas = density(temperature, pressure, WATER_GAS, densGas);
if (densGas <= 0.0) {
printf("error gas\n");
exit(-1);
throw Cantera::CanteraError("WaterPropsIAPWS::corr",
"Error occurred trying to find gas density at (T,P) = "
+ Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
}
setState_TR(temperature, densGas);
double gibbsGasRT = m_phi->gibbs_RT();
doublereal gibbsGasRT = m_phi->gibbs_RT();
delGRT = gibbsLiqRT - gibbsGasRT;
}
void WaterPropsIAPWS::
corr1(double temperature, double pressure, double &densLiq,
double &densGas, double &pcorr) {
corr1(doublereal temperature, doublereal pressure, doublereal &densLiq,
doublereal &densGas, doublereal &pcorr) {
densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
if (densLiq <= 0.0) {
throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
"Error occurred trying to find liquid density at (T,P) = "
+ Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
}
setState_TR(temperature, densLiq);
double prL = m_phi->phiR();
doublereal prL = m_phi->phiR();
densGas = density(temperature, pressure, WATER_GAS, densGas);
if (densGas <= 0.0) {
throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
"Error occurred trying to find gas density at (T,P) = "
+ Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
}
setState_TR(temperature, densGas);
double prG = m_phi->phiR();
doublereal prG = m_phi->phiR();
double rhs = (prL - prG) + log(densLiq/densGas);
doublereal rhs = (prL - prG) + log(densLiq/densGas);
rhs /= (1.0/densGas - 1.0/densLiq);
pcorr = rhs * Rgas * temperature / M_water;
@ -318,15 +406,15 @@ corr1(double temperature, double pressure, double &densLiq,
* p : Pascals : Newtons/m**2
*/
static int method = 1;
double WaterPropsIAPWS::psat(double temperature) {
double densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
double dp, pcorr;
double p = psat_est(temperature);
doublereal WaterPropsIAPWS::psat(doublereal temperature) {
doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
doublereal dp, pcorr;
doublereal p = psat_est(temperature);
bool conv = false;
for (int i = 0; i < 30; i++) {
if (method == 1) {
corr(temperature, p, densLiq, densGas, delGRT);
double delV = M_water * (1.0/densLiq - 1.0/densGas);
doublereal delV = M_water * (1.0/densLiq - 1.0/densGas);
dp = - delGRT * Rgas * temperature / delV;
} else {
corr1(temperature, p, densLiq, densGas, pcorr);
@ -347,15 +435,245 @@ double WaterPropsIAPWS::psat(double temperature) {
return p;
}
int WaterPropsIAPWS::phaseState() const {
int WaterPropsIAPWS::phaseState(bool checkState) const {
if (checkState) {
if (tau <= 1.0) {
iState = WATER_SUPERCRIT;
} else {
doublereal T = T_c / tau;
doublereal rho = delta * Rho_c;
//doublereal psatTable = psat_est(T);
doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3);
doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15);
int iStateGuess = WATER_LIQUID;
if (rho < rhoMid) {
iStateGuess = WATER_GAS;
}
doublereal kappa = isothermalCompressibility();
if (kappa >= 0.0) {
iState = iStateGuess;
} else {
// When we are here we are between the spinodal curves
doublereal rhoDel = rho * 1.000001;
//setState_TR(T, rhoDel);
doublereal deltaSave = delta;
doublereal deltaDel = rhoDel / Rho_c;
delta = deltaDel;
m_phi->tdpolycalc(tau, deltaDel);
doublereal kappaDel = isothermalCompressibility();
doublereal d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho);
if (d2rhodp2 > 0.0) {
iState = WATER_UNSTABLELIQUID;
} else {
iState = WATER_UNSTABLEGAS;
}
//setState_TR(T, rho);
delta = deltaSave;
m_phi->tdpolycalc(tau, delta);
}
}
}
return iState;
}
// Find the water spinodal density
doublereal WaterPropsIAPWS::densSpinodalWater() const {
doublereal temperature = T_c/tau;
doublereal delta_save = delta;
// return the critical density if we are above or even just a little below
// the critical temperature. We just don't want to worry about the critical
// point at this juncture.
if (temperature >= T_c - 0.001) {
return Rho_c;
}
doublereal p = psat_est(temperature);
doublereal rho_low = 0.0;
doublereal rho_high = 1000;
doublereal densSatLiq = density_const(p, WATER_LIQUID);
doublereal dens_old = densSatLiq;
delta = dens_old / Rho_c;
m_phi->tdpolycalc(tau, delta);
doublereal dpdrho_old = dpdrho();
if (dpdrho_old > 0.0) {
rho_high = MIN(dens_old, rho_high);
} else {
rho_low = MAX(rho_low, dens_old);
}
doublereal dens_new = densSatLiq* (1.0001);
delta = dens_new / Rho_c;
m_phi->tdpolycalc(tau, delta);
doublereal dpdrho_new = dpdrho();
if (dpdrho_new > 0.0) {
rho_high = MIN(dens_new, rho_high);
} else {
rho_low = MAX(rho_low, dens_new);
}
bool conv = false;
for (int it = 0; it < 50; it++) {
doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
if (slope >= 0.0) {
slope = MAX(slope, dpdrho_new *5.0/ dens_new);
} else {
slope = -dpdrho_new;
//slope = MIN(slope, dpdrho_new *5.0 / dens_new);
// shouldn't be here for liquid spinodal
}
doublereal delta_rho = - dpdrho_new / slope;
if (delta_rho > 0.0) {
delta_rho = MIN(delta_rho, dens_new * 0.1);
} else {
delta_rho = MAX(delta_rho, - dens_new * 0.1);
}
doublereal dens_est = dens_new + delta_rho;
if (dens_est < rho_low) {
dens_est = 0.5 * (rho_low + dens_new);
}
if (dens_est > rho_high) {
dens_est = 0.5 * (rho_high + dens_new);
}
dens_old = dens_new;
dpdrho_old = dpdrho_new;
dens_new = dens_est;
delta = dens_new / Rho_c;
m_phi->tdpolycalc(tau, delta);
dpdrho_new = dpdrho();
if (dpdrho_new > 0.0) {
rho_high = MIN(dens_new, rho_high);
} else if (dpdrho_new < 0.0) {
rho_low = MAX(rho_low, dens_new);
} else {
conv = true;
break;
}
if (fabs(dpdrho_new) < 1.0E-5) {
conv = true;
break;
}
}
if (!conv) {
throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalWater()",
" convergence failure");
}
// Restore the original delta
delta = delta_save;
m_phi->tdpolycalc(tau, delta);
return dens_new;
}
// Find the steam spinodal density
doublereal WaterPropsIAPWS::densSpinodalSteam() const {
doublereal temperature = T_c/tau;
doublereal delta_save = delta;
// return the critical density if we are above or even just a little below
// the critical temperature. We just don't want to worry about the critical
// point at this juncture.
if (temperature >= T_c - 0.001) {
return Rho_c;
}
doublereal p = psat_est(temperature);
doublereal rho_low = 0.0;
doublereal rho_high = 1000;
doublereal densSatGas = density_const(p, WATER_GAS);
doublereal dens_old = densSatGas;
delta = dens_old / Rho_c;
m_phi->tdpolycalc(tau, delta);
doublereal dpdrho_old = dpdrho();
if (dpdrho_old < 0.0) {
rho_high = MIN(dens_old, rho_high);
} else {
rho_low = MAX(rho_low, dens_old);
}
doublereal dens_new = densSatGas * (0.99);
delta = dens_new / Rho_c;
m_phi->tdpolycalc(tau, delta);
doublereal dpdrho_new = dpdrho();
if (dpdrho_new < 0.0) {
rho_high = MIN(dens_new, rho_high);
} else {
rho_low = MAX(rho_low, dens_new);
}
bool conv = false;
for (int it = 0; it < 50; it++) {
doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
if (slope >= 0.0) {
slope = dpdrho_new;
//slope = MAX(slope, dpdrho_new *5.0/ dens_new);
// shouldn't be here for gas spinodal
} else {
//slope = -dpdrho_new;
slope = MIN(slope, dpdrho_new *5.0 / dens_new);
}
doublereal delta_rho = - dpdrho_new / slope;
if (delta_rho > 0.0) {
delta_rho = MIN(delta_rho, dens_new * 0.1);
} else {
delta_rho = MAX(delta_rho, - dens_new * 0.1);
}
doublereal dens_est = dens_new + delta_rho;
if (dens_est < rho_low) {
dens_est = 0.5 * (rho_low + dens_new);
}
if (dens_est > rho_high) {
dens_est = 0.5 * (rho_high + dens_new);
}
dens_old = dens_new;
dpdrho_old = dpdrho_new;
dens_new = dens_est;
delta = dens_new / Rho_c;
m_phi->tdpolycalc(tau, delta);
dpdrho_new = dpdrho();
if (dpdrho_new < 0.0) {
rho_high = MIN(dens_new, rho_high);
} else if (dpdrho_new > 0.0) {
rho_low = MAX(rho_low, dens_new);
} else {
conv = true;
break;
}
if (fabs(dpdrho_new) < 1.0E-5) {
conv = true;
break;
}
}
if (!conv) {
throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalSteam()",
" convergence failure");
}
// Restore the original delta
delta = delta_save;
m_phi->tdpolycalc(tau, delta);
return dens_new;
}
/**
* Sets the internal state of the object to the
* specified temperature and density.
*/
void WaterPropsIAPWS::setState_TR(double temperature, double rho) {
void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho) {
calcDim(temperature, rho);
m_phi->tdpolycalc(tau, delta);
}
@ -365,9 +683,9 @@ void WaterPropsIAPWS::setState_TR(double temperature, double rho) {
* Calculate the enthalpy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::enthalpy() const {
double temperature = T_c/tau;
double hRT = m_phi->enthalpy_RT();
doublereal WaterPropsIAPWS::enthalpy() const {
doublereal temperature = T_c/tau;
doublereal hRT = m_phi->enthalpy_RT();
return (hRT * Rgas * temperature);
}
@ -376,9 +694,9 @@ double WaterPropsIAPWS::enthalpy() const {
* Calculate the internal Energy in mks units of
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::intEnergy() const{
double temperature = T_c / tau;
double uRT = m_phi->intEnergy_RT();
doublereal WaterPropsIAPWS::intEnergy() const {
doublereal temperature = T_c / tau;
doublereal uRT = m_phi->intEnergy_RT();
return (uRT * Rgas * temperature);
}
@ -386,8 +704,8 @@ double WaterPropsIAPWS::intEnergy() const{
* Calculate the enthalpy in mks units of356
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::entropy() const {
double sR = m_phi->entropy_R();
doublereal WaterPropsIAPWS::entropy() const {
doublereal sR = m_phi->entropy_R();
return (sR * Rgas);
}
@ -395,17 +713,17 @@ double WaterPropsIAPWS::entropy() const {
* Calculate heat capacity at constant volume
* J kmol-1 K-1.
*/
double WaterPropsIAPWS::cv() const {
double cvR = m_phi->cv_R();
doublereal WaterPropsIAPWS::cv() const {
doublereal cvR = m_phi->cv_R();
return (cvR * Rgas);
}
double WaterPropsIAPWS::cp() const {
double cpR = m_phi->cp_R();
doublereal WaterPropsIAPWS::cp() const {
doublereal cpR = m_phi->cp_R();
return (cpR * Rgas);
}
double WaterPropsIAPWS::molarVolume() const {
double rho = delta * Rho_c;
doublereal WaterPropsIAPWS::molarVolume() const {
doublereal rho = delta * Rho_c;
return (M_water / rho);
}

110
Cantera/src/thermo/WaterPropsIAPWS.h
View file

@ -17,7 +17,7 @@
#define WATERPROPSIAPWS_H
#include "WaterPropsIAPWSphi.h"
#include "config.h"
/**
* @name Names for the phase regions
*
@ -163,37 +163,37 @@ public:
* @param temperature temperature (kelvin)
* @param rho density (kg m-3)
*/
void setState_TR(double temperature, double rho);
void setState_TR(doublereal temperature, doublereal rho);
//! Calculate the Helmholtz free energy in mks units of J kmol-1 K-1,
//! using the last temperature and density
double helmholtzFE() const;
doublereal helmholtzFE() const;
//! Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
//! using the last temperature and density
double Gibbs() const;
doublereal Gibbs() const;
//! Calculate the enthalpy in mks units of J kmol-1
//! using the last temperature and density
double enthalpy() const;
doublereal enthalpy() const;
//! Calculate the internal energy in mks units of J kmol-1
double intEnergy() const;
doublereal intEnergy() const;
//! Calculate the entropy in mks units of J kmol-1 K-1
double entropy() const;
doublereal entropy() const;
//! Calculate the constant volume heat capacity in mks units of J kmol-1 K-1
//! at the last temperature and density
double cv() const;
doublereal cv() const;
//! Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
//! at the last temperature and density
double cp() const;
doublereal cp() const;
//! Calculate the molar volume (kmol m-3)
//! at the last temperature and density
double molarVolume() const;
doublereal molarVolume() const;
//! Calculates the pressure (Pascals), given the current value of the
//! temperature and density.
@ -203,7 +203,7 @@ public:
* @return
* returns the pressure (Pascal)
*/
double pressure() const;
doublereal pressure() const;
//! Calculates the density given the temperature and the pressure,
//! and a guess at the density. Sets the internal state.
@ -231,14 +231,41 @@ public:
* Returns the density. If an error is encountered in the calculation
* the value of -1.0 is returned.
*/
double density(double temperature, double pressure,
int phase = -1, double rhoguess = -1.0);
doublereal density(doublereal temperature, doublereal pressure,
int phase = -1, doublereal rhoguess = -1.0);
//! Calculates the density given the temperature and the pressure,
//! and a guess at the density, while not changing the internal state
/*!
* Note, below T_c, this is a multivalued function.
*
* The #density() function calculates the density that is consistent with
* a particular value of the temperature and pressure. It may therefore be
* multivalued or potentially there may be no answer from this function. It therefore
* takes a phase guess and a density guess as optional parameters. If no guesses are
* supplied to density(), a gas phase guess is assumed. This may or may not be what
* is wanted. Therefore, density() should usually at leat be supplied with a phase
* guess so that it may manufacture an appropriate density guess.
* #density() manufactures the initial density guess, nondimensionalizes everything,
* and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation
* to find the density condition that matches the desired input pressure.
*
* @param pressure : Pressure in Pascals (Newton/m**2)
* @param phase : guessed phase of water
* : -1: no guessed phase
* @param rhoguess : guessed density of the water
* : -1.0 no guessed density
* @return
* Returns the density. If an error is encountered in the calculation
* the value of -1.0 is returned.
*/
doublereal density_const(doublereal pressure, int phase = -1, doublereal rhoguess = -1.0) const;
//! Returns the density (kg m-3)
/*!
* The density is an independent variable in the underlying equation of state
*/
double density() const;
doublereal density() const;
//! Returns the coefficient of thermal expansion.
/*!
@ -249,7 +276,7 @@ public:
* @return
* Returns the coefficient of thermal expansion
*/
double coeffThermExp() const;
doublereal coeffThermExp() const;
//! Returns the isochoric pressure-temperature coefficient
/*!
@ -261,7 +288,7 @@ public:
* beta = delta (phi0_d() + phiR_d())
* - tau delta (phi0_dt() + phiR_dt())
*/
double coeffPresExp() const;
doublereal coeffPresExp() const;
//! Returns the coefficient of isothermal compressibility for the
//! state of the object
@ -273,7 +300,17 @@ public:
* @return
* returns the isothermal compressibility
*/
double isothermalCompressibility() const;
doublereal isothermalCompressibility() const;
//! Returns the value of dp / drho at constant T for the
//! state of the object
/*!
* units - Joules / kg
*
* @return
* returns dpdrho
*/
doublereal dpdrho() const;
//! This function returns an estimated value for the saturation pressure.
/*!
@ -285,7 +322,7 @@ public:
* @return
* Returns the estimated saturation pressure
*/
double psat_est(double temperature);
doublereal psat_est(doublereal temperature) const;
//! This function returns the saturation pressure given the
//! temperature as an input parameter.
@ -295,34 +332,43 @@ public:
* Returns the saturation pressure
* units = Pascal
*/
double psat(double temperature);
doublereal psat(doublereal temperature);
//! Return the value of the density at the water spinodal point
//! for the current temperature.
/*!
*
*/
doublereal densSpinodalWater() const;
doublereal densSpinodalSteam() const;
//! 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
* WATER_SUPERCRIT above the critical temperature
*/
int phaseState() const ;
int phaseState(bool checkState = false) const ;
//! Returns the critical temperature of water (Kelvin)
/*!
* This is hard coded to the value 647.096 Kelvin
*/
double Tcrit() { return 647.096;}
doublereal Tcrit() { return 647.096;}
//! Returns the critical pressure of water (22.064E6 Pa)
/*!
* This is hard coded to the value of 22.064E6 pascals
*/
double Pcrit() { return 22.064E6;}
doublereal Pcrit() { return 22.064E6;}
//! Return the critical density of water (kg m-3)
/*!
* This is equal to 322 kg m-3.
*/
double Rhocrit() { return 322.;}
doublereal Rhocrit() { return 322.;}
private:
/**
@ -331,7 +377,7 @@ private:
* @param temperature input temperature (kelvin)
* @param rho density in kg m-3
*/
void calcDim(double temperature, double rho);
void calcDim(doublereal temperature, doublereal rho);
//! Utility routine in the calculation of the saturation pressure
/*!
@ -341,8 +387,8 @@ private:
* @param densGas output Density of gas
* @param delGRT output delGRT
*/
void corr(double temperature, double pressure, double &densLiq,
double &densGas, double &delGRT);
void corr(doublereal temperature, doublereal pressure, doublereal &densLiq,
doublereal &densGas, doublereal &delGRT);
//! Utility routine in the calculation of the saturation pressure
/*!
@ -352,8 +398,8 @@ private:
* @param densGas output Density of gas
* @param pcorr output corrected pressure
*/
void corr1(double temperature, double pressure, double &densLiq,
double &densGas, double &pcorr);
void corr1(doublereal temperature, doublereal pressure, doublereal &densLiq,
doublereal &densGas, doublereal &pcorr);
private:
@ -364,15 +410,15 @@ private:
/*!
* tau = T_C / T
*/
double tau;
doublereal tau;
//! Dimensionless density
/*!
* delta = rho / rho_c
*/
double delta;
mutable doublereal delta;
//! Current state of the system
int iState;
mutable int iState;
};
#endif

4
Cantera/src/thermo/WaterPropsIAPWSphi.cpp
View file

@ -1111,7 +1111,7 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) {
bool conv = false;
double deldd = dd;
double pcheck = 1.0E-30 + 1.0E-8 * p_red;
for (int n = 0; n < 100; n++) {
for (int n = 0; n < 200; n++) {
/*
* Calculate the internal polynomials, and then calculate the
* phi deriv functions needed by this routine.
@ -1155,7 +1155,7 @@ double WaterPropsIAPWSphi::dfind(double p_red, double tau, double deltaGuess) {
if (n < 10) {
dpdx = dpddelta * 1.1;
}
if (dpdx < 0.1) dpdx = 0.1;
if (dpdx < 0.001) dpdx = 0.001;
/*
* Formulate the update to reduced density using

2
Cantera/src/thermo/WaterPropsIAPWSphi.h
View file

@ -17,6 +17,8 @@
#ifndef WATERPROPSIAPWSPHI_H
#define WATERPROPSIAPWSPHI_H
#include "config.h"
/*!
* the WaterPropsIAPSWSphi class support low level calls for
* the real description of water.