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

Iteration on unifying water treatments in the 2 objects.

Browse source
This commit is contained in:
Harry Moffat 2006-07-12 00:46:34 +00:00
parent ec04dbe060
commit 60ff5fb20e
4 changed files with 351 additions and 71 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

311
Cantera/src/thermo/DebyeHuckel.cpp
View file

@ -16,6 +16,7 @@
#include "DebyeHuckel.h"
#include "importCTML.h"
#include "WaterProps.h"
namespace Cantera {
@ -32,8 +33,10 @@ namespace Cantera {
m_maxIionicStrength(30.0),
m_useHelgesonFixedForm(false),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m
m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
m_waterProps(0)
{
m_npActCoeff.resize(3);
m_npActCoeff[0] = 0.1127;
@ -57,8 +60,10 @@ namespace Cantera {
m_maxIionicStrength(30.0),
m_useHelgesonFixedForm(false),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m
m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
m_waterProps(0)
{
m_npActCoeff.resize(3);
m_npActCoeff[0] = 0.1127;
@ -76,8 +81,10 @@ namespace Cantera {
m_maxIionicStrength(3.0),
m_useHelgesonFixedForm(false),
m_IionicMolalityStoich(0.0),
m_form_A_Debye(A_DEBYE_CONST),
m_A_Debye(1.172576), // units = sqrt(kg/gmol)
m_B_Debye(3.28640E9) // units = sqrt(kg/gmol) / m
m_B_Debye(3.28640E9), // units = sqrt(kg/gmol) / m
m_waterProps(0)
{
m_npActCoeff.resize(3);
m_npActCoeff[0] = 0.1127;
@ -121,9 +128,18 @@ namespace Cantera {
m_maxIionicStrength = b.m_maxIionicStrength;
m_useHelgesonFixedForm= b.m_useHelgesonFixedForm;
m_IionicMolalityStoich= b.m_IionicMolalityStoich;
m_form_A_Debye = b.m_form_A_Debye;
m_A_Debye = b.m_A_Debye;
m_B_Debye = b.m_B_Debye;
m_B_Dot = b.m_B_Dot;
m_npActCoeff = b.m_npActCoeff;
if (m_waterProps) {
delete m_waterProps;
m_waterProps = 0;
}
if (b.m_waterProps) {
m_waterProps = new WaterProps(*(b.m_waterProps));
}
m_expg0_RT = b.m_expg0_RT;
m_pe = b.m_pe;
m_pp = b.m_pp;
@ -411,7 +427,7 @@ namespace Cantera {
* Update the molality array, m_molalities()
* This requires an update due to mole fractions
*/
_updatelnMolalityActCoeff();
s_update_lnMolalityActCoeff();
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
@ -435,7 +451,9 @@ namespace Cantera {
*/
void DebyeHuckel::
getMolalityActivityCoefficients(doublereal* acMolality) const {
_updatelnMolalityActCoeff();
A_Debye_TP(-1.0, -1.0);
s_update_lnMolalityActCoeff();
copy(m_lnActCoeffMolal.begin(), m_lnActCoeffMolal.end(), acMolality);
for (int k = 0; k < m_kk; k++) {
acMolality[k] = exp(acMolality[k]);
@ -474,7 +492,7 @@ namespace Cantera {
* Update the activity coefficients
* This also updates the internal molality array.
*/
_updatelnMolalityActCoeff();
s_update_lnMolalityActCoeff();
/*
*
*/
@ -518,8 +536,8 @@ namespace Cantera {
* Update the activity coefficients, This also update the
* internally storred molalities.
*/
_updatelnMolalityActCoeff();
_updatedlnMolalityActCoeffdT();
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
double T = temperature();
double RTT = GasConstant * T * T;
for (int k = 0; k < m_kk; k++) {
@ -569,7 +587,7 @@ namespace Cantera {
* Update the activity coefficients, This also update the
* internally storred molalities.
*/
_updatelnMolalityActCoeff();
s_update_lnMolalityActCoeff();
doublereal R = GasConstant;
doublereal mm;
@ -594,7 +612,7 @@ namespace Cantera {
*/
double dAdT = dA_DebyedT_TP();
if (dAdT != 0.0) {
_updatedlnMolalityActCoeffdT();
s_update_dlnMolalityActCoeff_dT();
double RT = R * temperature();
for (k = 0; k < m_kk; k++) {
sbar[k] -= RT * m_dlnActCoeffMolaldT[k];
@ -654,9 +672,9 @@ namespace Cantera {
* Update the activity coefficients, This also update the
* internally storred molalities.
*/
_updatelnMolalityActCoeff();
_updatedlnMolalityActCoeffdT();
_updated2lnMolalityActCoeffdT2();
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
s_update_d2lnMolalityActCoeff_dT2();
double T = temperature();
double RT = GasConstant * T;
double RTT = RT * T;
@ -1496,10 +1514,28 @@ namespace Cantera {
*/
double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const {
double T = temperature();
double A;
if (tempArg != -1.0) {
T = tempArg;
}
return m_A_Debye;
double P = pressure();
if (presArg != -1.0) {
P = presArg;
}
switch (m_form_A_Debye) {
case A_DEBYE_CONST:
A = m_A_Debye;
break;
case A_DEBYE_WATER:
A = m_waterProps->ADebye(T, P, 0);
m_A_Debye = A;
break;
default:
printf("shouldn't be here\n");
exit(-1);
}
return A;
}
/**
@ -1509,16 +1545,32 @@ namespace Cantera {
* respect to temperature as a function of temperature
* and pressure.
*
* The default is to assume that it is equal to zero
* -> note, placeholder until a better formalism is
* put in place.
* units = A_Debye has units of sqrt(gmol kg-1).
* Temp has units of Kelvin.
*/
double DebyeHuckel::dA_DebyedT_TP(double tempArg, double presArg) const {
double T = temperature();
double T = temperature();
if (tempArg != -1.0) {
T = tempArg;
}
return 0.0;
double P = pressure();
if (presArg != -1.0) {
P = presArg;
}
double dAdT;
switch (m_form_A_Debye) {
case A_DEBYE_CONST:
dAdT = 0.0;
break;
case A_DEBYE_WATER:
dAdT = m_waterProps->ADebye(T, P, 1);
//dAdT = WaterProps::ADebye(T, P, 1);
break;
default:
printf("shouldn't be here\n");
exit(-1);
}
return dAdT;
}
/**
@ -1526,18 +1578,67 @@ namespace Cantera {
*
* Returns the 2nd derivative of the A_Debye parameter with
* respect to temperature as a function of temperature
* and pressure.
*
* The default is to assume that it is equal to zero
* -> note, placeholder until a better formalism is
* put in place.
* and pressure.
*
* units = A_Debye has units of sqrt(gmol kg-1).
* Temp has units of Kelvin.
*/
double DebyeHuckel::d2A_DebyedT2_TP(double tempArg, double presArg) const {
double T = temperature();
double T = temperature();
if (tempArg != -1.0) {
T = tempArg;
}
return 0.0;
double P = pressure();
if (presArg != -1.0) {
P = presArg;
}
double d2AdT2;
switch (m_form_A_Debye) {
case A_DEBYE_CONST:
d2AdT2 = 0.0;
break;
case A_DEBYE_WATER:
d2AdT2 = m_waterProps->ADebye(T, P, 2);
break;
default:
printf("shouldn't be here\n");
exit(-1);
}
return d2AdT2;
}
/**
* dA_DebyedP_TP() (virtual)
*
* Returns the derivative of the A_Debye parameter with
* respect to pressure, as a function of temperature
* and pressure.
*
* units = A_Debye has units of sqrt(gmol kg-1).
* Pressure has units of pascals.
*/
double DebyeHuckel::dA_DebyedP_TP(double tempArg, double presArg) const {
double T = temperature();
if (tempArg != -1.0) {
T = tempArg;
}
double P = pressure();
if (presArg != -1.0) {
P = presArg;
}
double dAdP;
switch (m_form_A_Debye) {
case A_DEBYE_CONST:
dAdP = 0.0;
break;
case A_DEBYE_WATER:
dAdP = m_waterProps->ADebye(T, P, 3);
break;
default:
printf("shouldn't be here\n");
exit(-1);
}
return dAdP;
}
/*
@ -1587,6 +1688,7 @@ namespace Cantera {
m_lnActCoeffMolal.resize(leng, 0.0);
m_dlnActCoeffMolaldT.resize(leng, 0.0);
m_d2lnActCoeffMolaldT2.resize(leng, 0.0);
m_dlnActCoeffMolaldP.resize(leng, 0.0);
m_B_Dot.resize(leng, 0.0);
m_expg0_RT.resize(leng, 0.0);
m_pe.resize(leng, 0.0);
@ -1671,7 +1773,7 @@ namespace Cantera {
}
/**
* _updatelnMolalityActCoeff():
* s_update_lnMolalityActCoeff():
*
* Using internally stored values, this function calculates
* the activity coefficients for all species.
@ -1683,7 +1785,7 @@ namespace Cantera {
* ( Note this is the main routine for implementing the
* activity coefficient formulation.)
*/
void DebyeHuckel::_updatelnMolalityActCoeff() const {
void DebyeHuckel::s_update_lnMolalityActCoeff() const {
double z_k, zs_k1, zs_k2;
/*
* Update the internally storred vector of molalities
@ -1941,7 +2043,7 @@ namespace Cantera {
}
/**
* _updatedMolalityActCoeffdT() (private, const )
* s_update_dMolalityActCoeff_dT() (private, const )
*
* Using internally stored values, this function calculates
* the temperature derivative of the logarithm of the
@ -1952,7 +2054,7 @@ namespace Cantera {
* solvent activity coefficient is on the molality
* scale. It's derivative is too.
*/
void DebyeHuckel::_updatedlnMolalityActCoeffdT() const {
void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const {
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
int k;
double dAdT = dA_DebyedT_TP();
@ -2074,7 +2176,7 @@ namespace Cantera {
}
/**
* _updated2lnMolalityActCoeffdT2() (private, const )
* s_update_d2lnMolalityActCoeff_dT2() (private, const )
*
* Using internally stored values, this function calculates
* the temperature 2nd derivative of the logarithm of the
@ -2087,17 +2189,17 @@ namespace Cantera {
* scale. It's derivatives are too.
*/
void DebyeHuckel::_updated2lnMolalityActCoeffdT2() const {
void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const {
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
int k;
double dAdT = dA_DebyedT_TP();
if (dAdT == 0.0) {
double d2AdT2 = d2A_DebyedT2_TP();
if (d2AdT2 == 0.0 && dAdT == 0.0) {
for (k = 0; k < m_kk; k++) {
m_dlnActCoeffMolaldT[k] = 0.0;
m_d2lnActCoeffMolaldT2[k] = 0.0;
}
return;
}
double d2AdT2 = d2A_DebyedT2_TP();
/*
* Calculate a safe value for the mole fraction
@ -2211,6 +2313,145 @@ namespace Cantera {
}
}
/**
* s_update_dlnMolalityActCoeff_dP() (private, const )
*
* Using internally stored values, this function calculates
* the pressure derivative of the logarithm of the
* activity coefficient
* for all species in the mechanism.
*
* We assume that the activity coefficients, molalities,
* and A_Debye are current.
*
* solvent activity coefficient is on the molality
* scale. It's derivatives are too.
*/
void DebyeHuckel::s_update_dlnMolalityActCoeff_dP() const {
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
int k, est;
double dAdP = dA_DebyedP_TP();
if (dAdP == 0.0) {
for (k = 0; k < m_kk; k++) {
m_dlnActCoeffMolaldP[k] = 0.0;
}
return;
}
/*
* Calculate a safe value for the mole fraction
* of the solvent
*/
double xmolSolvent = moleFraction(m_indexSolvent);
xmolSolvent = MAX(8.689E-3, xmolSolvent);
double sqrtI = sqrt(m_IionicMolality);
double numdAdPTmp = dAdP * sqrtI;
double denomTmp = m_B_Debye * sqrtI;
switch (m_formDH) {
case DHFORM_DILUTE_LIMIT:
for (int k = 0; k < m_kk; k++) {
m_dlnActCoeffMolaldP[k] =
m_lnActCoeffMolal[k] * dAdP / m_A_Debye;
}
break;
case DHFORM_BDOT_AK:
for (int k = 0; k < m_kk; k++) {
est = m_electrolyteSpeciesType[k];
if (est == cEST_nonpolarNeutral) {
m_lnActCoeffMolal[k] = 0.0;
} else {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp * m_Aionic[k]);
}
}
m_dlnActCoeffMolaldP[m_indexSolvent] = 0.0;
coeff = 2.0 / 3.0 * dAdP * m_Mnaught * sqrtI;
tmp = 0.0;
if (denomTmp > 0.0) {
for (int k = 0; k < m_kk; k++) {
y = denomTmp * m_Aionic[k];
yp1 = y + 1.0;
sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
z_k = m_speciesCharge[k];
tmp += m_molalities[k] * z_k * z_k * sigma / 2.0;
}
}
m_dlnActCoeffMolaldP[m_indexSolvent] += coeff * tmp;
break;
case DHFORM_BDOT_ACOMMON:
denomTmp *= m_Aionic[0];
for (int k = 0; k < m_kk; k++) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
}
if (denomTmp > 0.0) {
y = denomTmp;
yp1 = y + 1.0;
sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
2.0 /3.0 * dAdP * m_Mnaught * m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_BETAIJ:
denomTmp *= m_Aionic[0];
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp);
}
}
if (denomTmp > 0.0) {
y = denomTmp;
yp1 = y + 1.0;
sigma = 3.0 / (y * y * y) * (yp1 - 1.0/yp1 - 2.0*log(yp1));
} else {
sigma = 0.0;
}
m_dlnActCoeffMolaldP[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
case DHFORM_PITZER_BETAIJ:
denomTmp *= m_Aionic[0];
tmpLn = log(1.0 + denomTmp);
for (int k = 0; k < m_kk; k++) {
if (k != m_indexSolvent) {
z_k = m_speciesCharge[k];
m_dlnActCoeffMolaldP[k] =
- z_k * z_k * numdAdPTmp / (1.0 + denomTmp)
- 2.0 * z_k * z_k * dAdP * tmpLn
/ (m_B_Debye * m_Aionic[0]);
m_dlnActCoeffMolaldP[k] /= 3.0;
}
}
sigma = 1.0 / ( 1.0 + denomTmp);
m_dlnActCoeffMolaldP[m_indexSolvent] =
(xmolSolvent - 1.0)/xmolSolvent +
2.0 /3.0 * dAdP * m_Mnaught *
m_IionicMolality * sqrtI * sigma;
break;
default:
printf("ERROR\n");
exit(-1);
break;
}
}
}

71
Cantera/src/thermo/DebyeHuckel.h
View file

@ -90,6 +90,14 @@ namespace Cantera {
#define DHFORM_BETAIJ 3
#define DHFORM_PITZER_BETAIJ 4
/*
* Acceptable ways to calculate the value of A_Debye
*/
#define A_DEBYE_CONST 0
#define A_DEBYE_WATER 1
class WaterProps;
class WaterPDSS;
/**
* Definition of the DebyeHuckel object
@ -860,12 +868,42 @@ namespace Cantera {
virtual double A_Debye_TP(double temperature = -1.0,
double pressure = -1.0) const;
/**
* Value of the derivative of the Debye Huckel constant with
* respect to temperature as a function of temperature
* and pressure.
*
* A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
*
* Units = sqrt(kg/gmol)
*/
virtual double dA_DebyedT_TP(double temperature = -1.0,
double pressure = -1.0) const;
/**
* Value of the 2nd derivative of the Debye Huckel constant with
* respect to temperature as a function of temperature
* and pressure.
*
* A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
*
* Units = sqrt(kg/gmol)
*/
virtual double d2A_DebyedT2_TP(double temperature = -1.0,
double pressure = -1.0) const;
/**
* Value of the derivative of the Debye Huckel constant with
* respect to pressure, as a function of temperature
* and pressure.
*
* A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
*
* Units = sqrt(kg/gmol)
*/
virtual double dA_DebyedP_TP(double temperature = -1.0,
double pressure = -1.0) const;
/*
* AionicRadius()
*
@ -995,7 +1033,26 @@ namespace Cantera {
*/
mutable double m_IionicMolalityStoich;
public:
/**
* Form of the constant outside the Debye-Huckel term
* called A. It's normally a function of temperature
* and pressure. However, it can be set from the
* input file in order to aid in numerical comparisons.
* Acceptable forms:
*
* A_DEBYE_CONST 0
* A_DEBYE_WATER 1
*
* The A_DEBYE_WATER form may be used for water solvents
* with needs to cover varying temperatures and pressures.
* Note, the dielectric constant of water is a relatively
* strong function of T, and its variability must be
* accounted for,
*/
int m_form_A_Debye;
protected:
/**
* A_Debye -> this expression appears on the top of the
* ln actCoeff term in the general Debye-Huckel
@ -1059,6 +1116,12 @@ namespace Cantera {
*/
array_fp m_npActCoeff;
/**
* Pointer to the water property calculator
*/
WaterProps *m_waterProps;
/**
* Vector containing the species reference exp(-G/RT) functions
* at T = m_tlast
@ -1111,6 +1174,7 @@ namespace Cantera {
mutable array_fp m_lnActCoeffMolal;
mutable array_fp m_dlnActCoeffMolaldT;
mutable array_fp m_d2lnActCoeffMolaldT2;
mutable array_fp m_dlnActCoeffMolaldP;
private:
doublereal err(string msg) const;
@ -1122,10 +1186,11 @@ namespace Cantera {
* This function will be called to update the internally storred
* natural logarithm of the molality activity coefficients
*/
void _updatelnMolalityActCoeff() const;
void s_update_lnMolalityActCoeff() const;
void _updatedlnMolalityActCoeffdT() const;
void _updated2lnMolalityActCoeffdT2() const;
void s_update_dlnMolalityActCoeff_dT() const;
void s_update_d2lnMolalityActCoeff_dT2() const;
void s_update_dlnMolalityActCoeff_dP() const;
};
}

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

@ -146,10 +146,12 @@ namespace Cantera {
}
m_waterSS = b.m_waterSS;
m_densWaterSS = b.m_densWaterSS;
if (!m_waterProps) {
m_waterProps = new WaterProps(*b.m_waterProps);
} else {
m_waterProps = b.m_waterProps;
if (m_waterProps) {
delete m_waterProps;
m_waterProps = 0;
}
if (b.m_waterProps) {
m_waterProps = new WaterProps(*(b.m_waterProps));
}
m_waterSS = b.m_waterSS;
m_expg0_RT = b.m_expg0_RT;
@ -1323,7 +1325,6 @@ namespace Cantera {
break;
case A_DEBYE_WATER:
A = m_waterProps->ADebye(T, P, 0);
//A = WaterProps::ADebye(T, P, 0);
m_A_Debye = A;
break;
default:
@ -1482,10 +1483,6 @@ namespace Cantera {
* respect to temperature as a function of temperature
* and pressure.
*
* The default is to assume that it is equal to zero
* -> note, placeholder until a better formalism is
* put in place.
*
* units = A_Debye has units of sqrt(gmol kg-1).
* Temp has units of Kelvin.
*/

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

@ -1099,7 +1099,7 @@ namespace Cantera {
* strong function of T, and its variability must be
* accounted for,
*/
mutable int m_form_A_Debye;
int m_form_A_Debye;
protected:
/**
@ -1135,28 +1135,6 @@ namespace Cantera {
*/
mutable double m_A_Debye;
/**
* B_Debye -> this expression appears on the bottom of the
* ln actCoeff term in the general Debye-Huckel
* expression
* It depends on temperature
*
* B_Bebye = F / sqrt( epsilon R T / 2 )
*
* Units = sqrt(kg/gmol) / m
*
* Nominal value = 3.28640E9 sqrt(kg/gmol) / m
* based on:
* epsilon/epsilon_0 = 78.54
* (water at 25C)
* epsilon_0 = 8.854187817E12 C2 N-1 m-2
* e = 8.314472E3 kg m2 s-2 kmol-1 K-1
* F = 9.6485309E7 C kmol-1
* R = 8.314472E3 kg m2 s-2 kmol-1 K-1
* T = 298.15 K
*/
//double m_B_Debye;
/**
* Water standard state -> derived from the
* equation of state for water.
@ -1402,7 +1380,6 @@ namespace Cantera {
* This function will be called to update the internally storred
* natural logarithm of the molality activity coefficients
*/
//void s_updateDHlnMolalityActCoeff() const;
void s_update_lnMolalityActCoeff() const;
public:
void s_Pitzer_dlnMolalityActCoeff_dT() const;