From 08a855359a4e899a6cbcd6602c387aee2b39747a Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Tue, 30 May 2006 16:15:48 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/src/LatticePhase.cpp | 4 +- Cantera/src/LatticePhase.h | 15 +- Cantera/src/LatticeSolidPhase.cpp | 227 +++++++++++++++++++----------- Cantera/src/LatticeSolidPhase.h | 90 ++---------- Cantera/src/MultiPhaseEquil.cpp | 9 +- Cantera/src/ThermoFactory.cpp | 11 +- Cantera/src/mix_defs.h | 1 + 7 files changed, 174 insertions(+), 183 deletions(-) diff --git a/Cantera/src/LatticePhase.cpp b/Cantera/src/LatticePhase.cpp index bc5e528bd..5704b1a18 100644 --- a/Cantera/src/LatticePhase.cpp +++ b/Cantera/src/LatticePhase.cpp @@ -117,8 +117,8 @@ namespace Cantera { void LatticePhase::setParametersFromXML(const XML_Node& eosdata) { eosdata._require("model","Lattice"); - m_molar_density = getFloat(i, "site_density", "-"); - m_vacancy = getString(i, "vacancy_species"); + m_molar_density = getFloat(eosdata, "site_density", "-"); + m_vacancy = getString(eosdata, "vacancy_species"); } } diff --git a/Cantera/src/LatticePhase.h b/Cantera/src/LatticePhase.h index 04cfe779e..add856adf 100644 --- a/Cantera/src/LatticePhase.h +++ b/Cantera/src/LatticePhase.h @@ -94,10 +94,6 @@ namespace Cantera { // new methods defined here - double latticeMoleFraction(int k) { - return moleFraction(k)*molarDensity()/m_sitedens[m_lattice[k]]; - } - const array_fp& enthalpy_RT() const { _updateThermo(); return m_h0_RT; @@ -136,18 +132,15 @@ namespace Cantera { int m_mm; doublereal m_tmin, m_tmax, m_p0; - mutable doublereal m_tlast; mutable array_fp m_h0_RT; mutable array_fp m_cp0_R; mutable array_fp m_g0_RT; mutable array_fp m_s0_R; - doublereal m_press; - vector m_vac; - vector_fp m_sitedens; - doublereal m_molar_density; - vector m_lattice; - vector m_sp; + doublereal m_press; + string m_vacancy; + doublereal m_molar_density; + private: void _updateThermo() const; diff --git a/Cantera/src/LatticeSolidPhase.cpp b/Cantera/src/LatticeSolidPhase.cpp index 61bd55146..6c8fe3a39 100644 --- a/Cantera/src/LatticeSolidPhase.cpp +++ b/Cantera/src/LatticeSolidPhase.cpp @@ -13,40 +13,76 @@ #include "ct_defs.h" #include "mix_defs.h" #include "LatticeSolidPhase.h" +#include "LatticePhase.h" #include "SpeciesThermo.h" +#include "importCTML.h" namespace Cantera { doublereal LatticeSolidPhase:: enthalpy_mole() const { - doublereal p0 = m_spthermo->refPressure(); - return GasConstant * temperature() * - mean_X(&enthalpy_RT()[0]) - + (pressure() - p0)/molarDensity(); + _updateThermo(); + doublereal ndens, sum = 0.0; + int n; + for (n = 0; n < m_nlattice; n++) { + ndens = m_lattice[n]->molarDensity(); + sum += ndens * m_lattice[n]->enthalpy_mole(); + } + return sum/molarDensity(); } doublereal LatticeSolidPhase::intEnergy_mole() const { - doublereal p0 = m_spthermo->refPressure(); - return GasConstant * temperature() * - mean_X(&enthalpy_RT()[0]) - - p0/molarDensity(); + _updateThermo(); + doublereal ndens, sum = 0.0; + int n; + for (n = 0; n < m_nlattice; n++) { + ndens = m_lattice[n]->molarDensity(); + sum += ndens * m_lattice[n]->intEnergy_mole(); + } + return sum/molarDensity(); } doublereal LatticeSolidPhase::entropy_mole() const { - return GasConstant * (mean_X(&entropy_R()[0]) - - sum_xlogx()); + _updateThermo(); + doublereal ndens, sum = 0.0; + int n; + for (n = 0; n < m_nlattice; n++) { + ndens = m_lattice[n]->molarDensity(); + sum += ndens * m_lattice[n]->entropy_mole(); + } + return sum/molarDensity(); } doublereal LatticeSolidPhase::gibbs_mole() const { - return enthalpy_mole() - temperature() * entropy_mole(); + _updateThermo(); + doublereal ndens, sum = 0.0; + int n; + for (n = 0; n < m_nlattice; n++) { + ndens = m_lattice[n]->molarDensity(); + sum += ndens * m_lattice[n]->gibbs_mole(); + } + return sum/molarDensity(); } doublereal LatticeSolidPhase::cp_mole() const { - return GasConstant * mean_X(&cp_R()[0]); + _updateThermo(); + doublereal ndens, sum = 0.0; + int n; + for (n = 0; n < m_nlattice; n++) { + ndens = m_lattice[n]->molarDensity(); + sum += ndens * m_lattice[n]->cp_mole(); + } + return sum/molarDensity(); } void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const { - getMoleFractions(c); + _updateThermo(); + int n; + int strt = 0; + for (n = 0; n < m_nlattice; n++) { + m_lattice[n]->getMoleFractions(c+strt); + strt += m_lattice[n]->nSpecies(); + } } void LatticeSolidPhase::getActivityCoefficients(doublereal* ac) const { @@ -64,111 +100,130 @@ namespace Cantera { } void LatticeSolidPhase::getChemPotentials(doublereal* mu) const { - doublereal vdp = (pressure() - m_spthermo->refPressure())/ - molarDensity(); - doublereal xx; - doublereal rt = temperature() * GasConstant; - const array_fp& g_RT = gibbs_RT(); - for (int k = 0; k < m_kk; k++) { - xx = fmaxx(SmallNumber, moleFraction(k)); - mu[k] = rt*(g_RT[k] + log(xx)) + vdp; + _updateThermo(); + int n; + int strt = 0; + double dratio; + for (n = 0; n < m_nlattice; n++) { + dratio = m_lattice[n]->molarDensity()/molarDensity(); + m_lattice[n]->getChemPotentials(mu+strt); + scale(mu + strt, mu + strt + m_lattice[n]->nSpecies(), mu + strt, dratio); + strt += m_lattice[n]->nSpecies(); } } void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const { - getPureGibbs(mu0); + _updateThermo(); + int n; + int strt = 0; + double dratio; + for (n = 0; n < m_nlattice; n++) { + dratio = m_lattice[n]->molarDensity()/molarDensity(); + m_lattice[n]->getStandardChemPotentials(mu0+strt); + scale(mu0 + strt, mu0 + strt + m_lattice[n]->nSpecies(), mu0 + strt, dratio); + strt += m_lattice[n]->nSpecies(); + } } void LatticeSolidPhase::initThermo() { m_kk = nSpecies(); m_mm = nElements(); - doublereal tmin = m_spthermo->minTemp(); - doublereal tmax = m_spthermo->maxTemp(); - if (tmin > 0.0) m_tmin = tmin; - if (tmax > 0.0) m_tmax = tmax; - m_p0 = refPressure(); - - int leng = m_kk; - m_h0_RT.resize(leng); - m_g0_RT.resize(leng); - m_cp0_R.resize(leng); - m_s0_R.resize(leng); - setMolarDensity(m_molar_density); - - const vector& spnames = speciesNames(); - int n, k, kl, namesize; - int nl = m_sitedens.size(); - string s; - m_lattice.resize(m_kk,-1); - vector_fp conc(m_kk, 0.0); - - compositionMap xx; - for (n = 0; n < nl; n++) { - for (k = 0; k < m_kk; k++) { - xx[speciesName(k)] = -1.0; - } - parseCompString(m_sp[n], xx); - for (k = 0; k < m_kk; k++) { - if (xx[speciesName(k)] != -1.0) { - conc[k] = m_sitedens[n]*xx[speciesName(k)]; - m_lattice[k] = n; - } + m_x.resize(m_kk); + int n, nsp, k, loc = 0; + doublereal ndens; + m_molar_density = 0.0; + for (n = 0; n < m_nlattice; n++) { + nsp = m_lattice[n]->nSpecies(); + ndens = m_lattice[n]->molarDensity(); + for (k = 0; k < nsp; k++) { + m_x[loc] = ndens * m_lattice[n]->moleFraction(k); + loc++; } + m_molar_density += ndens; + } + setMoleFractions(DATA_PTR(m_x)); - } - for (k = 0; k < m_kk; k++) { - if (m_lattice[k] == -1) { - throw CanteraError("LatticeSolidPhase::" - "setParametersFromXML","Species "+speciesName(k) - +" not a member of any lattice."); - } - } - setMoleFractions(DATA_PTR(conc)); +// const vector& spnames = speciesNames(); +// int n, k, kl, namesize; +// int nl = m_sitedens.size(); +// string s; +// m_lattice.resize(m_kk,-1); +// vector_fp conc(m_kk, 0.0); + +// compositionMap xx; +// for (n = 0; n < nl; n++) { +// for (k = 0; k < m_kk; k++) { +// xx[speciesName(k)] = -1.0; +// } +// parseCompString(m_sp[n], xx); +// for (k = 0; k < m_kk; k++) { +// if (xx[speciesName(k)] != -1.0) { +// conc[k] = m_sitedens[n]*xx[speciesName(k)]; +// m_lattice[k] = n; +// } +// } + +// } +// for (k = 0; k < m_kk; k++) { +// if (m_lattice[k] == -1) { +// throw CanteraError("LatticeSolidPhase::" +// "setParametersFromXML","Species "+speciesName(k) +// +" not a member of any lattice."); +// } +// } +// setMoleFractions(DATA_PTR(conc)); } void LatticeSolidPhase::_updateThermo() const { doublereal tnow = temperature(); - if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) { - throw CanteraError("_updateThermo","molar density changed from " - +fp2str(m_molar_density)+" to "+fp2str(molarDensity())); - } + // if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) { + // throw CanteraError("_updateThermo","molar density changed from " + // +fp2str(m_molar_density)+" to "+fp2str(molarDensity())); + //} if (m_tlast != tnow) { - m_spthermo->update(tnow, &m_cp0_R[0], &m_h0_RT[0], - &m_s0_R[0]); - m_tlast = tnow; - int k; - for (k = 0; k < m_kk; k++) { - m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + int n; + getMoleFractions(DATA_PTR(m_x)); + int strt = 0; + for (n = 0; n < m_nlattice; n++) { + m_lattice[n]->setTemperature(tnow); + m_lattice[n]->setMoleFractions(DATA_PTR(m_x) + strt); + m_lattice[n]->setPressure(m_press); + strt += m_lattice[n]->nSpecies(); } m_tlast = tnow; } } + void LatticeSolidPhase::setLatticeMoleFractions(int nn, + string x) { + m_lattice[nn]->setMoleFractionsByName(x); + int n, k, loc=0, nsp; + doublereal ndens; + for (n = 0; n < m_nlattice; n++) { + nsp = m_lattice[n]->nSpecies(); + ndens = m_lattice[n]->molarDensity(); + for (k = 0; k < nsp; k++) { + m_x[loc] = ndens * m_lattice[n]->moleFraction(k); + loc++; + } + } + setMoleFractions(DATA_PTR(m_x)); + } + void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata) { eosdata._require("model","LatticeSolid"); XML_Node& la = eosdata.child("LatticeArray"); vector lattices; - la.getChildren("Lattice",lattices); + la.getChildren("phase",lattices); int n; int nl = lattices.size(); - doublereal site_density; - string vacancy; - doublereal sum = 0.0; - string s; + m_nlattice = nl; for (n = 0; n < nl; n++) { XML_Node& i = *lattices[n]; - site_density = getFloat(i, "site_density", "-"); - vacancy = getString(i, "vacancy_species"); - s = getString(i, "species"); - m_sp.push_back(s); - m_vac.push_back(vacancy); - m_sitedens.push_back(site_density); - sum += site_density; + m_lattice.push_back((LatticePhase*)newPhase(i)); } - m_molar_density = sum; } - } diff --git a/Cantera/src/LatticeSolidPhase.h b/Cantera/src/LatticeSolidPhase.h index bac65af91..760fc9dad 100644 --- a/Cantera/src/LatticeSolidPhase.h +++ b/Cantera/src/LatticeSolidPhase.h @@ -22,10 +22,8 @@ namespace Cantera { - /** - * Overloads the virtual methods of class Thermo to implement the - * incompressible equation of state. - */ + class LatticePhase; + class LatticeSolidPhase : public ThermoPhase { public: @@ -68,88 +66,24 @@ namespace Cantera { virtual doublereal standardConcentration(int k=0) const; virtual doublereal logStandardConc(int k=0) const; - virtual void getPureGibbs(doublereal* gpure) const { - const array_fp& gibbsrt = gibbs_RT(); - scale(gibbsrt.begin(), gibbsrt.end(), gpure, _RT()); - } - - void getEnthalpy_RT(doublereal* hrt) const { - const array_fp& _h = enthalpy_RT(); - copy(_h.begin(), _h.end(), hrt); - } - - void getEntropy_R(doublereal* sr) const { - const array_fp& _s = entropy_R(); - copy(_s.begin(), _s.end(), sr); - } - - virtual void getGibbs_RT(doublereal* grt) const { - const array_fp& gibbsrt = gibbs_RT(); - copy(gibbsrt.begin(), gibbsrt.end(), grt); - } - - void getCp_R(doublereal* cpr) const { - const array_fp& _cpr = cp_R(); - copy(_cpr.begin(), _cpr.end(), cpr); - } - - - // new methods defined here - - double latticeMoleFraction(int k) { - return moleFraction(k)*molarDensity()/m_sitedens[m_lattice[k]]; - } - - const array_fp& enthalpy_RT() const { - _updateThermo(); - return m_h0_RT; - } - - const array_fp& gibbs_RT() const { - _updateThermo(); - return m_g0_RT; - } - - const array_fp& entropy_R() const { - _updateThermo(); - return m_s0_R; - } - - const array_fp& cp_R() const { - _updateThermo(); - return m_cp0_R; - } - virtual void initThermo(); - // set the site density of sublattice n - virtual void setParameters(int n, doublereal* c) {} - - virtual void getParameters(int &n, doublereal * const c) { - double d = molarDensity(); - c[0] = d; - n = 1; - } - virtual void setParametersFromXML(const XML_Node& eosdata); + void LatticeSolidPhase::setLatticeMoleFractions(int n, + string x); protected: - + int m_mm; - doublereal m_tmin, m_tmax, m_p0; + int m_kk; + mutable doublereal m_tlast; + doublereal m_press; + doublereal m_molar_density; + int m_nlattice; + vector m_lattice; + mutable vector_fp m_x; - mutable doublereal m_tlast; - mutable array_fp m_h0_RT; - mutable array_fp m_cp0_R; - mutable array_fp m_g0_RT; - mutable array_fp m_s0_R; - doublereal m_press; - vector m_vac; - vector_fp m_sitedens; - doublereal m_molar_density; - vector m_lattice; - vector m_sp; private: void _updateThermo() const; diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp index 7c9e0685b..7927fed56 100644 --- a/Cantera/src/MultiPhaseEquil.cpp +++ b/Cantera/src/MultiPhaseEquil.cpp @@ -800,13 +800,16 @@ namespace Cantera { // don't require formation reactions for solution species // present in trace amounts to be equilibrated - if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) + if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) { err = 0.0; + } // for stoichiometric phase species, no error if not present and // delta G for the formation reaction is positive - else if (isStoichPhase(ik) && moles(ik) <= 0.0 && - m_deltaG_RT[j] >= 0.0) err = 0.0; + if (isStoichPhase(ik) && moles(ik) <= 0.0 && + m_deltaG_RT[j] >= 0.0) { + err = 0.0; + } else { err = fabs(m_deltaG_RT[j]); } diff --git a/Cantera/src/ThermoFactory.cpp b/Cantera/src/ThermoFactory.cpp index 8c32b5150..1fd2bc558 100644 --- a/Cantera/src/ThermoFactory.cpp +++ b/Cantera/src/ThermoFactory.cpp @@ -40,20 +40,21 @@ #ifdef WITH_LATTICE_SOLID #include "LatticeSolidPhase.h" +#include "LatticePhase.h" #endif namespace Cantera { ThermoFactory* ThermoFactory::s_factory = 0; - static int ntypes = 8; + static int ntypes = 9; static string _types[] = {"IdealGas", "Incompressible", "Surface", "Edge", "Metal", "StoichSubstance", - "PureFluid", "LatticeSolid"}; + "PureFluid", "LatticeSolid", "Lattice"}; static int _itypes[] = {cIdealGas, cIncompressible, cSurf, cEdge, cMetal, cStoichSubstance, - cPureFluid, cLatticeSolid}; + cPureFluid, cLatticeSolid, cLattice}; /** * This method returns a new instance of a subclass of ThermoPhase @@ -102,6 +103,10 @@ namespace Cantera { case cLatticeSolid: th = new LatticeSolidPhase; break; + + case cLattice: + th = new LatticePhase; + break; #endif #ifdef WITH_PURE_FLUIDS diff --git a/Cantera/src/mix_defs.h b/Cantera/src/mix_defs.h index df6d4f0d6..cbcb42718 100755 --- a/Cantera/src/mix_defs.h +++ b/Cantera/src/mix_defs.h @@ -56,6 +56,7 @@ namespace Cantera { const int cInterfaceKinetics = 4; const int cLineKinetics = 5; const int cEdgeKinetics = 6; + const int cSolidKinetics = 7; } #endif