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

*** empty log message ***

Browse source
This commit is contained in:
Dave Goodwin 2006-05-30 16:15:48 +00:00
parent c3d097140d
commit 08a855359a
7 changed files with 174 additions and 183 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/src/LatticePhase.cpp
View file

@ -117,8 +117,8 @@ namespace Cantera {
void LatticePhase::setParametersFromXML(const XML_Node& eosdata) { void LatticePhase::setParametersFromXML(const XML_Node& eosdata) {
eosdata._require("model","Lattice"); eosdata._require("model","Lattice");
m_molar_density = getFloat(i, "site_density", "-"); m_molar_density = getFloat(eosdata, "site_density", "-");
m_vacancy = getString(i, "vacancy_species"); m_vacancy = getString(eosdata, "vacancy_species");
} }
} }

15
Cantera/src/LatticePhase.h
View file

@ -94,10 +94,6 @@ namespace Cantera {
// new methods defined here // new methods defined here
double latticeMoleFraction(int k) {
return moleFraction(k)*molarDensity()/m_sitedens[m_lattice[k]];
}
const array_fp& enthalpy_RT() const { const array_fp& enthalpy_RT() const {
_updateThermo(); _updateThermo();
return m_h0_RT; return m_h0_RT;
@ -136,18 +132,15 @@ namespace Cantera {
int m_mm; int m_mm;
doublereal m_tmin, m_tmax, m_p0; doublereal m_tmin, m_tmax, m_p0;
mutable doublereal m_tlast; mutable doublereal m_tlast;
mutable array_fp m_h0_RT; mutable array_fp m_h0_RT;
mutable array_fp m_cp0_R; mutable array_fp m_cp0_R;
mutable array_fp m_g0_RT; mutable array_fp m_g0_RT;
mutable array_fp m_s0_R; mutable array_fp m_s0_R;
doublereal m_press; doublereal m_press;
vector<string> m_vac; string m_vacancy;
vector_fp m_sitedens; doublereal m_molar_density;
doublereal m_molar_density;
vector<int> m_lattice;
vector<string> m_sp;
private: private:
void _updateThermo() const; void _updateThermo() const;

227
Cantera/src/LatticeSolidPhase.cpp
View file

@ -13,40 +13,76 @@
#include "ct_defs.h" #include "ct_defs.h"
#include "mix_defs.h" #include "mix_defs.h"
#include "LatticeSolidPhase.h" #include "LatticeSolidPhase.h"
#include "LatticePhase.h"
#include "SpeciesThermo.h" #include "SpeciesThermo.h"
#include "importCTML.h"
namespace Cantera { namespace Cantera {
doublereal LatticeSolidPhase:: doublereal LatticeSolidPhase::
enthalpy_mole() const { enthalpy_mole() const {
doublereal p0 = m_spthermo->refPressure(); _updateThermo();
return GasConstant * temperature() * doublereal ndens, sum = 0.0;
mean_X(&enthalpy_RT()[0]) int n;
+ (pressure() - p0)/molarDensity(); 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 LatticeSolidPhase::intEnergy_mole() const {
doublereal p0 = m_spthermo->refPressure(); _updateThermo();
return GasConstant * temperature() * doublereal ndens, sum = 0.0;
mean_X(&enthalpy_RT()[0]) int n;
- p0/molarDensity(); 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 { doublereal LatticeSolidPhase::entropy_mole() const {
return GasConstant * (mean_X(&entropy_R()[0]) - _updateThermo();
sum_xlogx()); 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 { 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 { 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 { 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 { void LatticeSolidPhase::getActivityCoefficients(doublereal* ac) const {
@ -64,111 +100,130 @@ namespace Cantera {
} }
void LatticeSolidPhase::getChemPotentials(doublereal* mu) const { void LatticeSolidPhase::getChemPotentials(doublereal* mu) const {
doublereal vdp = (pressure() - m_spthermo->refPressure())/ _updateThermo();
molarDensity(); int n;
doublereal xx; int strt = 0;
doublereal rt = temperature() * GasConstant; double dratio;
const array_fp& g_RT = gibbs_RT(); for (n = 0; n < m_nlattice; n++) {
for (int k = 0; k < m_kk; k++) { dratio = m_lattice[n]->molarDensity()/molarDensity();
xx = fmaxx(SmallNumber, moleFraction(k)); m_lattice[n]->getChemPotentials(mu+strt);
mu[k] = rt*(g_RT[k] + log(xx)) + vdp; scale(mu + strt, mu + strt + m_lattice[n]->nSpecies(), mu + strt, dratio);
strt += m_lattice[n]->nSpecies();
} }
} }
void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const { 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() { void LatticeSolidPhase::initThermo() {
m_kk = nSpecies(); m_kk = nSpecies();
m_mm = nElements(); m_mm = nElements();
doublereal tmin = m_spthermo->minTemp(); m_x.resize(m_kk);
doublereal tmax = m_spthermo->maxTemp(); int n, nsp, k, loc = 0;
if (tmin > 0.0) m_tmin = tmin; doublereal ndens;
if (tmax > 0.0) m_tmax = tmax; m_molar_density = 0.0;
m_p0 = refPressure(); for (n = 0; n < m_nlattice; n++) {
nsp = m_lattice[n]->nSpecies();
int leng = m_kk; ndens = m_lattice[n]->molarDensity();
m_h0_RT.resize(leng); for (k = 0; k < nsp; k++) {
m_g0_RT.resize(leng); m_x[loc] = ndens * m_lattice[n]->moleFraction(k);
m_cp0_R.resize(leng); loc++;
m_s0_R.resize(leng);
setMolarDensity(m_molar_density);
const vector<string>& 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_molar_density += ndens;
}
setMoleFractions(DATA_PTR(m_x));
} // const vector<string>& spnames = speciesNames();
for (k = 0; k < m_kk; k++) { // int n, k, kl, namesize;
if (m_lattice[k] == -1) { // int nl = m_sitedens.size();
throw CanteraError("LatticeSolidPhase::" // string s;
"setParametersFromXML","Species "+speciesName(k) // m_lattice.resize(m_kk,-1);
+" not a member of any lattice."); // vector_fp conc(m_kk, 0.0);
}
} // compositionMap xx;
setMoleFractions(DATA_PTR(conc)); // 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 { void LatticeSolidPhase::_updateThermo() const {
doublereal tnow = temperature(); doublereal tnow = temperature();
if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) { // if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) {
throw CanteraError("_updateThermo","molar density changed from " // throw CanteraError("_updateThermo","molar density changed from "
+fp2str(m_molar_density)+" to "+fp2str(molarDensity())); // +fp2str(m_molar_density)+" to "+fp2str(molarDensity()));
} //}
if (m_tlast != tnow) { if (m_tlast != tnow) {
m_spthermo->update(tnow, &m_cp0_R[0], &m_h0_RT[0], int n;
&m_s0_R[0]); getMoleFractions(DATA_PTR(m_x));
m_tlast = tnow; int strt = 0;
int k; for (n = 0; n < m_nlattice; n++) {
for (k = 0; k < m_kk; k++) { m_lattice[n]->setTemperature(tnow);
m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; m_lattice[n]->setMoleFractions(DATA_PTR(m_x) + strt);
m_lattice[n]->setPressure(m_press);
strt += m_lattice[n]->nSpecies();
} }
m_tlast = tnow; 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) { void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata) {
eosdata._require("model","LatticeSolid"); eosdata._require("model","LatticeSolid");
XML_Node& la = eosdata.child("LatticeArray"); XML_Node& la = eosdata.child("LatticeArray");
vector<XML_Node*> lattices; vector<XML_Node*> lattices;
la.getChildren("Lattice",lattices); la.getChildren("phase",lattices);
int n; int n;
int nl = lattices.size(); int nl = lattices.size();
doublereal site_density; m_nlattice = nl;
string vacancy;
doublereal sum = 0.0;
string s;
for (n = 0; n < nl; n++) { for (n = 0; n < nl; n++) {
XML_Node& i = *lattices[n]; XML_Node& i = *lattices[n];
site_density = getFloat(i, "site_density", "-"); m_lattice.push_back((LatticePhase*)newPhase(i));
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_molar_density = sum;
} }
} }

90
Cantera/src/LatticeSolidPhase.h
View file

@ -22,10 +22,8 @@
namespace Cantera { namespace Cantera {
/** class LatticePhase;
* Overloads the virtual methods of class Thermo to implement the
* incompressible equation of state.
*/
class LatticeSolidPhase : public ThermoPhase { class LatticeSolidPhase : public ThermoPhase {
public: public:
@ -68,88 +66,24 @@ namespace Cantera {
virtual doublereal standardConcentration(int k=0) const; virtual doublereal standardConcentration(int k=0) const;
virtual doublereal logStandardConc(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(); 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); virtual void setParametersFromXML(const XML_Node& eosdata);
void LatticeSolidPhase::setLatticeMoleFractions(int n,
string x);
protected: protected:
int m_mm; 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<LatticePhase*> 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<string> m_vac;
vector_fp m_sitedens;
doublereal m_molar_density;
vector<int> m_lattice;
vector<string> m_sp;
private: private:
void _updateThermo() const; void _updateThermo() const;

9
Cantera/src/MultiPhaseEquil.cpp
View file

@ -800,13 +800,16 @@ namespace Cantera {
// don't require formation reactions for solution species // don't require formation reactions for solution species
// present in trace amounts to be equilibrated // present in trace amounts to be equilibrated
if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) if (!isStoichPhase(ik) && fabs(moles(ik)) <= SmallNumber) {
err = 0.0; err = 0.0;
}
// for stoichiometric phase species, no error if not present and // for stoichiometric phase species, no error if not present and
// delta G for the formation reaction is positive // delta G for the formation reaction is positive
else if (isStoichPhase(ik) && moles(ik) <= 0.0 && if (isStoichPhase(ik) && moles(ik) <= 0.0 &&
m_deltaG_RT[j] >= 0.0) err = 0.0; m_deltaG_RT[j] >= 0.0) {
err = 0.0;
}
else { else {
err = fabs(m_deltaG_RT[j]); err = fabs(m_deltaG_RT[j]);
} }

11
Cantera/src/ThermoFactory.cpp
View file

@ -40,20 +40,21 @@
#ifdef WITH_LATTICE_SOLID #ifdef WITH_LATTICE_SOLID
#include "LatticeSolidPhase.h" #include "LatticeSolidPhase.h"
#include "LatticePhase.h"
#endif #endif
namespace Cantera { namespace Cantera {
ThermoFactory* ThermoFactory::s_factory = 0; ThermoFactory* ThermoFactory::s_factory = 0;
static int ntypes = 8; static int ntypes = 9;
static string _types[] = {"IdealGas", "Incompressible", static string _types[] = {"IdealGas", "Incompressible",
"Surface", "Edge", "Metal", "StoichSubstance", "Surface", "Edge", "Metal", "StoichSubstance",
"PureFluid", "LatticeSolid"}; "PureFluid", "LatticeSolid", "Lattice"};
static int _itypes[] = {cIdealGas, cIncompressible, static int _itypes[] = {cIdealGas, cIncompressible,
cSurf, cEdge, cMetal, cStoichSubstance, cSurf, cEdge, cMetal, cStoichSubstance,
cPureFluid, cLatticeSolid}; cPureFluid, cLatticeSolid, cLattice};
/** /**
* This method returns a new instance of a subclass of ThermoPhase * This method returns a new instance of a subclass of ThermoPhase
@ -102,6 +103,10 @@ namespace Cantera {
case cLatticeSolid: case cLatticeSolid:
th = new LatticeSolidPhase; th = new LatticeSolidPhase;
break; break;
case cLattice:
th = new LatticePhase;
break;
#endif #endif
#ifdef WITH_PURE_FLUIDS #ifdef WITH_PURE_FLUIDS

1
Cantera/src/mix_defs.h
View file

@ -56,6 +56,7 @@ namespace Cantera {
const int cInterfaceKinetics = 4; const int cInterfaceKinetics = 4;
const int cLineKinetics = 5; const int cLineKinetics = 5;
const int cEdgeKinetics = 6; const int cEdgeKinetics = 6;
const int cSolidKinetics = 7;
} }
#endif #endif