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

This update adds in a Multicomponent Real Gas capability to Cantera.

Browse source 
The Object is called RedlichKwongMFTP, and it implements a multicomponent real gas approximation.
Most of the functionality has been verified against sample problems involving the CO2-H2O system.
This commit is contained in:
Harry Moffat 2011-07-23 00:30:35 +00:00
parent 1982d0c928
commit ff7e66e5a3
22 changed files with 7390 additions and 34 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

34
Cantera/src/base/checkFinite.cpp
View file

@ -50,11 +50,11 @@ namespace mdp {
void checkFinite(const double tmp) throw(std::range_error) {
if (_finite(tmp)) {
if(_isnan(tmp)) {
printf("ERROR: we have encountered a nan!\n");
printf("checkFinite() ERROR: we have encountered a nan!\n");
} else if (_fpclass(tmp) == _FPCLASS_PINF) {
printf("ERROR: we have encountered a pos inf!\n");
printf("checkFinite() ERROR: we have encountered a pos inf!\n");
} else {
printf("ERROR: we have encountered a neg inf!\n");
printf("checkFinite() ERROR: we have encountered a neg inf!\n");
}
const std::string s = "checkFinite()";
throw std::range_error(s);
@ -64,11 +64,11 @@ namespace mdp {
void checkFinite(const double tmp) throw(std::range_error) {
if (! finite(tmp)) {
if(isnan(tmp)) {
printf("ERROR: we have encountered a nan!\n");
printf("checkFinite() ERROR: we have encountered a nan!\n");
} else if (isinf(tmp) == 1) {
printf("ERROR: we have encountered a pos inf!\n");
printf("checkFinite() ERROR: we have encountered a pos inf!\n");
} else {
printf("ERROR: we have encountered a neg inf!\n");
printf("checkFinite() ERROR: we have encountered a neg inf!\n");
}
const std::string s = "checkFinite()";
throw std::range_error(s);
@ -102,7 +102,7 @@ namespace mdp {
checkFinite(tmp);
if (fabs(tmp) >= trigger) {
char sbuf[64];
sprintf(sbuf, "checkMagnitude: Trigger %g exceeded: %g\n", trigger,
sprintf(sbuf, "checkMagnitude() ERROR: Trigger %g exceeded: %g\n", trigger,
tmp);
throw std::range_error(sbuf);
}
@ -119,16 +119,16 @@ namespace mdp {
void checkZeroFinite(const double tmp) throw(std::range_error) {
if ((tmp == 0.0) || (! _finite(tmp))) {
if (tmp == 0.0) {
printf("ERROR: we have encountered a zero!\n");
printf("checkZeroFinite() ERROR: we have encountered a zero!\n");
} else if(_isnan(tmp)) {
printf("ERROR: we have encountered a nan!\n");
printf("checkZeroFinite() ERROR: we have encountered a nan!\n");
} else if (_fpclass(tmp) == _FPCLASS_PINF) {
printf("ERROR: we have encountered a pos inf!\n");
printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n");
} else {
printf("ERROR: we have encountered a neg inf!\n");
printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n");
}
char sbuf[64];
sprintf(sbuf, "checkZeroFinite: zero or indef exceeded: %g\n",
sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n",
tmp);
throw std::range_error(sbuf);
}
@ -137,16 +137,16 @@ void checkZeroFinite(const double tmp) throw(std::range_error) {
void checkZeroFinite(const double tmp) throw(std::range_error) {
if ((tmp == 0.0) || (! finite(tmp))) {
if (tmp == 0.0) {
printf("ERROR: we have encountered a zero!\n");
printf("checkZeroFinite() ERROR: we have encountered a zero!\n");
} else if(isnan(tmp)) {
printf("ERROR: we have encountered a nan!\n");
printf("checkZeroFinite() ERROR: we have encountered a nan!\n");
} else if (isinf(tmp) == 1) {
printf("ERROR: we have encountered a pos inf!\n");
printf("checkZeroFinite() ERROR: we have encountered a pos inf!\n");
} else {
printf("ERROR: we have encountered a neg inf!\n");
printf("checkZeroFinite() ERROR: we have encountered a neg inf!\n");
}
char sbuf[64];
sprintf(sbuf, "checkZeroFinite: zero or indef exceeded: %g\n",
sprintf(sbuf, "checkZeroFinite() ERROR: zero or indef exceeded: %g\n",
tmp);
throw std::range_error(sbuf);
}

1
Cantera/src/numerics/RootFind.h
View file

@ -137,7 +137,6 @@ namespace Cantera {
*/
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal &funcTargetValue, doublereal *xbest);
//! Return the function value
/*!
* This routine evaluates the following equation.

6
Cantera/src/thermo/Makefile.in
View file

@ -40,7 +40,7 @@ THERMO_OBJ = State.o Elements.o Constituents.o Phase.o \
ThermoFactory.o phasereport.o SpeciesThermoInterpType.o \
VPSSMgr.o VPSSMgrFactory.o VPSSMgr_General.o IdealSolnGasVPSS.o \
VPSSMgr_IdealGas.o VPSSMgr_ConstVol.o PDSS_ConstVol.o PDSS_IdealGas.o \
PDSS_SSVol.o @phase_object_files@
PDSS_SSVol.o MixtureFugacityTP.o RedlichKwongMFTP.o @phase_object_files@
THERMO_H = State.h Elements.h Constituents.h Phase.h mix_defs.h \
ThermoPhase.h IdealGasPhase.h ConstDensityThermo.h \
@ -54,7 +54,7 @@ THERMO_H = State.h Elements.h Constituents.h Phase.h mix_defs.h \
EdgePhase.h \
VPSSMgr.h VPSSMgrFactory.h VPSSMgr_General.h IdealSolnGasVPSS.h \
VPSSMgr_IdealGas.h VPSSMgr_ConstVol.h PDSS_ConstVol.h PDSS_IdealGas.h \
PDSS_SSVol.h @phase_header_files@
PDSS_SSVol.h MixtureFugacityTP.h RedlichKwongMFTP.h @phase_header_files@
# Extended Cantera Thermodynamics Object Files
@ -91,7 +91,7 @@ CATHERMO_OBJ = $(THERMO_OBJ) $(ELECTRO_OBJ) $(ISSP_OBJ)
CATHERMO_H = $(THERMO_H) $(ELECTRO_H) $(ISSP_H)
CXX_INCLUDES = -I../base @CXX_INCLUDES@
CXX_INCLUDES = -I../base -I../numerics @CXX_INCLUDES@
LIB = @buildlib@/libthermo.a
DEPENDS = $(CATHERMO_OBJ:.o=.d)

1214
Cantera/src/thermo/MixedSolventElectrolyte.cpp Normal file
View file

File diff suppressed because it is too large Load diff

1025
Cantera/src/thermo/MixedSolventElectrolyte.h Normal file
View file

File diff suppressed because it is too large Load diff

1367
Cantera/src/thermo/MixtureFugacityTP.cpp Normal file
View file

File diff suppressed because it is too large Load diff

968
Cantera/src/thermo/MixtureFugacityTP.h Normal file
View file

@ -0,0 +1,968 @@
/**
* @file MixtureFugacityTP.h
* 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).
*
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Date: 2010-08-04 10:06:14 -0600 (Wed, 04 Aug 2010) $
* $Revision: 547 $
*/
#ifndef CT_MIXTUREFUGACITYTP_H
#define CT_MIXTUREFUGACITYTP_H
#include "ThermoPhase.h"
#include "VPSSMgr.h"
#include "ResidEval.h"
namespace Cantera {
class XML_Node;
class PDSS;
//! Various states of the Fugacity object. In general there can be multiple liquid
//! objects for a single phase identified with each species.
#define FLUID_UNSTABLE -4
#define FLUID_UNDEFINED -3
#define FLUID_SUPERCRIT -2
#define FLUID_GAS -1
#define FLUID_LIQUID_0 0
#define FLUID_LIQUID_1 1
#define FLUID_LIQUID_2 2
#define FLUID_LIQUID_3 3
#define FLUID_LIQUID_4 4
#define FLUID_LIQUID_5 5
#define FLUID_LIQUID_6 6
#define FLUID_LIQUID_7 7
#define FLUID_LIQUID_8 8
#define FLUID_LIQUID_9 9
/**
* @ingroup thermoprops
*
* This is a filter class for ThermoPhase that implements some prepatory
* steps for efficiently handling mixture of gases that whose standard states
* are defined as ideal gases, but which describe also non-ideal solutions.
* In addition a multicomponent liquid phase below the critical temperature of the
* mixture is also allowed. The main subclass will be a mixture Redlich-kwong class.
*
* Several concepts are introduced. The first concept is there are temporary
* variables for holding the species standard state values
* of Cp, H, S, G, and V at the last temperature and pressure called. These functions are not recalculated
* if a new call is made using the previous temperature and pressure.
*
* The other concept is that the current state of the mixture is tracked.
* The state variable is either GAS, LIQUID, or SUPERCRIT fluid. Additionally,
* the variable LiquidContent is used and may vary between 0 and 1.
*
* To support the above functionality, pressure and temperature variables,
* m_Plast_ss and m_Tlast_ss, are kept which store the last pressure and temperature
* used in the evaluation of standard state properties.
*
* Typically, only one liquid phase is allowed to be formed within these classes.
* Additionally, there is an inherent contradiction between three phase models and
* the ThermoPhase class. The ThermoPhase class is really only meant to represent a
* single instanteation of a phase. The three phase models may be in equilibrium with
* multiple phases of the fluid in equilibrium with each other. This has yet to be resolved.
*
* This class is usually used for non-ideal gases.
*
*
* @nosubgrouping
*/
class MixtureFugacityTP : public ThermoPhase {
public:
/*!
*
* @name Constructors and Duplicators for %MixtureFugacityTP
*
*/
//! Constructor.
MixtureFugacityTP();
//! Copy Constructor.
/*!
* @param b Object to be copied
*/
MixtureFugacityTP(const MixtureFugacityTP &b);
//! Assignment operator
/*!
* @param b Object to be copied
*/
MixtureFugacityTP& operator=(const MixtureFugacityTP &b);
//! Destructor.
virtual ~MixtureFugacityTP();
//! Duplication routine
/*!
* @return Returns a duplication
*/
virtual ThermoPhase *duplMyselfAsThermoPhase() const;
//@}
/**
* @name Utilities (MixtureFugacityTP)
*/
//@{
/**
* Equation of state type flag. The base class returns
* zero. Subclasses should define this to return a unique
* non-zero value. Constants defined for this purpose are
* listed in mix_defs.h.
*/
virtual int eosType() const { return 0; }
//! 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
*/
virtual int standardStateConvention() const;
//! 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.
*/
virtual void setForcedSolutionBranch(int 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.
*/
virtual int forcedSolutionBranch() const;
//! 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.
*/
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.
* This quantity is to be used in conjunction with derivatives of
* that concentration-like variable when the derivative of the chemical
* potential is taken.
*
* units = dimensionless
*
* @param dlnActCoeffdlnN_diag Output vector of derivatives of the
* log Activity Coefficients. length = m_kk
*/
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const {
err("getdlnActCoeffdlnN_diag");
}
//@}
/// @name Partial Molar Properties of the Solution (MixtureFugacityTP)
//@{
//! 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. No need for child
* classes to handle.
*
* @param mu Output vector of non-dimensional species chemical potentials
* Length: m_kk.
*/
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(),
* which must be overloaded in inherited objects.
* The values are cached within this object, and are not recalculated unless
* the temperature or pressure changes.
*/
//@{
//! Get the array of chemical potentials at unit activity.
/*!
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
* \f$. The values are evaluated at the current temperature and pressure.
*
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* @param mu Output vector of standard state chemical potentials.
* length = m_kk. units are J / kmol.
*/
virtual void getStandardChemPotentials(doublereal* mu) const;
//! 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.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* @param hrt Output vector of standard state enthalpies.
* length = m_kk. units are unitless.
*/
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.
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* @param sr Output vector of nondimensional standard state
* entropies. length = m_kk.
*/
virtual void getEntropy_R(doublereal* sr) const;
//! Get the nondimensional Gibbs functions for the species
//! at their standard states of solution at the current T and P of the solution.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* @param grt Output vector of nondimensional standard state
* Gibbs free energies. length = m_kk.
*/
virtual void getGibbs_RT(doublereal* grt) const;
//! Get the nondimensional Gibbs functions for the standard
//! state of the species at the current T and P.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* @param gpure Output vector of standard state
* Gibbs free energies. length = m_kk.
* units are J/kmol.
*
* @todo This could be eliminated. It doesn't fit into the current
* naming convention.
*/
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.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution.
*
* \f[
* u^{ss}_k(T,P) = h^{ss}_k(T) - P * V^{ss}_k
* \f]
*
* @param urt Output vector of nondimensional standard state
* internal energies. length = m_kk.
*/
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.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure of the solution.
*
* @param cpr Output vector containing the
* the nondimensional Heat Capacities at constant
* pressure for the standard state of the species.
* Length: m_kk.
*/
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.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure of the solution.
*
* units = m^3 / kmol
*
* @param vol Output vector of species volumes. length = m_kk.
* units = m^3 / kmol
*/
virtual void getStandardVolumes(doublereal *vol) const;
//! Set the temperature of the phase
/*!
* Currently this passes down to setState_TP(). It does not
* make sense to calculate the standard state without first
* setting T and P.
*
* @param temp Temperature (kelvin)
*/
virtual void setTemperature(const doublereal temp);
//! Set the internally storred pressure (Pa) at constant
//! temperature and composition
/*!
* Currently this passes down to setState_TP(). It does not
* make sense to calculate the standard state without first
* setting T and P.
*
* @param p input Pressure (Pa)
*/
virtual void setPressure(doublereal p);
protected:
/**
* 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
* 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();
public:
//! Set the temperature and pressure at the same time
/*!
* Note this function triggers a reevalulation of the standard
* state quantities.
*
* @param T temperature (kelvin)
* @param pres pressure (pascal)
*/
virtual void setState_TP(doublereal T, doublereal pres);
//! Set the internally storred temperature (K) and density (kg/m^3)
/*!
* @param t Temperature in kelvin
* @param rho Density (kg/m^3)
*/
virtual void setState_TR(doublereal T, doublereal rho);
//! 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.
*/
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.
/*!
* @param y Array of unnormalized mass fraction values (input).
* Must have a length greater than or equal to the number of species.
*/
virtual void setMassFractions(const doublereal* const y);
//!Set the mass 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 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.
/*!
* @param x Array of unnormalized mole fraction values (input).
* Must have a length greater than or equal to the number of species.
*/
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 ofequations.
*
* @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
* units. For volumetric phases c[k] is the
* concentration of the kth species in kmol/m3.
* For surface phases, c[k] is the concentration
* in kmol/m2. The length of the vector is the number
* of species in the phase.
*/
virtual void setConcentrations(const doublereal* const c);
protected:
void setMoleFractions_NoState(const doublereal* const x);
public:
//! Returns the current pressure of the phase
/*!
* The pressure is an independent variable in this phase. Its current value
* is storred in the object MixtureFugacityTP.
*
* @return return the pressure in pascals.
*/
doublereal pressure() const {
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
* 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.
*
*/
virtual void _updateReferenceStateThermo() const;
public:
//@}
/// @name Thermodynamic Values for the Species Reference States (MixtureFugacityTP)
/*!
* There are also temporary
* variables for holding the species reference-state values of Cp, H, S, and V at the
* last temperature and reference pressure called. These functions are not recalculated
* if a new call is made using the previous temperature.
* All calculations are done within the routine _updateRefStateThermo().
*/
//@{
//! Returns the vector of nondimensional
//! enthalpies of the reference state at the current temperature
//! of the solution and the reference pressure for the species.
/*!
* @param hrt Output vector contains the nondimensional enthalpies
* of the reference state of the species
* length = m_kk, units = dimensionless.
*/
virtual void getEnthalpy_RT_ref(doublereal *hrt) const;
#ifdef H298MODIFY_CAPABILITY
//! Modify the value of the 298 K Heat of Formation of the standard state 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 Index of the species
* @param Hf298New Specify the new value of the Heat of Formation at 298K and 1 bar.
* units = J/kmol.
*/
void modifyOneHf298SS(const int k, const doublereal Hf298New);
#endif
//! 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 contains the nondimensional Gibbs free energies
* of the reference state of the species
* length = m_kk, units = dimensionless.
*/
virtual void getGibbs_RT_ref(doublereal *grt) const;
protected:
//! 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.
/*!
* @return Output vector contains the nondimensional Gibbs free energies
* of the reference state of the species
* length = m_kk, units = dimensionless.
*/
const vector_fp & gibbs_RT_ref() const;
public:
/*!
* 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 contain the Gibbs free energies
* of the reference state of the species
* length = m_kk, units = J/kmol.
*/
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 the species.
*
* @param er Output vector contain the nondimensional entropies
* of the species in their reference states
* length: m_kk, units: dimensionless.
*/
virtual void getEntropy_R_ref(doublereal *er) const;
/*!
* 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.
*
* @param cprt Output vector contains the nondimensional heat capacities
* of the species in their reference states
* length: m_kk, units: dimensionless.
*/
virtual void getCp_R_ref(doublereal *cprt) const;
//! 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
*
* @param vol Output vector containing the standard state volumes.
* Length: m_kk.
*/
virtual void getStandardVolumes_ref(doublereal *vol) const;
protected:
//@}
public:
//! @name Initialization Methods - For Internal use (VPStandardState)
/*!
* The following methods are used in the process of constructing
* the phase and setting its parameters from a specification in an
* input file. They are not normally used in application programs.
* To see how they are used, see files importCTML.cpp and
* ThermoFactory.cpp.
*/
//@{
//! 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
* the "state" entry for this phase in the input file.
*/
virtual void setStateFromXML(const XML_Node& state);
//! @internal Initialize the object
/*!
* 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 after calling installSpecies()
* for each species in the phase. It's called before calling
* initThermoXML() for the phase. Therefore, it's the correct
* place for initializing vectors which have lengths equal to the
* number of species.
*
* @see importCTML.cpp
*/
virtual void initThermo();
//! Initialize a ThermoPhase object, potentially reading activity
//! coefficient information from an XML database.
/*!
* This routine initializes the lengths in the current object and
* then calls the parent routine.
* 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 just prior to returning
* from function importPhase().
*
* @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, 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
//! Calculate the value of z
/*!
* \f[
* z = \frac{P v}{ R T}
* \f]
*
* returns the value of z
*/
doublereal z() const;
//! 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.
*/
virtual doublereal sresid() const;
//! 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.
*/
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
* accurate value for the saturation pressure.
*
* @param TKelvin temperature in kelvin
*
* @return returns the estimated saturation pressure at the given temperature
*/
virtual doublereal psatEst(doublereal TKelvin) const;
//! 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.
*/
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const;
protected:
//! 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 phaseRequested 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 accectable density at a
* different phase, we return a -2.
*/
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested,
doublereal rhoguess);
//! 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 delGRT output delGRT
*/
int corr0(doublereal TKelvin, doublereal pre, doublereal &densLiq,
doublereal &densGas, doublereal &liqGRT, doublereal &gasGRT);
public:
//! 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 paramters 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 phaseState(bool checkState = false) const ;
//! 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
*/
virtual doublereal densSpinodalLiquid() const;
//! 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
*/
virtual doublereal densSpinodalGas() const;
public:
//! Calculate the saturation pressure at the current mixture content for the given temperature
/*!
* @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
*/
doublereal calculatePsat(doublereal TKelvin, doublereal &molarVolGas,
doublereal &molarVolLiquid);
protected:
//! 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.
*/
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
*
* @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
*/
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 pressurees
/*!
* Because the pressure is now a calculation, we store the result of the calculation whenever
* it is recalculated.
*
* units = Pascals
*/
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
*/
int iState_;
//! Force the system to be on a particular side of the spinodal curve
int forcedState_;
//! The last temperature at which the reference state thermodynamic properties were calculated at.
mutable doublereal m_Tlast_ref;
//! Temporary storage for log of p/rt
mutable doublereal m_logc0;
//! Temporary storage for dimensionless reference state enthalpies
mutable array_fp m_h0_RT;
//! Temporary storage for dimensionless reference state heat capacities
mutable array_fp m_cp0_R;
//! Temporary storage for dimensionless reference state gibbs energies
mutable array_fp m_g0_RT;
//! Temporary storage for dimensionless reference state entropies
mutable array_fp m_s0_R;
spinodalFunc *fdpdv_;
private:
//! MixtureFugacityTP has its own err routine
/*!
* @param msg Error message string
*/
doublereal err(std::string msg) const;
};
}
#endif

5
Cantera/src/thermo/NasaThermo.h
View file

@ -209,11 +209,12 @@ namespace Cantera {
if (m_p0 < 0.0) {
m_p0 = refPressure;
} else if (fabs(m_p0 - refPressure) > 0.1) {
std::string logmsg = " WARNING NasaThermo: New Species, " + name + ", has a different reference pressure, "
std::string logmsg = " ERROR NasaThermo: New Species, " + name + ", has a different reference pressure, "
+ fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
writelog(logmsg);
logmsg = " This may become a fatal error in the future \n";
logmsg = " This is now a fatal error\n";
writelog(logmsg);
throw CanteraError("install()", "species have different reference pressures");
}
m_p0 = refPressure;
}

1895
Cantera/src/thermo/RedlichKwongMFTP.cpp Normal file
View file

File diff suppressed because it is too large Load diff

843
Cantera/src/thermo/RedlichKwongMFTP.h Normal file
View file

@ -0,0 +1,843 @@
/**
* @file RedlichKwongMFTP.h
* Definition file for a derived class of ThermoPhase that assumes either
* an ideal gas or ideal solution approximation and handles
* variable pressure standard state methods for calculating
* thermodynamic properties (see \ref thermoprops and
* class \link Cantera::RedlichKwongMFTP RedlichKwongMFTP\endlink).
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
/*
* $Author: hkmoffa $
* $Date: 2009-11-09 16:36:49 -0700 (Mon, 09 Nov 2009) $
* $Revision: 255 $
*/
#ifndef CT_REDLICHKWONGMFTP_H
#define CT_REDLICHKWONGMFTP_H
#include "MixtureFugacityTP.h"
namespace Cantera {
class XML_Node;
class PDSS;
/*!
* @name CONSTANTS - Models for the Standard State of IdealSolnPhase's
*/
//@{
#ifdef WITH_REAL_GASSES
/**
* @ingroup thermoprops
*
* 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 for %RedlichKwongMFTP
*
*/
//! Base constructor.
RedlichKwongMFTP();
//! Construct and initialize a RedlichKwongMFTP ThermoPhase object
//! directly from an asci 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 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 empty string.
*/
RedlichKwongMFTP(std::string infile, std::string id="");
//! Construct and initialize a RedlichKwongMFTP ThermoPhase 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)
*/
RedlichKwongMFTP(XML_Node& phaseRef, std::string id = "");
//! 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
*
* @param testProb Hard -coded test problem to instantiate.
* Current valid values are 1.
*/
RedlichKwongMFTP(int testProb);
//! Copy Constructor
/*!
* Copy constructor for the object. Constructed object will be a clone of this object, but will
* also own all of its data. This is a wrapper around the assignment operator
*
* @param right Object to be copied.
*/
RedlichKwongMFTP(const RedlichKwongMFTP &right);
//! Asignment 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.
*/
RedlichKwongMFTP& operator=(const RedlichKwongMFTP &right);
//! Destructor.
virtual ~RedlichKwongMFTP();
//! 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
*/
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
* non-zero value. Constants defined for this purpose are
* listed in mix_defs.h.
*/
virtual int eosType() const;
//@}
/// Molar enthalpy. Units: J/kmol.
virtual doublereal enthalpy_mole() const;
/// Molar internal energy. Units: J/kmol.
virtual doublereal intEnergy_mole() const;
/// Molar entropy. Units: J/kmol/K.
virtual doublereal entropy_mole() const;
/// Molar Gibbs function. Units: J/kmol.
virtual doublereal gibbs_mole() const;
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual doublereal cp_mole() const;
/// Molar heat capacity at constant volume. Units: J/kmol/K.
virtual doublereal cv_mole() const;
/**
* @}
* @name Mechanical Properties
* @{
*/
//! 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]
*
*/
virtual doublereal pressure() const;
//! Returns the isothermal compressibility. Units: 1/Pa.
/*!
* The isothermal compressibility is defined as
* \f[
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
* \f]
*/
virtual doublereal isothermalCompressibility() const;
protected:
/**
* 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
* 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();
//! Set the temperature (K)
/*!
* Overwritten setTemperature(double) from State.h. This
* function sets the temperature, and makes sure that
* the value propagates to underlying objects
*
* @todo Make State::setTemperature a virtual function
*
* @param temp Temperature in kelvin
*/
virtual void setTemperature(const doublereal temp);
//! Set the mass fractions to the specified values, and then
//! normalize them so that they sum to 1.0.
/*!
* @param y Array of unnormalized mass fraction values (input).
* Must have a length greater than or equal to the number of species.
*/
virtual void setMassFractions(const doublereal* const y);
//!Set the mass 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 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.
/*!
* @param x Array of unnormalized mole fraction values (input).
* Must have a length greater than or equal to the number of species.
*/
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 ofequations.
*
* @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
* units. For volumetric phases c[k] is the
* concentration of the kth species in kmol/m3.
* For surface phases, c[k] is the concentration
* in kmol/m2. The length of the vector is the number
* of species in the phase.
*/
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 /
* C^0_k, \f$ where \f$ C^0_k \f$ is a standard concentration
* defined below and \f$ a_k \f$ are activities used in the
* thermodynamic functions. These activity (or generalized)
* concentrations are used
* by kinetics manager classes to compute the forward and
* reverse rates of elementary reactions. Note that they may
* or may not have units of concentration --- they might be
* partial pressures, mole fractions, or surface coverages,
* for example.
*
* @param c Output array of generalized concentrations. The
* units depend upon the implementation of the
* reaction rate expressions within the phase.
*/
virtual void getActivityConcentrations(doublereal* c) const;
//! Returns the standard concentration \f$ C^0_k \f$, which is used to normalize
//! the generalized concentration.
/*!
* This is defined as the concentration by which the generalized
* concentration is normalized to produce the activity.
* In many cases, this quantity will be the same for all species in a phase.
* Since the activity for an ideal gas mixture is
* simply the mole fraction, for an ideal gas \f$ C^0_k = P/\hat R T \f$.
*
* @param k Optional parameter indicating the species. The default
* is to assume this refers to species 0.
* @return
* Returns the standard Concentration in units of m3 kmol-1.
*/
virtual doublereal standardConcentration(int k=0) const;
//! Returns the natural logarithm of the standard
//! concentration of the kth species
/*!
* @param k index of the species. (defaults to zero)
*/
virtual doublereal logStandardConc(int k=0) const;
//! Returns the units of the standard and generalized concentrations.
/*!
* Note they have the same units, as their
* ratio 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.
*
* The base %ThermoPhase class assigns the default quantities
* of (kmol/m3) for all species.
* 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 k species index. Defaults to 0.
* @param sizeUA output int containing the size of the vector.
* Currently, this is equal to 6.
*/
virtual void getUnitsStandardConc(double *uA, int k = 0, int sizeUA = 6) const;
//! Get the array of non-dimensional activity coefficients at
//! the current solution temperature, pressure, and solution concentration.
/*!
* For all objects with the Mixture Fugacity approximation, we define the
* standard state as an ideal gas at the current temperature and pressure
* of the solution. The activities are based on this standard state.
*
* @param ac Output vector of activity coefficients. Length: m_kk.
*/
virtual void getActivityCoefficients(doublereal* ac) const;
/// @name Partial Molar Properties of the Solution (RedlichKwongMFTP)
//@{
//! 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. No need for child
* classes to handle.
*
* @param mu Output vector of non-dimensional species chemical potentials
* Length: m_kk.
*/
void getChemPotentials_RT(doublereal* mu) const;
//! 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.
*
* @param mu Output vector of species chemical
* potentials. Length: m_kk. Units: J/kmol
*/
virtual void getChemPotentials(doublereal* mu) const;
//! Get the species partial molar enthalpies. Units: J/kmol.
/*!
* @param hbar Output vector of species partial molar enthalpies.
* Length: m_kk. units are J/kmol.
*/
virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
//! Get the species partial molar entropies. Units: J/kmol/K.
/*!
* @param sbar Output vector of species partial molar entropies.
* Length = m_kk. units are J/kmol/K.
*/
virtual void getPartialMolarEntropies(doublereal* sbar) const;
//! Get the species partial molar enthalpies. Units: J/kmol.
/*!
* @param ubar Output vector of speciar partial molar internal energies.
* Length = m_kk. units are J/kmol.
*/
virtual void getPartialMolarIntEnergies(doublereal* ubar) const;
//! Get the partial molar heat capacities 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.
*/
virtual void getPartialMolarCp(doublereal* cpbar) const;
//! Get the species partial molar volumes. Units: m^3/kmol.
/*!
* @param vbar Output vector of speciar partial molar volumes.
* Length = m_kk. units are m^3/kmol.
*/
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).
virtual doublereal critTemperature() const;
/// Critical pressure (Pa).
virtual doublereal critPressure() const;
/// Critical density (kg/m3).
virtual doublereal critDensity() const;
//@}
public:
//! @name Initialization Methods - For Internal use (VPStandardState)
/*!
* The following methods are used in the process of constructing
* the phase and setting its parameters from a specification in an
* input file. They are not normally used in application programs.
* To see how they are used, see files importCTML.cpp and
* ThermoFactory.cpp.
*/
//@{
//! 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
* any parameters that are specific to that particular phase
* model.
*
* @param thermoNode An XML_Node object corresponding to
* the "thermo" entry for this phase in the input file.
*/
virtual void setParametersFromXML(const XML_Node& thermoNode);
//! @internal Initialize the object
/*!
* 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 just prior to returning
* from function importPhase().
*
* @see importCTML.cpp
*/
virtual void initThermo();
//!This method is used by the ChemEquil equilibrium solver.
/*!
* It sets the state such that the chemical potentials satisfy
* \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
* \left(\frac{\lambda_m} {\hat R T}\right) \f] where
* \f$ \lambda_m \f$ is the element potential of element m. The
* temperature is unchanged. Any phase (ideal or not) that
* implements this method can be equilibrated by ChemEquil.
*
* @param lambda_RT Input vector of dimensionless element potentials
* The length is equal to nElements().
*/
void setToEquilState(const doublereal* lambda_RT);
//! Initialize a ThermoPhase object, potentially reading activity
//! coefficient information from an XML database.
/*!
*
* This routine initializes the lengths in the current object and
* then calls the parent routine.
* 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 just prior to returning
* from function importPhase().
*
* @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, std::string id);
private:
//! Read the pure species RedlichKwong input parameters
/*!
* @param pureFluidParam XML_Node for the pure fluid parameters
*/
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
*/
void initLengths();
//==============================================================================
// Special functions inherited from MixtureFugacityTP
protected:
//! 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.
*/
virtual doublereal sresid() const;
// 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 enthalpy in units of J kmol-1.
*/
virtual doublereal hresid() const;
//! 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.
*/
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const;
protected:
//! 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 phase 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 accectable density at a
* different phase, we return a -2.
*/
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.
/*!
* @return returns the density with units of kg m-3
*/
virtual doublereal densSpinodalLiquid() const;
//! 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
*/
virtual doublereal densSpinodalGas() const;
//! 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.
*/
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
*
* @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
*/
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const;
//! Calculate dpdV and dpdT at the current conditions
/*!
* These are storred 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.
* This function updates the internal numbers based on the state of the object.
*/
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 storred mole fractions in the object.
*
* @param temp Temperature (TKelvin)
*
* @param aCalc (output) Returns the a value
* @param bCalc (output) Returns the b value.
*/
void calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const;
//=========================================================================================
// 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;
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.
* - 0 = No, there are nonstaandard cross terms in the a coefficient matrices.
*/
int m_standardMixingRules;
//! Form of the temperature parameterization
/*!
* - 0 = There is no temperature parameterization of a or b
* - 1 = The a_ij parameter is a linear function of the temperature
*/
int m_formTempParam;
//! Value of b in the equation of state
/*!
* m_b is a function of the temperature and the mole fraction.
*/
doublereal m_b_current;
//! Value of a in the equation of state
/*!
* a_b is a function of the temperature and the mole fraction.
*/
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;
int NSolns_;
doublereal Vroot_[3];
//! Temporary storage - length = m_kk.
mutable vector_fp m_pp;
//! Temporary storage - length = m_kk.
mutable vector_fp m_tmpV;
// mutable vector_fp m_tmpV2;
// Partial molar volumes of the species
mutable vector_fp m_partialMolarVolumes;
//! The derivative of the pressure wrt the volume
/*!
* Calcualted at the current conditions
* temperature and mole number kept constant
*/
mutable doublereal dpdV_;
//! The derivative of the pressure wrt the temperature
/*!
* Calcualted at the current conditions
* Total volume and mole number kept constant
*/
mutable doublereal dpdT_;
//! Vector of derivatives of pressure wrt mole number
/*!
* Calcualted at the current conditions
* Total volume, temperature and other mole number kept constant
*/
mutable vector_fp dpdni_;
public:
//! Omega constant for a -> value of a in terms of critical properties
/*!
* this was calculated from a small nonlinear solve
*/
static const doublereal omega_a = 4.27480233540E-01;
//! Omega constant for b
static const doublereal omega_b = 8.66403499650E-02;
//! Omega constant for the critical molar volume
static const doublereal omega_vc = 3.33333333333333E-01;
};
#endif
}
#endif

5
Cantera/src/thermo/ShomateThermo.h
View file

@ -201,12 +201,13 @@ namespace Cantera {
if (m_p0 < 0.0) {
m_p0 = refPressure;
} else if (fabs(m_p0 - refPressure) > 0.1) {
std::string logmsg = " WARNING ShomateThermo: New Species, " + name
std::string logmsg = " ERROR ShomateThermo: New Species, " + name
+ ", has a different reference pressure, "
+ fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
writelog(logmsg);
logmsg = " This may become a fatal error in the future \n";
logmsg = " This is now a fatal error\n";
writelog(logmsg);
throw CanteraError("install()", "Species have different reference pressures");
}
m_p0 = refPressure;

3
Cantera/src/thermo/SimpleThermo.h
View file

@ -182,8 +182,9 @@ namespace Cantera {
", has a different reference pressure, "
+ fp2str(refPressure) + ", than existing reference pressure, " + fp2str(m_p0) + "\n";
writelog(logmsg);
logmsg = " This may become a fatal error in the future \n";
logmsg = " This is now a fatal error\n";
writelog(logmsg);
throw CanteraError("install()", "Species have different reference pressures");
}
m_p0 = refPressure;
}

11
Cantera/src/thermo/ThermoFactory.cpp
View file

@ -36,6 +36,10 @@
#include "PureFluidPhase.h"
#endif
#ifdef WITH_REAL_GASSES
#include "RedlichKwongMFTP.h"
#endif
#include "ConstDensityThermo.h"
#include "SurfPhase.h"
#include "EdgePhase.h"
@ -212,6 +216,13 @@ namespace Cantera {
th = new PureFluidPhase;
break;
#endif
#ifdef WITH_REAL_GASSES
case cRedlichKwongMFTP:
th = new RedlichKwongMFTP;
break;
#endif
#ifdef WITH_ELECTROLYTES
case cHMW:
th = new HMWSoln;

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

@ -1012,7 +1012,7 @@ namespace Cantera {
}
return m_speciesData;
}
//====================================================================================================================
/*
* Set the thermodynamic state.
*/
@ -1038,7 +1038,7 @@ namespace Cantera {
setDensity(rho);
}
}
//====================================================================================================================
/*
* Called by function 'equilibrate' in ChemEquil.h to transfer
* the element potentials to this object after every successful

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

@ -1544,7 +1544,7 @@ namespace Cantera {
* @param x Vector of mole fractions.
* Length is equal to m_kk.
*/
void setState_TPX(doublereal t, doublereal p, const doublereal* x);
virtual void setState_TPX(doublereal t, doublereal p, const doublereal* x);
//! Set the temperature (K), pressure (Pa), and mole fractions.
/*!

4
Cantera/src/thermo/mix_defs.h
View file

@ -73,6 +73,10 @@ namespace Cantera {
const int cIdealSolnGasVPSS = 500;
const int cIdealSolnGasVPSS_iscv = 501;
//! Fugacity Models
const int cMixtureFugacityTP = 700;
const int cRedlichKwongMFTP = 701;
const int cMargulesVPSSTP = 301;
const int cPhaseCombo_Interaction = 305;

3
config.h.in
View file

@ -183,6 +183,9 @@ typedef int ftnlen; // Fortran hidden string length type
// models for electrolyte solutions.
#undef WITH_ELECTROLYTES
// Enable Real Gasses
#undef WITH_REAL_GASSES
#undef WITH_PRIME
// Enable the VCS NonIdeal equilibrium solver. This is

16
configure vendored
View file

@ -2585,6 +2585,16 @@ _ACEOF
fi
if test "$WITH_REAL_GASSES" = "y"; then
cat >>confdefs.h <<\_ACEOF
#define WITH_REAL_GASSES 1
_ACEOF
hdrs=$hdrs' RedlichKwongMFTP.h'
objs=$objs' RedlichKwongMFTP.o'
fi
if test "$WITH_LATTICE_SOLID" = "y"; then
cat >>confdefs.h <<\_ACEOF
#define WITH_LATTICE_SOLID 1
@ -10058,7 +10068,7 @@ fi
# Provide some information about the compiler.
echo "$as_me:10061:" \
echo "$as_me:10071:" \
"checking for Fortran 77 compiler version" >&5
ac_compiler=`set X $ac_compile; echo $2`
{ (eval echo "$as_me:$LINENO: \"$ac_compiler --version </dev/null >&5\"") >&5
@ -10265,7 +10275,7 @@ _ACEOF
# flags.
ac_save_FFLAGS=$FFLAGS
FFLAGS="$FFLAGS $ac_verb"
(eval echo $as_me:10268: \"$ac_link\") >&5
(eval echo $as_me:10278: \"$ac_link\") >&5
ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'`
echo "$ac_f77_v_output" >&5
FFLAGS=$ac_save_FFLAGS
@ -10343,7 +10353,7 @@ _ACEOF
# flags.
ac_save_FFLAGS=$FFLAGS
FFLAGS="$FFLAGS $ac_cv_prog_f77_v"
(eval echo $as_me:10346: \"$ac_link\") >&5
(eval echo $as_me:10356: \"$ac_link\") >&5
ac_f77_v_output=`eval $ac_link 5>&1 2>&1 | grep -v 'Driving:'`
echo "$ac_f77_v_output" >&5
FFLAGS=$ac_save_FFLAGS

7
configure.in
View file

@ -486,6 +486,13 @@ if test "$WITH_PURE_FLUIDS" = "y"; then
fi
AC_SUBST(COMPILE_PURE_FLUIDS)
if test "$WITH_REAL_GASSES" = "y"; then
AC_DEFINE(WITH_REAL_GASSES)
hdrs=$hdrs' RedlichKwongMFTP.h'
objs=$objs' RedlichKwongMFTP.o'
fi
if test "$WITH_LATTICE_SOLID" = "y"; then
AC_DEFINE(WITH_LATTICE_SOLID)
hdrs=$hdrs' LatticeSolidPhase.h'

3
docs/install_examples/linux.64_sierra_gcc444_python264_numpy
View file

@ -34,6 +34,9 @@ export WITH_VCSNONIDEAL
WITH_H298MODIFY_CAPABILITY='y'
export WITH_H298MODIFY_CAPABILITY
WITH_REAL_GASSES='y'
export WITH_REAL_GASSES
BUILD_MATLAB_TOOLBOX="y"
export BUILD_MATLAB_TOOLBOX

4
preconfig
View file

@ -228,6 +228,10 @@ WITH_IDEAL_SOLUTIONS=${WITH_IDEAL_SOLUTIONS:="y"}
# models for electrolyte solutions
WITH_ELECTROLYTES=${WITH_ELECTROLYTES:="y"}
# Enable Real Gas Equations of State
# This supports the multicomponent Redlich-Kwong equation of state.
WITH_REAL_GASSES=${WITH_REAL_GASSES:="n"}
# Enable generating phase models from PrIMe models. For more
# information about PrIME, see http://www.primekinetics.org
# WARNING: Support for PrIMe is experimental!

4
test_problems/min_python/negATest/negATest_blessed.out
View file

@ -1,7 +1,7 @@
Number of species = 12
Number of reactions = 12
rop [ 0:O ] = 0.447058
rop [ 1:O2 ] = -0.0021443
rop [ 1:O2 ] = -0.00214428
rop [ 2:N ] = 0
rop [ 3:NO ] = -279.361
rop [ 4:NO2 ] = 0.00214319
@ -17,7 +17,7 @@ fwd_rop[ 0] = 479.304 rev_rop[ 0] = 97.94
fwd_rop[ 1] = -128.201 rev_rop[ 1] = -26.1964
fwd_rop[ 2] = 0 rev_rop[ 2] = 0
fwd_rop[ 3] = -0 rev_rop[ 3] = -0
fwd_rop[ 4] = 0 rev_rop[ 4] = 1.10334e-06
fwd_rop[ 4] = 0 rev_rop[ 4] = 1.08891e-06
fwd_rop[ 5] = 0 rev_rop[ 5] = 0
fwd_rop[ 6] = 0 rev_rop[ 6] = 0
fwd_rop[ 7] = 0 rev_rop[ 7] = 0