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

Cleaned up Doxygen docs for other descendants of ThermoPhase

Browse source
This commit is contained in:
Ray Speth 2013-02-14 01:02:58 +00:00
parent b15640a741
commit d59cac4ca0
15 changed files with 508 additions and 2577 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

44
include/cantera/thermo/ConstDensityThermo.h
View file

@ -22,31 +22,15 @@ namespace Cantera
//! Overloads the virtual methods of class ThermoPhase to implement the
//! incompressible equation of state.
/**
*
*
* <b> Specification of Species Standard State Properties </b>
*
*
* <b> Specification of Solution Thermodynamic Properties </b>
*
* The density is assumed to be constant, no matter what the concentration of the solution.
*
*
* <b> Application within %Kinetics Managers </b>
*
*
* <b> XML Example </b>
*
* An example of an XML Element named phase setting up a SurfPhase object named diamond_100
* is given below.
* The density is assumed to be constant, no matter what the concentration of the solution.
*
* @ingroup thermoprops
*/
class ConstDensityThermo : public ThermoPhase
{
public:
//! Constructor.
ConstDensityThermo();
@ -67,13 +51,13 @@ public:
//! Duplication routine for objects which inherit from %ThermoPhase
/*!
* This virtual routine can be used to duplicate %ThermoPhase objects
* inherited from %ThermoPhase even if the application only has
* This virtual routine can be used to duplicate objects
* derived from %ThermoPhase even if the application only has
* a pointer to %ThermoPhase to work with.
*/
virtual ThermoPhase* duplMyselfAsThermoPhase() const;
//! overloaded methods of class ThermoPhase
//! Returns a constant corresponding to this class's equation of state
virtual int eosType() const;
/// Molar enthalpy. Units: J/kmol.
@ -95,26 +79,11 @@ public:
virtual doublereal cv_mole() const;
//! Return the thermodynamic pressure (Pa).
/*!
* This method must be overloaded in derived classes. Since the
* mass density, temperature, and mass fractions are stored,
* this method should use these values to implement the
* mechanical equation of state \f$ P(T, \rho, Y_1, \dots,
* Y_K) \f$.
*/
virtual doublereal pressure() const;
//! Set the internally stored pressure (Pa) at constant
//! temperature and composition
/*!
* This method must be reimplemented in derived classes, where it
* may involve the solution of a nonlinear equation. Within %Cantera,
* the independent variable is the density. Therefore, this function
* solves for the density that will yield the desired input pressure.
* The temperature and composition iare held constant during this process.
*
* This base class function will print an error, if not overwritten.
*
* @param p input Pressure (Pa)
*/
virtual void setPressure(doublereal p);
@ -250,9 +219,6 @@ public:
std::copy(_cpr.begin(), _cpr.end(), cpr);
}
// new methods defined here
//! Returns a reference to the vector of nondimensional
//! enthalpies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
@ -319,7 +285,6 @@ public:
*/
virtual void setToEquilState(const doublereal* lambda_RT);
//! Set the equation of state parameters
/*!
* @internal
@ -384,7 +349,6 @@ protected:
doublereal m_press;
private:
//! Function to update the reference state thermo functions
void _updateThermo() const;
};

110
include/cantera/thermo/IdealGasPhase.h
View file

@ -97,10 +97,10 @@ namespace Cantera
* The molar volume of a species is given by the ideal gas law
*
* \f[
* V^o_k(T,P) = \frac{R T}{P} \mbox{\quad where}
* V^o_k(T,P) = \frac{R T}{P}
* \f]
*
* R is the molar gas constant. For a complete list of physical constants
* where R is the molar gas constant. For a complete list of physical constants
* used within %Cantera, see \ref physConstants .
*
* <HR>
@ -183,7 +183,7 @@ namespace Cantera
* \f]
* where
* \f[
* C_j^a = C^s a_j \mbox{\quad and \quad} C_k^a = C^s a_k
* C_j^a = C^s a_j \quad \mbox{and} \quad C_k^a = C^s a_k
* \f]
*
* \f$ C_j^a \f$ is the activity concentration of species j, and
@ -228,7 +228,8 @@ namespace Cantera
* For completeness, the pressure equilibrium constant may be obtained as well
*
* \f[
* \frac{P_j P_k}{ P_l P_{ref}} = K_p^1 = \exp(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} )
* \frac{P_j P_k}{ P_l P_{ref}} = K_p^1 =
\exp\left(\frac{\mu^{ref}_l - \mu^{ref}_j - \mu^{ref}_k}{R T} \right)
* \f]
*
* \f$ K_p \f$ is the simplest form of the equilibrium constant for ideal gases. However, it isn't
@ -255,7 +256,6 @@ namespace Cantera
* <H2> Instantiation of the Class </H2>
* <HR>
*
*
* The constructor for this phase is located in the default ThermoFactory
* for %Cantera. A new %IdealGasPhase may be created by the following code
* snippet:
@ -281,20 +281,20 @@ namespace Cantera
* An example of an XML Element named phase setting up a IdealGasPhase
* object named silane is given below.
*
* @verbatim
<!-- phase silane -->
<phase dim="3" id="silane">
<elementArray datasrc="elements.xml"> Si H He </elementArray>
<speciesArray datasrc="#species_data">
H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6
H2SISIH2 SI3H8 SI2 SI3
</speciesArray>
<reactionArray datasrc="#reaction_data"/>
<thermo model="IdealGas"/>
<kinetics model="GasKinetics"/>
<transport model="None"/>
</phase>
@endverbatim
* @code
* <!-- phase silane -->
* <phase dim="3" id="silane">
* <elementArray datasrc="elements.xml"> Si H He </elementArray>
* <speciesArray datasrc="#species_data">
* H2 H HE SIH4 SI SIH SIH2 SIH3 H3SISIH SI2H6
* H2SISIH2 SI3H8 SI2 SI3
* </speciesArray>
* <reactionArray datasrc="#reaction_data"/>
* <thermo model="IdealGas"/>
* <kinetics model="GasKinetics"/>
* <transport model="None"/>
* </phase>
* @endcode
*
* The model attribute "IdealGas" of the thermo XML element identifies the phase as
* being of the type handled by the IdealGasPhase object.
@ -304,22 +304,13 @@ namespace Cantera
*/
class IdealGasPhase: public ThermoPhase
{
public:
//! Default empty Constructor
IdealGasPhase();
//! Construct and initialize an IdealGasPhase ThermoPhase object
//! directly from an ASCII input file
/*!
* Working constructors
*
* The two constructors below are a direct way that
* the phase can initialize itself. They are shells that call
* the routine initThermo(), with a reference to the
* XML database to get the info for the phase.
*
* @param inputFile Name of the input file containing the phase XML data
* to set up the object
* @param id ID of the phase in the input file. Defaults to the
@ -360,10 +351,10 @@ public:
//! Destructor
virtual ~IdealGasPhase();
//! Duplicator from the %ThermoPhase parent class
//! Duplicator from the ThermoPhase parent class
/*!
* Given a pointer to a %ThermoPhase object, this function will
* duplicate the %ThermoPhase object and all underlying structures.
* Given a pointer to a ThermoPhase object, this function will
* duplicate the ThermoPhase object and all underlying structures.
* This is basically a wrapper around the inherited copy constructor.
*
* @return returns a pointer to a %ThermoPhase object, containing
@ -379,10 +370,8 @@ public:
return cIdealGas;
}
/**
* @name Molar Thermodynamic Properties of the Solution ------------------------------
* @{
*/
//! @name Molar Thermodynamic Properties of the Solution
//! @{
//! Return the Molar enthalpy. Units: J/kmol.
/*!
@ -496,19 +485,25 @@ public:
/**
* @returns species vibrational specific heat at
* constant volume. This is defined as
* constant volume,
* \f[
* C^{vib}_{v,s} \equiv \frac{\partial e^{vib}_s}{\partial T_V} = \frac{R_s \theta_{vs}^2 \exp\left(\theta_{vs}/T_V\right)}{\left[\left(\exp\left(\theta_{vs}/T_V\right)-1\right)T_V\right]^2}
* C^{vib}_{v,s} = \frac{\partial e^{vib}_{v,s} }{\partial T}
* \f]
* where the species vibration energy \f$ e^{vib}_{v,s} \f$ is
* - atom:
* 0
* - Diatomic:
* \f[ \frac{R_s \theta_{v,s}}{e^{\theta_{v,s}/T}-1} \f]
* - General Molecule:
* \f[
* \sum_i \frac{R_s \theta_{v,s,i}}{e^{\theta_{v,s,i}/T}-1}
* \f]
*/
virtual doublereal cv_vib(int k, doublereal T) const;
//@}
/**
* @name Mechanical Equation of State ------------------------------------------------
* @{
*/
//! @}
//! @name Mechanical Equation of State
//! @{
/**
* Pressure. Units: Pa.
@ -560,8 +555,7 @@ public:
//@}
/**
* @name Chemical Potentials and Activities ------------------------------------------
*
* @name Chemical Potentials and Activities
*
* The activity \f$a_k\f$ of a species in solution is
* related to the chemical potential by
@ -639,7 +633,7 @@ public:
virtual void getActivityCoefficients(doublereal* ac) const;
//@}
/// @name Partial Molar Properties of the Solution ----------------------------------
/// @name Partial Molar Properties of the Solution
//@{
//! Get the species chemical potentials. Units: J/kmol.
@ -689,7 +683,7 @@ public:
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution ----------
/// @name Properties of the Standard State of the Species in the Solution
//@{
//! Get the array of chemical potentials at unit activity for the
@ -765,7 +759,7 @@ public:
virtual void getStandardVolumes(doublereal* vol) const;
//@}
/// @name Thermodynamic Values for the Species Reference States ---------------------
/// @name Thermodynamic Values for the Species Reference States
//@{
//! Returns the vector of nondimensional
@ -906,11 +900,9 @@ public:
*/
virtual void initThermo();
//!This method is used by the ChemEquil equilibrium solver.
//! Method used by the ChemEquil equilibrium solver.
/*!
* @internal
* @name Chemical Equilibrium
* @{
*
* Set mixture to an equilibrium state consistent with specified
* element potentials and temperature.
@ -926,8 +918,6 @@ public:
*/
virtual void setToEquilState(const doublereal* lambda_RT);
//@}
protected:
//! Reference state pressure
/*!
@ -960,15 +950,19 @@ protected:
mutable vector_fp m_pp;
private:
//! Update the species reference state thermodynamic functions
/*!
* The polynomials for the standard state functions are only
* reevaluated if the temperature has changed.
*
* This method is called each time a thermodynamic property is requested,
* to check whether the internal species properties within the object
* need to be updated. Currently, this updates the species thermo
* polynomial values for the current value of the temperature. A check is
* made to see if the temperature has changed since the last evaluation.
* This object does not contain any persistent data that depends on the
* concentration, that needs to be updated. The state object modifies its
* concentration dependent information at the time the setMoleFractions()
* (or equivalent) call is made.
*/
void _updateThermo() const;
};
}

246
include/cantera/thermo/IdealSolidSolnPhase.h
View file

