Cantera  2.0
IdealMolalSoln.h
Go to the documentation of this file.
1 /**
2  * @file IdealMolalSoln.h
3  * ThermoPhase object for the ideal molal equation of
4  * state (see \ref thermoprops
5  * and class \link Cantera::IdealMolalSoln IdealMolalSoln\endlink).
6  *
7  * Header file for a derived class of ThermoPhase that handles
8  * variable pressure standard state methods for calculating
9  * thermodynamic properties that are further based upon
10  * activities on the molality scale. The Ideal molal
11  * solution assumes that all molality-based activity
12  * coefficients are equal to one. This turns out to be highly
13  * nonlinear in the limit of the solvent mole fraction going
14  * to zero.
15  */
16 /*
17  * Copyright (2006) Sandia Corporation. Under the terms of
18  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
19  * U.S. Government retains certain rights in this software.
20  */
21 #ifndef CT_IDEALMOLALSOLN_H
22 #define CT_IDEALMOLALSOLN_H
23 
24 #include "MolalityVPSSTP.h"
25 
26 namespace Cantera
27 {
28 
29 /** \addtogroup thermoprops */
30 /* @{
31  */
32 
33 
34 /**
35  * This phase is based upon the mixing-rule assumption that
36  * all molality-based activity coefficients are equal
37  * to one.
38  *
39  * This is a full instantiation of a ThermoPhase object.
40  * The assumption is that the molality-based activity
41  * coefficient is equal to one. This also implies that
42  * the osmotic coefficient is equal to one.
43  *
44  * Note, this does not mean that the solution is an
45  * ideal solution. In fact, there is a singularity in
46  * the formulation as
47  * the solvent concentration goes to zero.
48  *
49  * The mechanical equation of state is currently assumed to
50  * be that of an incompressible solution. This may change
51  * in the future. Each species has its own molar volume.
52  * The molar volume is a constant.
53  *
54  * Class IdealMolalSoln represents a condensed phase.
55  * The phase and the pure species phases which
56  * comprise the standard states of the species are assumed to have
57  * zero volume expansivity and zero isothermal compressibility.
58  * Each species does, however, have constant but distinct partial
59  * molar volumes equal to their pure species molar volumes.
60  * The class derives from class ThermoPhase,
61  * and overloads the virtual methods defined there with ones that
62  * use expressions appropriate for incompressible mixtures.
63  *
64  * The standard concentrations can have three different forms
65  * depending on the value of the member attribute m_formGC, which
66  * is supplied in the XML file.
67  *
68  * <TABLE>
69  * <TR><TD> m_formGC </TD><TD> ActivityConc </TD><TD> StandardConc </TD></TR>
70  * <TR><TD> 0 </TD><TD> \f$ {m_k}/ { m^{\Delta}}\f$ </TD><TD> \f$ 1.0 \f$ </TD></TR>
71  * <TR><TD> 1 </TD><TD> \f$ m_k / (m^{\Delta} V_k)\f$ </TD><TD> \f$ 1.0 / V_k \f$ </TD></TR>
72  * <TR><TD> 2 </TD><TD> \f$ m_k / (m^{\Delta} V^0_0)\f$</TD><TD> \f$ 1.0 / V^0_0\f$ </TD></TR>
73  * </TABLE>
74  *
75  * \f$ V^0_0 \f$ is the solvent standard molar volume. \f$ m^{\Delta} \f$ is a constant equal to a
76  * molality of \f$ 1.0 \quad\mbox{gm kmol}^{-1} \f$.
77  *
78  * The current default is to have mformGC = 2.
79  *
80  * The value and form of the activity concentration will affect
81  * reaction rate constants involving species in this phase.
82  *
83  * @verbatim
84  <thermo model="IdealMolalSoln">
85  <standardConc model="solvent_volume" />
86  <solvent> H2O(l) </solvent>
87 
88  <activityCoefficients model="IdealMolalSoln" >
89  <idealMolalSolnCutoff model="polyExp">
90  <gamma_O_limit> 1.0E-5 <gammaOlimit>
91  <gamma_k_limit> 1.0E-5 <gammaklimit>
92  <X_o_cutoff> 0.20 </X_o_cutoff>
93  <C_0_param> 0.05 </C_0_param>
94  <slope_f_limit> 0.6 </slopefLimit>
95  <slope_g_limit> 0.0 </slopegLimit>
96  </idealMolalSolnCutoff>
97  </activityCoefficients>
98 
99 
100 
101  </thermo>
102 
103 
104 
105  @endverbatim
106  *
107  */
109 {
110 
111 public:
112 
113  /// Constructors
114  IdealMolalSoln();
115 
116  //! Copy Constructor
118 
119  //! Assignment operator
121 
122  //! Constructor for phase initialization
123  /*!
124  * This constructor will initialize a phase, by reading the required
125  * information from an input file.
126  *
127  * @param inputFile Name of the Input file that contains information about the phase
128  * @param id id of the phase within the input file
129  */
130  IdealMolalSoln(std::string inputFile, std::string id = "");
131 
132  //! Constructor for phase initialization
133  /*!
134  * This constructor will initialize a phase, by reading the required
135  * information from XML_Node tree.
136  *
137  * @param phaseRef reference for an XML_Node tree that contains
138  * the information necessary to initialize the phase.
139  * @param id id of the phase within the input file
140  */
141  IdealMolalSoln(XML_Node& phaseRef, std::string id = "");
142 
143  /// Destructor.
144  virtual ~IdealMolalSoln();
145 
146  //! Duplication function
147  /*!
148  * This virtual function is used to create a duplicate of the
149  * current phase. It's used to duplicate the phase when given
150  * a ThermoPhase pointer to the phase.
151  *
152  * @return It returns a ThermoPhase pointer.
153  */
155 
156  /**
157  *
158  * @name Utilities
159  * @{
160  */
161 
162  /**
163  * Equation of state type flag. The base class returns
164  * zero. Subclasses should define this to return a unique
165  * non-zero value. Constants defined for this purpose are
166  * listed in mix_defs.h.
167  */
168  virtual int eosType() const {
169  return 0;
170  }
171 
172  /**
173  * @}
174  * @name Molar Thermodynamic Properties of the Solution ---------------
175  * @{
176  */
177 
178  //! Molar enthalpy of the solution. Units: J/kmol.
179  /*!
180  *
181  * Returns the amount of enthalpy per mole of solution.
182  * For an ideal molal solution,
183  * \f[
184  * \bar{h}(T, P, X_k) = \sum_k X_k \bar{h}_k(T)
185  * \f]
186  * The formula is written in terms of the partial molar enthalpies.
187  * \f$ \bar{h}_k(T, p, m_k) \f$.
188  * See the partial molar enthalpy function, getPartialMolarEnthalpies(),
189  * for details.
190  *
191  * Units: J/kmol
192  */
193  virtual doublereal enthalpy_mole() const;
194 
195  //! Molar internal energy of the solution: Units: J/kmol.
196  /*!
197  *
198  * Returns the amount of internal energy per mole of solution.
199  * For an ideal molal solution,
200  * \f[
201  * \bar{u}(T, P, X_k) = \sum_k X_k \bar{u}_k(T)
202  * \f]
203  * The formula is written in terms of the partial molar internal energy.
204  * \f$ \bar{u}_k(T, p, m_k) \f$.
205  */
206  virtual doublereal intEnergy_mole() const;
207 
208  //! Molar entropy of the solution. Units: J/kmol/K.
209  /*!
210  * Returns the amount of entropy per mole of solution.
211  * For an ideal molal solution,
212  * \f[
213  * \bar{s}(T, P, X_k) = \sum_k X_k \bar{s}_k(T)
214  * \f]
215  * The formula is written in terms of the partial molar entropies.
216  * \f$ \bar{s}_k(T, p, m_k) \f$.
217  * See the partial molar entropies function, getPartialMolarEntropies(),
218  * for details.
219  *
220  * Units: J/kmol/K.
221  */
222  virtual doublereal entropy_mole() const;
223 
224  //! Molar Gibbs function for the solution: Units J/kmol.
225  /*!
226  *
227  * Returns the gibbs free energy of the solution per mole
228  * of the solution.
229  *
230  * \f[
231  * \bar{g}(T, P, X_k) = \sum_k X_k \mu_k(T)
232  * \f]
233  *
234  * Units: J/kmol
235  */
236  virtual doublereal gibbs_mole() const;
237 
238  //! Molar heat capacity of the solution at constant pressure. Units: J/kmol/K.
239  /*!
240  * \f[
241  * \bar{c}_p(T, P, X_k) = \sum_k X_k \bar{c}_{p,k}(T)
242  * \f]
243  *
244  * Units: J/kmol/K
245  */
246  virtual doublereal cp_mole() const;
247 
248  //! Molar heat capacity of the solution at constant volume. Units: J/kmol/K.
249  /*!
250  * Molar heat capacity at constant volume: Units: J/kmol/K.
251  * NOT IMPLEMENTED.
252  * Units: J/kmol/K
253  */
254  virtual doublereal cv_mole() const;
255 
256  //@}
257  /** @name Mechanical Equation of State Properties -------------------------
258  //@{
259  *
260  * In this equation of state implementation, the density is a
261  * function only of the mole fractions. Therefore, it can't be
262  * an independent variable. Instead, the pressure is used as the
263  * independent variable. Functions which try to set the thermodynamic
264  * state by calling setDensity() may cause an exception to be
265  * thrown.
266  */
267 
268 
269 
270  /**
271  * Set the pressure at constant temperature. Units: Pa.
272  * This method sets a constant within the object.
273  * The mass density is not a function of pressure.
274  *
275  * @param p Input Pressure
276  */
277  virtual void setPressure(doublereal p);
278 
279 protected:
280  /**
281  * Calculate the density of the mixture using the partial
282  * molar volumes and mole fractions as input
283  *
284  * The formula for this is
285  *
286  * \f[
287  * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
288  * \f]
289  *
290  * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
291  * the molecular weights, and \f$V_k\f$ are the pure species
292  * molar volumes.
293  *
294  * Note, the basis behind this formula is that in an ideal
295  * solution the partial molar volumes are equal to the pure
296  * species molar volumes. We have additionally specified
297  * in this class that the pure species molar volumes are
298  * independent of temperature and pressure.
299  *
300  * NOTE: This is a non-virtual function, which is not a
301  * member of the ThermoPhase base class.
302  */
303  void calcDensity();
304 
305 public:
306  /**
307  * Overwritten setDensity() function is necessary because the
308  * density is not an independent variable.
309  *
310  * This function will now throw an error condition
311  *
312  * @internal May have to adjust the strategy here to make
313  * the eos for these materials slightly compressible, in order
314  * to create a condition where the density is a function of
315  * the pressure.
316  *
317  * This function will now throw an error condition.
318  *
319  * NOTE: This is an overwritten function from the State.h
320  * class
321  *
322  * @param rho Input Density
323  */
324  void setDensity(const doublereal rho);
325 
326  /**
327  * Overwritten setMolarDensity() function is necessary because the
328  * density is not an independent variable.
329  *
330  * This function will now throw an error condition.
331  *
332  * NOTE: This is an overwritten function from the State.h
333  * class
334  *
335  * @param rho Input Density
336  */
337  void setMolarDensity(const doublereal rho);
338 
339  //! Set the temperature (K) and pressure (Pa)
340  /*!
341  * Set the temperature and pressure.
342  *
343  * @param t Temperature (K)
344  * @param p Pressure (Pa)
345  */
346  virtual void setState_TP(doublereal t, doublereal p);
347 
348  //! The isothermal compressibility. Units: 1/Pa.
349  /*!
350  * The isothermal compressibility is defined as
351  * \f[
352  * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
353  * \f]
354  *
355  * It's equal to zero for this model, since the molar volume
356  * doesn't change with pressure or temperature.
357  */
358  virtual doublereal isothermalCompressibility() const;
359 
360  //! The thermal expansion coefficient. Units: 1/K.
361  /*!
362  * The thermal expansion coefficient is defined as
363  *
364  * \f[
365  * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
366  * \f]
367  *
368  * It's equal to zero for this model, since the molar volume
369  * doesn't change with pressure or temperature.
370  */
371  virtual doublereal thermalExpansionCoeff() const;
372 
373  /**
374  * @}
375  * @name Potential Energy
376  *
377  * Species may have an additional potential energy due to the
378  * presence of external gravitation or electric fields. These
379  * methods allow specifying a potential energy for individual
380  * species.
381  * @{
382  */
383 
384 
385  //!Set the potential energy of species k to pe.
386  /*!
387  * Units: J/kmol.
388  * This function must be reimplemented in inherited classes
389  * of ThermoPhase.
390  *
391  * @param k Species index
392  * @param pe Input potential energy.
393  */
394  virtual void setPotentialEnergy(int k, doublereal pe) {
395  err("setPotentialEnergy");
396  }
397 
398  /*
399  * Get the potential energy of species k.
400  * Units: J/kmol.
401  * This function must be reimplemented in inherited classes
402  * of ThermoPhase.
403  *
404  * @param k Species index
405  */
406  virtual doublereal potentialEnergy(int k) const {
407  return err("potentialEnergy");
408  }
409 
410  /*
411  * Set the electric potential of this phase (V).
412  * This is used by classes InterfaceKinetics and EdgeKinetics to
413  * compute the rates of charge-transfer reactions, and in computing
414  * the electrochemical potentials of the species.
415  *
416  * @param v input Electric Potential (volts).
417  */
418  void setElectricPotential(doublereal v) {
419  m_phi = v;
420  }
421 
422  //! Returns the electric potential of this phase (V).
423  doublereal electricPotential() const {
424  return m_phi;
425  }
426 
427 
428  /**
429  * @}
430  * @name Activities and Activity Concentrations
431  *
432  * The activity \f$a_k\f$ of a species in solution is
433  * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
434  * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T)\f$ is
435  * the chemical potential at unit activity, which depends only
436  * on temperature and the pressure.
437  * @{
438  */
439 
440  /*!
441  * This method returns an array of generalized concentrations
442  * \f$ C_k\f$ that are defined such that
443  * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
444  * is a standard concentration
445  * defined below. These generalized concentrations are used
446  * by kinetics manager classes to compute the forward and
447  * reverse rates of elementary reactions.
448  *
449  * @param c Array of generalized concentrations. The
450  * units depend upon the implementation of the
451  * reaction rate expressions within the phase.
452  */
453  virtual void getActivityConcentrations(doublereal* c) const;
454 
455  /**
456  * The standard concentration \f$ C^0_k \f$ used to normalize
457  * the generalized concentration. In many cases, this quantity
458  * will be the same for all species in a phase - for example,
459  * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
460  * reason, this method returns a single value, instead of an
461  * array. However, for phases in which the standard
462  * concentration is species-specific (e.g. surface species of
463  * different sizes), this method may be called with an
464  * optional parameter indicating the species.
465  *
466  * @param k Species index
467  */
468  virtual doublereal standardConcentration(size_t k=0) const;
469 
470  /*!
471  * Returns the natural logarithm of the standard
472  * concentration of the kth species
473  *
474  * @param k Species index
475  */
476  virtual doublereal logStandardConc(size_t k=0) const;
477 
478  /*!
479  * Returns the units of the standard and generalized
480  * concentrations Note they have the same units, as their
481  * ratio is defined to be equal to the activity of the kth
482  * species in the solution, which is unitless.
483  *
484  * This routine is used in print out applications where the
485  * units are needed. Usually, MKS units are assumed throughout
486  * the program and in the XML input files.
487  *
488  * @param uA Output vector containing the units
489  * uA[0] = kmol units - default = 1
490  * uA[1] = m units - default = -nDim(), the number of spatial
491  * dimensions in the Phase class.
492  * uA[2] = kg units - default = 0;
493  * uA[3] = Pa(pressure) units - default = 0;
494  * uA[4] = Temperature units - default = 0;
495  * uA[5] = time units - default = 0
496  * @param k species index. Defaults to 0.
497  * @param sizeUA output int containing the size of the vector.
498  * Currently, this is equal to 6.
499  */
500  virtual void getUnitsStandardConc(double* uA, int k = 0,
501  int sizeUA = 6) const;
502 
503  /*!
504  * Get the array of non-dimensional activities at
505  * the current solution temperature, pressure, and
506  * solution concentration.
507  *
508  * (note solvent is on molar scale)
509  *
510  * @param ac Output activity coefficients.
511  * Length: m_kk.
512  */
513  virtual void getActivities(doublereal* ac) const;
514 
515  /*!
516  * Get the array of non-dimensional molality-based
517  * activity coefficients at the current solution temperature,
518  * pressure, and solution concentration.
519  *
520  *
521  * (note solvent is on molar scale. The solvent molar
522  * based activity coefficient is returned).
523  *
524  * @param acMolality Output Molality-based activity coefficients.
525  * Length: m_kk.
526  */
527  virtual void
528  getMolalityActivityCoefficients(doublereal* acMolality) const;
529 
530  //@}
531  /// @name Partial Molar Properties of the Solution -----------------
532  //@{
533 
534 
535  //!Get the species chemical potentials: Units: J/kmol.
536  /*!
537  *
538  * This function returns a vector of chemical potentials of the
539  * species in solution.
540  *
541  * \f[
542  * \mu_k = \mu^{o}_k(T,P) + R T \ln(\frac{m_k}{m^\Delta})
543  * \f]
544  * \f[
545  * \mu_w = \mu^{o}_w(T,P) +
546  * R T ((X_w - 1.0) / X_w)
547  * \f]
548  *
549  * \f$ w \f$ refers to the solvent species.
550  * \f$ X_w \f$ is the mole fraction of the solvent.
551  * \f$ m_k \f$ is the molality of the kth solute.
552  * \f$ m^\Delta is 1 gmol solute per kg solvent. \f$
553  *
554  * Units: J/kmol.
555  *
556  * @param mu Output vector of species chemical potentials.
557  * Length: m_kk.
558  */
559  virtual void getChemPotentials(doublereal* mu) const;
560 
561  //! Returns an array of partial molar enthalpies for the species in the mixture.
562  /*!
563  * Units (J/kmol)
564  * For this phase, the partial molar enthalpies are equal to the
565  * species standard state enthalpies.
566  * \f[
567  * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k
568  * \f]
569  * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_k(T) \f$,
570  * at the reference pressure,\f$ P_{ref} \f$,
571  * are computed by the species thermodynamic
572  * property manager. They are polynomial functions of temperature.
573  * @see SpeciesThermo
574  *
575  * @param hbar Output vector of partial molar enthalpies.
576  * Length: m_kk.
577  */
578  virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
579 
580 
581  //! Returns an array of partial molar entropies of the species in the solution. Units: J/kmol.
582  /*!
583  *
584  * Maxwell's equations provide an insight in how to calculate this
585  * (p.215 Smith and Van Ness)
586  * \f[
587  * \frac{d(\mu_k)}{dT} = -\bar{s}_i
588  * \f]
589  * For this phase, the partial molar entropies are equal to the
590  * standard state species entropies plus the ideal molal solution contribution.
591  *
592  * \f[
593  * \bar{s}_k(T,P) = s^0_k(T) - R \ln( \frac{m_k}{m^{\triangle}} )
594  * \f]
595  * \f[
596  * \bar{s}_w(T,P) = s^0_w(T) - R ((X_w - 1.0) / X_w)
597  * \f]
598  *
599  * The subscript, w, refers to the solvent species. \f$ X_w \f$ is
600  * the mole fraction of solvent.
601  * The reference-state pure-species entropies,\f$ s^0_k(T) \f$,
602  * at the reference pressure, \f$ P_{ref} \f$, are computed by the
603  * species thermodynamic
604  * property manager. They are polynomial functions of temperature.
605  * @see SpeciesThermo
606  *
607  * @param sbar Output vector of partial molar entropies.
608  * Length: m_kk.
609  */
610  virtual void getPartialMolarEntropies(doublereal* sbar) const;
611 
612  // partial molar volumes of the species Units: m^3 kmol-1.
613  /*!
614  * For this solution, the partial molar volumes are equal to the
615  * constant species molar volumes.
616  *
617  * Units: m^3 kmol-1.
618  * @param vbar Output vector of partial molar volumes.
619  */
620  virtual void getPartialMolarVolumes(doublereal* vbar) const;
621 
622 
623  //! Partial molar heat capacity of the solution:. UnitsL J/kmol/K
624  /*!
625  * The kth partial molar heat capacity is equal to
626  * the temperature derivative of the partial molar
627  * enthalpy of the kth species in the solution at constant
628  * P and composition (p. 220 Smith and Van Ness).
629  * \f[
630  * \bar{Cp}_k(T,P) = {Cp}^0_k(T)
631  * \f]
632  *
633  * For this solution, this is equal to the reference state
634  * heat capacities.
635  *
636  * Units: J/kmol/K
637  *
638  * @param cpbar Output vector of partial molar heat capacities.
639  * Length: m_kk.
640  */
641  virtual void getPartialMolarCp(doublereal* cpbar) const;
642 
643  //@}
644  /// @name Properties of the Standard State of the Species
645  // in the Solution --
646  //@{
647 
648 
649 
650 
651  //@}
652  /// @name Thermodynamic Values for the Species Reference States ---
653  //@{
654 
655 
656  ///////////////////////////////////////////////////////
657  //
658  // The methods below are not virtual, and should not
659  // be overloaded.
660  //
661  //////////////////////////////////////////////////////
662 
663  /**
664  * @name Specific Properties
665  * @{
666  */
667 
668 
669  /**
670  * @name Setting the State
671  *
672  * These methods set all or part of the thermodynamic
673  * state.
674  * @{
675  */
676 
677  //@}
678 
679  /**
680  * @name Chemical Equilibrium
681  * Chemical equilibrium.
682  * @{
683  */
684 
685  /**
686  * This method is used by the ChemEquil equilibrium solver.
687  * It sets the state such that the chemical potentials satisfy
688  * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
689  * \left(\frac{\lambda_m} {\hat R T}\right) \f] where
690  * \f$ \lambda_m \f$ is the element potential of element m. The
691  * temperature is unchanged. Any phase (ideal or not) that
692  * implements this method can be equilibrated by ChemEquil.
693  *
694  * Not implemented.
695  *
696  * @param lambda_RT vector of Nondimensional element potentials.
697  */
698  virtual void setToEquilState(const doublereal* lambda_RT) {
699  err("setToEquilState");
700  }
701 
702  //@}
703 
704  /**
705  * @internal
706  * Set equation of state parameters. The number and meaning of
707  * these depends on the subclass.
708  * @param n number of parameters
709  * @param c array of <I>n</I> coefficients
710  *
711  */
712  virtual void setParameters(int n, doublereal* const c);
713 
714  /*!
715  * @internal
716  * Get the parameters used to initialize the phase.
717  *
718  * @param n number of parameters (output)
719  * @param c array of <I>n</I> coefficients
720  */
721  virtual void getParameters(int& n, doublereal* const c) const;
722 
723  /*!
724  * Set equation of state parameter values from XML
725  * entries. This method is called by function importPhase in
726  * file importCTML.cpp when processing a phase definition in
727  * an input file. It should be overloaded in subclasses to set
728  * any parameters that are specific to that particular phase
729  * model.
730  *
731  * @param eosdata An XML_Node object corresponding to
732  * the "thermo" entry for this phase in the input file.
733  */
734  virtual void setParametersFromXML(const XML_Node& eosdata);
735 
736  //---------------------------------------------------------
737  /// @name Critical state properties.
738  /// These methods are only implemented by some subclasses.
739 
740  //@{
741 
742  /**
743  * Critical temperature (K).
744  * Not implemented for this phase type.
745  */
746  virtual doublereal critTemperature() const {
747  err("critTemperature");
748  return -1.0;
749  }
750 
751  /**
752  * Critical pressure (Pa).
753  *
754  * Not implemented for this phase type.
755  */
756  virtual doublereal critPressure() const {
757  err("critPressure");
758  return -1.0;
759  }
760 
761  /**
762  * Critical density (kg/m3).
763  * Not implemented for this phase type.
764  */
765  virtual doublereal critDensity() const {
766  err("critDensity");
767  return -1.0;
768  }
769 
770  //@}
771 
772  /// @name Saturation properties.
773  /// These methods are only implemented by subclasses that
774  /// implement full liquid-vapor equations of state.
775  ///
776 
777  //@}
778 
779 
780  /*
781  * -------------- Utilities -------------------------------
782  */
783 
784  //! Initialization routine for an IdealMolalSoln phase.
785  /*!
786  * This internal routine is responsible for setting up
787  * the internal storage. This is reimplemented from the ThermoPhase
788  * class.
789  */
790  virtual void initThermo();
791 
792  //! Import and initialize an IdealMolalSoln phase
793  //! specification in an XML tree into the current object.
794  /*!
795  * Here we read an XML description of the phase.
796  * We import descriptions of the elements that make up the
797  * species in a phase.
798  * We import information about the species, including their
799  * reference state thermodynamic polynomials. We then freeze
800  * the state of the species.
801  *
802  * Then, we read the species molar volumes from the xml
803  * tree to finish the initialization.
804  *
805  * This routine is a precursor to constructPhaseXML(XML_Node*)
806  * routine, which does most of the work.
807  *
808  * This is a virtual routine, first used here.
809  *
810  * @param infile XML file containing the description of the
811  * phase
812  *
813  * @param id Optional parameter identifying the name of the
814  * phase. If none is given, the first XML
815  * phase element will be used.
816  */
817  void constructPhaseFile(std::string infile, std::string id="");
818 
819  //! Import and initialize an IdealMolalSoln phase
820  //! specification in an XML tree into the current object.
821  /*!
822  * This is the main routine for constructing the phase.
823  * It processes the XML file, and then it calls importPhase().
824  * Then, initThermoXML() is called after importPhase().
825  *
826  * Here we read an XML description of the phase.
827  * We import descriptions of the elements that make up the
828  * species in a phase.
829  * We import information about the species, including their
830  * reference state thermodynamic polynomials. We then freeze
831  * the state of the species.
832  *
833  * Then, we read the species molar volumes from the xml
834  * tree to finish the initialization.
835  *
836  * This is a virtual routine, first used in this class.
837  *
838  * @param phaseNode This object must be the phase node of a
839  * complete XML tree
840  * description of the phase, including all of the
841  * species data. In other words while "phase" must
842  * point to an XML phase object, it must have
843  * sibling nodes "speciesData" that describe
844  * the species in the phase.
845  * @param id ID of the phase. If nonnull, a check is done
846  * to see if phaseNode is pointing to the phase
847  * with the correct id.
848  */
849  void constructPhaseXML(XML_Node& phaseNode, std::string id);
850 
851  //! Import and initialize an IdealMolalSoln phase
852  //! specification in an XML tree into the current object.
853  /*!
854  * This routine is called from importPhase() to finish
855  * up the initialization of the thermo object. It reads in the
856  * species molar volumes.
857  *
858  * @param phaseNode This object must be the phase node of a
859  * complete XML tree
860  * description of the phase, including all of the
861  * species data. In other words while "phase" must
862  * point to an XML phase object, it must have
863  * sibling nodes "speciesData" that describe
864  * the species in the phase.
865  * @param id ID of the phase. If nonnull, a check is done
866  * to see if phaseNode is pointing to the phase
867  * with the correct id.
868  */
869  virtual void initThermoXML(XML_Node& phaseNode, std::string id="");
870 
871  //! Report the molar volume of species k
872  /*!
873  *
874  * units - \f$ m^3 kmol^-1 \f$
875  *
876  * @param k Species index.
877  */
878  double speciesMolarVolume(int k) const;
879 
880  /*!
881  * Fill in a return vector containing the species molar volumes
882  * units - \f$ m^3 kmol^-1 \f$
883  *
884  * @param smv Output vector of species molar volumes.
885  */
886  void getSpeciesMolarVolumes(double* smv) const;
887  //@}
888 
889 protected:
890  /**
891  * Species molar volume \f$ m^3 kmol^-1 \f$
892  */
894 
895  /**
896  * The standard concentrations can have three different forms
897  * depending on the value of the member attribute m_formGC, which
898  * is supplied in the XML file.
899  *
900  * <TABLE>
901  * <TR><TD> m_formGC </TD><TD> ActivityConc </TD><TD> StandardConc </TD></TR>
902  * <TR><TD> 0 </TD><TD> \f$ {m_k}/ { m^{\Delta}}\f$ </TD><TD> \f$ 1.0 \f$ </TD></TR>
903  * <TR><TD> 1 </TD><TD> \f$ m_k / (m^{\Delta} V_k)\f$ </TD><TD> \f$ 1.0 / V_k \f$ </TD></TR>
904  * <TR><TD> 2 </TD><TD> \f$ m_k / (m^{\Delta} V^0_0)\f$</TD><TD> \f$ 1.0 / V^0_0\f$ </TD></TR>
905  * </TABLE>
906  */
907  int m_formGC;
908 
909 public:
910  //! Cutoff type
912 
913 private:
914 
915  /**
916  * Vector containing the species reference exp(-G/RT) functions
917  * at T = m_tlast
918  */
920 
921  /**
922  * Vector of potential energies for the species.
923  */
924  mutable vector_fp m_pe;
925 
926  /**
927  * Temporary array used in equilibrium calculations
928  */
929  mutable vector_fp m_pp;
930 
931  /**
932  * vector of size m_kk, used as a temporary holding area.
933  */
934  mutable vector_fp m_tmpV;
935 
936  //! Logarithm of the molal activity coefficients
937  /*!
938  * Normally these are all one. However, stability schemes will change that
939  */
941 public:
942  //! value of the solute mole fraction that centers the cutoff polynomials
943  //! for the cutoff =1 process;
944  doublereal IMS_X_o_cutoff_;
945 
946  //! gamma_o value for the cutoff process at the zero solvent point
947  doublereal IMS_gamma_o_min_;
948 
949  //! gamma_k minimum for the cutoff process at the zero solvent point
950  doublereal IMS_gamma_k_min_;
951 
952  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
953  doublereal IMS_cCut_;
954 
955  //! Parameter in the polyExp cutoff treatment
956  /*!
957  * This is the slope of the f function at the zero solvent point
958  * Default value is 0.6
959  */
960  doublereal IMS_slopefCut_;
961 
962  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
963  doublereal IMS_dfCut_;
964 
965  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
966  doublereal IMS_efCut_;
967 
968  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
969  doublereal IMS_afCut_;
970 
971  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
972  doublereal IMS_bfCut_;
973 
974  //! Parameter in the polyExp cutoff treatment
975  /*!
976  * This is the slope of the g function at the zero solvent point
977  * Default value is 0.0
978  */
979  doublereal IMS_slopegCut_;
980 
981  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
982  doublereal IMS_dgCut_;
983 
984  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
985  doublereal IMS_egCut_;
986 
987  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
988  doublereal IMS_agCut_;
989 
990  //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
991  doublereal IMS_bgCut_;
992 
993 private:
994 
995  //! Internal error message
996  /*!
997  * @param msg message to be printed
998  */
999  doublereal err(std::string msg) const;
1000 
1001  //! This function will be called to update the internally stored
1002  //! natural logarithm of the molality activity coefficients
1003  /*!
1004  * Normally the solutes are all zero. However, sometimes they are not,
1005  * due to stability schemes
1006  */
1007  void s_updateIMS_lnMolalityActCoeff() const;
1008 
1009  //! This internal function adjusts the lengths of arrays.
1010  /*!
1011  * This function is not virtual nor is it inherited
1012  */
1013  void initLengths();
1014 
1015  //! Calculate parameters for cutoff treatments of activity coefficients
1016  /*!
1017  * Some cutoff treatments for the activity coefficients
1018  * actually require some calculations to create a consistent treatment.
1019  *
1020  * This routine is called during the setup to calculate these parameters
1021  */
1022  void calcIMSCutoffParams_();
1023 };
1024 
1025 /* @} */
1026 }
1027 
1028 #endif
1029 
1030 
1031 
1032 
1033