From 67d7e3f290f231494bce4c1d2415764e77737e2f Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 18 Jul 2007 19:27:51 +0000 Subject: [PATCH] Cropping Update -Did a better job at coming up with cropping strategy for molalalities which are used for the calculation of the activity coefficents. Without such a strategy, act coeffs may become unbounded and cause the equilibrium solver to fail. --- Cantera/src/thermo/HMWSoln.cpp | 182 ++++++++++++++++++++++++++++----- Cantera/src/thermo/HMWSoln.h | 13 +++ 2 files changed, 171 insertions(+), 24 deletions(-) diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index ad6fbc8ef..002342100 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -46,6 +46,7 @@ namespace Cantera { m_waterSS(0), m_densWaterSS(1000.), m_waterProps(0), + m_molalitiesAreCropped(false), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -75,6 +76,7 @@ namespace Cantera { m_waterSS(0), m_densWaterSS(1000.), m_waterProps(0), + m_molalitiesAreCropped(false), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -98,6 +100,7 @@ namespace Cantera { m_waterSS(0), m_densWaterSS(1000.), m_waterProps(0), + m_molalitiesAreCropped(false), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -127,6 +130,7 @@ namespace Cantera { m_waterSS(0), m_densWaterSS(1000.), m_waterProps(0), + m_molalitiesAreCropped(false), m_debugCalc(0) { /* @@ -243,6 +247,8 @@ namespace Cantera { m_gamma = b.m_gamma; m_CounterIJ = b.m_CounterIJ; + m_molalitiesCropped = b.m_molalitiesCropped; + m_molalitiesAreCropped= b.m_molalitiesAreCropped; m_debugCalc = b.m_debugCalc; } return *this; @@ -708,7 +714,7 @@ namespace Cantera { // ------- 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$ @@ -729,7 +735,7 @@ namespace Cantera { } } - /** + /* * 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, @@ -751,7 +757,7 @@ namespace Cantera { return 1.0 / mvSolvent; } - /** + /* * Returns the natural logarithm of the standard * concentration of the kth species */ @@ -760,7 +766,7 @@ namespace Cantera { 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 @@ -794,7 +800,7 @@ namespace Cantera { } - /** + /* * Get the array of non-dimensional activities at * the current solution temperature, pressure, and * solution concentration. @@ -821,7 +827,7 @@ namespace Cantera { exp(m_lnActCoeffMolal[m_indexSolvent]) * xmolSolvent; } - /** + /* * getMolalityActivityCoefficients() (virtual, const) * * Get the array of non-dimensional Molality based @@ -1149,9 +1155,6 @@ namespace Cantera { } /* - * - * getPureGibbs() - * * Get the Gibbs functions for the pure species * at the current T and P of the solution. * We assume an incompressible constant partial molar @@ -1168,7 +1171,6 @@ namespace Cantera { } /* - * * getEnthalpy_RT() (virtual, const) * * Get the array of nondimensional Enthalpy functions for the ss @@ -1747,7 +1749,7 @@ namespace Cantera { m_pe.resize(leng, 0.0); m_pp.resize(leng, 0.0); m_tmpV.resize(leng, 0.0); - + m_molalitiesCropped.resize(leng, 0.0); int maxCounterIJlen = 1 + (leng-1) * (leng-2) / 2; @@ -1838,6 +1840,7 @@ namespace Cantera { m_gamma.resize(leng, 0.0); + counterIJ_setup(); } @@ -1854,6 +1857,11 @@ namespace Cantera { * State objects' data. */ calcMolalities(); + /* + * Calculate a cropped set of molalities that will be used + * in all activity coefficent calculations. + */ + calcMolalitiesCropped(); /* * Calculate the stoichiometric ionic charge. This isn't used in the * Pitzer formulation. @@ -1887,9 +1895,134 @@ namespace Cantera { s_updatePitzerSublnMolalityActCoeff(); } + + /* + * Calculate cropped molalities + */ + void HMWSoln::calcMolalitiesCropped() const { + int i, j, k; + doublereal Imax = 0.0, Itmp; + doublereal Iac_max; + m_molalitiesAreCropped = false; + + for (k = 0; k < m_kk; k++) { + m_molalitiesCropped[k] = m_molalities[k]; + double charge = m_speciesCharge[k]; + Itmp = m_molalities[k] * charge * charge; + if (Itmp > Imax) { + Imax = Itmp; + } + } + /* + * Quick return + */ + if (Imax < m_maxIionicStrength) { + return; + } + + m_molalitiesAreCropped = true; + + for (i = 1; i < (m_kk - 1); i++) { + double charge_i = m_speciesCharge[i]; + double abs_charge_i = fabs(charge_i); + if (charge_i == 0.0) { + continue; + } + for (j = (i+1); j < m_kk; j++) { + double charge_j = m_speciesCharge[j]; + double abs_charge_j = fabs(charge_j); + /* + * Find the counterIJ for the symmetric binary interaction + */ + //n = m_kk*i + j; + //counterIJ = m_CounterIJ[n]; + /* + * Only loop over oppositely charge species + */ + if (charge_i * charge_j < 0) { + Iac_max = m_maxIionicStrength; + + if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) { + Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i; + if (Imax > Iac_max) { + m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i); + } + Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i); + if (Imax > Iac_max) { + m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i); + } + } else { + Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j; + if (Imax > Iac_max) { + m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j); + } + Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i; + if (Imax > Iac_max) { + m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i); + } + } + } + } + } + + /* + * Do this loop 10 times until we have achieved charge neutrality + * in the cropped molalities + */ + for (int times = 0; times< 10; times++) { + double anion_charge = 0.0; + double cation_charge = 0.0; + int anion_contrib_max_i = -1; + double anion_contrib_max = -1.0; + int cation_contrib_max_i = -1; + double cation_contrib_max = -1.0; + for (i = 0; i < m_kk; i++) { + double charge_i = m_speciesCharge[i]; + if (charge_i < 0.0) { + double anion_contrib = - m_molalitiesCropped[i] * charge_i; + anion_charge += anion_contrib ; + if (anion_contrib > anion_contrib_max) { + anion_contrib_max = anion_contrib; + anion_contrib_max_i = i; + } + } else if (charge_i > 0.0) { + double cation_contrib = m_molalitiesCropped[i] * charge_i; + cation_charge += cation_contrib ; + if (cation_contrib > cation_contrib_max) { + cation_contrib_max = cation_contrib; + cation_contrib_max_i = i; + } + } + } + double total_charge = cation_charge - anion_charge; + if (total_charge > 1.0E-8) { + double desiredCrop = total_charge/m_speciesCharge[cation_contrib_max_i]; + double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i]; + if (desiredCrop < maxCrop) { + m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop; + break; + } else { + m_molalitiesCropped[cation_contrib_max_i] -= maxCrop; + } + } else if (total_charge < -1.0E-8) { + double desiredCrop = total_charge/m_speciesCharge[anion_contrib_max_i]; + double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i]; + if (desiredCrop < maxCrop) { + m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop; + break; + } else { + m_molalitiesCropped[anion_contrib_max_i] -= maxCrop; + } + } else { + break; + } + } + + } + /* * Set up a counter variable for keeping track of symmetric binary - * interactactions amongst the solute species. + * interactions amongst the solute species. * * n = m_kk*i + j * m_Counter[n] = counter @@ -2053,7 +2186,8 @@ namespace Cantera { } } - /** + + /* * Calculate the Pitzer portion of the activity coefficients. * * This is the main routine in the whole module. It calculates the @@ -2082,9 +2216,9 @@ namespace Cantera { std::string sni, snj, snk; /* - * This is the molality of the species in solution. + * Use the CROPPED molality of the species in solution. */ - const double *molality = DATA_PTR(m_molalities); + const double *molality = DATA_PTR(m_molalitiesCropped); /* * These are the charges of the species accessed from Constituents.h */ @@ -2878,7 +3012,7 @@ namespace Cantera { /*************************************************************************************/ - /** + /* * Calculate the Pitzer portion of the temperature * derivative of the log activity coefficients. * This is an internal routine. @@ -2906,7 +3040,7 @@ namespace Cantera { double d_wateract_dT; std::string sni, snj, snk; - const double *molality = DATA_PTR(m_molalities); + const double *molality = DATA_PTR(m_molalitiesCropped); const double *charge = DATA_PTR(m_speciesCharge); const double *beta0MX_L = DATA_PTR(m_Beta0MX_ij_L); const double *beta1MX_L = DATA_PTR(m_Beta1MX_ij_L); @@ -3650,7 +3784,7 @@ namespace Cantera { /*************************************************************************************/ - /** + /* * s_update_d2lnMolalityActCoeff_dT2() (private, const ) * * Using internally stored values, this function calculates @@ -3685,7 +3819,7 @@ namespace Cantera { std::string sni, snj, snk; - const double *molality = DATA_PTR(m_molalities); + const double *molality = DATA_PTR(m_molalitiesCropped); const double *charge = DATA_PTR(m_speciesCharge); const double *beta0MX_LL= DATA_PTR(m_Beta0MX_ij_LL); const double *beta1MX_LL= DATA_PTR(m_Beta1MX_ij_LL); @@ -4436,7 +4570,7 @@ namespace Cantera { /********************************************************************************************/ - /** + /* * s_Pitzer_dlnMolalityActCoeff_dP() (private, const ) * * Using internally stored values, this function calculates @@ -4456,7 +4590,7 @@ namespace Cantera { s_update_dlnMolalityActCoeff_dP(); } - /** + /* * s_update_dlnMolalityActCoeff_dP() (private, const ) * * Using internally stored values, this function calculates @@ -4492,7 +4626,7 @@ namespace Cantera { double d_wateract_dP; std::string sni, snj, snk; - const double *molality = DATA_PTR(m_molalities); + const double *molality = DATA_PTR(m_molalitiesCropped); const double *charge = DATA_PTR(m_speciesCharge); const double *beta0MX_P = DATA_PTR(m_Beta0MX_ij_P); const double *beta1MX_P = DATA_PTR(m_Beta1MX_ij_P); @@ -5381,7 +5515,7 @@ namespace Cantera { std::string sni, snj; calcMolalities(); const double *charge = DATA_PTR(m_speciesCharge); - double *molality = DATA_PTR(m_molalities); + double *molality = DATA_PTR(m_molalitiesCropped); double *moleF = DATA_PTR(m_tmpV); /* * Update the coefficients wrt Temperature @@ -5390,7 +5524,7 @@ namespace Cantera { s_updatePitzerCoeffWRTemp(2); getMoleFractions(moleF); - printf("Index Name MoleF Molality Charge\n"); + printf("Index Name MoleF MolalityCropped Charge\n"); for (k = 0; k < m_kk; k++) { sni = speciesName(k); printf("%2d %-16s %14.7le %14.7le %5.1f \n", diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index f26a5382c..e852b8a30 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -2852,6 +2852,12 @@ namespace Cantera { * -------- Temporary Variables Used in the Activity Coeff Calc */ + //! Cropped values of the molalities used in activity + //! coefficient calculations + mutable vector_fp m_molalitiesCropped; + + //! Boolean indicating whether the molalities are cropped + mutable bool m_molalitiesAreCropped; //! a counter variable for keeping track of symmetric binary //! interactions amongst the solute species. @@ -3170,6 +3176,13 @@ namespace Cantera { */ void counterIJ_setup() const; + //! Calculate the cropped molalities + /*! + * This is an internal routine that calculates values + * of m_molalitiesCropped from m_molalities + */ + void calcMolalitiesCropped() const; + //! Process an XML node called "binarySaltParameters" /*! * This node contains all of the parameters necessary to describe