@ -22,7 +22,6 @@
#include "ThermoFactory.h"
#include "SpeciesThermo.h"
namespace Cantera
{
@ -45,17 +44,10 @@ const int cIdealSolidSolnPhase2 = 5012;
* The class derives from class ThermoPhase,
* and overloads the virtual methods defined there with ones that
* use expressions appropriate for ideal solution mixtures.
* File name for the XML datafile containing information
* for this phase
*
* The generalized concentrations can have three different forms
* depending on the value of the member attribute m_formGC, which
* depending on the value of the member attribute #m_formGC, which
* is supplied in the constructor and in the XML file.
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
* The value and form of the generalized concentration will affect
* reaction rate constants involving species in this phase.
*
@ -63,49 +55,30 @@ const int cIdealSolidSolnPhase2 = 5012;
*/
class IdealSolidSolnPhase : public ThermoPhase
{
public:
/**
* Constructor for IdealSolidSolnPhase.
* The generalized concentrations can have three different forms
* depending on the value of the member attribute m_formGC, which
* depending on the value of the member attribute #m_formGC, which
* is supplied in the constructor or read from the xml data file.
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
*
* @param formCG This parameter initializes the m_formGC variable. The default
* is a value of 0.
* @param formCG This parameter initializes the #m_formGC variable.
*/
IdealSolidSolnPhase(int formCG=0);
//! Construct and initialize an IdealSolidSolnPhase ThermoPhase object
//! directly from an ASCII input file
/*!
*
* This constructor will also fully initialize the object.
* The generalized concentrations can have three different forms
* depending on the value of the member attribute m_formGC, which
* depending on the value of the member attribute #m_formGC, which
* is supplied in the constructor or read from the xml data file.
*
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
*
* @param infile File name for the XML datafile containing information
* for this phase
* @param id The name of this phase. This is used to look up
* the phase in the XML datafile.
* @param formCG This parameter initializes the m_formGC variable. The default
* is a value of 0.
* @param formCG This parameter initializes the #m_formGC variable.
*/
IdealSolidSolnPhase(const std::string& infile, const std::string& id="", int formCG=0);
@ -113,42 +86,28 @@ public:
//! directly from an XML database
/*!
* The generalized concentrations can have three different forms
* depending on the value of the member attribute m_formGC, which
* depending on the value of the member attribute #m_formGC, which
* is supplied in the constructor and/or read from the data file.
*
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
*
* @param root XML tree containing a description of the phase.
* The tree must be positioned at the XML element
* named phase with id, "id", on input to this routine.
* @param id The name of this phase. This is used to look up
* the phase in the XML datafile.
* @param formCG This parameter initializes the m_formGC variable. The default
* is a value of 0.
* @param formCG This parameter initializes the #m_formGC variable.
*/
IdealSolidSolnPhase(XML_Node& root, const std::string& id="", int formCG=0);
/*!
* Copy Constructor
*/
//! Copy Constructor
IdealSolidSolnPhase(const IdealSolidSolnPhase&);
/*!
* Assignment operator
*/
//! Assignment operator
IdealSolidSolnPhase& operator=(const IdealSolidSolnPhase&);
/*!
* Base Class Duplication Function
* -> given a pointer to ThermoPhase, this function can
* duplicate the object. (note has to be a separate function
* not the copy constructor, because it has to be
* a virtual function)
*
* Given a pointer to ThermoPhase, this function can duplicate the object.
*/
virtual ThermoPhase* duplMyselfAsThermoPhase() const;
@ -157,14 +116,12 @@ public:
/**
* Equation of state flag. Returns a value depending upon the value of
* m_formGC, which is defined at instantiation.
* #m_formGC, which is defined at instantiation.
*/
virtual int eosType() const;
/**
* @name Molar Thermodynamic Properties of the Solution ------------------------
* @{
*/
//! @name Molar Thermodynamic Properties of the Solution
//! @{
/**
* Molar enthalpy of the solution. Units: J/kmol.
@ -257,7 +214,7 @@ public:
}
//@}
/** @name Mechanical Equation of State Properties ------------------------------------
/** @name Mechanical Equation of State Properties
*
* In this equation of state implementation, the density is a
* function only of the mole fractions. Therefore, it can't be
@ -305,9 +262,6 @@ public:
* species molar volumes. We have additionally specified
* in this class that the pure species molar volumes are
* independent of temperature and pressure.
*
* NOTE: This is a non-virtual function, which is not a
* member of the ThermoPhase base class.
*/
void calcDensity();
@ -322,11 +276,6 @@ public:
* to create a condition where the density is a function of
* the pressure.
*
* This function will now throw an error condition.
*
* NOTE: This is a virtual function that overwrites the State.h
* class
*
* @param rho Input density
*/
virtual void setDensity(const doublereal rho);
@ -337,9 +286,6 @@ public:
*
* This function will now throw an error condition.
*
* NOTE: This is virtual function that overwrites the State.h
* class
*
* @param rho Input Density
*/
virtual void setMolarDensity(const doublereal rho);
@ -379,11 +325,10 @@ public:
*/
virtual void setConcentrations(const doublereal* const c);
//@}
/**
* @name Chemical Potentials and Activities -----------------------------------------
* @name Chemical Potentials and Activities
*
* The activity \f$a_k\f$ of a species in solution is
* related to the chemical potential by
@ -443,14 +388,7 @@ public:
* standard concentration to be independent of the mole fractions.
*
* In this implementation the form of the generalized concentrations
* depend upon the member attribute, m_formGC:
*
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
* depend upon the member attribute, #m_formGC.
*
* HKM Note: We have absorbed the pressure dependence of the pure species
* state into the thermodynamics functions. Therefore the
@ -464,14 +402,12 @@ public:
virtual void getActivityConcentrations(doublereal* c) const;
/**
* The standard concentration \f$ C^0_k \f$ used to normalize
* the generalized concentration.
* In many cases, this quantity
* will be the same for all species in a phase.
* However, for this case, we will return a distinct concentration
* for each species. This is the inverse of the species molar
* volume. Units for the standard concentration are
* kmol m<SUP>-3</SUP>.
* The standard concentration \f$ C^0_k \f$ used to normalize the
* generalized concentration. In many cases, this quantity will be the
* same for all species in a phase. However, for this case, we will return
* a distinct concentration for each species. This is the inverse of the
* species molar volume. Units for the standard concentration are kmol
* m<SUP>-3</SUP>.
*
* @param k Species number: this is a require parameter,
* a change from the ThermoPhase base class, where it was
@ -511,14 +447,16 @@ public:
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* @param uA Output vector containing the units
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
* @param uA Output vector containing the units:
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*
* @param k species index. Defaults to 0.
* @param sizeUA output int containing the size of the vector.
* Currently, this is equal to 6.
@ -531,7 +469,6 @@ public:
virtual void getUnitsStandardConc(double* uA, int k = 0,
int sizeUA = 6) const;
//! Get the array of species activity coefficients
/*!
* @param ac output vector of activity coefficients. Length: m_kk
@ -572,10 +509,9 @@ public:
virtual void getChemPotentials_RT(doublereal* mu) const;
//@}
/// @name Partial Molar Properties of the Solution -----------------------------
/// @name Partial Molar Properties of the Solution
//@{
//! Returns an array of partial molar enthalpies for the species in the mixture.
/*!
* Units (J/kmol)
@ -637,10 +573,9 @@ public:
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution -------------------------------------
/// @name Properties of the Standard State of the Species in the Solution
//@{
/**
* Get the standard state chemical potentials of the species.
* This is the array of chemical potentials at unit activity
@ -660,7 +595,6 @@ public:
getPureGibbs(mu0);
}
//! Get the array of nondimensional Enthalpy functions for the standard state species
//! at the current <I>T</I> and <I>P</I> of the solution.
/*!
@ -679,7 +613,6 @@ public:
*/
void getEnthalpy_RT(doublereal* hrt) const;
//! Get the nondimensional Entropies for the species
//! standard states at the current T and P of the solution.
/*!
@ -727,21 +660,18 @@ public:
*/
virtual void getPureGibbs(doublereal* gpure) const;
//! Returns the vector of nondimensional
//! internal Energies of the standard state at the current
//! temperature and pressure of the solution for each species.
/*!
*
* @param urt Output vector of standard state nondimensional internal energies.
* Length: m_kk.
*/
virtual void getIntEnergy_RT(doublereal* urt) const;
/**
* Get the nondimensional heat capacity at constant pressure
* function for the species
* standard states at the current T and P of the solution.
* Get the nondimensional heat capacity at constant pressure function for
* the species standard states at the current T and P of the solution.
* \f[
* Cp^0_k(T,P) = Cp^{ref}_k(T)
* \f]
@ -757,8 +687,7 @@ public:
/**
* Get the molar volumes of each species in their standard
* states at the current
* <I>T</I> and <I>P</I> of the solution.
* states at the current <I>T</I> and <I>P</I> of the solution.
* units = m^3 / kmol
*
* @param vol Output vector of standard state volumes.
@ -766,12 +695,10 @@ public:
*/
virtual void getStandardVolumes(doublereal* vol) const;
//@}
/// @name Thermodynamic Values for the Species Reference States ------
/// @name Thermodynamic Values for the Species Reference States
//@{
/**
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
@ -884,8 +811,9 @@ public:
virtual doublereal potentialEnergy(int k) const {
return m_pe[k];
}
//@}
/// @name Utility Functions -----------------------------------------------
/// @name Utility Functions
//@{
/**
@ -896,45 +824,34 @@ public:
virtual void initThermo();
/**
* @internal
* Import and initialize a ThermoPhase object
* using an XML tree.
* Here we read extra information about the XML description
* of a phase. Regular information about elements and species
* and their reference state thermodynamic information
* have already been read at this point.
* For example, we do not need to call this function for
* ideal gas equations of state.
* This function is called from importPhase()
* after the elements and the
* species are initialized with default ideal solution
* level data.
* @internal Import and initialize a ThermoPhase object using an XML
* tree. Here we read extra information about the XML description of a
* phase. Regular information about elements and species and their
* reference state thermodynamic information have already been read at
* this point. For example, we do not need to call this function for
* ideal gas equations of state. This function is called from
* importPhase() after the elements and the species are initialized
* with default ideal solution level data.
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
* @param phaseNode This object must be the phase node of a complete XML
* tree description of the phase, including all of the
* species data. In other words while "phase" must point to
* an XML phase object, it must have sibling nodes
* "speciesData" that describe the species in the phase.
* @param id ID of the phase. If nonnull, a check is done to see if
* phaseNode is pointing to the phase with the correct id.
*/
virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
/**
* Set mixture to an equilibrium state consistent with specified
* element potentials and the temperature.
*
* @param lambda_RT vector of non-dimensional element potentials
* \f$ \lambda_m/RT \f$.
*
*/
virtual void setToEquilState(const doublereal* lambda_RT);
/**
* Report the molar volume of species k
*
@ -957,12 +874,18 @@ public:
//@}
protected:
/**
* Format for the generalized concentrations
* 0 = C_k = X_k. (default)
* 1 = C_k = X_k / V_k
* 2 = C_k = X_k / V_N
* Format for the generalized concentrations.
*
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 (default) </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
*
* The value and form of the generalized concentration will affect
* reaction rate constants involving species in this phase.
*/
int m_formGC;
@ -995,9 +918,7 @@ protected:
*/
mutable doublereal m_tlast;
/**
* Vector containing the species reference enthalpies at T = m_tlast
*/
//! Vector containing the species reference enthalpies at T = m_tlast
mutable vector_fp m_h0_RT;
/**
@ -1006,16 +927,10 @@ protected:
*/
mutable vector_fp m_cp0_R;
/**
* Vector containing the species reference Gibbs functions
* at T = m_tlast
*/
//! Vector containing the species reference Gibbs functions at T = m_tlast
mutable vector_fp m_g0_RT;
/**
* Vector containing the species reference entropies
* at T = m_tlast
*/
//! Vector containing the species reference entropies at T = m_tlast
mutable vector_fp m_s0_R;
/**
@ -1024,18 +939,14 @@ protected:
*/
mutable vector_fp m_expg0_RT;
/**
* Vector of potential energies for the species.
*/
//! Vector of potential energies for the species.
mutable vector_fp m_pe;
/**
* Temporary array used in equilibrium calculations
*/
//! Temporary array used in equilibrium calculations
mutable vector_fp m_pp;
private:
/// @name Utility Functions ------------------------------------------
/// @name Utility Functions
//@{
/**
* This function gets called for every call to functions in this
@ -1047,9 +958,7 @@ private:
*/
void _updateThermo() const;
/**
* This internal function adjusts the lengths of arrays
*/
//! This internal function adjusts the lengths of arrays
void initLengths();
//@}
@ -1057,8 +966,3 @@ private:
}
#endif

112
include/cantera/thermo/LatticePhase.h
View file

@ -24,17 +24,16 @@ namespace Cantera
//! A simple thermodynamic model for a bulk phase,
//! assuming a lattice of solid atoms
/*!
* The bulk consists of a matrix of equivalent sites whose molar density
* does not vary with temperature or pressure. The thermodynamics
* obeys the ideal solution laws. The phase and the pure species phases which
* The bulk consists of a matrix of equivalent sites whose molar density
* does not vary with temperature or pressure. The thermodynamics
* obeys the ideal solution laws. The phase and the pure species phases which
* comprise the standard states of the species are assumed to have
* zero volume expansivity and zero isothermal compressibility.
*
* The density of matrix sites is given by the variable \f$ C_o \f$,
* which has SI units of kmol m-3.
*
*
* <b> Specification of Species Standard %State Properties </b>
* <b> Specification of Species Standard State Properties </b>
*
* It is assumed that the reference state thermodynamics may be
* obtained by a pointer to a populated species thermodynamic property
@ -48,7 +47,7 @@ namespace Cantera
* which has a weak dependence on the system pressure, \f$P\f$.
*
* \f[
* \raggedright h^o_k(T,P) =
* h^o_k(T,P) =
* h^{ref}_k(T) + \left( \frac{P - P_{ref}}{C_o} \right)
* \f]
*
@ -73,7 +72,6 @@ namespace Cantera
* V^o_k(T,P) = \frac{1.0}{C_o}
* \f]
*
*
* <HR>
* <H2> Specification of Solution Thermodynamic Properties </H2>
* <HR>
@ -179,24 +177,23 @@ namespace Cantera
* \mu_l(T,P) = \mu^o_l(T, P) + R T \log(a_l)
* \f]
*
* The concentration equilibrium constant, \f$ K_c \f$, may be obtained by changing over
* to activity concentrations. When this is done:
* The concentration equilibrium constant, \f$ K_c \f$, may be obtained by changing over
* to activity concentrations. When this is done:
*
* \f[
* \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 =
* \exp(\frac{\mu^{o}_l - \mu^{o}_j - \mu^{o}_k}{R T} )
* \f]
* \f[
* \frac{C^a_j C^a_k}{ C^a_l} = C^o K_a^{o,1} = K_c^1 =
* \exp(\frac{\mu^{o}_l - \mu^{o}_j - \mu^{o}_k}{R T} )
* \f]
*
*
* %Kinetics managers will calculate the concentration equilibrium constant, \f$ K_c \f$,
* using the second and third part of the above expression as a definition for the concentration
* equilibrium constant.
* %Kinetics managers will calculate the concentration equilibrium constant, \f$ K_c \f$,
* using the second and third part of the above expression as a definition for the concentration
* equilibrium constant.
*
* <HR>
* <H2> Instantiation of the Class </H2>
* <HR>
*
*
* The constructor for this phase is located in the default ThermoFactory
* for %Cantera. A new %LatticePhase object may be created by the following code snippet:
*
@ -224,34 +221,31 @@ namespace Cantera
* An example of an XML Element named phase setting up a LatticePhase object named "O_lattice_SiO2"
* is given below.
*
* @verbatim
<!-- phase O_lattice_SiO2 -->
<phase dim="3" id="O_lattice_SiO2">
<elementArray datasrc="elements.xml"> Si H He </elementArray>
<speciesArray datasrc="#species_data">
O_O Vac_O
</speciesArray>
<reactionArray datasrc="#reaction_data"/>
<thermo model="Lattice">
<site_density> 73.159 </site_density>
<vacancy_species> Vac_O </vacancy_species>
</thermo>
<kinetics model="BulkKinetics"/>
<transport model="None"/>
</phase>
@endverbatim
* @code
* <!-- phase O_lattice_SiO2 -->
* <phase dim="3" id="O_lattice_SiO2">
* <elementArray datasrc="elements.xml"> Si H He </elementArray>
* <speciesArray datasrc="#species_data">
* O_O Vac_O
* </speciesArray>
* <reactionArray datasrc="#reaction_data"/>
* <thermo model="Lattice">
* <site_density> 73.159 </site_density>
* <vacancy_species> Vac_O </vacancy_species>
* </thermo>
* <kinetics model="BulkKinetics"/>
* <transport model="None"/>
* </phase>
* @endcode
*
* The model attribute "Lattice" of the thermo XML element identifies the phase as
* being of the type handled by the LatticePhase object.
*
* @ingroup thermoprops
*
*/
class LatticePhase : public ThermoPhase
{
public:
//! Base Empty constructor
LatticePhase();
@ -289,8 +283,6 @@ public:
* This virtual function is used to create a duplicate of the
* current phase. It's used to duplicate the phase when given
* a ThermoPhase pointer to the phase.
*
* @return It returns a ThermoPhase pointer.
*/
ThermoPhase* duplMyselfAsThermoPhase() const;
@ -299,10 +291,8 @@ public:
return cLattice;
}
/**
* @name Molar Thermodynamic Properties of the Solution ------------------------
* @{
*/
//! @name Molar Thermodynamic Properties of the Solution
//! @{
//! Return the Molar Enthalpy. Units: J/kmol.
/*!
@ -405,8 +395,7 @@ public:
virtual doublereal cv_mole() const;
//@}
/// @name Mechanical Equation of State Properties ------------------------------------
//@{
/// @name Mechanical Equation of State Properties
/**
* In this equation of state implementation, the density is a
* function only of the mole fractions. Therefore, it can't be
@ -417,7 +406,6 @@ public:
*/
//@{
//! Pressure. Units: Pa.
/*!
* For this incompressible system, we return the internally stored
@ -455,9 +443,6 @@ public:
* species molar volumes. We have additionally specified
* in this class that the pure species molar volumes are
* independent of temperature and pressure.
*
* NOTE: This is a non-virtual function, which is not a
* member of the ThermoPhase base class.
*/
doublereal calcDensity();
@ -496,11 +481,9 @@ public:
*/
virtual void setConcentrations(const doublereal* const c);
//@}
/// @name Activities, Standard States, and Activity Concentrations
/**
*
* The activity \f$a_k\f$ of a species in solution is
* related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
* + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
@ -564,7 +547,6 @@ public:
//@}
/// @name Partial Molar Properties of the Solution
///
//@{
//! Get the species chemical potentials. Units: J/kmol.
@ -578,12 +560,10 @@ public:
*/
virtual void getChemPotentials(doublereal* mu) const;
//@}
/// @name Partial Molar Properties of the Solution -----------------------------
/// @name Partial Molar Properties of the Solution
//@{
/**
* Returns an array of partial molar enthalpies for the species
* in the mixture.
@ -642,7 +622,6 @@ public:
*/
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//! Get the array of chemical potentials at unit activity for the
//! species standard states at the current <I>T</I> and <I>P</I> of the solution.
/*!
@ -664,7 +643,6 @@ public:
*/
virtual void getPureGibbs(doublereal* gpure) const;
//@}
/// @name Properties of the Standard State of the Species in the Solution
//@{
@ -844,7 +822,6 @@ public:
*/
virtual void initThermo();
//! Import and initialize a ThermoPhase object using an XML tree.
/*!
* Here we read extra information about the XML description
@ -912,17 +889,16 @@ public:
*
* eosdata points to the thermo block, and looks like this:
*
* @verbatim
<phase id="O_lattice_SiO2" >
<thermo model="Lattice">
<site_density units="kmol/m^3"> 73.159 </site_density>
<vacancy_species> "O_vacancy" </vacancy_species>
</thermo>
</phase> @endverbatim
*
* @code
* <phase id="O_lattice_SiO2" >
* <thermo model="Lattice">
* <site_density units="kmol/m^3"> 73.159 </site_density>
* <vacancy_species> "O_vacancy" </vacancy_species>
* </thermo>
* </phase>
* @endcode
*/
virtual void setParametersFromXML(const XML_Node& eosdata);
//@}
protected:
@ -953,7 +929,6 @@ protected:
//! Temporary storage for the reference state entropies at the current temperature
mutable vector_fp m_s0_R;
//! String name for the species which represents a vacancy
//! in the lattice
/*!
@ -975,10 +950,7 @@ protected:
*/
doublereal m_site_density;
// doublereal m_molar_lattice_volume;
private:
//! Update the species reference state thermodynamic functions
/*!
* The polynomials for the standard state functions are only

24
include/cantera/thermo/LatticeSolidPhase.h
View file

@ -32,7 +32,6 @@ namespace Cantera
* lattice phase and the molar densities of the sublattice and the molar density within the LatticeSolidPhase
* have the same values.
*
*
* The mole fraction vector is redefined witin the the LatticeSolidPhase object. Each of the mole
* fractions sum to one on each of the sublattices. The routine getMoleFraction() and setMoleFraction()
* have been redefined to use this convention.
@ -65,7 +64,6 @@ namespace Cantera
* solid. \f$ \theta_i \f$ is a fixed weighting factor for the ith lattice representing the lattice
* stoichiometric coefficient. For this object the \f$ \theta_i \f$ values are fixed.
*
*
* Let's take FeS2 as an example, which may be thought of as a combination of two lattices: Fe and S lattice.
* The Fe sublattice has a molar density of 1 gmol cm-3. The S sublattice has a molar density of 2 gmol cm-3.
* We then define the LatticeSolidPhase object as having a nominal composition of FeS2, and having a
@ -73,7 +71,6 @@ namespace Cantera
* associated with the sublattices. The Fe sublattice will have a weight of 1.0 associated with it. The
* S sublattice will have a weight of 2.0 associated with it.
*
*
* <HR>
* <H3> Specification of Solution Density Properties </H3>
* <HR>
@ -108,9 +105,7 @@ namespace Cantera
*/
class LatticeSolidPhase : public ThermoPhase
{
public:
//! Base empty constructor
LatticeSolidPhase();
@ -141,7 +136,7 @@ public:
//! Equation of state type flag.
/*!
* Redefine this to return cLatticeSolid, listed in mix_defs.h.
* Returns cLatticeSolid, listed in mix_defs.h.
*/
virtual int eosType() const {
return cLatticeSolid;
@ -175,7 +170,6 @@ public:
*/
virtual doublereal maxTemp(size_t k = npos) const;
//! Returns the reference pressure in Pa. This function is a wrapper
//! that calls the species thermo refPressure function.
virtual doublereal refPressure() const ;
@ -205,7 +199,6 @@ public:
*/
virtual doublereal enthalpy_mole() const;
//! Return the Molar Internal Energy. Units: J/kmol.
/*!
* The molar enthalpy is determined by the following formula, where \f$ \theta_n \f$ is the
@ -309,8 +302,6 @@ public:
* \f]
*
* where \f$ \rho_n \f$ is the density of the nth sublattice
*
* Note this is a nonvirtual function.
*/
doublereal calcDensity();
@ -370,7 +361,6 @@ public:
return err("not implemented");
}
//! Set the mass fractions to the specified values, and then
//! normalize them so that they sum to 1.0.
/*!
@ -386,7 +376,6 @@ public:
err("not implemented");
}
//! Set the mass fractions to the specified values without normalizing.
/*!
* This is useful when the normalization
@ -452,7 +441,6 @@ public:
*/
virtual void getChemPotentials(doublereal* mu) const;
//! Returns an array of partial molar enthalpies for the species in the mixture.
/*!
* Units (J/kmol)
@ -520,7 +508,6 @@ public:
* \f$. The values are evaluated at the current
* temperature and pressure of the solution.
*
*
* This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of
* the sublattice species.
*
@ -555,7 +542,7 @@ public:
*/
virtual doublereal logStandardConc(size_t k=0) const;
//@}
/// @name Thermodynamic Values for the Species Reference States --------------------
/// @name Thermodynamic Values for the Species Reference States
//@{
//! Returns the vector of nondimensional enthalpies of the reference state at the current
@ -571,7 +558,6 @@ public:
*/
virtual void getGibbs_RT_ref(doublereal* grt) const;
//! Returns the vector of the gibbs function of the reference state at the current
//! temperatureof the solution and the reference pressure for the species.
/*!
@ -586,7 +572,6 @@ public:
*/
virtual void getGibbs_ref(doublereal* g) const;
//! Initialize the ThermoPhase object after all species have been set up
/*!
* @internal Initialize.
@ -617,10 +602,8 @@ public:
*/
virtual void installSlavePhases(Cantera::XML_Node* phaseNode);
//! Set equation of state parameter values from XML entries.
/*!
*
* This method is called by function importPhase() in
* file importCTML.cpp when processing a phase definition in
* an input file. It should be overloaded in subclasses to set
@ -633,10 +616,8 @@ public:
*/
virtual void setParametersFromXML(const XML_Node& eosdata);
//! Set the Lattice mole fractions using a string
/*!
*
* @param n Integer value of the lattice whose mole fractions are being set
* @param x string containing Name:value pairs that will specify the mole fractions
* of species on a particular lattice
@ -697,7 +678,6 @@ protected:
std::vector<size_t> lkstart_;
private:
//! Update the reference thermodynamic functions
void _updateThermo() const;
};

185
include/cantera/thermo/MixtureFugacityTP.h
View file

@ -3,7 +3,6 @@
* Header file for a derived class of ThermoPhase that handles
* non-ideal mixtures based on the fugacity models (see \ref thermoprops and
* class \link Cantera::MixtureFugacityTP MixtureFugacityTP\endlink).
*
*/
/*
* Copyright (2005) Sandia Corporation. Under the terms of
@ -72,19 +71,14 @@ class PDSS;
*
* This class is usually used for non-ideal gases.
*
*
* @nosubgrouping
*/
class MixtureFugacityTP : public ThermoPhase
{
public:
//! @name Constructors and Duplicators for %MixtureFugacityTP
//! @{
/*!
*
* @name Constructors and Duplicators for %MixtureFugacityTP
*
*/
//! Constructor.
MixtureFugacityTP();
@ -103,19 +97,15 @@ public:
//! Destructor.
virtual ~MixtureFugacityTP();
//! Duplication routine
/*!
* @return Returns a duplication
*/
virtual ThermoPhase* duplMyselfAsThermoPhase() const;
//@}
/**
* @name Utilities (MixtureFugacityTP)
*/
//@{
//! @}
//! @name Utilities
//! @{
/**
* Equation of state type flag. The base class returns
* zero. Subclasses should define this to return a unique
@ -131,12 +121,10 @@ public:
//! temperature based, and variable pressure based.
/*!
* Currently, there are two standard state conventions:
* - Temperature-based activities
* cSS_CONVENTION_TEMPERATURE 0
* - default
*
* - Variable Pressure and Temperature -based activities
* cSS_CONVENTION_VPSS 1
* - Temperature-based activities,
* `cSS_CONVENTION_TEMPERATURE 0` (default)
* - Variable Pressure and Temperature based activities,
* `cSS_CONVENTION_VPSS 1`
*/
virtual int standardStateConvention() const;
@ -164,17 +152,15 @@ public:
*/
virtual int reportSolnBranchActual() const;
//! Get the array of log concentration-like derivatives of the
//! log activity coefficients
/*!
* This function is a virtual method. For ideal mixtures
* (unity activity coefficients), this can return zero.
* Implementations should take the derivative of the
* logarithm of the activity coefficient with respect to the
* logarithm of the concentration-like variable (i.e. moles)
* that represents the standard state.
* For ideal mixtures (unity activity coefficients), this can return zero.
* Implementations should take the derivative of the logarithm of the
* activity coefficient with respect to the logarithm of the
* concentration-like variable (i.e. moles) that represents the standard
* state.
*
* This quantity is to be used in conjunction with derivatives of
* that concentration-like variable when the derivative of the chemical
* potential is taken.
@ -188,12 +174,10 @@ public:
err("getdlnActCoeffdlnN_diag");
}
//@}
/// @name Partial Molar Properties of the Solution (MixtureFugacityTP)
/// @name Partial Molar Properties of the Solution
//@{
//! Get the array of non-dimensional species chemical potentials
//! These are partial molar Gibbs free energies.
/*!
@ -210,10 +194,8 @@ public:
void getChemPotentials_RT(doublereal* mu) const;
//@}
/*!
* @name Properties of the Standard State of the Species in the Solution
* (MixtureFugacityTP)
*
* Within MixtureFugacityTP, these properties are calculated via a common routine,
* _updateStandardStateThermo(),
@ -249,7 +231,6 @@ public:
*/
virtual void getEnthalpy_RT(doublereal* hrt) const;
//! Get the array of nondimensional Enthalpy functions for the standard state species
/*!
* at the current <I>T</I> and <I>P</I> of the solution.
@ -297,7 +278,6 @@ public:
*/
virtual void getIntEnergy_RT(doublereal* urt) const;
//! Get the nondimensional Heat Capacities at constant
//! pressure for the standard state of the species at the current T and P.
/*!
@ -311,7 +291,6 @@ public:
*/
virtual void getCp_R(doublereal* cpr) const;
//! Get the molar volumes of each species in their standard
//! states at the current <I>T</I> and <I>P</I> of the solution.
/*!
@ -324,7 +303,7 @@ public:
* units = m^3 / kmol
*/
virtual void getStandardVolumes(doublereal* vol) const;
// @}
//! Set the temperature of the phase
/*!
@ -336,7 +315,6 @@ public:
*/
virtual void setTemperature(const doublereal temp);
//! Set the internally stored pressure (Pa) at constant
//! temperature and composition
/*!
@ -348,8 +326,6 @@ public:
*/
virtual void setPressure(doublereal p);
protected:
/**
* Calculate the density of the mixture using the partial
@ -370,9 +346,6 @@ protected:
* species molar volumes. We have additionally specified
* in this class that the pure species molar volumes are
* independent of temperature and pressure.
*
* NOTE: This is a non-virtual function, which is not a
* member of the ThermoPhase base class.
*/
virtual void calcDensity();
@ -401,12 +374,10 @@ public:
*
* @param t Temperature (K)
* @param p Pressure (Pa)
* @param x Vector of mole fractions.
* Length is equal to m_kk.
* @param x Vector of mole fractions. Length is equal to m_kk.
*/
virtual void setState_TPX(doublereal t, doublereal p, const doublereal* x);
//! Set the mass fractions to the specified values, and then
//! normalize them so that they sum to 1.0.
/*!
@ -415,7 +386,6 @@ public:
*/
virtual void setMassFractions(const doublereal* const y);
//!Set the mass fractions to the specified values without normalizing.
/*!
* This is useful when the normalization
@ -423,13 +393,10 @@ public:
* by a constraint equation as part of a larger set of
* equations.
*
* @param y Input vector of mass fractions.
* Length is m_kk.
* @param y Input vector of mass fractions. Length is m_kk.
*/
virtual void setMassFractions_NoNorm(const doublereal* const y);
//! Set the mole fractions to the specified values, and then
//! normalize them so that they sum to 1.0.
/*!
@ -438,19 +405,16 @@ public:
*/
virtual void setMoleFractions(const doublereal* const x);
//! Set the mole fractions to the specified values without normalizing.
/*!
* This is useful when the normalization
* condition is being handled by some other means, for example
* by a constraint equation as part of a larger set of equations.
*
* @param x Input vector of mole fractions.
* Length is m_kk.
* @param x Input vector of mole fractions. Length is m_kk.
*/
virtual void setMoleFractions_NoNorm(const doublereal* const x);
//! Set the concentrations to the specified values within the phase.
/*!
* @param c The input vector to this routine is in dimensional
@ -465,7 +429,6 @@ public:
protected:
void setMoleFractions_NoState(const doublereal* const x);
public:
//! Returns the current pressure of the phase
/*!
@ -478,38 +441,24 @@ public:
return m_Pcurrent;
}
protected:
//! Updates the reference state thermodynamic functions at the current T of the solution.
/*!
*
* If m_useTmpStandardStateStorage is true,
* this function must be called for every call to functions in this
* class. It checks to see whether the temperature or pressure has changed and
* This function must be called for every call to functions in this
* class. It checks to see whether the temperature has changed and
* thus the ss thermodynamics functions for all of the species
* must be recalculated.
*
* This function is responsible for updating the following internal members,
* when m_useTmpStandardStateStorage is true.
*
* - m_hss_RT;
* - m_cpss_R;
* - m_gss_RT;
* - m_sss_R;
* - m_Vss
*
* If m_useTmpStandardStateStorage is not true, this function may be
* required to be called by child classes to update internal member data.
* This function is responsible for updating the following internal members:
*
* - m_h0_RT;
* - m_cp0_R;
* - m_g0_RT;
* - m_s0_R;
*/
virtual void _updateReferenceStateThermo() const;
public:
//@}
/// @name Thermodynamic Values for the Species Reference States (MixtureFugacityTP)
/*!
* There are also temporary
@ -520,7 +469,6 @@ public:
*/
//@{
//! Returns the vector of nondimensional
//! enthalpies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
@ -549,7 +497,6 @@ public:
//! Gibbs free energies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
/*!
*
* @param grt Output vector contains the nondimensional Gibbs free energies
* of the reference state of the species
* length = m_kk, units = dimensionless.
@ -613,16 +560,8 @@ public:
*/
virtual void getStandardVolumes_ref(doublereal* vol) const;
protected:
//@}
public:
//! @name Initialization Methods - For Internal use (VPStandardState)
//! @name Initialization Methods - For Internal use
/*!
* The following methods are used in the process of constructing
* the phase and setting its parameters from a specification in an
@ -632,13 +571,11 @@ public:
*/
//@{
//! Set the initial state of the phase to the conditions specified in the state XML element.
/*!
*
* This method sets the temperature, pressure, and mole fraction vector to a set default value.
*
* @param state AN XML_Node object corresponding to
* @param state An XML_Node object corresponding to
* the "state" entry for this phase in the input file.
*/
virtual void setStateFromXML(const XML_Node& state);
@ -690,17 +627,14 @@ public:
*/
virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
private:
//! @internal Initialize the internal lengths in this object.
/*!
* Note this is not a virtual function.
*/
void initLengths();
//@}
protected:
// Special Functions for fugacity classes
//! @name Special Functions for fugacity classes
//! @{
//! Calculate the value of z
/*!
@ -729,7 +663,6 @@ protected:
*/
virtual doublereal hresid() const;
//! Estimate for the saturation pressure
/*!
* Note: this is only used as a starting guess for later routines that actually calculate an
@ -740,6 +673,7 @@ protected:
* @return returns the estimated saturation pressure at the given temperature
*/
virtual doublereal psatEst(doublereal TKelvin) const;
public:
//! Estimate for the molar volume of the liquid
/*!
@ -757,7 +691,6 @@ public:
*/
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal& pres) const;
public:
//! Calculates the density given the temperature and the pressure and a guess at the density.
/*!
* Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this.
@ -777,26 +710,26 @@ public:
* @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there
* is no guessed density
*
*
* @return We return the density of the fluid at the requested phase. If we have not found any
* acceptable density we return a -1. If we have found an acceptable density at a
* different phase, we return a -2.
*/
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested,
doublereal rhoguess);
protected:
//! Utility routine in the calculation of the saturation pressure
/*!
* Private routine
*
* @param TKelvin temperature (kelvin)
* @param pres pressure (Pascal)
* @param densLiq Output density of liquid
* @param densGas output density of gas
* @param gasGRT output delGRT
* @param[out] densLiq density of liquid
* @param[out] densGas density of gas
* @param[out] liqGRT deltaG/RT of liquid
* @param[out] gasGRT deltaG/RT of gas
*/
int corr0(doublereal TKelvin, doublereal pres, doublereal& densLiq,
doublereal& densGas, doublereal& liqGRT, doublereal& gasGRT);
public:
//! Returns the Phase State flag for the current state of the object
/*!
@ -804,9 +737,9 @@ public:
* in parameters space we are
*
* 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_SUPERCRIT above the critical temperature
* - WATER_GAS below the critical temperature but below the critical density
* - WATER_LIQUID below the critical temperature but above the critical density
* - WATER_SUPERCRIT above the critical temperature
*/
int phaseState(bool checkState = false) const ;
@ -817,7 +750,6 @@ public:
*/
virtual doublereal densSpinodalLiquid() const;
//! Return the value of the density at the gas spinodal point (on the gas side)
//! for the current temperature.
/*!
@ -825,9 +757,6 @@ public:
*/
virtual doublereal densSpinodalGas() const;
public:
//! Calculate the saturation pressure at the current mixture content for the given temperature
/*!
@ -839,6 +768,7 @@ public:
*/
doublereal calculatePsat(doublereal TKelvin, doublereal& molarVolGas,
doublereal& molarVolLiquid);
protected:
//! Calculate the pressure given the temperature and the molar volume
/*!
@ -851,7 +781,6 @@ protected:
*/
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const;
//! Calculate the pressure and the pressure derivative given the temperature and the molar volume
/*!
* Temperature and mole number are held constant
@ -865,28 +794,19 @@ protected:
*/
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const;
virtual void updateMixingExpressions();
//@}
class spinodalFunc : public Cantera::ResidEval
{
public:
spinodalFunc(MixtureFugacityTP* tp);
virtual int evalSS(const doublereal t, const doublereal* const y, doublereal* const r);
MixtureFugacityTP* m_tp;
};
protected:
//! Current value of the pressures
/*!
* Because the pressure is now a calculation, we store the result of the calculation whenever
@ -896,29 +816,21 @@ protected:
*/
doublereal m_Pcurrent;
//! Storage for the current values of the mole fractions of the species
/*!
* This vector is kept up-to-date when some the setState functions are called.
*
* The State object is allowed to com
*
* Therefore, it may be considered to be an independent variable.
*
*/
std::vector<doublereal> moleFractions_;
//! Current state of the fluid
/*!
* There are three possible states of the fluid
* FLUID_GAS
* FLUID_LIQUID
* FLUID_SUPERCRIT
* There are three possible states of the fluid:
* - FLUID_GAS
* - FLUID_LIQUID
* - FLUID_SUPERCRIT
*/
int iState_;
//! Force the system to be on a particular side of the spinodal curve
int forcedState_;
@ -941,14 +853,13 @@ protected:
mutable vector_fp m_s0_R;
spinodalFunc* fdpdv_;
private:
private:
//! MixtureFugacityTP has its own err routine
/*!
* @param msg Error message string
*/
doublereal err(const std::string& msg) const;
};
}

