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

Simplify treatment of units within the equilibrium solver

Browse source 
-> always assumes numbers are in kmol
  -> always assumes pressures are in pascal
  -> always assumes volumes are in m**3
This commit is contained in:
Harry Moffat 2008-04-28 16:08:35 +00:00
parent 72d3f7b7f7
commit 7b4e4ec1c4
16 changed files with 141 additions and 146 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
Cantera/src/equil/vcs_Gibbs.cpp
View file

@ -45,11 +45,7 @@ namespace VCSnonideal {
g += TPhInertMoles[iph] *
log(TPhInertMoles[iph] / tPhMoles[iph]);
if (Vphase->GasPhase) {
if (m_VCS_UnitsFormat == 3) {
g += TPhInertMoles[iph] * log(m_pressure/(Cantera::OneAtm));
} else {
g += TPhInertMoles[iph] * log(m_pressure);
}
g += TPhInertMoles[iph] * log(m_pressurePA/(1.01325E5));
}
}
}
@ -86,11 +82,7 @@ namespace VCSnonideal {
g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols);
vcs_VolPhase *Vphase = VPhaseList[iphase];
if (Vphase->GasPhase == iphase) {
if (m_VCS_UnitsFormat == 3) {
g += TPhInertMoles[iphase] * log(m_pressure/(Cantera::OneAtm));
} else {
g += TPhInertMoles[iphase] * log(m_pressure);
}
g += TPhInertMoles[iphase] * log(m_pressurePA/1.01325E5);
}
}

35
Cantera/src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -668,12 +668,7 @@ namespace VCSnonideal {
}
plogf("\n");
plogf("Temperature = %g Kelvin\n", m_vprob->T);
plogf("Pressure = %g ", m_vprob->Pres);
if (m_vprob->m_VCS_UnitsFormat == VCS_UNITS_MKS) {
plogf("Pa\n");
} else {
plogf("atm\n");
}
plogf("Pressure = %g Pa\n", m_vprob->PresPA);
plogf("\n");
plogf("----------------------------------------"
"---------------------\n");
@ -1099,7 +1094,7 @@ namespace VCSnonideal {
// We will work out the details later.
vprob->iest = -1;
vprob->T = mphase->temperature();
vprob->Pres = mphase->pressure();
vprob->PresPA = mphase->pressure();
vprob->Vol = mphase->volume();
vprob->Title = "MultiPhase Object";
@ -1225,7 +1220,7 @@ namespace VCSnonideal {
*/
vprob->addPhaseElements(VolPhase);
VolPhase->setState_TP(vprob->T, vprob->Pres);
VolPhase->setState_TP(vprob->T, vprob->PresPA);
vector<double> muPhase(tPhase->nSpecies(),0.0);
tPhase->getChemPotentials(&muPhase[0]);
double tMoles = 0.0;
@ -1426,14 +1421,14 @@ namespace VCSnonideal {
print_char('=', 80); plogf("\n");
plogf(" Phase IDs of species\n");
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_gMols\n");
plogf(" Initial_Estimated_kMols\n");
for (int i = 0; i < vprob->nspecies; i++) {
int iphase = vprob->PhaseID[i];
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
VolPhase->PhaseName.c_str());
plogf(" %-10.5g\n", 1.0E3 * vprob->w[i]);
plogf(" %-10.5g\n", vprob->w[i]);
}
/*
@ -1442,7 +1437,7 @@ namespace VCSnonideal {
plogf("\n"); print_char('-', 80); plogf("\n");
plogf(" Information about phases\n");
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
plogf(" TMolesInert Tmoles(gmol)\n");
plogf(" TMolesInert Tmoles(kmol)\n");
for (int iphase = 0; iphase < vprob->NPhase; iphase++) {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
@ -1450,8 +1445,8 @@ namespace VCSnonideal {
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->SingleSpecies,
VolPhase->GasPhase, sEOS.c_str(),
VolPhase->NVolSpecies, VolPhase->TMolesInert * 1.0E3);
plogf("%16e\n", VolPhase->TotalMoles() * 1.0E3);
VolPhase->NVolSpecies, VolPhase->TMolesInert );
plogf("%16e\n", VolPhase->TotalMoles());
}
plogf("\n"); print_char('=', 80); plogf("\n");
@ -1479,7 +1474,7 @@ namespace VCSnonideal {
// the call to the equilibrium solver.
vprob->iest = -1;
vprob->T = mphase->temperature();
vprob->Pres = mphase->pressure();
vprob->PresPA = mphase->pressure();
vprob->Vol = mphase->volume();
Cantera::ThermoPhase *tPhase = 0;
@ -1492,7 +1487,7 @@ namespace VCSnonideal {
*/
volPhase->setElectricPotential(tPhase->electricPotential());
volPhase->setState_TP(vprob->T, vprob->Pres);
volPhase->setState_TP(vprob->T, vprob->PresPA);
vector<double> muPhase(tPhase->nSpecies(),0.0);
tPhase->getChemPotentials(&muPhase[0]);
/*
@ -1544,14 +1539,14 @@ namespace VCSnonideal {
print_char('=', 80); plogf("\n\n");
plogf(" Phase IDs of species\n");
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_gMols\n");
plogf(" Initial_Estimated_kMols\n");
for (int i = 0; i < vprob->nspecies; i++) {
int iphase = vprob->PhaseID[i];
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
plogf("%16s %5d %16s", vprob->SpName[i].c_str(), iphase,
VolPhase->PhaseName.c_str());
plogf(" %-10.5g\n", 1.0E3 * vprob->w[i]);
plogf(" %-10.5g\n", vprob->w[i]);
}
/*
@ -1560,7 +1555,7 @@ namespace VCSnonideal {
plogf("\n"); print_char('-', 80); plogf("\n");
plogf(" Information about phases\n");
plogf(" PhaseName PhaseNum SingSpec GasPhase EqnState NumSpec");
plogf(" TMolesInert Tmoles(gmol)\n");
plogf(" TMolesInert Tmoles(kmol)\n");
for (int iphase = 0; iphase < vprob->NPhase; iphase++) {
vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
@ -1568,8 +1563,8 @@ namespace VCSnonideal {
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID, VolPhase->SingleSpecies,
VolPhase->GasPhase, sEOS.c_str(),
VolPhase->NVolSpecies, VolPhase->TMolesInert * 1.0E3);
plogf("%16e\n", VolPhase->TotalMoles() * 1.0E3);
VolPhase->NVolSpecies, VolPhase->TMolesInert );
plogf("%16e\n", VolPhase->TotalMoles() );
}
plogf("\n"); print_char('=', 80); plogf("\n");

8
Cantera/src/equil/vcs_TP.cpp
View file

@ -36,7 +36,7 @@ namespace VCSnonideal {
* ip1 = 1 -> Print intermediate results.
* maxit -> Maximum number of iterations for the algorithm
* T = Temperature (Kelvin)
* pres = Pressure (units given by if__ variable)
* pres = Pressure (pascal)
*
* Return Codes
* ------------------
@ -55,7 +55,7 @@ namespace VCSnonideal {
* Store the temperature and pressure in the private global variables
*/
m_temperature = T_arg;
m_pressure = pres_arg;
m_pressurePA = pres_arg;
/*
* Evaluate the standard state free energies
* at the current temperatures and pressures.
@ -113,7 +113,7 @@ namespace VCSnonideal {
* 0 -> don't report on anything
* IP1 = 1 -> Print intermediate results.
* T = Temperature (Kelvin)
* Pres = Pressure (units of if__ variable)
* Pres = Pressure (Pascal)
*
* Evaluate the standard state free energies at the current temperature
* and pressure. Ideal gas pressure contribution is added in here.
@ -154,7 +154,7 @@ namespace VCSnonideal {
for (int iph = 0; iph < NPhase; iph++) {
vcs_VolPhase* vph = VPhaseList[iph];
vph->setState_TP(m_temperature, m_pressure);
vph->setState_TP(m_temperature, m_pressurePA);
vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies));
}

26
Cantera/src/equil/vcs_VolPhase.cpp
View file

@ -40,6 +40,7 @@ vcs_VolPhase::vcs_VolPhase() :
ElGlobalIndex(0),
NVolSpecies(0),
TMolesInert(0.0),
m_molarVolInert(1000.),
ActivityConvention(0),
Existence(0),
IndexSpecialSpecies(-1),
@ -402,7 +403,7 @@ double vcs_VolPhase::G0_calc_one(int kspec, double tkelvin) {
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
* @param pres Current pressure (pascal)
*/
void vcs_VolPhase::GStar_calc(double tkelvin, double pres) {
setState_TP(tkelvin, pres);
@ -431,8 +432,8 @@ void vcs_VolPhase::GStar_calc(double tkelvin, double pres) {
* The kth species standard state G is returned
*
* @param kspec Species number (within the phase)
* @param TKelvin Current temperature
* @param pres Current pressure
* @param TKelvin Current temperature (kelvin)
* @param pres Current pressure (pascal)
*
* @return Gstar[kspec] returns the gibbs free energy for the
* standard state of the kth species.
@ -694,7 +695,9 @@ void vcs_VolPhase::setState_TP(double temp, double pres)
* The results are held internally within the object.
*
* @param TKelvin Current temperature
* @param pres Current pressure
* @param pres Current pressure (pascal)
*
* Calculations are in m**3/kmol
*/
void vcs_VolPhase::VolStar_calc(double tkelvin, double pres) {
setState_TP(tkelvin, pres);
@ -722,7 +725,7 @@ void vcs_VolPhase::VolStar_calc(double tkelvin, double pres) {
*
* @param kspec Species number (within the phase)
* @param TKelvin Current temperature
* @param pres Current pressure
* @param pres Current pressure (pascal)
*
* @return molar volume of the kspec species's standard
* state
@ -747,8 +750,7 @@ double vcs_VolPhase::VolPM_calc() const {
kglob = IndSpecies[k];
vcs_SpeciesProperties *sProp = ListSpeciesPtr[k];
VCS_SPECIES_THERMO *sTherm = sProp->SpeciesThermo;
StarMolarVol[k] =
(sTherm->VolStar_calc(kglob, Temp, Pres));
StarMolarVol[k] = (sTherm->VolStar_calc(kglob, Temp, Pres));
}
for (k = 0; k < NVolSpecies; k++) {
PartialMolarVol[k] = StarMolarVol[k];
@ -760,6 +762,16 @@ double vcs_VolPhase::VolPM_calc() const {
Vol += PartialMolarVol[k] * Xmol[k];
}
Vol *= TMoles;
if (TMolesInert > 0.0) {
if (GasPhase) {
double volI = TMolesInert * 8314.47215 * Temp / Pres;
Vol += volI;
} else {
printf("unknown situation\n");
exit(-1);
}
}
}
m_UpToDate_VolPM = true;
return Vol;

24
Cantera/src/equil/vcs_VolPhase.h
View file

@ -215,6 +215,12 @@ public:
//! Total moles of inert in the phase
double TMolesInert;
//! molar volume of the inert species if present
/*!
* units m**3 / kg
*/
double m_molarVolInert;
//! Convention for the activity formulation
/*!
* 0 = molar based activities (default)
@ -285,10 +291,10 @@ public:
*
* Chem_Pot Pres vol moles
* ----------------------------------------------------------------------
* -1 VCS_UNITS_KCALMOL = kcal/mol atm cm**3 gmol
* 0 VCS_UNITS_UNITLESS = MU / RT -> no units atm cm**3 gmol
* 1 VCS_UNITS_KJMOL = kJ / mol atm cm**3 gmol
* 2 VCS_UNITS_KELVIN = KELVIN -> MU / R atm cm**3 gmol
* -1 VCS_UNITS_KCALMOL = kcal/gmol Pa m**3 kmol
* 0 VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol
* 1 VCS_UNITS_KJMOL = kJ / gmol Pa m**3 kmol
* 2 VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol
* 3 VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol
* ----------------------------------------------------------------------
*
@ -358,15 +364,13 @@ public:
//! Vector of the Star molar Volumes of the species.
/*!
* units depends on VCS_UnitsFormat variable
* Cantera -> m3 / kmol
* units m3 / kmol
*/
mutable std::vector<double> StarMolarVol;
//! Vector of the Partial molar Volumes of the species.
/*!
* units depends on VCS_UnitsFormat variable
* Cantera -> m3 / kmol
* units m3 / kmol
*/
mutable std::vector<double> PartialMolarVol;
@ -591,6 +595,8 @@ public:
*
* @param TKelvin Current temperature
* @param pres Current pressure
*
* Units are in m**3/kmol
*/
void VolStar_calc(double TKelvin, double pres);
@ -605,7 +611,7 @@ public:
* @param pres Current pressure
*
* @return molar volume of the kspec species's standard
* state
* state (m**3/kmol)
*/
double VolStar_calc_one(int kglob, double TKelvin, double pres);

16
Cantera/src/equil/vcs_defs.h
View file

@ -191,10 +191,10 @@ namespace VCSnonideal {
* @verbatim
Chem_Pot Pres vol moles
-------------------------------------------------
VCS_UNITS_KCALMOL = kcal/mol atm cm**3 gmol
VCS_UNITS_UNITLESS = MU / RT -> no units atm cm**3 gmol
VCS_UNITS_KJMOL = kJ / mol atm cm**3 gmol
VCS_UNITS_KELVIN = KELVIN -> MU / R atm cm**3 gmol
VCS_UNITS_KCALMOL = kcal/mol Pa m**3 kmol
VCS_UNITS_UNITLESS = MU / RT -> no units Pa m**3 kmol
VCS_UNITS_KJMOL = kJ / mol Pa m**3 kmol
VCS_UNITS_KELVIN = KELVIN -> MU / R Pa m**3 kmol
VCS_UNITS_MKS = Joules / Kmol (Cantera) Pa m**3 kmol
Energy:
@ -205,10 +205,10 @@ namespace VCSnonideal {
VCS_UNITS_MKS = J / kmol
Pressure: (Pref and Pres)
VCS_UNITS_KCALMOL = atm
VCS_UNITS_UNITLESS = no units
VCS_UNITS_KJMOL = atm
VCS_UNITS_KELVIN = atm
VCS_UNITS_KCALMOL = Pa
VCS_UNITS_UNITLESS = Pa
VCS_UNITS_KJMOL = Pa
VCS_UNITS_KELVIN = Pa
VCS_UNITS_MKS = Pa = kg / m s2
@endverbatim
* @{

12
Cantera/src/equil/vcs_nondim.cpp
View file

@ -113,11 +113,13 @@ void VCS_SOLVE::vcs_nondim_TP(void) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= 1.0E3;
//m_molNumSpecies_old[i] *= 1.0E3;
m_molNumSpecies_old[i] *= 1.0;
}
}
for (i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] *= 1.0E3;
//m_elemAbundancesGoal[i] *= 1.0E3;
m_elemAbundancesGoal[i] *= 1.0;
}
}
}
@ -153,11 +155,13 @@ void VCS_SOLVE::vcs_redim_TP(void)
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] /= 1.0E3;
//m_molNumSpecies_old[i] /= 1.0E3;
m_molNumSpecies_old[i] /= 1.0;
}
}
for (i = 0; i < m_numElemConstraints; ++i) {
m_elemAbundancesGoal[i] /= 1.0E3;
//m_elemAbundancesGoal[i] /= 1.0E3;
m_elemAbundancesGoal[i] /= 1.0;
}
}

