WaterPropsIAPWS.cpp

Go to the documentation of this file.
00001 /**
00002  * @file WaterPropsIAPWS.cpp
00003  * Definitions for a class for calculating the equation of state of water
00004  * from the IAPWS 1995 Formulation based on the steam tables thermodynamic
00005  * basis (See class \link Cantera::WaterPropsIAPWS WaterPropsIAPWS\endlink).
00006  */
00007 /*
00008  * Copywrite (2006) Sandia Corporation. Under the terms of
00009  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
00010  * U.S. Government retains certain rights in this software.
00011  */
00012 /*
00013  * $Id: WaterPropsIAPWS.cpp 398 2010-02-09 20:24:11Z hkmoffa $
00014  */
00015 
00016 #include "WaterPropsIAPWS.h"
00017 #include "ctexceptions.h"
00018 #include "stringUtils.h"
00019 #include <cmath>
00020 #include <cstdio>
00021 #include <cstdlib>
00022 
00023 namespace Cantera {
00024 /*
00025  * Critical Point values of water in mks units
00026  */
00027 //! Critical Temperature value (kelvin)
00028 const doublereal T_c = 647.096;  
00029 //! Critical Pressure (Pascals)
00030 static const doublereal P_c = 22.064E6; 
00031 //! Value of the Density at the critical point (kg  m-3)
00032 const doublereal Rho_c = 322.;
00033 //! Molecular Weight of water that is consistent with the paper (kg kmol-1)
00034 static const doublereal M_water = 18.015268; 
00035 
00036 /*
00037  * Note, this is the Rgas value quoted in the paper. For consistency
00038  * we have to use that value and not the updated value
00039  *
00040  * The Ratio of R/M = 0.46151805 kJ kg-1 K-1 , which is Eqn. (6.3) in the paper.
00041  */
00042 //static const doublereal Rgas = 8.314472E3;   // Joules kmol-1 K-1
00043 static const doublereal Rgas = 8.314371E3;   // Joules kmol-1 K-1
00044 //@{
00045 #ifndef MAX
00046 # define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
00047 #endif
00048 
00049 #ifndef MIN
00050 # define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
00051 #endif
00052 //@}
00053 
00054 // Base constructor
00055 WaterPropsIAPWS:: WaterPropsIAPWS() :
00056   m_phi(0),
00057   tau(-1.0),
00058   delta(-1.0),
00059   iState(-30000)
00060 {
00061   m_phi = new WaterPropsIAPWSphi();
00062 }
00063 
00064 // Copy constructor
00065 /*
00066  * @param b Object to be copied
00067  */
00068 WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS &b) :
00069   m_phi(0),
00070   tau(b.tau),
00071   delta(b.delta),
00072   iState(b.iState)
00073 {
00074   m_phi = new WaterPropsIAPWSphi();
00075   m_phi->tdpolycalc(tau, delta);
00076 }
00077 
00078 // assignment constructor
00079 /*
00080  * @param right Object to be copied
00081  */
00082 WaterPropsIAPWS & WaterPropsIAPWS::operator=(const WaterPropsIAPWS &b) {
00083   if (this == &b) return *this;
00084   tau = b.tau;
00085   delta = b.delta;
00086   iState = b.iState;
00087   m_phi->tdpolycalc(tau, delta);
00088   return *this;
00089 }
00090 
00091 // destructor
00092 WaterPropsIAPWS::~WaterPropsIAPWS() {
00093   delete (m_phi);
00094   m_phi = 0;
00095 }
00096 
00097 /*
00098  * Calculate the dimensionless temp and rho and store internally.
00099  *
00100  * @param temperature   input temperature (kelvin)
00101  *  @param rho          density in kg m-3
00102  *
00103  *  this is a private function
00104  */
00105 void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) {
00106   tau = T_c / temperature;
00107   delta = rho / Rho_c;
00108   /*
00109    * Determine the internal state
00110    */
00111   if (temperature > T_c) {
00112       iState = WATER_SUPERCRIT;
00113   } else {
00114     if (delta < 1.0) {
00115       iState = WATER_GAS;
00116     } else {
00117       iState = WATER_LIQUID;
00118     }
00119   }
00120 }
00121 
00122 // Calculate the Helmholtz free energy in mks units of J kmol-1 K-1,
00123 // using the last temperature and density
00124 doublereal  WaterPropsIAPWS::helmholtzFE() const {
00125   doublereal retn = m_phi->phi(tau, delta);
00126   doublereal temperature = T_c/tau;
00127   doublereal RT = Rgas * temperature;
00128   return (retn * RT);
00129 }
00130 
00131 /*
00132  * Calculate the pressure (Pascals), using the 
00133  * current internally storred temperature and density
00134  *  Temperature: kelvin
00135  *  rho: density in kg m-3 
00136  */
00137 doublereal  WaterPropsIAPWS::pressure() const {
00138   doublereal retn = m_phi->pressureM_rhoRT(tau, delta);
00139   doublereal rho = delta * Rho_c;
00140   doublereal temperature = T_c / tau;
00141   return (retn * rho * Rgas * temperature/M_water);
00142 }
00143 
00144 /*
00145  * Calculates the density given the temperature and the pressure,
00146  * and a guess at the density. Note, below T_c, this is a
00147  * multivalued function. 
00148  *
00149  * parameters:
00150  *    temperature: Kelvin
00151  *    pressure   : Pressure in Pascals (Newton/m**2)
00152  *    phase      : guessed phase of water 
00153  *               : -1: no guessed phase
00154  *    rhoguess   : guessed density of the water
00155  *               : -1.0 no guessed density
00156  *
00157  * If a problem is encountered, a negative 1 is returned.
00158  */
00159 doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
00160                                     int phase, doublereal rhoguess) {
00161 
00162   doublereal deltaGuess = 0.0;
00163   if (rhoguess == -1.0) {
00164     if (phase != -1) {
00165       if (temperature > T_c) {
00166         rhoguess = pressure * M_water / (Rgas * temperature);
00167       } else {
00168         if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
00169           rhoguess = pressure * M_water / (Rgas * temperature);
00170         } else if (phase == WATER_LIQUID) {
00171           /*
00172            * Provide a guess about the liquid density that is 
00173            * relatively high -> convergnce from above seems robust.
00174            */
00175           rhoguess = 1000.;
00176         } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
00177           throw Cantera::CanteraError("WaterPropsIAPWS::density", 
00178                                       "Unstable Branch finder is untested");
00179         } else {
00180           throw Cantera::CanteraError("WaterPropsIAPWS::density", 
00181                                       "unknown state: " + Cantera::int2str(phase));
00182         }
00183       }
00184     } else {
00185       /*
00186        * Assume the Gas phase initial guess, if nothing is
00187        * specified to the routine
00188        */
00189       rhoguess = pressure * M_water / (Rgas * temperature);
00190     }
00191 
00192   }
00193   doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
00194   deltaGuess = rhoguess / Rho_c;
00195   setState_TR(temperature, rhoguess);
00196   doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
00197   doublereal density_retn;
00198   if (delta_retn >0.0) {
00199     delta = delta_retn;
00200       
00201     /*
00202      * Dimensionalize the density before returning
00203      */
00204     density_retn = delta_retn * Rho_c;
00205     /*
00206      * Set the internal state -> this may be
00207      * a duplication. However, let's just be sure.
00208      */
00209     setState_TR(temperature, density_retn);
00210     
00211 
00212   } else {
00213     density_retn = -1.0;
00214   }
00215   return density_retn;
00216 }
00217 
00218 // Calculates the density given the temperature and the pressure,
00219 // and a guess at the density, while not changing the internal state
00220 /*
00221  *  Note, below T_c, this is a multivalued function.
00222  *
00223  * The #density() function calculates the density that is consistent with
00224  * a particular value of the temperature and pressure. It may therefore be
00225  * multivalued or potentially there may be no answer from this function. It therefore
00226  * takes a phase guess and a density guess as optional parameters. If no guesses are
00227  *
00228  * supplied to density(), a gas phase guess is assumed. This may or may not be what
00229  * is wanted. Therefore, density() should usually at leat be supplied with a phase
00230  * guess so that it may manufacture an appropriate density guess.
00231  * #density() manufactures the initial density guess, nondimensionalizes everything,
00232  * and then calls #WaterPropsIAPWSphi::dfind(), which does the iterative calculation
00233  * to find the density condition that matches the desired input pressure.
00234  *
00235  *  @param  pressure   : Pressure in Pascals (Newton/m**2)
00236  *  @param  phase      : guessed phase of water
00237  *                     : -1: no guessed phase
00238  *  @param rhoguess    : guessed density of the water
00239  *                     : -1.0 no guessed density
00240  *  @return
00241  *     Returns the density. If an error is encountered in the calculation
00242  *     the value of -1.0 is returned.
00243  */
00244 doublereal WaterPropsIAPWS::density_const(doublereal pressure,
00245                                           int phase, doublereal rhoguess) const {
00246   doublereal temperature = T_c / tau;
00247   doublereal deltaGuess = 0.0;
00248   doublereal deltaSave = delta;
00249   if (rhoguess == -1.0) {
00250     if (phase != -1) {
00251       if (temperature > T_c) {
00252         rhoguess = pressure * M_water / (Rgas * temperature);
00253       } else {
00254         if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
00255           rhoguess = pressure * M_water / (Rgas * temperature);
00256         } else if (phase == WATER_LIQUID) {
00257           /*
00258            * Provide a guess about the liquid density that is 
00259            * relatively high -> convergnce from above seems robust.
00260            */
00261           rhoguess = 1000.;
00262         } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
00263           throw Cantera::CanteraError("WaterPropsIAPWS::density", 
00264                                       "Unstable Branch finder is untested");
00265         } else {
00266           throw Cantera::CanteraError("WaterPropsIAPWS::density", 
00267                                       "unknown state: " + Cantera::int2str(phase));
00268         }
00269       }
00270     } else {
00271       /*
00272        * Assume the Gas phase initial guess, if nothing is
00273        * specified to the routine
00274        */
00275       rhoguess = pressure * M_water / (Rgas * temperature);
00276     }
00277 
00278   }
00279   doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
00280   deltaGuess = rhoguess / Rho_c;
00281 
00282   delta = deltaGuess;
00283   m_phi->tdpolycalc(tau, delta);
00284   //  setState_TR(temperature, rhoguess);
00285 
00286   doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
00287   doublereal density_retn;
00288   if (delta_retn > 0.0) {
00289     delta = delta_retn;
00290       
00291     /*
00292      * Dimensionalize the density before returning
00293      */
00294     density_retn = delta_retn * Rho_c;
00295   
00296   } else {
00297     density_retn = -1.0;
00298   }
00299 
00300   delta = deltaSave;
00301   m_phi->tdpolycalc(tau, delta);
00302   return density_retn;
00303 }
00304 
00305 // Returns the density (kg m-3)
00306 /*
00307  * The density is an independent variable in the underlying equation of state
00308  *
00309  * @return  Returns the density (kg m-3)
00310  */
00311 doublereal WaterPropsIAPWS::density() const {
00312   return (delta * Rho_c);
00313 }
00314 
00315 // Returns the temperature (Kelvin)
00316 /*
00317  * @return  Returns the internally storred temperature
00318  */
00319 doublereal WaterPropsIAPWS::temperature() const {
00320   return (T_c / tau);
00321 }
00322 
00323 /*
00324  * psat_est provides a rough estimate of the saturation
00325  * pressure given the temperature. This is used as an initial
00326  * guess for refining the pressure.
00327  *
00328  * Input
00329  *   temperature (kelvin)
00330  *
00331  * return:
00332  *   psat (Pascals)
00333  */
00334 doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const {
00335     
00336   static const doublereal A[8] = {
00337     -7.8889166E0,
00338     2.5514255E0,
00339     -6.716169E0,
00340     33.2239495E0,
00341     -105.38479E0,
00342     174.35319E0,
00343     -148.39348E0,
00344     48.631602E0
00345   };
00346   doublereal ps;
00347   if (temperature < 314.) {
00348     doublereal pl = 6.3573118E0 - 8858.843E0 / temperature
00349       + 607.56335E0 * pow(temperature, -0.6);
00350     ps = 0.1 * exp(pl);
00351   } else {
00352     doublereal v = temperature / 647.25;
00353     doublereal w = fabs(1.0-v);
00354     doublereal b = 0.0;
00355     for (int i = 0; i < 8; i++) {
00356       doublereal z = i + 1;
00357       b += A[i] * pow(w, ((z+1.0)/2.0));
00358     }
00359     doublereal q = b / v;
00360     ps = 22.093*exp(q);
00361   }
00362   /*
00363    * Original correlation was in cgs. Convert to mks
00364    */
00365   ps *= 1.0E6;
00366   return ps;
00367 }
00368 
00369 /*
00370  * Returns the coefficient of isothermal compressibility
00371  * of temperature and pressure.
00372  *          kappa = - d (ln V) / dP at constant T.
00373  */
00374 doublereal WaterPropsIAPWS::isothermalCompressibility() const {
00375   doublereal dpdrho_val = dpdrho();
00376   doublereal dens = delta * Rho_c;
00377   return (1.0 / (dens * dpdrho_val));
00378 }
00379 
00380 // Returns the value of dp / drho at constant T at  the current
00381 // state of the object
00382 /*
00383  *  units - Joules / kg
00384  *
00385  * @return  returns dpdrho
00386  */
00387 doublereal WaterPropsIAPWS::dpdrho() const {
00388   doublereal retn = m_phi->dimdpdrho(tau, delta);
00389   doublereal temperature = T_c/tau;
00390   doublereal val = retn * Rgas * temperature / M_water;
00391   return val;
00392 }
00393 
00394 // Returns the isochoric pressure derivative wrt temperature
00395 /*
00396  *     beta = M / (rho * Rgas) (d (pressure) / dT) at constant rho
00397  *
00398  *  Note for ideal gases this is equal to one.
00399  *
00400  *    beta = delta (phi0_d() + phiR_d())
00401  *            - tau delta (phi0_dt() + phiR_dt())
00402  */
00403 doublereal WaterPropsIAPWS:: coeffPresExp() const {
00404   doublereal retn = m_phi->dimdpdT(tau, delta);
00405   return (retn);
00406 }
00407 
00408 // Returns the coefficient of thermal expansion.
00409 /*
00410  *           alpha = d (ln V) / dT at constant P.
00411  *
00412  * @return  Returns the coefficient of thermal expansion
00413  */
00414 doublereal WaterPropsIAPWS:: coeffThermExp() const {
00415   doublereal kappa = isothermalCompressibility();
00416   doublereal beta = coeffPresExp();
00417   doublereal dens = delta * Rho_c;
00418   return (kappa * dens * Rgas * beta / M_water);
00419 }
00420 
00421 // Calculate the Gibbs free energy in mks units of J kmol-1 K-1.
00422 // using the last temperature and density
00423 doublereal WaterPropsIAPWS::Gibbs() const {
00424   doublereal gRT = m_phi->gibbs_RT();
00425   doublereal temperature = T_c/tau;
00426   return (gRT * Rgas * temperature);
00427 }
00428 
00429 
00430 // Utility routine in the calculation of the saturation pressure
00431 /*
00432  *  Private routine
00433  *
00434  * Calculate the Gibbs free energy in mks units of
00435  * J kmol-1 K-1.
00436  *
00437  * @param temperature    temperature (kelvin)
00438  * @param pressure       pressure (Pascal)
00439  * @param densLiq        Output density of liquid
00440  * @param densGas        output Density of gas
00441  * @param delGRT         output delGRT
00442  */
00443 void  WaterPropsIAPWS::
00444 corr(doublereal temperature, doublereal pressure, doublereal &densLiq, 
00445      doublereal &densGas, doublereal &delGRT) {
00446     
00447   densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
00448   if (densLiq <= 0.0) {
00449     throw Cantera::CanteraError("WaterPropsIAPWS::corr",
00450                                 "Error occurred trying to find liquid density at (T,P) = "
00451                                 + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
00452   }
00453   setState_TR(temperature, densLiq);
00454   doublereal gibbsLiqRT =  m_phi->gibbs_RT();
00455 
00456   densGas = density(temperature, pressure, WATER_GAS, densGas);
00457   if (densGas <= 0.0) {
00458     throw Cantera::CanteraError("WaterPropsIAPWS::corr",
00459                                 "Error occurred trying to find gas density at (T,P) = "
00460                                 + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
00461   }
00462   setState_TR(temperature, densGas);
00463   doublereal gibbsGasRT = m_phi->gibbs_RT();
00464     
00465   delGRT = gibbsLiqRT - gibbsGasRT;
00466 }
00467 
00468 // Utility routine in the calculation of the saturation pressure
00469 /*
00470  *  Private routine
00471  *
00472  * @param temperature    temperature (kelvin)
00473  * @param pressure       pressure (Pascal)
00474  * @param densLiq        Output density of liquid
00475  * @param densGas        output Density of gas
00476  * @param pcorr          output corrected pressure
00477  */
00478 void WaterPropsIAPWS::
00479 corr1(doublereal temperature, doublereal pressure, doublereal &densLiq, 
00480       doublereal &densGas, doublereal &pcorr) {
00481     
00482   densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
00483   if (densLiq <= 0.0) {
00484     throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
00485                                 "Error occurred trying to find liquid density at (T,P) = "
00486                                 + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
00487   }
00488   setState_TR(temperature, densLiq);
00489   doublereal prL = m_phi->phiR();
00490 
00491   densGas = density(temperature, pressure, WATER_GAS, densGas);
00492   if (densGas <= 0.0) {
00493     throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
00494                                 "Error occurred trying to find gas density at (T,P) = "
00495                                 + Cantera::fp2str(temperature) + "  " + Cantera::fp2str(pressure));
00496   }
00497   setState_TR(temperature, densGas);
00498   doublereal prG = m_phi->phiR();
00499     
00500   doublereal rhs = (prL - prG) + log(densLiq/densGas);
00501   rhs /= (1.0/densGas - 1.0/densLiq);
00502 
00503   pcorr = rhs * Rgas * temperature / M_water;
00504 }
00505 
00506 
00507 // This function returns the saturation pressure given the
00508 // temperature as an input parameter, and sets the internal state to the saturated
00509 // conditions.
00510 /*
00511  *  Note this function will return the saturation pressure, given the temperature.
00512  *  It will then set the state of the system to the saturation condition. The input
00513  *  parameter waterState is used to either specify the liquid state or the
00514  *  gas state at the desired temperatue and saturated pressure.
00515  *
00516  *  If the input temperature, T, is above T_c, this routine will set the internal
00517  *  state to T and the pressure to P_c. Then, return P_c.
00518  *
00519  * @param temperature   input temperature (kelvin)
00520  * @param waterState    integer specifying the water state
00521  *
00522  * @return Returns the saturation pressure
00523  *                units = Pascal
00524  */
00525 doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState) {
00526   static int method = 1;
00527   doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
00528   doublereal dp, pcorr;
00529   if (temperature >= T_c) {
00530     densGas = density(temperature, P_c, WATER_SUPERCRIT);
00531     setState_TR(temperature, densGas);
00532     return P_c;
00533   }
00534   doublereal p = psat_est(temperature);
00535   bool conv = false;
00536   for (int i = 0; i < 30; i++) {
00537     if (method == 1) {
00538       corr(temperature, p, densLiq, densGas, delGRT);
00539       doublereal delV = M_water * (1.0/densLiq - 1.0/densGas);
00540       dp = - delGRT * Rgas * temperature / delV;
00541     } else {
00542       corr1(temperature, p, densLiq, densGas, pcorr);
00543       dp = pcorr - p;
00544     }
00545     p += dp;
00546       
00547     if ((method == 1) && delGRT < 1.0E-8) {
00548       conv = true;
00549       break;
00550     } else {
00551       if (fabs(dp/p) < 1.0E-9) {
00552         conv = true;
00553         break;
00554       }
00555     }
00556   }
00557   // Put the fluid in the desired end condition
00558   if (waterState == WATER_LIQUID) {
00559     setState_TR(temperature, densLiq);
00560   } else if (waterState == WATER_GAS) {
00561     setState_TR(temperature, densGas);
00562   } else {
00563     throw Cantera::CanteraError("WaterPropsIAPWS::psat",
00564                                 "unknown water state input: " + Cantera::int2str(waterState));
00565   }
00566   return p;
00567 }
00568 
00569 // Returns the Phase State flag for the current state of the object
00570 /*
00571  * @param checkState If true, this function does a complete check to see where
00572  *        in paramters space we are
00573  *
00574  *  There are three values:
00575  *     WATER_GAS   below the critical temperature but below the critical density
00576  *     WATER_LIQUID  below the critical temperature but above the critical density
00577  *     WATER_SUPERCRIT   above the critical temperature
00578  */
00579 int WaterPropsIAPWS::phaseState(bool checkState) const {
00580   if (checkState) {
00581     if (tau <= 1.0) {
00582       iState = WATER_SUPERCRIT;
00583     } else {
00584       doublereal T = T_c / tau;
00585       doublereal rho = delta * Rho_c;
00586       //doublereal psatTable = psat_est(T);
00587       doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3);
00588       doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15);
00589       int iStateGuess = WATER_LIQUID;
00590       if (rho < rhoMid) {
00591         iStateGuess = WATER_GAS;
00592       }
00593       doublereal kappa = isothermalCompressibility();
00594       if (kappa >= 0.0) {
00595         iState = iStateGuess;
00596       } else {
00597         // When we are here we are between the spinodal curves
00598         doublereal rhoDel = rho * 1.000001;
00599 
00600         //setState_TR(T, rhoDel);
00601         doublereal deltaSave = delta;
00602         doublereal deltaDel = rhoDel / Rho_c;
00603         delta = deltaDel;
00604         m_phi->tdpolycalc(tau, deltaDel);
00605 
00606         doublereal kappaDel = isothermalCompressibility();
00607         doublereal d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho);
00608         if  (d2rhodp2 > 0.0) {
00609           iState = WATER_UNSTABLELIQUID;
00610         } else {
00611           iState = WATER_UNSTABLEGAS;
00612         }
00613         //setState_TR(T, rho);
00614         delta = deltaSave;
00615 
00616         m_phi->tdpolycalc(tau, delta);
00617       }
00618     }
00619   }
00620   return iState;
00621 }
00622 
00623 // Return the value of the density at the water spinodal point (on the liquid side)
00624 // for the current temperature.
00625 /*
00626  * @return returns the density with units of kg m-3
00627  */
00628 doublereal WaterPropsIAPWS::densSpinodalWater() const {
00629   doublereal temperature = T_c/tau;
00630   doublereal delta_save = delta;
00631   // return the critical density if we are above or even just a little below
00632   // the critical temperature. We just don't want to worry about the critical
00633   // point at this juncture.
00634   if (temperature >= T_c - 0.001) {
00635     return Rho_c;
00636   }
00637   doublereal p = psat_est(temperature);
00638   doublereal rho_low = 0.0;
00639   doublereal rho_high = 1000;
00640   
00641   doublereal densSatLiq = density_const(p, WATER_LIQUID);
00642   doublereal dens_old = densSatLiq;
00643   delta = dens_old / Rho_c;
00644   m_phi->tdpolycalc(tau, delta);
00645   doublereal dpdrho_old = dpdrho();
00646   if (dpdrho_old > 0.0) {
00647     rho_high = MIN(dens_old, rho_high);
00648   } else {
00649     rho_low = MAX(rho_low, dens_old);
00650   }
00651   doublereal dens_new = densSatLiq* (1.0001);
00652   delta = dens_new / Rho_c;
00653   m_phi->tdpolycalc(tau, delta);
00654   doublereal dpdrho_new = dpdrho();
00655   if (dpdrho_new > 0.0) {
00656     rho_high = MIN(dens_new, rho_high);
00657   } else {
00658     rho_low = MAX(rho_low, dens_new);
00659   }
00660   bool conv = false;
00661  
00662   for (int it = 0; it < 50; it++) {
00663     doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
00664     if (slope >= 0.0) {
00665       slope = MAX(slope, dpdrho_new *5.0/ dens_new);
00666     } else {  
00667       slope = -dpdrho_new;
00668       //slope = MIN(slope, dpdrho_new *5.0 / dens_new);
00669       // shouldn't be here for liquid spinodal
00670     }
00671     doublereal delta_rho =  - dpdrho_new / slope;
00672     if (delta_rho > 0.0) {
00673       delta_rho = MIN(delta_rho, dens_new * 0.1);
00674     } else {
00675       delta_rho = MAX(delta_rho, - dens_new * 0.1);
00676     }
00677     doublereal dens_est =  dens_new + delta_rho;
00678     if (dens_est < rho_low) {
00679       dens_est = 0.5 * (rho_low + dens_new); 
00680     }
00681     if (dens_est > rho_high) {
00682       dens_est = 0.5 * (rho_high + dens_new); 
00683     }
00684 
00685     
00686     dens_old = dens_new;
00687     dpdrho_old = dpdrho_new;
00688     dens_new = dens_est;
00689 
00690     delta = dens_new / Rho_c;
00691     m_phi->tdpolycalc(tau, delta);
00692     dpdrho_new = dpdrho();
00693     if (dpdrho_new > 0.0) {
00694       rho_high = MIN(dens_new, rho_high);
00695     } else  if (dpdrho_new < 0.0) {
00696       rho_low = MAX(rho_low, dens_new);
00697     } else {
00698       conv = true;
00699       break;
00700     }
00701     
00702     if (fabs(dpdrho_new) < 1.0E-5) {
00703       conv = true;
00704       break;
00705     }
00706   }
00707 
00708   if (!conv) {
00709     throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalWater()",
00710                                 " convergence failure");
00711   }
00712   // Restore the original delta
00713   delta = delta_save;
00714   m_phi->tdpolycalc(tau, delta);
00715 
00716   return dens_new;
00717 }
00718 
00719 // Return the value of the density at the water spinodal point (on the gas side)
00720 // for the current temperature.
00721 /*
00722  * @return returns the density with units of kg m-3
00723  */
00724 doublereal WaterPropsIAPWS::densSpinodalSteam() const {
00725   doublereal temperature = T_c/tau;
00726   doublereal delta_save = delta;
00727   // return the critical density if we are above or even just a little below
00728   // the critical temperature. We just don't want to worry about the critical
00729   // point at this juncture.
00730   if (temperature >= T_c - 0.001) {
00731     return Rho_c;
00732   }
00733   doublereal p = psat_est(temperature);
00734   doublereal rho_low = 0.0;
00735   doublereal rho_high = 1000;
00736   
00737   doublereal densSatGas = density_const(p, WATER_GAS);
00738   doublereal dens_old = densSatGas;
00739   delta = dens_old / Rho_c;
00740   m_phi->tdpolycalc(tau, delta);
00741   doublereal dpdrho_old = dpdrho();
00742   if (dpdrho_old < 0.0) {
00743     rho_high = MIN(dens_old, rho_high);
00744   } else {
00745     rho_low = MAX(rho_low, dens_old);
00746   }
00747   doublereal dens_new = densSatGas * (0.99);
00748   delta = dens_new / Rho_c;
00749   m_phi->tdpolycalc(tau, delta);
00750   doublereal dpdrho_new = dpdrho();
00751   if (dpdrho_new < 0.0) {
00752     rho_high = MIN(dens_new, rho_high);
00753   } else {
00754     rho_low = MAX(rho_low, dens_new);
00755   }
00756   bool conv = false;
00757  
00758   for (int it = 0; it < 50; it++) {
00759     doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
00760     if (slope >= 0.0) {
00761       slope = dpdrho_new;
00762       //slope = MAX(slope, dpdrho_new *5.0/ dens_new);
00763       // shouldn't be here for gas spinodal
00764     } else {  
00765       //slope = -dpdrho_new;
00766       slope = MIN(slope, dpdrho_new *5.0 / dens_new);
00767     
00768     }
00769     doublereal delta_rho = - dpdrho_new / slope;
00770     if (delta_rho > 0.0) {
00771       delta_rho = MIN(delta_rho, dens_new * 0.1);
00772     } else {
00773       delta_rho = MAX(delta_rho, - dens_new * 0.1);
00774     }
00775     doublereal dens_est =  dens_new + delta_rho;
00776     if (dens_est < rho_low) {
00777       dens_est = 0.5 * (rho_low + dens_new); 
00778     }
00779     if (dens_est > rho_high) {
00780       dens_est = 0.5 * (rho_high + dens_new); 
00781     }
00782 
00783     
00784     dens_old = dens_new;
00785     dpdrho_old = dpdrho_new;
00786     dens_new = dens_est;
00787 
00788     delta = dens_new / Rho_c;
00789     m_phi->tdpolycalc(tau, delta);
00790     dpdrho_new = dpdrho();
00791     if (dpdrho_new < 0.0) {
00792       rho_high = MIN(dens_new, rho_high);
00793     } else  if (dpdrho_new > 0.0) {
00794       rho_low = MAX(rho_low, dens_new);
00795     } else {
00796       conv = true;
00797       break;
00798     }
00799     
00800     if (fabs(dpdrho_new) < 1.0E-5) {
00801       conv = true;
00802       break;
00803     }
00804   }
00805 
00806   if (!conv) {
00807     throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalSteam()",
00808                                 " convergence failure");
00809   }
00810   // Restore the original delta
00811   delta = delta_save;
00812   m_phi->tdpolycalc(tau, delta);
00813 
00814   return dens_new;
00815 }
00816 
00817 /*
00818  * Sets the internal state of the object to the
00819  * specified temperature and density.
00820  */
00821 void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho) {
00822   calcDim(temperature, rho);
00823   m_phi->tdpolycalc(tau, delta);
00824 }
00825 
00826 /*
00827  * Calculate the enthalpy in mks units of
00828  * J kmol-1 K-1.
00829  */
00830 doublereal WaterPropsIAPWS::enthalpy() const {
00831   doublereal temperature = T_c/tau;
00832   doublereal hRT =  m_phi->enthalpy_RT();
00833   return (hRT * Rgas * temperature);
00834 }
00835 
00836 /*
00837  * Calculate the internal Energy in mks units of
00838  * J kmol-1 K-1.
00839  */
00840 doublereal WaterPropsIAPWS::intEnergy() const {
00841   doublereal temperature = T_c / tau;
00842   doublereal uRT = m_phi->intEnergy_RT();
00843   return (uRT * Rgas * temperature);
00844 }
00845 
00846 /*
00847  * Calculate the enthalpy in mks units of356
00848  * J kmol-1 K-1.
00849  */
00850 doublereal WaterPropsIAPWS::entropy() const {
00851   doublereal sR = m_phi->entropy_R();
00852   return (sR * Rgas);
00853 }
00854 
00855 /*
00856  * Calculate heat capacity at constant volume
00857  * J kmol-1 K-1.
00858  */
00859 doublereal WaterPropsIAPWS::cv() const {
00860   doublereal cvR = m_phi->cv_R();
00861   return (cvR * Rgas);
00862 }
00863 
00864 // Calculate the constant pressure heat capacity in mks units of J kmol-1 K-1
00865 // at the last temperature and density
00866 doublereal  WaterPropsIAPWS::cp() const {
00867   doublereal cpR = m_phi->cp_R();
00868   return (cpR * Rgas);
00869 }
00870 
00871 // Calculate the molar volume (kmol m-3)
00872 // at the last temperature and density
00873 doublereal WaterPropsIAPWS::molarVolume() const {
00874   doublereal rho = delta * Rho_c;
00875   return (M_water / rho);
00876 }
00877 
00878 }
Generated by  doxygen 1.6.3