35
include/cantera/thermo/PureFluidPhase.h
View file

@ -19,7 +19,6 @@
namespace Cantera
{
//! This phase object consists of a single component that can be a
//! gas, a liquid, a mixed gas-liquid fluid, or a fluid beyond its
//! critical point
@ -27,23 +26,10 @@ namespace Cantera
* The object inherits from ThermoPhase. However, it's built on top
* of the tpx package.
*
*
* <H2> Specification of Species Standard State Properties </H2>
*
*
* <H2> Application within %Kinetics Managers </H2>
*
*
* <H2> XML Example </H2>
*
*
* <H2> Instantiation of the Class </H2>
*
* @ingroup thermoprops
*/
class PureFluidPhase : public ThermoPhase
{
public:
//! Empty Base Constructor
@ -126,8 +112,6 @@ public:
mu[0] = gibbs_mole();
}
//! Get the species electrochemical potentials.
/*!
* These are partial molar quantities. This method adds a term \f$ F z_k
@ -159,7 +143,6 @@ public:
*/
virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
//! Returns an array of partial molar entropies of the species in the
//! solution. Units: J/kmol/K.
/*!
@ -344,9 +327,6 @@ public:
//! Returns the vector of nondimensional enthalpies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
/*!
* This base function will throw a Cantera exception unless
* it is overwritten in a derived class.
*
* @param hrt Output vector containing the nondimensional reference state enthalpies
* Length: m_kk.
*/
@ -370,8 +350,6 @@ public:
*/
virtual void getGibbs_ref(doublereal* g) const;
//! Returns the vector of nondimensional entropies of the reference state at the current temperature
//! of the solution and the reference pressure for each species.
/*!
@ -435,13 +413,9 @@ public:
*/
virtual void setState_SP(doublereal s, doublereal p,
doublereal tol = 1.e-8);
//@}
//! @name Critical State Properties
/*!
* Critical properties for the pure fluid
*/
//@{
//! critical temperature
@ -456,14 +430,8 @@ public:
//@}
//! @name Saturation properties.
/*!
* These methods are only implemented by subclasses that
* implement full liquid-vapor equations of state. They may be
* moved out of ThermoPhase at a later date.
*/
//@{
//! saturation temperature
/*!
* @param p Pressure (Pa)
@ -515,7 +483,6 @@ public:
//! Set equation of state parameter values from XML entries.
/*!
*
* This method is called by function importPhase() in
* file importCTML.cpp when processing a phase definition in
* an input file. It should be overloaded in subclasses to set
@ -528,7 +495,6 @@ public:
*/
virtual void setParametersFromXML(const XML_Node& eosdata);
//! returns a summary of the state of the phase as a string
/*!
* @param show_thermo If true, extra information is printed out
@ -550,7 +516,6 @@ protected:
void setTPXState() const;
private:
//! Pointer to the underlying tpx object Substance that does the work
mutable tpx::Substance* m_sub;

174
include/cantera/thermo/RedlichKwongMFTP.h
View file

@ -21,12 +21,6 @@ namespace Cantera
{
class XML_Node;
class PDSS;
/*!
* @name CONSTANTS - Models for the Standard State of IdealSolnPhase's
*/
//@{
/**
* @ingroup thermoprops
@ -34,39 +28,28 @@ class PDSS;
* This class can handle either an ideal solution or an ideal gas approximation
* of a phase.
*
*
* @nosubgrouping
*/
class RedlichKwongMFTP : public MixtureFugacityTP
{
public:
//! @name Constructors and Duplicators
//! @{
/*!
*
* @name Constructors and Duplicators for %RedlichKwongMFTP
*
*/
//! Base constructor.
RedlichKwongMFTP();
//! Construct and initialize a RedlichKwongMFTP ThermoPhase object
//! directly from an ASCII input file
//! Construct and initialize a RedlichKwongMFTP object directly from an
//! ASCII input file
/*!
* Working constructors
*
* The two constructors below are the normal way the phase initializes itself. They are shells that call
* the routine initThermo(), with a reference to the
* XML database to get the info for the phase.
*
* @param infile Name of the input file containing the phase XML data
* to set up the object
* @param id ID of the phase in the input file. Defaults to the empty string.
*/
RedlichKwongMFTP(const std::string& infile, std::string id="");
//! Construct and initialize a RedlichKwongMFTP ThermoPhase object
//! directly from an XML database
//! Construct and initialize a RedlichKwongMFTP object directly from an
//! XML database
/*!
* @param phaseRef XML phase node containing the description of the phase
* @param id id attribute containing the name of the phase. (default is the empty string)
@ -76,7 +59,6 @@ public:
//! This is a special constructor, used to replicate test problems
//! during the initial verification of the object
/*!
*
* test problems:
* 1: Pure CO2 problem
* input file = CO2_RedlickKwongMFTP.xml
@ -107,7 +89,6 @@ public:
//! Destructor.
virtual ~RedlichKwongMFTP();
//! Duplicator from the ThermoPhase parent class
/*!
* Given a pointer to a ThermoPhase object, this function will
@ -118,12 +99,6 @@ public:
*/
virtual ThermoPhase* duplMyselfAsThermoPhase() const;
//@}
/**
* @name Utilities (RedlichKwongMFTP)
*/
//@{
/**
* Equation of state type flag. The base class returns
* zero. Subclasses should define this to return a unique
@ -132,7 +107,8 @@ public:
*/
virtual int eosType() const;
//@}
//! @name Molar Thermodynamic properties
//! @{
/// Molar enthalpy. Units: J/kmol.
virtual doublereal enthalpy_mole() const;
@ -152,11 +128,9 @@ public:
/// Molar heat capacity at constant volume. Units: J/kmol/K.
virtual doublereal cv_mole() const;
/**
* @}
* @name Mechanical Properties
* @{
*/
//! @}
//! @name Mechanical Properties
//! @{
//! Return the thermodynamic pressure (Pa).
/*!
@ -167,7 +141,6 @@ public:
* \f[
* P = \frac{RT}{v-b_{mix}} - \frac{a_{mix}}{T^{0.5} v \left( v + b_{mix} \right) }
* \f]
*
*/
virtual doublereal pressure() const;
@ -179,6 +152,7 @@ public:
* \f]
*/
virtual doublereal isothermalCompressibility() const;
// @}
protected:
/**
@ -200,9 +174,6 @@ protected:
* species standard state molar volumes.
* The species molar volumes may be functions
* of temperature and pressure.
*
* NOTE: This is a non-virtual function, which is not a
* member of the ThermoPhase base class.
*/
virtual void calcDensity();
@ -213,8 +184,6 @@ protected:
* function sets the temperature, and makes sure that
* the value propagates to underlying objects
*
* @todo Make Phase::setTemperature a virtual function
*
* @param temp Temperature in kelvin
*/
virtual void setTemperature(const doublereal temp);
@ -234,8 +203,7 @@ protected:
* by a constraint equation as part of a larger set of
* equations.
*
* @param y Input vector of mass fractions.
* Length is m_kk.
* @param y Input vector of mass fractions. Length is m_kk.
*/
virtual void setMassFractions_NoNorm(const doublereal* const y);
@ -253,12 +221,10 @@ protected:
* condition is being handled by some other means, for example
* by a constraint equation as part of a larger set ofequations.
*
* @param x Input vector of mole fractions.
* Length is m_kk.
* @param x Input vector of mole fractions. Length is m_kk.
*/
virtual void setMoleFractions_NoNorm(const doublereal* const x);
//! Set the concentrations to the specified values within the phase.
/*!
* @param c The input vector to this routine is in dimensional
@ -270,9 +236,7 @@ protected:
*/
virtual void setConcentrations(const doublereal* const c);
public:
//! This method returns an array of generalized concentrations
/*!
* \f$ C^a_k\f$ are defined such that \f$ a_k = C^a_k /
@ -330,14 +294,16 @@ public:
* Inherited classes are responsible for overriding the default
* values if necessary.
*
* @param uA Output vector containing the units
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
* @param uA Output vector containing the units:
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*
* @param k species index. Defaults to 0.
* @param sizeUA output int containing the size of the vector.
* Currently, this is equal to 6.
@ -355,8 +321,7 @@ public:
*/
virtual void getActivityCoefficients(doublereal* ac) const;
/// @name Partial Molar Properties of the Solution (RedlichKwongMFTP)
/// @name Partial Molar Properties of the Solution
//@{
//! Get the array of non-dimensional species chemical potentials.
@ -422,35 +387,9 @@ public:
virtual void getPartialMolarVolumes(doublereal* vbar) const;
//@}
/*!
* @name Properties of the Standard State of the Species in the Solution
*
* Properties of the standard states are delegated to the VPSSMgr object.
* The values are cached within this object, and are not recalculated unless
* the temperature or pressure changes.
*/
//@{
//@}
/// @name Thermodynamic Values for the Species Reference States (RedlichKwongMFTP)
/*!
* Properties of the reference states are delegated to the VPSSMgr object.
* The values are cached within this object, and are not recalculated unless
* the temperature or pressure changes.
*/
//@{
//@}
//---------------------------------------------------------
/// @name Critical State Properties.
/// These methods are only implemented by some subclasses, and may
/// be moved out of ThermoPhase at a later date.
//@{
/// Critical temperature (K).
@ -461,14 +400,11 @@ public:
/// Critical density (kg/m3).
virtual doublereal critDensity() const;
//@}
public:
//! @name Initialization Methods - For Internal use (VPStandardState)
//@}
//! @name Initialization Methods - For Internal use
/*!
* The following methods are used in the process of constructing
* the phase and setting its parameters from a specification in an
@ -478,9 +414,7 @@ public:
*/
//@{
//! Set equation of state parameter values from XML
//! entries.
//! Set equation of state parameter values from XML entries.
/*!
* This method is called by function importPhase in
* file importCTML.cpp when processing a phase definition in
@ -509,7 +443,6 @@ public:
*/
virtual void initThermo();
//!This method is used by the ChemEquil equilibrium solver.
/*!
* It sets the state such that the chemical potentials satisfy
@ -560,32 +493,25 @@ private:
*/
void readXMLPureFluid(XML_Node& pureFluidParam);
//! Apply mixing rules for a coefficients
void applyStandardMixingRules();
//! Read the cross species RedlichKwong input parameters
/*!
* @param pureFluidParam XML_Node for the cross fluid parameters
*/
void readXMLCrossFluid(XML_Node& pureFluidParam);
//==============================================================================
private:
//! @internal Initialize the internal lengths in this object.
/*!
* Note this is not a virtual function and only handles
* this object
* Note this is not a virtual function and only handles this object
*/
void initLengths();
//==============================================================================
// Special functions inherited from MixtureFugacityTP
// @}
protected:
// Special functions inherited from MixtureFugacityTP
//! Calculate the deviation terms for the total entropy of the mixture from the
//! ideal gas mixture
@ -647,7 +573,6 @@ public:
*/
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phase, doublereal rhoguess);
public:
//! Return the value of the density at the liquid spinodal point (on the liquid side)
//! for the current temperature.
/*!
@ -655,7 +580,6 @@ public:
*/
virtual doublereal densSpinodalLiquid() const;
//! Return the value of the density at the gas spinodal point (on the gas side)
//! for the current temperature.
/*!
@ -663,8 +587,6 @@ public:
*/
virtual doublereal densSpinodalGas() const;
//! Calculate the pressure given the temperature and the molar volume
/*!
* Calculate the pressure given the temperature and the molar volume
@ -676,7 +598,6 @@ public:
*/
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const;
//! Calculate the pressure and the pressure derivative given the temperature and the molar volume
/*!
* Temperature and mole number are held constant
@ -690,17 +611,14 @@ public:
*/
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const;
//! Calculate dpdV and dpdT at the current conditions
/*!
* These are stored internally.
*/
void pressureDerivatives() const;
virtual void updateMixingExpressions();
//! Update the a and b parameters
/*!
* The a and the b parameters depend on the mole fraction and the temperature.
@ -708,10 +626,8 @@ public:
*/
void updateAB();
//! Calculate the a and the b parameters given the temperature
/*!
*
* This function doesn't change the internal state of the object, so it is a const
* function. It does use the stored mole fractions in the object.
*
@ -722,24 +638,27 @@ public:
*/
void calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const;
//=========================================================================================
// Special functions not inherited from MixtureFugacityTP
// Special functions not inherited from MixtureFugacityTP
doublereal da_dt() const;
void calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
doublereal& pc, doublereal& tc, doublereal& vc) const;
//! Solve the cubic equation of state
/*!
* The R-K equation of state may be solved via the following formula:
*
* V**3 - V**2(RT/P) - V(RTb/P - a/(P T**.5) + b*b) - (a b / (P T**.5)) = 0
*
* Returns the number of solutions found. If it only finds the liquid
* branch solution, it will return a -1 or a -2 instead of 1 or 2. If it
* returns 0, then there is an error.
*/
int NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
doublereal Vroot[3]) const;
//@}
//==============================================================================
protected:
//! boolean indicating whether standard mixing rules are applied
/*!
* - 1 = Yes, there are standard cross terms in the a coefficient matrices.
@ -754,7 +673,6 @@ protected:
*/
int m_formTempParam;
//! Value of b in the equation of state
/*!
* m_b is a function of the temperature and the mole fraction.
@ -767,13 +685,11 @@ protected:
*/
doublereal m_a_current;
vector_fp a_vec_Curr_;
vector_fp b_vec_Curr_;
Array2D a_coeff_vec;
vector_fp m_pc_Species;
vector_fp m_tc_Species;
vector_fp m_vc_Species;
@ -782,8 +698,6 @@ protected:
doublereal Vroot_[3];
//! Temporary storage - length = m_kk.
mutable vector_fp m_pp;
@ -795,8 +709,6 @@ protected:
// Partial molar volumes of the species
mutable vector_fp m_partialMolarVolumes;
//! The derivative of the pressure wrt the volume
/*!
* Calculated at the current conditions
@ -830,8 +742,6 @@ public:
//! Omega constant for the critical molar volume
static const doublereal omega_vc;
};
}

220
src/thermo/IdealGasPhase.cpp
View file

@ -14,7 +14,7 @@ using namespace std;
namespace Cantera
{
// Default empty Constructor
IdealGasPhase::IdealGasPhase() :
m_p0(-1.0),
m_tlast(0.0),
@ -42,7 +42,6 @@ IdealGasPhase::~IdealGasPhase()
{
}
// Copy Constructor
IdealGasPhase::IdealGasPhase(const IdealGasPhase& right) :
m_p0(right.m_p0),
m_tlast(right.m_tlast),
@ -55,14 +54,6 @@ IdealGasPhase::IdealGasPhase(const IdealGasPhase& right) :
*this = right;
}
// Assignment operator
/*
* Assignment operator for the object. Constructed
* object will be a clone of this object, but will
* also own all of its data.
*
* @param right Object to be copied.
*/
IdealGasPhase& IdealGasPhase::operator=(const IdealGasPhase& right)
{
if (&right != this) {
@ -80,14 +71,6 @@ IdealGasPhase& IdealGasPhase::operator=(const IdealGasPhase& right)
return *this;
}
// Duplicator from the %ThermoPhase parent class
/*
* Given a pointer to a %ThermoPhase object, this function will
* duplicate the %ThermoPhase object and all underlying structures.
* This is basically a wrapper around the copy constructor.
*
* @return returns a pointer to a %ThermoPhase
*/
ThermoPhase* IdealGasPhase::duplMyselfAsThermoPhase() const
{
return new IdealGasPhase(*this);
@ -95,80 +78,31 @@ ThermoPhase* IdealGasPhase::duplMyselfAsThermoPhase() const
// Molar Thermodynamic Properties of the Solution ------------------
/*
* Molar internal energy. J/kmol. For an ideal gas mixture,
* \f[
* \hat u(T) = \sum_k X_k \hat h^0_k(T) - \hat R T,
* \f]
* and is a function only of temperature.
* The reference-state pure-species enthalpies
* \f$ \hat h^0_k(T) \f$ are computed by the species thermodynamic
* property manager.
* @see SpeciesThermo
*/
doublereal IdealGasPhase::intEnergy_mole() const
{
return GasConstant * temperature() * (mean_X(&enthalpy_RT_ref()[0]) - 1.0);
}
/*
* Molar entropy. Units: J/kmol/K.
* For an ideal gas mixture,
* \f[
* \hat s(T, P) = \sum_k X_k \hat s^0_k(T) - \hat R \log (P/P^0).
* \f]
* The reference-state pure-species entropies
* \f$ \hat s^0_k(T) \f$ are computed by the species thermodynamic
* property manager.
* @see SpeciesThermo
*/
doublereal IdealGasPhase::entropy_mole() const
{
return GasConstant * (mean_X(&entropy_R_ref()[0]) - sum_xlogx() - std::log(pressure() / m_spthermo->refPressure()));
}
/*
* Molar Gibbs free Energy for an ideal gas.
* Units = J/kmol.
*/
doublereal IdealGasPhase::gibbs_mole() const
{
return enthalpy_mole() - temperature() * entropy_mole();
}
/*
* Molar heat capacity at constant pressure. Units: J/kmol/K.
* For an ideal gas mixture,
* \f[
* \hat c_p(t) = \sum_k \hat c^0_{p,k}(T).
* \f]
* The reference-state pure-species heat capacities
* \f$ \hat c^0_{p,k}(T) \f$ are computed by the species thermodynamic
* property manager.
* @see SpeciesThermo
*/
doublereal IdealGasPhase::cp_mole() const
{
return GasConstant * mean_X(&cp_R_ref()[0]);
}
/*
* Molar heat capacity at constant volume. Units: J/kmol/K.
* For an ideal gas mixture,
* \f[ \hat c_v = \hat c_p - \hat R. \f]
*/
doublereal IdealGasPhase::cv_mole() const
{
return cp_mole() - GasConstant;
}
/**
* @returns species translational/rotational specific heat at
* constant volume.
*
* Either: $5/2 R_s$ or $3/2 R_s$ for molecules/atoms.
*
*/
doublereal IdealGasPhase::cv_tr(doublereal atomicity) const
{
// k is the species number
@ -189,44 +123,16 @@ doublereal IdealGasPhase::cv_tr(doublereal atomicity) const
return c[3];
}
/**
* @returns species translational specific heat at constant volume.
*/
doublereal IdealGasPhase::cv_trans() const
{
return 1.5 * GasConstant;
}
/**
* @returns species rotational specific heat at constant volume.
*
*/
doublereal IdealGasPhase::cv_rot(double atom) const
{
return std::max(cv_tr(atom) - cv_trans(), 0.);
}
/**
* @returns species vibrational specific heat at
* constant volume.
*
* C^{vib}_{v,s} = \frac{\partial e^{vib}_{v,s} }{\partial T}
*
* The species vibration energy ($e^{vib}_{v,s}$) is:
*
* 0: atom
*
* Diatomic:
* \f[
* \frac{R_s \theta_{v,s}}{e^{\theta_{v,s}/T}-1}
* \f]
*
* General Molecules:
* \f[
* \sum_i \frac{R_s \theta_{v,s,i}}{e^{\theta_{v,s,i}/T}-1}
* \f]
*
*/
doublereal IdealGasPhase::cv_vib(const int k, const doublereal T) const
{
@ -253,23 +159,12 @@ doublereal IdealGasPhase::cv_vib(const int k, const doublereal T) const
}
// Mechanical Equation of State ----------------------------
// Chemical Potentials and Activities ----------------------
/*
* Returns the standard concentration \f$ C^0_k \f$, which is used to normalize
* the generalized concentration.
*/
doublereal IdealGasPhase::standardConcentration(size_t k) const
{
double p = pressure();
return p / (GasConstant * temperature());
}
/*
* Returns the natural logarithm of the standard
* concentration of the kth species
*/
doublereal IdealGasPhase::logStandardConc(size_t k) const
{
_updateThermo();
@ -278,9 +173,6 @@ doublereal IdealGasPhase::logStandardConc(size_t k) const
return lc;
}
/*
* Get the array of non-dimensional activity coefficients
*/
void IdealGasPhase::getActivityCoefficients(doublereal* ac) const
{
for (size_t k = 0; k < m_kk; k++) {
@ -288,10 +180,6 @@ void IdealGasPhase::getActivityCoefficients(doublereal* ac) const
}
}
/*
* Get the array of chemical potentials at unit activity \f$
* \mu^0_k(T,P) \f$.
*/
void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
@ -318,10 +206,6 @@ void IdealGasPhase::getChemPotentials(doublereal* mu) const
}
}
/*
* Get the array of partial molar enthalpies of the species
* units = J / kmol
*/
void IdealGasPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -329,10 +213,6 @@ void IdealGasPhase::getPartialMolarEnthalpies(doublereal* hbar) const
scale(_h.begin(), _h.end(), hbar, rt);
}
/*
* Get the array of partial molar entropies of the species
* units = J / kmol / K
*/
void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const
{
const vector_fp& _s = entropy_R_ref();
@ -345,10 +225,6 @@ void IdealGasPhase::getPartialMolarEntropies(doublereal* sbar) const
}
}
/*
* Get the array of partial molar internal energies of the species
* units = J / kmol
*/
void IdealGasPhase::getPartialMolarIntEnergies(doublereal* ubar) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -358,19 +234,12 @@ void IdealGasPhase::getPartialMolarIntEnergies(doublereal* ubar) const
}
}
/*
* Get the array of partial molar heat capacities
*/
void IdealGasPhase::getPartialMolarCp(doublereal* cpbar) const
{
const vector_fp& _cp = cp_R_ref();
scale(_cp.begin(), _cp.end(), cpbar, GasConstant);
}
/*
* Get the array of partial molar volumes
* units = m^3 / kmol
*/
void IdealGasPhase::getPartialMolarVolumes(doublereal* vbar) const
{
double vol = 1.0 / molarDensity();
@ -381,22 +250,12 @@ void IdealGasPhase::getPartialMolarVolumes(doublereal* vbar) const
// Properties of the Standard State of the Species in the Solution --
/*
* Get the nondimensional Enthalpy functions for the species
* at their standard states at the current T and P of the
* solution
*/
void IdealGasPhase::getEnthalpy_RT(doublereal* hrt) const
{
const vector_fp& _h = enthalpy_RT_ref();
copy(_h.begin(), _h.end(), hrt);
}
/*
* Get the array of nondimensional entropy functions for the
* standard state species
* at the current <I>T</I> and <I>P</I> of the solution.
*/
void IdealGasPhase::getEntropy_R(doublereal* sr) const
{
const vector_fp& _s = entropy_R_ref();
@ -407,10 +266,6 @@ void IdealGasPhase::getEntropy_R(doublereal* sr) const
}
}
/*
* Get the nondimensional gibbs function for the species
* standard states at the current T and P of the solution.
*/
void IdealGasPhase::getGibbs_RT(doublereal* grt) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
@ -421,11 +276,6 @@ void IdealGasPhase::getGibbs_RT(doublereal* grt) const
}
}
/*
* get the pure Gibbs free energies of each species assuming
* it is in its standard state. This is the same as
* getStandardChemPotentials().
*/
void IdealGasPhase::getPureGibbs(doublereal* gpure) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
@ -437,11 +287,6 @@ void IdealGasPhase::getPureGibbs(doublereal* gpure) const
}
}
/*
* Returns the vector of nondimensional
* internal Energies of the standard state at the current temperature
* and pressure of the solution for each species.
*/
void IdealGasPhase::getIntEnergy_RT(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -450,25 +295,12 @@ void IdealGasPhase::getIntEnergy_RT(doublereal* urt) const
}
}
/*
* Get the nondimensional heat capacity at constant pressure
* function for the species
* standard states at the current T and P of the solution.
*/
void IdealGasPhase::getCp_R(doublereal* cpr) const
{
const vector_fp& _cpr = cp_R_ref();
copy(_cpr.begin(), _cpr.end(), cpr);
}
/*
* Get the molar volumes of the species standard states at the current
* <I>T</I> and <I>P</I> of the solution.
* units = m^3 / kmol
*
* @param vol Output vector containing the standard state volumes.
* Length: m_kk.
*/
void IdealGasPhase::getStandardVolumes(doublereal* vol) const
{
double tmp = 1.0 / molarDensity();
@ -479,56 +311,30 @@ void IdealGasPhase::getStandardVolumes(doublereal* vol) const
// Thermodynamic Values for the Species Reference States ---------
/*
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* and reference pressure.
*/
void IdealGasPhase::getEnthalpy_RT_ref(doublereal* hrt) const
{
const vector_fp& _h = enthalpy_RT_ref();
copy(_h.begin(), _h.end(), hrt);
}
/*
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* and reference pressure.
*/
void IdealGasPhase::getGibbs_RT_ref(doublereal* grt) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
copy(gibbsrt.begin(), gibbsrt.end(), grt);
}
/*
* Returns the vector of the
* gibbs function of the reference state at the current temperature
* and reference pressure.
* units = J/kmol
*/
void IdealGasPhase::getGibbs_ref(doublereal* g) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), g, _RT());
}
/*
* Returns the vector of nondimensional
* entropies of the reference state at the current temperature
* and reference pressure.
*/
void IdealGasPhase::getEntropy_R_ref(doublereal* er) const
{
const vector_fp& _s = entropy_R_ref();
copy(_s.begin(), _s.end(), er);
}
/*
* Returns the vector of nondimensional
* internal Energies of the reference state at the current temperature
* of the solution and the reference pressure for each species.
*/
void IdealGasPhase::getIntEnergy_RT_ref(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -537,11 +343,6 @@ void IdealGasPhase::getIntEnergy_RT_ref(doublereal* urt) const
}
}
/*
* Returns the vector of nondimensional
* constant pressure heat capacities of the reference state
* at the current temperature and reference pressure.
*/
void IdealGasPhase::getCp_R_ref(doublereal* cprt) const
{
const vector_fp& _cpr = cp_R_ref();
@ -556,8 +357,6 @@ void IdealGasPhase::getStandardVolumes_ref(doublereal* vol) const
}
}
// new methods defined here -------------------------------
void IdealGasPhase::initThermo()
{
m_p0 = refPressure();
@ -569,11 +368,6 @@ void IdealGasPhase::initThermo()
m_pp.resize(m_kk);
}
/*
* Set mixture to an equilibrium state consistent with specified
* chemical potentials and temperature. This method is needed by
* the ChemEquil equilibrium solver.
*/
void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
{
double tmp, tmp2;
@ -605,17 +399,6 @@ void IdealGasPhase::setToEquilState(const doublereal* mu_RT)
setState_PX(pres, &m_pp[0]);
}
/// This method is called each time a thermodynamic property is
/// requested, to check whether the internal species properties
/// within the object need to be updated.
/// Currently, this updates the species thermo polynomial values
/// for the current value of the temperature. A check is made
/// to see if the temperature has changed since the last
/// evaluation. This object does not contain any persistent
/// data that depends on the concentration, that needs to be
/// updated. The state object modifies its concentration
/// dependent information at the time the setMoleFractions()
/// (or equivalent) call is made.
void IdealGasPhase::_updateThermo() const
{
doublereal tnow = temperature();
@ -635,4 +418,3 @@ void IdealGasPhase::_updateThermo() const
}
}
}

