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

Reworking of cropping algorithm to adhere to the sand report.

Browse source
This commit is contained in:
Harry Moffat 2009-01-17 04:40:28 +00:00
parent 6e447bf49c
commit 793daa0194
4 changed files with 284 additions and 110 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

304
Cantera/src/thermo/HMWSoln.cpp
View file

@ -56,7 +56,7 @@ namespace Cantera {
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(1.0E-5),
IMS_gamma_k_min_(10.0),
IMS_cCut_(0.5),
IMS_cCut_(0.05),
IMS_slopefCut_(0.6),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -67,6 +67,14 @@ namespace Cantera {
IMS_egCut_(0.0),
IMS_agCut_(0.0),
IMS_bgCut_(0.0),
MC_X_o_cutoff_(0.0),
MC_X_o_min_(0.0),
MC_slopepCut_(0.0),
MC_dpCut_(0.0),
MC_epCut_(0.0),
MC_apCut_(0.0),
MC_bpCut_(0.0),
MC_cpCut_(0.0),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -101,7 +109,7 @@ namespace Cantera {
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(1.0E-5),
IMS_gamma_k_min_(10.0),
IMS_cCut_(0.5),
IMS_cCut_(0.05),
IMS_slopefCut_(0.6),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -112,6 +120,14 @@ namespace Cantera {
IMS_egCut_(0.0),
IMS_agCut_(0.0),
IMS_bgCut_(0.0),
MC_X_o_cutoff_(0.0),
MC_X_o_min_(0.0),
MC_slopepCut_(0.0),
MC_dpCut_(0.0),
MC_epCut_(0.0),
MC_apCut_(0.0),
MC_bpCut_(0.0),
MC_cpCut_(0.0),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -140,7 +156,7 @@ namespace Cantera {
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(1.0E-5),
IMS_gamma_k_min_(10.0),
IMS_cCut_(0.5),
IMS_cCut_(0.05),
IMS_slopefCut_(0.6),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -151,6 +167,14 @@ namespace Cantera {
IMS_egCut_(0.0),
IMS_agCut_(0.0),
IMS_bgCut_(0.0),
MC_X_o_cutoff_(0.0),
MC_X_o_min_(0.0),
MC_slopepCut_(0.0),
MC_dpCut_(0.0),
MC_epCut_(0.0),
MC_apCut_(0.0),
MC_bpCut_(0.0),
MC_cpCut_(0.0),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -185,7 +209,7 @@ namespace Cantera {
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(1.0E-5),
IMS_gamma_k_min_(10.0),
IMS_cCut_(0.5),
IMS_cCut_(0.05),
IMS_slopefCut_(0.6),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -196,6 +220,14 @@ namespace Cantera {
IMS_egCut_(0.0),
IMS_agCut_(0.0),
IMS_bgCut_(0.0),
MC_X_o_cutoff_(0.0),
MC_X_o_min_(0.0),
MC_slopepCut_(0.0),
MC_dpCut_(0.0),
MC_epCut_(0.0),
MC_apCut_(0.0),
MC_bpCut_(0.0),
MC_cpCut_(0.0),
m_debugCalc(0)
{
/*
@ -340,6 +372,14 @@ namespace Cantera {
IMS_egCut_ = b.IMS_egCut_;
IMS_agCut_ = b.IMS_agCut_;
IMS_bgCut_ = b.IMS_bgCut_;
MC_X_o_cutoff_ = b.MC_X_o_cutoff_;
MC_X_o_min_ = b.MC_X_o_min_;
MC_slopepCut_ = b.MC_slopepCut_;
MC_dpCut_ = b.MC_dpCut_;
MC_epCut_ = b.MC_epCut_;
MC_apCut_ = b.MC_apCut_;
MC_bpCut_ = b.MC_bpCut_;
MC_cpCut_ = b.MC_cpCut_;
m_CounterIJ = b.m_CounterIJ;
m_molalitiesCropped = b.m_molalitiesCropped;
@ -390,6 +430,30 @@ namespace Cantera {
m_waterSS(0),
m_densWaterSS(1000.),
m_waterProps(0),
m_molalitiesAreCropped(false),
IMS_typeCutoff_(0),
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(1.0E-5),
IMS_gamma_k_min_(10.0),
IMS_cCut_(0.05),
IMS_slopefCut_(0.6),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
IMS_afCut_(0.0),
IMS_bfCut_(0.0),
IMS_slopegCut_(0.0),
IMS_dgCut_(0.0),
IMS_egCut_(0.0),
IMS_agCut_(0.0),
IMS_bgCut_(0.0),
MC_X_o_cutoff_(0.0),
MC_X_o_min_(0.0),
MC_slopepCut_(0.0),
MC_dpCut_(0.0),
MC_epCut_(0.0),
MC_apCut_(0.0),
MC_bpCut_(0.0),
MC_cpCut_(0.0),
m_debugCalc(0)
{
if (testProb != 1) {
@ -1733,6 +1797,17 @@ namespace Cantera {
*/
s_updatePitzer_lnMolalityActCoeff();
double xmolSolvent = moleFraction(m_indexSolvent);
double xx = MAX(m_xmolSolventMIN, xmolSolvent);
double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
for (int k = 1; k < m_kk; k++) {
m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
}
m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
/*
* Now do the pH Scaling
*/
@ -1757,111 +1832,135 @@ namespace Cantera {
Imax = Itmp;
}
}
/*
* Quick return
*/
if (Imax < m_maxIionicStrength) {
return;
}
m_molalitiesAreCropped = true;
int cropMethod = 1;
for (i = 1; i < (m_kk - 1); i++) {
double charge_i = m_speciesCharge[i];
double abs_charge_i = fabs(charge_i);
if (charge_i == 0.0) {
continue;
if (cropMethod == 0) {
/*
* Quick return
*/
if (Imax < m_maxIionicStrength) {
return;
}
for (j = (i+1); j < m_kk; j++) {
double charge_j = m_speciesCharge[j];
double abs_charge_j = fabs(charge_j);
/*
* Find the counterIJ for the symmetric binary interaction
*/
//n = m_kk*i + j;
//counterIJ = m_CounterIJ[n];
/*
* Only loop over oppositely charge species
*/
if (charge_i * charge_j < 0) {
Iac_max = m_maxIionicStrength;
if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) {
Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
if (Imax > Iac_max) {
m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
}
Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
if (Imax > Iac_max) {
m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
}
} else {
Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
if (Imax > Iac_max) {
m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
}
Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
if (Imax > Iac_max) {
m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
}
}
}
}
}
m_molalitiesAreCropped = true;
/*
* Do this loop 10 times until we have achieved charge neutrality
* in the cropped molalities
*/
for (int times = 0; times< 10; times++) {
double anion_charge = 0.0;
double cation_charge = 0.0;
int anion_contrib_max_i = -1;
double anion_contrib_max = -1.0;
int cation_contrib_max_i = -1;
double cation_contrib_max = -1.0;
for (i = 0; i < m_kk; i++) {
for (i = 1; i < (m_kk - 1); i++) {
double charge_i = m_speciesCharge[i];
if (charge_i < 0.0) {
double anion_contrib = - m_molalitiesCropped[i] * charge_i;
anion_charge += anion_contrib ;
if (anion_contrib > anion_contrib_max) {
anion_contrib_max = anion_contrib;
anion_contrib_max_i = i;
}
} else if (charge_i > 0.0) {
double cation_contrib = m_molalitiesCropped[i] * charge_i;
cation_charge += cation_contrib ;
if (cation_contrib > cation_contrib_max) {
cation_contrib_max = cation_contrib;
cation_contrib_max_i = i;
double abs_charge_i = fabs(charge_i);
if (charge_i == 0.0) {
continue;
}
for (j = (i+1); j < m_kk; j++) {
double charge_j = m_speciesCharge[j];
double abs_charge_j = fabs(charge_j);
/*
* Find the counterIJ for the symmetric binary interaction
*/
//n = m_kk*i + j;
//counterIJ = m_CounterIJ[n];
/*
* Only loop over oppositely charge species
*/
if (charge_i * charge_j < 0) {
Iac_max = m_maxIionicStrength;
if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) {
Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
if (Imax > Iac_max) {
m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
}
Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
if (Imax > Iac_max) {
m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
}
} else {
Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
if (Imax > Iac_max) {
m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
}
Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
if (Imax > Iac_max) {
m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
}
}
}
}
}
double total_charge = cation_charge - anion_charge;
if (total_charge > 1.0E-8) {
double desiredCrop = total_charge/m_speciesCharge[cation_contrib_max_i];
double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
if (desiredCrop < maxCrop) {
m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
break;
} else {
m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
/*
* Do this loop 10 times until we have achieved charge neutrality
* in the cropped molalities
*/
for (int times = 0; times< 10; times++) {
double anion_charge = 0.0;
double cation_charge = 0.0;
int anion_contrib_max_i = -1;
double anion_contrib_max = -1.0;
int cation_contrib_max_i = -1;
double cation_contrib_max = -1.0;
for (i = 0; i < m_kk; i++) {
double charge_i = m_speciesCharge[i];
if (charge_i < 0.0) {
double anion_contrib = - m_molalitiesCropped[i] * charge_i;
anion_charge += anion_contrib ;
if (anion_contrib > anion_contrib_max) {
anion_contrib_max = anion_contrib;
anion_contrib_max_i = i;
}
} else if (charge_i > 0.0) {
double cation_contrib = m_molalitiesCropped[i] * charge_i;
cation_charge += cation_contrib ;
if (cation_contrib > cation_contrib_max) {
cation_contrib_max = cation_contrib;
cation_contrib_max_i = i;
}
}
}
} else if (total_charge < -1.0E-8) {
double desiredCrop = total_charge/m_speciesCharge[anion_contrib_max_i];
double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
if (desiredCrop < maxCrop) {
m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
break;
double total_charge = cation_charge - anion_charge;
if (total_charge > 1.0E-8) {
double desiredCrop = total_charge/m_speciesCharge[cation_contrib_max_i];
double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
if (desiredCrop < maxCrop) {
m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
break;
} else {
m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
}
} else if (total_charge < -1.0E-8) {
double desiredCrop = total_charge/m_speciesCharge[anion_contrib_max_i];
double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
if (desiredCrop < maxCrop) {
m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
break;
} else {
m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
}
} else {
m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
break;
}
} else {
break;
}
}
if (cropMethod == 1) {
double *molF = DATA_PTR(m_gamma_tmp);
getMoleFractions(molF);
double xmolSolvent = molF[m_indexSolvent];
if (xmolSolvent >= MC_X_o_cutoff_) {
return;
}
double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
double denomInv = 1.0/ (m_Mnaught * p);
for (int k = 0; k < m_kk; k++) {
m_molalitiesCropped[k] = molF[k] * denomInv ;
}
}
}
/*
@ -2247,7 +2346,7 @@ namespace Cantera {
* molalitysum is the sum of the molalities over all solutes,
* even those with zero charge.
*/
double molalitysum = 0.0;
double molalitysumUncropped = 0.0;
double *gfunc = DATA_PTR(m_gfunc_IJ);
double *g2func = DATA_PTR(m_g2func_IJ);
@ -2261,7 +2360,7 @@ namespace Cantera {
double *Phiphi = DATA_PTR(m_PhiPhi_IJ);
double *CMX = DATA_PTR(m_CMX_IJ);
double x1, x2;
double Aphi, F, zsqF;
double sum1, sum2, sum3, sum4, sum5, term1;
@ -2288,7 +2387,7 @@ namespace Cantera {
Is += charge[n] * charge[n] * molality[n];
// total molar charge
molarcharge += fabs(charge[n]) * molality[n];
molalitysum += molality[n];
molalitysumUncropped += m_molalities[n];
}
Is *= 0.5;
if (Is > m_maxIionicStrength) {
@ -3178,8 +3277,8 @@ namespace Cantera {
* Calculate the osmotic coefficient from
* osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
*/
if (molalitysum > 1.0E-150) {
osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysum);
if (molalitysumUncropped > 1.0E-150) {
osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
} else {
osmotic_coef = 1.0;
}
@ -3201,7 +3300,7 @@ namespace Cantera {
printf(" Step 10: \n");
}
#endif
lnwateract = -(m_weightSolvent/1000.0) * molalitysum * osmotic_coef;
lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
wateract = exp(lnwateract);
/*
@ -3213,11 +3312,12 @@ namespace Cantera {
* ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
*/
double xmolSolvent = moleFraction(m_indexSolvent);
m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xmolSolvent);
double xx = MAX(m_xmolSolventMIN, xmolSolvent);
m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
#ifdef DEBUG_MODE
if (m_debugCalc) {
printf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
printf(" molalitySum = %16.7g\n", molalitysum);
printf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
printf(" ln_a_water=%10.6f a_water=%10.6f\n\n",
lnwateract, wateract);
}

30
Cantera/src/thermo/HMWSoln.h
View file

@ -3165,6 +3165,34 @@ namespace Cantera {
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_bgCut_;
//! value of the solvent mole fraction that centers the cutoff polynomials
//! for the cutoff =1 process;
doublereal MC_X_o_cutoff_;
//! gamma_o value for the cutoff process at the zero solvent point
doublereal MC_X_o_min_;
//! Parameter in the Molality Exp cutoff treatment
/*!
* This is the slope of the p function at the zero solvent point
* Default value is 0.0
*/
doublereal MC_slopepCut_;
//! Parameter in the Molality Exp cutoff treatment
doublereal MC_dpCut_;
//! Parameter in the Molality Exp cutoff treatment
doublereal MC_epCut_;
//! Parameter in the Molality Exp cutoff treatment
doublereal MC_apCut_;
//! Parameter in the Molality Exp cutoff treatment
doublereal MC_bpCut_;
//! Parameter in the Molality Exp cutoff treatment
doublereal MC_cpCut_;
//! Local error routine
/*!
@ -3416,7 +3444,7 @@ namespace Cantera {
//! Precalculate the IMS Cutoff parameters for typeCutoff = 2
void calcIMSCutoffParams_();
void calcMCCutoffParams_();
//! Utility function to assign an integer value from a string
//! for the ElectrolyteSpeciesType field.
/*!

52
Cantera/src/thermo/HMWSoln_input.cpp
View file

@ -1615,7 +1615,14 @@ namespace Cantera {
}
}
VPStandardStateTP::initThermoXML(phaseNode, id);
IMS_typeCutoff_ = 2;
if (IMS_typeCutoff_ == 2) {
calcIMSCutoffParams_();
}
calcMCCutoffParams_();
setMoleFSolventMin(1.0E-5);
MolalityVPSSTP::initThermoXML(phaseNode, id);
/*
* Lastly set the state
*/
@ -1648,7 +1655,7 @@ namespace Cantera {
}
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
" failed to converge on the f polynomial");
}
converged = false;
@ -1671,8 +1678,45 @@ namespace Cantera {
}
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
throw CanteraError("HMWSoln::calcIMSCutoffParams_()",
" failed to converge on the g polynomial");
}
}
// Precalculate the MC Cutoff parameters
void HMWSoln::calcMCCutoffParams_() {
MC_X_o_min_ = 0.3;
MC_X_o_cutoff_ = 0.6;
MC_slopepCut_ = 0.05;
MC_cpCut_ = 0.25;
// Initial starting values
MC_apCut_ = MC_X_o_min_;
MC_epCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
double damp = 0.5;
for (its = 0; its < 500 && !converged; its++) {
oldV = MC_epCut_;
MC_apCut_ = damp *(MC_X_o_min_ - MC_epCut_) + (1-damp) * MC_apCut_;
double MC_bpCutNew = MC_apCut_ / MC_cpCut_ + MC_slopepCut_ - 1.0;
MC_bpCut_ = damp * MC_bpCutNew + (1-damp) * MC_bpCut_;
double MC_dpCutNew = ((- MC_apCut_/MC_cpCut_ + MC_bpCut_ - MC_bpCut_ * MC_X_o_cutoff_/MC_cpCut_)
/
(MC_X_o_cutoff_ * MC_X_o_cutoff_/MC_cpCut_ - 2.0 * MC_X_o_cutoff_));
MC_dpCut_ = damp * MC_dpCutNew + (1-damp) * MC_dpCut_;
double tmp = MC_apCut_ + MC_X_o_cutoff_*( MC_bpCut_ + MC_dpCut_ * MC_X_o_cutoff_);
double eterm = std::exp(- MC_X_o_cutoff_ / MC_cpCut_);
MC_epCut_ = - eterm * (tmp);
double diff = MC_epCut_ - oldV;
if (fabs(diff) < 1.0E-14) {
converged = true;
}
}
if (!converged) {
throw CanteraError("HMWSoln::calcMCCutoffParams_()",
" failed to converge on the p polynomial");
}
}

8
Cantera/src/thermo/IdealMolalSoln.cpp
View file

@ -1,5 +1,4 @@
/**
*
* @file IdealMolalSoln.cpp
* ThermoPhase object for the ideal molal equation of
* state (see \ref thermoprops
@ -19,7 +18,7 @@
* U.S. Government retains certain rights in this software.
*/
/*
* $Author$
*
* $Date$
* $Revision$
*/
@ -1091,6 +1090,9 @@ namespace Cantera {
}
MolalityVPSSTP::initThermoXML(phaseNode, id);
setMoleFSolventMin(1.0E-5);
/*
* Set the state
*/
@ -1340,7 +1342,7 @@ namespace Cantera {
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
" failed to converge on the g polynomial");
}
}