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

[Thermo] Clean up HMWSoln variables used in cutoff calculation

Browse source
This commit is contained in:
Ray Speth 2017-08-05 13:18:35 -04:00
parent 1bd950fef5
commit d07908f9c9
3 changed files with 29 additions and 131 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

26
include/cantera/thermo/HMWSoln.h
View file

@ -2216,29 +2216,13 @@ private:
//! one. However, stability schemes will change that
mutable vector_fp IMS_lnActCoeffMolal_;
//! IMS Cutoff type
int IMS_typeCutoff_;
//! value of the solute mole fraction that centers the cutoff polynomials
//! for the cutoff =1 process;
doublereal IMS_X_o_cutoff_;
//! gamma_o value for the cutoff process at the zero solvent point
doublereal IMS_gamma_o_min_;
//! gamma_k minimum for the cutoff process at the zero solvent point
doublereal IMS_gamma_k_min_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_cCut_;
//! Parameter in the polyExp cutoff treatment
/*!
* This is the slope of the f function at the zero solvent point
* Default value is 0.6
*/
doublereal IMS_slopefCut_;
//! Parameter in the polyExp cutoff treatment
/*!
* This is the slope of the g function at the zero solvent point
@ -2262,16 +2246,6 @@ private:
//! 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_;
//! @name Parameters in the Molality Exp cutoff treatment
//! @{
doublereal MC_dpCut_;

121
src/thermo/HMWSoln.cpp
View file

@ -34,12 +34,8 @@ HMWSoln::HMWSoln() :
m_waterSS(0),
m_densWaterSS(1000.),
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_slopegCut_(0.0),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -50,8 +46,6 @@ HMWSoln::HMWSoln() :
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),
@ -80,12 +74,8 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
m_waterSS(0),
m_densWaterSS(1000.),
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_slopegCut_(0.0),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -96,8 +86,6 @@ HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
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),
@ -123,12 +111,8 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
m_waterSS(0),
m_densWaterSS(1000.),
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_slopegCut_(0.0),
IMS_dfCut_(0.0),
IMS_efCut_(0.0),
@ -139,8 +123,6 @@ HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
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),
@ -3774,97 +3756,38 @@ void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
calcMolalities();
double xmolSolvent = moleFraction(0);
double xx = std::max(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
// Exponentials - trial 2
if (xmolSolvent > IMS_X_o_cutoff_) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
double tmp = log(xx * IMS_gamma_k_min_);
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
return;
} else {
// If we are in the middle region, calculate the connecting polynomials
double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
double xminus2 = xminus * xminus;
double xminus3 = xminus2 * xminus;
double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
double eterm = std::exp(-xoverc);
double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
+ 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
double f_prime = 1.0 + eterm*fptmp;
double f = xmolSolvent + IMS_efCut_
+ eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
+ 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
double g_prime = 1.0 + eterm*gptmp;
double g = xmolSolvent + IMS_egCut_
+ eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
double h1_g = h1 / IMS_gamma_o_min_;
double h1_g_prime = h1_prime / IMS_gamma_o_min_;
double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
double h1_f = h1 * alpha;
double h1_f_prime = h1_prime * alpha;
double f = h2 + h1_f;
double f_prime = h2_prime + h1_f_prime;
double g = h2 + h1_g;
double g_prime = h2_prime + h1_g_prime;
double tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
tmp = log(xmolSolvent) + lngammak;
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[0] = lngammao;
}
} else if (IMS_typeCutoff_ == 2) {
// Exponentials - trial 2
if (xmolSolvent > IMS_X_o_cutoff_) {
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/IMS_cCut_;
double eterm = std::exp(-xoverc);
double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
+ 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
double f_prime = 1.0 + eterm*fptmp;
double f = xmolSolvent + IMS_efCut_
+ eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
+ 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
double g_prime = 1.0 + eterm*gptmp;
double g = xmolSolvent + IMS_egCut_
+ eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
tmp = log(xx) + lngammak;
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[0] = lngammao;
tmp = log(xx) + lngammak;
for (size_t k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[0] = lngammao;
}
return;
}

13
src/thermo/HMWSoln_input.cpp
View file

@ -1087,10 +1087,7 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
readXMLCroppingCoefficients(acNode);
}
IMS_typeCutoff_ = 2;
if (IMS_typeCutoff_ == 2) {
calcIMSCutoffParams_();
}
calcIMSCutoffParams_();
calcMCCutoffParams_();
setMoleFSolventMin(1.0E-5);
@ -1163,6 +1160,10 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_)
void HMWSoln::calcIMSCutoffParams_()
{
double IMS_gamma_o_min_ = 1.0E-5; // value at the zero solvent point
double IMS_gamma_k_min_ = 10.0; // minimum at the zero solvent point
double IMS_slopefCut_ = 0.6; // slope of the f function at the zero solvent point
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
IMS_efCut_ = 0.0;
bool converged = false;
@ -1212,9 +1213,9 @@ void HMWSoln::calcIMSCutoffParams_()
void HMWSoln::calcMCCutoffParams_()
{
MC_X_o_min_ = 0.35;
double MC_X_o_min_ = 0.35; // value at the zero solvent point
MC_X_o_cutoff_ = 0.6;
MC_slopepCut_ = 0.02;
double MC_slopepCut_ = 0.02; // slope of the p function at the zero solvent point
MC_cpCut_ = 0.25;
// Initial starting values