592
src/thermo/IdealSolidSolnPhase.cpp
View file

@ -21,11 +21,6 @@ using namespace std;
namespace Cantera
{
/*
* Constructor for IdealSolidSolnPhase class:
* The default form for the generalized concentrations is 0
* i.e., unity.
*/
IdealSolidSolnPhase::IdealSolidSolnPhase(int formGC) :
ThermoPhase(),
m_formGC(formGC),
@ -53,7 +48,7 @@ IdealSolidSolnPhase::IdealSolidSolnPhase(const std::string& inputFile,
}
initThermoFile(inputFile, id);
}
//====================================================================================================================
IdealSolidSolnPhase::IdealSolidSolnPhase(XML_Node& root, const std::string& id,
int formGC) :
ThermoPhase(),
@ -68,12 +63,11 @@ IdealSolidSolnPhase::IdealSolidSolnPhase(XML_Node& root, const std::string& id,
}
importPhase(*findXMLPhase(&root, id), this);
}
//====================================================================================================================
IdealSolidSolnPhase::IdealSolidSolnPhase(const IdealSolidSolnPhase& b)
{
*this = b;
}
//====================================================================================================================
IdealSolidSolnPhase& IdealSolidSolnPhase::
operator=(const IdealSolidSolnPhase& b)
@ -97,22 +91,11 @@ operator=(const IdealSolidSolnPhase& b)
return *this;
}
/*
* Base Class Duplication Function
* -> given a pointer to ThermoPhase, this function can
* duplicate the object. (note has to be a separate function
* not the copy constructor, because it has to be
* a virtual function)
*/
ThermoPhase* IdealSolidSolnPhase::duplMyselfAsThermoPhase() const
{
return new IdealSolidSolnPhase(*this);
}
//====================================================================================================================
/**
* Equation of state flag. Returns the value cIdealGas, defined
* in mix_defs.h.
*/
int IdealSolidSolnPhase::eosType() const
{
integer res;
@ -136,19 +119,7 @@ int IdealSolidSolnPhase::eosType() const
/********************************************************************
* Molar Thermodynamic Properties of the Solution
********************************************************************/
/**
* Molar enthalpy of the solution. Units: J/kmol.
* For an ideal, constant partial molar volume solution mixture with
* pure species phases which exhibit zero volume expansivity and
* zero isothermal compressibility:
* \f[
* \hat h(T,P) = \sum_k X_k \hat h^0_k(T) + (P - P_{ref}) (\sum_k X_k \hat V^0_k)
* \f]
* The reference-state pure-species enthalpies at the reference pressure Pref
* \f$ \hat h^0_k(T) \f$, are computed by the species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*/
doublereal IdealSolidSolnPhase::
enthalpy_mole() const
{
@ -157,21 +128,6 @@ enthalpy_mole() const
return (htp + (pressure() - m_Pref)/molarDensity());
}
/**
* Molar internal energy of the solution. J/kmol.
* For an ideal, constant partial molar volume solution mixture with
* pure species phases which exhibit zero volume expansivity and
* zero isothermal compressibility:
* \f[
* \hat u(T) = \hat h(T,P) - p \hat V = \sum_k X_k \hat h^0_k(T)
* - P_{ref} (\sum_k X_k \hat V^0_k)
* \f]
* and is a function only of temperature.
* The reference-state pure-species enthalpies
* \f$ \hat h^0_k(T) \f$ are computed by the species thermodynamic
* property manager.
* @see SpeciesThermo
*/
doublereal IdealSolidSolnPhase::intEnergy_mole() const
{
const double* eptr = DATA_PTR(enthalpy_RT_ref().begin());
@ -180,39 +136,12 @@ doublereal IdealSolidSolnPhase::intEnergy_mole() const
return (htp - m_Pref / molarDensity());
}
/**
* Molar entropy of the solution. Units: J/kmol/K.
* For an ideal, constant partial molar volume solution mixture with
* pure species phases which exhibit zero volume expansivity:
* \f[
* \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T)
* - \hat R \sum_k X_k log(X_k)
* \f]
* The reference-state pure-species entropies
* \f$ \hat s^0_k(T,p_{ref}) \f$ are computed by the species thermodynamic
* property manager. The pure species entropies are independent of
* temperature since the volume expansivities are equal to zero.
* @see SpeciesThermo
*/
doublereal IdealSolidSolnPhase::entropy_mole() const
{
const double* dptr = DATA_PTR(entropy_R_ref());
return GasConstant * (mean_X(dptr) - sum_xlogx());
}
/**
* Molar gibbs free energy of the solution. Units: J/kmol.
* For an ideal, constant partial molar volume solution mixture with
* pure species phases which exhibit zero volume expansivity:
* \f[
* \hat g(T, P) = \sum_k X_k \hat g^0_k(T,P) + \hat R T \sum_k X_k log(X_k)
* \f]
* The reference-state pure-species gibbs free energies
* \f$ \hat g^0_k(T) \f$ are computed by the species thermodynamic
* property manager, while the standard state gibbs free energies
* \f$ \hat g^0_k(T,P) \f$ are computed by the member function, gibbs_RT().
* @see SpeciesThermo
*/
doublereal IdealSolidSolnPhase::gibbs_mole() const
{
const double* dptr = DATA_PTR(gibbs_RT_ref());
@ -220,20 +149,6 @@ doublereal IdealSolidSolnPhase::gibbs_mole() const
return (GasConstant * temperature() * (g + sum_xlogx()));
}
/**
* Molar heat capacity at constant pressure of the solution.
* Units: J/kmol/K.
* For an ideal, constant partial molar volume solution mixture with
* pure species phases which exhibit zero volume expansivity:
* \f[
* \hat c_p(T,P) = \sum_k X_k \hat c^0_{p,k}(T) .
* \f]
* The heat capacity is independent of pressure.
* The reference-state pure-species heat capacities
* \f$ \hat c^0_{p,k}(T) \f$ are computed by the species thermodynamic
* property manager.
* @see SpeciesThermo
*/
doublereal IdealSolidSolnPhase::cp_mole() const
{
const double* dptr = DATA_PTR(cp_R_ref());
@ -243,27 +158,7 @@ doublereal IdealSolidSolnPhase::cp_mole() const
/********************************************************************
* Mechanical Equation of State
********************************************************************/
/**
*
* Calculate the density of the mixture using the partial
* molar volumes and mole fractions as input
*
* The formula for this is
*
* \f[
* \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}}
* \f]
*
* where \f$ X_k \f$ are the mole fractions, \f$W_k\f$ are
* the molecular weights, and \f$V_k\f$ are the pure species
* molar volumes.
*
* Note, the basis behind this formula is that in an ideal
* solution the partial molar volumes are equal to the pure
* species molar volumes. We have additionally specified that
* in this class that the pure species molar volumes are
* independent of temperature and pressure.
*/
void IdealSolidSolnPhase::calcDensity()
{
/*
@ -280,22 +175,6 @@ void IdealSolidSolnPhase::calcDensity()
Phase::setDensity(dens);
}
/**
* Overwritten setDensity() function is necessary because the
* density is not an independent variable.
*
* This function will now throw an error condition
*
* @internal May have to adjust the strategy here to make
* the eos for these materials slightly compressible, in order
* to create a condition where the density is a function of
* the pressure.
*
* This function will now throw an error condition.
*
* NOTE: This is a virtual function that overwrites the State.h
* class
*/
void IdealSolidSolnPhase::
setDensity(const doublereal rho)
{
@ -312,90 +191,42 @@ setDensity(const doublereal rho)
}
}
/*
* setPressure(double) (virtual from ThermoPhase)
*
* Set the pressure at constant temperature. Units: Pa.
* This method sets a constant within the object.
* The mass density is not a function of pressure.
* Note: This function overrides the setPressure() function
* in the ThermoPhase object.
* We calculate the density and store it in the
* State object, because this density is supposed to
* be current after setting the pressure, and is now
* a dependent variable.
*/
void IdealSolidSolnPhase::setPressure(doublereal p)
{
m_Pcurrent = p;
calcDensity();
}
/*
* setMolarDensity() (virtual from State)
* Overwritten setMolarDensity() function is necessary because the
* density is not an independent variable.
*
* This function will now throw an error condition.
*
* NOTE: This is a virtual function that overrides the State.h
* class
*/
void IdealSolidSolnPhase::setMolarDensity(const doublereal n)
{
throw CanteraError("IdealSolidSolnPhase::setMolarDensity",
"Density is not an independent variable");
}
/*
* setMoleFractions() (virtual from State)
*
* Sets the mole fractions and adjusts the internal density.
*/
void IdealSolidSolnPhase::setMoleFractions(const doublereal* const x)
{
Phase::setMoleFractions(x);
calcDensity();
}
/**
* setMoleFractions_NoNorm() (virtual from State)
*
* Sets the mole fractions and adjusts the internal density.
*/
void IdealSolidSolnPhase::setMoleFractions_NoNorm(const doublereal* const x)
{
Phase::setMoleFractions(x);
calcDensity();
}
/*
* setMassFractions() (virtual from State)
*
* Sets the mass fractions and adjusts the internal density.
*/
void IdealSolidSolnPhase::setMassFractions(const doublereal* const y)
{
Phase::setMassFractions(y);
calcDensity();
}
/*
* setMassFractions_NoNorm() (virtual from State)
*
* Sets the mass fractions and adjusts the internal density.
*/
void IdealSolidSolnPhase::setMassFractions_NoNorm(const doublereal* const y)
{
Phase::setMassFractions_NoNorm(y);
calcDensity();
}
/*
* setConcentrations (virtual from State)
*
* Sets the concentrations and adjusts the internal density
*/
void IdealSolidSolnPhase::setConcentrations(const doublereal* const c)
{
Phase::setConcentrations(c);
@ -406,54 +237,6 @@ void IdealSolidSolnPhase::setConcentrations(const doublereal* const c)
* Chemical Potentials and Activities
********************************************************************/
/********************************************************************
*
* getActivitConcentrations():
*
* This method returns the array of generalized
* concentrations. The generalized concentrations are used
* in the evaluation of the rates of progress for reactions
* involving species in this phase. The generalized
* concentration divided by the standard concentration is also
* equal to the activity of species.
*
* For this implementation the activity is defined to be the
* mole fraction of the species. The generalized concentration
* is defined to be equal to the mole fraction divided by
* the partial molar volume. The generalized concentrations
* for species in this phase therefore have units of
* kmol m<SUP>-3</SUP>. Rate constants must reflect this fact.
*
* On a general note, the following must be true.
* For an ideal solution, the generalized concentration must consist
* of the mole fraction multiplied by a constant. The constant may be
* fairly arbitrarily chosen, with differences adsorbed into the
* reaction rate expression. 1/V_N, 1/V_k, or 1 are equally good,
* as long as the standard concentration is adjusted accordingly.
* However, it must be a constant (and not the concentration, btw,
* which is a function of the mole fractions) in order for the
* ideal solution properties to hold at the same time having the
* standard concentration to be independent of the mole fractions.
*
* In this implementation the form of the generalized concentrations
* depend upon the member attribute, m_formGC:
*
* <TABLE>
* <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
* <TR><TD> 0 </TD><TD> X_k </TD><TD> 1.0 </TD></TR>
* <TR><TD> 1 </TD><TD> X_k / V_k </TD><TD> 1.0 / V_k </TD></TR>
* <TR><TD> 2 </TD><TD> X_k / V_N </TD><TD> 1.0 / V_N </TD></TR>
* </TABLE>
*
* HKM Note: We have absorbed the pressure dependence of the pure species
* state into the thermodynamics functions. Therefore the
* standard state on which the activities are based depend
* on both temperature and pressure. If we hadn't, it would have
* appeared in this function in a very awkward exp[] format.
*
* @param c[] Pointer to array of doubles of length m_kk, which on exit
* will contain the generalized concentrations.
*/
void IdealSolidSolnPhase::
getActivityConcentrations(doublereal* c) const
{
@ -479,23 +262,6 @@ getActivityConcentrations(doublereal* c) const
}
}
/*********************************************************************
*
* standardConcentration()
*
* The standard concentration \f$ C^0_k \f$ used to normalize
* the generalized concentration.
* In many cases, this quantity
* will be the same for all species in a phase.
* However, for this case, we will return a distinct concentration
* for each species. This is the inverse of the species molar
* volume. Units are m<SUP>3</SUP> kmol<SUP>-1</SUP>.
*
*
* @param k Species number: this is a require parameter,
* a change from the ThermoPhase base class, where it was
* an optional parameter.
*/
doublereal IdealSolidSolnPhase::
standardConcentration(size_t k) const
{
@ -523,16 +289,6 @@ referenceConcentration(int k) const
return 0.0;
}
/*********************************************************************
*
* logStandardConc()
*
* Returns the log of the standard concentration
*
* @param k Species number: this is a require parameter,
* a change from the ThermoPhase base class, where it was
* an optional parameter.
*/
doublereal IdealSolidSolnPhase::
logStandardConc(size_t k) const
{
@ -555,32 +311,6 @@ logStandardConc(size_t k) const
return res;
}
/***********************************************************************
*
* getUnitsStandardConcentration()
*
* Returns the units of the standard and general concentrations
* Note they have the same units, as their divisor is
* defined to be equal to the activity of the kth species
* in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*
* For EOS types other than cIdealSolidSolnPhase1, the default
* kmol/m3 holds for standard concentration units. For
* cIdealSolidSolnPhase0 type, the standard concentration is
* unitless.
*/
void IdealSolidSolnPhase::
getUnitsStandardConc(double* uA, int, int sizeUA) const
{
@ -613,10 +343,6 @@ getUnitsStandardConc(double* uA, int, int sizeUA) const
}
}
/*
* getActivityCoefficients():
*
*/
void IdealSolidSolnPhase::
getActivityCoefficients(doublereal* ac) const
{
@ -624,23 +350,7 @@ getActivityCoefficients(doublereal* ac) const
ac[k] = 1.0;
}
}
//================================================================================================
/*
*
* getChemPotentials():
*
* This function returns a vector of chemical potentials of the
* species.
* \f[
* \mu_k = \mu^o_k(T) + V_k * (p - p_o) + R T ln(X_k)
* \f]
* or another way to phrase this is
* \f[
* \mu_k = \mu^o_k(T,p) + R T ln(X_k)
* \f]
* where \f$ \mu^o_k(T,p) = \mu^o_k(T) + V_k * (p - p_o)\f$
*
*/
void IdealSolidSolnPhase::
getChemPotentials(doublereal* mu) const
{
@ -654,24 +364,7 @@ getChemPotentials(doublereal* mu) const
+ delta_p * m_speciesMolarVolume[k];
}
}
//================================================================================================
/*
*
* getChemPotentials_RT()
*
* Get the array of non-dimensional chemical potentials \f$
* \mu_k / \hat R T \f$, where
*
* \f[
* \mu_k = \mu^o_k(T) + V_k * (p - p_o) + R T ln(X_k)
* \f]
* or another way to phrase this is
* \f[
* \mu_k = \mu^o_k(T,p) + R T ln(X_k)
* \f]
* where \f$ \mu^o_k(T,p) = \mu^o_k(T) + V_k * (p - p_o)\f$
*
*/
void IdealSolidSolnPhase::
getChemPotentials_RT(doublereal* mu) const
{
@ -690,21 +383,6 @@ getChemPotentials_RT(doublereal* mu) const
* Partial Molar Properties
********************************************************************/
/********************************************************************
*
* getPartialMolarEnthalpies()
*
* For this phase, the partial molar enthalpies are equal to the
* pure species enthalpies.
* \f[
* \hat h_k(T,P) = \sum_k X_k \hat h^0_k(T) + (p - p_{ref}) (\sum_k X_k \hat V^0_k)
* \f]
* The reference-state pure-species enthalpies at the reference
* pressure p_ref
* \f$ \hat h^0_k(T) \f$, are computed by the species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*/
void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -712,23 +390,6 @@ void IdealSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const
scale(_h.begin(), _h.end(), hbar, rt);
}
/********************************************************************
*
* getPartialMolarEntropies()
*
* Returns an array of partial molar entropies of the species in the
* solution. Units: J/kmol.
* For this phase, the partial molar entropies are equal to the
* pure species entropies plus the ideal solution contribution.
* \f[
* \bar s_k(T,P) = \hat s^0_k(T) - R log(X_k)
* \f]
* The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$,
* at the reference pressure, \f$ P_{ref} \f$, are computed by the
* species thermodynamic
* property manager. They are polynomial functions of temperature.
* @see SpeciesThermo
*/
void IdealSolidSolnPhase::
getPartialMolarEntropies(doublereal* sbar) const
{
@ -741,14 +402,6 @@ getPartialMolarEntropies(doublereal* sbar) const
}
}
/********************************************************************
*
* getPartialMolarCp()
*
* For this phase, the partial molar heat capacities are equal
* to the standard state heat capacities.
*/
void IdealSolidSolnPhase::
getPartialMolarCp(doublereal* cpbar) const
{
@ -758,16 +411,6 @@ getPartialMolarCp(doublereal* cpbar) const
}
}
/******************************************************************
*
* getPartialMolarVolumes()
*
* returns an array of partial molar volumes of the species
* in the solution. Units: m^3 kmol-1.
*
* For this solution, thepartial molar volumes are equal to the
* constant species molar volumes.
*/
void IdealSolidSolnPhase::
getPartialMolarVolumes(doublereal* vbar) const
{
@ -775,25 +418,9 @@ getPartialMolarVolumes(doublereal* vbar) const
}
/*****************************************************************
* Properties of the Standard State of the Species
* in the Solution
* Properties of the Standard State of the Species in the Solution
*****************************************************************/
/******************************************************************
*
* getPureGibbs()
*
* Get the Gibbs functions for the pure species
* at the current <I>T</I> and <I>P</I> of the solution.
* We assume an incompressible constant partial molar
* volume here:
* \f[
* \mu^0_k(T,p) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k
* \f]
* where \f$V_k\f$ is the molar volume of pure species <I>k<\I>.
* \f$ u^{ref}_k(T)\f$ is the chemical potential of pure
* species <I>k<\I> at the reference pressure, \f$P_{ref}\f$.
*/
void IdealSolidSolnPhase::
getPureGibbs(doublereal* gpure) const
{
@ -806,21 +433,6 @@ getPureGibbs(doublereal* gpure) const
}
}
/**
* Get the nondimensional gibbs function for the species
* standard states at the current T and P of the solution.
*
* \f[
* \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k
* \f]
* where \f$V_k\f$ is the molar volume of pure species <I>k</I>.
* \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure
* species <I>k</I> at the reference pressure, \f$P_{ref}\f$.
*
* @param grt Vector of length m_kk, which on return sr[k]
* will contain the nondimensional
* standard state gibbs function for species k.
*/
void IdealSolidSolnPhase::
getGibbs_RT(doublereal* grt) const
{
@ -833,21 +445,6 @@ getGibbs_RT(doublereal* grt) const
}
}
/********************************************************************
*
* getEnthalpy_RT()
*
* Get the array of nondimensional Enthalpy functions for the ss
* species at the current <I>T</I> and <I>P</I> of the solution.
* We assume an incompressible constant partial molar
* volume here:
* \f[
* h^0_k(T,P) = h^{ref}_k(T) + (P - P_{ref}) * V_k
* \f]
* where \f$V_k\f$ is the molar volume of pure species <I>k<\I>.
* \f$ h^{ref}_k(T)\f$ is the enthalpy of the pure
* species <I>k<\I> at the reference pressure, \f$P_{ref}\f$.
*/
void IdealSolidSolnPhase::
getEnthalpy_RT(doublereal* hrt) const
{
@ -859,36 +456,12 @@ getEnthalpy_RT(doublereal* hrt) const
}
}
/**
* Get the nondimensional Entropies for the species
* standard states at the current T and P of the solution.
*
* Note, this is equal to the reference state entropies
* due to the zero volume expansivity:
* i.e., (dS/dp)_T = (dV/dT)_P = 0.0
*
* @param sr Vector of length m_kk, which on return sr[k]
* will contain the nondimensional
* standard state entropy of species k.
*/
void IdealSolidSolnPhase::getEntropy_R(doublereal* sr) const
{
const vector_fp& _s = entropy_R_ref();
copy(_s.begin(), _s.end(), sr);
}
/*
* Returns the vector of nondimensional
* internal Energies of the standard state at the current temperature
* of the solution and current pressure for each species.
* \f[
* u^0_k(T,P) = h^{ref}_k(T) - P_{ref} * V_k
* \f]
*
* The standard state internal energy is independent of
* pressure in this equation of state.
* (inherited from ThermoPhase.h)
*/
void IdealSolidSolnPhase::getIntEnergy_RT(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -898,50 +471,21 @@ void IdealSolidSolnPhase::getIntEnergy_RT(doublereal* urt) const
}
}
/*
* Get the nondimensional heat capacity at constant pressure
* function for the species
* standard states at the current T and P of the solution.
*
* \f[
* Cp^0_k(T,P) = Cp^{ref}_k(T)
* \f]
* where \f$V_k\f$ is the molar volume of pure species <I>k<\I>.
* \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity
* of species <I>k<\I> at the reference pressure, \f$P_{ref}\f$.
*
* @param cpr Vector of length m_kk, which on return cpr[k]
* will contain the nondimensional
* constant pressure heat capacity for species k.
*/
void IdealSolidSolnPhase::getCp_R(doublereal* cpr) const
{
const vector_fp& _cpr = cp_R_ref();
copy(_cpr.begin(), _cpr.end(), cpr);
}
/*
* Get the molar volumes of each species in their standard
* states at the current
* <I>T</I> and <I>P</I> of the solution.
* units = m^3 / kmol
*/
void IdealSolidSolnPhase::getStandardVolumes(doublereal* vol) const
{
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vol);
}
/*********************************************************************
* Thermodynamic Values for the Species Reference States
*********************************************************************/
/*
* Returns the vector of non-dimensional Enthalpy function
* of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Units = unitless
*/
void IdealSolidSolnPhase::getEnthalpy_RT_ref(doublereal* hrt) const
{
_updateThermo();
@ -950,12 +494,6 @@ void IdealSolidSolnPhase::getEnthalpy_RT_ref(doublereal* hrt) const
}
}
/*
* Returns the vector of non-dimensional Gibbs function
* of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Units = unitless
*/
void IdealSolidSolnPhase::getGibbs_RT_ref(doublereal* grt) const
{
_updateThermo();
@ -964,12 +502,6 @@ void IdealSolidSolnPhase::getGibbs_RT_ref(doublereal* grt) const
}
}
/*
* Returns the vector of Gibbs function
* of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Units = J / kmol
*/
void IdealSolidSolnPhase::getGibbs_ref(doublereal* g) const
{
_updateThermo();
@ -979,12 +511,6 @@ void IdealSolidSolnPhase::getGibbs_ref(doublereal* g) const
}
}
/*
* Returns the vector of nondimensional
* internal Energies of the standard state at the current temperature
* of the solution and current pressure for each species.
* (inherited from ThermoPhase.h)
*/
void IdealSolidSolnPhase::getIntEnergy_RT_ref(doublereal* urt) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -994,12 +520,6 @@ void IdealSolidSolnPhase::getIntEnergy_RT_ref(doublereal* urt) const
}
}
/*
* Returns the vector of non-dimensional Entropy function
* of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Units = unitless
*/
void IdealSolidSolnPhase::getEntropy_R_ref(doublereal* er) const
{
_updateThermo();
@ -1008,12 +528,6 @@ void IdealSolidSolnPhase::getEntropy_R_ref(doublereal* er) const
}
}
/*
* Returns the vector of non-dimensional Entropy function
* of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* Units = unitless
*/
void IdealSolidSolnPhase::getCp_R_ref(doublereal* cpr) const
{
_updateThermo();
@ -1022,26 +536,12 @@ void IdealSolidSolnPhase::getCp_R_ref(doublereal* cpr) const
}
}
/*
* Returns a reference to the vector of nondimensional
* enthalpies of the reference state at the current temperature.
* Real reason for its existence is that it also checks
* to see if a recalculation of the reference thermodynamics
* functions needs to be done.
*/
const vector_fp& IdealSolidSolnPhase::enthalpy_RT_ref() const
{
_updateThermo();
return m_h0_RT;
}
/*
* Returns a reference to the vector of nondimensional
* enthalpies of the reference state at the current temperature.
* Real reason for its existence is that it also checks
* to see if a recalculation of the reference thermodynamics
* functions needs to be done.
*/
const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const
{
_updateThermo();
@ -1051,41 +551,11 @@ const vector_fp& IdealSolidSolnPhase::entropy_R_ref() const
/*********************************************************************
* Utility Functions
*********************************************************************/
/*
* initThermo() function initializes the object for use.
*
* Before its invocation, the class isn't ready for calculation.
*/
void IdealSolidSolnPhase::initThermo()
{
}
/*
* @internal
* Import and initialize a ThermoPhase object
* using an XML tree.
* Here we read extra information about the XML description
* of a phase. Regular information about elements and species
* and their reference state thermodynamic information
* have already been read at this point.
* For example, we do not need to call this function for
* ideal gas equations of state.
* This function is called from importPhase()
* after the elements and the
* species are initialized with default ideal solution
* level data.
*
* @param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* @param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*/
void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
string subname = "IdealSolidSolnPhase::initThermoXML";
@ -1166,9 +636,6 @@ void IdealSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string&
ThermoPhase::initThermoXML(phaseNode, id);
}
/*
* This internal function adjusts the lengths of arrays
*/
void IdealSolidSolnPhase::
initLengths()
{
@ -1190,14 +657,6 @@ initLengths()
m_speciesMolarVolume.resize(m_kk);
}
/*
* Set mixture to an equilibrium state consistent with specified
* element potentials and temperature.
*
* @param lambda_RT vector of non-dimensional element potentials
* \f$ \lambda_m/RT \f$.
*
*/
void IdealSolidSolnPhase::
setToEquilState(const doublereal* lambda_RT)
{
@ -1217,45 +676,19 @@ setToEquilState(const doublereal* lambda_RT)
doublereal* dptr = DATA_PTR(m_pp);
setState_PX(pres, dptr);
}
//================================================================================================
/*
*
* speciesMolarVolume()
*
* Report the molar volume of species k
*
* units - \f$ m^3 kmol^-1 \f$
*/
double IdealSolidSolnPhase::
speciesMolarVolume(int k) const
{
return m_speciesMolarVolume[k];
}
/*
*
* getSpeciesMolarVolumes():
*
* Fill in a return vector containing the species molar volumes
* units - \f$ m^3 kmol^-1 \f$
*/
void IdealSolidSolnPhase::
getSpeciesMolarVolumes(doublereal* smv) const
{
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), smv);
}
//================================================================================================
/*
*
* _updateThermo()
*
* This function gets called for every call to functions in this
* class. It checks to see whether the temperature has changed and
* thus the reference thermodynamics functions for all of the species
* must be recalculated.
* If the temperature has changed, the species thermo manager is called
* to recalculate G, Cp, H, and S at the current temperature.
*/
void IdealSolidSolnPhase::
_updateThermo() const
{
@ -1277,6 +710,5 @@ _updateThermo() const
m_tlast = tnow;
}
}
//================================================================================================
} // end namespace Cantera
//==================================================================================================

148
src/thermo/LatticePhase.cpp
View file

@ -21,7 +21,6 @@
namespace Cantera
{
// Base Empty constructor
LatticePhase::LatticePhase() :
m_Pref(OneAtm),
m_Pcurrent(OneAtm),
@ -31,10 +30,6 @@ LatticePhase::LatticePhase() :
{
}
// Copy Constructor
/*
* @param right Object to be copied
*/
LatticePhase::LatticePhase(const LatticePhase& right) :
m_Pref(OneAtm),
m_Pcurrent(OneAtm),
@ -45,10 +40,6 @@ LatticePhase::LatticePhase(const LatticePhase& right) :
*this = operator=(right);
}
// Assignment operator
/*
* @param right Object to be copied
*/
LatticePhase& LatticePhase::operator=(const LatticePhase& right)
{
if (&right != this) {
@ -67,41 +58,20 @@ LatticePhase& LatticePhase::operator=(const LatticePhase& right)
return *this;
}
// Destructor
LatticePhase::~LatticePhase()
{
}
// Full constructor for a lattice phase
/*
* @param inputFile String name of the input file
* @param id string id of the phase name
*/
LatticePhase::LatticePhase(const std::string& inputFile, const std::string& id)
{
initThermoFile(inputFile, id);
}
// Full constructor for a water phase
/*
* @param phaseRef XML node referencing the lattice phase.
* @param id string id of the phase name
*/
LatticePhase::LatticePhase(XML_Node& phaseRef, const std::string& id)
{
importPhase(*findXMLPhase(&phaseRef, id), this);
}
// Duplication function
/*
* This virtual function is used to create a duplicate of the
* current phase. It's used to duplicate the phase when given
* a ThermoPhase pointer to the phase.
*
* @return It returns a ThermoPhase pointer.
*/
ThermoPhase* LatticePhase::duplMyselfAsThermoPhase() const
{
return new LatticePhase(*this);
@ -128,22 +98,22 @@ doublereal LatticePhase::entropy_mole() const
return GasConstant * (mean_X(&entropy_R_ref()[0]) -
sum_xlogx());
}
//====================================================================================================================
doublereal LatticePhase::gibbs_mole() const
{
return enthalpy_mole() - temperature() * entropy_mole();
}
//====================================================================================================================
doublereal LatticePhase::cp_mole() const
{
return GasConstant * mean_X(&cp_R_ref()[0]);
}
//====================================================================================================================
doublereal LatticePhase::cv_mole() const
{
return cp_mole();
}
//====================================================================================================================
doublereal LatticePhase::calcDensity()
{
setMolarDensity(m_site_density);
@ -162,65 +132,65 @@ doublereal LatticePhase::calcDensity()
// Phase::setDensity(dens);
return dens;
}
//====================================================================================================================
void LatticePhase::setPressure(doublereal p)
{
m_Pcurrent = p;
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMoleFractions(const doublereal* const x)
{
Phase::setMoleFractions(x);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMoleFractions_NoNorm(const doublereal* const x)
{
Phase::setMoleFractions(x);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMassFractions(const doublereal* const y)
{
Phase::setMassFractions(y);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setMassFractions_NoNorm(const doublereal* const y)
{
Phase::setMassFractions_NoNorm(y);
calcDensity();
}
//====================================================================================================================
void LatticePhase::setConcentrations(const doublereal* const c)
{
Phase::setConcentrations(c);
calcDensity();
}
//====================================================================================================================
void LatticePhase::getActivityConcentrations(doublereal* c) const
{
getMoleFractions(c);
}
//====================================================================================================================
void LatticePhase::getActivityCoefficients(doublereal* ac) const
{
for (size_t k = 0; k < m_kk; k++) {
ac[k] = 1.0;
}
}
//====================================================================================================================
doublereal LatticePhase::standardConcentration(size_t k) const
{
return 1.0;
}
//====================================================================================================================
doublereal LatticePhase::logStandardConc(size_t k) const
{
return 0.0;
}
//====================================================================================================================
void LatticePhase::getChemPotentials(doublereal* mu) const
{
doublereal delta_p = m_Pcurrent - m_Pref;
@ -234,14 +204,14 @@ void LatticePhase::getChemPotentials(doublereal* mu) const
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
const vector_fp& _h = enthalpy_RT_ref();
doublereal rt = GasConstant * temperature();
scale(_h.begin(), _h.end(), hbar, rt);
}
//====================================================================================================================
void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const
{
const vector_fp& _s = entropy_R_ref();
@ -252,7 +222,7 @@ void LatticePhase::getPartialMolarEntropies(doublereal* sbar) const
sbar[k] = r * (_s[k] - log(xx));
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarCp(doublereal* cpbar) const
{
getCp_R(cpbar);
@ -260,18 +230,18 @@ void LatticePhase::getPartialMolarCp(doublereal* cpbar) const
cpbar[k] *= GasConstant;
}
}
//====================================================================================================================
void LatticePhase::getPartialMolarVolumes(doublereal* vbar) const
{
getStandardVolumes(vbar);
}
//====================================================================================================================
void LatticePhase::getStandardChemPotentials(doublereal* mu0) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), mu0, _RT());
}
//====================================================================================================================
void LatticePhase::getPureGibbs(doublereal* gpure) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
@ -281,7 +251,7 @@ void LatticePhase::getPureGibbs(doublereal* gpure) const
gpure[k] = RT * gibbsrt[k] + delta_p * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getEnthalpy_RT(doublereal* hrt) const
{
const vector_fp& _h = enthalpy_RT_ref();
@ -290,13 +260,13 @@ void LatticePhase::getEnthalpy_RT(doublereal* hrt) const
hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getEntropy_R(doublereal* sr) const
{
const vector_fp& _s = entropy_R_ref();
std::copy(_s.begin(), _s.end(), sr);
}
//====================================================================================================================
void LatticePhase::getGibbs_RT(doublereal* grt) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
@ -306,7 +276,7 @@ void LatticePhase::getGibbs_RT(doublereal* grt) const
grt[k] = gibbsrt[k] + delta_prt * m_speciesMolarVolume[k];
}
}
//====================================================================================================================
void LatticePhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
@ -314,42 +284,30 @@ void LatticePhase::getGibbs_ref(doublereal* g) const
g[k] *= GasConstant * temperature();
}
}
//===================================================================================================================
void LatticePhase::getCp_R(doublereal* cpr) const
{
const vector_fp& _cpr = cp_R_ref();
std::copy(_cpr.begin(), _cpr.end(), cpr);
}
//===================================================================================================================
void LatticePhase::getStandardVolumes(doublereal* vbar) const
{
copy(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), vbar);
}
//=======================================================================================================
// Returns the vector of nondimensional Enthalpies of the reference state at the current temperature
// of the solution and the reference pressure for the phase.
/*
* @return Output vector of nondimensional reference state Enthalpies of the species.
* Length: m_kk
*/
const vector_fp& LatticePhase::enthalpy_RT_ref() const
{
_updateThermo();
return m_h0_RT;
}
//=======================================================================================================
// Returns a reference to the dimensionless reference state Gibbs free energy vector.
/*
* This function is part of the layer that checks/recalculates the reference
* state thermo functions.
*/
const vector_fp& LatticePhase::gibbs_RT_ref() const
{
_updateThermo();
return m_g0_RT;
}
//====================================================================================================================
void LatticePhase::getGibbs_RT_ref(doublereal* grt) const
{
_updateThermo();
@ -357,43 +315,19 @@ void LatticePhase::getGibbs_RT_ref(doublereal* grt) const
grt[k] = m_g0_RT[k];
}
}
//=======================================================================================================
// Returns a reference to the dimensionless reference state Entropy vector.
/*
* This function is part of the layer that checks/recalculates the reference
* state thermo functions.
*/
const vector_fp& LatticePhase::entropy_R_ref() const
{
_updateThermo();
return m_s0_R;
}
//=======================================================================================================
// Returns a reference to the dimensionless reference state Heat Capacity vector.
/*
* This function is part of the layer that checks/recalculates the reference
* state thermo functions.
*/
const vector_fp& LatticePhase::cp_R_ref() const
{
_updateThermo();
return m_cp0_R;
}
//====================================================================================================================
// Initialize the ThermoPhase object after all species have been set up
/*
* @internal Initialize.
*
* This method performs any initialization required after all
* species have been added. For example, it is used to
* resize internal work arrays that must have an entry for
* each species.
* This method is called from ThermoPhase::initThermoXML(),
* which is called from importPhase(),
* just prior to returning from the function, importPhase().
*
* @see importCTML.cpp
*/
void LatticePhase::initThermo()
{
m_Pref = refPressure();
@ -406,7 +340,7 @@ void LatticePhase::initThermo()
ThermoPhase::initThermo();
}
//====================================================================================================================
void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
std::string idattrib = phaseNode.id();
@ -459,12 +393,7 @@ void LatticePhase::initThermoXML(XML_Node& phaseNode, const std::string& id)
*/
ThermoPhase::initThermoXML(phaseNode, id);
}
//=====================================================================================================
// Update the species reference state thermodynamic functions
/*
* The polynomials for the standard state functions are only
* reevaluated if the temperature has changed.
*/
void LatticePhase::_updateThermo() const
{
doublereal tnow = temperature();
@ -477,26 +406,25 @@ void LatticePhase::_updateThermo() const
m_tlast = tnow;
}
}
//=====================================================================================================
void LatticePhase::setParameters(int n, doublereal* const c)
{
m_site_density = c[0];
setMolarDensity(m_site_density);
}
//=====================================================================================================
void LatticePhase::getParameters(int& n, doublereal* const c) const
{
double d = molarDensity();
c[0] = d;
n = 1;
}
//=====================================================================================================
void LatticePhase::setParametersFromXML(const XML_Node& eosdata)
{
eosdata._require("model", "Lattice");
m_site_density = ctml::getFloat(eosdata, "site_density", "toSI");
m_vacancy = ctml::getChildValue(eosdata, "vacancy_species");
}
//=====================================================================================================
}
//=======================================================================================================

