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

Started transfering over some of the cropping capabilities from

Browse source 
IdealMolalSoln to HNWSoln.
This commit is contained in:
Harry Moffat 2008-11-07 01:52:46 +00:00
parent 022f3da1da
commit dda3f02e71
5 changed files with 497 additions and 187 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -51,6 +51,21 @@ namespace Cantera {
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.5),
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),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -81,6 +96,21 @@ namespace Cantera {
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.5),
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),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -105,6 +135,21 @@ namespace Cantera {
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.5),
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),
m_debugCalc(0)
{
for (int i = 0; i < 17; i++) {
@ -135,6 +180,21 @@ namespace Cantera {
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.5),
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),
m_debugCalc(0)
{
/*
@ -257,6 +317,23 @@ namespace Cantera {
m_CMX_IJ_P = b.m_CMX_IJ_P;
m_gamma = b.m_gamma;
IMS_lnActCoeffMolal_ = b.IMS_lnActCoeffMolal_;
IMS_typeCutoff_ = b.IMS_typeCutoff_;
IMS_X_o_cutoff_ = b.IMS_X_o_cutoff_;
IMS_gamma_o_min_ = b.IMS_gamma_o_min_;
IMS_gamma_k_min_ = b.IMS_gamma_k_min_;
IMS_cCut_ = b.IMS_cCut_;
IMS_slopefCut_ = b.IMS_slopefCut_;
IMS_dfCut_ = b.IMS_dfCut_;
IMS_efCut_ = b.IMS_efCut_;
IMS_afCut_ = b.IMS_afCut_;
IMS_bfCut_ = b.IMS_bfCut_;
IMS_slopegCut_ = b.IMS_slopegCut_;
IMS_dgCut_ = b.IMS_dgCut_;
IMS_egCut_ = b.IMS_egCut_;
IMS_agCut_ = b.IMS_agCut_;
IMS_bgCut_ = b.IMS_bgCut_;
m_CounterIJ = b.m_CounterIJ;
m_molalitiesCropped = b.m_molalitiesCropped;
m_molalitiesAreCropped= b.m_molalitiesAreCropped;
@ -1561,6 +1638,7 @@ namespace Cantera {
m_gamma.resize(leng, 0.0);
IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
counterIJ_setup();
}
@ -1610,6 +1688,11 @@ namespace Cantera {
*/
s_updatePitzerCoeffWRTemp();
/*
* Calculate the IMS cutoff factors
*/
s_updateIMS_lnMolalityActCoeff();
/*
* Now do the main calculation.
*/
@ -5297,6 +5380,128 @@ namespace Cantera {
*etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
}
}
// This function will be called to update the internally storred
// natural logarithm of the molality activity coefficients
/*
* Normally they are all one. However, sometimes they are not,
* due to stability schemes
*
* gamma_k_molar = gamma_k_molal / Xmol_solvent
*
* gamma_o_molar = gamma_o_molal
*/
void HMWSoln::s_updateIMS_lnMolalityActCoeff() const {
int k;
double tmp;
/*
* Calculate the molalities. Currently, the molalities
* may not be current with respect to the contents of the
* State objects' data.
*/
calcMolalities();
double xmolSolvent = moleFraction(m_indexSolvent);
double xx = MAX(m_xmolSolventMIN, xmolSolvent);
if (IMS_typeCutoff_ == 0) {
for (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0 ) {
for (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
tmp = log(xx * IMS_gamma_k_min_);
for (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = 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_;
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 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 h1_g = h1 / IMS_gamma_o_min_;
double h1_g_prime = h1_prime / IMS_gamma_o_min_;
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;
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 (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
}
}
// Exponentials - trial 2
else if (IMS_typeCutoff_ == 2) {
if (xmolSolvent > IMS_X_o_cutoff_) {
for (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= 0.0;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = - 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));
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 (k = 1; k < m_kk; k++) {
IMS_lnActCoeffMolal_[k]= tmp;
}
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
}
}
return;
}
/**
* This routine prints out the input pitzer coefficients for the

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

@ -2196,7 +2196,7 @@ namespace Cantera {
//@}
protected:
private:
/**
* This is the form of the Pitzer parameterization
@ -2294,7 +2294,6 @@ namespace Cantera {
*/
double m_TempPitzerRef;
protected:
/**
* Stoichiometric ionic strength on the molality scale.
* This differs from m_IionicMolality in the sense that
@ -2322,7 +2321,7 @@ namespace Cantera {
*/
int m_form_A_Debye;
protected:
private:
/**
* A_Debye -> this expression appears on the top of the
* ln actCoeff term in the general Debye-Huckel
@ -2356,13 +2355,11 @@ namespace Cantera {
*/
mutable double m_A_Debye;
//! Water standard state calculator
/*!
* derived from the equation of state for water.
*/
PDSS *m_waterSS;
//PDSS *m_waterSS;
//! density of standard-state water
/*!
@ -2719,11 +2716,11 @@ namespace Cantera {
* -------- Temporary Variables Used in the Activity Coeff Calc
*/
//! Cropped values of the molalities used in activity
//! Cropped and modified values of the molalities used in activity
//! coefficient calculations
mutable vector_fp m_molalitiesCropped;
//! Boolean indicating whether the molalities are cropped
//! Boolean indicating whether the molalities are cropped or are modified
mutable bool m_molalitiesAreCropped;
//! a counter variable for keeping track of symmetric binary
@ -2934,7 +2931,66 @@ namespace Cantera {
*/
mutable vector_fp m_gamma;
private:
//! Logarithm of the molal activity coefficients
/*!
* Normally these are all 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 minimun 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 having to do with rate of exp decay
doublereal IMS_dfCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_efCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_afCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_bfCut_;
//! Parameter in the polyExp cutoff treatment
/*!
* This is the slope of the g function at the zero solvent point
* Default value is 0.0
*/
doublereal IMS_slopegCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_dgCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_egCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_agCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal IMS_bgCut_;
//! Local error routine
/*!
@ -2951,6 +3007,18 @@ namespace Cantera {
*/
void s_update_lnMolalityActCoeff() const;
//! This function will be called to update the internally storred
//! natural logarithm of the molality activity coefficients
/*
* Normally they are all one. However, sometimes they are not,
* due to stability schemes
*
* gamma_k_molar = gamma_k_molal / Xmol_solvent
*
* gamma_o_molar = gamma_o_molal
*/
void s_updateIMS_lnMolalityActCoeff() const;
public:
//! Calculates the temperature derivative of the
@ -3132,7 +3200,11 @@ namespace Cantera {
*/
void readXMLLambdaNeutral(XML_Node &BinSalt);
//! utility function to assign an integer value from a string
//! Precalculate the IMS Cutoff parameters for typeCutoff = 2
void HMWSoln::calcIMSCutoffParams_();
//! Utility function to assign an integer value from a string
//! for the ElectrolyteSpeciesType field.
/*!
* @param estString string name of the electrolyte species type

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

@ -1303,4 +1303,55 @@ namespace Cantera {
//}
}
// Precalculate the IMS Cutoff parameters for typeCutoff = 2
void HMWSoln::calcIMSCutoffParams_() {
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
IMS_efCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
for (its = 0; its < 100 && !converged; its++) {
oldV = IMS_efCut_;
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) -IMS_efCut_;
IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
/
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
double tmp = IMS_afCut_ + IMS_X_o_cutoff_*( IMS_bfCut_ + IMS_dfCut_ *IMS_X_o_cutoff_);
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
IMS_efCut_ = - eterm * (tmp);
if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
converged = true;
}
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
}
converged = false;
double f_0 = IMS_afCut_ + IMS_efCut_;
double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
IMS_egCut_ = 0.0;
for (its = 0; its < 100 && !converged; its++) {
oldV = IMS_egCut_;
double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
IMS_agCut_ = exp(lng_0) - IMS_egCut_;
IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
/
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
double tmp = IMS_agCut_ + IMS_X_o_cutoff_*( IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
IMS_egCut_ = - eterm * (tmp);
if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
converged = true;
}
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
}
}
}

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

@ -42,21 +42,21 @@ namespace Cantera {
IdealMolalSoln::IdealMolalSoln() :
MolalityVPSSTP(),
m_formGC(2),
typeCutoff_(0),
X_o_cutoff_(0.20),
gamma_o_min_(0.00001),
gamma_k_min_(10.0),
cCut_(.05),
slopefCut_(0.6),
dfCut_(0.0),
efCut_(0.0),
afCut_(0.0),
bfCut_(0.0),
slopegCut_(0.0),
dgCut_(0.0),
egCut_(0.0),
agCut_(0.0),
bgCut_(0.0)
IMS_typeCutoff_(0),
IMS_X_o_cutoff_(0.20),
IMS_gamma_o_min_(0.00001),
IMS_gamma_k_min_(10.0),
IMS_cCut_(.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)
{
}
@ -88,26 +88,26 @@ namespace Cantera {
MolalityVPSSTP::operator=(b);
m_speciesMolarVolume = b.m_speciesMolarVolume;
m_formGC = b.m_formGC;
typeCutoff_ = b.typeCutoff_;
X_o_cutoff_ = b.X_o_cutoff_;
gamma_o_min_ = b.gamma_o_min_;
gamma_k_min_ = b.gamma_k_min_;
cCut_ = b.cCut_;
slopefCut_ = b.slopefCut_;
dfCut_ = b.dfCut_;
efCut_ = b.efCut_;
afCut_ = b.afCut_;
bfCut_ = b.bfCut_;
slopegCut_ = b.slopegCut_;
dgCut_ = b.dgCut_;
egCut_ = b.egCut_;
agCut_ = b.agCut_;
bgCut_ = b.bgCut_;
IMS_typeCutoff_ = b.IMS_typeCutoff_;
IMS_X_o_cutoff_ = b.IMS_X_o_cutoff_;
IMS_gamma_o_min_ = b.IMS_gamma_o_min_;
IMS_gamma_k_min_ = b.IMS_gamma_k_min_;
IMS_cCut_ = b.IMS_cCut_;
IMS_slopefCut_ = b.IMS_slopefCut_;
IMS_dfCut_ = b.IMS_dfCut_;
IMS_efCut_ = b.IMS_efCut_;
IMS_afCut_ = b.IMS_afCut_;
IMS_bfCut_ = b.IMS_bfCut_;
IMS_slopegCut_ = b.IMS_slopegCut_;
IMS_dgCut_ = b.IMS_dgCut_;
IMS_egCut_ = b.IMS_egCut_;
IMS_agCut_ = b.IMS_agCut_;
IMS_bgCut_ = b.IMS_bgCut_;
m_expg0_RT = b.m_expg0_RT;
m_pe = b.m_pe;
m_pp = b.m_pp;
m_tmpV = b.m_tmpV;
m_lnActCoeffMolal = b.m_lnActCoeffMolal;
IMS_lnActCoeffMolal_ = b.IMS_lnActCoeffMolal_;
}
return *this;
}
@ -115,21 +115,21 @@ namespace Cantera {
IdealMolalSoln::IdealMolalSoln(std::string inputFile, std::string id) :
MolalityVPSSTP(),
m_formGC(2),
typeCutoff_(0),
X_o_cutoff_(0.2),
gamma_o_min_(0.00001),
gamma_k_min_(10.0),
cCut_(.05),
slopefCut_(0.6),
dfCut_(0.0),
efCut_(0.0),
afCut_(0.0),
bfCut_(0.0),
slopegCut_(0.0),
dgCut_(0.0),
egCut_(0.0),
agCut_(0.0),
bgCut_(0.0)
IMS_typeCutoff_(0),
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(0.00001),
IMS_gamma_k_min_(10.0),
IMS_cCut_(.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)
{
constructPhaseFile(inputFile, id);
}
@ -137,21 +137,21 @@ namespace Cantera {
IdealMolalSoln::IdealMolalSoln(XML_Node& root, std::string id) :
MolalityVPSSTP(),
m_formGC(2),
typeCutoff_(0),
X_o_cutoff_(0.2),
gamma_o_min_(0.00001),
gamma_k_min_(10.0),
cCut_(.05),
slopefCut_(0.6),
dfCut_(0.0),
efCut_(0.0),
afCut_(0.0),
bfCut_(0.0),
slopegCut_(0.0),
dgCut_(0.0),
egCut_(0.0),
agCut_(0.0),
bgCut_(0.0)
IMS_typeCutoff_(0),
IMS_X_o_cutoff_(0.2),
IMS_gamma_o_min_(0.00001),
IMS_gamma_k_min_(10.0),
IMS_cCut_(.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)
{
constructPhaseXML(root, id);
}
@ -532,7 +532,7 @@ namespace Cantera {
* Update the molality array, m_molalities()
* This requires an update due to mole fractions
*/
if (typeCutoff_ == 0) {
if (IMS_typeCutoff_ == 0) {
calcMolalities();
for (int k = 0; k < m_kk; k++) {
ac[k] = m_molalities[k];
@ -548,11 +548,11 @@ namespace Cantera {
* Now calculate the array of activities.
*/
for (int k = 1; k < m_kk; k++) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
}
double xmolSolvent = moleFraction(m_indexSolvent);
ac[m_indexSolvent] =
exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent;
exp(IMS_lnActCoeffMolal_[m_indexSolvent]) * xmolSolvent;
}
}
@ -570,7 +570,7 @@ namespace Cantera {
*/
void IdealMolalSoln::
getMolalityActivityCoefficients(doublereal* acMolality) const {
if (typeCutoff_ == 0) {
if (IMS_typeCutoff_ == 0) {
for (int k = 0; k < m_kk; k++) {
acMolality[k] = 1.0;
}
@ -580,7 +580,7 @@ namespace Cantera {
exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
} else {
s_updateIMS_lnMolalityActCoeff();
std::copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality);
std::copy(IMS_lnActCoeffMolal_.begin(), IMS_lnActCoeffMolal_.end(), acMolality);
for (int k = 0; k < m_kk; k++) {
acMolality[k] = exp(acMolality[k]);
}
@ -640,7 +640,7 @@ namespace Cantera {
*/
doublereal RT = GasConstant * temperature();
if (typeCutoff_ == 0 || xmolSolvent > 3.* X_o_cutoff_/2.0) {
if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
for (int k = 1; k < m_kk; k++) {
xx = fmaxx(m_molalities[k], xxSmall);
@ -664,11 +664,11 @@ namespace Cantera {
for (int k = 1; k < m_kk; k++) {
xx = MAX(m_molalities[k], xxSmall);
mu[k] += RT * (log(xx) + m_lnActCoeffMolal[k]);
mu[k] += RT * (log(xx) + IMS_lnActCoeffMolal_[k]);
}
xx = MAX(xmolSolvent, xxSmall);
mu[m_indexSolvent] +=
RT * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]);
RT * (log(xx) + IMS_lnActCoeffMolal_[m_indexSolvent]);
}
}
@ -718,7 +718,7 @@ namespace Cantera {
doublereal R = GasConstant;
doublereal mm;
calcMolalities();
if (typeCutoff_ == 0) {
if (IMS_typeCutoff_ == 0) {
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
mm = fmaxx(SmallNumber, m_molalities[k]);
@ -741,12 +741,12 @@ namespace Cantera {
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
mm = fmaxx(SmallNumber, m_molalities[k]);
sbar[k] -= R * (log(mm) + m_lnActCoeffMolal[k]);
sbar[k] -= R * (log(mm) + IMS_lnActCoeffMolal_[k]);
}
}
double xmolSolvent = moleFraction(m_indexSolvent);
mm = fmaxx(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
sbar[m_indexSolvent] -= R *(log(mm) + IMS_lnActCoeffMolal_[m_indexSolvent]);
}
}
@ -1003,7 +1003,7 @@ namespace Cantera {
if (thermoNode.hasChild("activityCoefficients")) {
XML_Node& acNode = thermoNode.child("activityCoefficients");
std::string modelString = acNode.attrib("model");
typeCutoff_ = 0;
IMS_typeCutoff_ = 0;
if (modelString != "IdealMolalSoln") {
throw CanteraError("IdealMolalSoln::initThermoXML",
"unknown ActivityCoefficient model: " + modelString);
@ -1013,31 +1013,31 @@ namespace Cantera {
modelString = ccNode.attrib("model");
if (modelString != "") {
if (modelString == "polyExp") {
typeCutoff_ = 2;
IMS_typeCutoff_ = 2;
} else if (modelString == "poly") {
typeCutoff_ = 1;
IMS_typeCutoff_ = 1;
} else {
throw CanteraError("IdealMolalSoln::initThermoXML",
"Unknown idealMolalSolnCutoff form: " + modelString);
}
if (ccNode.hasChild("gamma_o_limit")) {
gamma_o_min_ = getFloat(ccNode, "gamma_o_limit");
IMS_gamma_o_min_ = getFloat(ccNode, "gamma_o_limit");
}
if (ccNode.hasChild("gamma_k_limit")) {
gamma_k_min_ = getFloat(ccNode, "gamma_k_limit");
IMS_gamma_k_min_ = getFloat(ccNode, "gamma_k_limit");
}
if (ccNode.hasChild("X_o_cutoff")) {
X_o_cutoff_ = getFloat(ccNode, "X_o_cutoff");
IMS_X_o_cutoff_ = getFloat(ccNode, "X_o_cutoff");
}
if (ccNode.hasChild("c_0_param")) {
cCut_ = getFloat(ccNode, "c_0_param");
IMS_cCut_ = getFloat(ccNode, "c_0_param");
}
if (ccNode.hasChild("slope_f_limit")) {
slopefCut_ = getFloat(ccNode, "slope_f_limit");
IMS_slopefCut_ = getFloat(ccNode, "slope_f_limit");
}
if (ccNode.hasChild("slope_g_limit")) {
slopegCut_ = getFloat(ccNode, "slope_g_limit");
IMS_slopegCut_ = getFloat(ccNode, "slope_g_limit");
}
}
@ -1083,8 +1083,8 @@ namespace Cantera {
m_speciesMolarVolume[k] = getFloat(*ss, "molarVolume", "toSI");
}
typeCutoff_ = 2;
if (typeCutoff_ == 2) {
IMS_typeCutoff_ = 2;
if (IMS_typeCutoff_ == 2) {
calcIMSCutoffParams_();
}
@ -1173,47 +1173,47 @@ namespace Cantera {
double xmolSolvent = moleFraction(m_indexSolvent);
double xx = MAX(m_xmolSolventMIN, xmolSolvent);
if (typeCutoff_ == 0) {
if (IMS_typeCutoff_ == 0) {
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= 0.0;
IMS_lnActCoeffMolal_[k]= 0.0;
}
m_lnActCoeffMolal[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * X_o_cutoff_/2.0 ) {
} else if (IMS_typeCutoff_ == 1) {
if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0 ) {
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= 0.0;
IMS_lnActCoeffMolal_[k]= 0.0;
}
m_lnActCoeffMolal[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else if (xmolSolvent < X_o_cutoff_/2.0) {
tmp = log(xx * gamma_k_min_);
} else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
tmp = log(xx * IMS_gamma_k_min_);
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= tmp;
IMS_lnActCoeffMolal_[k]= tmp;
}
m_lnActCoeffMolal[m_indexSolvent] = log(gamma_o_min_);
IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
return;
} else {
/*
* If we are in the middle region, calculate the connecting polynomials
*/
double xminus = xmolSolvent - X_o_cutoff_/2.0;
double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
double xminus2 = xminus * xminus;
double xminus3 = xminus2 * xminus;
double x_o_cut2 = X_o_cutoff_ * X_o_cutoff_;
double x_o_cut3 = x_o_cut2 * X_o_cutoff_;
double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
double h2 = 3.5 * xminus2 / X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
double h2_prime = 7.0 * xminus / X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
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 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 h1_g = h1 / gamma_o_min_;
double h1_g_prime = h1_prime / gamma_o_min_;
double h1_g = h1 / IMS_gamma_o_min_;
double h1_g_prime = h1_prime / IMS_gamma_o_min_;
double alpha = 1.0 / ( exp(1.0) * gamma_k_min_);
double alpha = 1.0 / ( exp(1.0) * IMS_gamma_k_min_);
double h1_f = h1 * alpha;
double h1_f_prime = h1_prime * alpha;
@ -1229,34 +1229,34 @@ namespace Cantera {
tmp = log(xmolSolvent) + lngammak;
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= tmp;
IMS_lnActCoeffMolal_[k]= tmp;
}
m_lnActCoeffMolal[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
}
}
// Exponentials - trial 2
else if (typeCutoff_ == 2) {
if (xmolSolvent > X_o_cutoff_) {
else if (IMS_typeCutoff_ == 2) {
if (xmolSolvent > IMS_X_o_cutoff_) {
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= 0.0;
IMS_lnActCoeffMolal_[k]= 0.0;
}
m_lnActCoeffMolal[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
return;
} else {
double xoverc = xmolSolvent/cCut_;
double xoverc = xmolSolvent/IMS_cCut_;
double eterm = std::exp(-xoverc);
double fptmp = bfCut_ - afCut_ / cCut_ - bfCut_*xoverc
+ 2.0*dfCut_*xmolSolvent - dfCut_*xmolSolvent*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 + efCut_ + eterm * (afCut_ + xmolSolvent * (bfCut_ + dfCut_*xmolSolvent));
double f = xmolSolvent + IMS_efCut_ + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
double gptmp = bgCut_ - agCut_ / cCut_ - bgCut_*xoverc
+ 2.0*dgCut_*xmolSolvent - dgCut_*xmolSolvent*xoverc;
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 + egCut_ + eterm * (agCut_ + xmolSolvent * (bgCut_ + dgCut_*xmolSolvent));
double g = xmolSolvent + IMS_egCut_ + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
@ -1264,9 +1264,9 @@ namespace Cantera {
tmp = log(xx) + lngammak;
for (k = 1; k < m_kk; k++) {
m_lnActCoeffMolal[k]= tmp;
IMS_lnActCoeffMolal_[k]= tmp;
}
m_lnActCoeffMolal[m_indexSolvent] = lngammao;
IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
}
}
return;
@ -1289,67 +1289,50 @@ namespace Cantera {
m_pp.resize(leng);
m_speciesMolarVolume.resize(leng);
m_tmpV.resize(leng);
m_lnActCoeffMolal.resize(leng);
IMS_lnActCoeffMolal_.resize(leng);
}
void IdealMolalSoln::calcIMSCutoffParams_() {
afCut_ = 1.0 / (std::exp(1.0) * gamma_k_min_);
efCut_ = 0.0;
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
IMS_efCut_ = 0.0;
bool converged = false;
double oldV = 0.0;
int its;
for (its = 0; its < 100 && !converged; its++) {
oldV = efCut_;
afCut_ = 1.0 / (std::exp(1.0) * gamma_k_min_) -efCut_;
bfCut_ = afCut_ / cCut_ + slopefCut_ - 1.0;
dfCut_ = ((- afCut_/cCut_ + bfCut_ - bfCut_*X_o_cutoff_/cCut_)
/
(X_o_cutoff_*X_o_cutoff_/cCut_ - 2.0 * X_o_cutoff_));
double tmp = afCut_ + X_o_cutoff_*( bfCut_ + dfCut_ *X_o_cutoff_);
double eterm = std::exp(-X_o_cutoff_/cCut_);
efCut_ = - eterm * (tmp);
if (fabs(efCut_ - oldV) < 1.0E-14) {
oldV = IMS_efCut_;
IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) - IMS_efCut_;
IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
/
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
double tmp = IMS_afCut_ + IMS_X_o_cutoff_*( IMS_bfCut_ + IMS_dfCut_ * IMS_X_o_cutoff_);
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
IMS_efCut_ = - eterm * (tmp);
if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
converged = true;
}
}
if (!converged) {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
}
converged = false;
double f_0 = afCut_ + efCut_;
double f_prime_0 = 1.0 - afCut_ / cCut_ + bfCut_;
egCut_ = 0.0;
double f_0 = IMS_afCut_ + IMS_efCut_;
double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
IMS_egCut_ = 0.0;
for (its = 0; its < 100 && !converged; its++) {
oldV = egCut_;
double lng_0 = -log(gamma_o_min_) - f_prime_0 / f_0;
agCut_ = exp(lng_0) - egCut_;
bgCut_ = agCut_ / cCut_ + slopegCut_ - 1.0;
dgCut_ = ((- agCut_/cCut_ + bgCut_ - bgCut_*X_o_cutoff_/cCut_)
/
(X_o_cutoff_*X_o_cutoff_/cCut_ - 2.0 * X_o_cutoff_));
double tmp = agCut_ + X_o_cutoff_*( bgCut_ + dgCut_ *X_o_cutoff_);
double eterm = std::exp(-X_o_cutoff_/cCut_);
egCut_ = - eterm * (tmp);
if (fabs(egCut_ - oldV) < 1.0E-14) {
oldV = IMS_egCut_;
double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
IMS_agCut_ = exp(lng_0) - IMS_egCut_;
IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
/
(IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
double tmp = IMS_agCut_ + IMS_X_o_cutoff_*( IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
IMS_egCut_ = - eterm * (tmp);
if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
converged = true;
}
}
@ -1357,7 +1340,6 @@ namespace Cantera {
throw CanteraError(" IdealMolalSoln::calcCutoffParams_()",
" failed to converge on the f polynomial");
}
}
}

32
Cantera/src/thermo/IdealMolalSoln.h
View file

@ -902,7 +902,7 @@ namespace Cantera {
public:
//! Cutoff type
int typeCutoff_;
int IMS_typeCutoff_;
private:
@ -931,58 +931,58 @@ namespace Cantera {
/*!
* Normally these are all one. However, stability schemes will change that
*/
mutable vector_fp m_lnActCoeffMolal;
mutable vector_fp IMS_lnActCoeffMolal_;
public:
//! value of the solute mole fraction that centers the cutoff polynomials
//! for the cutoff =1 process;
doublereal X_o_cutoff_;
doublereal IMS_X_o_cutoff_;
//! gamma_o value for the cutoff process at the zero solvent point
doublereal gamma_o_min_;
doublereal IMS_gamma_o_min_;
//! gamma_k minimun for the cutoff process at the zero solvent point
doublereal gamma_k_min_;
doublereal IMS_gamma_k_min_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal cCut_;
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 slopefCut_;
doublereal IMS_slopefCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal dfCut_;
doublereal IMS_dfCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal efCut_;
doublereal IMS_efCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal afCut_;
doublereal IMS_afCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal bfCut_;
doublereal IMS_bfCut_;
//! Parameter in the polyExp cutoff treatment
/*!
* This is the slope of the g function at the zero solvent point
* Default value is 0.0
*/
doublereal slopegCut_;
doublereal IMS_slopegCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal dgCut_;
doublereal IMS_dgCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal egCut_;
doublereal IMS_egCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal agCut_;
doublereal IMS_agCut_;
//! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
doublereal bgCut_;
doublereal IMS_bgCut_;
private: