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.
This commit is contained in:
Victor Brunini
parent
a313873b03
commit
93fcf181b1
7 changed files with 833 additions and 3 deletions
318
include/cantera/thermo/MaskellSolidSolnPhase.h
Normal file
318
include/cantera/thermo/MaskellSolidSolnPhase.h
Normal file
|
|
@ -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 <I>T</I> and <I>P</I> 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 <I>T</I> and <I>P</I> of the solution.
|
||||
/*!
|
||||
* These are the standard state chemical potentials \f$ \mu^0_k(T,P)
|
||||
* \f$. The values are evaluated at the current
|
||||
* temperature and pressure of the solution
|
||||
*
|
||||
* @param mu Output vector of chemical potentials.
|
||||
* Length: m_kk.
|
||||
*/
|
||||
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
|
||||
|
|
@ -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;
|
||||
|
||||
|
|
|
|||
284
src/thermo/MaskellSolidSolnPhase.cpp
Normal file
284
src/thermo/MaskellSolidSolnPhase.cpp
Normal file
|
|
@ -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 <cassert>
|
||||
#include <iostream>
|
||||
|
||||
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:
|
||||
* <thermo model="MaskellSolidSolution" />
|
||||
*/
|
||||
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
|
||||
|
|
@ -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);
|
||||
|
|
|
|||
57
test/data/MaskellSolidSolnPhase_nohmix.xml
Normal file
57
test/data/MaskellSolidSolnPhase_nohmix.xml
Normal file
|
|
@ -0,0 +1,57 @@
|
|||
|
||||
<?xml version="1.0"?>
|
||||
<ctml>
|
||||
<validate reactions="yes" species="yes"/>
|
||||
|
||||
<phase dim="3" id="Soln_Maskell9">
|
||||
<elementArray datasrc="elements.xml">
|
||||
H He
|
||||
</elementArray>
|
||||
<speciesArray datasrc="#species_Soln">
|
||||
H(s) He(s)
|
||||
</speciesArray>
|
||||
<thermo model="MaskellSolidSolnPhase">
|
||||
<density units="kg/m3"> 1.0 </density>
|
||||
</thermo>
|
||||
<kinetics model="none"/>
|
||||
<state>
|
||||
<temperature> 298.15 </temperature>
|
||||
<pressure units="atm"> 1.0 </pressure>
|
||||
<moleFractions>
|
||||
H(s):0.90 He(s):0.10
|
||||
</moleFractions>
|
||||
</state>
|
||||
</phase>
|
||||
|
||||
<!-- species definitions -->
|
||||
<speciesData id="species_Soln">
|
||||
<species name="H(s)">
|
||||
<atomArray> H:1 He:2 </atomArray>
|
||||
<thermo>
|
||||
<Shomate Pref="1 bar" Tmax="800." Tmin="250.0">
|
||||
<floatArray size="1">
|
||||
1.
|
||||
</floatArray>
|
||||
</Shomate>
|
||||
</thermo>
|
||||
<standardState model="constant_incompressible">
|
||||
<molarVolume units="m3/kmol"> 0.005 </molarVolume>
|
||||
</standardState>
|
||||
</species>
|
||||
|
||||
<species name="He(s)">
|
||||
<atomArray> H:0 He:1 </atomArray>
|
||||
<thermo>
|
||||
<Shomate Pref="1 bar" Tmax="800." Tmin="250.0">
|
||||
<floatArray size="1">
|
||||
6.94544000E+01,
|
||||
</floatArray>
|
||||
</Shomate>
|
||||
</thermo>
|
||||
<standardState model="constant_incompressible">
|
||||
<molarVolume units="m3/kmol"> 0.005 </molarVolume>
|
||||
</standardState>
|
||||
</species>
|
||||
</speciesData>
|
||||
|
||||
</ctml>
|
||||
58
test/data/MaskellSolidSolnPhase_valid.xml
Normal file
58
test/data/MaskellSolidSolnPhase_valid.xml
Normal file
|
|
@ -0,0 +1,58 @@
|
|||
|
||||
<?xml version="1.0"?>
|
||||
<ctml>
|
||||
<validate reactions="yes" species="yes"/>
|
||||
|
||||
<phase dim="3" id="Soln_Maskell9">
|
||||
<elementArray datasrc="elements.xml">
|
||||
H He
|
||||
</elementArray>
|
||||
<speciesArray datasrc="#species_Soln">
|
||||
H(s) He(s)
|
||||
</speciesArray>
|
||||
<thermo model="MaskellSolidSolnPhase">
|
||||
<h_mix>-1000.</h_mix>
|
||||
<density units="kg/m3"> 1.0 </density>
|
||||
</thermo>
|
||||
<kinetics model="none"/>
|
||||
<state>
|
||||
<temperature> 298.15 </temperature>
|
||||
<pressure units="atm"> 1.0 </pressure>
|
||||
<moleFractions>
|
||||
H(s):0.90 He(s):0.10
|
||||
</moleFractions>
|
||||
</state>
|
||||
</phase>
|
||||
|
||||
<!-- species definitions -->
|
||||
<speciesData id="species_Soln">
|
||||
<species name="H(s)">
|
||||
<atomArray> H:1 He:2 </atomArray>
|
||||
<thermo>
|
||||
<Shomate Pref="1 bar" Tmax="800." Tmin="250.0">
|
||||
<floatArray size="1">
|
||||
1.
|
||||
</floatArray>
|
||||
</Shomate>
|
||||
</thermo>
|
||||
<standardState model="constant_incompressible">
|
||||
<molarVolume units="m3/kmol"> 0.005 </molarVolume>
|
||||
</standardState>
|
||||
</species>
|
||||
|
||||
<species name="He(s)">
|
||||
<atomArray> H:0 He:1 </atomArray>
|
||||
<thermo>
|
||||
<Shomate Pref="1 bar" Tmax="800." Tmin="250.0">
|
||||
<floatArray size="1">
|
||||
6.94544000E+01,
|
||||
</floatArray>
|
||||
</Shomate>
|
||||
</thermo>
|
||||
<standardState model="constant_incompressible">
|
||||
<molarVolume units="m3/kmol"> 0.005 </molarVolume>
|
||||
</standardState>
|
||||
</species>
|
||||
</speciesData>
|
||||
|
||||
</ctml>
|
||||
106
test/thermo/MaskellSolidSolnPhase_Test.cpp
Normal file
106
test/thermo/MaskellSolidSolnPhase_Test.cpp
Normal file
|
|
@ -0,0 +1,106 @@
|
|||
#include "gtest/gtest.h"
|
||||
#include "cantera/thermo/MaskellSolidSolnPhase.h"
|
||||
#include "cantera/thermo/SimpleThermo.h"
|
||||
#include "cantera/thermo/ThermoFactory.h"
|
||||
#include <iostream>
|
||||
|
||||
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<double> 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<double> 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<double> moleFracs(2);
|
||||
moleFracs[0] = r;
|
||||
moleFracs[1] = 1-r;
|
||||
test_phase->setMoleFractions(&moleFracs[0]);
|
||||
}
|
||||
|
||||
void check_chemPotentials(const double expected_result[9])
|
||||
{
|
||||
std::vector<double> 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<MaskellSolidSolnPhase *>(test_phase);
|
||||
}
|
||||
|
||||
TEST_F(MaskellSolidSolnPhase_Test, chem_potentials)
|
||||
{
|
||||
initializeTestPhaseWithSimpleThermo();
|
||||
|
||||
MaskellSolidSolnPhase * maskell_phase = dynamic_cast<MaskellSolidSolnPhase *>(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);
|
||||
}
|
||||
|
||||
};
|
||||
Loading…
Add table
Reference in a new issue