222
src/thermo/LatticeSolidPhase.cpp
View file

@ -18,12 +18,9 @@
#include <string>
using namespace std;
//======================================================================================================================
namespace Cantera
{
//====================================================================================================================
// Base empty constructor
LatticeSolidPhase::LatticeSolidPhase() :
m_tlast(0.0),
m_press(-1.0),
@ -35,11 +32,7 @@ LatticeSolidPhase::LatticeSolidPhase() :
tmpV_(0)
{
}
//====================================================================================================================
// Copy Constructor
/*
* @param right Object to be copied
*/
LatticeSolidPhase::LatticeSolidPhase(const LatticeSolidPhase& right) :
m_tlast(0.0),
m_press(-1.0),
@ -52,11 +45,7 @@ LatticeSolidPhase::LatticeSolidPhase(const LatticeSolidPhase& right) :
{
*this = operator=(right);
}
//====================================================================================================================
// Assignment operator
/*
* @param right Object to be copied
*/
LatticeSolidPhase&
LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
{
@ -73,8 +62,7 @@ LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
}
return *this;
}
//====================================================================================================================
// Destructor
LatticeSolidPhase::~LatticeSolidPhase()
{
// We own the sublattices. So we have to delete the sublattices
@ -83,33 +71,12 @@ LatticeSolidPhase::~LatticeSolidPhase()
m_lattice[n] = 0;
}
}
//====================================================================================================================
// Duplication function
/*
* This virtual function is used to create a duplicate of the
* current phase. It's used to duplicate the phase when given
* a ThermoPhase pointer to the phase.
*
* @return It returns a %ThermoPhase pointer.
*/
ThermoPhase* LatticeSolidPhase::duplMyselfAsThermoPhase() const
{
return new LatticeSolidPhase(*this);
}
//====================================================================================================================
// Minimum temperature for which the thermodynamic data for the species
// or phase are valid.
/*
* If no argument is supplied, the
* value returned will be the lowest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo minTemp function.
*
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
doublereal LatticeSolidPhase::minTemp(size_t k) const
{
if (k != npos) {
@ -127,19 +94,7 @@ doublereal LatticeSolidPhase::minTemp(size_t k) const
}
return mm;
}
//====================================================================================================================
// Maximum temperature for which the thermodynamic data for the species
// or phase are valid.
/*
* If no argument is supplied, the
* value returned will be the lowest temperature at which the
* data for \e all species are valid. Otherwise, the value
* will be only for species \a k. This function is a wrapper
* that calls the species thermo minTemp function.
*
* @param k index of the species. Default is -1, which will return the max of the min value
* over all species.
*/
doublereal LatticeSolidPhase::maxTemp(size_t k) const
{
if (k != npos) {
@ -157,18 +112,13 @@ doublereal LatticeSolidPhase::maxTemp(size_t k) const
}
return mm;
}
//====================================================================================================================
/*
* Returns the reference pressure in Pa. This function is a wrapper
* that calls the species thermo refPressure function.
*/
doublereal LatticeSolidPhase::refPressure() const
{
return m_lattice[0]->refPressure();
}
//====================================================================================================================
doublereal LatticeSolidPhase::
enthalpy_mole() const
doublereal LatticeSolidPhase::enthalpy_mole() const
{
_updateThermo();
doublereal sum = 0.0;
@ -177,7 +127,7 @@ enthalpy_mole() const
}
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::intEnergy_mole() const
{
_updateThermo();
@ -187,7 +137,7 @@ doublereal LatticeSolidPhase::intEnergy_mole() const
}
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::entropy_mole() const
{
_updateThermo();
@ -197,7 +147,7 @@ doublereal LatticeSolidPhase::entropy_mole() const
}
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::gibbs_mole() const
{
_updateThermo();
@ -207,7 +157,7 @@ doublereal LatticeSolidPhase::gibbs_mole() const
}
return sum;
}
//====================================================================================================================
doublereal LatticeSolidPhase::cp_mole() const
{
_updateThermo();
@ -217,7 +167,7 @@ doublereal LatticeSolidPhase::cp_mole() const
}
return sum;
}
//====================================================================================================================
void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const
{
_updateThermo();
@ -227,30 +177,24 @@ void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const
strt += m_lattice[n]->nSpecies();
}
}
//====================================================================================================================
void LatticeSolidPhase::getActivityCoefficients(doublereal* ac) const
{
for (size_t k = 0; k < m_kk; k++) {
ac[k] = 1.0;
}
}
//====================================================================================================================
doublereal LatticeSolidPhase::standardConcentration(size_t k) const
{
return 1.0;
}
//====================================================================================================================
doublereal LatticeSolidPhase::logStandardConc(size_t k) const
{
return 0.0;
}
//====================================================================================================================
// Set the pressure at constant temperature. Units: Pa.
/*
*
* @param p Pressure (units - Pa)
*/
void LatticeSolidPhase::setPressure(doublereal p)
{
m_press = p;
@ -259,19 +203,7 @@ void LatticeSolidPhase::setPressure(doublereal p)
}
calcDensity();
}
//====================================================================================================================
// Calculate the density of the solid mixture
/*
* The formula for this is
*
* \f[
* \rho = \sum_n{ \rho_n \theta_n }
* \f]
*
* where \f$ \rho_n \f$ is the density of the nth sublattice
*
* Note this is a nonvirtual function.
*/
doublereal LatticeSolidPhase::calcDensity()
{
double sum = 0.0;
@ -281,20 +213,7 @@ doublereal LatticeSolidPhase::calcDensity()
Phase::setDensity(sum);
return sum;
}
//====================================================================================================================
// Set the mole fractions to the specified values, and then
// normalize them so that they sum to 1.0 for each of the subphases
/*
* On input, the mole fraction vector is assumed to sum to one for each of the sublattices. The sublattices
* are updated with this mole fraction vector. The mole fractions are also stored within this object, after
* they are normalized to one by dividing by the number of sublattices.
*
* @param x Input vector of mole fractions. There is no restriction
* on the sum of the mole fraction vector. Internally,
* this object will pass portions of this vector to the sublattices which assume that the portions
* individually sum to one.
* Length is m_kk.
*/
void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
{
size_t nsp, strt = 0;
@ -309,14 +228,7 @@ void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
Phase::setMoleFractions(DATA_PTR(m_x));
calcDensity();
}
//====================================================================================================================
// Get the species mole fraction vector.
/*
* On output the mole fraction vector will sum to one for each of the subphases which make up this phase.
*
* @param x On return, x contains the mole fractions. Must have a
* length greater than or equal to the number of species.
*/
void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
{
size_t nsp, strt = 0;
@ -348,18 +260,7 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
strt += nsp;
}
}
//====================================================================================================================
// Get the species chemical potentials. Units: J/kmol.
/*
* This function returns a vector of chemical potentials of the
* species in solution at the current temperature, pressure
* and mole fraction of the solution.
*
* This returns the underlying lattice chemical potentials
*
* @param mu Output vector of species chemical
* potentials. Length: m_kk. Units: J/kmol
*/
void LatticeSolidPhase::getChemPotentials(doublereal* mu) const
{
_updateThermo();
@ -370,7 +271,7 @@ void LatticeSolidPhase::getChemPotentials(doublereal* mu) const
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
_updateThermo();
@ -381,7 +282,7 @@ void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const
{
_updateThermo();
@ -392,7 +293,7 @@ void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const
{
_updateThermo();
@ -403,7 +304,7 @@ void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const
strt += nlsp;
}
}
//====================================================================================================================
void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const
{
_updateThermo();
@ -414,20 +315,7 @@ void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const
strt += nlsp;
}
}
//====================================================================================================================
// Get the array of standard state chemical potentials at unit activity for the species
// at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current
* temperature and pressure of the solution.
*
* This returns the underlying lattice standard chemical potentials, as the units are kmol-1 of
* the sublattice species.
*
* @param mu0 Output vector of chemical potentials.
* Length: m_kk. Units: J/kmol
*/
void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
{
_updateThermo();
@ -437,7 +325,7 @@ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
strt += m_lattice[n]->nSpecies();
}
}
//====================================================================================================================
void LatticeSolidPhase::getGibbs_RT_ref(doublereal* grt) const
{
_updateThermo();
@ -445,7 +333,7 @@ void LatticeSolidPhase::getGibbs_RT_ref(doublereal* grt) const
m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
}
}
//====================================================================================================================
void LatticeSolidPhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
@ -453,13 +341,7 @@ void LatticeSolidPhase::getGibbs_ref(doublereal* g) const
g[k] *= GasConstant * temperature();
}
}
//====================================================================================================================
// Add in species from Slave phases
/*
* This hook is used for cSS_CONVENTION_SLAVE phases
*
* @param phaseNode XML_Node for the current phase
*/
void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
{
size_t kk = 0;
@ -537,22 +419,6 @@ void LatticeSolidPhase::installSlavePhases(Cantera::XML_Node* phaseNode)
}
}
//====================================================================================================================
// Initialize the ThermoPhase object after all species have been set up
/*
* @internal Initialize.
*
* This method is provided to allow subclasses to perform any initialization required after all
* species have been added. For example, it might be used to
* resize internal work arrays that must have an entry for
* each species. The base class implementation does nothing,
* and subclasses that do not require initialization do not
* need to overload this method. When importing a CTML phase
* description, this method is called from ThermoPhase::initThermoXML(),
* which is called from importPhase(), just prior to returning from function importPhase().
*
* @see importCTML.cpp
*/
void LatticeSolidPhase::initThermo()
{
initLengths();
@ -569,11 +435,7 @@ void LatticeSolidPhase::initThermo()
setMoleFractions(DATA_PTR(m_x));
ThermoPhase::initThermo();
}
//====================================================================================================================
// Initialize vectors that depend on the number of species and sublattices
/*
*
*/
void LatticeSolidPhase::initLengths()
{
theta_.resize(m_nlattice,0);
@ -581,7 +443,7 @@ void LatticeSolidPhase::initLengths()
m_x.resize(m_kk, 0.0);
tmpV_.resize(m_kk, 0.0);
}
//====================================================================================================================
void LatticeSolidPhase::_updateThermo() const
{
doublereal tnow = temperature();
@ -601,7 +463,7 @@ void LatticeSolidPhase::_updateThermo() const
m_tlast = tnow;
}
}
//====================================================================================================================
void LatticeSolidPhase::setLatticeMoleFractionsByName(int nn, const std::string& x)
{
m_lattice[nn]->setMoleFractionsByName(x);
@ -617,15 +479,7 @@ void LatticeSolidPhase::setLatticeMoleFractionsByName(int nn, const std::string&
}
setMoleFractions(DATA_PTR(m_x));
}
//====================================================================================================================
//====================================================================================================================
// Set the parameters from the XML file
/*!
* Currently, this is the spot that we read in all of the sublattice phases.
* The SetParametersFromXML() call is carried out at
*/
void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
{
eosdata._require("model","LatticeSolid");
@ -661,18 +515,8 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
}
}
//====================================================================================================================
#ifdef H298MODIFY_CAPABILITY
//! Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
/*!
* The 298K heat of formation is defined as the enthalpy change to create the standard state
* of the species from its constituent elements in their standard states at 298 K and 1 bar.
*
* @param k Species k
* @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar
*/
void LatticeSolidPhase::modifyOneHf298SS(const size_t& k, const doublereal Hf298New)
{
for (size_t n = 0; n < m_nlattice; n++) {
@ -686,7 +530,6 @@ void LatticeSolidPhase::modifyOneHf298SS(const size_t& k, const doublereal Hf298
_updateThermo();
}
#endif
//====================================================================================================================
doublereal LatticeSolidPhase::err(const std::string& msg) const
{
@ -695,4 +538,3 @@ doublereal LatticeSolidPhase::err(const std::string& msg) const
}
} // End namespace Cantera
//======================================================================================================================

441
src/thermo/MixtureFugacityTP.cpp
View file

@ -3,7 +3,6 @@
* Methods file for a derived class of ThermoPhase that handles
* non-ideal mixtures based on the fugacity models (see \ref thermoprops and
* class \link Cantera::MixtureFugacityTP MixtureFugacityTP\endlink).
*
*/
/*
* Copyright (2005) Sandia Corporation. Under the terms of
@ -20,10 +19,7 @@ using namespace std;
namespace Cantera
{
//====================================================================================================================
/*
* Default constructor
*/
MixtureFugacityTP::MixtureFugacityTP() :
ThermoPhase(),
m_Pcurrent(-1.0),
@ -38,16 +34,7 @@ MixtureFugacityTP::MixtureFugacityTP() :
m_s0_R(0)
{
}
//====================================================================================================================
/*
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor.
*
* The copy constructor just calls the assignment operator
* to do the heavy lifting.
*/
MixtureFugacityTP::MixtureFugacityTP(const MixtureFugacityTP& b) :
ThermoPhase(),
m_Pcurrent(-1.0),
@ -63,13 +50,7 @@ MixtureFugacityTP::MixtureFugacityTP(const MixtureFugacityTP& b) :
{
MixtureFugacityTP::operator=(b);
}
//====================================================================================================================
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
MixtureFugacityTP&
MixtureFugacityTP::operator=(const MixtureFugacityTP& b)
{
@ -92,102 +73,43 @@ MixtureFugacityTP::operator=(const MixtureFugacityTP& b)
m_cp0_R = b.m_cp0_R;
m_g0_RT = b.m_g0_RT;
m_s0_R = b.m_s0_R;
/*
* The VPSSMgr object contains shallow pointers. Whenever you have shallow
* pointers, they have to be fixed up to point to the correct objects referring
* back to this ThermoPhase's properties.
*/
//m_VPSS_ptr->initAllPtrs(this, m_spthermo);
/*
* The PDSS objects contains shallow pointers. Whenever you have shallow
* pointers, they have to be fixed up to point to the correct objects referring
* back to this ThermoPhase's properties. This function also sets m_VPSS_ptr
* so it occurs after m_VPSS_ptr is set.
*/
/*
* Ok, the VPSSMgr object is ready for business.
* We need to resync the temperature and the pressure of the new standard states
* with what is stored in this object.
*/
// m_VPSS_ptr->setState_TP(m_Tlast_ss, m_Plast_ss);
}
return *this;
}
//====================================================================================================================
/*
* ~MixtureFugacityTP(): (virtual)
*
*/
MixtureFugacityTP::~MixtureFugacityTP()
{
}
/*
* Duplication function.
* This calls the copy constructor for this object.
*/
ThermoPhase* MixtureFugacityTP::duplMyselfAsThermoPhase() const
{
return new MixtureFugacityTP(*this);
}
//====================================================================================================================
// This method returns the convention used in specification
// of the standard state, of which there are currently two,
// temperature based, and variable pressure based.
/*
* Currently, there are two standard state conventions:
* - Temperature-based activities
* cSS_CONVENTION_TEMPERATURE 0
* - default
*
* - Variable Pressure and Temperature -based activities
* cSS_CONVENTION_VPSS 1
*/
int MixtureFugacityTP::standardStateConvention() const
{
return cSS_CONVENTION_TEMPERATURE;
}
//====================================================================================================================
// Set the solution branch to force the ThermoPhase to exist on one branch or another
/*
* @param solnBranch Branch that the solution is restricted to.
* the value -1 means gas. The value -2 means unrestricted.
* Values of zero or greater refer to species dominated condensed phases.
*/
void MixtureFugacityTP::setForcedSolutionBranch(int solnBranch)
{
forcedState_ = solnBranch;
}
//====================================================================================================================
// Report the solution branch which the solution is restricted to
/*
* @return Branch that the solution is restricted to.
* the value -1 means gas. The value -2 means unrestricted.
* Values of zero or greater refer to species dominated condensed phases.
*/
int MixtureFugacityTP::forcedSolutionBranch() const
{
return forcedState_;
}
//====================================================================================================================
// Report the solution branch which the solution is actually on
/*
* @return Branch that the solution is restricted to.
* the value -1 means gas. The value -2 means superfluid..
* Values of zero or greater refer to species dominated condensed phases.
*/
int MixtureFugacityTP::reportSolnBranchActual() const
{
return iState_;
}
//====================================================================================================================
/*
* ------------Molar Thermodynamic Properties -------------------------
*/
//====================================================================================================================
doublereal MixtureFugacityTP::err(const std::string& msg) const
{
@ -195,20 +117,11 @@ doublereal MixtureFugacityTP::err(const std::string& msg) const
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
//====================================================================================================================
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
//====================================================================================================================
/*
* Get the array of non-dimensional species chemical potentials
* These are partial molar Gibbs free energies.
* \f$ \mu_k / \hat R T \f$.
* Units: unitless
*
* We close the loop on this function, here, calling
* getChemPotentials() and then dividing by RT.
*/
void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
{
getChemPotentials(muRT);
@ -217,10 +130,11 @@ void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
muRT[k] *= invRT;
}
}
//====================================================================================================================
/*
* ----- Thermodynamic Values for the Species Standard States States ----
*/
void MixtureFugacityTP::getStandardChemPotentials(doublereal* g) const
{
_updateReferenceStateThermo();
@ -231,33 +145,20 @@ void MixtureFugacityTP::getStandardChemPotentials(doublereal* g) const
g[k] = RT * (g[k] + tmp);
}
}
//====================================================================================================================
void MixtureFugacityTP::getEnthalpy_RT(doublereal* hrt) const
{
getEnthalpy_RT_ref(hrt);
}
//================================================================================================
#ifdef H298MODIFY_CAPABILITY
// Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1)
/*
* The 298K heat of formation is defined as the enthalpy change to create the standard state
* of the species from its constituent elements in their standard states at 298 K and 1 bar.
*
* @param k Species k
* @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar
*/
void MixtureFugacityTP::modifyOneHf298SS(const int k, const doublereal Hf298New)
{
m_spthermo->modifyOneHf298(k, Hf298New);
m_Tlast_ref += 0.0001234;
}
#endif
//====================================================================================================================
/*
* Get the array of nondimensional entropy functions for the
* standard state species
* at the current <I>T</I> and <I>P</I> of the solution.
*/
void MixtureFugacityTP::getEntropy_R(doublereal* sr) const
{
_updateReferenceStateThermo();
@ -267,11 +168,7 @@ void MixtureFugacityTP::getEntropy_R(doublereal* sr) const
sr[k] -= tmp;
}
}
//====================================================================================================================
/*
* Get the nondimensional gibbs function for the species
* standard states at the current T and P of the solution.
*/
void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const
{
_updateReferenceStateThermo();
@ -281,12 +178,7 @@ void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const
grt[k] += tmp;
}
}
//====================================================================================================================
/*
* get the pure Gibbs free energies of each species assuming
* it is in its standard state. This is the same as
* getStandardChemPotentials().
*/
void MixtureFugacityTP::getPureGibbs(doublereal* g) const
{
_updateReferenceStateThermo();
@ -297,12 +189,7 @@ void MixtureFugacityTP::getPureGibbs(doublereal* g) const
g[k] += tmp;
}
}
//====================================================================================================================
/*
* Returns the vector of nondimensional
* internal Energies of the standard state at the current temperature
* and pressure of the solution for each species.
*/
void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const
{
_updateReferenceStateThermo();
@ -314,26 +201,13 @@ void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const
urt[i] -= tmp * v0;
}
}
//====================================================================================================================
/*
* Get the nondimensional heat capacity at constant pressure
* function for the species
* standard states at the current T and P of the solution.
*/
void MixtureFugacityTP::getCp_R(doublereal* cpr) const
{
_updateReferenceStateThermo();
copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
}
//====================================================================================================================
/*
* Get the molar volumes of the species standard states at the current
* <I>T</I> and <I>P</I> of the solution.
* units = m^3 / kmol
*
* @param vol Output vector containing the standard state volumes.
* Length: m_kk.
*/
void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const
{
_updateReferenceStateThermo();
@ -342,84 +216,49 @@ void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const
vol[i]= v0;
}
}
//====================================================================================================================
/*
* ----- Thermodynamic Values for the Species Reference States ----
*/
/*
* Returns the vector of nondimensional enthalpies of the
* reference state at the current temperature of the solution and
* the reference pressure for the species.
*/
void MixtureFugacityTP::getEnthalpy_RT_ref(doublereal* hrt) const
{
_updateReferenceStateThermo();
copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
}
//====================================================================================================================
/*
* Returns the vector of nondimensional
* enthalpies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
void MixtureFugacityTP::getGibbs_RT_ref(doublereal* grt) const
{
_updateReferenceStateThermo();
copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
}
//====================================================================================================================
/*
* Returns the vector of the
* gibbs function of the reference state at the current temperature
* of the solution and the reference pressure for the species.
* units = J/kmol
*
* This is filled in here so that derived classes don't have to
* take care of it.
*/
void MixtureFugacityTP::getGibbs_ref(doublereal* g) const
{
const vector_fp& gibbsrt = gibbs_RT_ref();
scale(gibbsrt.begin(), gibbsrt.end(), g, _RT());
}
//====================================================================================================================
const vector_fp& MixtureFugacityTP::gibbs_RT_ref() const
{
_updateReferenceStateThermo();
return m_g0_RT;
}
//====================================================================================================================
/*
* Returns the vector of nondimensional
* entropies of the reference state at the current temperature
* of the solution and the reference pressure for the species.
*/
void MixtureFugacityTP::getEntropy_R_ref(doublereal* er) const
{
_updateReferenceStateThermo();
copy(m_s0_R.begin(), m_s0_R.end(), er);
return;
}
//====================================================================================================================
/*
* Returns the vector of nondimensional
* constant pressure heat capacities of the reference state
* at the current temperature of the solution
* and reference pressure for the species.
*/
void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const
{
_updateReferenceStateThermo();
copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
}
//====================================================================================================================
/*
* Get the molar volumes of the species reference states at the current
* <I>T</I> and reference pressure of the solution.
*
* units = m^3 / kmol
*/
void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const
{
_updateReferenceStateThermo();
@ -429,17 +268,7 @@ void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const
vol[i]= v0;
}
}
//====================================================================================================================
// Set the initial state of the phase to the conditions specified in the state XML element.
/*
*
* This method sets the temperature, pressure, and mole fraction vector to a set default value.
* We modify the default behavior here so that TP is evaluated at the same time.
*
* @param state AN XML_Node object corresponding to
* the "state" entry for this phase in the
* input file.
*/
void MixtureFugacityTP::setStateFromXML(const XML_Node& state)
{
int doTP = 0;
@ -472,21 +301,13 @@ void MixtureFugacityTP::setStateFromXML(const XML_Node& state)
setState_TR(t, rho);
}
}
//====================================================================================================================
/*
* Perform initializations after all species have been
* added.
*/
void MixtureFugacityTP::initThermo()
{
initLengths();
ThermoPhase::initThermo();
}
//====================================================================================================================
/*
* Initialize the internal lengths.
* (this is not a virtual function)
*/
void MixtureFugacityTP::initLengths()
{
m_kk = nSpecies();
@ -497,13 +318,13 @@ void MixtureFugacityTP::initLengths()
m_g0_RT.resize(m_kk, 0.0);
m_s0_R.resize(m_kk, 0.0);
}
//====================================================================================================================
void MixtureFugacityTP::setTemperature(const doublereal temp)
{
_updateReferenceStateThermo();
setState_TR(temperature(), density());
}
//====================================================================================================================
void MixtureFugacityTP::setPressure(doublereal p)
{
setState_TP(temperature(), p);
@ -514,49 +335,48 @@ void MixtureFugacityTP::setPressure(doublereal p)
// printf(" MixFug:setPres: mu(%d = %g) = %18.8g\n", i, mf[i], chemPot[i]);
// }
}
//====================================================================================================================
void MixtureFugacityTP::setMassFractions(const doublereal* const y)
{
Phase::setMassFractions(y);
getMoleFractions(DATA_PTR(moleFractions_));
}
//====================================================================================================================
void MixtureFugacityTP::setMassFractions_NoNorm(const doublereal* const y)
{
Phase::setMassFractions_NoNorm(y);
getMoleFractions(DATA_PTR(moleFractions_));
}
//====================================================================================================================
void MixtureFugacityTP::setMoleFractions(const doublereal* const x)
{
Phase::setMoleFractions(x);
getMoleFractions(DATA_PTR(moleFractions_));
}
//====================================================================================================================
void MixtureFugacityTP::setMoleFractions_NoNorm(const doublereal* const x)
{
Phase::setMoleFractions_NoNorm(x);
getMoleFractions(DATA_PTR(moleFractions_));
}
//====================================================================================================================
void MixtureFugacityTP::setConcentrations(const doublereal* const c)
{
Phase::setConcentrations(c);
getMoleFractions(DATA_PTR(moleFractions_));
}
//====================================================================================================================
void MixtureFugacityTP::setMoleFractions_NoState(const doublereal* const x)
{
Phase::setMoleFractions(x);
getMoleFractions(DATA_PTR(moleFractions_));
updateMixingExpressions();
}
//====================================================================================================================
void MixtureFugacityTP::calcDensity()
{
err("MixtureFugacityTP::calcDensity() called, but EOS for phase is not known");
}
//====================================================================================================================
void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres)
{
@ -647,15 +467,7 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres)
//setPressure(pres);
//calcDensity();
}
//====================================================================================================================
// Set the internally stored temperature (K) and density (kg/m^3)
/*
* This overrides the default behavior. In addition to just storing the state in the object, we need to do
* an equation of state calculation and figure out what phase state we are in.
*
* @param t Temperature in kelvin
* @param rho Density (kg/m^3)
*/
void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho)
{
getMoleFractions(DATA_PTR(moleFractions_));
@ -672,40 +484,12 @@ void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho)
// printf("setState_TR: state at T = %g, rho = %g, mv = %g, P = %20.13g, iState = %d\n", T, rho, mv, m_Pcurrent, iState_);
}
//====================================================================================================================
// Set the temperature (K), pressure (Pa), and mole fractions.
/*
* Note, the mole fractions are set first before the pressure is set.
* Setting the pressure may involve the solution of a nonlinear equation.
*
* @param t Temperature (K)
* @param p Pressure (Pa)
* @param x Vector of mole fractions.
* Length is equal to m_kk.
*/
void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublereal* x)
{
setMoleFractions_NoState(x);
setState_TP(t,p);
}
//====================================================================================================================
/*
* Import and initialize a ThermoPhase object
*
* param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*
* This routine initializes the lengths in the current object and
* then calls the parent routine.
*/
void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
MixtureFugacityTP::initLengths();
@ -715,7 +499,7 @@ void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id
// m_VPSS_ptr->initThermoXML(phaseNode, id);
ThermoPhase::initThermoXML(phaseNode, id);
}
//====================================================================================================================
doublereal MixtureFugacityTP::z() const
{
doublereal p = pressure();
@ -726,19 +510,19 @@ doublereal MixtureFugacityTP::z() const
doublereal zz = p * molarV / rt;
return zz;
}
//====================================================================================================================
doublereal MixtureFugacityTP::sresid() const
{
throw CanteraError("MixtureFugacityTP::sresid()", "Base Class: not implemented");
return 0.0;
}
//====================================================================================================================
doublereal MixtureFugacityTP::hresid() const
{
throw CanteraError("MixtureFugacityTP::hresid()", "Base Class: not implemented");
return 0.0;
}
//====================================================================================================================
doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
{
doublereal tcrit = critTemperature();
@ -750,32 +534,13 @@ doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
doublereal lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
return pcrit*exp(lpr);
}
//====================================================================================================================
doublereal MixtureFugacityTP::liquidVolEst(doublereal TKelvin, doublereal& pres) const
{
throw CanteraError("MixtureFugacityTP::liquidVolEst()", "unimplemented");
return 0.0;
}
//====================================================================================================================
/*
* Calculates the density given the temperature and the pressure,
* and a guess at the density. Note, below T_c, this is a
* multivalued function. This function assumes that the phase is on one side of the vapor dome
* or the other. It does not allow for crosses of the vapor dome.
*
* parameters:
* temperature: Kelvin
* pressure : Pressure in Pascals (Newton/m**2)
* phase : guessed phase of water
* : -1: no guessed phase
* rhoguess : guessed density of the water
*
* -1.0 no guessed density
*
* If a problem is encountered, a negative 1 is returned.
*
* @TODO make this a const function
*/
doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
int phase, doublereal rhoguess)
{
@ -956,18 +721,17 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
}
return densBase;
}
//====================================================================================================================
void MixtureFugacityTP::updateMixingExpressions()
{
}
//====================================================================================================================
MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
ResidEval(),
m_tp(tp)
{
}
//====================================================================================================================
int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal* const y,
doublereal* const r)
{
@ -979,7 +743,6 @@ int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal
r[0] = val;
return status;
}
//====================================================================================================================
int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
@ -1017,17 +780,7 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de
// delGRT = gibbsLiqRT - gibbsGasRT;
return retn;
}
//====================================================================================================================
// Returns the Phase State flag for the current state of the object
/*
* @param checkState If true, this function does a complete check to see where
* in parameter space we are
*
* 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_SUPERCRIT above the critical temperature
*/
int MixtureFugacityTP::phaseState(bool checkState) const
{
int state = iState_;
@ -1070,62 +823,40 @@ int MixtureFugacityTP::phaseState(bool checkState) const
}
return state;
}
//====================================================================================================================
// Return the value of the density at the liquid spinodal point (on the liquid side)
// for the current temperature.
/*
* @return returns the density with units of kg m-3
*/
doublereal MixtureFugacityTP::densSpinodalLiquid() const
{
throw CanteraError("", "unimplmented");
return 0.0;
}
//====================================================================================================================
// Return the value of the density at the gas spinodal point (on the gas side)
// for the current temperature.
/*
* @return returns the density with units of kg m-3
*/
doublereal MixtureFugacityTP::densSpinodalGas() const
{
throw CanteraError("", "unimplmented");
return 0.0;
}
//====================================================================================================================
// Calculate the saturation pressure at the current mixture content for the given temperature
/*
* This is a non-const routine that is public.
*
* The algorithm for this routine has undergone quite a bit of work. It probably needs more work.
* However, it seems now to be fairly robust.
* The key requirement is to find an initial pressure where both the liquid and the gas exist. This
* is not as easy as it sounds, and it gets exceedingly hard as the critical temperature is approached
* from below.
* Once we have this initial state, then we seek to equilibrate the gibbs free energies of the
* gas and liquid and use the formula
*
* dp = VdG
*
* to create an update condition for deltaP using
*
* - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
*
*
*
* @param TKelvin (input) Temperature (Kelvin)
* @param molarVolGas (return) Molar volume of the gas
* @param molarVolLiquid (return) Molar volume of the liquid
*
* @return Returns the saturation pressure at the given temperature
*
* @TODO Suggestions for the future would be to switch it to an algorithm that uses the gas molar volume
* and the liquid molar volumes as the fundamental unknowns.
*
*/
doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& molarVolGas,
doublereal& molarVolLiquid)
{
/*
* The algorithm for this routine has undergone quite a bit of work. It probably needs more work.
* However, it seems now to be fairly robust.
* The key requirement is to find an initial pressure where both the liquid and the gas exist. This
* is not as easy as it sounds, and it gets exceedingly hard as the critical temperature is approached
* from below.
* Once we have this initial state, then we seek to equilibrate the gibbs free energies of the
* gas and liquid and use the formula
*
* dp = VdG
*
* to create an update condition for deltaP using
*
* - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
*
* @TODO Suggestions for the future would be to switch it to an algorithm that uses the gas molar volume
* and the liquid molar volumes as the fundamental unknowns.
*/
// we need this because this is a non-const routine that is public
setTemperature(TKelvin);
double tcrit = critTemperature();
@ -1316,36 +1047,18 @@ doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& mola
return pres;
}
//====================================================================================================================
// Calculate the pressure given the temperature and the molar volume
doublereal MixtureFugacityTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
{
throw CanteraError("MixtureFugacityTP::pressureCalc", "unimplemented");
return 0.0;
}
//====================================================================================================================
// Calculate the pressure given the temperature and the molar volume
doublereal MixtureFugacityTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
{
throw CanteraError("MixtureFugacityTP::dpdVCalc", "unimplemented");
return 0.0;
}
//====================================================================================================================
/*
* void _updateStandardStateThermo() (protected, virtual, const)
*
* If m_useTmpStandardStateStorage is true,
* This function must be called for every call to functions in this
* class that need standard state properties.
* Child classes may require that it be called even if m_useTmpStandardStateStorage
* is not true.
* It checks to see whether the temperature has changed and
* thus the ss thermodynamics functions for all of the species
* must be recalculated.
*
* This
*/
void MixtureFugacityTP::_updateReferenceStateThermo() const
{
double Tnow = temperature();
@ -1367,9 +1080,5 @@ void MixtureFugacityTP::_updateReferenceStateThermo() const
m_logc0 = log(pref/(GasConstant * Tnow));
}
}
//====================================================================================================================
}

