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

Fixed an error in the phase stability calculation

Browse source
This commit is contained in:
Harry Moffat 2010-12-18 18:26:16 +00:00
parent fc7b7e8dc6
commit 8c4d47f59c
6 changed files with 242 additions and 18 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -507,6 +507,10 @@ namespace VCSnonideal {
const std::vector<double> & vcs_VolPhase::moleFractions() const {
return Xmol_;
}
double vcs_VolPhase::moleFraction(int k) const {
return Xmol_[k];
}
/***************************************************************************/
// Set the moles and/or mole fractions within the phase

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

@ -419,6 +419,8 @@ namespace VCSnonideal {
//! object.
const std::vector<double> & moleFractions() const;
double moleFraction(int klocal) const;
//! Sets the creationMoleNum's within the phase object
/*!
* @param F_k Pointer to a vector of n_k's

13
Cantera/src/equil/vcs_phaseStability.cpp
View file

@ -671,7 +671,7 @@ namespace VCSnonideal {
vector<doublereal> fracDelta_old(Vphase->nSpecies(), 0.0);
vector<doublereal> fracDelta_raw(Vphase->nSpecies(), 0.0);
vector<int> creationGlobalRxnNumbers(Vphase->nSpecies(), -1);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
vector<doublereal> m_feSpecies_Deficient(m_numComponents, 0.0);
doublereal damp = 1.0;
@ -779,18 +779,16 @@ namespace VCSnonideal {
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
/*
* first Calculate altered chemical potentials for component species
* First calculate altered chemical potentials for component species
* belonging to this phase.
*/
for (i = 0; i < (int) componentList.size(); i++) {
kc = componentList[i];
kc_spec = Vphase->spGlobalIndexVCS(kc);
if ( X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
+ log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
} else {
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
+ log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
}
}
@ -807,8 +805,7 @@ namespace VCSnonideal {
}
double *dtmp_ptr = m_stoichCoeffRxnMatrix[irxn];
if (dtmp_ptr[kc_spec] != 0.0) {
m_deltaGRxn_Deficient[irxn] +=
dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
m_deltaGRxn_Deficient[irxn] += dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
}
}

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

@ -495,13 +495,20 @@ public:
*/
void vcs_dfe(const int stateCalc, const int ll, const int lbot, const int ltop);
//! Print out a table of chemical potentials
/*!
* @param vcsState Determines where to get the mole numbers from.
* - VCS_STATECALC_OLD -> from m_molNumSpecies_old
* - VCS_STATECALC_NEW -> from m_molNumSpecies_new
*/
void vcs_printSpeciesChemPot(const int stateCalc) const;
//! This routine uploads the state of the system into all of the
//! vcs_VolumePhase objects in the current problem.
/*!
* @param vcsState Determines where to get the mole numbers from.
* - VCS_STATECALC_OLD -> from m_molNumSpecies_old
* - VCS_STATECALC_NEW -> from m_molNumSpecies_new
*
*/
void vcs_updateVP(const int stateCalc);
@ -621,6 +628,8 @@ public:
void vcs_deltag(const int l, const bool doDeleted, const int vcsState,
const bool alterZeroedPhases = true);
void vcs_printDeltaG(const int stateCalc);
//! Calculate deltag of formation for all species in a single phase.
/*!
* Calculate deltag of formation for all species in a single
@ -1618,8 +1627,10 @@ public:
//! Last deltag[irxn] from the previous step
std::vector<double> m_deltaGRxn_old;
//! Last deltag[irxn] from the previous step with additions for
//! possible births of zeroed phases.
//! Last deltag[irxn] from the previous step with additions for possible births of zeroed phases for component species
/*!
*
*/
std::vector<double> m_deltaGRxn_Deficient;
//! Temporary vector of Rxn DeltaG's
@ -1684,10 +1695,10 @@ public:
std::vector<double> m_tPhaseMoles_new;
//! Temporary vector of length NPhase
std::vector<double> m_TmpPhase;
mutable std::vector<double> m_TmpPhase;
//! Temporary vector of length NPhase
std::vector<double> m_TmpPhase2;
mutable std::vector<double> m_TmpPhase2;
//! Change in the total moles in each phase
/*!

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

@ -4371,10 +4371,10 @@ namespace VCSnonideal {
feSpecies[kspec] = m_SSfeSpecies[kspec]
+ log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[m_phaseID[kspec]] - m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
} else {
feSpecies[kspec] = m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase]; ;
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iphase];
}
} else {
feSpecies[kspec] = m_SSfeSpecies[kspec]
@ -4438,7 +4438,93 @@ namespace VCSnonideal {
}
}
}
}
//====================================================================================================================
// Print out a table of chemical potentials
/*
* @param vcsState Determines where to get the mole numbers from.
* - VCS_STATECALC_OLD -> from m_molNumSpecies_old
* - VCS_STATECALC_NEW -> from m_molNumSpecies_new
*/
void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const {
double mfValue = 1.0;
bool zeroedPhase = false;
int kspec;
const double * molNum = VCS_DATA_PTR(m_molNumSpecies_old);
const double * tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
const double * actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
if (stateCalc == VCS_STATECALC_NEW) {
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
molNum = VCS_DATA_PTR(m_molNumSpecies_new);
}
double * tMoles = VCS_DATA_PTR(m_TmpPhase);
const double *tPhInertMoles = VCS_DATA_PTR(TPhInertMoles);
for (int iph = 0; iph < m_numPhases; iph++) {
tMoles[iph] = tPhInertMoles[iph];
}
for (kspec = 0; kspec < m_numSpeciesTot; kspec++) {
if(m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
int iph = m_phaseID[kspec];
tMoles[iph] += molNum[kspec];
}
}
double RT = m_temperature * Cantera::GasConstant;
printf(" --- CHEMICAL POT TABLE (J/kmol) Name PhID MolFR ChemoSS "
" logMF Gamma Elect extra ElectrChem\n");
printf(" ");
vcs_print_line("-", 132);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
mfValue = 1.0;
int iphase = m_phaseID[kspec];
const vcs_VolPhase * Vphase = m_VolPhaseList[iphase];
if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS) ) {
zeroedPhase = true;
} else {
zeroedPhase = false;
}
if (tMoles[iphase] > 0.0) {
if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tMoles[iphase];
} else {
mfValue = molNum[kspec]/tMoles[iphase];
}
} else {
int klocal = m_speciesLocalPhaseIndex[kspec];
mfValue = Vphase->moleFraction(klocal);
}
double volts = Vphase->electricPotential();
double elect = m_chargeSpecies[kspec] * m_Faraday_dim * volts;
double comb = - m_lnMnaughtSpecies[kspec];
double total = (m_SSfeSpecies[kspec] + log(mfValue) + elect + log(actCoeff_ptr[kspec]) + comb);
if (zeroedPhase) {
printf(" --- ** zp *** ");
} else {
printf(" --- ");
}
printf("%-24.24s", m_speciesName[kspec].c_str());
printf(" %-3d", iphase);
printf(" % -12.4e", mfValue);
printf(" % -12.4e", m_SSfeSpecies[kspec] * RT);
printf(" % -12.4e", log(mfValue) * RT);
printf(" % -12.4e", log(actCoeff_ptr[kspec]) * RT);
printf(" % -12.4e", elect * RT);
printf(" % -12.4e", comb * RT);
printf(" % -12.4e\n", total *RT);
}
printf(" ");
vcs_print_line("-", 132);
}
/*****************************************************************************/
#ifdef DEBUG_MODE
@ -4990,8 +5076,125 @@ namespace VCSnonideal {
}
#endif
}
/*****************************************************************************/
//====================================================================================================================
void VCS_SOLVE::vcs_printDeltaG( const int stateCalc) {
int j;
double * deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
double * feSpecies = VCS_DATA_PTR(m_feSpecies_old);
double * molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
const double * tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
const double * actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
if (stateCalc == VCS_STATECALC_NEW) {
deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
feSpecies = VCS_DATA_PTR(m_feSpecies_new);
molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
}
double RT = m_temperature * Cantera::GasConstant;
bool zeroedPhase = false;
if (m_debug_print_lvl >= 2) {
plogf(" --- DELTA_G TABLE Components:");
for (j = 0; j < m_numComponents; j++) {
plogf(" %3d ", j);
}
plogf("\n --- Components Moles:");
for (j = 0; j < m_numComponents; j++) {
plogf("%10.3g", m_molNumSpecies_old[j]);
}
plogf("\n --- NonComponent| Moles | ");
for (j = 0; j < m_numComponents; j++) {
plogf("%-10.10s", m_speciesName[j].c_str());
}
//plogf("| m_scSize");
plogf("\n");
for (int i = 0; i < m_numRxnTot; i++) {
plogf(" --- %3d ", m_indexRxnToSpecies[i]);
plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[i]].c_str());
plogf("|%10.3g|", m_molNumSpecies_old[m_indexRxnToSpecies[i]]);
for (j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]);
}
//plogf(" | %6.2f", m_scSize[i]);
plogf("\n");
}
plogf(" "); for(int i=0; i<77; i++) plogf("-"); plogf("\n");
}
printf(" --- DeltaG Table (J/kmol) Name PhID MoleNum MolFR "
" ElectrChemStar ElectrChem DeltaGStar DeltaG(Pred) Stability\n");
printf(" ");
vcs_print_line("-", 132);
for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) {
int irxn = kspec - m_numComponents;
double mfValue = 1.0;
int iphase = m_phaseID[kspec];
const vcs_VolPhase * Vphase = m_VolPhaseList[iphase];
if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS) ) {
zeroedPhase = true;
} else {
zeroedPhase = false;
}
if (tPhMoles_ptr[iphase] > 0.0) {
if (molNumSpecies[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
mfValue = VCS_DELETE_MINORSPECIES_CUTOFF / tPhMoles_ptr[iphase];
} else {
mfValue = molNumSpecies[kspec] / tPhMoles_ptr[iphase];
}
} else {
int klocal = m_speciesLocalPhaseIndex[kspec];
mfValue = Vphase->moleFraction(klocal);
}
if (zeroedPhase) {
printf(" --- ** zp *** ");
} else {
printf(" --- ");
}
double feFull = feSpecies[kspec];
if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ) {
feFull += log(actCoeff_ptr[kspec]) + log(mfValue);
}
printf("%-24.24s", m_speciesName[kspec].c_str());
printf(" %-3d", iphase);
printf(" % -12.4e", molNumSpecies[kspec]);
printf(" % -12.4e", mfValue);
printf(" % -12.4e", feSpecies[kspec] * RT);
printf(" % -12.4e", feFull * RT);
if (irxn >= 0) {
printf(" % -12.4e", deltaGRxn[irxn] * RT);
printf(" % -12.4e", (deltaGRxn[irxn] + feFull - feSpecies[kspec]) * RT);
if (deltaGRxn[irxn] < 0.0) {
if ( molNumSpecies[kspec] > 0.0) {
printf(" growing");
} else {
printf(" stable");
}
} else if (deltaGRxn[irxn] > 0.0) {
if ( molNumSpecies[kspec] > 0.0) {
printf(" shrinking");
} else {
printf(" unstable");
}
} else {
printf(" balanced");
}
}
printf(" \n");
}
printf(" ");
vcs_print_line("-", 132);
}
//====================================================================================================================
// Calculate deltag of formation for all species in a single phase.
/*
* Calculate deltag of formation for all species in a single

11
Cantera/src/equil/vcs_solve_phaseStability.cpp
View file

@ -203,7 +203,7 @@ namespace VCSnonideal {
*/
int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc,
double &funcVal,
int printLv) {
int printLvl) {
int retn = 0;
double test = -1.0E-10;
int usedZeroedSpecies;
@ -222,11 +222,18 @@ namespace VCSnonideal {
retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa),
VCS_DATA_PTR(sm), VCS_DATA_PTR(ss),
test, &usedZeroedSpecies);
vcs_evaluate_speciesType();
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
if (printLvl > 3) {
vcs_printSpeciesChemPot(VCS_STATECALC_OLD);
}
vcs_deltag(0, true, VCS_STATECALC_OLD);
if (printLvl > 3) {
vcs_printDeltaG(VCS_STATECALC_OLD);
}
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
phasePopPhaseIDs.clear();
iphasePop = vcs_popPhaseID(phasePopPhaseIDs);
funcVal = vcs_phaseStabilityTest(iph);