PDSS_Water.cpp

Go to the documentation of this file.
00001 /**
00002  * @file PDSS_Water.cpp
00003  *
00004  */
00005 /*
00006  * Copywrite (2006) Sandia Corporation. Under the terms of
00007  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008  * U.S. Government retains certain rights in this software.
00009  */
00010 /*
00011  * $Id: PDSS_Water.cpp 385 2010-01-17 17:05:46Z hkmoffa $
00012  */
00013 #include "ct_defs.h"
00014 
00015 #include "xml.h"
00016 #include "ctml.h"
00017 #include "PDSS_Water.h"
00018 
00019 #include "WaterPropsIAPWS.h"
00020 #include "ThermoFactory.h"
00021 #include "WaterProps.h"
00022 #include "VPStandardStateTP.h"
00023 
00024 #include <cmath>
00025 
00026 namespace Cantera {
00027   /**
00028    * Basic list of constructors and duplicators
00029    */
00030   PDSS_Water::PDSS_Water() :
00031     PDSS(),
00032     m_sub(0),
00033     m_waterProps(0),
00034     m_dens(1000.0),
00035     m_iState(WATER_LIQUID),
00036     EW_Offset(0.0),
00037     SW_Offset(0.0),
00038     m_verbose(0),
00039     m_allowGasPhase(false)
00040   {
00041     m_pdssType = cPDSS_WATER;
00042     m_sub = new WaterPropsIAPWS();
00043     m_waterProps =  new WaterProps(m_sub);
00044     m_spthermo = 0;
00045     constructSet();
00046     m_minTemp = 200.;
00047     m_maxTemp = 10000.;
00048   }
00049 
00050   PDSS_Water::PDSS_Water(VPStandardStateTP *tp, int spindex) :
00051     PDSS(tp, spindex),
00052     m_sub(0), 
00053     m_waterProps(0),
00054     m_dens(1000.0),
00055     m_iState(WATER_LIQUID),
00056     EW_Offset(0.0),
00057     SW_Offset(0.0),
00058     m_verbose(0),
00059     m_allowGasPhase(false)
00060   {
00061     m_pdssType = cPDSS_WATER;
00062     m_sub = new WaterPropsIAPWS();
00063     m_waterProps =  new WaterProps(m_sub);
00064     m_spthermo = 0;
00065     constructSet();
00066     m_minTemp = 200.;
00067     m_maxTemp = 10000.;
00068   }
00069 
00070 
00071   PDSS_Water::PDSS_Water(VPStandardStateTP *tp, int spindex, 
00072                          std::string inputFile, std::string id) :
00073     PDSS(tp, spindex),
00074     m_sub(0),
00075     m_waterProps(0),
00076     m_dens(1000.0),
00077     m_iState(WATER_LIQUID),
00078     EW_Offset(0.0),
00079     SW_Offset(0.0),
00080     m_verbose(0),
00081     m_allowGasPhase(false)
00082   {
00083     m_pdssType = cPDSS_WATER;
00084     m_sub = new WaterPropsIAPWS();
00085     m_waterProps =  new WaterProps(m_sub); 
00086     constructPDSSFile(tp, spindex, inputFile, id);
00087     m_spthermo = 0;
00088     m_minTemp = 200.;
00089     m_maxTemp = 10000.;
00090   }
00091 
00092   PDSS_Water::PDSS_Water(VPStandardStateTP *tp, int spindex,
00093                          const XML_Node& speciesNode, 
00094                          const XML_Node& phaseRoot, bool spInstalled) :
00095     PDSS(tp, spindex),
00096     m_sub(0),
00097     m_waterProps(0),
00098     m_dens(1000.0),
00099     m_iState(WATER_LIQUID),
00100     EW_Offset(0.0),
00101     SW_Offset(0.0),
00102     m_verbose(0),
00103     m_allowGasPhase(false)
00104   {
00105     m_pdssType = cPDSS_WATER;
00106     m_sub = new WaterPropsIAPWS();
00107     m_waterProps =  new WaterProps(m_sub);
00108     std::string id= "";
00109     constructPDSSXML(tp, spindex, phaseRoot, id) ;
00110     initThermo();
00111     m_spthermo = 0;
00112     m_minTemp = 200.;
00113     m_maxTemp = 10000.;
00114   }
00115 
00116 
00117 
00118   PDSS_Water::PDSS_Water(const PDSS_Water &b) :
00119     PDSS(),
00120     m_sub(0),
00121     m_waterProps(0),
00122     m_dens(1000.0),
00123     m_iState(WATER_LIQUID),
00124     EW_Offset(b.EW_Offset),
00125     SW_Offset(b.SW_Offset),
00126     m_verbose(b.m_verbose),
00127     m_allowGasPhase(b.m_allowGasPhase)
00128   {
00129     m_sub = new WaterPropsIAPWS();
00130     /*
00131      * Use the assignment operator to do the brunt
00132      * of the work for the copy construtor.
00133      */
00134     *this = b;
00135   }
00136 
00137   /**
00138    * Assignment operator
00139    */
00140   PDSS_Water& PDSS_Water::operator=(const PDSS_Water&b) {
00141     if (&b == this) return *this;
00142     /*
00143      * Call the base class operator
00144      */
00145     PDSS::operator=(b);
00146 
00147     if (!m_sub) {
00148       m_sub = new WaterPropsIAPWS();
00149     }
00150     m_sub->operator=(*(b.m_sub));
00151 
00152     if (!m_waterProps) {
00153       m_waterProps = new WaterProps(m_sub);
00154     }
00155     m_waterProps->operator=(*(b.m_waterProps));
00156 
00157     m_dens          = b.m_dens;
00158     m_iState        = b.m_iState;
00159     EW_Offset       = b.EW_Offset;
00160     SW_Offset       = b.SW_Offset;
00161     m_verbose       = b.m_verbose;
00162     m_allowGasPhase = b.m_allowGasPhase;
00163 
00164     return *this;
00165   }
00166 
00167   PDSS_Water::~PDSS_Water() {
00168     delete m_waterProps;
00169     delete m_sub;
00170   }
00171 
00172   PDSS *PDSS_Water::duplMyselfAsPDSS() const {
00173     PDSS_Water *kPDSS = new PDSS_Water(*this);
00174     return (PDSS *) kPDSS;
00175   }
00176 
00177   /*
00178    * constructPDSSXML:
00179    *
00180    * Initialization of a Debye-Huckel phase using an
00181    * xml file.
00182    *
00183    * This routine is a precursor to  constructSet
00184    * routine, which does most of the work.
00185    *
00186    * @param infile XML file containing the description of the
00187    *        phase
00188    *
00189    * @param id  Optional parameter identifying the name of the
00190    *            phase. If none is given, the first XML
00191    *            phase element will be used.
00192    */
00193   void PDSS_Water::constructPDSSXML(VPStandardStateTP *tp, int spindex,
00194                                     const XML_Node& phaseNode, std::string id) {
00195     constructSet();
00196   }
00197    
00198   /*
00199    * constructPDSSFile():
00200    *
00201    * Initialization of a Debye-Huckel phase using an
00202    * xml file.
00203    *
00204    * This routine is a precursor to constructPDSSXML(XML_Node*)
00205    * routine, which does most of the work.
00206    *
00207    * @param infile XML file containing the description of the
00208    *        phase
00209    *
00210    * @param id  Optional parameter identifying the name of the
00211    *            phase. If none is given, the first XML
00212    *            phase element will be used.
00213    */
00214   void PDSS_Water::constructPDSSFile(VPStandardStateTP *tp, int spindex,
00215                                      std::string inputFile, std::string id) {
00216 
00217     if (inputFile.size() == 0) {
00218       throw CanteraError("PDSS_Water::constructPDSSFile",
00219                          "input file is null");
00220     }
00221     std::string path = findInputFile(inputFile);
00222     std::ifstream fin(path.c_str());
00223     if (!fin) {
00224       throw CanteraError("PDSS_Water::initThermo","could not open "
00225                          +path+" for reading.");
00226     }
00227     /*
00228      * The phase object automatically constructs an XML object.
00229      * Use this object to store information.
00230      */
00231 
00232     XML_Node *fxml = new XML_Node();
00233     fxml->build(fin);
00234     XML_Node *fxml_phase = findXMLPhase(fxml, id);
00235     if (!fxml_phase) {
00236       throw CanteraError("PDSS_Water::initThermo",
00237                          "ERROR: Can not find phase named " +
00238                          id + " in file named " + inputFile);
00239     }   
00240     constructPDSSXML(tp, spindex, *fxml_phase, id);
00241     delete fxml;
00242   }
00243 
00244 
00245 
00246   void PDSS_Water::constructSet() {
00247     if (m_sub) delete m_sub;
00248     m_sub = new WaterPropsIAPWS();
00249     if (m_sub == 0) {
00250       throw CanteraError("PDSS_Water::initThermo",
00251                          "could not create new substance object.");
00252     }
00253     /*
00254      * Calculate the molecular weight. 
00255      *  hard coded to Cantera's elements and Water.
00256      */
00257     m_mw = 2 * 1.00794 + 15.9994;
00258 
00259     /*
00260      * Set the baseline 
00261      */
00262     doublereal T = 298.15;
00263 
00264     m_p0 = OneAtm;
00265 
00266     doublereal presLow = 1.0E-2;
00267     doublereal oneBar = 1.0E5;
00268     doublereal dens = 1.0E-9;
00269     m_dens = m_sub->density(T, presLow, WATER_GAS, dens);
00270     m_pres = presLow;
00271     SW_Offset = 0.0;
00272     doublereal s = entropy_mole();
00273     s -=  GasConstant * log(oneBar/presLow);
00274     if (s != 188.835E3) {
00275       SW_Offset = 188.835E3 - s;
00276     }
00277     s = entropy_mole();
00278     s -=  GasConstant * log(oneBar/presLow);
00279     //printf("s = %g\n", s);
00280 
00281     doublereal h = enthalpy_mole();
00282     if (h != -241.826E6) {
00283       EW_Offset = -241.826E6 - h;
00284     }
00285     h = enthalpy_mole();
00286     //printf("h = %g\n", h);
00287 
00288     /*
00289      * Set the initial state of the system to 298.15 K and 
00290      * 1 bar.
00291      */
00292     setTemperature(298.15);
00293     m_dens = m_sub->density(298.15, OneAtm, WATER_LIQUID);
00294     m_pres = OneAtm;
00295   }
00296 
00297   void PDSS_Water::initThermo() {
00298     PDSS::initThermo();
00299   }
00300 
00301   void PDSS_Water::initThermoXML(const XML_Node& phaseNode, std::string& id) {
00302     PDSS::initThermoXML(phaseNode, id);
00303   }
00304 
00305   doublereal PDSS_Water::enthalpy_mole() const {
00306     doublereal h = m_sub->enthalpy();
00307     return (h + EW_Offset);
00308   }
00309 
00310   doublereal PDSS_Water::intEnergy_mole() const {
00311     doublereal u = m_sub->intEnergy();
00312     return (u + EW_Offset);            
00313   }
00314 
00315   doublereal PDSS_Water::entropy_mole() const {
00316     doublereal s = m_sub->entropy();
00317     return (s + SW_Offset); 
00318   }
00319 
00320   doublereal PDSS_Water::gibbs_mole() const {
00321     doublereal g = m_sub->Gibbs();
00322     return (g + EW_Offset - SW_Offset*m_temp);
00323   }
00324 
00325   doublereal PDSS_Water::cp_mole() const {
00326     doublereal cp = m_sub->cp();
00327     return cp; 
00328   }
00329 
00330   doublereal PDSS_Water::cv_mole() const {
00331     doublereal cv = m_sub->cv();
00332     return cv;
00333   }
00334 
00335   doublereal  PDSS_Water::molarVolume() const {
00336     doublereal mv = m_sub->molarVolume();
00337     return (mv);
00338   }
00339 
00340   doublereal PDSS_Water::gibbs_RT_ref() const {
00341     doublereal T = m_temp;
00342     m_sub->density(T, m_p0);
00343     doublereal h = m_sub->enthalpy();
00344     m_sub->setState_TR(m_temp, m_dens);
00345     return ((h + EW_Offset - SW_Offset*T)/(T * GasConstant));
00346   }
00347 
00348   doublereal PDSS_Water::enthalpy_RT_ref() const {
00349     doublereal T = m_temp;
00350     m_sub->density(T, m_p0);
00351     doublereal h = m_sub->enthalpy();
00352     m_sub->setState_TR(m_temp, m_dens);
00353     return ((h + EW_Offset)/(T * GasConstant));
00354   }
00355 
00356   doublereal PDSS_Water::entropy_R_ref() const {
00357     doublereal T = m_temp;
00358     m_sub->density(T, m_p0);
00359     doublereal s = m_sub->entropy();
00360     m_sub->setState_TR(m_temp, m_dens);
00361     return ((s + SW_Offset)/GasConstant); 
00362   }
00363 
00364   doublereal PDSS_Water::cp_R_ref() const {
00365     doublereal T = m_temp;
00366     m_sub->density(T, m_p0);
00367     doublereal cp = m_sub->cp();
00368     m_sub->setState_TR(m_temp, m_dens);
00369     return (cp/GasConstant); 
00370   }
00371 
00372   doublereal PDSS_Water::molarVolume_ref() const {
00373     doublereal T = m_temp;
00374     m_sub->density(T, m_p0);
00375     doublereal mv = m_sub->molarVolume();
00376     m_sub->setState_TR(m_temp, m_dens);
00377     return (mv); 
00378   }
00379 
00380 
00381   /**
00382    * Calculate the pressure (Pascals), given the temperature and density
00383    *  Temperature: kelvin
00384    *  rho: density in kg m-3
00385    */
00386   doublereal PDSS_Water::pressure() const {
00387     doublereal p = m_sub->pressure();
00388     m_pres = p;
00389     return p;
00390   }
00391 
00392 
00393   // In this routine we must be sure to only find the water branch of the
00394   // curve and not the gas branch
00395   void PDSS_Water::setPressure(doublereal p) {
00396     doublereal T = m_temp;
00397     doublereal dens = m_dens;
00398     int waterState = WATER_LIQUID;
00399     if (T > m_sub->Tcrit()) {
00400       waterState = WATER_SUPERCRIT;
00401     }
00402  
00403 
00404 #ifdef DEBUG_HKM
00405     //printf("waterPDSS: set pres = %g t = %g, waterState = %d\n",
00406     //      p, T, waterState);
00407 #endif
00408     doublereal dd = m_sub->density(T, p, waterState, dens);
00409     if (dd <= 0.0) {
00410       std::string stateString = "T = " +
00411         fp2str(T) + " K and p = " + fp2str(p) + " Pa";
00412       throw CanteraError("PDSS_Water:setPressure()", 
00413                          "Failed to set water SS state: " + stateString);
00414     }
00415     m_dens = dd;
00416     m_pres = p;
00417 
00418     // We are only putting the phase check here because of speed considerations.
00419     m_iState = m_sub->phaseState(true);
00420     if (! m_allowGasPhase) {
00421       if (m_iState != WATER_SUPERCRIT && m_iState != WATER_LIQUID && m_iState != WATER_UNSTABLELIQUID) {
00422         throw CanteraError("PDSS_Water::setPressure",
00423                            "Water State isn't liquid or crit");
00424       }
00425     }
00426   }
00427  
00428   // Return the volumetric thermal expansion coefficient. Units: 1/K.
00429   /*
00430    * The thermal expansion coefficient is defined as
00431    * \f[
00432    * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
00433    * \f]
00434    */
00435   doublereal PDSS_Water::thermalExpansionCoeff() const {
00436     doublereal val = m_sub->coeffThermExp();
00437     return val; 
00438   }
00439 
00440   doublereal PDSS_Water::dthermalExpansionCoeffdT() const {
00441     doublereal pres = pressure();
00442     doublereal dens_save = m_dens;
00443     doublereal tt = m_temp - 0.04;
00444     doublereal dd = m_sub->density(tt, pres, m_iState, m_dens);
00445     if (dd < 0.0) {
00446       throw CanteraError("PDSS_Water::dthermalExpansionCoeffdT", 
00447                          "unable to solve for the density at T = " + fp2str(tt) + ", P = " + fp2str(pres));
00448     }
00449     doublereal vald = m_sub->coeffThermExp();
00450     m_sub->setState_TR(m_temp, dens_save);
00451     doublereal val2 = m_sub->coeffThermExp();
00452     doublereal val = (val2 - vald) / 0.04;
00453     return val; 
00454   }
00455    
00456   doublereal PDSS_Water::isothermalCompressibility() const {
00457     doublereal val = m_sub->isothermalCompressibility();
00458     return val; 
00459   }
00460 
00461   /// critical temperature 
00462   doublereal PDSS_Water::critTemperature() const { return m_sub->Tcrit(); }
00463         
00464   /// critical pressure
00465   doublereal PDSS_Water::critPressure() const { return m_sub->Pcrit(); }
00466         
00467   /// critical density
00468   doublereal PDSS_Water::critDensity() const { return m_sub->Rhocrit(); }
00469         
00470   void PDSS_Water::setDensity(doublereal dens) {
00471     m_dens = dens;
00472     m_sub->setState_TR(m_temp, m_dens);
00473   }
00474 
00475   doublereal PDSS_Water::density() const {
00476     return m_dens;
00477   }
00478 
00479   void PDSS_Water::setTemperature(doublereal temp) {
00480     m_temp = temp;
00481     doublereal dd = m_dens;
00482     m_sub->setState_TR(temp, dd);
00483   }
00484 
00485   void PDSS_Water::setState_TP(doublereal temp, doublereal pres) {
00486     m_temp = temp;
00487     setPressure(pres);
00488   }
00489 
00490   void PDSS_Water::setState_TR(doublereal temp, doublereal dens) {
00491     m_temp = temp;
00492     m_dens = dens;
00493     m_sub->setState_TR(m_temp, m_dens);
00494   }
00495 
00496   doublereal PDSS_Water::pref_safe(doublereal temp) const {
00497     if (temp < m_sub->Tcrit()) {
00498       doublereal pp = m_sub->psat_est(temp);
00499       if (pp > OneAtm) {
00500         return pp;
00501       }
00502     } else  {
00503       return m_sub->Pcrit();
00504     }
00505     return OneAtm;
00506   }
00507 
00508   // saturation pressure
00509   doublereal PDSS_Water::satPressure(doublereal t){
00510     doublereal pp = m_sub->psat(t, WATER_LIQUID);
00511     m_dens = m_sub->density();
00512     m_temp = t;
00513     return pp;
00514   }
00515   
00516 }
Generated by  doxygen 1.6.3