ca4f9cc31d
does, putting in doxygen documentation. Tried to fix some holes in DebyeHuckel wherein calculations may be returned without updating the underlying water standard state.
2726 lines
77 KiB
C++
2726 lines
77 KiB
C++
/**
|
|
* @file DebyeHuckel.cpp
|
|
*
|
|
* Definitions of the DebyeHuckel object.
|
|
*/
|
|
/*
|
|
* Copywrite (2006) Sandia Corporation. Under the terms of
|
|
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
|
|
* U.S. Government retains certain rights in this software.
|
|
*/
|
|
/*
|
|
* $Id$
|
|
*/
|
|
//! Max function
|
|
#ifndef MAX
|
|
#define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
|
|
#endif
|
|
|
|
#include "DebyeHuckel.h"
|
|
#include "importCTML.h"
|
|
#include "WaterProps.h"
|
|
#include "WaterPDSS.h"
|
|
|
|
using namespace std;
|
|
|
|
namespace Cantera {
|
|
|
|
/*
|
|
* Default constructor
|
|
*/
|
|
DebyeHuckel::DebyeHuckel() :
|
|
MolalityVPSSTP(),
|
|
m_formDH(DHFORM_DILUTE_LIMIT),
|
|
m_formGC(2),
|
|
m_IionicMolality(0.0),
|
|
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_waterSS(0),
|
|
m_densWaterSS(1000.),
|
|
m_waterProps(0)
|
|
{
|
|
m_npActCoeff.resize(3);
|
|
m_npActCoeff[0] = 0.1127;
|
|
m_npActCoeff[1] = -0.01049;
|
|
m_npActCoeff[2] = 1.545E-3;
|
|
}
|
|
|
|
/*
|
|
* Working constructors
|
|
*
|
|
* The two constructors below are the normal way
|
|
* the phase initializes itself. They are shells that call
|
|
* the routine initThermo(), with a reference to the
|
|
* XML database to get the info for the phase.
|
|
*/
|
|
DebyeHuckel::DebyeHuckel(std::string inputFile, std::string id) :
|
|
MolalityVPSSTP(),
|
|
m_formDH(DHFORM_DILUTE_LIMIT),
|
|
m_formGC(2),
|
|
m_IionicMolality(0.0),
|
|
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_waterSS(0),
|
|
m_densWaterSS(1000.),
|
|
m_waterProps(0)
|
|
{
|
|
m_npActCoeff.resize(3);
|
|
m_npActCoeff[0] = 0.1127;
|
|
m_npActCoeff[1] = -0.01049;
|
|
m_npActCoeff[2] = 1.545E-3;
|
|
constructPhaseFile(inputFile, id);
|
|
}
|
|
|
|
DebyeHuckel::DebyeHuckel(XML_Node& phaseRoot, std::string id) :
|
|
MolalityVPSSTP(),
|
|
m_formDH(DHFORM_DILUTE_LIMIT),
|
|
m_formGC(2),
|
|
m_IionicMolality(0.0),
|
|
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_waterSS(0),
|
|
m_densWaterSS(1000.),
|
|
m_waterProps(0)
|
|
{
|
|
m_npActCoeff.resize(3);
|
|
m_npActCoeff[0] = 0.1127;
|
|
m_npActCoeff[1] = -0.01049;
|
|
m_npActCoeff[2] = 1.545E-3;
|
|
constructPhaseXML(phaseRoot, id);
|
|
}
|
|
|
|
/*
|
|
* Copy Constructor:
|
|
*
|
|
* Note this stuff will not work until the underlying phase
|
|
* has a working copy constructor
|
|
*/
|
|
DebyeHuckel::DebyeHuckel(const DebyeHuckel &b) :
|
|
MolalityVPSSTP(b)
|
|
{
|
|
/*
|
|
* Use the assignment operator to do the brunt
|
|
* of the work for the copy construtor.
|
|
*/
|
|
*this = b;
|
|
}
|
|
|
|
/*
|
|
* operator=()
|
|
*
|
|
* Note this stuff will not work until the underlying phase
|
|
* has a working assignment operator
|
|
*/
|
|
DebyeHuckel& DebyeHuckel::
|
|
operator=(const DebyeHuckel &b) {
|
|
if (&b != this) {
|
|
MolalityVPSSTP::operator=(b);
|
|
m_formDH = b.m_formDH;
|
|
m_formGC = b.m_formGC;
|
|
m_Aionic = b.m_Aionic;
|
|
m_npActCoeff = b.m_npActCoeff;
|
|
m_IionicMolality = b.m_IionicMolality;
|
|
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_waterSS) {
|
|
delete m_waterSS;
|
|
m_waterSS = 0;
|
|
}
|
|
if (b.m_waterSS) {
|
|
m_waterSS = new WaterPDSS(*(b.m_waterSS));
|
|
}
|
|
m_densWaterSS = b.m_densWaterSS;
|
|
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;
|
|
m_tmpV = b.m_tmpV;
|
|
m_speciesCharge_Stoich= b.m_speciesCharge_Stoich;
|
|
m_Beta_ij = b.m_Beta_ij;
|
|
m_lnActCoeffMolal = b.m_lnActCoeffMolal;
|
|
m_d2lnActCoeffMolaldT2= b.m_d2lnActCoeffMolaldT2;
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
|
|
/*
|
|
* ~DebyeHuckel(): (virtual)
|
|
*
|
|
* Destructor for DebyeHuckel. Release objects that
|
|
* it owns.
|
|
*/
|
|
DebyeHuckel::~DebyeHuckel() {
|
|
if (m_waterProps) {
|
|
delete m_waterProps; m_waterProps = 0;
|
|
}
|
|
if (m_waterSS) {
|
|
delete m_waterSS; m_waterSS = 0;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* duplMyselfAsThermoPhase():
|
|
*
|
|
* This routine operates at the ThermoPhase level to
|
|
* duplicate the current object. It uses the copy constructor
|
|
* defined above.
|
|
*/
|
|
ThermoPhase* DebyeHuckel::duplMyselfAsThermoPhase() {
|
|
DebyeHuckel* mtp = new DebyeHuckel(*this);
|
|
return (ThermoPhase *) mtp;
|
|
}
|
|
|
|
/*
|
|
* Equation of state type flag. The base class returns
|
|
* zero. Subclasses should define this to return a unique
|
|
* non-zero value. Constants defined for this purpose are
|
|
* listed in mix_defs.h.
|
|
*/
|
|
int DebyeHuckel::eosType() const {
|
|
int res;
|
|
switch (m_formGC) {
|
|
case 0:
|
|
res = cDebyeHuckel0;
|
|
break;
|
|
case 1:
|
|
res = cDebyeHuckel1;
|
|
break;
|
|
case 2:
|
|
res = cDebyeHuckel2;
|
|
break;
|
|
default:
|
|
throw CanteraError("eosType", "Unknown type");
|
|
break;
|
|
}
|
|
return res;
|
|
}
|
|
|
|
//
|
|
// -------- Molar Thermodynamic Properties of the Solution ---------------
|
|
//
|
|
/*
|
|
* Molar enthalpy of the solution. Units: J/kmol.
|
|
*/
|
|
doublereal DebyeHuckel::enthalpy_mole() const {
|
|
getPartialMolarEnthalpies(DATA_PTR(m_tmpV));
|
|
return mean_X(DATA_PTR(m_tmpV));
|
|
}
|
|
|
|
/*
|
|
* Molar internal energy of the solution. Units: J/kmol.
|
|
*
|
|
* This is calculated from the soln enthalpy and then
|
|
* subtracting pV.
|
|
*/
|
|
doublereal DebyeHuckel::intEnergy_mole() const {
|
|
double hh = enthalpy_mole();
|
|
double pres = pressure();
|
|
double molarV = 1.0/molarDensity();
|
|
double uu = hh - pres * molarV;
|
|
return uu;
|
|
}
|
|
|
|
/*
|
|
* Molar soln entropy at constant pressure. Units: J/kmol/K.
|
|
*
|
|
* This is calculated from the partial molar entropies.
|
|
*/
|
|
doublereal DebyeHuckel::entropy_mole() const {
|
|
getPartialMolarEntropies(DATA_PTR(m_tmpV));
|
|
return mean_X(DATA_PTR(m_tmpV));
|
|
}
|
|
|
|
// Molar Gibbs function. Units: J/kmol.
|
|
doublereal DebyeHuckel::gibbs_mole() const {
|
|
getChemPotentials(DATA_PTR(m_tmpV));
|
|
return mean_X(DATA_PTR(m_tmpV));
|
|
}
|
|
|
|
/*
|
|
* Molar heat capacity at constant pressure. Units: J/kmol/K.
|
|
*
|
|
* Returns the solution heat capacition at constant pressure.
|
|
* This is calculated from the partial molar heat capacities.
|
|
*/
|
|
doublereal DebyeHuckel::cp_mole() const {
|
|
getPartialMolarCp(DATA_PTR(m_tmpV));
|
|
double val = mean_X(DATA_PTR(m_tmpV));
|
|
return val;
|
|
}
|
|
|
|
/// Molar heat capacity at constant volume. Units: J/kmol/K.
|
|
doublereal DebyeHuckel::cv_mole() const {
|
|
//getPartialMolarCv(m_tmpV.begin());
|
|
//return mean_X(m_tmpV.begin());
|
|
err("not implemented");
|
|
return 0.0;
|
|
}
|
|
|
|
//
|
|
// ------- Mechanical Equation of State Properties ------------------------
|
|
//
|
|
|
|
/*
|
|
* Pressure. Units: Pa.
|
|
* For this incompressible system, we return the internally storred
|
|
* independent value of the pressure.
|
|
*/
|
|
doublereal DebyeHuckel::pressure() const {
|
|
return m_Pcurrent;
|
|
}
|
|
|
|
/*
|
|
* Set the pressure at constant temperature. Units: Pa.
|
|
* This method sets a constant within the object.
|
|
* The mass density is not a function of pressure.
|
|
*/
|
|
void DebyeHuckel::setPressure(doublereal p) {
|
|
|
|
#ifdef DEBUG_MODE
|
|
//printf("setPressure: %g\n", p);
|
|
#endif
|
|
/*
|
|
* Store the current pressure
|
|
*/
|
|
m_Pcurrent = p;
|
|
|
|
/*
|
|
* update the standard state thermo
|
|
* -> This involves calling the water function and setting the pressure
|
|
*/
|
|
_updateStandardStateThermo();
|
|
|
|
if (m_waterSS) {
|
|
|
|
/*
|
|
* Store the internal density of the water SS.
|
|
* Note, we would have to do this for all other
|
|
* species if they had pressure dependent properties.
|
|
*/
|
|
m_densWaterSS = m_waterSS->density();
|
|
}
|
|
/*
|
|
* Calculate all of the other standard volumes
|
|
* -> note these are constant for now
|
|
*/
|
|
/*
|
|
* Get the partial molar volumes of all of the
|
|
* species. -> note this is a lookup for
|
|
* water, here since it was done above.
|
|
*/
|
|
double *vbar = &m_pp[0];
|
|
getPartialMolarVolumes(vbar);
|
|
|
|
/*
|
|
* Get mole fractions of all species.
|
|
*/
|
|
double *x = &m_tmpV[0];
|
|
getMoleFractions(x);
|
|
|
|
/*
|
|
* Calculate the solution molar volume and the
|
|
* solution density.
|
|
*/
|
|
doublereal vtotal = 0.0;
|
|
for (int i = 0; i < m_kk; i++) {
|
|
vtotal += vbar[i] * x[i];
|
|
}
|
|
doublereal dd = meanMolecularWeight() / vtotal;
|
|
|
|
/*
|
|
* Now, update the State class with the results. This
|
|
* stores the density.
|
|
*/
|
|
State::setDensity(dd);
|
|
|
|
}
|
|
|
|
/*
|
|
* The isothermal compressibility. Units: 1/Pa.
|
|
* The isothermal compressibility is defined as
|
|
* \f[
|
|
* \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
|
|
* \f]
|
|
*
|
|
* It's equal to zero for this model, since the molar volume
|
|
* doesn't change with pressure or temperature.
|
|
*/
|
|
doublereal DebyeHuckel::isothermalCompressibility() const {
|
|
throw CanteraError("DebyeHuckel::isothermalCompressibility",
|
|
"unimplemented");
|
|
return 0.0;
|
|
}
|
|
|
|
/*
|
|
* The thermal expansion coefficient. Units: 1/K.
|
|
* The thermal expansion coefficient is defined as
|
|
*
|
|
* \f[
|
|
* \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
|
|
* \f]
|
|
*
|
|
* It's equal to zero for this model, since the molar volume
|
|
* doesn't change with pressure or temperature.
|
|
*/
|
|
doublereal DebyeHuckel::thermalExpansionCoeff() const {
|
|
throw CanteraError("DebyeHuckel::thermalExpansionCoeff",
|
|
"unimplemented");
|
|
return 0.0;
|
|
}
|
|
|
|
/*
|
|
* Overwritten setDensity() function is necessary because the
|
|
* density is not an indendent variable.
|
|
*
|
|
* This function will now throw an error condition
|
|
*
|
|
* @internal May have to adjust the strategy here to make
|
|
* the eos for these materials slightly compressible, in order
|
|
* to create a condition where the density is a function of
|
|
* the pressure.
|
|
*
|
|
* This function will now throw an error condition.
|
|
*
|
|
* NOTE: This is an overwritten function from the State.h
|
|
* class
|
|
*/
|
|
void DebyeHuckel::setDensity(doublereal rho) {
|
|
double dens = density();
|
|
if (rho != dens) {
|
|
throw CanteraError("Idea;MolalSoln::setDensity",
|
|
"Density is not an independent variable");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Overwritten setMolarDensity() function is necessary because the
|
|
* density is not an indendent variable.
|
|
*
|
|
* This function will now throw an error condition.
|
|
*
|
|
* NOTE: This is a virtual function, overwritten function from the State.h
|
|
* class
|
|
*/
|
|
void DebyeHuckel::setMolarDensity(doublereal conc) {
|
|
double concI = molarDensity();
|
|
if (conc != concI) {
|
|
throw CanteraError("Idea;MolalSoln::setMolarDensity",
|
|
"molarDensity/density is not an independent variable");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Overwritten setTemperature(double) from State.h. This
|
|
* function sets the temperature, and makes sure that
|
|
* the value propagates to underlying objects.
|
|
*/
|
|
void DebyeHuckel::setTemperature(doublereal temp) {
|
|
_updateStandardStateThermo();
|
|
State::setTemperature(temp);
|
|
}
|
|
|
|
|
|
//
|
|
// ------- Activities and Activity Concentrations
|
|
//
|
|
|
|
/*
|
|
* This method returns an array of generalized concentrations
|
|
* \f$ C_k\f$ that are defined such that
|
|
* \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
|
|
* is a standard concentration
|
|
* defined below. These generalized concentrations are used
|
|
* by kinetics manager classes to compute the forward and
|
|
* reverse rates of elementary reactions.
|
|
*
|
|
* @param c Array of generalized concentrations. The
|
|
* units depend upon the implementation of the
|
|
* reaction rate expressions within the phase.
|
|
*/
|
|
void DebyeHuckel::getActivityConcentrations(doublereal* c) const {
|
|
double c_solvent = standardConcentration();
|
|
getActivities(c);
|
|
for (int k = 0; k < m_kk; k++) {
|
|
c[k] *= c_solvent;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* The standard concentration \f$ C^0_k \f$ used to normalize
|
|
* the generalized concentration. In many cases, this quantity
|
|
* will be the same for all species in a phase - for example,
|
|
* for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
|
|
* reason, this method returns a single value, instead of an
|
|
* array. However, for phases in which the standard
|
|
* concentration is species-specific (e.g. surface species of
|
|
* different sizes), this method may be called with an
|
|
* optional parameter indicating the species.
|
|
*
|
|
* For the time being we will use the concentration of pure
|
|
* solvent for the the standard concentration of all species.
|
|
* This has the effect of making reaction rates
|
|
* based on the molality of species proportional to the
|
|
* molality of the species.
|
|
*/
|
|
doublereal DebyeHuckel::standardConcentration(int k) const {
|
|
double mvSolvent = m_speciesSize[m_indexSolvent];
|
|
return 1.0 / mvSolvent;
|
|
}
|
|
|
|
/*
|
|
* Returns the natural logarithm of the standard
|
|
* concentration of the kth species
|
|
*/
|
|
doublereal DebyeHuckel::logStandardConc(int k) const {
|
|
double c_solvent = standardConcentration(k);
|
|
return log(c_solvent);
|
|
}
|
|
|
|
/*
|
|
* Returns the units of the standard and general concentrations
|
|
* Note they have the same units, as their divisor is
|
|
* defined to be equal to the activity of the kth species
|
|
* in the solution, which is unitless.
|
|
*
|
|
* This routine is used in print out applications where the
|
|
* units are needed. Usually, MKS units are assumed throughout
|
|
* the program and in the XML input files.
|
|
*
|
|
* On return uA contains the powers of the units (MKS assumed)
|
|
* of the standard concentrations and generalized concentrations
|
|
* for the kth species.
|
|
*
|
|
* uA[0] = kmol units - default = 1
|
|
* uA[1] = m units - default = -nDim(), the number of spatial
|
|
* dimensions in the Phase class.
|
|
* uA[2] = kg units - default = 0;
|
|
* uA[3] = Pa(pressure) units - default = 0;
|
|
* uA[4] = Temperature units - default = 0;
|
|
* uA[5] = time units - default = 0
|
|
*/
|
|
void DebyeHuckel::getUnitsStandardConc(double *uA, int k, int sizeUA) {
|
|
for (int i = 0; i < sizeUA; i++) {
|
|
if (i == 0) uA[0] = 1.0;
|
|
if (i == 1) uA[1] = -nDim();
|
|
if (i == 2) uA[2] = 0.0;
|
|
if (i == 3) uA[3] = 0.0;
|
|
if (i == 4) uA[4] = 0.0;
|
|
if (i == 5) uA[5] = 0.0;
|
|
}
|
|
}
|
|
|
|
|
|
/*
|
|
* Get the array of non-dimensional activities at
|
|
* the current solution temperature, pressure, and
|
|
* solution concentration.
|
|
* (note solvent activity coefficient is on the molar scale).
|
|
*
|
|
*/
|
|
void DebyeHuckel::getActivities(doublereal* ac) const {
|
|
_updateStandardStateThermo();
|
|
/*
|
|
* Update the molality array, m_molalities()
|
|
* This requires an update due to mole fractions
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent) {
|
|
ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal[k]);
|
|
}
|
|
}
|
|
double xmolSolvent = moleFraction(m_indexSolvent);
|
|
ac[m_indexSolvent] =
|
|
exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent;
|
|
}
|
|
|
|
/*
|
|
* getMolalityActivityCoefficients() (virtual, const)
|
|
*
|
|
* Get the array of non-dimensional Molality based
|
|
* activity coefficients at
|
|
* the current solution temperature, pressure, and
|
|
* solution concentration.
|
|
* (note solvent activity coefficient is on the molar scale).
|
|
*
|
|
* Note, most of the work is done in an internal private routine
|
|
*/
|
|
void DebyeHuckel::
|
|
getMolalityActivityCoefficients(doublereal* acMolality) const {
|
|
_updateStandardStateThermo();
|
|
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]);
|
|
}
|
|
}
|
|
|
|
//
|
|
// ------ Partial Molar Properties of the Solution -----------------
|
|
//
|
|
/*
|
|
* Get the species chemical potentials. Units: J/kmol.
|
|
*
|
|
* This function returns a vector of chemical potentials of the
|
|
* species in solution.
|
|
*
|
|
* \f[
|
|
* \mu_k = \mu^{o}_k(T,P) + R T ln(m_k)
|
|
* \f]
|
|
*
|
|
* \f[
|
|
* \mu_solvent = \mu^{o}_solvent(T,P) +
|
|
* R T ((X_solvent - 1.0) / X_solvent)
|
|
* \f]
|
|
*/
|
|
void DebyeHuckel::getChemPotentials(doublereal* mu) const{
|
|
double xx;
|
|
const double xxSmall = 1.0E-150;
|
|
/*
|
|
* First get the standard chemical potentials in
|
|
* molar form.
|
|
* -> this requires updates of standard state as a function
|
|
* of T and P
|
|
*/
|
|
getStandardChemPotentials(mu);
|
|
/*
|
|
* Update the activity coefficients
|
|
* This also updates the internal molality array.
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
/*
|
|
*
|
|
*/
|
|
doublereal RT = GasConstant * temperature();
|
|
double xmolSolvent = moleFraction(m_indexSolvent);
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (m_indexSolvent != k) {
|
|
xx = MAX(m_molalities[k], xxSmall);
|
|
mu[k] += RT * (log(xx) + m_lnActCoeffMolal[k]);
|
|
}
|
|
}
|
|
xx = MAX(xmolSolvent, xxSmall);
|
|
mu[m_indexSolvent] +=
|
|
RT * (log(xx) + m_lnActCoeffMolal[m_indexSolvent]);
|
|
}
|
|
|
|
|
|
/*
|
|
* Returns an array of partial molar enthalpies for the species
|
|
* in the mixture.
|
|
* Units (J/kmol)
|
|
*
|
|
* We calculate this quantity partially from the relation and
|
|
* partially by calling the standard state enthalpy function.
|
|
*
|
|
* hbar_i = - T**2 * d(chemPot_i/T)/dT
|
|
*
|
|
* We calculate
|
|
*/
|
|
void DebyeHuckel::getPartialMolarEnthalpies(doublereal* hbar) const {
|
|
/*
|
|
* Get the nondimensional standard state enthalpies
|
|
*/
|
|
getEnthalpy_RT(hbar);
|
|
/*
|
|
* Dimensionalize it.
|
|
*/
|
|
double T = temperature();
|
|
double RT = GasConstant * T;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
hbar[k] *= RT;
|
|
}
|
|
/*
|
|
* Check to see whether activity coefficients are temperature
|
|
* dependent. If they are, then calculate the their temperature
|
|
* derivatives and add them into the result.
|
|
*/
|
|
double dAdT = dA_DebyedT_TP();
|
|
if (dAdT != 0.0) {
|
|
/*
|
|
* Update the activity coefficients, This also update the
|
|
* internally storred molalities.
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
s_update_dlnMolalityActCoeff_dT();
|
|
double RTT = GasConstant * T * T;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
hbar[k] -= RTT * m_dlnActCoeffMolaldT[k];
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
*
|
|
* getPartialMolarEntropies() (virtual, const)
|
|
*
|
|
* Returns an array of partial molar entropies of the species in the
|
|
* solution. Units: J/kmol.
|
|
*
|
|
* Maxwell's equations provide an insight in how to calculate this
|
|
* (p.215 Smith and Van Ness)
|
|
*
|
|
* d(chemPot_i)/dT = -sbar_i
|
|
*
|
|
* For this phase, the partial molar entropies are equal to the
|
|
* SS species entropies plus the ideal solution contribution.following
|
|
* contribution:
|
|
* \f[
|
|
* \bar s_k(T,P) = \hat s^0_k(T) - R log(M0 * molality[k])
|
|
* \f]
|
|
* \f[
|
|
* \bar s_solvent(T,P) = \hat s^0_solvent(T)
|
|
* - R ((xmolSolvent - 1.0) / xmolSolvent)
|
|
* \f]
|
|
*
|
|
* The reference-state pure-species entropies,\f$ \hat s^0_k(T) \f$,
|
|
* at the reference pressure, \f$ P_{ref} \f$, are computed by the
|
|
* species thermodynamic
|
|
* property manager. They are polynomial functions of temperature.
|
|
* @see SpeciesThermo
|
|
*/
|
|
void DebyeHuckel::
|
|
getPartialMolarEntropies(doublereal* sbar) const {
|
|
int k;
|
|
/*
|
|
* Get the standard state entropies at the temperature
|
|
* and pressure of the solution.
|
|
*/
|
|
getEntropy_R(sbar);
|
|
/*
|
|
* Dimensionalize the entropies
|
|
*/
|
|
doublereal R = GasConstant;
|
|
for (k = 0; k < m_kk; k++) {
|
|
sbar[k] *= R;
|
|
}
|
|
/*
|
|
* Update the activity coefficients, This also update the
|
|
* internally storred molalities.
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
/*
|
|
* First we will add in the obvious dependence on the T
|
|
* term out front of the log activity term
|
|
*/
|
|
doublereal mm;
|
|
for (k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent) {
|
|
mm = fmaxx(SmallNumber, m_molalities[k]);
|
|
sbar[k] -= R * (log(mm) + m_lnActCoeffMolal[k]);
|
|
}
|
|
}
|
|
double xmolSolvent = moleFraction(m_indexSolvent);
|
|
mm = fmaxx(SmallNumber, xmolSolvent);
|
|
sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
|
|
/*
|
|
* Check to see whether activity coefficients are temperature
|
|
* dependent. If they are, then calculate the their temperature
|
|
* derivatives and add them into the result.
|
|
*/
|
|
double dAdT = dA_DebyedT_TP();
|
|
if (dAdT != 0.0) {
|
|
s_update_dlnMolalityActCoeff_dT();
|
|
double RT = R * temperature();
|
|
for (k = 0; k < m_kk; k++) {
|
|
sbar[k] -= RT * m_dlnActCoeffMolaldT[k];
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* getPartialMolarVolumes() (virtual, const)
|
|
*
|
|
* returns an array of partial molar volumes of the species
|
|
* in the solution. Units: m^3 kmol-1.
|
|
*
|
|
* For this solution, the partial molar volumes are normally
|
|
* equal to theconstant species molar volumes, except
|
|
* when the activity coefficients depend on pressure.
|
|
*
|
|
* The general relation is
|
|
*
|
|
* vbar_i = d(chemPot_i)/dP at const T, n
|
|
*
|
|
* = V0_i + d(Gex)/dP)_T,M
|
|
*
|
|
* = V0_i + RT d(lnActCoeffi)dP _T,M
|
|
*
|
|
*/
|
|
void DebyeHuckel::getPartialMolarVolumes(doublereal* vbar) const {
|
|
getStandardVolumes(vbar);
|
|
/*
|
|
* Update the derivatives wrt the activity coefficients.
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
s_update_dlnMolalityActCoeff_dP();
|
|
double T = temperature();
|
|
double RT = GasConstant * T;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
vbar[k] += RT * m_dlnActCoeffMolaldP[k];
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Partial molar heat capacity of the solution:
|
|
* The kth partial molar heat capacity is equal to
|
|
* the temperature derivative of the partial molar
|
|
* enthalpy of the kth species in the solution at constant
|
|
* P and composition (p. 220 Smith and Van Ness).
|
|
*
|
|
* Cp = -T d2(chemPot_i)/dT2
|
|
*/
|
|
void DebyeHuckel::getPartialMolarCp(doublereal* cpbar) const {
|
|
/*
|
|
* Get the nondimensional gibbs standard state of the
|
|
* species at the T and P of the solution.
|
|
*/
|
|
getCp_R(cpbar);
|
|
|
|
for (int k = 0; k < m_kk; k++) {
|
|
cpbar[k] *= GasConstant;
|
|
}
|
|
|
|
/*
|
|
* Check to see whether activity coefficients are temperature
|
|
* dependent. If they are, then calculate the their temperature
|
|
* derivatives and add them into the result.
|
|
*/
|
|
double dAdT = dA_DebyedT_TP();
|
|
if (dAdT != 0.0) {
|
|
/*
|
|
* Update the activity coefficients, This also update the
|
|
* internally storred molalities.
|
|
*/
|
|
s_update_lnMolalityActCoeff();
|
|
s_update_dlnMolalityActCoeff_dT();
|
|
s_update_d2lnMolalityActCoeff_dT2();
|
|
double T = temperature();
|
|
double RT = GasConstant * T;
|
|
double RTT = RT * T;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
cpbar[k] -= (2.0 * RT * m_dlnActCoeffMolaldT[k] +
|
|
RTT * m_d2lnActCoeffMolaldT2[k]);
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* -------- Properties of the Standard State of the Species
|
|
* in the Solution ------------------
|
|
*/
|
|
|
|
/*
|
|
* getStandardChemPotentials() (virtual, const)
|
|
*
|
|
*
|
|
* Get the standard state chemical potentials of the species.
|
|
* This is the array of chemical potentials at unit activity
|
|
* (Mole fraction scale)
|
|
* \f$ \mu^0_k(T,P) \f$.
|
|
* We define these here as the chemical potentials of the pure
|
|
* species at the temperature and pressure of the solution.
|
|
* This function is used in the evaluation of the
|
|
* equilibrium constant Kc. Therefore, Kc will also depend
|
|
* on T and P. This is the norm for liquid and solid systems.
|
|
*
|
|
* units = J / kmol
|
|
*/
|
|
void DebyeHuckel::getStandardChemPotentials(doublereal* mu) const {
|
|
_updateStandardStateThermo();
|
|
getGibbs_ref(mu);
|
|
doublereal pref;
|
|
doublereal delta_p;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
pref = m_spthermo->refPressure(k);
|
|
delta_p = m_Pcurrent - pref;
|
|
mu[k] += delta_p * m_speciesSize[k];
|
|
}
|
|
if (m_waterSS) {
|
|
mu[0] = m_waterSS->gibbs_mole();
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Get the nondimensional gibbs function for the species
|
|
* standard states at the current T and P of the solution.
|
|
*
|
|
* \f[
|
|
* \mu^0_k(T,P) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k
|
|
* \f]
|
|
* where \f$V_k\f$ is the molar volume of pure species <I>k</I>.
|
|
* \f$ \mu^{ref}_k(T)\f$ is the chemical potential of pure
|
|
* species <I>k</I> at the reference pressure, \f$P_{ref}\f$.
|
|
*
|
|
* @param grt Vector of length m_kk, which on return sr[k]
|
|
* will contain the nondimensional
|
|
* standard state gibbs function for species k.
|
|
*/
|
|
void DebyeHuckel::getGibbs_RT(doublereal* grt) const {
|
|
_updateStandardStateThermo();
|
|
getPureGibbs(grt);
|
|
doublereal invRT = 1.0 / _RT();
|
|
for (int k = 0; k < m_kk; k++) {
|
|
grt[k] *= invRT;
|
|
}
|
|
}
|
|
|
|
/*
|
|
*
|
|
* getPureGibbs()
|
|
*
|
|
* Get the Gibbs functions for the pure species
|
|
* at the current <I>T</I> and <I>P</I> of the solution.
|
|
* We assume an incompressible constant partial molar
|
|
* volume here:
|
|
* \f[
|
|
* \mu^0_k(T,p) = \mu^{ref}_k(T) + (P - P_{ref}) * V_k
|
|
* \f]
|
|
* where \f$V_k\f$ is the molar volume of pure species <I>k<\I>.
|
|
* \f$ u^{ref}_k(T)\f$ is the chemical potential of pure
|
|
* species <I>k<\I> at the reference pressure, \f$P_{ref}\f$.
|
|
*/
|
|
void DebyeHuckel::getPureGibbs(doublereal* gpure) const {
|
|
getStandardChemPotentials(gpure);
|
|
}
|
|
|
|
/*
|
|
* Get the array of nondimensional Enthalpy functions for the ss
|
|
* species at the current <I>T</I> and <I>P</I> of the solution.
|
|
* We assume an incompressible constant partial molar
|
|
* volume here:
|
|
* \f[
|
|
* h^0_k(T,P) = h^{ref}_k(T) + (P - P_{ref}) * V_k
|
|
* \f]
|
|
* where \f$V_k\f$ is the molar volume of SS species <I>k<\I>.
|
|
* \f$ h^{ref}_k(T)\f$ is the enthalpy of the SS
|
|
* species <I>k<\I> at the reference pressure, \f$P_{ref}\f$.
|
|
*/
|
|
void DebyeHuckel::
|
|
getEnthalpy_RT(doublereal* hrt) const {
|
|
_updateStandardStateThermo();
|
|
getEnthalpy_RT_ref(hrt);
|
|
doublereal pref;
|
|
doublereal delta_p;
|
|
double RT = _RT();
|
|
for (int k = 0; k < m_kk; k++) {
|
|
pref = m_spthermo->refPressure(k);
|
|
delta_p = m_Pcurrent - pref;
|
|
hrt[k] += delta_p/ RT * m_speciesSize[k];
|
|
}
|
|
if (m_waterSS) {
|
|
hrt[0] = m_waterSS->enthalpy_mole();
|
|
hrt[0] /= RT;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Get the nondimensional Entropies for the species
|
|
* standard states at the current T and P of the solution.
|
|
*
|
|
* Note, this is equal to the reference state entropies
|
|
* due to the zero volume expansivity:
|
|
* i.e., (dS/dp)_T = (dV/dT)_P = 0.0
|
|
*
|
|
* The solvent water entropy is obtained from a pure water
|
|
* equation of state model.
|
|
*
|
|
* @param sr Vector of length m_kk, which on return sr[k]
|
|
* will contain the nondimensional
|
|
* standard state entropy of species k.
|
|
*/
|
|
void DebyeHuckel::
|
|
getEntropy_R(doublereal* sr) const {
|
|
_updateStandardStateThermo();
|
|
getEntropy_R_ref(sr);
|
|
if (m_waterSS) {
|
|
sr[0] = m_waterSS->entropy_mole();
|
|
sr[0] /= GasConstant;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Get the nondimensional heat capacity at constant pressure
|
|
* function for the species
|
|
* standard states at the current T and P of the solution.
|
|
* \f[
|
|
* Cp^0_k(T,P) = Cp^{ref}_k(T)
|
|
* \f]
|
|
* where \f$V_k\f$ is the molar volume of pure species <I>k</I>.
|
|
* \f$ Cp^{ref}_k(T)\f$ is the constant pressure heat capacity
|
|
* of species <I>k</I> at the reference pressure, \f$p_{ref}\f$.
|
|
*
|
|
* The solvent water heat capacity is obtained from a pure water
|
|
* equation of state model.
|
|
*
|
|
* @param cpr Vector of length m_kk, which on return cpr[k]
|
|
* will contain the nondimensional
|
|
* constant pressure heat capacity for species k.
|
|
*/
|
|
void DebyeHuckel::getCp_R(doublereal* cpr) const {
|
|
_updateStandardStateThermo();
|
|
getCp_R_ref(cpr);
|
|
if (m_waterSS) {
|
|
cpr[0] = m_waterSS->cp_mole();
|
|
cpr[0] /= GasConstant;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Get the molar volumes of each species in their standard
|
|
* states at the current
|
|
* <I>T</I> and <I>P</I> of the solution.
|
|
* units = m^3 / kmol
|
|
*/
|
|
void DebyeHuckel::getStandardVolumes(doublereal *vol) const {
|
|
_updateStandardStateThermo();
|
|
copy(m_speciesSize.begin(),
|
|
m_speciesSize.end(), vol);
|
|
if (m_waterSS) {
|
|
double dd = m_waterSS->density();
|
|
vol[0] = molecularWeight(0)/dd;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* ------ Thermodynamic Values for the Species Reference States ---
|
|
*/
|
|
|
|
// -> This is handled by VPStandardStatesTP
|
|
|
|
/*
|
|
* -------------- Utilities -------------------------------
|
|
*/
|
|
|
|
/*
|
|
* Initialization routine for a DebyeHuckel phase.
|
|
*
|
|
* This is a virtual routine. This routine will call initThermo()
|
|
* for the parent class as well.
|
|
*/
|
|
void DebyeHuckel::initThermo() {
|
|
MolalityVPSSTP::initThermo();
|
|
initLengths();
|
|
}
|
|
|
|
/*
|
|
* constructPhaseFile
|
|
*
|
|
* Initialization of a Debye-Huckel phase using an
|
|
* xml file.
|
|
*
|
|
* This routine is a precursor to initThermo(XML_Node*)
|
|
* routine, which does most of the work.
|
|
*
|
|
* @param infile XML file containing the description of the
|
|
* phase
|
|
*
|
|
* @param id Optional parameter identifying the name of the
|
|
* phase. If none is given, the first XML
|
|
* phase element will be used.
|
|
*/
|
|
void DebyeHuckel::constructPhaseFile(std::string inputFile, std::string id) {
|
|
|
|
if (inputFile.size() == 0) {
|
|
throw CanteraError("DebyeHuckel::initThermo",
|
|
"input file is null");
|
|
}
|
|
std::string path = findInputFile(inputFile);
|
|
ifstream fin(path.c_str());
|
|
if (!fin) {
|
|
throw CanteraError("DebyeHuckel::initThermo","could not open "
|
|
+path+" for reading.");
|
|
}
|
|
/*
|
|
* The phase object automatically constructs an XML object.
|
|
* Use this object to store information.
|
|
*/
|
|
XML_Node &phaseNode_XML = xml();
|
|
XML_Node *fxml = new XML_Node();
|
|
fxml->build(fin);
|
|
XML_Node *fxml_phase = findXMLPhase(fxml, id);
|
|
if (!fxml_phase) {
|
|
throw CanteraError("DebyeHuckel::initThermo",
|
|
"ERROR: Can not find phase named " +
|
|
id + " in file named " + inputFile);
|
|
}
|
|
fxml_phase->copy(&phaseNode_XML);
|
|
constructPhaseXML(*fxml_phase, id);
|
|
delete fxml;
|
|
}
|
|
|
|
/**
|
|
* interp_est() (static)
|
|
*
|
|
* utility function to assign an integer value from a string
|
|
* for the ElectrolyteSpeciesType field.
|
|
*/
|
|
static int interp_est(std::string estString) {
|
|
const char *cc = estString.c_str();
|
|
if (!strcasecmp(cc, "solvent")) {
|
|
return cEST_solvent;
|
|
} else if (!strcasecmp(cc, "chargedspecies")) {
|
|
return cEST_chargedSpecies;
|
|
} else if (!strcasecmp(cc, "weakAcidAssociated")) {
|
|
return cEST_weakAcidAssociated;
|
|
} else if (!strcasecmp(cc, "strongAcidAssociated")) {
|
|
return cEST_strongAcidAssociated;
|
|
} else if (!strcasecmp(cc, "polarNeutral")) {
|
|
return cEST_polarNeutral;
|
|
} else if (!strcasecmp(cc, "nonpolarNeutral")) {
|
|
return cEST_nonpolarNeutral;
|
|
}
|
|
int retn, rval;
|
|
if ((retn = sscanf(cc, "%d", &rval)) != 1) {
|
|
return -1;
|
|
}
|
|
return rval;
|
|
}
|
|
|
|
/*
|
|
* Import and initialize a DebyeHuckel phase
|
|
* specification in an XML tree into the current object.
|
|
* Here we read an XML description of the phase.
|
|
* We import descriptions of the elements that make up the
|
|
* species in a phase.
|
|
* We import information about the species, including their
|
|
* reference state thermodynamic polynomials. We then freeze
|
|
* the state of the species.
|
|
*
|
|
* Then, we read the species molar volumes from the xml
|
|
* tree to finish the initialization.
|
|
*
|
|
* @param phaseNode This object must be the phase node of a
|
|
* complete XML tree
|
|
* description of the phase, including all of the
|
|
* species data. In other words while "phase" must
|
|
* point to an XML phase object, it must have
|
|
* sibling nodes "speciesData" that describe
|
|
* the species in the phase.
|
|
* @param id ID of the phase. If nonnull, a check is done
|
|
* to see if phaseNode is pointing to the phase
|
|
* with the correct id.
|
|
*/
|
|
void DebyeHuckel::constructPhaseXML(XML_Node& phaseNode, std::string id) {
|
|
|
|
if (id.size() > 0) {
|
|
std::string idp = phaseNode.id();
|
|
if (idp != id) {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"phasenode and Id are incompatible");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Find the Thermo XML node
|
|
*/
|
|
if (!phaseNode.hasChild("thermo")) {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"no thermo XML node");
|
|
}
|
|
XML_Node& thermoNode = phaseNode.child("thermo");
|
|
|
|
/*
|
|
* Possibly change the form of the standard concentrations
|
|
*/
|
|
if (thermoNode.hasChild("standardConc")) {
|
|
XML_Node& scNode = thermoNode.child("standardConc");
|
|
m_formGC = 2;
|
|
std::string formString = scNode.attrib("model");
|
|
if (formString != "") {
|
|
if (formString == "unity") {
|
|
m_formGC = 0;
|
|
printf("exit standardConc = unity not done\n");
|
|
exit(-1);
|
|
} else if (formString == "molar_volume") {
|
|
m_formGC = 1;
|
|
printf("exit standardConc = molar_volume not done\n");
|
|
exit(-1);
|
|
} else if (formString == "solvent_volume") {
|
|
m_formGC = 2;
|
|
} else {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"Unknown standardConc model: " + formString);
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Get the Name of the Solvent:
|
|
* <solvent> solventName </solvent>
|
|
*/
|
|
std::string solventName = "";
|
|
if (thermoNode.hasChild("solvent")) {
|
|
XML_Node& scNode = thermoNode.child("solvent");
|
|
vector<std::string> nameSolventa;
|
|
getStringArray(scNode, nameSolventa);
|
|
int nsp = static_cast<int>(nameSolventa.size());
|
|
if (nsp != 1) {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"badly formed solvent XML node");
|
|
}
|
|
solventName = nameSolventa[0];
|
|
}
|
|
|
|
/*
|
|
* Determine the form of the Debye-Huckel model,
|
|
* m_formDH. We will use this information to size arrays below.
|
|
*/
|
|
if (thermoNode.hasChild("activityCoefficients")) {
|
|
XML_Node& scNode = thermoNode.child("activityCoefficients");
|
|
m_formDH = DHFORM_DILUTE_LIMIT;
|
|
std::string formString = scNode.attrib("model");
|
|
if (formString != "") {
|
|
if (formString == "Dilute_limit") {
|
|
m_formDH = DHFORM_DILUTE_LIMIT;
|
|
} else if (formString == "Bdot_with_variable_a") {
|
|
m_formDH = DHFORM_BDOT_AK ;
|
|
} else if (formString == "Bdot_with_common_a") {
|
|
m_formDH = DHFORM_BDOT_ACOMMON;
|
|
} else if (formString == "Beta_ij") {
|
|
m_formDH = DHFORM_BETAIJ;
|
|
} else if (formString == "Pitzer_with_Beta_ij") {
|
|
m_formDH = DHFORM_PITZER_BETAIJ;
|
|
} else {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"Unknown standardConc model: " + formString);
|
|
}
|
|
}
|
|
} else {
|
|
/*
|
|
* If there is no XML node named "activityCoefficients", assume
|
|
* that we are doing the extreme dilute limit assumption
|
|
*/
|
|
m_formDH = DHFORM_DILUTE_LIMIT;
|
|
}
|
|
|
|
/*
|
|
* Call the Cantera importPhase() function. This will import
|
|
* all of the species into the phase. This will also handle
|
|
* all of the solvent and solute standard states
|
|
*/
|
|
bool m_ok = importPhase(phaseNode, this);
|
|
if (!m_ok) {
|
|
throw CanteraError("DebyeHuckel::constructPhaseXML",
|
|
"importPhase failed ");
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Process the XML file after species are set up.
|
|
*
|
|
* This gets called from importPhase(). It processes the XML file
|
|
* after the species are set up. This is the main routine for
|
|
* reading in activity coefficient parameters.
|
|
*
|
|
* @param phaseNode This object must be the phase node of a
|
|
* complete XML tree
|
|
* description of the phase, including all of the
|
|
* species data. In other words while "phase" must
|
|
* point to an XML phase object, it must have
|
|
* sibling nodes "speciesData" that describe
|
|
* the species in the phase.
|
|
* @param id ID of the phase. If nonnull, a check is done
|
|
* to see if phaseNode is pointing to the phase
|
|
* with the correct id.
|
|
*/
|
|
void DebyeHuckel::
|
|
initThermoXML(XML_Node& phaseNode, std::string id) {
|
|
int k;
|
|
std::string stemp;
|
|
/*
|
|
* Find the Thermo XML node
|
|
*/
|
|
if (!phaseNode.hasChild("thermo")) {
|
|
throw CanteraError("HMWSoln::initThermoXML",
|
|
"no thermo XML node");
|
|
}
|
|
XML_Node& thermoNode = phaseNode.child("thermo");
|
|
|
|
/*
|
|
* Initialize all of the lengths of arrays in the object
|
|
* now that we know what species are in the phase.
|
|
*/
|
|
initThermo();
|
|
|
|
/*
|
|
* Reconcile the solvent name and index.
|
|
*/
|
|
/*
|
|
* Get the Name of the Solvent:
|
|
* <solvent> solventName </solvent>
|
|
*/
|
|
std::string solventName = "";
|
|
if (thermoNode.hasChild("solvent")) {
|
|
XML_Node& scNode = thermoNode.child("solvent");
|
|
vector<std::string> nameSolventa;
|
|
getStringArray(scNode, nameSolventa);
|
|
int nsp = static_cast<int>(nameSolventa.size());
|
|
if (nsp != 1) {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"badly formed solvent XML node");
|
|
}
|
|
solventName = nameSolventa[0];
|
|
}
|
|
for (k = 0; k < m_kk; k++) {
|
|
std::string sname = speciesName(k);
|
|
if (solventName == sname) {
|
|
m_indexSolvent = k;
|
|
break;
|
|
}
|
|
}
|
|
if (m_indexSolvent == -1) {
|
|
cout << "DebyeHuckel::initThermoXML: Solvent Name not found"
|
|
<< endl;
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Solvent name not found");
|
|
}
|
|
if (m_indexSolvent != 0) {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Solvent " + solventName +
|
|
" should be first species");
|
|
}
|
|
|
|
/*
|
|
* Now go get the specification of the standard states for
|
|
* species in the solution. This includes the molar volumes
|
|
* data blocks for incompressible species.
|
|
*/
|
|
XML_Node& speciesList = phaseNode.child("speciesArray");
|
|
XML_Node* speciesDB =
|
|
get_XML_NameID("speciesData", speciesList["datasrc"],
|
|
&phaseNode.root());
|
|
const vector<string>&sss = speciesNames();
|
|
|
|
for (k = 0; k < m_kk; k++) {
|
|
XML_Node* s = speciesDB->findByAttr("name", sss[k]);
|
|
if (!s) {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Species Data Base " + sss[k] + " not found");
|
|
}
|
|
XML_Node *ss = s->findByName("standardState");
|
|
if (!ss) {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Species " + sss[k] +
|
|
" standardState XML block not found");
|
|
}
|
|
std::string modelStringa = ss->attrib("model");
|
|
if (modelStringa == "") {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Species " + sss[k] +
|
|
" standardState XML block model attribute not found");
|
|
}
|
|
std::string modelString = lowercase(modelStringa);
|
|
|
|
if (k == 0) {
|
|
if (modelString == "wateriapws" || modelString == "real_water" ||
|
|
modelString == "waterpdss") {
|
|
/*
|
|
* Initialize the water standard state model
|
|
*/
|
|
if (m_waterSS) delete m_waterSS;
|
|
m_waterSS = new WaterPDSS(this, 0);
|
|
/*
|
|
* Fill in the molar volume of water (m3/kmol)
|
|
* at standard conditions to fill in the m_speciesSize entry
|
|
* with something reasonable.
|
|
*/
|
|
m_waterSS->setState_TP(300., OneAtm);
|
|
double dens = m_waterSS->density();
|
|
double mw = m_waterSS->molecularWeight();
|
|
m_speciesSize[0] = mw / dens;
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "Solvent species " << sss[k] << " has volume " <<
|
|
m_speciesSize[k] << endl;
|
|
#endif
|
|
} else if (modelString == "constant_incompressible") {
|
|
m_speciesSize[k] = getFloat(*ss, "molarVolume", "-");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "species " << sss[k] << " has volume " <<
|
|
m_speciesSize[k] << endl;
|
|
#endif
|
|
} else {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Solvent SS Model \"" + modelStringa +
|
|
"\" is not known");
|
|
}
|
|
} else {
|
|
if (modelString != "constant_incompressible") {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"Solute SS Model \"" + modelStringa +
|
|
"\" is not known");
|
|
}
|
|
m_speciesSize[k] = getFloat(*ss, "molarVolume", "-");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "species " << sss[k] << " has volume " <<
|
|
m_speciesSize[k] << endl;
|
|
#endif
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* Go get all of the coefficients and factors in the
|
|
* activityCoefficients XML block
|
|
*/
|
|
XML_Node *acNodePtr = 0;
|
|
if (thermoNode.hasChild("activityCoefficients")) {
|
|
XML_Node& acNode = thermoNode.child("activityCoefficients");
|
|
acNodePtr = &acNode;
|
|
/*
|
|
* Look for parameters for A_Debye
|
|
*/
|
|
if (acNode.hasChild("A_Debye")) {
|
|
XML_Node *ss = acNode.findByName("A_Debye");
|
|
string modelStringa = ss->attrib("model");
|
|
string modelString = lowercase(modelStringa);
|
|
if (modelString != "") {
|
|
if (modelString == "water") {
|
|
m_form_A_Debye = A_DEBYE_WATER;
|
|
} else {
|
|
throw CanteraError("DebyeHuckel::initThermoXML",
|
|
"A_Debye Model \"" + modelStringa +
|
|
"\" is not known");
|
|
}
|
|
} else {
|
|
m_A_Debye = getFloat(acNode, "A_Debye");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "A_Debye = " << m_A_Debye << endl;
|
|
#endif
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Initialize the water property calculator. It will share
|
|
* the internal eos water calculator.
|
|
*/
|
|
if (m_form_A_Debye == A_DEBYE_WATER) {
|
|
if (m_waterProps) delete m_waterProps;
|
|
m_waterProps = new WaterProps(m_waterSS);
|
|
}
|
|
|
|
/*
|
|
* Look for parameters for B_Debye
|
|
*/
|
|
if (acNode.hasChild("B_Debye")) {
|
|
m_B_Debye = getFloat(acNode, "B_Debye");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "B_Debye = " << m_B_Debye << endl;
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
* Look for parameters for B_dot
|
|
*/
|
|
if (acNode.hasChild("B_dot")) {
|
|
if (m_formDH == DHFORM_BETAIJ ||
|
|
m_formDH == DHFORM_DILUTE_LIMIT ||
|
|
m_formDH == DHFORM_PITZER_BETAIJ) {
|
|
throw CanteraError("DebyeHuckel:init",
|
|
"B_dot entry in the wrong DH form");
|
|
}
|
|
double bdot_common = getFloat(acNode, "B_dot");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "B_dot = " << bdot_common << endl;
|
|
#endif
|
|
/*
|
|
* Set B_dot parameters for charged species
|
|
*/
|
|
for (int k = 0; k < m_kk; k++) {
|
|
double z_k = charge(k);
|
|
if (fabs (z_k) > 0.0001) {
|
|
m_B_Dot[k] = bdot_common;
|
|
} else {
|
|
m_B_Dot[k] = 0.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Look for Parameters for the Maximum Ionic Strength
|
|
*/
|
|
if (acNode.hasChild("maxIonicStrength")) {
|
|
m_maxIionicStrength = getFloat(acNode, "maxIonicStrength");
|
|
#ifdef DEBUG_HKM_NOT
|
|
cout << "m_maxIionicStrength = "
|
|
<<m_maxIionicStrength << endl;
|
|
#endif
|
|
}
|
|
|
|
/*
|
|
* Look for Helgeson Parameters
|
|
*/
|
|
if (acNode.hasChild("UseHelgesonFixedForm")) {
|
|
m_useHelgesonFixedForm = true;
|
|
} else {
|
|
m_useHelgesonFixedForm = false;
|
|
}
|
|
|
|
/*
|
|
* Look for parameters for the Ionic radius
|
|
*/
|
|
if (acNode.hasChild("ionicRadius")) {
|
|
XML_Node& irNode = acNode.child("ionicRadius");
|
|
|
|
std::string Aunits = "";
|
|
double Afactor = 1.0;
|
|
if (irNode.hasAttrib("units")) {
|
|
std::string Aunits = irNode.attrib("units");
|
|
Afactor = toSI(Aunits);
|
|
}
|
|
|
|
if (irNode.hasAttrib("default")) {
|
|
std::string ads = irNode.attrib("default");
|
|
double ad = fpValue(ads);
|
|
for (int k = 0; k < m_kk; k++) {
|
|
m_Aionic[k] = ad * Afactor;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* If the Debye-Huckel form is BDOT_AK, we can
|
|
* have separate values for the denominator's ionic
|
|
* size. -> That's how the activity coefficient is
|
|
* parameterized. In this case only do we allow the
|
|
* code to read in these parameters.
|
|
*/
|
|
if (m_formDH == DHFORM_BDOT_AK) {
|
|
/*
|
|
* Define a string-string map, and interpret the
|
|
* value of the xml element as binary pairs separated
|
|
* by colons, e.g.:
|
|
* Na+:3.0
|
|
* Cl-:4.0
|
|
* H+:9.0
|
|
* OH-:3.5
|
|
* Read them into the map.
|
|
*/
|
|
map<string, string> m;
|
|
getMap(irNode, m);
|
|
/*
|
|
* Iterate over the map pairs, interpreting the
|
|
* first string as a species in the current phase.
|
|
* If no match is made, silently ignore the
|
|
* lack of agreement (HKM -> may be changed in the
|
|
* future).
|
|
*/
|
|
map<std::string,std::string>::const_iterator _b = m.begin();
|
|
for (; _b != m.end(); ++_b) {
|
|
int kk = speciesIndex(_b->first);
|
|
if (kk < 0) {
|
|
//throw CanteraError(
|
|
// "DebyeHuckel::initThermoXML error",
|
|
// "no species match was found"
|
|
// );
|
|
} else {
|
|
m_Aionic[kk] = fpValue(_b->second) * Afactor;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Get the matrix of coefficients for the Beta
|
|
* binary interaction parameters. We assume here that
|
|
* this matrix is symmetric, so that we only have to
|
|
* input 1/2 of the values.
|
|
*/
|
|
if (acNode.hasChild("DHBetaMatrix")) {
|
|
if (m_formDH == DHFORM_BETAIJ ||
|
|
m_formDH == DHFORM_PITZER_BETAIJ) {
|
|
XML_Node& irNode = acNode.child("DHBetaMatrix");
|
|
const vector<string>& sn = speciesNames();
|
|
getMatrixValues(irNode, sn, sn, m_Beta_ij, true, true);
|
|
} else {
|
|
throw CanteraError("DebyeHuckel::initThermoXML:",
|
|
"DHBetaMatrix found for wrong type");
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Fill in parameters for the calculation of the
|
|
* stoichiometric Ionic Strength
|
|
*
|
|
* The default is that stoich charge is the same as the
|
|
* regular charge.
|
|
*/
|
|
m_speciesCharge_Stoich.resize(m_kk, 0.0);
|
|
for (k = 0; k < m_kk; k++) {
|
|
m_speciesCharge_Stoich[k] = m_speciesCharge[k];
|
|
}
|
|
/*
|
|
* First look at the species database.
|
|
* -> Look for the subelement "stoichIsMods"
|
|
* in each of the species SS databases.
|
|
*/
|
|
const XML_Node *phaseSpecies = speciesData();
|
|
if (phaseSpecies) {
|
|
std::string kname, jname;
|
|
vector<XML_Node*> xspecies;
|
|
phaseSpecies->getChildren("species",xspecies);
|
|
int jj = xspecies.size();
|
|
for (k = 0; k < m_kk; k++) {
|
|
int jmap = -1;
|
|
kname = speciesName(k);
|
|
for (int j = 0; j < jj; j++) {
|
|
const XML_Node& sp = *xspecies[j];
|
|
jname = sp["name"];
|
|
if (jname == kname) {
|
|
jmap = j;
|
|
break;
|
|
}
|
|
}
|
|
if (jmap > -1) {
|
|
const XML_Node& sp = *xspecies[jmap];
|
|
if (sp.hasChild("stoichIsMods")) {
|
|
double val = getFloat(sp, "stoichIsMods");
|
|
m_speciesCharge_Stoich[k] = val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Now look at the activity coefficient database
|
|
*/
|
|
if (acNodePtr) {
|
|
if (acNodePtr->hasChild("stoichIsMods")) {
|
|
XML_Node& sIsNode = acNodePtr->child("stoichIsMods");
|
|
|
|
map<std::string, std::string> msIs;
|
|
getMap(sIsNode, msIs);
|
|
map<std::string,std::string>::const_iterator _b = msIs.begin();
|
|
for (; _b != msIs.end(); ++_b) {
|
|
int kk = speciesIndex(_b->first);
|
|
if (kk < 0) {
|
|
//throw CanteraError(
|
|
// "DebyeHuckel::initThermoXML error",
|
|
// "no species match was found"
|
|
// );
|
|
} else {
|
|
double val = fpValue(_b->second);
|
|
m_speciesCharge_Stoich[kk] = val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Fill in the vector specifying the electrolyte species
|
|
* type
|
|
*
|
|
* First fill in default values. Everthing is either
|
|
* a charge species, a nonpolar neutral, or the solvent.
|
|
*/
|
|
for (k = 0; k < m_kk; k++) {
|
|
if (fabs(m_speciesCharge[k]) > 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_chargedSpecies;
|
|
if (fabs(m_speciesCharge_Stoich[k] - m_speciesCharge[k])
|
|
> 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
|
|
}
|
|
} else if (fabs(m_speciesCharge_Stoich[k]) > 0.0001) {
|
|
m_electrolyteSpeciesType[k] = cEST_weakAcidAssociated;
|
|
} else {
|
|
m_electrolyteSpeciesType[k] = cEST_nonpolarNeutral;
|
|
}
|
|
}
|
|
m_electrolyteSpeciesType[m_indexSolvent] = cEST_solvent;
|
|
/*
|
|
* First look at the species database.
|
|
* -> Look for the subelement "stoichIsMods"
|
|
* in each of the species SS databases.
|
|
*/
|
|
const XML_Node *phaseSpecies = speciesData();
|
|
const XML_Node *spPtr = 0;
|
|
if (phaseSpecies) {
|
|
std::string kname;
|
|
for (k = 0; k < m_kk; k++) {
|
|
kname = speciesName(k);
|
|
spPtr = speciesXML_Node(kname, phaseSpecies);
|
|
if (!spPtr) {
|
|
if (spPtr->hasChild("electrolyteSpeciesType")) {
|
|
std::string est = getString(*spPtr, "electrolyteSpeciesType");
|
|
if ((m_electrolyteSpeciesType[k] = interp_est(est)) == -1) {
|
|
throw CanteraError("DebyeHuckel:initThermoXML",
|
|
"Bad electrolyte type: " + est);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
* Then look at the phase thermo specification
|
|
*/
|
|
if (acNodePtr) {
|
|
if (acNodePtr->hasChild("electrolyteSpeciesType")) {
|
|
XML_Node& ESTNode = acNodePtr->child("electrolyteSpeciesType");
|
|
map<std::string, std::string> msEST;
|
|
getMap(ESTNode, msEST);
|
|
map<std::string,std::string>::const_iterator _b = msEST.begin();
|
|
for (; _b != msEST.end(); ++_b) {
|
|
int kk = speciesIndex(_b->first);
|
|
if (kk < 0) {
|
|
} else {
|
|
std::string est = _b->second;
|
|
if ((m_electrolyteSpeciesType[kk] = interp_est(est)) == -1) {
|
|
throw CanteraError("DebyeHuckel:initThermoXML",
|
|
"Bad electrolyte type: " + est);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Lastly set the state
|
|
*/
|
|
if (phaseNode.hasChild("state")) {
|
|
XML_Node& stateNode = phaseNode.child("state");
|
|
setStateFromXML(stateNode);
|
|
}
|
|
|
|
}
|
|
|
|
/*
|
|
* @internal
|
|
* Set equation of state parameters. The number and meaning of
|
|
* these depends on the subclass.
|
|
* @param n number of parameters
|
|
* @param c array of \i n coefficients
|
|
*
|
|
*/
|
|
void DebyeHuckel::setParameters(int n, doublereal* c) {
|
|
}
|
|
void DebyeHuckel::getParameters(int &n, doublereal * const c) {
|
|
}
|
|
|
|
/*
|
|
* Set equation of state parameter values from XML
|
|
* entries. This method is called by function importPhase in
|
|
* file importCTML.cpp when processing a phase definition in
|
|
* an input file. It should be overloaded in subclasses to set
|
|
* any parameters that are specific to that particular phase
|
|
* model.
|
|
*
|
|
* @param eosdata An XML_Node object corresponding to
|
|
* the "thermo" entry for this phase in the input file.
|
|
*
|
|
* HKM -> Right now, the parameters are set elsewhere (initThermoXML)
|
|
* It just didn't seem to fit.
|
|
*/
|
|
void DebyeHuckel::setParametersFromXML(const XML_Node& eosdata) {
|
|
}
|
|
|
|
/*
|
|
* Report the molar volume of species k
|
|
*
|
|
* units - \f$ m^3 kmol^-1 \f$
|
|
*/
|
|
// double DebyeHuckel::speciesMolarVolume(int k) const {
|
|
// return m_speciesSize[k];
|
|
//}
|
|
|
|
|
|
/*
|
|
* A_Debye_TP() (virtual)
|
|
*
|
|
* Returns the A_Debye parameter as a function of temperature
|
|
* and pressure.
|
|
*
|
|
* The default is to assume that it is constant, given
|
|
* in the initialization process and storred in the
|
|
* member double, m_A_Debye
|
|
*/
|
|
double DebyeHuckel::A_Debye_TP(double tempArg, double presArg) const {
|
|
double T = temperature();
|
|
double A;
|
|
if (tempArg != -1.0) {
|
|
T = tempArg;
|
|
}
|
|
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;
|
|
}
|
|
|
|
/*
|
|
* dA_DebyedT_TP() (virtual)
|
|
*
|
|
* Returns the derivative of the A_Debye parameter with
|
|
* respect to temperature as a function of temperature
|
|
* and pressure.
|
|
*
|
|
* 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();
|
|
if (tempArg != -1.0) {
|
|
T = tempArg;
|
|
}
|
|
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;
|
|
}
|
|
|
|
/*
|
|
* d2A_DebyedT2_TP() (virtual)
|
|
*
|
|
* Returns the 2nd derivative of the A_Debye parameter with
|
|
* respect to temperature as a function of temperature
|
|
* 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();
|
|
if (tempArg != -1.0) {
|
|
T = tempArg;
|
|
}
|
|
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;
|
|
}
|
|
|
|
/*
|
|
* ----------- Critical State Properties --------------------------
|
|
*/
|
|
|
|
/*
|
|
* ---------- Other Property Functions
|
|
*/
|
|
double DebyeHuckel::AionicRadius(int k) const {
|
|
return m_Aionic[k];
|
|
}
|
|
|
|
/*
|
|
* ------------ Private and Restricted Functions ------------------
|
|
*/
|
|
|
|
/**
|
|
* Bail out of functions with an error exit if they are not
|
|
* implemented.
|
|
*/
|
|
doublereal DebyeHuckel::err(std::string msg) const {
|
|
throw CanteraError("DebyeHuckel",
|
|
"Unfinished func called: " + msg );
|
|
return 0.0;
|
|
}
|
|
|
|
|
|
/*
|
|
* initLengths():
|
|
*
|
|
* This internal function adjusts the lengths of arrays based on
|
|
* the number of species
|
|
*/
|
|
void DebyeHuckel::initLengths() {
|
|
m_kk = nSpecies();
|
|
|
|
/*
|
|
* Obtain the limits of the temperature from the species
|
|
* thermo handler's limits.
|
|
*/
|
|
int leng = m_kk;
|
|
m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
|
|
m_speciesSize.resize(leng);
|
|
m_Aionic.resize(leng, 0.0);
|
|
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);
|
|
m_pp.resize(leng, 0.0);
|
|
m_tmpV.resize(leng, 0.0);
|
|
if (m_formDH == DHFORM_BETAIJ ||
|
|
m_formDH == DHFORM_PITZER_BETAIJ) {
|
|
m_Beta_ij.resize(leng, leng, 0.0);
|
|
}
|
|
}
|
|
|
|
/**
|
|
* nonpolarActCoeff() (private)
|
|
*
|
|
* Static function that implements the non-polar species
|
|
* salt-out modifications.
|
|
* Returns the calculated activity coefficients.
|
|
*/
|
|
double DebyeHuckel::_nonpolarActCoeff(double IionicMolality) const {
|
|
double I2 = IionicMolality * IionicMolality;
|
|
double l10actCoeff =
|
|
m_npActCoeff[0] * IionicMolality +
|
|
m_npActCoeff[1] * I2 +
|
|
m_npActCoeff[2] * I2 * IionicMolality;
|
|
return pow(10.0 , l10actCoeff);
|
|
}
|
|
|
|
|
|
/**
|
|
* _osmoticCoeffHelgesonFixedForm()
|
|
*
|
|
* Formula for the osmotic coefficient that occurs in
|
|
* the GWB. It is originally from Helgeson for a variable
|
|
* NaCl brine. It's to be used with extreme caution.
|
|
*/
|
|
double DebyeHuckel::
|
|
_osmoticCoeffHelgesonFixedForm() const {
|
|
const double a0 = 1.454;
|
|
const double b0 = 0.02236;
|
|
const double c0 = 9.380E-3;
|
|
const double d0 = -5.362E-4;
|
|
double Is = m_IionicMolalityStoich;
|
|
if (Is <= 0.0) {
|
|
return 0.0;
|
|
}
|
|
double Is2 = Is * Is;
|
|
double bhat = 1.0 + a0 * sqrt(Is);
|
|
double func = bhat - 2.0 * log(bhat) - 1.0/bhat;
|
|
double v1 = m_A_Debye / (a0 * a0 * a0 * Is) * func;
|
|
double oc = 1.0 - v1 + b0 * Is / 2.0 + 2.0 * c0 * Is2 / 3.0
|
|
+ 3.0 * d0 * Is2 * Is / 4.0;
|
|
return oc;
|
|
}
|
|
|
|
|
|
/**
|
|
* _activityWaterHelgesonFixedForm()
|
|
*
|
|
* Formula for the log of the activity of the water
|
|
* solvent that occurs in
|
|
* the GWB. It is originally from Helgeson for a variable
|
|
* NaCl brine. It's to be used with extreme caution.
|
|
*/
|
|
double DebyeHuckel::
|
|
_lnactivityWaterHelgesonFixedForm() const {
|
|
/*
|
|
* Update the internally storred vector of molalities
|
|
*/
|
|
calcMolalities();
|
|
double oc = _osmoticCoeffHelgesonFixedForm();
|
|
double sum = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent) {
|
|
sum += MAX(m_molalities[k], 0.0);
|
|
}
|
|
}
|
|
if (sum > 2.0 * m_maxIionicStrength) {
|
|
sum = 2.0 * m_maxIionicStrength;
|
|
};
|
|
double lac = - m_Mnaught * sum * oc;
|
|
return lac;
|
|
}
|
|
|
|
/*
|
|
* s_update_lnMolalityActCoeff():
|
|
*
|
|
* Using internally stored values, this function calculates
|
|
* the activity coefficients for all species.
|
|
*
|
|
* The ln(activity_solvent) is first calculated for the
|
|
* solvent. Then the molar based activity coefficient
|
|
* is calculated and returned.
|
|
*
|
|
* ( Note this is the main routine for implementing the
|
|
* activity coefficient formulation.)
|
|
*/
|
|
void DebyeHuckel::s_update_lnMolalityActCoeff() const {
|
|
double z_k, zs_k1, zs_k2;
|
|
/*
|
|
* Update the internally storred vector of molalities
|
|
*/
|
|
calcMolalities();
|
|
/*
|
|
* Calculate the apparent (real) ionic strength.
|
|
*
|
|
* Note this is not the stoichiometric ionic strengh,
|
|
* where reactions of ions forming neutral salts
|
|
* are ignorred in calculating the ionic strength.
|
|
*/
|
|
m_IionicMolality = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
m_IionicMolality += m_molalities[k] * z_k * z_k;
|
|
}
|
|
m_IionicMolality /= 2.0;
|
|
|
|
if (m_IionicMolality > m_maxIionicStrength) {
|
|
m_IionicMolality = m_maxIionicStrength;
|
|
}
|
|
|
|
/*
|
|
* Calculate the stoichiometric ionic charge
|
|
*/
|
|
m_IionicMolalityStoich = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
zs_k1 = m_speciesCharge_Stoich[k];
|
|
if (z_k == zs_k1) {
|
|
m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
|
|
} else {
|
|
zs_k2 = z_k - zs_k1;
|
|
m_IionicMolalityStoich
|
|
+= m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
|
|
}
|
|
}
|
|
m_IionicMolalityStoich /= 2.0;
|
|
|
|
if (m_IionicMolalityStoich > m_maxIionicStrength) {
|
|
m_IionicMolalityStoich = m_maxIionicStrength;
|
|
}
|
|
|
|
/*
|
|
* Possibly update the storred value of the
|
|
* Debye-Huckel parameter A_Debye
|
|
* This parameter appears on the top of the activity
|
|
* coefficient expression.
|
|
* It depends on T (and P), as it depends explicity
|
|
* on the temperature. Also, the dielectric constant
|
|
* is usually a fairly strong function of T, also.
|
|
*/
|
|
m_A_Debye = A_Debye_TP();
|
|
|
|
/*
|
|
* Calculate a safe value for the mole fraction
|
|
* of the solvent
|
|
*/
|
|
double xmolSolvent = moleFraction(m_indexSolvent);
|
|
xmolSolvent = MAX(8.689E-3, xmolSolvent);
|
|
|
|
int est;
|
|
double ac_nonPolar = 1.0;
|
|
double numTmp = m_A_Debye * sqrt(m_IionicMolality);
|
|
double denomTmp = m_B_Debye * sqrt(m_IionicMolality);
|
|
double coeff;
|
|
double lnActivitySolvent = 0.0;
|
|
double tmp;
|
|
double tmpLn;
|
|
double y, yp1, sigma;
|
|
switch (m_formDH) {
|
|
case DHFORM_DILUTE_LIMIT:
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
m_lnActCoeffMolal[k] = - z_k * z_k * numTmp;
|
|
}
|
|
lnActivitySolvent =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 / 3.0 * m_A_Debye * m_Mnaught *
|
|
m_IionicMolality * sqrt(m_IionicMolality);
|
|
break;
|
|
|
|
case DHFORM_BDOT_AK:
|
|
ac_nonPolar = _nonpolarActCoeff(m_IionicMolality);
|
|
for (int k = 0; k < m_kk; k++) {
|
|
est = m_electrolyteSpeciesType[k];
|
|
if (est == cEST_nonpolarNeutral) {
|
|
m_lnActCoeffMolal[k] = log(ac_nonPolar);
|
|
} else {
|
|
z_k = m_speciesCharge[k];
|
|
m_lnActCoeffMolal[k] =
|
|
- z_k * z_k * numTmp / (1.0 + denomTmp * m_Aionic[k])
|
|
+ log(10.0) * m_B_Dot[k] * m_IionicMolality;
|
|
}
|
|
}
|
|
|
|
lnActivitySolvent = (xmolSolvent - 1.0)/xmolSolvent;
|
|
coeff = 2.0 / 3.0 * m_A_Debye * m_Mnaught
|
|
* sqrt(m_IionicMolality);
|
|
tmp = 0.0;
|
|
if (denomTmp > 0.0) {
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent || m_Aionic[k] != 0.0) {
|
|
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;
|
|
}
|
|
}
|
|
}
|
|
lnActivitySolvent += coeff * tmp;
|
|
tmp = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
if ((k != m_indexSolvent) && (z_k != 0.0)) {
|
|
tmp += m_B_Dot[k] * m_molalities[k];
|
|
}
|
|
}
|
|
lnActivitySolvent -=
|
|
m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
|
|
|
|
/*
|
|
* Special section to implement the Helgeson fixed form
|
|
* for the water brine activity coefficient.
|
|
*/
|
|
if (m_useHelgesonFixedForm) {
|
|
lnActivitySolvent = _lnactivityWaterHelgesonFixedForm();
|
|
}
|
|
break;
|
|
|
|
case DHFORM_BDOT_ACOMMON:
|
|
denomTmp *= m_Aionic[0];
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
m_lnActCoeffMolal[k] =
|
|
- z_k * z_k * numTmp / (1.0 + denomTmp)
|
|
+ log(10.0) * m_B_Dot[k] * m_IionicMolality;
|
|
}
|
|
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;
|
|
}
|
|
lnActivitySolvent =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * m_A_Debye * m_Mnaught *
|
|
m_IionicMolality * sqrt(m_IionicMolality) * sigma;
|
|
tmp = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
if ((k != m_indexSolvent) && (z_k != 0.0)) {
|
|
tmp += m_B_Dot[k] * m_molalities[k];
|
|
}
|
|
}
|
|
lnActivitySolvent -=
|
|
m_Mnaught * log(10.0) * m_IionicMolality * tmp / 2.0;
|
|
|
|
break;
|
|
|
|
case DHFORM_BETAIJ:
|
|
denomTmp = m_B_Debye * m_Aionic[0];
|
|
denomTmp *= sqrt(m_IionicMolality);
|
|
lnActivitySolvent =
|
|
(xmolSolvent - 1.0)/xmolSolvent;
|
|
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent) {
|
|
z_k = m_speciesCharge[k];
|
|
m_lnActCoeffMolal[k] =
|
|
- z_k * z_k * numTmp / (1.0 + denomTmp);
|
|
for (int j = 0; j < m_kk; j++) {
|
|
double beta = m_Beta_ij.value(k, j);
|
|
#ifdef DEBUG_HKM_NOT
|
|
if (beta != 0.0) {
|
|
printf("b: k = %d, j = %d, betakj = %g\n",
|
|
k, j, beta);
|
|
}
|
|
#endif
|
|
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] * beta;
|
|
}
|
|
}
|
|
}
|
|
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;
|
|
}
|
|
lnActivitySolvent =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * m_A_Debye * m_Mnaught *
|
|
m_IionicMolality * sqrt(m_IionicMolality) * sigma;
|
|
tmp = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
for (int j = 0; j < m_kk; j++) {
|
|
tmp +=
|
|
m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
|
|
}
|
|
}
|
|
lnActivitySolvent -= m_Mnaught * tmp;
|
|
break;
|
|
|
|
case DHFORM_PITZER_BETAIJ:
|
|
denomTmp = m_B_Debye * sqrt(m_IionicMolality);
|
|
denomTmp *= m_Aionic[0];
|
|
numTmp = m_A_Debye * sqrt(m_IionicMolality);
|
|
tmpLn = log(1.0 + denomTmp);
|
|
for (int k = 0; k < m_kk; k++) {
|
|
if (k != m_indexSolvent) {
|
|
z_k = m_speciesCharge[k];
|
|
m_lnActCoeffMolal[k] =
|
|
- z_k * z_k * numTmp / 3.0 / (1.0 + denomTmp);
|
|
m_lnActCoeffMolal[k] +=
|
|
- 2.0 * z_k * z_k * m_A_Debye * tmpLn /
|
|
(3.0 * m_B_Debye * m_Aionic[0]);
|
|
for (int j = 0; j < m_kk; j++) {
|
|
m_lnActCoeffMolal[k] += 2.0 * m_molalities[j] *
|
|
m_Beta_ij.value(k, j);
|
|
}
|
|
}
|
|
}
|
|
sigma = 1.0 / (1.0 + denomTmp);
|
|
lnActivitySolvent =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * m_A_Debye * m_Mnaught *
|
|
m_IionicMolality * sqrt(m_IionicMolality) * sigma;
|
|
tmp = 0.0;
|
|
for (int k = 0; k < m_kk; k++) {
|
|
for (int j = 0; j < m_kk; j++) {
|
|
tmp +=
|
|
m_Beta_ij.value(k, j) * m_molalities[k] * m_molalities[j];
|
|
}
|
|
}
|
|
lnActivitySolvent -= m_Mnaught * tmp;
|
|
break;
|
|
|
|
default:
|
|
printf("ERROR\n");
|
|
exit(-1);
|
|
}
|
|
/*
|
|
* Above, we calculated the ln(activitySolvent). Translate that
|
|
* into the molar-based activity coefficient by dividing by
|
|
* the solvent mole fraction. Solvents are not on the molality
|
|
* scale.
|
|
*/
|
|
xmolSolvent = moleFraction(m_indexSolvent);
|
|
m_lnActCoeffMolal[m_indexSolvent] =
|
|
lnActivitySolvent - log(xmolSolvent);
|
|
}
|
|
|
|
/*
|
|
* s_update_dMolalityActCoeff_dT() (private, const )
|
|
*
|
|
* Using internally stored values, this function calculates
|
|
* the temperature derivative of the logarithm of the
|
|
* activity coefficient for all species in the mechanism.
|
|
*
|
|
* We assume that the activity coefficients are current.
|
|
*
|
|
* solvent activity coefficient is on the molality
|
|
* scale. It's derivative is too.
|
|
*/
|
|
void DebyeHuckel::s_update_dlnMolalityActCoeff_dT() const {
|
|
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
|
|
int k;
|
|
double dAdT = dA_DebyedT_TP();
|
|
if (dAdT == 0.0) {
|
|
for (k = 0; k < m_kk; k++) {
|
|
m_dlnActCoeffMolaldT[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 numdAdTTmp = dAdT * sqrtI;
|
|
double denomTmp = m_B_Debye * sqrtI;
|
|
|
|
switch (m_formDH) {
|
|
case DHFORM_DILUTE_LIMIT:
|
|
for (int k = 0; k < m_kk; k++) {
|
|
m_dlnActCoeffMolaldT[k] =
|
|
m_lnActCoeffMolal[k] * dAdT / m_A_Debye;
|
|
}
|
|
break;
|
|
|
|
case DHFORM_BDOT_AK:
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
m_dlnActCoeffMolaldT[k] =
|
|
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp * m_Aionic[k]);
|
|
}
|
|
|
|
m_dlnActCoeffMolaldT[m_indexSolvent] = 0.0;
|
|
|
|
coeff = 2.0 / 3.0 * dAdT * 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_dlnActCoeffMolaldT[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_dlnActCoeffMolaldT[k] =
|
|
- z_k * z_k * numdAdTTmp / (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_dlnActCoeffMolaldT[m_indexSolvent] =
|
|
2.0 /3.0 * dAdT * 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_dlnActCoeffMolaldT[k] =
|
|
- z_k * z_k * numdAdTTmp / (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_dlnActCoeffMolaldT[m_indexSolvent] =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * dAdT * 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_dlnActCoeffMolaldT[k] =
|
|
- z_k * z_k * numdAdTTmp / (1.0 + denomTmp)
|
|
- 2.0 * z_k * z_k * dAdT * tmpLn
|
|
/ (m_B_Debye * m_Aionic[0]);
|
|
m_dlnActCoeffMolaldT[k] /= 3.0;
|
|
}
|
|
}
|
|
|
|
sigma = 1.0 / ( 1.0 + denomTmp);
|
|
m_dlnActCoeffMolaldT[m_indexSolvent] =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * dAdT * m_Mnaught *
|
|
m_IionicMolality * sqrtI * sigma;
|
|
break;
|
|
|
|
default:
|
|
printf("ERROR\n");
|
|
exit(-1);
|
|
break;
|
|
}
|
|
}
|
|
|
|
/**
|
|
* s_update_d2lnMolalityActCoeff_dT2() (private, const )
|
|
*
|
|
* Using internally stored values, this function calculates
|
|
* the temperature 2nd derivative of the logarithm of the
|
|
* activity coefficient
|
|
* for all species in the mechanism.
|
|
*
|
|
* We assume that the activity coefficients are current.
|
|
*
|
|
* solvent activity coefficient is on the molality
|
|
* scale. It's derivatives are too.
|
|
*/
|
|
void DebyeHuckel::s_update_d2lnMolalityActCoeff_dT2() const {
|
|
double z_k, coeff, tmp, y, yp1, sigma, tmpLn;
|
|
int k;
|
|
double dAdT = dA_DebyedT_TP();
|
|
double d2AdT2 = d2A_DebyedT2_TP();
|
|
if (d2AdT2 == 0.0 && dAdT == 0.0) {
|
|
for (k = 0; k < m_kk; k++) {
|
|
m_d2lnActCoeffMolaldT2[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 numd2AdT2Tmp = d2AdT2 * sqrtI;
|
|
double denomTmp = m_B_Debye * sqrtI;
|
|
|
|
switch (m_formDH) {
|
|
case DHFORM_DILUTE_LIMIT:
|
|
for (int k = 0; k < m_kk; k++) {
|
|
m_d2lnActCoeffMolaldT2[k] =
|
|
m_lnActCoeffMolal[k] * d2AdT2 / m_A_Debye;
|
|
}
|
|
break;
|
|
|
|
case DHFORM_BDOT_AK:
|
|
for (int k = 0; k < m_kk; k++) {
|
|
z_k = m_speciesCharge[k];
|
|
m_d2lnActCoeffMolaldT2[k] =
|
|
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp * m_Aionic[k]);
|
|
}
|
|
|
|
m_d2lnActCoeffMolaldT2[m_indexSolvent] = 0.0;
|
|
|
|
coeff = 2.0 / 3.0 * d2AdT2 * 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_d2lnActCoeffMolaldT2[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_d2lnActCoeffMolaldT2[k] =
|
|
- z_k * z_k * numd2AdT2Tmp / (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_d2lnActCoeffMolaldT2[m_indexSolvent] =
|
|
2.0 /3.0 * d2AdT2 * 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_d2lnActCoeffMolaldT2[k] =
|
|
- z_k * z_k * numd2AdT2Tmp / (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_d2lnActCoeffMolaldT2[m_indexSolvent] =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * d2AdT2 * 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_dlnActCoeffMolaldT[k] =
|
|
- z_k * z_k * numd2AdT2Tmp / (1.0 + denomTmp)
|
|
- 2.0 * z_k * z_k * d2AdT2 * tmpLn
|
|
/ (m_B_Debye * m_Aionic[0]);
|
|
m_dlnActCoeffMolaldT[k] /= 3.0;
|
|
}
|
|
}
|
|
|
|
sigma = 1.0 / ( 1.0 + denomTmp);
|
|
m_dlnActCoeffMolaldT[m_indexSolvent] =
|
|
(xmolSolvent - 1.0)/xmolSolvent +
|
|
2.0 /3.0 * d2AdT2 * m_Mnaught *
|
|
m_IionicMolality * sqrtI * sigma;
|
|
break;
|
|
|
|
default:
|
|
printf("ERROR\n");
|
|
exit(-1);
|
|
break;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* 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;
|
|
}
|
|
}
|
|
|
|
/*
|
|
* Updates the standard state thermodynamic functions at the current T and P of the solution.
|
|
*
|
|
* @internal
|
|
*
|
|
* This function gets called for every call to functions in this
|
|
* class. It checks to see whether the temperature or pressure has changed and
|
|
* thus the ss thermodynamics functions for all of the species
|
|
* must be recalculated.
|
|
*
|
|
*
|
|
* Note, this will throw an error. It must be reimplemented in derived classes.
|
|
*/
|
|
void DebyeHuckel::_updateStandardStateThermo(doublereal pnow) const {
|
|
_updateRefStateThermo();
|
|
doublereal tnow = temperature();
|
|
if (pnow == -1.0) {
|
|
pnow = m_Pcurrent;
|
|
}
|
|
if (m_tlast != tnow || m_plast != pnow) {
|
|
if (m_waterSS) {
|
|
m_waterSS->setTempPressure(tnow, pnow);
|
|
}
|
|
m_tlast = tnow;
|
|
m_plast = pnow;
|
|
}
|
|
}
|
|
|
|
}
|