From 7971238d37697ca0f106ee3cad9b2befdfa92ae9 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 13 Jun 2006 16:02:41 +0000 Subject: [PATCH] Added an IdealSolidSolnPhase capability. Changed existing thermo files to use initThermoXML(). --- Cantera/src/thermo/IdealSolidSolnPhase.cpp | 1345 ++++++++++++++++++++ Cantera/src/thermo/IdealSolidSolnPhase.h | 1038 +++++++++++++++ Cantera/src/thermo/Makefile.in | 16 +- Cantera/src/thermo/MolalityVPSSTP.cpp | 135 +- Cantera/src/thermo/MolalityVPSSTP.h | 36 +- Cantera/src/thermo/VPStandardStateTP.cpp | 44 +- Cantera/src/thermo/VPStandardStateTP.h | 21 +- 7 files changed, 2587 insertions(+), 48 deletions(-) create mode 100644 Cantera/src/thermo/IdealSolidSolnPhase.cpp create mode 100644 Cantera/src/thermo/IdealSolidSolnPhase.h diff --git a/Cantera/src/thermo/IdealSolidSolnPhase.cpp b/Cantera/src/thermo/IdealSolidSolnPhase.cpp new file mode 100644 index 000000000..7cdf830e3 --- /dev/null +++ b/Cantera/src/thermo/IdealSolidSolnPhase.cpp @@ -0,0 +1,1345 @@ +/** + * @file IdealSolidSolnPhase.cpp + * Code for the class IdealSolidSolnPhase, + * which implements an ideal solid solution model + * with incompressible thermodynamics. + */ +/* + * $Id$ + */ +/* + * Copywrite 2006 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000, with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + */ + +#include "IdealSolidSolnPhase.h" +#include +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), + m_mm(0), + m_tmin(0.0), + m_tmax(1000000.), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(0.0) + { + if (formGC < 0 || formGC > 2) { + throw CanteraError(" IdealSolidSolnPhase Constructor", + " Illegal value of formGC"); + } + } + + IdealSolidSolnPhase::IdealSolidSolnPhase(string inputFile, string id, + int formGC) : + ThermoPhase(), + m_formGC(formGC), + m_mm(0), + m_tmin(0.0), + m_tmax(1000000.), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(0.0) + { + if (formGC < 0 || formGC > 2) { + throw CanteraError(" IdealSolidSolnPhase Constructor", + " Illegal value of formGC"); + } + constructPhaseFile(inputFile, id); + } + + IdealSolidSolnPhase::IdealSolidSolnPhase(XML_Node& root, string id, + int formGC) : + ThermoPhase(), + m_formGC(formGC), + m_mm(0), + m_tmin(0.0), + m_tmax(1000000.), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(0.0) + { + if (formGC < 0 || formGC > 2) { + throw CanteraError(" IdealSolidSolnPhase Constructor", + " Illegal value of formGC"); + } + constructPhaseXML(root, id); + } + + IdealSolidSolnPhase::IdealSolidSolnPhase(const IdealSolidSolnPhase &b) + { + *this = b; + } + + + IdealSolidSolnPhase& IdealSolidSolnPhase:: + operator=(const IdealSolidSolnPhase &b) { + if (this != &b) { + //ThermoPhase::operator=(b); + // m_spthermo = dupMyselfAsSpeciesThermo(b.m_spthermo); + m_formGC = b.m_formGC; + m_mm = b.m_mm; + m_tmin = b.m_tmin; + m_tmax = b.m_tmax; + m_Pref = b.m_Pref; + m_Pcurrent = b.m_Pcurrent; + m_speciesMolarVolume = b.m_speciesMolarVolume; + m_tlast = b.m_tlast; + m_h0_RT = b.m_h0_RT; + m_cp0_R = b.m_cp0_R; + m_g0_RT = b.m_g0_RT; + m_s0_R = b.m_s0_R; + m_expg0_RT = b.m_expg0_RT; + m_pe = b.m_pe; + m_pp = b.m_pp; + } + 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() { + IdealSolidSolnPhase *ii = new IdealSolidSolnPhase(*this); + return (ThermoPhase*) ii; + } + + /** + * Equation of state flag. Returns the value cIdealGas, defined + * in mix_defs.h. + */ + int IdealSolidSolnPhase::eosType() const { + integer res; + switch (m_formGC) { + case 0: + res = cIdealSolidSolnPhase0; + break; + case 1: + res = cIdealSolidSolnPhase1; + break; + case 2: + res = cIdealSolidSolnPhase2; + break; + default: + throw CanteraError("eosType", "Unknown type"); + break; + } + return res; + } + + /******************************************************************** + * 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 { + const double *eptr = &(enthalpy_RT_ref()[0]); + doublereal htp = (GasConstant * temperature() * mean_X(eptr)); + 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()); + doublereal htp = (GasConstant * temperature() * + mean_X(eptr)); + 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()); + doublereal g = mean_X(dptr); + 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()); + return GasConstant * mean_X(dptr); + } + + /******************************************************************** + * 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() { + /* + * Calculate the molarVolume of the solution (m**3 kmol-1) + */ + const doublereal * const dtmp = moleFractdivMMW(); + double invDens = dot(m_speciesMolarVolume.begin(), + m_speciesMolarVolume.end(), dtmp); + /* + * Set the density in the parent State object directly, + * by calling the State::setDensity() function. + */ + double dens = 1.0/invDens; + State::setDensity(dens); + } + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent 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(doublereal rho) { + /* + * Unless the input density is exactly equal to the density + * calculated and storred in the State object, we throw an + * exception. This is because the density is NOT an + * independent variable. + */ + double dens = density(); + if (rho != dens) { + throw CanteraError("IdealSolidSolnPhase::setDensity", + "Density is not an independent variable"); + } + } + + /** + * 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 indendent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is a virtual function that overrides the State.h + * class + */ + void IdealSolidSolnPhase::setMolarDensity(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(doublereal *x) { + State::setMoleFractions(x); + calcDensity(); + } + + /** + * setMoleFractions_NoNorm() (virtual from State) + * + * Sets the mole fractions and adjusts the internal density. + */ + void IdealSolidSolnPhase::setMoleFractions_NoNorm(doublereal *x) { + State::setMoleFractions(x); + calcDensity(); + } + + /** + * setMassFractions() (virtual from State) + * + * Sets the mass fractions and adjusts the internal density. + */ + void IdealSolidSolnPhase::setMassFractions(doublereal *y) { + State::setMassFractions(y); + calcDensity(); + } + + /** + * setMassFractions_NoNorm() (virtual from State) + * + * Sets the mass fractions and adjusts the internal density. + */ + void IdealSolidSolnPhase::setMassFractions_NoNorm(doublereal *y) { + State::setMassFractions_NoNorm(y); + calcDensity(); + } + + /** + * setConcentrations (virtual from State) + * + * Sets the concentrations and adjusts the internal density + */ + void IdealSolidSolnPhase::setConcentrations(doublereal *c) { + State::setConcentrations(c); + calcDensity(); + } + + /******************************************************************** + * 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 dividied by the standard concentration is also + * equal to the activity of species. + * + * For this implentation 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-3. 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: + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * 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 awkwards 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 { + const doublereal * const dtmp = moleFractdivMMW(); + const double mmw = meanMolecularWeight(); + switch (m_formGC) { + case 0: + for (int k = 0; k < m_kk; k++) { + c[k] = dtmp[k] * mmw; + } + break; + case 1: + for (int k = 0; k < m_kk; k++) { + c[k] = dtmp[k] * mmw / m_speciesMolarVolume[k]; + } + break; + case 2: + double atmp = mmw / m_speciesMolarVolume[m_kk-1]; + for (int k = 0; k < m_kk; k++) { + c[k] = dtmp[k] * atmp; + } + break; + } + } + + /********************************************************************* + * + * 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 m3 kmol-1. + * + * + * @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(int k) const { + switch (m_formGC) { + case 0: + return 1.0; + break; + case 1: + return 1.0 / m_speciesMolarVolume[k]; + break; + case 2: + return 1.0/m_speciesMolarVolume[m_kk-1]; + break; + } + return 0.0; + } + doublereal IdealSolidSolnPhase:: + referenceConcentration(int k) const { + switch (m_formGC) { + case 0: + return 1.0; + break; + case 1: + return 1.0 / m_speciesMolarVolume[k]; + break; + case 2: + return 1.0 / m_speciesMolarVolume[m_kk-1]; + break; + } + 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(int k) const { + _updateThermo(); + double res; + switch (m_formGC) { + case 0: + res = 0.0; + break; + case 1: + res = log(1.0/m_speciesMolarVolume[k]); + break; + case 2: + res = log(1.0/m_speciesMolarVolume[m_kk-1]); + break; + default: + throw CanteraError("eosType", "Unknown type"); + break; + } + 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 concentrtion is + * unitless. + */ + void IdealSolidSolnPhase:: + getUnitsStandardConc(double *uA, int, int sizeUA) { + int eos = eosType(); + if (eos == cIdealSolidSolnPhase0) { + for (int i = 0; i < sizeUA; i++) { + uA[i] = 0.0; + } + } else { + for (int i = 0; i < sizeUA; i++) { + if (i == 0) uA[0] = 1.0; + if (i == 1) uA[1] = -nDim(); + if (i == 2) uA[2] = 0.0; + if (i == 3) uA[3] = 0.0; + if (i == 4) uA[4] = 0.0; + if (i == 5) uA[5] = 0.0; + } + } + } + + /******************************************************************** + * getActivityCoefficients(): + * + */ + void IdealSolidSolnPhase:: + getActivityCoefficients(doublereal *ac) const { + for (int k = 0; k < m_kk; k++) { + 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 { + doublereal delta_p = m_Pcurrent - m_Pref; + doublereal xx; + doublereal RT = temperature() * GasConstant; + const array_fp& g_RT = gibbs_RT_ref(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(SmallNumber, moleFraction(k)); + mu[k] = RT * (g_RT[k] + log(xx)) + + 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 { + doublereal RT = temperature() * GasConstant; + doublereal delta_pdRT = (m_Pcurrent - m_Pref) / RT; + doublereal xx; + const array_fp& g_RT = gibbs_RT_ref(); + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(SmallNumber, moleFraction(k)); + mu[k] = (g_RT[k] + log(xx)) + + delta_pdRT * m_speciesMolarVolume[k]; + } + } + + /******************************************************************** + * 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 array_fp& _h = enthalpy_RT_ref(); + doublereal rt = GasConstant * temperature(); + 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 { + const array_fp& _s = entropy_R_ref(); + doublereal r = GasConstant; + doublereal xx; + for (int k = 0; k < m_kk; k++) { + xx = fmaxx(SmallNumber, moleFraction(k)); + sbar[k] = r * (_s[k] - log(xx)); + } + } + + /******************************************************************** + * + * getPartialMolarCp() + * + * For this phase, the partial molar heat capacities are equal + * to the standard state heat capacities. + + */ + void IdealSolidSolnPhase:: + getPartialMolarCp(doublereal* cpbar) const { + getCp_R(cpbar); + for (int k = 0; k < m_kk; k++) { + cpbar[k] *= GasConstant; + } + } + + /****************************************************************** + * + * 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 { + getStandardVolumes(vbar); + } + + /***************************************************************** + * Properties of the Standard State of the Species + * in the Solution + *****************************************************************/ + + /****************************************************************** + * + * getPureGibbs() + * + * Get the Gibbs functions for the pure species + * at the current T and P 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 k<\I>. + * \f$ u^{ref}_k(T)\f$ is the chemical potential of pure + * species k<\I> at the reference pressure, \f$P_{ref}\f$. + */ + void IdealSolidSolnPhase:: + getPureGibbs(doublereal* gpure) const { + const array_fp& gibbsrt = gibbs_RT_ref(); + doublereal RT = _RT(); + const doublereal * const gk = DATA_PTR(gibbsrt); + doublereal delta_p = (m_Pcurrent - m_Pref); + for (int k = 0; k < m_kk; k++) { + gpure[k] = RT * gk[k] + delta_p * m_speciesMolarVolume[k]; + } + } + + /** + * 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 k<\I>. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species 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 { + const array_fp& gibbsrt = gibbs_RT_ref(); + doublereal RT = _RT(); + const doublereal * const gk = DATA_PTR(gibbsrt); + doublereal delta_prt = (m_Pcurrent - m_Pref)/ RT; + for (int k = 0; k < m_kk; k++) { + grt[k] = gk[k] + delta_prt * m_speciesMolarVolume[k]; + } + } + + /******************************************************************** + * + * getEnthalpy_RT() + * + * Get the array of nondimensional Enthalpy functions for the ss + * species at the current T and P 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 k<\I>. + * \f$ h^{ref}_k(T)\f$ is the enthalpy of the pure + * species k<\I> at the reference pressure, \f$P_{ref}\f$. + */ + void IdealSolidSolnPhase:: + getEnthalpy_RT(doublereal* hrt) const { + const array_fp& _h = enthalpy_RT_ref(); + doublereal delta_prt = ((m_Pcurrent - m_Pref) / + (GasConstant * temperature())); + for (int k = 0; k < m_kk; k++) { + hrt[k] = _h[k] + delta_prt * m_speciesMolarVolume[k]; + } + } + + /** + * 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 array_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 array_fp& _h = enthalpy_RT_ref(); + doublereal prefrt = m_Pref / (GasConstant * temperature()); + for (int k = 0; k < m_kk; k++) { + urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k]; + } + } + + /** + * 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 k<\I>. + * \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity + * of species 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 array_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 + * T and P 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(); + for (int k = 0; k != m_kk; k++) { + hrt[k] = m_h0_RT[k]; + } + } + + /** + * 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(); + for (int k = 0; k != m_kk; k++) { + grt[k] = m_g0_RT[k]; + } + } + + /** + * 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(); + double tmp = GasConstant * temperature(); + for (int k = 0; k != m_kk; k++) { + g[k] = tmp * m_g0_RT[k]; + } + } + + /** + * 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 array_fp& _h = enthalpy_RT_ref(); + doublereal prefrt = m_Pref / (GasConstant * temperature()); + for (int k = 0; k < m_kk; k++) { + urt[k] = _h[k] - prefrt * m_speciesMolarVolume[k]; + } + } + + /** + * 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(); + for (int k = 0; k != m_kk; k++) { + er[k] = m_s0_R[k]; + } + } + + /** + * 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(); + for (int k = 0; k != m_kk; k++) { + cpr[k] = m_cp0_R[k]; + } + } + + /** + * 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 array_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 array_fp& IdealSolidSolnPhase::expGibbs_RT_ref() const { + _updateThermo(); + int k; + for (k = 0; k != m_kk; k++) m_expg0_RT[k] = exp(m_g0_RT[k]); + return m_expg0_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 array_fp& IdealSolidSolnPhase::entropy_R_ref() const { + _updateThermo(); + return m_s0_R; + } + + /********************************************************************* + * Utility Functions + *********************************************************************/ + /** + * initThermo() function initializes the object for use. + * + * Before its invokation, the class isn't ready for calculation. + */ + void IdealSolidSolnPhase::initThermo() { + } + + /** + * Import and initialize an IdealSolidSolnPhase phase + * specification in an XML tree into the current object. + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * + * This routine calls importPhase() to do most of its work. + * Then, importPhase() calls initThermoXML() to finish + * off the work. + * + * @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:: + constructPhaseXML(XML_Node& phaseNode, string id) { + string subname = "IdealSolidSolnPhase::constructPhaseXML"; + if (id.size() > 0) { + string idp = phaseNode.id(); + if (idp != id) { + throw CanteraError(subname.c_str(), + "phasenode and Id are incompatible"); + } + } + + /* + * Check on the thermo field. Must have: + * + */ + if (phaseNode.hasChild("thermo")) { + XML_Node& thNode = phaseNode.child("thermo"); + string mStringa = thNode.attrib("model"); + string mString = lowercase(mStringa); + if (mString != "idealsolidsolution") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa); + } + } else { + throw CanteraError(subname.c_str(), + "Unspecified thermo model"); + } + + /* + * Form of the standard concentrations. Must have one of: + * + * + * + * + */ + if (phaseNode.hasChild("standardConc")) { + XML_Node& scNode = phaseNode.child("standardConc"); + string formStringa = scNode.attrib("model"); + string formString = lowercase(formStringa); + if (formString == "unity") { + m_formGC = 0; + } else if (formString == "molar_volume") { + m_formGC = 1; + } else if (formString == "solvent_volume") { + m_formGC = 2; + } else { + throw CanteraError(subname.c_str(), + "Unknown standardConc model: " + formStringa); + } + } else { + throw CanteraError(subname.c_str(), + "Unspecified standardConc model"); + } + + bool m_ok = importPhase(phaseNode, this); + if (!m_ok) { + throw CanteraError(subname.c_str(),"importPhase failed "); + } + } + + /** + * Initialization of an IdealSolidSolnPhase phase using an + * xml file + * + * This routine is a precursor to constructPhaseFile(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + void IdealSolidSolnPhase:: + constructPhaseFile(string inputFile, string id) { + if (inputFile.size() == 0) { + throw CanteraError("IdealSolidSolnPhase::constructPhaseFile", + "input file is null"); + } + string path = findInputFile(inputFile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("IdealSolidSolnPhase::constructPhaseFile","could not open " + +path+" for reading."); + } + /* + * The phase object automatically constructs an XML object. + * Use this object to store information. + */ + XML_Node &phaseNode_XML = xml(); + XML_Node *fxml = new XML_Node(); + fxml->build(fin); + XML_Node *fxml_phase = findXMLPhase(fxml, id); + if (!fxml_phase) { + throw CanteraError("IdealSolidSolnPhase::constructPhaseFile", + "ERROR: Can not find phase named " + + id + " in file named " + inputFile); + } + fxml_phase->copy(&phaseNode_XML); + + constructPhaseXML(*fxml_phase, id); + delete fxml; + } + + /** + * @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, string id) { + /* + * Initialize all of the lengths now that we know how many species + * there are in the phase. + */ + initLengths(); + /* + * Now go get the molar volumes + */ + XML_Node& speciesList = phaseNode.child("speciesArray"); + XML_Node* speciesDB = get_XML_NameID("speciesData", speciesList["datasrc"], + &phaseNode.root()); + const vector&sss = speciesNames(); + + for (int k = 0; k < m_kk; k++) { + XML_Node* s = speciesDB->findByAttr("name", sss[k]); + XML_Node *ss = s->findByName("standardState"); + m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "-"); + } + + /* + * Call the base initThermo, which handles setting the initial + * state. + */ + ThermoPhase::initThermoXML(phaseNode, id); + } + + /** + * This internal function adjusts the lengths of arrays + */ + void IdealSolidSolnPhase:: + initLengths() { + m_kk = nSpecies(); + m_mm = nElements(); + + /* + * Obtain the limits of the temperature from the species + * thermo handler's limits. + */ + doublereal tmin = m_spthermo->minTemp(); + doublereal tmax = m_spthermo->maxTemp(); + if (tmin > 0.0) m_tmin = tmin; + if (tmax > 0.0) m_tmax = tmax; + + /* + * Obtain the reference pressure by calling the ThermoPhase + * function refPressure, which in turm calls the + * species thermo reference pressure function of the + * same name. + */ + m_Pref = refPressure(); + + int leng = m_kk; + m_h0_RT.resize(leng); + m_g0_RT.resize(leng); + m_expg0_RT.resize(leng); + m_cp0_R.resize(leng); + m_s0_R.resize(leng); + m_pe.resize(leng, 0.0); + m_pp.resize(leng); + m_speciesMolarVolume.resize(leng); + } + + /** + * 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) + { + const array_fp& grt = gibbs_RT_ref(); + + // set the pressure and composition to be consistent with + // the temperature, + doublereal pres = 0.0; + for (int k = 0; k < m_kk; k++) { + m_pp[k] = -grt[k]; + for (int m = 0; m < m_mm; m++) { + m_pp[k] += nAtoms(k,m)*lambda_RT[m]; + } + m_pp[k] = m_Pref * exp(m_pp[k]); + pres += m_pp[k]; + } + 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 { + doublereal tnow = temperature(); + if (m_tlast != tnow) { + /* + * Update the thermodynamic functions of the reference state. + */ + m_spthermo->update(tnow, DATA_PTR(m_cp0_R), DATA_PTR(m_h0_RT), + DATA_PTR(m_s0_R)); + m_tlast = tnow; + doublereal rrt = 1.0 / (GasConstant * tnow); + int k; + doublereal deltaE; + for (k = 0; k < m_kk; k++) { + deltaE = rrt * m_pe[k]; + m_h0_RT[k] += deltaE; + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } + m_tlast = tnow; + } + } +} diff --git a/Cantera/src/thermo/IdealSolidSolnPhase.h b/Cantera/src/thermo/IdealSolidSolnPhase.h new file mode 100644 index 000000000..74bf25b7a --- /dev/null +++ b/Cantera/src/thermo/IdealSolidSolnPhase.h @@ -0,0 +1,1038 @@ +/** + * @file IdealSolidSolnPhase.h + * Header file for the class IdealSolidSolnPhase + * This class implements an ideal solid solution model + * with incompressible thermodynamics. + * + * This class inherits from the Cantera class ThermoPhase + * and implements an ideal solid solution model with incompressible + * thermodynamics. + * + * The concept of a monomer unit is mapped onto a condensed + * phase species. + */ + +/* + * $Author$ + * $Date$ + * $Revision$ + */ +/* + * Copywrite 2006 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000, with Sandia Corporation, the U.S. Government + * retains certain rights in this software. + */ + +#ifndef CT_IDEALSOLIDSOLNPHASE_H +#define CT_IDEALSOLIDSOLNPHASE_H + +#include "mix_defs.h" +#include "ThermoPhase.h" +#include "importCTML.h" +#include "SpeciesThermo.h" + + +namespace Cantera { + + const int cIdealSolidSolnPhase0 = 5010; + const int cIdealSolidSolnPhase1 = 5011; + const int cIdealSolidSolnPhase2 = 5012; + + /** + * Class IdealSolidSolnPhase represents a condensed phase ideal + * solution compound. 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. + * Each species does, however, have constant but distinct partial + * molar volumes equal to their pure species molar volumes. + * 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 + * is supplied in the constructor and in the XML file. + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * The value and form of the generalized concentration will affect + * reaction rate constants involving species in this phase. + * + * @ingroup thermoprops + */ + 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 + * is supplied in the constructor or read from the xml data file. + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + */ + IdealSolidSolnPhase(int formCG=0); + + /** + * Constructor for IdealSolidSolnPhase. + * + * 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 + * is supplied in the constructor or read from the xml data file. + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * @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. + */ + IdealSolidSolnPhase(string infile, string id="", int formCG=0); + + + /** + * Constructor for IdealSolidSolnPhase. + * 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 + * is supplied in the constructor and/or read from the data file. + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * @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. + * + */ + IdealSolidSolnPhase(XML_Node& root, string id="", int formCG=0); + + /* + * Copy Constructor + */ + IdealSolidSolnPhase(const IdealSolidSolnPhase &); + + /* + * 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) + */ + virtual ThermoPhase* duplMyselfAsThermoPhase(); + + + virtual ~IdealSolidSolnPhase() {} + + /** + * Equation of state flag. Returns a value depending upon the value of + * m_formGC, which is defined at instantiation. + */ + virtual int eosType() const; + + /** + * @name 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 + */ + virtual doublereal enthalpy_mole() const; + + /** + * Molar internal 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 and + * zero isothermal compressibility: + * \f[ + * \hat u(T,X) = \hat h(T,P,X) - 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 + */ + virtual doublereal intEnergy_mole() const; + + /** + * 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 + */ + virtual doublereal entropy_mole() const; + + /** + * 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 + */ + virtual doublereal gibbs_mole() const; + + /** + * 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 + */ + virtual doublereal cp_mole() const; + + /** + * Molar heat capacity at constant volume 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_v(T,P) = \hat c_p(T,P) \f] + * The two heat capacities are equal. + */ + virtual doublereal cv_mole() const { + return cp_mole(); + } + + //@} + /** @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 + * an independent variable. Instead, the pressure is used as the + * independent variable. Functions which try to set the thermodynamic + * state by calling setDensity() may cause an exception to be + * thrown. + */ + //@{ + + /** + * Pressure. Units: Pa. + * For this incompressible system, we return the internally storred + * independent value of the pressure. + */ + virtual doublereal pressure() const { + return m_Pcurrent; + } + + /** + * 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. + */ + virtual void setPressure(doublereal p); + + /** + * 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. + */ + void calcDensity(); + + /** + * Overwritten setDensity() function is necessary because the + * density is not an indendent 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 + */ + virtual void setDensity(doublereal rho); + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an independent variable. + * + * This function will now throw an error condition. + * + * NOTE: This is virtual function that overwrites the State.h + * class + */ + virtual void setMolarDensity(doublereal rho); + + virtual void setMoleFractions(doublereal *x); + virtual void setMoleFractions_NoNorm(doublereal *x); + virtual void setMassFractions(doublereal *y); + virtual void setMassFractions_NoNorm(doublereal *y); + virtual void setConcentrations(doublereal *c); + + + + //@} + + /** + * @name Chemical Potentials and Activities ----------------------------------------- + * + * The activity \f$a_k\f$ of a species in solution is + * related to the chemical potential by + * \f[ + * \mu_k(T,P,X_k) = \mu_k^0(T,P) + * + \hat R T \log a_k. + * \f] + * The quantity \f$\mu_k^0(T,P)\f$ is + * the standard state chemical potential at unit activity. + * It may depend on the pressure and the temperature. However, + * it may not depend on the mole fractions of the species + * in the solid solution. + * + * The activities are related to the generalized + * concentrations, \f$\tilde C_k\f$, and standard + * concentrations, \f$C^0_k\f$, by the following formula: + * + * \f[ + * a_k = \frac{\tilde C_k}{C^0_k} + * \f] + * The generalized concentrations are used in the kinetics classes + * to describe the rates of progress of reactions involving the + * species. Their formulation depends upons the specification + * of the rate constants for reaction, especially the units used + * in specifying the rate constants. The bridge between the + * thermodynamic equilibrium expressions that use a_k and the + * kinetics expressions which use the generalized concentrations + * is provided by the multiplicative factor of the + * standard concentrations. + * @{ + */ + + /** + * 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 dividied by the standard concentration is also + * equal to the activity of species. + * + * For this implentation 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-3. 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: + * + * + * + * + * + * + *
m_formGC GeneralizedConc StandardConc
0 X_k 1.0
1 X_k / V_k 1.0 / V_k
2 X_k / V_N 1.0 / V_N
+ * + * HKM Note: We have absorbed the pressure dependence of the pures 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 awkwards exp[] format. + * + * @param c[] Pointer to array of doubles of length m_kk, which on exit + * will contain the generalized concentrations. + */ + 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-3. + * + * @param k Species number: this is a require parameter, + * a change from the ThermoPhase base class, where it was + * an optional parameter. + */ + virtual doublereal standardConcentration(int k) const; + + /** + * The reference (ie 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. (clone of the standard concentration -> + * suggest changing the name). This is the inverse of the species molar + * volume. + */ + virtual doublereal referenceConcentration(int k) const; + + /** + * Returns the log of the standard concentration of the kth species + * + * @param k Species number: this is a require parameter, + * a change from the ThermoPhase base class, where it was + * an optional parameter. + */ + virtual doublereal logStandardConc(int k) const; + + /** + * 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 cIdealSolidSolnPhase0, the default + * kmol/m3 holds for standard concentration units. For + * cIdealSolidSolnPhase0 type, the standard concentrtion is + * unitless. + */ + virtual void getUnitsStandardConc(double *uA, int k = 0, + int sizeUA = 6); + + /** + * Get the array of species activity coefficients + */ + virtual void getActivityCoefficients(doublereal * ac) const; + + /** + * Get the species chemical potentials. Units: J/kmol. + * + * This function returns a vector of chemical potentials of the + * species in solution. + * \f[ + * \mu_k = \mu^{ref}_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^{ref}_k(T) + V_k * (p - p_o)\f$ + */ + virtual void getChemPotentials(doublereal* mu) const; + + /** + * Get the array of non-dimensional species solution + * chemical potentials at the current T and P + * \f$\mu_k / \hat R T \f$. + * \f[ + * \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k + RT ln(X_k) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + */ + virtual void getChemPotentials_RT(doublereal* mu) const; + + //@} + /// @name Partial Molar Properties of the Solution ----------------------------- + //@{ + + /** + * Returns an array of partial molar enthalpies for the species + * in the mixture. + * Units (J/kmol) + * For this phase, the partial molar enthalpies are equal to the + * pure species enthalpies + * \f[ + * \bar h_k(T,P) = \hat h^{ref}_k(T) + (P - P_{ref}) \hat V^0_k + * \f] + * The reference-state pure-species enthalpies, \f$ \hat h^{ref}_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 + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol/K. + * 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^{ref}_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 + */ + virtual void getPartialMolarEntropies(doublereal* sbar) const; + + /** + * Returns an array of partial molar Heat Capacities at constant + * pressure of the species in the + * solution. Units: J/kmol/K. + * For this phase, the partial molar heat capacities are equal + * to the standard state heat capacities. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + /** + * 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. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //@} + /// @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 + * \f$ \mu^0_k(T,P) \f$. + * We define these here as the chemical potentials of the pure + * species at the temperature and pressure of the solution. + * This function is used in the evaluation of the + * equilibrium constant Kc. Therefore, Kc will also depend + * on T and P. This is the norm for liquid and solid systems. + * + * units = J / kmol + */ + virtual void getStandardChemPotentials(doublereal* mu0) const { + getPureGibbs(mu0); + } + + /** + * Get the array of nondimensional Enthalpy functions for the + * standard state species + * at the current T and P 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 k. + * \f$ h^{ref}_k(T)\f$ is the enthalpy of the pure + * species k at the reference pressure, \f$P_{ref}\f$. + * + * @param hrt Vector of length m_kk, which on return hrt[k] + * will contain the nondimensional + * standard state enthalpy of species k. + */ + void 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 for species k. + */ + void getEntropy_R(doublereal* sr) 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 k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k 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. + */ + virtual void getGibbs_RT(doublereal* grt) const; + + /** + * Get the Gibbs functions for the pure species + * at the current T and P 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 k. + * \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure + * species k at the reference pressure, \f$P_{ref}\f$. + */ + 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. + */ + 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. + * \f[ + * Cp^0_k(T,P) = Cp^{ref}_k(T) + * \f] + * where \f$V_k\f$ is the molar volume of pure species k. + * \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity + * of species k 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 getCp_R(doublereal* cpr) const; + + /** + * Get the molar volumes of each species in their standard + * states at the current + * T and P of the solution. + * units = m^3 / kmol + */ + virtual void getStandardVolumes(doublereal *vol) const; + + + //@} + /// @name 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. + */ + virtual void getEnthalpy_RT_ref(doublereal *hrt) const; + + /** + * Returns the vector of nondimensional + * enthalpies of the reference state at the current temperature + * of the solution and the reference pressure for the species. + */ + virtual void getGibbs_RT_ref(doublereal *grt) const; + + /** + * 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 + */ + 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. + */ + virtual void getEntropy_R_ref(doublereal *er) const; + + /** + * 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. + */ + virtual void getIntEnergy_RT_ref(doublereal *urt) 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. + */ + virtual void getCp_R_ref(doublereal *cprt) 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 array_fp& enthalpy_RT_ref() 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 array_fp& gibbs_RT_ref() const { + _updateThermo(); + return m_g0_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 array_fp& expGibbs_RT_ref() 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 array_fp& entropy_R_ref() 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 array_fp& cp_R_ref() const { + _updateThermo(); + return m_cp0_R; + } + + virtual void setPotentialEnergy(int k, doublereal pe) { + m_pe[k] = pe; + _updateThermo(); + } + + virtual doublereal potentialEnergy(int k) const { + return m_pe[k]; + } + //@} + /// @name Utility Functions ----------------------------------------------- + //@{ + + + + /** + * Initialization of an IdealSolidSolnPhase phase using an + * xml file + * + * This routine is a precursor to constructPhaseXML(XML_Node*) + * routine, which does most of the work. + * + * @param infile XML file containing the description of the + * phase + * + * @param id Optional parameter identifying the name of the + * phase. If none is given, the first XML + * phase element will be used. + */ + virtual void constructPhaseFile(string infile, string id=""); + + /** + * Import and initialize an IdealSolidSolnPhase phase + * specification in an XML tree into the current object. + * Here we read an XML description of the phase. + * We import descriptions of the elements that make up the + * species in a phase. + * We import information about the species, including their + * reference state thermodynamic polynomials. We then freeze + * the state of the species. + * This routine calls importPhase() to do most of its work. + * Then, importPhase() calls initThermoXML() to finish + * off the work. + * + + * @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 constructPhaseXML(XML_Node& phaseNode, string id=""); + + /** + * Initialization of an IdealSolidSolnPhase phase: + * Note this function is pretty much useless because it doesn't + * get the xml tree passed to it. Suggest a change. + */ + 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. + * + * @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, 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$. + * @param t temperature in K. + * @param work. Temporary work space. Must be dimensioned at least + * as large as the number of species. + * + */ + virtual void setToEquilState(const doublereal* lambda_RT); + + + /** + * Report the molar volume of species k + * + * units - \f$ m^3 kmol^-1 \f$ + */ + double speciesMolarVolume(int k) const; + + /** + * Fill in a return vector containing the species molar volumes + * units - \f$ m^3 kmol^-1 \f$ + */ + void getSpeciesMolarVolumes(double *smv) const; + //@} + + 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 + */ + int m_formGC; + /** + * m_mm = Number of distinct elements defined in species in this + * phase + */ + int m_mm; + + /** + * Maximum temperature that this phase can accurately describe + * the thermodynamics. + */ + doublereal m_tmin; + + /** + * Minimum temperature that this phase can accurately describe + * the thermodynamics. + */ + doublereal m_tmax; + /** + * Value of the reference pressure for all species in this phase. + * The T dependent polynomials are evaluated at the reference + * pressure. Note, because this is a single value, all species + * are required to have the same reference pressure. + */ + doublereal m_Pref; + + /** + * m_Pcurrent = The current pressure + * Since the density isn't a function of pressure, but only of the + * mole fractions, we need to independently specify the pressure. + * The density variable which is inherited as part of the State class, + * m_dens, is always kept current whenever T, P, or X[] change. + */ + doublereal m_Pcurrent; + + /** + * Species molar volume \f$ m^3 kmol^-1 \f$ + */ + array_fp m_speciesMolarVolume; + + /** + * Value of the temperature at which the thermodynamics functions + * for the reference state of the species were last evaluated. + */ + mutable doublereal m_tlast; + + /** + * Vector containing the species reference enthalpies at T = m_tlast + */ + mutable array_fp m_h0_RT; + + /** + * Vector containing the species reference constant pressure + * heat capacities at T = m_tlast + */ + mutable array_fp m_cp0_R; + + /** + * Vector containing the species reference Gibbs functions + * at T = m_tlast + */ + mutable array_fp m_g0_RT; + + /** + * Vector containing the species reference entropies + * at T = m_tlast + */ + mutable array_fp m_s0_R; + + /** + * Vector containing the species reference exp(-G/RT) functions + * at T = m_tlast + */ + mutable array_fp m_expg0_RT; + + /** + * Vector of potential energies for the species. + */ + mutable array_fp m_pe; + + /** + * Temporary array used in equilibrium calculations + */ + mutable array_fp m_pp; + + private: + /// @name Utility Functions ------------------------------------------ + //@{ + /** + * 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 _updateThermo() const; + + /** + * This internal function adjusts the lengths of arrays + */ + void initLengths(); + + //@} + }; +} + +#endif + + + + + diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in index bea93ec37..eb93ebbf4 100644 --- a/Cantera/src/thermo/Makefile.in +++ b/Cantera/src/thermo/Makefile.in @@ -19,16 +19,18 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) # Extended Cantera Thermodynamics Object Files CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \ - MolalityVPSSTP.o VPStandardStateTP.o + MolalityVPSSTP.o VPStandardStateTP.o \ + IdealSolidSolnPhase.o CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \ - MolalityVPSSTP.h VPStandardStateTP.h + MolalityVPSSTP.h VPStandardStateTP.h \ + IdealSolidSolnPhase.h CXX_INCLUDES = -I.. @CXX_INCLUDES@ LIB = @buildlib@/libcaThermo.a DEPENDS = $(CATHERMO_OBJ:.o=.d) -all: $(LIB) +all: $(LIB) .depends @(@INSTALL@ -d $(INCDIR)) @(for lh in $(CATHERMO_H) ; do \ $(INSTALL_TSC) "$${lh}" $(INCDIR) ; \ @@ -57,14 +59,16 @@ clean: $(RM) $(LIB) ; \ echo "$(RM) $(LIB)" ; \ fi) - $(RM) *.o *~ .depends + $(RM) *.o *~ .depends *.d (if test -d SunWS_cache ; then \ $(RM) -rf SunWS_cache ; \ fi ) -depends: $(DEPENDS) +depends: + @MAKE@ .depends + +.depends: $(DEPENDS) cat *.d > .depends - $(RM) $(DEPENDS) TAGS: etags *.h *.cpp diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp index aef04292a..4ccac77bf 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.cpp +++ b/Cantera/src/thermo/MolalityVPSSTP.cpp @@ -26,14 +26,18 @@ namespace Cantera { * Default constructor. * * This doesn't do much more than initialize constants with - * default values for water at 25C. + * default values for water at 25C. Water molecular weight + * comes from the default elements.xml file. It actually + * differs slightly from the IAPWS95 value of 18.015268. However, + * density conservation and therefore element conservation + * is the more important principle to follow. */ MolalityVPSSTP::MolalityVPSSTP() : VPStandardStateTP(), m_indexSolvent(0), - m_weightSolvent(18.0), + m_weightSolvent(18.01528), m_xmolSolventMIN(0.01), - m_Mnaught(18.0E-3) + m_Mnaught(18.01528E-3) { } @@ -109,7 +113,8 @@ namespace Cantera { */ void MolalityVPSSTP::setSolvent(int k) { if (k < 0 || k >= m_kk) { - throw CanteraError("MolalityVPSSTP::setSolute ", "trouble"); + throw CanteraError("MolalityVPSSTP::setSolute ", + "bad value"); } m_indexSolvent = k; m_weightSolvent = molecularWeight(k); @@ -144,6 +149,33 @@ namespace Cantera { return m_xmolSolventMIN; } + /** + * calcMolalities(): + * We calculate the vector of molalities of the species + * in the phase and store the result internally: + * \f[ + * m_i = (n_i) / (1000 * M_o * n_{o,p}) + * \f] + * where + * - \f$ M_o \f$ is the molecular weight of the solvent + * - \f$ n_o \f$ is the mole fraction of the solvent + * - \f$ n_i \f$ is the mole fraction of the solute. + * - \f$ n_{o,p} = max (n_{o, min}, n_o) \f$ + * - \f$ n_{o,min} \f$ = minimum mole fraction of solvent allowed + * in the denominator. + */ + void MolalityVPSSTP::calcMolalities() const { + getMoleFractions(DATA_PTR(m_molalities)); + double xmolSolvent = m_molalities[m_indexSolvent]; + if (xmolSolvent < m_xmolSolventMIN) { + xmolSolvent = m_xmolSolventMIN; + } + double denomInv = 1.0/ (m_Mnaught * xmolSolvent); + for (int k = 0; k < m_kk; k++) { + m_molalities[k] *= denomInv; + } + } + /** * getMolalities(): * We calculate the vector of molalities of the species @@ -160,19 +192,10 @@ namespace Cantera { * in the denominator. */ void MolalityVPSSTP::getMolalities(doublereal * const molal) const { - getMoleFractions(molal); - double xmolSolvent = molal[m_indexSolvent]; - if (xmolSolvent < m_xmolSolventMIN) { - xmolSolvent = m_xmolSolventMIN; - } - double denomInv = 1.0/ - (m_Mnaught * xmolSolvent); - for (int k = 0; k < m_kk; k++) { - molal[k] *= denomInv; - } - for (int k = 0; k < m_kk; k++) { - m_molalities[k] = molal[k]; - } + calcMolalities(); + for (int k = 0; k < m_kk; k++) { + molal[k] = m_molalities[k]; + } } /** @@ -220,7 +243,7 @@ namespace Cantera { * the molalities from the mole fractions that we * just obtained. */ - getMolalities(DATA_PTR(m_molalities)); + calcMolalities(); } /* @@ -309,7 +332,7 @@ namespace Cantera { * calculate the molalities again and store it in * this object. */ - getMolalities(DATA_PTR(m_molalities)); + calcMolalities(); } /* @@ -327,16 +350,6 @@ namespace Cantera { setMolalitiesByName(xx); } - - /* - * Update the internal array that contains the molalities of the - * species. - */ - void MolalityVPSSTP::updateMolalities() const { - getMolalities(DATA_PTR(m_molalities)); - } - - /* * ------------ Molar Thermodynamic Properties ---------------------- @@ -401,14 +414,22 @@ namespace Cantera { * * Calculate the osmotic coefficient of the solvent. Note there * are lots of definitions of the osmotic coefficient floating - * around. We use the one defined in the Pitzer paper: + * around. We use the one defined in the Pitzer's book: + * (Activity Coeff in Electrolyte Solutions, K. S. Pitzer + * CRC Press, Boca Raton, 1991, p. 85, Eqn. 28). * * Definition: - * - sum(m_i) * M0 * oc = ln(activity_solvent) + * - sum(m_i) * Mnaught * oc = ln(activity_solvent) */ doublereal MolalityVPSSTP::osmoticCoefficient() const { + /* + * First, we calculate the activities all over again + */ vector_fp act(m_kk); getActivities(DATA_PTR(act)); + /* + * Then, we calculate the sum of the solvent molalities + */ double sum = 0; for (int k = 0; k < m_kk; k++) { if (k != m_indexSolvent) { @@ -483,6 +504,17 @@ namespace Cantera { } } + /** + * Set the temperature (K), pressure (Pa), and molalities + * (gmol kg-1) of the solutes + */ + void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, + const doublereal * const molalities) { + setMolalities(molalities); + setTemperature(t); + setPressure(p); + } + /** * @internal Initialize. This method is provided to allow * subclasses to perform any initialization required after all @@ -497,10 +529,47 @@ namespace Cantera { * @see importCTML.cpp */ void MolalityVPSSTP::initThermo() { - VPStandardStateTP::initThermo(); - m_molalities.resize(m_kk); + initLengths(); + VPStandardStateTP::initThermo(); + + /* + * The solvent defaults to species 0 + */ + setSolvent(0); } + void MolalityVPSSTP::initLengths() { + int m_kk = nSpecies(); + m_molalities.resize(m_kk); + } + + /** + * initThermoXML() (virtual from ThermoPhase) + * 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. + */ + void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, string id) { + + initLengths(); + /* + * The solvent defaults to species 0 + */ + setSolvent(0); + + VPStandardStateTP::initThermoXML(phaseNode, id); + } + + } diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h index fccf6772f..dcdc2ad57 100644 --- a/Cantera/src/thermo/MolalityVPSSTP.h +++ b/Cantera/src/thermo/MolalityVPSSTP.h @@ -119,6 +119,12 @@ namespace Cantera { */ doublereal moleFSolventMin() const; + /** + * Calculates the molality of all species and + * stores the result internally. + */ + void calcMolalities() const; + /** * getMolalities() * This function will return the molalities of the @@ -131,7 +137,7 @@ namespace Cantera { void setMolalities(const doublereal * const molal); void setMolalitiesByName(compositionMap& xMap); void setMolalitiesByName(const string &); - void updateMolalities() const; + /** * @} * @name Mechanical Properties @@ -415,6 +421,34 @@ namespace Cantera { * @see importCTML.cpp */ virtual void initThermo(); + + + /** + * 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. + */ + void initThermoXML(XML_Node& phaseNode, string id); + + /** + * Set the temperature (K), pressure (Pa), and molalities + * (gmol kg-1) of the solutes + */ + void setState_TPM(doublereal t, doublereal p, + const doublereal * const molalities); + + + private: + void initLengths(); protected: diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp index 870817d2b..490771cfa 100644 --- a/Cantera/src/thermo/VPStandardStateTP.cpp +++ b/Cantera/src/thermo/VPStandardStateTP.cpp @@ -263,15 +263,45 @@ namespace Cantera { * added. */ void VPStandardStateTP::initThermo() { - ThermoPhase::initThermo(); - m_kk = nSpecies(); - int leng = m_kk; - m_h0_RT.resize(leng); - m_g0_RT.resize(leng); - m_cp0_R.resize(leng); - m_s0_R.resize(leng); + initLengths(); + ThermoPhase::initThermo(); } + /** + * Initialize the internal lengths. + * (this is not a virtual function) + */ + void VPStandardStateTP::initLengths() { + m_kk = nSpecies(); + int leng = m_kk; + m_h0_RT.resize(leng); + m_g0_RT.resize(leng); + m_cp0_R.resize(leng); + m_s0_R.resize(leng); + } + + /** + * 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 VPStandardStateTP::initThermoXML(XML_Node& phaseNode, string id) { + VPStandardStateTP::initLengths(); + ThermoPhase::initThermoXML(phaseNode, id); + } + /** * void _updateRefStateThermo() (private, const) * diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h index 90331ee8d..19c5b26bf 100644 --- a/Cantera/src/thermo/VPStandardStateTP.h +++ b/Cantera/src/thermo/VPStandardStateTP.h @@ -380,7 +380,26 @@ namespace Cantera { * @see importCTML.cpp */ virtual void initThermo(); - + + /** + * 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. + */ + void initThermoXML(XML_Node& phaseNode, string id); + + private: + void initLengths(); + protected: /* * The last temperature at which the reference thermodynamic