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

Unifying the water treatment between objects.

Browse source 
Found a bug in getPartialMolarEntropy in HMWSoln that wasn't covered
in the test suite.
This commit is contained in:
Harry Moffat 2006-07-12 15:41:38 +00:00
parent 60ff5fb20e
commit 07d5b6a6aa
4 changed files with 230 additions and 81 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -17,7 +17,7 @@
#include "DebyeHuckel.h"
#include "importCTML.h"
#include "WaterProps.h"
#include "WaterPDSS.h"
namespace Cantera {
@ -36,6 +36,8 @@ namespace Cantera {
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);
@ -63,6 +65,8 @@ namespace Cantera {
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);
@ -84,6 +88,8 @@ namespace Cantera {
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);
@ -133,6 +139,14 @@ namespace Cantera {
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;
@ -156,9 +170,16 @@ namespace Cantera {
/**
* ~DebyeHuckel(): (virtual)
*
* Destructor: does nothing:
* 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;
}
}
/**
@ -205,17 +226,29 @@ namespace Cantera {
* Molar enthalpy of the solution. Units: J/kmol.
*/
doublereal DebyeHuckel::enthalpy_mole() const {
return err("not implemented");
}
/**
* Molar internal energy of the solution. Units: J/kmol.
*/
doublereal DebyeHuckel::intEnergy_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));
@ -227,10 +260,16 @@ namespace Cantera {
return mean_X(DATA_PTR(m_tmpV));
}
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
/**
* 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));
return mean_X(DATA_PTR(m_tmpV));
double val = mean_X(DATA_PTR(m_tmpV));
return val;
}
/// Molar heat capacity at constant volume. Units: J/kmol/K.
@ -254,6 +293,69 @@ namespace Cantera {
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
double temp = temperature();
if (m_waterSS) {
/*
* Call the water SS and set it's internal state
*/
m_waterSS->setTempPressure(temp, p);
/*
* 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();
}
/*
* Store the current pressure
*/
m_Pcurrent = p;
/*
* 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
* store the denisty.
*/
State::setDensity(dd);
}
/**
* The isothermal compressibility. Units: 1/Pa.
* The isothermal compressibility is defined as
@ -265,6 +367,8 @@ namespace Cantera {
* doesn't change with pressure or temperature.
*/
doublereal DebyeHuckel::isothermalCompressibility() const {
throw CanteraError("DebyeHuckel::isothermalCompressibility",
"unimplemented");
return 0.0;
}
@ -280,6 +384,8 @@ namespace Cantera {
* doesn't change with pressure or temperature.
*/
doublereal DebyeHuckel::thermalExpansionCoeff() const {
throw CanteraError("DebyeHuckel::thermalExpansionCoeff",
"unimplemented");
return 0.0;
}
@ -523,8 +629,18 @@ namespace Cantera {
* 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
@ -538,7 +654,6 @@ namespace Cantera {
*/
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
double T = temperature();
double RTT = GasConstant * T * T;
for (int k = 0; k < m_kk; k++) {
hbar[k] -= RTT * m_dlnActCoeffMolaldT[k];
@ -583,18 +698,23 @@ namespace Cantera {
* 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();
doublereal R = GasConstant;
doublereal mm;
/*
* 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]);
@ -603,8 +723,7 @@ namespace Cantera {
}
double xmolSolvent = moleFraction(m_indexSolvent);
mm = fmaxx(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -=
R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
/*
* Check to see whether activity coefficients are temperature
* dependent. If they are, then calculate the their temperature
@ -626,21 +745,33 @@ namespace Cantera {
* 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 equal to the
* constant species molar volumes.
* 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
*
* So, if the activity coefficients depended on pressure this
* function would be nontrivial.
* = 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
@ -685,7 +816,6 @@ namespace Cantera {
}
}
/*
* -------- Properties of the Standard State of the Species
* in the Solution ------------------
@ -716,6 +846,9 @@ namespace Cantera {
delta_p = m_Pcurrent - pref;
mu[k] += delta_p * m_speciesSize[k];
}
if (m_waterSS) {
mu[0] = m_waterSS->gibbs_mole();
}
}
/**
@ -2318,8 +2451,7 @@ namespace Cantera {
*
* Using internally stored values, this function calculates
* the pressure derivative of the logarithm of the
* activity coefficient
* for all species in the mechanism.
* activity coefficient for all species in the mechanism.
*
* We assume that the activity coefficients, molalities,
* and A_Debye are current.

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

@ -182,7 +182,7 @@ namespace Cantera {
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
/*
* (HKM -> Bump up to Parent object)
*
*/
virtual doublereal cp_mole() const;
@ -216,10 +216,7 @@ namespace Cantera {
* This method sets a constant within the object.
* The mass density is not a function of pressure.
*/
virtual void setPressure(doublereal p) {
m_Pcurrent = p;
}
virtual void setPressure(doublereal p);
/**
* Calculate the density of the mixture using the partial
@ -1116,6 +1113,12 @@ namespace Cantera {
*/
array_fp m_npActCoeff;
/**
* Water standard state -> derived from the
* equation of state for water.
*/
WaterPDSS *m_waterSS;
double m_densWaterSS;
/**
* Pointer to the water property calculator

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

@ -141,10 +141,13 @@ namespace Cantera {
m_IionicMolalityStoich= b.m_IionicMolalityStoich;
m_form_A_Debye = b.m_form_A_Debye;
m_A_Debye = b.m_A_Debye;
if (!m_waterSS) {
m_waterSS = new WaterPDSS(this, 0);
if (m_waterSS) {
delete m_waterSS;
m_waterSS = 0;
}
if (b.m_waterSS) {
m_waterSS = new WaterPDSS(*(b.m_waterSS));
}
m_waterSS = b.m_waterSS;
m_densWaterSS = b.m_densWaterSS;
if (m_waterProps) {
delete m_waterProps;
@ -153,7 +156,6 @@ namespace Cantera {
if (b.m_waterProps) {
m_waterProps = new WaterProps(*(b.m_waterProps));
}
m_waterSS = b.m_waterSS;
m_expg0_RT = b.m_expg0_RT;
m_pe = b.m_pe;
m_pp = b.m_pp;
@ -366,8 +368,12 @@ namespace Cantera {
* Destructor: does nothing:
*/
HMWSoln::~HMWSoln() {
delete m_waterProps;
delete m_waterSS;
if (m_waterProps) {
delete m_waterProps; m_waterProps = 0;
}
if (m_waterSS) {
delete m_waterSS; m_waterSS = 0;
}
}
/**
@ -488,6 +494,9 @@ namespace Cantera {
/**
* Molar internal energy of the solution. Units: J/kmol.
*
* This is calculated from the soln enthalpy and then
* subtracting pV.
*/
doublereal HMWSoln::intEnergy_mole() const {
double hh = enthalpy_mole();
@ -497,6 +506,11 @@ namespace Cantera {
return uu;
}
/**
* Molar soln entropy at constant pressure. Units: J/kmol/K.
*
* This is calculated from the partial molar entropies.
*/
doublereal HMWSoln::entropy_mole() const {
getPartialMolarEntropies(DATA_PTR(m_tmpV));
return mean_X(DATA_PTR(m_tmpV));
@ -508,7 +522,11 @@ namespace Cantera {
return mean_X(DATA_PTR(m_tmpV));
}
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
/** 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 HMWSoln::cp_mole() const {
getPartialMolarCp(DATA_PTR(m_tmpV));
double val = mean_X(DATA_PTR(m_tmpV));
@ -594,7 +612,6 @@ namespace Cantera {
* store the denisty.
*/
State::setDensity(dd);
}
/**
@ -608,6 +625,8 @@ namespace Cantera {
* doesn't change with pressure or temperature.
*/
doublereal HMWSoln::isothermalCompressibility() const {
throw CanteraError("HMWSoln::isothermalCompressibility",
"unimplemented");
return 0.0;
}
@ -882,19 +901,24 @@ namespace Cantera {
* We calculate
*/
void HMWSoln::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;
}
/*
* Update the activity coefficients, This also update the
* internally storred molalities.
*/
s_update_lnMolalityActCoeff();
s_update_dlnMolalityActCoeff_dT();
double T = temperature();
double RT = GasConstant * T;
for (int k = 0; k < m_kk; k++) {
hbar[k] *= RT;
}
double RTT = RT * T;
for (int k = 0; k < m_kk; k++) {
hbar[k] -= RTT * m_dlnActCoeffMolaldT[k];
@ -935,18 +959,23 @@ namespace Cantera {
* 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 stored molalities.
*/
s_update_lnMolalityActCoeff();
doublereal R = GasConstant;
doublereal mm;
/*
* 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]);
@ -955,8 +984,7 @@ namespace Cantera {
}
double xmolSolvent = moleFraction(m_indexSolvent);
mm = fmaxx(SmallNumber, xmolSolvent);
sbar[m_indexSolvent] -=
R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
sbar[m_indexSolvent] -= R *(log(mm) + m_lnActCoeffMolal[m_indexSolvent]);
/*
* Check to see whether activity coefficients are temperature
* dependent. If they are, then calculate the their temperature
@ -967,18 +995,16 @@ namespace Cantera {
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
* 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 equal to the
* species standard state molar volumes. However, extensions
* to this will be implemented in the future.
* For this solution, the partial molar volumes are a
* complex function of pressure.
*
* The general relation is
*
@ -988,12 +1014,10 @@ namespace Cantera {
*
* = V0_i + RT d(lnActCoeffi)dP _T,M
*
* So, if the activity coefficients depended on pressure this
* function would be nontrivial.
*/
void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const {
/*
* Get the standard state values
* Get the standard state values in m^3 kmol-1
*/
getStandardVolumes(vbar);
/*
@ -1027,32 +1051,22 @@ namespace Cantera {
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.
* Update the activity coefficients, This also update the
* internally storred molalities.
*/
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]);
}
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 ------------------
@ -4250,7 +4264,7 @@ namespace Cantera {
#endif
}
/***********************************************************************************************/
/********************************************************************************************/
/**
* s_Pitzer_dlnMolalityActCoeff_dP() (private, const )

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

@ -190,7 +190,7 @@ namespace Cantera {
/// Molar heat capacity at constant pressure. Units: J/kmol/K.
/*
* (HKM -> Bump up to Parent object)
*
*/
virtual doublereal cp_mole() const;