4
Cantera/src/equil/vcs_prep.cpp
View file

@ -191,8 +191,8 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
}
if (fabs(sum) < 1.0E-6) {
modifiedSoln = true;
if (m_pressure <= 0.0) pres = 1.0;
else pres = m_pressure;
if (m_pressurePA <= 0.0) pres = 1.01325E5;
else pres = m_pressurePA;
retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {

20
Cantera/src/equil/vcs_prob.cpp
View file

@ -42,7 +42,7 @@ VCS_PROB::VCS_PROB(int nsp, int nel, int nph) :
NPhase(nph),
NPHASE0(0),
T(298.15),
Pres(1.0),
PresPA(1.0),
Vol(0.0),
m_VCS_UnitsFormat(VCS_UNITS_UNITLESS),
/* Set the units for the chemical potential data to be
@ -243,10 +243,8 @@ void VCS_PROB::prob_report(int print_lvl) {
if (prob_type == 0) {
plogf("\tSolve a constant T, P problem:\n");
plogf("\t\tT = %g K\n", T);
double pres_atm = Pres;
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
pres_atm = Pres / 1.0133E5;
}
double pres_atm = PresPA / 1.01325E5;
plogf("\t\tPres = %g atm\n", pres_atm);
} else {
plogf("\tUnknown problem type\n");
@ -279,7 +277,7 @@ void VCS_PROB::prob_report(int print_lvl) {
plogf(" Information about phases\n");
plogf(" PhaseName PhaseNum SingSpec GasPhase "
" EqnState NumSpec");
plogf(" TMolesInert Tmoles\n");
plogf(" TMolesInert TKmoles\n");
for (iphase = 0; iphase < NPhase; iphase++) {
Vphase = VPhaseList[iphase];
@ -293,10 +291,11 @@ void VCS_PROB::prob_report(int print_lvl) {
}
plogf("\nElemental Abundances: ");
plogf(" Target_gmol ElemType ElActive\n");
plogf(" Target_kmol ElemType ElActive\n");
double fac = 1.0;
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
fac = 1.0E3;
//fac = 1.0E3;
fac = 1.0;
}
for (i = 0; i < ne; ++i) {
print_space(26); plogf("%-2.2s", ElName[i].c_str());
@ -322,7 +321,7 @@ void VCS_PROB::prob_report(int print_lvl) {
for (iphase = 0; iphase < NPhase; iphase++) {
Vphase = VPhaseList[iphase];
Vphase->G0_calc(T);
Vphase->GStar_calc(T, Pres);
Vphase->GStar_calc(T, PresPA);
for (int kindex = 0; kindex < Vphase->NVolSpecies; kindex++) {
int kglob = Vphase->IndSpecies[kindex];
plogf("%16s ", SpName[kglob].c_str());
@ -475,7 +474,6 @@ void VCS_PROB::reportCSV(const std::string &reportFile) {
exit(-1);
}
double Temp = T;
double pres = Pres;
std::vector<double> volPM(nspecies, 0.0);
std::vector<double> activity(nspecies, 0.0);;
@ -508,7 +506,7 @@ void VCS_PROB::reportCSV(const std::string &reportFile) {
fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
" -----------------------------\n");
fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
fprintf(FP,"Pressure = %11.5g Pascal\n", PresPA);
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations);
fprintf(FP,"Number VCS iterations = %d\n", m_Iterations);

3
Cantera/src/equil/vcs_prob.h
View file

@ -136,8 +136,9 @@ namespace VCSnonideal {
//! Pressure
/*!
* units given by m_VCS_UnitsFormat
* -> are now PA
*/
double Pres;
double PresPA;
//! Volume of the entire system
/*!

24
Cantera/src/equil/vcs_report.cpp
View file

@ -105,23 +105,19 @@ int VCS_SOLVE::vcs_report(int iconv)
* Calculate some quantities that may need updating
*/
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressure,
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
std::string punits = "Atm";
if (m_VCS_UnitsFormat == 3) {
punits = "Pa ";
}
plogf("\t\tPressure = %15.5g %3s\n", m_pressure, punits.c_str());
plogf("\t\tVolume = %15.5g cm**3\n", Vol);
plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
plogf("\t\tVolume = %15.5g m**3\n", Vol);
/*
* -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS --------------
*/
plogf("\n\n");
print_line("-", 80);
plogf(" Species Equilibrium moles ");
plogf(" Species Equilibrium kmoles ");
plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
print_line("-", 80);
for (i = 0; i < m_numComponents; ++i) {
@ -138,7 +134,7 @@ int VCS_SOLVE::vcs_report(int iconv)
if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) {
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l], m_molNumSpecies_new[l], m_feSpecies_curr[l]);
plogf(" MolNum ");
plogf(" KMolNum ");
} else if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_curr[l]);
plogf(" Voltage = %14.7E", m_molNumSpecies_old[l]);
@ -162,13 +158,13 @@ int VCS_SOLVE::vcs_report(int iconv)
}
}
if (m_numSpeciesRdc != nspecies) {
plogf("\n SPECIES WITH LESS THAN 1.0E-32 MOLES:\n\n");
plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
plogf(" %-12.12s", SpName[kspec].c_str());
plogf(" %14.7E %14.7E %12.4E",
m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" Mol_Num");
plogf(" KMol_Num");
} else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Voltage");
} else {
@ -237,7 +233,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %10.10s", (ElName[j]).c_str());
}
plogf(" | |\n");
plogf(" PhaseName | MolTarget |");
plogf(" PhaseName |KMolTarget |");
for (j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", m_elemAbundancesGoal[j]);
}
@ -291,7 +287,7 @@ int VCS_SOLVE::vcs_report(int iconv)
if (inertYes)
plogf("\t\t(Inert species have standard free energy of zero)\n");
plogf("\nElemental Abundances: ");
plogf("\nElemental Abundances (kmol): ");
plogf(" Actual Target Type ElActive\n");
for (i = 0; i < m_numElemConstraints; ++i) {
print_space(26); plogf("%-2.2s", (ElName[i]).c_str());
@ -310,7 +306,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\t\t(RT = %g ", rt);
vcs_printChemPotUnits(m_VCS_UnitsFormat);
plogf(")\n");
plogf(" Name TMoles StandStateChemPot "
plogf(" Name TKMoles StandStateChemPot "
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)\n");
print_line("-", 115);
for (i = 0; i < nspecies; ++i) {

14
Cantera/src/equil/vcs_solve.cpp
View file

@ -44,7 +44,7 @@ namespace VCSnonideal {
iest(0),
TMoles(0.0),
m_temperature(0.0),
m_pressure(0.0),
m_pressurePA(0.0),
tolmaj(0.0),
tolmin(0.0),
tolmaj2(0.0),
@ -389,7 +389,7 @@ namespace VCSnonideal {
* calculate the residual and Jacobian)
*/
iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->Pres);
iconv = vcs_TP(ipr, ip1, maxit, vprob->T, vprob->PresPA);
/*
@ -582,8 +582,8 @@ namespace VCSnonideal {
*/
if (pub->T > 0.0) m_temperature = pub->T;
else m_temperature = 293.15;
if (pub->Pres > 0.0) m_pressure = pub->Pres;
else m_pressure = 1.0;
if (pub->PresPA > 0.0) m_pressurePA = pub->PresPA;
else m_pressurePA = Cantera::OneAtm;
/*
* TPhInertMoles[] -> must be copied over here
*/
@ -773,7 +773,7 @@ namespace VCSnonideal {
bool status_change = false;
m_temperature = pub->T;
m_pressure = pub->Pres;
m_pressurePA = pub->PresPA;
m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
iest = pub->iest;
@ -902,7 +902,7 @@ namespace VCSnonideal {
int k1 = 0;
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressure,
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
for (i = 0; i < m_numSpeciesTot; ++i) {
@ -930,7 +930,7 @@ namespace VCSnonideal {
}
pub->T = m_temperature;
pub->Pres = m_pressure;
pub->PresPA = m_pressurePA;
pub->Vol = Vol;
int kT = 0;
for (int iph = 0; iph < pub->NPhase; iph++) {

25
Cantera/src/equil/vcs_solve.h
View file

@ -621,7 +621,7 @@ public:
* Note other constraint conditions are added to this vector.
* This is input from the input file and
* is considered a constant from thereon.
* units = gmoles
* units = kmoles
*/
std::vector<double> m_elemAbundances;
@ -632,7 +632,7 @@ public:
* Note other constraint conditions are added to this vector.
* This is input from the input file and
* is considered a constant from thereon.
* units = gmoles
* units = kmoles
*/
std::vector<double> m_elemAbundancesGoal;
@ -641,7 +641,7 @@ public:
* -> Don't use this except for scaling
* purposes only */
//! total gmols of species in each phase
//! total kmols of species in each phase
/*!
* This contains the total number of moles of species in each phase
*
@ -649,7 +649,7 @@ public:
*/
std::vector<double> m_tPhaseMoles_old;
//! total gmols of species in each phase in the tentative soln vector
//! total kmols of species in each phase in the tentative soln vector
/*!
* This contains the total number of moles of species in each phase
* in the tentative solution vector
@ -681,12 +681,13 @@ public:
* 1: atm
* 2: atm
* 3: Pa
* Units being changed to Pa
*/
double m_pressure;
double m_pressurePA;
//! Total kmoles of inert to add to each phase
/*!
* TPhInertMoles[iph] = Total gmoles of inert to add to each phase
* TPhInertMoles[iph] = Total kmoles of inert to add to each phase
* length = number of phases
*/
std::vector<double> TPhInertMoles;
@ -864,7 +865,7 @@ public:
//! Molecular weight of each species
/*!
* units = gm/gmol
* units = kg/kmol
* length = number of species
*/
std::vector<double> WtSpecies;
@ -891,7 +892,7 @@ public:
*/
int UseActCoeffJac;
double Vol; /* Vol = Volume (cm^3) */
double Vol; /* Vol = Volume (m^3) */
//! Partialm molar volumes of the species
/*!
@ -928,10 +929,10 @@ public:
* 2: Kelvin
* 3: J / kmol
* and pressure data:
* -1: atm
* 0: atm
* 1: atm
* 2: atm
* -1: Pa
* 0: Pa
* 1: Pa
* 2: pa
* 3: Pa
*/
int m_VCS_UnitsFormat;

36
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -237,11 +237,11 @@ namespace VCSnonideal {
plogf("%10d PHASE2 SPECIES%8d SINGLE SPECIES PHASES\n\n",
numSpecliquid,
m_numSpeciesTot - (VPhaseList[0])->NVolSpecies - numSpecliquid);
string punits = "atm";
if (m_VCS_UnitsFormat == 3) {
punits = "Pa ";
}
plogf(" PRESSURE%22.8g %3s\n", m_pressure, punits.c_str());
//string punits = "atm";
//if (m_VCS_UnitsFormat == 3) {
// punits = "Pa ";
//}
plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa ");
plogf(" TEMPERATURE%19.3f K\n", m_temperature);
Vphase = VPhaseList[0];
if (Vphase->NVolSpecies > 0) {
@ -532,7 +532,7 @@ namespace VCSnonideal {
if (iti == 0) plogf("- Full Calculation:\n");
else plogf("- Major Components Calculation:\n");
plogf(" --- Species IC ");
plogf(" Moles Tent_Moles Rxn_Adj | Comment \n");
plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n");
}
#endif
@ -785,7 +785,7 @@ namespace VCSnonideal {
*/
if (m_molNumSpecies_new[kspec] <= 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "initial nonpos moles= %11.3E",
sprintf(ANOTE, "initial nonpos kmoles= %11.3E",
m_molNumSpecies_new[kspec]);
#endif
/* ************************************************* */
@ -801,7 +801,7 @@ namespace VCSnonideal {
if (! (SSPhase[kspec])) {
/*
* Section for multispecies phases:
* - Cut reaction adjustment for positive moles of
* - Cut reaction adjustment for positive kmoles of
* major species in multispecies phases.
* Decrease its concentration by a factor of 10.
*/
@ -811,7 +811,7 @@ namespace VCSnonideal {
/*
* Change major to minor if the current species
* has a mole number that is less than 1/100 of the
* total moles in the problem.
* total kmoles in the problem.
* However, it also has to be a small species within its
* own phase as well.
* we can't call vcs_species_type() because the phase moles
@ -954,7 +954,7 @@ namespace VCSnonideal {
} /* End of Loop on ic[irxn] -> the type of species */
/***********************************************************************/
/****** CALCULATE MOLE NUMBER CHANGE FOR THE COMPONENT BASIS ***********/
/****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/
/***********************************************************************/
if (dx != 0.0 && (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
/*
@ -1143,7 +1143,7 @@ namespace VCSnonideal {
plogf(" -----------------------------------------------------\n");
plogf(" --- FORCER SUBROUTINE changed the solution:\n");
plogf(" --- SPECIES Status INIT MOLES TENT_MOLES");
plogf(" FINAL MOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
plogf(" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
for (i = 0; i < m_numComponents; ++i) {
plogf(" --- %-12.12s", SpName[i].c_str());
plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i], m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]);
@ -1160,11 +1160,11 @@ namespace VCSnonideal {
plogf("Norms of Delta G():%14.6E%14.6E\n",
l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)),
l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
plogf(" Total moles of gas = %15.7E\n", m_tPhaseMoles_old[0]);
plogf(" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]);
if ((NPhase > 1) && (! (VPhaseList[1])->SingleSpecies)) {
plogf(" Total moles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
plogf(" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]);
} else {
plogf(" Total moles of liquid = %15.7E\n", 0.0);
plogf(" Total kmoles of liquid = %15.7E\n", 0.0);
}
plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
@ -1185,7 +1185,7 @@ namespace VCSnonideal {
} else {
plogf(" (only major species):\n");
}
plogf(" --- Species Status Initial_Moles Final_Moles Initial_Mu/RT");
plogf(" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n");
for (i = 0; i < m_numComponents; ++i) {
plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" ");
@ -1214,7 +1214,7 @@ namespace VCSnonideal {
l2normdg(VCS_DATA_PTR(m_deltaGRxn_new)));
plogendl();
plogf(" --- Phase_Name Moles(after update)\n");
plogf(" --- Phase_Name KMoles(after update)\n");
plogf(" --- "); vcs_print_line("-", 50);
for (iph = 0; iph < NPhase; iph++) {
Vphase = VPhaseList[iph];
@ -2129,7 +2129,7 @@ namespace VCSnonideal {
*
* Change the concentration of a species by delta moles.
* Make sure to conserve
* elements and keep track of the total moles in all phases.
* elements and keep track of the total kmoles in all phases.
*
* return:
* 1: succeeded
@ -2831,7 +2831,7 @@ namespace VCSnonideal {
plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n");
plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n");
plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n");
plogf(" --- Species Moles Rxn_Adjustment DeltaG"
plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
" | Comment\n");
}
#endif

27
Cantera/src/equil/vcs_species_thermo.cpp
View file

@ -43,7 +43,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase,
SS0_H0(0.0),
SS0_S0(0.0),
SS0_Cp0(0.0),
SS0_Pref(1.0),
SS0_Pref(1.01325E5),
SS0_Params(0),
SSStar_Model(VCS_SSSTAR_CONSTANT),
SSStar_Params(0),
@ -55,15 +55,7 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase,
UseCanteraCalls(false),
m_VCS_UnitsFormat(VCS_UNITS_UNITLESS)
{
/*
* Set up the numerical value for P_reference, based on the current
* global units choice.
*/
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
SS0_Pref = 1.0133E5;
} else {
SS0_Pref = 1.0;
}
SS0_Pref = 1.01325E5;
}
@ -224,11 +216,11 @@ double VCS_SPECIES_THERMO::GStar_R_calc(int kglob, double TKelvin,
* Input
*
* Output
* return value = standard state volume in cm**3 per mol.
* (VCS_UNITS_MKS) m**3 / kmol
* return value = standard state volume in m**3 per kmol.
* (VCS_UNITS_MKS)
*/
double VCS_SPECIES_THERMO::
VolStar_calc(int kglob, double TKelvin, double pres)
VolStar_calc(int kglob, double TKelvin, double presPA)
{
char yo[] = "VCS_SPECIES_THERMO::VStar_calc ";
double vol, T;
@ -237,18 +229,15 @@ VolStar_calc(int kglob, double TKelvin, double pres)
if (UseCanteraCalls) {
AssertThrowVCS(m_VCS_UnitsFormat == VCS_UNITS_MKS, "Possible inconsistency");
int kspec = IndexSpeciesPhase;
vol = OwningPhase->VolStar_calc_one(kspec, TKelvin, pres);
vol = OwningPhase->VolStar_calc_one(kspec, TKelvin, presPA);
} else {
switch(SSStar_Vol_Model) {
case VCS_SSVOL_CONSTANT:
vol = SSStar_Vol0;
break;
case VCS_SSVOL_IDEALGAS:
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
vol = 8.31451E3 * T / pres;
} else {
vol= 83.14510 / 1.01325 * T / pres;
}
// R J/kmol/K (2006 CODATA value)
vol= 8314.47215 * T / presPA;
break;
default:
plogf("%sERROR: unknown SSVol model\n", yo);

1
Cantera/src/equil/vcs_species_thermo.h
View file

@ -112,6 +112,7 @@ public:
/**
* Value of the pressure for the reference state.
* defaults to 1.01325E5 = 1 atm
*/
double SS0_Pref;
/**