185
src/thermo/PureFluidPhase.cpp
View file

@ -22,7 +22,6 @@ using std::setw;
namespace Cantera
{
// Base Constructor
PureFluidPhase::PureFluidPhase() :
ThermoPhase(),
m_sub(0),
@ -32,7 +31,6 @@ PureFluidPhase::PureFluidPhase() :
{
}
// CopyConstructor
PureFluidPhase::PureFluidPhase(const PureFluidPhase& right) :
ThermoPhase(),
m_sub(0),
@ -43,10 +41,6 @@ PureFluidPhase::PureFluidPhase(const PureFluidPhase& right) :
*this = right;
}
//! Assignment operator
/*!
* @param right Object to be copied
*/
PureFluidPhase& PureFluidPhase::operator=(const PureFluidPhase& right)
{
if (&right != this) {
@ -62,22 +56,11 @@ PureFluidPhase& PureFluidPhase::operator=(const PureFluidPhase& right)
return *this;
}
// Duplicator from the %ThermoPhase parent class
/*
* Given a pointer to a %ThermoPhase object, this function will
* duplicate the %ThermoPhase object and all underlying structures.
* This is basically a wrapper around the copy constructor.
*
* @return returns a pointer to a %ThermoPhase
*/
ThermoPhase* PureFluidPhase::duplMyselfAsThermoPhase() const
{
return new PureFluidPhase(*this);
}
PureFluidPhase::~PureFluidPhase()
{
delete m_sub;
@ -175,178 +158,108 @@ pressure() const
setTPXState();
return m_sub->P();
}
//====================================================================================================================
void PureFluidPhase::
setPressure(doublereal p)
void PureFluidPhase::setPressure(doublereal p)
{
Set(tpx::PropertyPair::TP, temperature(), p);
setDensity(1.0/m_sub->v());
}
//====================================================================================================================
void PureFluidPhase::Set(tpx::PropertyPair::type n, double x, double y) const
{
m_sub->Set(n, x, y);
}
//====================================================================================================================
void PureFluidPhase::setTPXState() const
{
Set(tpx::PropertyPair::TV, temperature(), 1.0/density());
}
//====================================================================================================================
doublereal PureFluidPhase::isothermalCompressibility() const
{
return m_sub->isothermalCompressibility();
}
//====================================================================================================================
doublereal PureFluidPhase::thermalExpansionCoeff() const
{
return m_sub->thermalExpansionCoeff();
}
//====================================================================================================================
tpx::Substance& PureFluidPhase::TPX_Substance()
{
return *m_sub;
}
//====================================================================================================================
// Returns an array of partial molar enthalpies for the species
// in the mixture. Units (J/kmol)
/*
* @param hbar Output vector of species partial molar enthalpies.
* Length: m_kk. units are J/kmol.
*/
void PureFluidPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
hbar[0] = enthalpy_mole();
}
//====================================================================================================================
// Returns an array of partial molar entropies of the species in the
// solution. Units: J/kmol/K.
/*
* @param sbar Output vector of species partial molar entropies.
* Length = m_kk. units are J/kmol/K.
*/
void PureFluidPhase::getPartialMolarEntropies(doublereal* sbar) const
{
sbar[0] = entropy_mole();
}
//====================================================================================================================
// Return an array of partial molar internal energies for the
// species in the mixture. Units: J/kmol.
/*
* @param ubar Output vector of species partial molar internal energies.
* Length = m_kk. units are J/kmol.
*/
void PureFluidPhase::getPartialMolarIntEnergies(doublereal* ubar) const
{
ubar[0] = intEnergy_mole();
}
//====================================================================================================================
// Return an array of partial molar heat capacities for the
// species in the mixture. Units: J/kmol/K
/*
* @param cpbar Output vector of species partial molar heat
* capacities at constant pressure.
* Length = m_kk. units are J/kmol/K.
*/
void PureFluidPhase::getPartialMolarCp(doublereal* cpbar) const
{
cpbar[0] = cp_mole();
}
//====================================================================================================================
// Return an array of partial molar volumes for the
// species in the mixture. Units: m^3/kmol.
/*
* @param vbar Output vector of species partial molar volumes.
* Length = m_kk. units are m^3/kmol.
*/
void PureFluidPhase::getPartialMolarVolumes(doublereal* vbar) const
{
vbar[0] = 1.0 / molarDensity();
}
//====================================================================================================================
int PureFluidPhase::standardStateConvention() const
{
return cSS_CONVENTION_TEMPERATURE;
}
//====================================================================================================================
void PureFluidPhase::getActivityConcentrations(doublereal* c) const
{
c[0] = 1.0;
}
//====================================================================================================================
doublereal PureFluidPhase::standardConcentration(size_t k) const
{
return 1.0;
}
//====================================================================================================================
void PureFluidPhase::getActivities(doublereal* a) const
{
a[0] = 1.0;
}
//====================================================================================================================
// Get the array of chemical potentials at unit activity for the species
// at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current
* temperature and pressure of the solution
*
* @param mu Output vector of chemical potentials.
* Length: m_kk.
*/
void PureFluidPhase::getStandardChemPotentials(doublereal* mu) const
{
mu[0] = gibbs_mole();
}
//====================================================================================================================
// Get the nondimensional Enthalpy functions for the species
// at their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*
* @param hrt Output vector of nondimensional standard state enthalpies.
* Length: m_kk.
*/
void PureFluidPhase::getEnthalpy_RT(doublereal* hrt) const
{
doublereal rt = _RT();
doublereal h = enthalpy_mole();
hrt[0] = h / rt;
}
//====================================================================================================================
// Get the array of nondimensional Entropy functions for the
// standard state species at the current <I>T</I> and <I>P</I> of the solution.
/*
* @param sr Output vector of nondimensional standard state entropies.
* Length: m_kk.
*/
void PureFluidPhase::getEntropy_R(doublereal* sr) const
{
doublereal s = entropy_mole();
sr[0] = s / GasConstant;
}
//====================================================================================================================
// Get the nondimensional Gibbs functions for the species
// in their standard states at the current <I>T</I> and <I>P</I> of the solution.
/*
* @param grt Output vector of nondimensional standard state gibbs free energies
* Length: m_kk.
*/
void PureFluidPhase::getGibbs_RT(doublereal* grt) const
{
doublereal rt = _RT();
doublereal g = gibbs_mole();
grt[0] = g / rt;
}
//====================================================================================================================
// Returns the vector of nondimensional enthalpies of the reference state at the current temperature
// of the solution and the reference pressure for the species.
/*
* This base function will throw a CanteraException unless
* it is overwritten in a derived class.
*
* @param hrt Output vector containing the nondimensional reference state enthalpies
* Length: m_kk.
*/
void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const
{
double psave = pressure();
@ -358,13 +271,7 @@ void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const
Set(tpx::PropertyPair::TP, t, psave);
}
//====================================================================================================================
// Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temperature
// of the solution and the reference pressure for the species.
/*
* @param grt Output vector containing the nondimensional reference state
* Gibbs Free energies. Length: m_kk.
*/
void PureFluidPhase::getGibbs_RT_ref(doublereal* grt) const
{
double psave = pressure();
@ -376,27 +283,13 @@ void PureFluidPhase::getGibbs_RT_ref(doublereal* grt) const
grt[0] += log(pref/plow);
Set(tpx::PropertyPair::TP, t, psave);
}
//====================================================================================================================
// Returns the vector of the gibbs function of the reference state at the current temperature
// of the solution and the reference pressure for the species.
/*
* units = J/kmol
*
* @param g Output vector containing the reference state
* Gibbs Free energies. Length: m_kk. Units: J/kmol.
*/
void PureFluidPhase::getGibbs_ref(doublereal* g) const
{
getGibbs_RT_ref(g);
g[0] *= (GasConstant * temperature());
}
//====================================================================================================================
// Returns the vector of nondimensional entropies of the reference state at the current temperature
// of the solution and the reference pressure for each species.
/*
* @param er Output vector containing the nondimensional reference state
* entropies. Length: m_kk.
*/
void PureFluidPhase::getEntropy_R_ref(doublereal* er) const
{
double psave = pressure();
@ -408,62 +301,56 @@ void PureFluidPhase::getEntropy_R_ref(doublereal* er) const
er[0] -= log(pref/plow);
Set(tpx::PropertyPair::TP, t, psave);
}
//====================================================================================================================
// critical temperature
doublereal PureFluidPhase::critTemperature() const
{
return m_sub->Tcrit();
}
//====================================================================================================================
/// critical pressure
doublereal PureFluidPhase::critPressure() const
{
return m_sub->Pcrit();
}
//====================================================================================================================
/// critical density
doublereal PureFluidPhase::critDensity() const
{
return 1.0/m_sub->Vcrit();
}
//====================================================================================================================
/// saturation temperature
doublereal PureFluidPhase::satTemperature(doublereal p) const
{
doublereal ts = m_sub->Tsat(p);
return ts;
}
//====================================================================================================================
void PureFluidPhase::setState_HP(doublereal h, doublereal p,
doublereal tol)
{
Set(tpx::PropertyPair::HP, h, p);
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
}
//====================================================================================================================
void PureFluidPhase::setState_UV(doublereal u, doublereal v,
doublereal tol)
{
Set(tpx::PropertyPair::UV, u, v);
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
}
//====================================================================================================================
void PureFluidPhase::setState_SV(doublereal s, doublereal v,
doublereal tol)
{
Set(tpx::PropertyPair::SV, s, v);
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
}
//====================================================================================================================
void PureFluidPhase::setState_SP(doublereal s, doublereal p,
doublereal tol)
{
Set(tpx::PropertyPair::SP, s, p);
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
}
//====================================================================================================================
// saturation pressure
doublereal PureFluidPhase::satPressure(doublereal t) const
{
doublereal vsv = m_sub->v();
@ -471,13 +358,13 @@ doublereal PureFluidPhase::satPressure(doublereal t) const
doublereal ps = m_sub->Ps();
return ps;
}
//====================================================================================================================
doublereal PureFluidPhase::vaporFraction() const
{
setTPXState();
return m_sub->x();
}
//====================================================================================================================
void PureFluidPhase::setState_Tsat(doublereal t, doublereal x)
{
setTemperature(t);
@ -485,7 +372,7 @@ void PureFluidPhase::setState_Tsat(doublereal t, doublereal x)
Set(tpx::PropertyPair::TX, t, x);
setDensity(1.0/m_sub->v());
}
//====================================================================================================================
void PureFluidPhase::setState_Psat(doublereal p, doublereal x)
{
setTPXState();
@ -494,10 +381,6 @@ void PureFluidPhase::setState_Psat(doublereal p, doublereal x)
setDensity(1.0/m_sub->v());
}
//====================================================================================================================
/**
* Format a summary of the mixture state for output.
*/
std::string PureFluidPhase::report(bool show_thermo) const
{
char p[800];

347
src/thermo/RedlichKwongMFTP.cpp
View file

@ -28,10 +28,6 @@ const doublereal RedlichKwongMFTP::omega_a = 4.27480233540E-01;
const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
//====================================================================================================================
/*
* Default constructor
*/
RedlichKwongMFTP::RedlichKwongMFTP() :
MixtureFugacityTP(),
m_standardMixingRules(0),
@ -56,7 +52,7 @@ RedlichKwongMFTP::RedlichKwongMFTP() :
Vroot_[1] = 0.0;
Vroot_[2] = 0.0;
}
//====================================================================================================================
RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, std::string id) :
MixtureFugacityTP(),
m_standardMixingRules(0),
@ -91,7 +87,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, std::string id) :
}
importPhase(*xphase, this);
}
//====================================================================================================================
RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id) :
MixtureFugacityTP(),
m_standardMixingRules(0),
@ -122,7 +118,6 @@ RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id
importPhase(*xphase, this);
}
//====================================================================================================================
RedlichKwongMFTP::RedlichKwongMFTP(int testProb) :
MixtureFugacityTP(),
m_standardMixingRules(0),
@ -161,16 +156,7 @@ RedlichKwongMFTP::RedlichKwongMFTP(int testProb) :
}
importPhase(*xphase, this);
}
//====================================================================================================================
/*
* Copy Constructor:
*
* Note this stuff will not work until the underlying phase
* has a working copy constructor.
*
* The copy constructor just calls the assignment operator
* to do the heavy lifting.
*/
RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP& b) :
MixtureFugacityTP(),
m_standardMixingRules(0),
@ -194,13 +180,6 @@ RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP& b) :
*this = b;
}
//====================================================================================================================
/*
* operator=()
*
* Note this stuff will not work until the underlying phase
* has a working assignment operator
*/
RedlichKwongMFTP& RedlichKwongMFTP::
operator=(const RedlichKwongMFTP& b)
{
@ -237,35 +216,25 @@ operator=(const RedlichKwongMFTP& b)
}
return *this;
}
//====================================================================================================================
/*
* ~RedlichKwongMFTP(): (virtual)
*
*/
RedlichKwongMFTP::~RedlichKwongMFTP()
{
}
//====================================================================================================================
/*
* Duplication function.
* This calls the copy constructor for this object.
*/
ThermoPhase* RedlichKwongMFTP::duplMyselfAsThermoPhase() const
{
return new RedlichKwongMFTP(*this);
}
//====================================================================================================================
int RedlichKwongMFTP::eosType() const
{
return cRedlichKwongMFTP;
}
//====================================================================================================================
/*
* ------------Molar Thermodynamic Properties -------------------------
*/
//====================================================================================================================
// Molar enthalpy. Units: J/kmol.
doublereal RedlichKwongMFTP::enthalpy_mole() const
{
_updateReferenceStateThermo();
@ -274,16 +243,14 @@ doublereal RedlichKwongMFTP::enthalpy_mole() const
doublereal h_nonideal = hresid();
return (h_ideal + h_nonideal);
}
//====================================================================================================================
// Molar internal energy. Units: J/kmol.
doublereal RedlichKwongMFTP::intEnergy_mole() const
{
doublereal p0 = pressure();
doublereal md = molarDensity();
return (enthalpy_mole() - p0 / md);
}
//====================================================================================================================
// Molar entropy. Units: J/kmol/K.
doublereal RedlichKwongMFTP::entropy_mole() const
{
_updateReferenceStateThermo();
@ -292,14 +259,12 @@ doublereal RedlichKwongMFTP::entropy_mole() const
doublereal sr_nonideal = sresid();
return (sr_ideal + sr_nonideal);
}
//====================================================================================================================
// Molar Gibbs function. Units: J/kmol.
doublereal RedlichKwongMFTP::gibbs_mole() const
{
return enthalpy_mole() - temperature() * entropy_mole();
}
//====================================================================================================================
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
doublereal RedlichKwongMFTP::cp_mole() const
{
_updateReferenceStateThermo();
@ -316,29 +281,15 @@ doublereal RedlichKwongMFTP::cp_mole() const
double cp = dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
return cp;
}
//====================================================================================================================
/// Molar heat capacity at constant volume. Units: J/kmol/K.
doublereal RedlichKwongMFTP::cv_mole() const
{
throw CanteraError("", "unimplemented");
return cp_mole() - GasConstant;
}
//====================================================================================================================
// Return the thermodynamic pressure (Pa).
/*
* Since the mass density, temperature, and mass fractions are stored,
* this method uses these values to implement the
* mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$.
*
* \f[
* P = \frac{RT}{v-b_{mix}} - \frac{a_{mix}}{T^{0.5} v \left( v + b_{mix} \right) }
* \f]
*
*/
doublereal RedlichKwongMFTP::pressure() const
{
#ifdef DEBUG_MODE
_updateReferenceStateThermo();
@ -354,10 +305,9 @@ doublereal RedlichKwongMFTP::pressure() const
throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe");
}
#endif
return m_Pcurrent;
}
//====================================================================================================================
void RedlichKwongMFTP::calcDensity()
{
/*
@ -373,58 +323,52 @@ void RedlichKwongMFTP::calcDensity()
*/
double dens = 1.0/invDens;
Phase::setDensity(dens);
}
//====================================================================================================================
void RedlichKwongMFTP::setTemperature(const doublereal temp)
{
Phase::setTemperature(temp);
_updateReferenceStateThermo();
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::setMassFractions(const doublereal* const x)
{
MixtureFugacityTP::setMassFractions(x);
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::setMassFractions_NoNorm(const doublereal* const x)
{
MixtureFugacityTP::setMassFractions_NoNorm(x);
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::setMoleFractions(const doublereal* const x)
{
MixtureFugacityTP::setMoleFractions(x);
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::setMoleFractions_NoNorm(const doublereal* const x)
{
MixtureFugacityTP::setMoleFractions(x);
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::setConcentrations(const doublereal* const c)
{
MixtureFugacityTP::setConcentrations(c);
updateAB();
}
//====================================================================================================================
doublereal RedlichKwongMFTP::isothermalCompressibility() const
{
throw CanteraError("RedlichKwongMFTP::isothermalCompressibility() ",
"not implemented");
return 0.0;
}
//====================================================================================================================
void RedlichKwongMFTP::getActivityConcentrations(doublereal* c) const
{
getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes));
@ -432,59 +376,20 @@ void RedlichKwongMFTP::getActivityConcentrations(doublereal* c) const
c[k] = moleFraction(k) / m_partialMolarVolumes[k];
}
}
//====================================================================================================================
/*
* Returns the standard concentration \f$ C^0_k \f$, which is used to normalize
* the generalized concentration.
*/
doublereal RedlichKwongMFTP::standardConcentration(size_t k) const
{
getStandardVolumes(DATA_PTR(m_tmpV));
return 1.0 / m_tmpV[k];
}
//====================================================================================================================
/*
* Returns the natural logarithm of the standard
* concentration of the kth species
*/
doublereal RedlichKwongMFTP::logStandardConc(size_t k) const
{
double c = standardConcentration(k);
double lc = std::log(c);
return lc;
}
//====================================================================================================================
/*
*
* getUnitsStandardConcentration()
*
* Returns the units of the standard and general concentrations
* Note they have the same units, as their divisor is
* defined to be equal to the activity of the kth species
* in the solution, which is unitless.
*
* This routine is used in print out applications where the
* units are needed. Usually, MKS units are assumed throughout
* the program and in the XML input files.
*
* uA[0] = kmol units - default = 1
* uA[1] = m units - default = -nDim(), the number of spatial
* dimensions in the Phase class.
* uA[2] = kg units - default = 0;
* uA[3] = Pa(pressure) units - default = 0;
* uA[4] = Temperature units - default = 0;
* uA[5] = time units - default = 0
*
* For EOS types other than cIdealSolidSolnPhase1, the default
* kmol/m3 holds for standard concentration units. For
* cIdealSolidSolnPhase0 type, the standard concentration is
* unitless.
*/
void RedlichKwongMFTP::getUnitsStandardConc(double* uA, int, int sizeUA) const
{
//int eos = eosType();
@ -512,14 +417,6 @@ void RedlichKwongMFTP::getUnitsStandardConc(double* uA, int, int sizeUA) const
}
//====================================================================================================================
//! Get the array of non-dimensional activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* For ideal gases, the activity coefficients are all equal to one.
*
* @param ac Output vector of activity coefficients. Length: m_kk.
*/
void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const
{
doublereal TKelvin = temperature();
@ -551,20 +448,11 @@ void RedlichKwongMFTP::getActivityCoefficients(doublereal* ac) const
ac[k] = exp(ac[k]/rt);
}
}
//====================================================================================================================
/*
* ---- Partial Molar Properties of the Solution -----------------
*/
//====================================================================================================================
/*
* Get the array of non-dimensional species chemical potentials
* These are partial molar Gibbs free energies.
* \f$ \mu_k / \hat R T \f$.
* Units: unitless
*
* We close the loop on this function, here, calling
* getChemPotentials() and then dividing by RT.
*/
void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
{
getChemPotentials(muRT);
@ -573,7 +461,7 @@ void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
muRT[k] *= invRT;
}
}
//====================================================================================================================
void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
{
getGibbs_ref(mu);
@ -610,7 +498,7 @@ void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
);
}
}
//====================================================================================================================
void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
@ -669,7 +557,7 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(doublereal* hbar) const
}
}
//====================================================================================================================
void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
{
getEntropy_R_ref(sbar);
@ -726,21 +614,21 @@ void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
}
}
//====================================================================================================================
void RedlichKwongMFTP::getPartialMolarIntEnergies(doublereal* ubar) const
{
getIntEnergy_RT(ubar);
doublereal rt = GasConstant * temperature();
scale(ubar, ubar+m_kk, ubar, rt);
}
//====================================================================================================================
void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const
{
getCp_R(cpbar);
doublereal r = GasConstant;
scale(cpbar, cpbar+m_kk, cpbar, r);
}
//====================================================================================================================
void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
{
// getStandardVolumes(vbar);
@ -786,7 +674,7 @@ void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
}
}
//====================================================================================================================
doublereal RedlichKwongMFTP::critTemperature() const
{
double pc, tc, vc;
@ -802,7 +690,7 @@ doublereal RedlichKwongMFTP::critTemperature() const
calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
return tc;
}
//====================================================================================================================
doublereal RedlichKwongMFTP::critPressure() const
{
double pc, tc, vc;
@ -819,7 +707,7 @@ doublereal RedlichKwongMFTP::critPressure() const
return pc;
}
//====================================================================================================================
doublereal RedlichKwongMFTP::critDensity() const
{
double pc, tc, vc;
@ -837,26 +725,13 @@ doublereal RedlichKwongMFTP::critDensity() const
double mmw = meanMolecularWeight();
return mmw / vc;
}
//====================================================================================================================
/*
* ----- Thermodynamic Values for the Species Reference States ----
*/
//====================================================================================================================
/*
* Perform initializations after all species have been
* added.
*/
void RedlichKwongMFTP::initThermo()
{
initLengths();
MixtureFugacityTP::initThermo();
}
//====================================================================================================================
void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
{
double tmp, tmp2;
@ -891,11 +766,7 @@ void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
// set state
setState_PX(pres, &m_pp[0]);
}
//====================================================================================================================
/*
* Initialize the internal lengths.
* (this is not a virtual function)
*/
void RedlichKwongMFTP::initLengths()
{
@ -916,24 +787,7 @@ void RedlichKwongMFTP::initLengths()
m_partialMolarVolumes.resize(m_kk, 0.0);
dpdni_.resize(m_kk, 0.0);
}
//====================================================================================================================
/*
* Import and initialize a ThermoPhase object
*
* param phaseNode This object must be the phase node of a
* complete XML tree
* description of the phase, including all of the
* species data. In other words while "phase" must
* point to an XML phase object, it must have
* sibling nodes "speciesData" that describe
* the species in the phase.
* param id ID of the phase. If nonnull, a check is done
* to see if phaseNode is pointing to the phase
* with the correct id.
*
* This routine initializes the lengths in the current object and
* then calls the parent routine.
*/
void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
{
RedlichKwongMFTP::initLengths();
@ -1018,7 +872,6 @@ void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
MixtureFugacityTP::initThermoXML(phaseNode, id);
}
//====================================================================================================================
void RedlichKwongMFTP::readXMLPureFluid(XML_Node& pureFluidParam)
{
@ -1083,7 +936,7 @@ void RedlichKwongMFTP::readXMLPureFluid(XML_Node& pureFluidParam)
}
}
}
//====================================================================================================================
void RedlichKwongMFTP::applyStandardMixingRules()
{
int nParam = 2;
@ -1100,7 +953,6 @@ void RedlichKwongMFTP::applyStandardMixingRules()
}
}
}
//====================================================================================================================
void RedlichKwongMFTP::readXMLCrossFluid(XML_Node& CrossFluidParam)
{
@ -1168,22 +1020,13 @@ void RedlichKwongMFTP::readXMLCrossFluid(XML_Node& CrossFluidParam)
}
}
}
//====================================================================================================================
void RedlichKwongMFTP::setParametersFromXML(const XML_Node& thermoNode)
{
MixtureFugacityTP::setParametersFromXML(thermoNode);
std::string model = thermoNode["model"];
}
//====================================================================================================================
// Calculate the deviation terms for the total entropy of the mixture from the
// ideal gas mixture
/*
* Here we use the current state conditions
*
* @return Returns the change in entropy in units of J kmol-1 K-1.
*/
doublereal RedlichKwongMFTP::sresid() const
{
// note this agrees with tpx
@ -1200,14 +1043,7 @@ doublereal RedlichKwongMFTP::sresid() const
double sp = GasConstant * sresid_mol_R;
return sp;
}
//====================================================================================================================
// Calculate the deviation terms for the total enthalpy of the mixture from the
// ideal gas mixture
/*
* Here we use the current state conditions
*
* @return Returns the change in entropy in units of J kmol-1.
*/
doublereal RedlichKwongMFTP::hresid() const
{
// note this agrees with tpx
@ -1223,21 +1059,7 @@ doublereal RedlichKwongMFTP::hresid() const
double hresid_mol = GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
return hresid_mol;
}
//====================================================================================================================
// Estimate for the molar volume of the liquid
/*
* Note: this is only used as a starting guess for later routines that actually calculate an
* accurate value for the liquid molar volume.
* This routine doesn't change the state of the system.
*
* @param TKelvin temperature in kelvin
* @param pres Pressure in Pa. This is used as an initial guess. If the routine
* needs to change the pressure to find a stable liquid state, the
* new pressure is returned in this variable.
*
* @return Returns the estimate of the liquid volume. If the liquid can't be found, this
* routine returns -1.
*/
doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
{
double v = m_b_current * 1.1;
@ -1282,31 +1104,7 @@ doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGu
//printf (" RedlichKwongMFTP::liquidVolEst %g %g converged in %d its\n", TKelvin, pres, i);
return v;
}
//====================================================================================================================
// Calculates the density given the temperature and the pressure and a guess at the density.
/*
* Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this.
* This is protected because it is called during setState_TP() routines. Infinite loops would result
* if it were not protected.
*
* -> why is this not const?
*
* parameters:
* @param TKelvin Temperature in Kelvin
* @param pressure Pressure in Pascals (Newton/m**2)
* @param phaseReqested int representing the phase whose density we are requesting. If we put
* a gas or liquid phase here, we will attempt to find a volume in that
* part of the volume space, only, in this routine. A value of FLUID_UNDEFINED
* means that we will accept anything.
*
* @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there
* is no guessed density
*
*
* @return We return the density of the fluid at the requested phase. If we have not found any
* acceptable density we return a -1. If we have found an acceptable density at a
* different phase, we return a -2.
*/
doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
{
@ -1382,12 +1180,7 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa,
densBase = mmw / molarVolLast;
return densBase;
}
//====================================================================================================================
// Return the value of the density at the liquid spinodal point (on the liquid side)
// for the current temperature.
/*
* @return returns the density with units of kg m-3
*/
doublereal RedlichKwongMFTP::densSpinodalLiquid() const
{
if (NSolns_ != 3) {
@ -1412,12 +1205,7 @@ doublereal RedlichKwongMFTP::densSpinodalLiquid() const
doublereal rho = mmw / vbest;
return rho;
}
//====================================================================================================================
// Return the value of the density at the gas spinodal point (on the gas side)
// for the current temperature.
/*
* @return returns the density with units of kg m-3
*/
doublereal RedlichKwongMFTP::densSpinodalGas() const
{
if (NSolns_ != 3) {
@ -1442,16 +1230,7 @@ doublereal RedlichKwongMFTP::densSpinodalGas() const
doublereal rho = mmw / vbest;
return rho;
}
//====================================================================================================================
// Calculate the pressure given the temperature and the molar volume
/*
* Calculate the pressure given the temperature and the molar volume
*
* @param TKelvin temperature in kelvin
* @param molarVol molar volume ( m3/kmol)
*
* @return Returns the pressure.
*/
doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
{
doublereal sqt = sqrt(TKelvin);
@ -1459,18 +1238,7 @@ doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVo
- m_a_current / (sqt * molarVol * (molarVol + m_b_current));
return pres;
}
//====================================================================================================================
// Calculate the pressure and the pressure derivative given the temperature and the molar volume
/*
* Temperature and mole number are held constant
*
* @param TKelvin temperature in kelvin
* @param molarVol molar volume ( m3/kmol)
*
* @param presCalc Returns the pressure.
*
* @return Returns the derivative of the pressure wrt the molar volume
*/
doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
{
doublereal sqt = sqrt(TKelvin);
@ -1483,7 +1251,6 @@ doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol,
+ m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
return dpdv;
}
//====================================================================================================================
void RedlichKwongMFTP::pressureDerivatives() const
{
@ -1501,12 +1268,12 @@ void RedlichKwongMFTP::pressureDerivatives() const
dpdT_ = (GasConstant / (vmb) - fac / (sqt * mv * vpb));
}
//====================================================================================================================
void RedlichKwongMFTP::updateMixingExpressions()
{
updateAB();
}
//====================================================================================================================
void RedlichKwongMFTP::updateAB()
{
double temp = temperature();
@ -1528,7 +1295,7 @@ void RedlichKwongMFTP::updateAB()
}
}
}
//====================================================================================================================
void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
{
bCalc = 0.0;
@ -1553,7 +1320,7 @@ void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublerea
}
}
}
//====================================================================================================================
doublereal RedlichKwongMFTP::da_dt() const
{
@ -1568,7 +1335,7 @@ doublereal RedlichKwongMFTP::da_dt() const
}
return dadT;
}
//====================================================================================================================
void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
doublereal& pc, doublereal& tc, doublereal& vc) const
{
@ -1612,18 +1379,6 @@ void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, double
vc = omega_vc * GasConstant * tc / pc;
}
//====================================================================================================================
// Solve the cubic equation of state
/*
* The R-K equation of state may be solved via the following formula
*
* V**3 - V**2(RT/P) - V(RTb/P - a/(P T**.5) + b*b) - (a b / (P T**.5)) = 0
*
* Returns the number of solutions found. If it only finds the liquid branch solution, it will return a -1 or a -2
* instead of 1 or 2. If it returns 0, then there is an error.
*
*/
int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
doublereal Vroot[3]) const
{