From 93fcf181b1ed4e186c74432915142022d0283e7e Mon Sep 17 00:00:00 2001 From: Victor Brunini Date: Fri, 28 Feb 2014 00:38:09 +0000 Subject: [PATCH] Initial pass at implementing Maskell solid solution ThermoPhase. Based on Maskell, Shaw, and Tye, Electrochimica Acta 28(2) 225-230 1983. Includes unit tests checking calculation of the chemical potential values and parsing of the mixing enthalpy parameter from an XML file. --- .../cantera/thermo/MaskellSolidSolnPhase.h | 318 ++++++++++++++++++ include/cantera/thermo/mix_defs.h | 2 + src/thermo/MaskellSolidSolnPhase.cpp | 284 ++++++++++++++++ src/thermo/ThermoFactory.cpp | 11 +- test/data/MaskellSolidSolnPhase_nohmix.xml | 57 ++++ test/data/MaskellSolidSolnPhase_valid.xml | 58 ++++ test/thermo/MaskellSolidSolnPhase_Test.cpp | 106 ++++++ 7 files changed, 833 insertions(+), 3 deletions(-) create mode 100644 include/cantera/thermo/MaskellSolidSolnPhase.h create mode 100644 src/thermo/MaskellSolidSolnPhase.cpp create mode 100644 test/data/MaskellSolidSolnPhase_nohmix.xml create mode 100644 test/data/MaskellSolidSolnPhase_valid.xml create mode 100644 test/thermo/MaskellSolidSolnPhase_Test.cpp diff --git a/include/cantera/thermo/MaskellSolidSolnPhase.h b/include/cantera/thermo/MaskellSolidSolnPhase.h new file mode 100644 index 000000000..007a38c41 --- /dev/null +++ b/include/cantera/thermo/MaskellSolidSolnPhase.h @@ -0,0 +1,318 @@ +/** + * @file MaskellSolidSolnPhase.h Header file for a solid solution model + * following Maskell, Shaw, and Tye. Electrochimica Acta 1982 + * + * This class inherits from the Cantera class ThermoPhase and implements a + * non-ideal solid solution model with incompressible thermodynamics. + */ +/* + * Copyright 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_MASKELLSOLIDSOLNPHASE_H +#define CT_MASKELLSOLIDSOLNPHASE_H + +#include "mix_defs.h" +#include "ThermoPhase.h" +#include "ThermoFactory.h" +#include "SpeciesThermo.h" + +namespace Cantera +{ + +/** + * Class MaskellSolidSolnPhase represents a condensed phase + * non-ideal solution with 2 species following the thermodynamic + * model described in Maskell, Shaw, and Tye, Manganese Dioxide Electrode -- IX, + * Electrochimica Acta 28(2) pp 231-235, 1983. + * + * @ingroup thermoprops + */ +class MaskellSolidSolnPhase : public ThermoPhase +{ +public: + /** + * 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. + * + * @param formCG This parameter initializes the #m_formGC variable. + */ + MaskellSolidSolnPhase(); + + //! Construct and initialize an MaskellSolidSolnPhase ThermoPhase object + //! directly from an XML database + /*! + * @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. + */ + //MaskellSolidSolnPhase(XML_Node& root, const std::string& id=""); + + //! Copy Constructor + MaskellSolidSolnPhase(const MaskellSolidSolnPhase&); + + //! Assignment operator + MaskellSolidSolnPhase& operator=(const MaskellSolidSolnPhase&); + + /*! + * Base Class Duplication Function + * + * Given a pointer to ThermoPhase, this function can duplicate the object. + */ + virtual ThermoPhase* duplMyselfAsThermoPhase() const; + + //! @name Molar Thermodynamic Properties of the Solution + //! @{ + /** + * Molar enthalpy of the solution. Units: J/kmol. + */ + virtual doublereal enthalpy_mole() const; + + /** + * Molar entropy of the solution. Units: J/kmol/K. + */ + virtual doublereal entropy_mole() const; + + /** + * Molar heat capacity at constant pressure of the solution. + * Units: J/kmol/K. + */ + //virtual doublereal cp_mole() const; + + /** + * Molar heat capacity at constant volume of the solution. + * Units: J/kmol/K. + */ + //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 stored + * 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. + * + * @param p Input Pressure (Pa) + */ + virtual void setPressure(doublereal p); + + /** + * Overwritten setDensity() function is necessary because the + * density is not an independent variable. + * + * This function will now throw an error condition + * + * @internal May have to adjust the strategy here to make + * the eos for these materials slightly compressible, in order + * to create a condition where the density is a function of + * the pressure. + * + * @param rho Input density + */ + virtual void setDensity(const doublereal rho); + + /** + * Overwritten setMolarDensity() function is necessary because the + * density is not an independent variable. + * + * This function will now throw an error condition. + * + * @param rho Input Density + */ + virtual void setMolarDensity(const doublereal rho); + + //@} + + /** + * @name Chemical Potentials and Activities + * @{ + */ + + //! Get the array of species activity coefficients + /*! + * @param ac output vector of activity coefficients. Length: m_kk + */ + virtual void getActivityCoefficients(doublereal* ac) const; + + /** + * Get the species chemical potentials. Units: J/kmol. + * + * @param mu Output vector of chemical potentials. + */ + virtual void getChemPotentials(doublereal* mu) const; + + /** + * Get the array of non-dimensional species solution + * chemical potentials at the current T and P + * + * @param mu Output vector of dimensionless chemical potentials. Length = m_kk. + */ + 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) + * + * @param hbar Output vector containing partial molar enthalpies. + * Length: m_kk. + */ + virtual void getPartialMolarEnthalpies(doublereal* hbar) const; + + /** + * Returns an array of partial molar entropies of the species in the + * solution. Units: J/kmol/K. + * + * @param sbar Output vector containing partial molar entropies. + * Length: m_kk. + */ + 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. + * + * @param cpbar Output vector of partial heat capacities. Length: m_kk. + */ + virtual void getPartialMolarCp(doublereal* cpbar) const; + + /** + * returns an array of partial molar volumes of the species + * in the solution. Units: m^3 kmol-1. + * + * @param vbar Output vector of partial molar volumes. Length: m_kk. + */ + virtual void getPartialMolarVolumes(doublereal* vbar) const; + + //! Get the Gibbs functions for the standard + //! state of the species at the current T and P of the solution + /*! + * Units are Joules/kmol + * @param gpure Output vector of standard state gibbs free energies + * Length: m_kk. + */ + virtual void getPureGibbs(doublereal* gpure) const; + + //! Get the array of chemical potentials at unit activity for the species + //! at their standard states at the current T and P of the solution. + /*! + * These are the standard state chemical potentials \f$ \mu^0_k(T,P) + * \f$. The values are evaluated at the current + * temperature and pressure of the solution + * + * @param mu Output vector of chemical potentials. + * Length: m_kk. + */ + virtual void getStandardChemPotentials(doublereal* mu) const; + + //@} + /// @name Utility Functions + //@{ + + /** + * @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, const std::string& id); + + + void set_h_mix(const doublereal hmix) { h_mixing = hmix; } +protected: + /** + * 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. + */ + doublereal m_Pcurrent; + + /** + * Function to call through to m_spthermo->update and fill m_h0_RT, + * m_cp0_R, m_g0_RT, m_s0_R. + */ + void _updateThermo() const; + + /** + * 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 vector_fp m_h0_RT; + + /** + * Vector containing the species reference constant pressure + * heat capacities at T = m_tlast + */ + mutable vector_fp m_cp0_R; + + //! Vector containing the species reference Gibbs functions at T = m_tlast + mutable vector_fp m_g0_RT; + + //! Vector containing the species reference entropies at T = m_tlast + mutable vector_fp m_s0_R; + + //! Value of the enthalpy change on mixing due to protons changing from type B to type A configurations. + doublereal h_mixing; + +private: + // Functions to calculate some of the pieces of the mixing terms. + doublereal s() const; + doublereal fm(const doublereal r) const; + doublereal p(const doublereal r) const; +}; +} + +#endif diff --git a/include/cantera/thermo/mix_defs.h b/include/cantera/thermo/mix_defs.h index 807649524..45dfbf3cb 100644 --- a/include/cantera/thermo/mix_defs.h +++ b/include/cantera/thermo/mix_defs.h @@ -63,6 +63,8 @@ const int cFixedChemPot = 70; /// Constant partial molar volume solution IdealSolidSolnPhase.h const int cIdealSolidSolnPhase = 5009; +const int cMaskellSolidSolnPhase = 5010; + //! HMW - Strong electrolyte using the Pitzer formulation const int cHMW = 40; diff --git a/src/thermo/MaskellSolidSolnPhase.cpp b/src/thermo/MaskellSolidSolnPhase.cpp new file mode 100644 index 000000000..032f102ff --- /dev/null +++ b/src/thermo/MaskellSolidSolnPhase.cpp @@ -0,0 +1,284 @@ +/** + * @file MaskellSolidSolnPhase.cpp Implementation file for an ideal solid + * solution model with incompressible thermodynamics (see \ref + * thermoprops and \link Cantera::MaskellSolidSolnPhase + * MaskellSolidSolnPhase\endlink). + */ +/* + * Copyright 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 "cantera/thermo/MaskellSolidSolnPhase.h" + +#include "cantera/base/stringUtils.h" + +#include +#include + +namespace Cantera +{ + +MaskellSolidSolnPhase::MaskellSolidSolnPhase() : + ThermoPhase(), + m_Pref(OneAtm), + m_Pcurrent(OneAtm), + m_tlast(-1), + m_h0_RT(2), + m_cp0_R(2), + m_g0_RT(2), + m_s0_R(2), + h_mixing(0.) +{ +} + +MaskellSolidSolnPhase::MaskellSolidSolnPhase(const MaskellSolidSolnPhase& b) +{ + *this = b; +} + +MaskellSolidSolnPhase& MaskellSolidSolnPhase:: +operator=(const MaskellSolidSolnPhase& b) +{ + if (this != &b) { + ThermoPhase::operator=(b); + } + return *this; +} + +ThermoPhase* MaskellSolidSolnPhase::duplMyselfAsThermoPhase() const +{ + return new MaskellSolidSolnPhase(*this); +} + +/******************************************************************** + * Molar Thermodynamic Properties of the Solution + ********************************************************************/ + +doublereal MaskellSolidSolnPhase:: +enthalpy_mole() const +{ + _updateThermo(); + const doublereal h0 = GasConstant * temperature() * mean_X(&m_h0_RT[0]); + const doublereal r = moleFraction(0); + const doublereal fmval = fm(r); + return h0 + r * fmval * h_mixing; +} + +doublereal xlogx(doublereal x) +{ + return x * std::log(x); +} + +doublereal MaskellSolidSolnPhase::entropy_mole() const +{ + _updateThermo(); + const doublereal s0 = GasConstant * mean_X(&m_s0_R[0]); + const doublereal r = moleFraction(0); + const doublereal fmval = fm(r); + const doublereal rfm = r * fmval; + return s0 + GasConstant * (xlogx(1-rfm) - xlogx(rfm) - xlogx(1-r-rfm) - xlogx((1-fmval)*r) - xlogx(1-r) - xlogx(r)); +} + +/******************************************************************** + * Mechanical Equation of State + ********************************************************************/ + +void MaskellSolidSolnPhase:: +setDensity(const doublereal rho) +{ + /* + * Unless the input density is exactly equal to the density + * calculated and stored 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("MaskellSolidSolnPhase::setDensity", + "Density is not an independent variable"); + } +} + +void MaskellSolidSolnPhase::setPressure(doublereal p) +{ + m_Pcurrent = p; +} + +void MaskellSolidSolnPhase::setMolarDensity(const doublereal n) +{ + throw CanteraError("MaskellSolidSolnPhase::setMolarDensity", + "Density is not an independent variable"); +} + +/******************************************************************** + * Chemical Potentials and Activities + ********************************************************************/ + +void MaskellSolidSolnPhase:: +getActivityCoefficients(doublereal* ac) const +{ +} + +void MaskellSolidSolnPhase:: +getChemPotentials(doublereal* mu) const +{ + _updateThermo(); + const doublereal r = moleFraction(0); + const doublereal pval = p(r); + const doublereal fmval = fm(r); + const doublereal rfm = r * fmval; + const doublereal RT = GasConstant * temperature(); + const doublereal muDelta = -pval * h_mixing - GasConstant * temperature() + * std::log( (std::pow(1 - rfm, pval) * std::pow(rfm, pval) * std::pow(r - rfm, 1 - pval) * r) / + (std::pow(1 - r - rfm, 1 + pval) * (1 - r)) ); + mu[0] = RT * m_g0_RT[0] + muDelta; + mu[1] = RT * m_g0_RT[1] - muDelta; +} + +void MaskellSolidSolnPhase:: +getChemPotentials_RT(doublereal* mu) const +{ + const doublereal invRT = 1.0 / (GasConstant * temperature()); + getChemPotentials(mu); + for(unsigned sp=0; sp < m_kk; ++sp) + { + mu[sp] *= invRT; + } +} + +/******************************************************************** + * Partial Molar Properties + ********************************************************************/ + +void MaskellSolidSolnPhase::getPartialMolarEnthalpies(doublereal* hbar) const +{ +} + +void MaskellSolidSolnPhase:: +getPartialMolarEntropies(doublereal* sbar) const +{ +} + +void MaskellSolidSolnPhase:: +getPartialMolarCp(doublereal* cpbar) const +{ +} + +void MaskellSolidSolnPhase:: +getPartialMolarVolumes(doublereal* vbar) const +{ +} + +void MaskellSolidSolnPhase:: +getPureGibbs(doublereal* gpure) const +{ + _updateThermo(); + const doublereal RT = GasConstant * temperature(); + for(unsigned sp=0; sp < m_kk; ++sp) + { + gpure[sp] = RT * m_g0_RT[sp]; + } +} + +void MaskellSolidSolnPhase:: +getStandardChemPotentials(doublereal* mu) const +{ + // What is the difference between this and getPureGibbs? IdealSolidSolnPhase gives the same for both + getPureGibbs(mu); +} + +/********************************************************************* + * Utility Functions + *********************************************************************/ +void MaskellSolidSolnPhase::initThermoXML(XML_Node& phaseNode, const std::string& id_) +{ + std::string subname = "MaskellSolidSolnPhase::initThermoXML"; + if (id_.size() > 0) { + std::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"); + std::string mStringa = thNode.attrib("model"); + std::string mString = lowercase(mStringa); + if (mString != "maskellsolidsolnphase") { + throw CanteraError(subname.c_str(), + "Unknown thermo model: " + mStringa); + } + + + /* + * Parse the enthalpy of mixing constant + */ + if (thNode.hasChild("h_mix")) { + XML_Node& scNode = thNode.child("h_mix"); + set_h_mix(fpValue(scNode.value())); + } else { + throw CanteraError(subname.c_str(), + "Mixing enthalpy parameter not specified."); + } + } else { + throw CanteraError(subname.c_str(), + "Unspecified thermo model"); + } + + + // Confirm that the phase only contains 2 species + if( m_kk != 2 ) + { + throw CanteraError( subname.c_str(), "MaskellSolidSolution model requires exactly 2 species."); + } + + /* + * Call the base initThermo, which handles setting the initial + * state. + */ + ThermoPhase::initThermoXML(phaseNode, id_); +} + +void MaskellSolidSolnPhase::_updateThermo() const +{ + assert(m_kk == 2); + 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; + for (size_t k = 0; k < m_kk; k++) { + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } + m_tlast = tnow; + } +} + +doublereal MaskellSolidSolnPhase::s() const +{ + return 1 + std::exp(h_mixing / (GasConstant * temperature())); +} + +doublereal MaskellSolidSolnPhase::fm(const doublereal r) const +{ + const doublereal sval = s(); + return (1 - std::sqrt(1 - 4*r*(1-r)/sval)) / (2*r); +} + +doublereal MaskellSolidSolnPhase::p(const doublereal r) const +{ + const doublereal sval = s(); + return (1 - 2*r) / std::sqrt(sval*sval - 4 * sval * r + 4 * sval * r * r); +} + +} // end namespace Cantera diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 0c86e8c99..dca149f69 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -15,6 +15,7 @@ #include "VPSSMgrFactory.h" #include "cantera/thermo/IdealSolidSolnPhase.h" +#include "cantera/thermo/MaskellSolidSolnPhase.h" #include "cantera/thermo/MargulesVPSSTP.h" #include "cantera/thermo/RedlichKisterVPSSTP.h" #include "cantera/thermo/IonsFromNeutralVPSSTP.h" @@ -63,7 +64,7 @@ ThermoFactory* ThermoFactory::s_factory = 0; mutex_t ThermoFactory::thermo_mutex; //! Define the number of %ThermoPhase types for use in this factory routine -static int ntypes = 26; +static int ntypes = 27; //! Define the string name of the %ThermoPhase types that are handled by this factory routine static string _types[] = {"IdealGas", "Incompressible", @@ -74,7 +75,7 @@ static string _types[] = {"IdealGas", "Incompressible", "MineralEQ3", "MetalSHEelectrons", "Margules", "PhaseCombo_Interaction", "IonsFromNeutralMolecule", "FixedChemPot", "MolarityIonicVPSSTP", "MixedSolventElectrolyte", "Redlich-Kister", "RedlichKwong", - "RedlichKwongMFTP" + "RedlichKwongMFTP", "MaskellSolidSolnPhase" }; //! Define the integer id of the %ThermoPhase types that are handled by this factory routine @@ -86,7 +87,7 @@ static int _itypes[] = {cIdealGas, cIncompressible, cMineralEQ3, cMetalSHEelectrons, cMargulesVPSSTP, cPhaseCombo_Interaction, cIonsFromNeutral, cFixedChemPot, cMolarityIonicVPSSTP, cMixedSolventElectrolyte, cRedlichKisterVPSSTP, - cRedlichKwongMFTP, cRedlichKwongMFTP + cRedlichKwongMFTP, cRedlichKwongMFTP, cMaskellSolidSolnPhase }; ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model) @@ -203,6 +204,10 @@ ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model) th = new IdealSolnGasVPSS; break; + case cMaskellSolidSolnPhase: + th = new MaskellSolidSolnPhase; + break; + default: throw UnknownThermoPhaseModel("ThermoFactory::newThermoPhase", model); diff --git a/test/data/MaskellSolidSolnPhase_nohmix.xml b/test/data/MaskellSolidSolnPhase_nohmix.xml new file mode 100644 index 000000000..8516b09cf --- /dev/null +++ b/test/data/MaskellSolidSolnPhase_nohmix.xml @@ -0,0 +1,57 @@ + + + + + + + + H He + + + H(s) He(s) + + + 1.0 + + + + 298.15 + 1.0 + + H(s):0.90 He(s):0.10 + + + + + + + + H:1 He:2 + + + + 1. + + + + + 0.005 + + + + + H:0 He:1 + + + + 6.94544000E+01, + + + + + 0.005 + + + + + diff --git a/test/data/MaskellSolidSolnPhase_valid.xml b/test/data/MaskellSolidSolnPhase_valid.xml new file mode 100644 index 000000000..3848a51f6 --- /dev/null +++ b/test/data/MaskellSolidSolnPhase_valid.xml @@ -0,0 +1,58 @@ + + + + + + + + H He + + + H(s) He(s) + + + -1000. + 1.0 + + + + 298.15 + 1.0 + + H(s):0.90 He(s):0.10 + + + + + + + + H:1 He:2 + + + + 1. + + + + + 0.005 + + + + + H:0 He:1 + + + + 6.94544000E+01, + + + + + 0.005 + + + + + diff --git a/test/thermo/MaskellSolidSolnPhase_Test.cpp b/test/thermo/MaskellSolidSolnPhase_Test.cpp new file mode 100644 index 000000000..206ab20c0 --- /dev/null +++ b/test/thermo/MaskellSolidSolnPhase_Test.cpp @@ -0,0 +1,106 @@ +#include "gtest/gtest.h" +#include "cantera/thermo/MaskellSolidSolnPhase.h" +#include "cantera/thermo/SimpleThermo.h" +#include "cantera/thermo/ThermoFactory.h" +#include + +namespace Cantera +{ + +class MaskellSolidSolnPhase_Test : public testing::Test +{ +protected: + ThermoPhase *test_phase; +public: + MaskellSolidSolnPhase_Test() : test_phase(NULL) {} + + ~MaskellSolidSolnPhase_Test() { delete test_phase; } + + void initializeTestPhaseWithSimpleThermo() + { + test_phase = new MaskellSolidSolnPhase(); + test_phase->addElement("A", 1.); + test_phase->addElement("B", 2.); + std::vector comp(2); + comp[0] = 1.; + comp[1] = 0.; + test_phase->addSpecies("A", &comp[0], 0., 1.); + comp[0] = 0.; + comp[1] = 1.; + test_phase->addSpecies("B", &comp[0], 0., 1.); + + // Setup simple thermo so that the standard state enthalpy and + // gibbs free energies are always 0 so that we can just test the + // additional contribution from the Maskell model + SimpleThermo * spec_thermo = new SimpleThermo(); + std::vector coeffs(4); + coeffs[0] = 1; + coeffs[1] = 0; + coeffs[2] = 0; + coeffs[3] = 0; + spec_thermo->install("A", 0, 0, &coeffs[0], 0., 1000., 1.); + coeffs[1] = 1000; + spec_thermo->install("B", 1, 0, &coeffs[0], 0., 1000., 1.); + test_phase->setSpeciesThermo(spec_thermo); + + test_phase->setState_TP(298., 1.); + set_r(0.5); + } + + void initializeTestPhaseWithXML(const std::string & filename) + { + test_phase = newPhase(filename.c_str(), ""); + } + + void set_r(const double r) { + std::vector moleFracs(2); + moleFracs[0] = r; + moleFracs[1] = 1-r; + test_phase->setMoleFractions(&moleFracs[0]); + } + + void check_chemPotentials(const double expected_result[9]) + { + std::vector chemPotentials(2); + for(int i=0; i < 9; ++i) + { + const double r = 0.1 * (i+1); + set_r(r); + test_phase->getChemPotentials(&chemPotentials[0]); + EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6); + EXPECT_NEAR(1000.-expected_result[i], chemPotentials[1], 1.e-6); + } + } +}; + +TEST_F(MaskellSolidSolnPhase_Test, construct_from_xml) +{ + const std::string invalid_file("../data/MaskellSolidSolnPhase_nohmix.xml"); + EXPECT_THROW(initializeTestPhaseWithXML(invalid_file), CanteraError); + delete test_phase; + + const std::string valid_file("../data/MaskellSolidSolnPhase_valid.xml"); + initializeTestPhaseWithXML(valid_file); + MaskellSolidSolnPhase * maskell_phase = dynamic_cast(test_phase); +} + +TEST_F(MaskellSolidSolnPhase_Test, chem_potentials) +{ + initializeTestPhaseWithSimpleThermo(); + + MaskellSolidSolnPhase * maskell_phase = dynamic_cast(test_phase); + + maskell_phase->set_h_mix(0.); + const double expected_result_0[9] = {1.2338461168724738e7, 8.011774549216799e6, 4.990989640314685e6, 2.415973128783114e6, 0., -2.415973128783114e6, -4.99098964031469e6, -8.0117745492168e6, -1.2338461168724738e7}; + check_chemPotentials(expected_result_0); + + maskell_phase->set_h_mix(5000.); + const double expected_result_5000[9] = { 1.233625377465302e7, 8.00995666545047e6, 4.989677478024063e6, 2.41528026460977e6, 0., -2.415280264609771e6, -4.989677478024068e6, -8.00995666545047e6, -1.233625377465302e7 }; + check_chemPotentials(expected_result_5000); + + maskell_phase->set_h_mix(-5000.); + const double expected_result_minus_5000[9] = { 1.2340671035887627e7, 8.013594700219031e6, 4.992303607179179e6, 2.4166670154679064e6, 0., -2.4166670154679064e6, -4.9923036071791835e6, -8.013594700219034e6, -1.2340671035887627e7}; + check_chemPotentials(expected_result_minus_5000); +} + +};