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

Changed the algorithms a little bit to accomodate PhaseCombo, which

Browse source 
doesn't have an ideal mixing term. This caused problems in places.
This commit is contained in:
Harry Moffat 2011-04-14 18:25:24 +00:00
parent 59a800218f
commit 5f73e33c9e
3 changed files with 104 additions and 38 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -123,7 +123,7 @@ namespace Cantera {
//! Add a phase to the mixture.
/*!
* This function must be called befure the init() function is called,
* This function must be called before the init() function is called,
* which serves to freeze the MultiPhase.
*
* @param p pointer to the phase object

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

@ -984,15 +984,30 @@ namespace VCSnonideal {
*/
void vcs_VolPhase::_updateLnActCoeffJac() {
int k, j;
double deltaMoles_j = 0.0;
/*
* Evaluate the current base activity coefficients if necessary
*/
if (!m_UpToDate_AC) {
_updateActCoeff();
}
#ifndef NOOLD
if (!TP_ptr) return;
TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &dLnActCoeffdMolNumber[0][0]);
for (j = 0; j < m_numSpecies; j++) {
double moles_j_base = v_totalMoles * Xmol_[j];
double * const lnActCoeffCol = dLnActCoeffdMolNumber[j];
if (moles_j_base < 1.0E-200) {
moles_j_base = 1.0E-7 * moles_j_base + 1.0E-20 * v_totalMoles + 1.0E-150;
}
for (k = 0; k < m_numSpecies; k++) {
lnActCoeffCol[k] /= moles_j_base;
}
}
#endif
double deltaMoles_j = 0.0;
// Make copies of ActCoeff and Xmol_ for use in taking differences
std::vector<double> ActCoeff_Base(ActCoeff);
std::vector<double> Xmol_Base(Xmol_);
@ -1008,7 +1023,7 @@ namespace VCSnonideal {
* quantities.
*/
double moles_j_base = v_totalMoles * Xmol_Base[j];
deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-20 * v_totalMoles + 1.0E-150;
deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * v_totalMoles + 1.0E-150;
/*
* Now, update the total moles in the phase and all of the
* mole fractions based on this.
@ -1030,8 +1045,15 @@ namespace VCSnonideal {
*/
double * const lnActCoeffCol = dLnActCoeffdMolNumber[j];
for (k = 0; k < m_numSpecies; k++) {
lnActCoeffCol[k] = (ActCoeff[k] - ActCoeff_Base[k]) /
double tmp;
tmp = (ActCoeff[k] - ActCoeff_Base[k]) /
((ActCoeff[k] + ActCoeff_Base[k]) * 0.5 * deltaMoles_j);
if (fabs(tmp - lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(lnActCoeffCol[k])) {
// printf(" we have an error\n");
}
//tmp = lnActCoeffCol[k];
}
/*
* Revert to the base case Xmol_, v_totalMoles
@ -1048,6 +1070,7 @@ namespace VCSnonideal {
setMoleFractions(VCS_DATA_PTR(Xmol_Base));
_updateMoleFractionDependencies();
_updateActCoeff();
}
/***************************************************************************/

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

@ -1243,11 +1243,12 @@ namespace VCSnonideal {
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new),
VCS_DATA_PTR(m_tPhaseMoles_new)));
plogendl();
#ifdef DEBUG_MODE
if (m_VCount->Its > 550) {
plogf(" --- Troublesome solve");
plogendl();
}
#endif
}
/*************************************************************************/
@ -2056,12 +2057,15 @@ namespace VCSnonideal {
, char *ANOTE
#endif
) const {
double dx = 0.0;
double dx = 0.0, a;
double w_kspec = m_molNumSpecies_old[kspec];
double molNum_kspec_new;
double wTrial;
double dg_irxn = m_deltaGRxn_old[irxn];
double dg_irxn = m_deltaGRxn_old[irxn];
doublereal s;
vcs_VolPhase * Vphase = 0;
int iph = m_phaseID[kspec];
*do_delete = FALSE;
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (w_kspec <= 0.0) {
@ -2086,8 +2090,36 @@ namespace VCSnonideal {
return 0.0;
}
}
wTrial = w_kspec * exp(-dg_irxn);
/*
* get the diagonal of the activity coefficent jacobian
*/
Vphase = m_VolPhaseList[iph];
s = m_dLnActCoeffdMolNum[kspec][kspec];
// s *= (m_tPhaseMoles_old[iph]);
/*
* We fit it to a power law approximation of the activity coefficient
*
* gamma = gamma_0 * ( x / x0)**a
*
* where a is forced to be a little bit greater than -1.
* We do this so that the resulting expression is always nonnegative
*
* We then solve the resulting calculation:
*
* gamma * x = gamma_0 * x0 exp (-deltaG/RT);
*
*
*/
a = w_kspec * s;
if (a < (-1.0 + 0.01)) {
a = -1.0 + 0.01;
} else if (a > 1.0) {
a = 1.0;
}
wTrial = w_kspec * (exp( -dg_irxn / (1.0 + a)));
// wTrial = w_kspec * exp(-dg_irxn);
molNum_kspec_new = wTrial;
if (wTrial > 100. * w_kspec) {
@ -2785,10 +2817,8 @@ namespace VCSnonideal {
}
return false;
}
/*************************************************************************************/
// Provide an estimate for the deleted species in phases that
// are not zeroed out
//====================================================================================================================
// Provide an estimate for the deleted species in phases that are not zeroed out
/*
* Try to add back in all deleted species. An estimate of the kmol numbers
* are obtained and the species is added back into the equation system,
@ -2798,38 +2828,53 @@ namespace VCSnonideal {
int iph, kspec, retn;
if (m_numSpeciesRdc == m_numSpeciesTot) return 0;
/*
* Use the standard chemical potentials for the chemical potentials
* of deleted species. Then, calculate Delta G for
* Use the standard chemical potentials for the chemical potentials of deleted species. Then, calculate Delta G for
* for formation reactions.
* We are relying here on a old saved value of m_actCoeffSpecies_old[kspec]
* being sufficiently good. Note, we will recalculate everything at the
* end of the routine.
* being sufficiently good. Note, we will recalculate everything at the end of the routine.
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
iph = m_phaseID[kspec];
m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
- m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), m_numSpeciesTot);
for (int cits = 0; cits < 3; cits++) {
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
iph = m_phaseID[kspec];
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
if (m_molNumSpecies_new[kspec] == 0.0) {
m_molNumSpecies_new[kspec] = VCS_DELETE_MINORSPECIES_CUTOFF * 1.0E-10;
}
if (!Vphase->m_singleSpecies) {
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_NEW, VCS_DATA_PTR(m_actCoeffSpecies_new));
}
m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_new[kspec]) - m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
}
/*
* Recalculate the DeltaG's of the formation reactions for the deleted species in the mechanism
*/
vcs_deltag(0, true, VCS_STATECALC_NEW);
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 690.0);
double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
m_molNumSpecies_new[kspec] = dx;
if (m_molNumSpecies_new[kspec] > 2 *VCS_DELETE_MINORSPECIES_CUTOFF) {
m_molNumSpecies_new[kspec] = 2 * VCS_DELETE_MINORSPECIES_CUTOFF;
}
}
}
}
/*
* Recalculate the DeltaG's of the formation reactions for the
* deleted species in the mechanism
*/
vcs_deltag(0, true, VCS_STATECALC_NEW);
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 690.0);
double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
double dx = m_molNumSpecies_new[kspec];
retn = delta_species(kspec, &dx);
if (retn == 0) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for "
"species %s (%d) with mol number %g\n",
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec].c_str(), kspec, dx);
}
#endif
@ -2839,8 +2884,7 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (retn == 0) {
if (m_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for "
"species %s (%d) with mol number %g\n",
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec].c_str(), kspec, dx);
}
}
@ -2877,8 +2921,7 @@ namespace VCSnonideal {
retn++;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- add_deleted(): species %s "
"with mol number %g not converged: DG = %g",
plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
m_deltaGRxn_old[irxn]);
plogendl();