Cantera  2.3.0
HMWSoln.cpp
Go to the documentation of this file.
1 /**
2  * @file HMWSoln.cpp
3  * Definitions for the HMWSoln ThermoPhase object, which
4  * models concentrated electrolyte solutions
5  * (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
6  *
7  * Class HMWSoln represents a concentrated liquid electrolyte phase which obeys
8  * the Pitzer formulation for nonideality using molality-based standard states.
9  *
10  * This version of the code was modified to have the binary Beta2 Pitzer
11  * parameter consistent with the temperature expansions used for Beta0,
12  * Beta1, and Cphi.(CFJC, SNL)
13  */
14 
15 // This file is part of Cantera. See License.txt in the top-level directory or
16 // at http://www.cantera.org/license.txt for license and copyright information.
17 
18 #include "cantera/thermo/HMWSoln.h"
19 #include "cantera/thermo/ThermoFactory.h"
20 #include "cantera/thermo/PDSS_Water.h"
21 #include "cantera/thermo/electrolytes.h"
22 #include "cantera/base/stringUtils.h"
23 
24 namespace Cantera
25 {
26 
27 HMWSoln::HMWSoln() :
28  m_formPitzer(PITZERFORM_BASE),
29  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
30  m_formGC(2),
31  m_IionicMolality(0.0),
32  m_maxIionicStrength(100.0),
33  m_TempPitzerRef(298.15),
34  m_IionicMolalityStoich(0.0),
35  m_form_A_Debye(A_DEBYE_WATER),
36  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
37  m_waterSS(0),
38  m_densWaterSS(1000.),
39  m_molalitiesAreCropped(false),
40  IMS_typeCutoff_(0),
41  IMS_X_o_cutoff_(0.2),
42  IMS_gamma_o_min_(1.0E-5),
43  IMS_gamma_k_min_(10.0),
44  IMS_cCut_(0.05),
45  IMS_slopefCut_(0.6),
46  IMS_slopegCut_(0.0),
47  IMS_dfCut_(0.0),
48  IMS_efCut_(0.0),
49  IMS_afCut_(0.0),
50  IMS_bfCut_(0.0),
51  IMS_dgCut_(0.0),
52  IMS_egCut_(0.0),
53  IMS_agCut_(0.0),
54  IMS_bgCut_(0.0),
55  MC_X_o_cutoff_(0.0),
56  MC_X_o_min_(0.0),
57  MC_slopepCut_(0.0),
58  MC_dpCut_(0.0),
59  MC_epCut_(0.0),
60  MC_apCut_(0.0),
61  MC_bpCut_(0.0),
62  MC_cpCut_(0.0),
63  CROP_ln_gamma_o_min(-6.0),
64  CROP_ln_gamma_o_max(3.0),
65  CROP_ln_gamma_k_min(-5.0),
66  CROP_ln_gamma_k_max(15.0),
67  m_last_is(-1.0),
68  m_debugCalc(0)
69 {
70 }
71 
72 HMWSoln::HMWSoln(const std::string& inputFile, const std::string& id_) :
73  m_formPitzer(PITZERFORM_BASE),
74  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
75  m_formGC(2),
76  m_IionicMolality(0.0),
77  m_maxIionicStrength(100.0),
78  m_TempPitzerRef(298.15),
79  m_IionicMolalityStoich(0.0),
80  m_form_A_Debye(A_DEBYE_WATER),
81  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
82  m_waterSS(0),
83  m_densWaterSS(1000.),
84  m_molalitiesAreCropped(false),
85  IMS_typeCutoff_(0),
86  IMS_X_o_cutoff_(0.2),
87  IMS_gamma_o_min_(1.0E-5),
88  IMS_gamma_k_min_(10.0),
89  IMS_cCut_(0.05),
90  IMS_slopefCut_(0.6),
91  IMS_slopegCut_(0.0),
92  IMS_dfCut_(0.0),
93  IMS_efCut_(0.0),
94  IMS_afCut_(0.0),
95  IMS_bfCut_(0.0),
96  IMS_dgCut_(0.0),
97  IMS_egCut_(0.0),
98  IMS_agCut_(0.0),
99  IMS_bgCut_(0.0),
100  MC_X_o_cutoff_(0.0),
101  MC_X_o_min_(0.0),
102  MC_slopepCut_(0.0),
103  MC_dpCut_(0.0),
104  MC_epCut_(0.0),
105  MC_apCut_(0.0),
106  MC_bpCut_(0.0),
107  MC_cpCut_(0.0),
108  CROP_ln_gamma_o_min(-6.0),
109  CROP_ln_gamma_o_max(3.0),
110  CROP_ln_gamma_k_min(-5.0),
111  CROP_ln_gamma_k_max(15.0),
112  m_last_is(-1.0),
113  m_debugCalc(0)
114 {
115  initThermoFile(inputFile, id_);
116 }
117 
118 HMWSoln::HMWSoln(XML_Node& phaseRoot, const std::string& id_) :
119  m_formPitzer(PITZERFORM_BASE),
120  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
121  m_formGC(2),
122  m_IionicMolality(0.0),
123  m_maxIionicStrength(100.0),
124  m_TempPitzerRef(298.15),
125  m_IionicMolalityStoich(0.0),
126  m_form_A_Debye(A_DEBYE_WATER),
127  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
128  m_waterSS(0),
129  m_densWaterSS(1000.),
130  m_molalitiesAreCropped(false),
131  IMS_typeCutoff_(0),
132  IMS_X_o_cutoff_(0.2),
133  IMS_gamma_o_min_(1.0E-5),
134  IMS_gamma_k_min_(10.0),
135  IMS_cCut_(0.05),
136  IMS_slopefCut_(0.6),
137  IMS_slopegCut_(0.0),
138  IMS_dfCut_(0.0),
139  IMS_efCut_(0.0),
140  IMS_afCut_(0.0),
141  IMS_bfCut_(0.0),
142  IMS_dgCut_(0.0),
143  IMS_egCut_(0.0),
144  IMS_agCut_(0.0),
145  IMS_bgCut_(0.0),
146  MC_X_o_cutoff_(0.0),
147  MC_X_o_min_(0.0),
148  MC_slopepCut_(0.0),
149  MC_dpCut_(0.0),
150  MC_epCut_(0.0),
151  MC_apCut_(0.0),
152  MC_bpCut_(0.0),
153  MC_cpCut_(0.0),
154  CROP_ln_gamma_o_min(-6.0),
155  CROP_ln_gamma_o_max(3.0),
156  CROP_ln_gamma_k_min(-5.0),
157  CROP_ln_gamma_k_max(15.0),
158  m_last_is(-1.0),
159  m_debugCalc(0)
160 {
161  importPhase(phaseRoot, this);
162 }
163 
164 HMWSoln::HMWSoln(const HMWSoln& b) :
165  m_formPitzer(PITZERFORM_BASE),
166  m_formPitzerTemp(PITZER_TEMP_CONSTANT),
167  m_formGC(2),
168  m_IionicMolality(0.0),
169  m_maxIionicStrength(100.0),
170  m_TempPitzerRef(298.15),
171  m_IionicMolalityStoich(0.0),
172  m_form_A_Debye(A_DEBYE_WATER),
173  m_A_Debye(1.172576), // units = sqrt(kg/gmol)
174  m_waterSS(0),
175  m_densWaterSS(1000.),
176  m_molalitiesAreCropped(false),
177  IMS_typeCutoff_(0),
178  IMS_X_o_cutoff_(0.2),
179  IMS_gamma_o_min_(1.0E-5),
180  IMS_gamma_k_min_(10.0),
181  IMS_cCut_(0.05),
182  IMS_slopefCut_(0.6),
183  IMS_slopegCut_(0.0),
184  IMS_dfCut_(0.0),
185  IMS_efCut_(0.0),
186  IMS_afCut_(0.0),
187  IMS_bfCut_(0.0),
188  IMS_dgCut_(0.0),
189  IMS_egCut_(0.0),
190  IMS_agCut_(0.0),
191  IMS_bgCut_(0.0),
192  MC_X_o_cutoff_(0.0),
193  MC_X_o_min_(0.0),
194  MC_slopepCut_(0.0),
195  MC_dpCut_(0.0),
196  MC_epCut_(0.0),
197  MC_apCut_(0.0),
198  MC_bpCut_(0.0),
199  MC_cpCut_(0.0),
200  CROP_ln_gamma_o_min(-6.0),
201  CROP_ln_gamma_o_max(3.0),
202  CROP_ln_gamma_k_min(-5.0),
203  CROP_ln_gamma_k_max(15.0),
204  m_last_is(-1.0),
205  m_debugCalc(0)
206 {
207  // Use the assignment operator to do the brunt of the work for the copy
208  // constructor.
209  *this = b;
210 }
211 
212 HMWSoln& HMWSoln::operator=(const HMWSoln& b)
213 {
214  if (&b != this) {
215  MolalityVPSSTP::operator=(b);
216  m_formPitzer = b.m_formPitzer;
217  m_formPitzerTemp = b.m_formPitzerTemp;
218  m_formGC = b.m_formGC;
219  m_Aionic = b.m_Aionic;
220  m_IionicMolality = b.m_IionicMolality;
221  m_maxIionicStrength = b.m_maxIionicStrength;
222  m_TempPitzerRef = b.m_TempPitzerRef;
223  m_IionicMolalityStoich= b.m_IionicMolalityStoich;
224  m_form_A_Debye = b.m_form_A_Debye;
225  m_A_Debye = b.m_A_Debye;
226 
227  // This is an internal shallow copy of the PDSS_Water pointer
228  m_waterSS = providePDSS(0);
229  if (!m_waterSS) {
230  throw CanteraError("HMWSoln::operator=()", "Dynamic cast to PDSS_Water failed");
231  }
232 
233  m_densWaterSS = b.m_densWaterSS;
234 
235  m_waterProps = 0;
236  if (b.m_waterProps) {
237  m_waterProps.reset(new WaterProps(&dynamic_cast<PDSS_Water&>(*m_waterSS)));
238  }
239 
240  m_tmpV = b.m_tmpV;
241  m_speciesCharge_Stoich= b.m_speciesCharge_Stoich;
242  m_Beta0MX_ij = b.m_Beta0MX_ij;
243  m_Beta0MX_ij_L = b.m_Beta0MX_ij_L;
244  m_Beta0MX_ij_LL = b.m_Beta0MX_ij_LL;
245  m_Beta0MX_ij_P = b.m_Beta0MX_ij_P;
246  m_Beta0MX_ij_coeff = b.m_Beta0MX_ij_coeff;
247  m_Beta1MX_ij = b.m_Beta1MX_ij;
248  m_Beta1MX_ij_L = b.m_Beta1MX_ij_L;
249  m_Beta1MX_ij_LL = b.m_Beta1MX_ij_LL;
250  m_Beta1MX_ij_P = b.m_Beta1MX_ij_P;
251  m_Beta1MX_ij_coeff = b.m_Beta1MX_ij_coeff;
252  m_Beta2MX_ij = b.m_Beta2MX_ij;
253  m_Beta2MX_ij_L = b.m_Beta2MX_ij_L;
254  m_Beta2MX_ij_LL = b.m_Beta2MX_ij_LL;
255  m_Beta2MX_ij_P = b.m_Beta2MX_ij_P;
256  m_Beta2MX_ij_coeff = b.m_Beta2MX_ij_coeff;
257  m_Alpha1MX_ij = b.m_Alpha1MX_ij;
258  m_Alpha2MX_ij = b.m_Alpha2MX_ij;
259  m_CphiMX_ij = b.m_CphiMX_ij;
260  m_CphiMX_ij_L = b.m_CphiMX_ij_L;
261  m_CphiMX_ij_LL = b.m_CphiMX_ij_LL;
262  m_CphiMX_ij_P = b.m_CphiMX_ij_P;
263  m_CphiMX_ij_coeff = b.m_CphiMX_ij_coeff;
264  m_Theta_ij = b.m_Theta_ij;
265  m_Theta_ij_L = b.m_Theta_ij_L;
266  m_Theta_ij_LL = b.m_Theta_ij_LL;
267  m_Theta_ij_P = b.m_Theta_ij_P;
268  m_Theta_ij_coeff = b.m_Theta_ij_coeff;
269  m_Psi_ijk = b.m_Psi_ijk;
270  m_Psi_ijk_L = b.m_Psi_ijk_L;
271  m_Psi_ijk_LL = b.m_Psi_ijk_LL;
272  m_Psi_ijk_P = b.m_Psi_ijk_P;
273  m_Psi_ijk_coeff = b.m_Psi_ijk_coeff;
274  m_Lambda_nj = b.m_Lambda_nj;
275  m_Lambda_nj_L = b.m_Lambda_nj_L;
276  m_Lambda_nj_LL = b.m_Lambda_nj_LL;
277  m_Lambda_nj_P = b.m_Lambda_nj_P;
278  m_Lambda_nj_coeff = b.m_Lambda_nj_coeff;
279 
280  m_Mu_nnn = b.m_Mu_nnn;
281  m_Mu_nnn_L = b.m_Mu_nnn_L;
282  m_Mu_nnn_LL = b.m_Mu_nnn_LL;
283  m_Mu_nnn_P = b.m_Mu_nnn_P;
284  m_Mu_nnn_coeff = b.m_Mu_nnn_coeff;
285 
286  m_lnActCoeffMolal_Scaled = b.m_lnActCoeffMolal_Scaled;
287  m_lnActCoeffMolal_Unscaled = b.m_lnActCoeffMolal_Unscaled;
288  m_dlnActCoeffMolaldT_Scaled = b.m_dlnActCoeffMolaldT_Scaled;
289  m_dlnActCoeffMolaldT_Unscaled = b.m_dlnActCoeffMolaldT_Unscaled;
290  m_d2lnActCoeffMolaldT2_Scaled = b.m_d2lnActCoeffMolaldT2_Scaled;
291  m_d2lnActCoeffMolaldT2_Unscaled = b.m_d2lnActCoeffMolaldT2_Unscaled;
292  m_dlnActCoeffMolaldP_Scaled = b.m_dlnActCoeffMolaldP_Scaled;
293  m_dlnActCoeffMolaldP_Unscaled = b.m_dlnActCoeffMolaldP_Unscaled;
294 
295  m_molalitiesCropped = b.m_molalitiesCropped;
296  m_molalitiesAreCropped = b.m_molalitiesAreCropped;
297  m_CounterIJ = b.m_CounterIJ;
298  m_gfunc_IJ = b.m_gfunc_IJ;
299  m_g2func_IJ = b.m_g2func_IJ;
300  m_hfunc_IJ = b.m_hfunc_IJ;
301  m_h2func_IJ = b.m_h2func_IJ;
302  m_BMX_IJ = b.m_BMX_IJ;
303  m_BMX_IJ_L = b.m_BMX_IJ_L;
304  m_BMX_IJ_LL = b.m_BMX_IJ_LL;
305  m_BMX_IJ_P = b.m_BMX_IJ_P;
306  m_BprimeMX_IJ = b.m_BprimeMX_IJ;
307  m_BprimeMX_IJ_L = b.m_BprimeMX_IJ_L;
308  m_BprimeMX_IJ_LL = b.m_BprimeMX_IJ_LL;
309  m_BprimeMX_IJ_P = b.m_BprimeMX_IJ_P;
310  m_BphiMX_IJ = b.m_BphiMX_IJ;
311  m_BphiMX_IJ_L = b.m_BphiMX_IJ_L;
312  m_BphiMX_IJ_LL = b.m_BphiMX_IJ_LL;
313  m_BphiMX_IJ_P = b.m_BphiMX_IJ_P;
314  m_Phi_IJ = b.m_Phi_IJ;
315  m_Phi_IJ_L = b.m_Phi_IJ_L;
316  m_Phi_IJ_LL = b.m_Phi_IJ_LL;
317  m_Phi_IJ_P = b.m_Phi_IJ_P;
318  m_Phiprime_IJ = b.m_Phiprime_IJ;
319  m_PhiPhi_IJ = b.m_PhiPhi_IJ;
320  m_PhiPhi_IJ_L = b.m_PhiPhi_IJ_L;
321  m_PhiPhi_IJ_LL = b.m_PhiPhi_IJ_LL;
322  m_PhiPhi_IJ_P = b.m_PhiPhi_IJ_P;
323  m_CMX_IJ = b.m_CMX_IJ;
324  m_CMX_IJ_L = b.m_CMX_IJ_L;
325  m_CMX_IJ_LL = b.m_CMX_IJ_LL;
326  m_CMX_IJ_P = b.m_CMX_IJ_P;
327  m_gamma_tmp = b.m_gamma_tmp;
328 
329  IMS_lnActCoeffMolal_ = b.IMS_lnActCoeffMolal_;
330  IMS_typeCutoff_ = b.IMS_typeCutoff_;
331  IMS_X_o_cutoff_ = b.IMS_X_o_cutoff_;
332  IMS_gamma_o_min_ = b.IMS_gamma_o_min_;
333  IMS_gamma_k_min_ = b.IMS_gamma_k_min_;
334  IMS_cCut_ = b.IMS_cCut_;
335  IMS_slopefCut_ = b.IMS_slopefCut_;
336  IMS_dfCut_ = b.IMS_dfCut_;
337  IMS_efCut_ = b.IMS_efCut_;
338  IMS_afCut_ = b.IMS_afCut_;
339  IMS_bfCut_ = b.IMS_bfCut_;
340  IMS_slopegCut_ = b.IMS_slopegCut_;
341  IMS_dgCut_ = b.IMS_dgCut_;
342  IMS_egCut_ = b.IMS_egCut_;
343  IMS_agCut_ = b.IMS_agCut_;
344  IMS_bgCut_ = b.IMS_bgCut_;
345  MC_X_o_cutoff_ = b.MC_X_o_cutoff_;
346  MC_X_o_min_ = b.MC_X_o_min_;
347  MC_slopepCut_ = b.MC_slopepCut_;
348  MC_dpCut_ = b.MC_dpCut_;
349  MC_epCut_ = b.MC_epCut_;
350  MC_apCut_ = b.MC_apCut_;
351  MC_bpCut_ = b.MC_bpCut_;
352  MC_cpCut_ = b.MC_cpCut_;
353  CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min;
354  CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max;
355  CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min;
356  CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max;
357  CROP_speciesCropped_ = b.CROP_speciesCropped_;
358  m_debugCalc = b.m_debugCalc;
359  }
360  return *this;
361 }
362 
363 HMWSoln::~HMWSoln()
364 {
365 }
366 
367 ThermoPhase* HMWSoln::duplMyselfAsThermoPhase() const
368 {
369  return new HMWSoln(*this);
370 }
371 
372 int HMWSoln::eosType() const
373 {
374  warn_deprecated("HMWSoln::eosType",
375  "To be removed after Cantera 2.3.");
376  int res;
377  switch (m_formGC) {
378  case 0:
379  res = cHMWSoln0;
380  break;
381  case 1:
382  res = cHMWSoln1;
383  break;
384  case 2:
385  res = cHMWSoln2;
386  break;
387  default:
388  throw CanteraError("eosType", "Unknown type");
389  }
390  return res;
391 }
392 
393 // -------- Molar Thermodynamic Properties of the Solution ---------------
394 
395 doublereal HMWSoln::enthalpy_mole() const
396 {
397  getPartialMolarEnthalpies(m_tmpV.data());
398  return mean_X(m_tmpV);
399 }
400 
401 doublereal HMWSoln::relative_enthalpy() const
402 {
403  getPartialMolarEnthalpies(m_tmpV.data());
404  double hbar = mean_X(m_tmpV);
405  getEnthalpy_RT(m_gamma_tmp.data());
406  for (size_t k = 0; k < m_kk; k++) {
407  m_gamma_tmp[k] *= RT();
408  }
409  double h0bar = mean_X(m_gamma_tmp);
410  return hbar - h0bar;
411 }
412 
413 doublereal HMWSoln::relative_molal_enthalpy() const
414 {
415  double L = relative_enthalpy();
416  getMoleFractions(m_tmpV.data());
417  double xanion = 0.0;
418  size_t kcation = npos;
419  double xcation = 0.0;
420  size_t kanion = npos;
421  for (size_t k = 0; k < m_kk; k++) {
422  if (charge(k) > 0.0) {
423  if (m_tmpV[k] > xanion) {
424  xanion = m_tmpV[k];
425  kanion = k;
426  }
427  } else if (charge(k) < 0.0) {
428  if (m_tmpV[k] > xcation) {
429  xcation = m_tmpV[k];
430  kcation = k;
431  }
432  }
433  }
434  if (kcation == npos || kanion == npos) {
435  return L;
436  }
437  double xuse = xcation;
438  double factor = 1;
439  if (xanion < xcation) {
440  xuse = xanion;
441  if (charge(kcation) != 1.0) {
442  factor = charge(kcation);
443  }
444  } else {
445  if (charge(kanion) != 1.0) {
446  factor = charge(kanion);
447  }
448  }
449  xuse = xuse / factor;
450  return L / xuse;
451 }
452 
453 doublereal HMWSoln::entropy_mole() const
454 {
455  getPartialMolarEntropies(m_tmpV.data());
456  return mean_X(m_tmpV);
457 }
458 
459 doublereal HMWSoln::gibbs_mole() const
460 {
461  getChemPotentials(m_tmpV.data());
462  return mean_X(m_tmpV);
463 }
464 
465 doublereal HMWSoln::cp_mole() const
466 {
467  getPartialMolarCp(m_tmpV.data());
468  return mean_X(m_tmpV);
469 }
470 
471 doublereal HMWSoln::cv_mole() const
472 {
473  double kappa_t = isothermalCompressibility();
474  double beta = thermalExpansionCoeff();
475  double cp = cp_mole();
476  double tt = temperature();
477  double molarV = molarVolume();
478  return cp - beta * beta * tt * molarV / kappa_t;
479 }
480 
481 // ------- Mechanical Equation of State Properties ------------------------
482 
483 void HMWSoln::calcDensity()
484 {
485  static const int cacheId = m_cache.getId();
486  CachedScalar cached = m_cache.getScalar(cacheId);
487  if(cached.validate(temperature(), pressure(), stateMFNumber())) {
488  return;
489  }
490 
491  // Store the internal density of the water SS. Note, we would have to do
492  // this for all other species if they had pressure dependent properties.
493  m_densWaterSS = m_waterSS->density();
494 
495  // Calculate all of the other standard volumes. Note these are constant for
496  // now
497  getPartialMolarVolumes(m_tmpV.data());
498  double dd = meanMolecularWeight() / mean_X(m_tmpV);
499  Phase::setDensity(dd);
500 }
501 
502 void HMWSoln::setDensity(const doublereal rho)
503 {
504  double dens_old = density();
505  if (rho != dens_old) {
506  throw CanteraError("HMWSoln::setDensity",
507  "Density is not an independent variable");
508  }
509 }
510 
511 void HMWSoln::setMolarDensity(const doublereal rho)
512 {
513  throw CanteraError("HMWSoln::setMolarDensity",
514  "Density is not an independent variable");
515 }
516 
517 // ------- Activities and Activity Concentrations
518 
519 void HMWSoln::getActivityConcentrations(doublereal* c) const
520 {
521  double cs_solvent = standardConcentration();
522  getActivities(c);
523  c[0] *= cs_solvent;
524  if (m_kk > 1) {
525  double cs_solute = standardConcentration(1);
526  for (size_t k = 1; k < m_kk; k++) {
527  c[k] *= cs_solute;
528  }
529  }
530 }
531 
532 doublereal HMWSoln::standardConcentration(size_t k) const
533 {
534  getStandardVolumes(m_tmpV.data());
535  double mvSolvent = m_tmpV[m_indexSolvent];
536  if (k > 0) {
537  return m_Mnaught / mvSolvent;
538  }
539  return 1.0 / mvSolvent;
540 }
541 
542 void HMWSoln::getActivities(doublereal* ac) const
543 {
544  updateStandardStateThermo();
545 
546  // Update the molality array, m_molalities(). This requires an update due to
547  // mole fractions
548  s_update_lnMolalityActCoeff();
549 
550  // Now calculate the array of activities.
551  for (size_t k = 0; k < m_kk; k++) {
552  if (k != m_indexSolvent) {
553  ac[k] = m_molalities[k] * exp(m_lnActCoeffMolal_Scaled[k]);
554  }
555  }
556  double xmolSolvent = moleFraction(m_indexSolvent);
557  ac[m_indexSolvent] =
558  exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent;
559 }
560 
561 void HMWSoln::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
562 {
563  updateStandardStateThermo();
564  A_Debye_TP(-1.0, -1.0);
565  s_update_lnMolalityActCoeff();
566  std::copy(m_lnActCoeffMolal_Unscaled.begin(), m_lnActCoeffMolal_Unscaled.end(), acMolality);
567  for (size_t k = 0; k < m_kk; k++) {
568  acMolality[k] = exp(acMolality[k]);
569  }
570 }
571 
572 // ------ Partial Molar Properties of the Solution -----------------
573 
574 void HMWSoln::getChemPotentials(doublereal* mu) const
575 {
576  double xx;
577 
578  // First get the standard chemical potentials in molar form. This requires
579  // updates of standard state as a function of T and P
580  getStandardChemPotentials(mu);
581 
582  // Update the activity coefficients. This also updates the internal molality
583  // array.
584  s_update_lnMolalityActCoeff();
585  double xmolSolvent = moleFraction(m_indexSolvent);
586  for (size_t k = 0; k < m_kk; k++) {
587  if (m_indexSolvent != k) {
588  xx = std::max(m_molalities[k], SmallNumber);
589  mu[k] += RT() * (log(xx) + m_lnActCoeffMolal_Scaled[k]);
590  }
591  }
592  xx = std::max(xmolSolvent, SmallNumber);
593  mu[m_indexSolvent] +=
594  RT() * (log(xx) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
595 }
596 
597 void HMWSoln::getPartialMolarEnthalpies(doublereal* hbar) const
598 {
599  // Get the nondimensional standard state enthalpies
600  getEnthalpy_RT(hbar);
601 
602  // dimensionalize it.
603  for (size_t k = 0; k < m_kk; k++) {
604  hbar[k] *= RT();
605  }
606 
607  // Update the activity coefficients, This also update the internally stored
608  // molalities.
609  s_update_lnMolalityActCoeff();
610  s_update_dlnMolalityActCoeff_dT();
611  for (size_t k = 0; k < m_kk; k++) {
612  hbar[k] -= RT() * temperature() * m_dlnActCoeffMolaldT_Scaled[k];
613  }
614 }
615 
616 void HMWSoln::getPartialMolarEntropies(doublereal* sbar) const
617 {
618  // Get the standard state entropies at the temperature and pressure of the
619  // solution.
620  getEntropy_R(sbar);
621 
622  // Dimensionalize the entropies
623  for (size_t k = 0; k < m_kk; k++) {
624  sbar[k] *= GasConstant;
625  }
626 
627  // Update the activity coefficients, This also update the internally stored
628  // molalities.
629  s_update_lnMolalityActCoeff();
630 
631  // First we will add in the obvious dependence on the T term out front of
632  // the log activity term
633  doublereal mm;
634  for (size_t k = 0; k < m_kk; k++) {
635  if (k != m_indexSolvent) {
636  mm = std::max(SmallNumber, m_molalities[k]);
637  sbar[k] -= GasConstant * (log(mm) + m_lnActCoeffMolal_Scaled[k]);
638  }
639  }
640  double xmolSolvent = moleFraction(m_indexSolvent);
641  mm = std::max(SmallNumber, xmolSolvent);
642  sbar[m_indexSolvent] -= GasConstant *(log(mm) + m_lnActCoeffMolal_Scaled[m_indexSolvent]);
643 
644  // Check to see whether activity coefficients are temperature dependent. If
645  // they are, then calculate the their temperature derivatives and add them
646  // into the result.
647  s_update_dlnMolalityActCoeff_dT();
648  for (size_t k = 0; k < m_kk; k++) {
649  sbar[k] -= RT() * m_dlnActCoeffMolaldT_Scaled[k];
650  }
651 }
652 
653 void HMWSoln::getPartialMolarVolumes(doublereal* vbar) const
654 {
655  // Get the standard state values in m^3 kmol-1
656  getStandardVolumes(vbar);
657 
658  // Update the derivatives wrt the activity coefficients.
659  s_update_lnMolalityActCoeff();
660  s_update_dlnMolalityActCoeff_dP();
661  for (size_t k = 0; k < m_kk; k++) {
662  vbar[k] += RT() * m_dlnActCoeffMolaldP_Scaled[k];
663  }
664 }
665 
666 void HMWSoln::getPartialMolarCp(doublereal* cpbar) const
667 {
668  getCp_R(cpbar);
669  for (size_t k = 0; k < m_kk; k++) {
670  cpbar[k] *= GasConstant;
671  }
672 
673  // Update the activity coefficients, This also update the internally stored
674  // molalities.
675  s_update_lnMolalityActCoeff();
676  s_update_dlnMolalityActCoeff_dT();
677  s_update_d2lnMolalityActCoeff_dT2();
678  for (size_t k = 0; k < m_kk; k++) {
679  cpbar[k] -= (2.0 * RT() * m_dlnActCoeffMolaldT_Scaled[k] +
680  RT() * temperature() * m_d2lnActCoeffMolaldT2_Scaled[k]);
681  }
682 }
683 
684 // -------------- Utilities -------------------------------
685 
686 doublereal HMWSoln::satPressure(doublereal t) {
687  double p_old = pressure();
688  double t_old = temperature();
689  double pres = m_waterSS->satPressure(t);
690 
691  // Set the underlying object back to its original state.
692  m_waterSS->setState_TP(t_old, p_old);
693  return pres;
694 }
695 
696 double HMWSoln::A_Debye_TP(double tempArg, double presArg) const
697 {
698  double T = temperature();
699  double A;
700  if (tempArg != -1.0) {
701  T = tempArg;
702  }
703  double P = pressure();
704  if (presArg != -1.0) {
705  P = presArg;
706  }
707 
708  static const int cacheId = m_cache.getId();
709  CachedScalar cached = m_cache.getScalar(cacheId);
710  if(cached.validate(T, P)) {
711  return m_A_Debye;
712  }
713 
714  switch (m_form_A_Debye) {
715  case A_DEBYE_CONST:
716  A = m_A_Debye;
717  break;
718  case A_DEBYE_WATER:
719  A = m_waterProps->ADebye(T, P, 0);
720  m_A_Debye = A;
721  break;
722  default:
723  throw CanteraError("HMWSoln::A_Debye_TP", "shouldn't be here");
724  }
725  return A;
726 }
727 
728 double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const
729 {
730  doublereal T = temperature();
731  if (tempArg != -1.0) {
732  T = tempArg;
733  }
734  doublereal P = pressure();
735  if (presArg != -1.0) {
736  P = presArg;
737  }
738  doublereal dAdT;
739  switch (m_form_A_Debye) {
740  case A_DEBYE_CONST:
741  dAdT = 0.0;
742  break;
743  case A_DEBYE_WATER:
744  dAdT = m_waterProps->ADebye(T, P, 1);
745  break;
746  default:
747  throw CanteraError("HMWSoln::dA_DebyedT_TP", "shouldn't be here");
748  }
749  return dAdT;
750 }
751 
752 double HMWSoln::dA_DebyedP_TP(double tempArg, double presArg) const
753 {
754  double T = temperature();
755  if (tempArg != -1.0) {
756  T = tempArg;
757  }
758  double P = pressure();
759  if (presArg != -1.0) {
760  P = presArg;
761  }
762 
763  double dAdP;
764  static const int cacheId = m_cache.getId();
765  CachedScalar cached = m_cache.getScalar(cacheId);
766  switch (m_form_A_Debye) {
767  case A_DEBYE_CONST:
768  dAdP = 0.0;
769  break;
770  case A_DEBYE_WATER:
771  if(cached.validate(T, P)) {
772  dAdP = cached.value;
773  } else {
774  dAdP = m_waterProps->ADebye(T, P, 3);
775  cached.value = dAdP;
776  }
777  break;
778  default:
779  throw CanteraError("HMWSoln::dA_DebyedP_TP", "shouldn't be here");
780  }
781  return dAdP;
782 }
783 
784 double HMWSoln::ADebye_L(double tempArg, double presArg) const
785 {
786  double dAdT = dA_DebyedT_TP();
787  double dAphidT = dAdT /3.0;
788  double T = temperature();
789  if (tempArg != -1.0) {
790  T = tempArg;
791  }
792  return dAphidT * (4.0 * GasConstant * T * T);
793 }
794 
795 double HMWSoln::ADebye_V(double tempArg, double presArg) const
796 {
797  double dAdP = dA_DebyedP_TP();
798  double dAphidP = dAdP /3.0;
799  double T = temperature();
800  if (tempArg != -1.0) {
801  T = tempArg;
802  }
803  return - dAphidP * (4.0 * GasConstant * T);
804 }
805 
806 double HMWSoln::ADebye_J(double tempArg, double presArg) const
807 {
808  double T = temperature();
809  if (tempArg != -1.0) {
810  T = tempArg;
811  }
812  double A_L = ADebye_L(T, presArg);
813  double d2 = d2A_DebyedT2_TP(T, presArg);
814  double d2Aphi = d2 / 3.0;
815  return 2.0 * A_L / T + 4.0 * GasConstant * T * T *d2Aphi;
816 }
817 
818 double HMWSoln::d2A_DebyedT2_TP(double tempArg, double presArg) const
819 {
820  double T = temperature();
821  if (tempArg != -1.0) {
822  T = tempArg;
823  }
824  double P = pressure();
825  if (presArg != -1.0) {
826  P = presArg;
827  }
828  double d2AdT2;
829  switch (m_form_A_Debye) {
830  case A_DEBYE_CONST:
831  d2AdT2 = 0.0;
832  break;
833  case A_DEBYE_WATER:
834  d2AdT2 = m_waterProps->ADebye(T, P, 2);
835  break;
836  default:
837  throw CanteraError("HMWSoln::d2A_DebyedT2_TP", "shouldn't be here");
838  }
839  return d2AdT2;
840 }
841 
842 // ---------- Other Property Functions
843 
844 double HMWSoln::AionicRadius(int k) const
845 {
846  return m_Aionic[k];
847 }
848 
849 // ------------ Private and Restricted Functions ------------------
850 
851 void HMWSoln::initLengths()
852 {
853  // Resize lengths equal to the number of species in the phase.
854  m_electrolyteSpeciesType.resize(m_kk, cEST_polarNeutral);
855  m_speciesSize.resize(m_kk);
856  m_speciesCharge_Stoich.resize(m_kk, 0.0);
857  m_Aionic.resize(m_kk, 0.0);
858  m_tmpV.resize(m_kk, 0.0);
859  m_molalitiesCropped.resize(m_kk, 0.0);
860 
861  size_t maxCounterIJlen = 1 + (m_kk-1) * (m_kk-2) / 2;
862 
863  // Figure out the size of the temperature coefficient arrays
864  int TCoeffLength = 1;
865  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
866  TCoeffLength = 2;
867  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
868  TCoeffLength = 5;
869  }
870 
871  m_Beta0MX_ij.resize(maxCounterIJlen, 0.0);
872  m_Beta0MX_ij_L.resize(maxCounterIJlen, 0.0);
873  m_Beta0MX_ij_LL.resize(maxCounterIJlen, 0.0);
874  m_Beta0MX_ij_P.resize(maxCounterIJlen, 0.0);
875  m_Beta0MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
876 
877  m_Beta1MX_ij.resize(maxCounterIJlen, 0.0);
878  m_Beta1MX_ij_L.resize(maxCounterIJlen, 0.0);
879  m_Beta1MX_ij_LL.resize(maxCounterIJlen, 0.0);
880  m_Beta1MX_ij_P.resize(maxCounterIJlen, 0.0);
881  m_Beta1MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
882 
883  m_Beta2MX_ij.resize(maxCounterIJlen, 0.0);
884  m_Beta2MX_ij_L.resize(maxCounterIJlen, 0.0);
885  m_Beta2MX_ij_LL.resize(maxCounterIJlen, 0.0);
886  m_Beta2MX_ij_P.resize(maxCounterIJlen, 0.0);
887  m_Beta2MX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
888 
889  m_CphiMX_ij.resize(maxCounterIJlen, 0.0);
890  m_CphiMX_ij_L.resize(maxCounterIJlen, 0.0);
891  m_CphiMX_ij_LL.resize(maxCounterIJlen, 0.0);
892  m_CphiMX_ij_P.resize(maxCounterIJlen, 0.0);
893  m_CphiMX_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
894 
895  m_Alpha1MX_ij.resize(maxCounterIJlen, 2.0);
896  m_Alpha2MX_ij.resize(maxCounterIJlen, 12.0);
897  m_Theta_ij.resize(maxCounterIJlen, 0.0);
898  m_Theta_ij_L.resize(maxCounterIJlen, 0.0);
899  m_Theta_ij_LL.resize(maxCounterIJlen, 0.0);
900  m_Theta_ij_P.resize(maxCounterIJlen, 0.0);
901  m_Theta_ij_coeff.resize(TCoeffLength, maxCounterIJlen, 0.0);
902 
903  size_t n = m_kk*m_kk*m_kk;
904  m_Psi_ijk.resize(n, 0.0);
905  m_Psi_ijk_L.resize(n, 0.0);
906  m_Psi_ijk_LL.resize(n, 0.0);
907  m_Psi_ijk_P.resize(n, 0.0);
908  m_Psi_ijk_coeff.resize(TCoeffLength, n, 0.0);
909 
910  m_Lambda_nj.resize(m_kk, m_kk, 0.0);
911  m_Lambda_nj_L.resize(m_kk, m_kk, 0.0);
912  m_Lambda_nj_LL.resize(m_kk, m_kk, 0.0);
913  m_Lambda_nj_P.resize(m_kk, m_kk, 0.0);
914  m_Lambda_nj_coeff.resize(TCoeffLength, m_kk * m_kk, 0.0);
915 
916  m_Mu_nnn.resize(m_kk, 0.0);
917  m_Mu_nnn_L.resize(m_kk, 0.0);
918  m_Mu_nnn_LL.resize(m_kk, 0.0);
919  m_Mu_nnn_P.resize(m_kk, 0.0);
920  m_Mu_nnn_coeff.resize(TCoeffLength, m_kk, 0.0);
921 
922  m_lnActCoeffMolal_Scaled.resize(m_kk, 0.0);
923  m_dlnActCoeffMolaldT_Scaled.resize(m_kk, 0.0);
924  m_d2lnActCoeffMolaldT2_Scaled.resize(m_kk, 0.0);
925  m_dlnActCoeffMolaldP_Scaled.resize(m_kk, 0.0);
926  m_lnActCoeffMolal_Unscaled.resize(m_kk, 0.0);
927  m_dlnActCoeffMolaldT_Unscaled.resize(m_kk, 0.0);
928  m_d2lnActCoeffMolaldT2_Unscaled.resize(m_kk, 0.0);
929  m_dlnActCoeffMolaldP_Unscaled.resize(m_kk, 0.0);
930 
931  m_CounterIJ.resize(m_kk*m_kk, 0);
932  m_gfunc_IJ.resize(maxCounterIJlen, 0.0);
933  m_g2func_IJ.resize(maxCounterIJlen, 0.0);
934  m_hfunc_IJ.resize(maxCounterIJlen, 0.0);
935  m_h2func_IJ.resize(maxCounterIJlen, 0.0);
936  m_BMX_IJ.resize(maxCounterIJlen, 0.0);
937  m_BMX_IJ_L.resize(maxCounterIJlen, 0.0);
938  m_BMX_IJ_LL.resize(maxCounterIJlen, 0.0);
939  m_BMX_IJ_P.resize(maxCounterIJlen, 0.0);
940  m_BprimeMX_IJ.resize(maxCounterIJlen, 0.0);
941  m_BprimeMX_IJ_L.resize(maxCounterIJlen, 0.0);
942  m_BprimeMX_IJ_LL.resize(maxCounterIJlen, 0.0);
943  m_BprimeMX_IJ_P.resize(maxCounterIJlen, 0.0);
944  m_BphiMX_IJ.resize(maxCounterIJlen, 0.0);
945  m_BphiMX_IJ_L.resize(maxCounterIJlen, 0.0);
946  m_BphiMX_IJ_LL.resize(maxCounterIJlen, 0.0);
947  m_BphiMX_IJ_P.resize(maxCounterIJlen, 0.0);
948  m_Phi_IJ.resize(maxCounterIJlen, 0.0);
949  m_Phi_IJ_L.resize(maxCounterIJlen, 0.0);
950  m_Phi_IJ_LL.resize(maxCounterIJlen, 0.0);
951  m_Phi_IJ_P.resize(maxCounterIJlen, 0.0);
952  m_Phiprime_IJ.resize(maxCounterIJlen, 0.0);
953  m_PhiPhi_IJ.resize(maxCounterIJlen, 0.0);
954  m_PhiPhi_IJ_L.resize(maxCounterIJlen, 0.0);
955  m_PhiPhi_IJ_LL.resize(maxCounterIJlen, 0.0);
956  m_PhiPhi_IJ_P.resize(maxCounterIJlen, 0.0);
957  m_CMX_IJ.resize(maxCounterIJlen, 0.0);
958  m_CMX_IJ_L.resize(maxCounterIJlen, 0.0);
959  m_CMX_IJ_LL.resize(maxCounterIJlen, 0.0);
960  m_CMX_IJ_P.resize(maxCounterIJlen, 0.0);
961 
962  m_gamma_tmp.resize(m_kk, 0.0);
963  IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
964  CROP_speciesCropped_.resize(m_kk, 0);
965 
966  counterIJ_setup();
967 }
968 
969 void HMWSoln::s_update_lnMolalityActCoeff() const
970 {
971  static const int cacheId = m_cache.getId();
972  CachedScalar cached = m_cache.getScalar(cacheId);
973  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
974  return;
975  }
976 
977  // Calculate the molalities. Currently, the molalities may not be current
978  // with respect to the contents of the State objects' data.
979  calcMolalities();
980 
981  // Calculate a cropped set of molalities that will be used in all activity
982  // coefficient calculations.
983  calcMolalitiesCropped();
984 
985  // Calculate the stoichiometric ionic charge. This isn't used in the Pitzer
986  // formulation.
987  m_IionicMolalityStoich = 0.0;
988  for (size_t k = 0; k < m_kk; k++) {
989  double z_k = charge(k);
990  double zs_k1 = m_speciesCharge_Stoich[k];
991  if (z_k == zs_k1) {
992  m_IionicMolalityStoich += m_molalities[k] * z_k * z_k;
993  } else {
994  double zs_k2 = z_k - zs_k1;
995  m_IionicMolalityStoich += m_molalities[k] * (zs_k1 * zs_k1 + zs_k2 * zs_k2);
996  }
997  }
998 
999  // Update the temperature dependence of the pitzer coefficients and their
1000  // derivatives
1001  s_updatePitzer_CoeffWRTemp();
1002 
1003  // Calculate the IMS cutoff factors
1004  s_updateIMS_lnMolalityActCoeff();
1005 
1006  // Now do the main calculation.
1007  s_updatePitzer_lnMolalityActCoeff();
1008 
1009  double xmolSolvent = moleFraction(m_indexSolvent);
1010  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
1011  double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
1012  double lnxs = log(xx);
1013 
1014  for (size_t k = 1; k < m_kk; k++) {
1015  CROP_speciesCropped_[k] = 0;
1016  m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
1017  if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
1018  CROP_speciesCropped_[k] = 2;
1019  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
1020  }
1021  if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
1022  // -1.0 and -1.5 caused multiple solutions
1023  CROP_speciesCropped_[k] = 2;
1024  m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
1025  }
1026  }
1027  CROP_speciesCropped_[0] = 0;
1028  m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
1029  if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
1030  CROP_speciesCropped_[0] = 2;
1031  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
1032  }
1033  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
1034  CROP_speciesCropped_[0] = 2;
1035  // -0.5 caused multiple solutions
1036  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
1037  }
1038  if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
1039  CROP_speciesCropped_[0] = 2;
1040  m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
1041  }
1042 
1043  // Now do the pH Scaling
1044  s_updateScaling_pHScaling();
1045 }
1046 
1047 void HMWSoln::calcMolalitiesCropped() const
1048 {
1049  doublereal Imax = 0.0;
1050  m_molalitiesAreCropped = false;
1051 
1052  for (size_t k = 0; k < m_kk; k++) {
1053  m_molalitiesCropped[k] = m_molalities[k];
1054  Imax = std::max(m_molalities[k] * charge(k) * charge(k), Imax);
1055  }
1056 
1057  int cropMethod = 1;
1058  if (cropMethod == 0) {
1059  // Quick return
1060  if (Imax < m_maxIionicStrength) {
1061  return;
1062  }
1063 
1064  m_molalitiesAreCropped = true;
1065  for (size_t i = 1; i < (m_kk - 1); i++) {
1066  double charge_i = charge(i);
1067  double abs_charge_i = fabs(charge_i);
1068  if (charge_i == 0.0) {
1069  continue;
1070  }
1071  for (size_t j = (i+1); j < m_kk; j++) {
1072  double charge_j = charge(j);
1073  double abs_charge_j = fabs(charge_j);
1074 
1075  // Only loop over oppositely charge species
1076  if (charge_i * charge_j < 0) {
1077  double Iac_max = m_maxIionicStrength;
1078 
1079  if (m_molalitiesCropped[i] > m_molalitiesCropped[j]) {
1080  Imax = m_molalitiesCropped[i] * abs_charge_i * abs_charge_i;
1081  if (Imax > Iac_max) {
1082  m_molalitiesCropped[i] = Iac_max / (abs_charge_i * abs_charge_i);
1083  }
1084  Imax = m_molalitiesCropped[j] * fabs(abs_charge_j * abs_charge_i);
1085  if (Imax > Iac_max) {
1086  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_i);
1087  }
1088  } else {
1089  Imax = m_molalitiesCropped[j] * abs_charge_j * abs_charge_j;
1090  if (Imax > Iac_max) {
1091  m_molalitiesCropped[j] = Iac_max / (abs_charge_j * abs_charge_j);
1092  }
1093  Imax = m_molalitiesCropped[i] * abs_charge_j * abs_charge_i;
1094  if (Imax > Iac_max) {
1095  m_molalitiesCropped[i] = Iac_max / (abs_charge_j * abs_charge_i);
1096  }
1097  }
1098  }
1099  }
1100  }
1101 
1102  // Do this loop 10 times until we have achieved charge neutrality in the
1103  // cropped molalities
1104  for (int times = 0; times< 10; times++) {
1105  double anion_charge = 0.0;
1106  double cation_charge = 0.0;
1107  size_t anion_contrib_max_i = npos;
1108  double anion_contrib_max = -1.0;
1109  size_t cation_contrib_max_i = npos;
1110  double cation_contrib_max = -1.0;
1111  for (size_t i = 0; i < m_kk; i++) {
1112  double charge_i = charge(i);
1113  if (charge_i < 0.0) {
1114  double anion_contrib = - m_molalitiesCropped[i] * charge_i;
1115  anion_charge += anion_contrib;
1116  if (anion_contrib > anion_contrib_max) {
1117  anion_contrib_max = anion_contrib;
1118  anion_contrib_max_i = i;
1119  }
1120  } else if (charge_i > 0.0) {
1121  double cation_contrib = m_molalitiesCropped[i] * charge_i;
1122  cation_charge += cation_contrib;
1123  if (cation_contrib > cation_contrib_max) {
1124  cation_contrib_max = cation_contrib;
1125  cation_contrib_max_i = i;
1126  }
1127  }
1128  }
1129  double total_charge = cation_charge - anion_charge;
1130  if (total_charge > 1.0E-8) {
1131  double desiredCrop = total_charge/charge(cation_contrib_max_i);
1132  double maxCrop = 0.66 * m_molalitiesCropped[cation_contrib_max_i];
1133  if (desiredCrop < maxCrop) {
1134  m_molalitiesCropped[cation_contrib_max_i] -= desiredCrop;
1135  break;
1136  } else {
1137  m_molalitiesCropped[cation_contrib_max_i] -= maxCrop;
1138  }
1139  } else if (total_charge < -1.0E-8) {
1140  double desiredCrop = total_charge/charge(anion_contrib_max_i);
1141  double maxCrop = 0.66 * m_molalitiesCropped[anion_contrib_max_i];
1142  if (desiredCrop < maxCrop) {
1143  m_molalitiesCropped[anion_contrib_max_i] -= desiredCrop;
1144  break;
1145  } else {
1146  m_molalitiesCropped[anion_contrib_max_i] -= maxCrop;
1147  }
1148  } else {
1149  break;
1150  }
1151  }
1152  }
1153 
1154  if (cropMethod == 1) {
1155  double* molF = m_gamma_tmp.data();
1156  getMoleFractions(molF);
1157  double xmolSolvent = molF[m_indexSolvent];
1158  if (xmolSolvent >= MC_X_o_cutoff_) {
1159  return;
1160  }
1161 
1162  m_molalitiesAreCropped = true;
1163  double poly = MC_apCut_ + MC_bpCut_ * xmolSolvent + MC_dpCut_* xmolSolvent * xmolSolvent;
1164  double p = xmolSolvent + MC_epCut_ + exp(- xmolSolvent/ MC_cpCut_) * poly;
1165  double denomInv = 1.0/ (m_Mnaught * p);
1166  for (size_t k = 0; k < m_kk; k++) {
1167  m_molalitiesCropped[k] = molF[k] * denomInv;
1168  }
1169 
1170  // Do a further check to see if the Ionic strength is below a max value
1171  // Reduce the molalities to enforce this. Note, this algorithm preserves
1172  // the charge neutrality of the solution after cropping.
1173  double Itmp = 0.0;
1174  for (size_t k = 0; k < m_kk; k++) {
1175  Itmp += m_molalitiesCropped[k] * charge(k) * charge(k);
1176  }
1177  if (Itmp > m_maxIionicStrength) {
1178  double ratio = Itmp / m_maxIionicStrength;
1179  for (size_t k = 0; k < m_kk; k++) {
1180  if (charge(k) != 0.0) {
1181  m_molalitiesCropped[k] *= ratio;
1182  }
1183  }
1184  }
1185  }
1186 }
1187 
1188 void HMWSoln::counterIJ_setup() const
1189 {
1190  m_CounterIJ.resize(m_kk * m_kk);
1191  int counter = 0;
1192  for (size_t i = 0; i < m_kk; i++) {
1193  size_t n = i;
1194  size_t nc = m_kk * i;
1195  m_CounterIJ[n] = 0;
1196  m_CounterIJ[nc] = 0;
1197  }
1198  for (size_t i = 1; i < (m_kk - 1); i++) {
1199  size_t n = m_kk * i + i;
1200  m_CounterIJ[n] = 0;
1201  for (size_t j = (i+1); j < m_kk; j++) {
1202  n = m_kk * j + i;
1203  size_t nc = m_kk * i + j;
1204  counter++;
1205  m_CounterIJ[n] = counter;
1206  m_CounterIJ[nc] = counter;
1207  }
1208  }
1209 }
1210 
1211 void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const
1212 {
1213  double T = temperature();
1214  const double twoT = 2.0 * T;
1215  const double invT = 1.0 / T;
1216  const double invT2 = invT * invT;
1217  const double twoinvT3 = 2.0 * invT * invT2;
1218  double tinv = 0.0, tln = 0.0, tlin = 0.0, tquad = 0.0;
1219  if (m_formPitzerTemp == PITZER_TEMP_LINEAR) {
1220  tlin = T - m_TempPitzerRef;
1221  } else if (m_formPitzerTemp == PITZER_TEMP_COMPLEX1) {
1222  tlin = T - m_TempPitzerRef;
1223  tquad = T * T - m_TempPitzerRef * m_TempPitzerRef;
1224  tln = log(T/ m_TempPitzerRef);
1225  tinv = 1.0/T - 1.0/m_TempPitzerRef;
1226  }
1227 
1228  for (size_t i = 1; i < (m_kk - 1); i++) {
1229  for (size_t j = (i+1); j < m_kk; j++) {
1230  // Find the counterIJ for the symmetric binary interaction
1231  size_t n = m_kk*i + j;
1232  size_t counterIJ = m_CounterIJ[n];
1233 
1234  const double* beta0MX_coeff = m_Beta0MX_ij_coeff.ptrColumn(counterIJ);
1235  const double* beta1MX_coeff = m_Beta1MX_ij_coeff.ptrColumn(counterIJ);
1236  const double* beta2MX_coeff = m_Beta2MX_ij_coeff.ptrColumn(counterIJ);
1237  const double* CphiMX_coeff = m_CphiMX_ij_coeff.ptrColumn(counterIJ);
1238  const double* Theta_coeff = m_Theta_ij_coeff.ptrColumn(counterIJ);
1239 
1240  switch (m_formPitzerTemp) {
1241  case PITZER_TEMP_CONSTANT:
1242  break;
1243  case PITZER_TEMP_LINEAR:
1244 
1245  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1246  + beta0MX_coeff[1]*tlin;
1247  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1];
1248  m_Beta0MX_ij_LL[counterIJ] = 0.0;
1249  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1250  + beta1MX_coeff[1]*tlin;
1251  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1];
1252  m_Beta1MX_ij_LL[counterIJ] = 0.0;
1253  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1254  + beta2MX_coeff[1]*tlin;
1255  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1];
1256  m_Beta2MX_ij_LL[counterIJ] = 0.0;
1257  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1258  + CphiMX_coeff[1]*tlin;
1259  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1];
1260  m_CphiMX_ij_LL[counterIJ] = 0.0;
1261  m_Theta_ij[counterIJ] = Theta_coeff[0] + Theta_coeff[1]*tlin;
1262  m_Theta_ij_L[counterIJ] = Theta_coeff[1];
1263  m_Theta_ij_LL[counterIJ] = 0.0;
1264  break;
1265 
1266  case PITZER_TEMP_COMPLEX1:
1267  m_Beta0MX_ij[counterIJ] = beta0MX_coeff[0]
1268  + beta0MX_coeff[1]*tlin
1269  + beta0MX_coeff[2]*tquad
1270  + beta0MX_coeff[3]*tinv
1271  + beta0MX_coeff[4]*tln;
1272  m_Beta1MX_ij[counterIJ] = beta1MX_coeff[0]
1273  + beta1MX_coeff[1]*tlin
1274  + beta1MX_coeff[2]*tquad
1275  + beta1MX_coeff[3]*tinv
1276  + beta1MX_coeff[4]*tln;
1277  m_Beta2MX_ij[counterIJ] = beta2MX_coeff[0]
1278  + beta2MX_coeff[1]*tlin
1279  + beta2MX_coeff[2]*tquad
1280  + beta2MX_coeff[3]*tinv
1281  + beta2MX_coeff[4]*tln;
1282  m_CphiMX_ij[counterIJ] = CphiMX_coeff[0]
1283  + CphiMX_coeff[1]*tlin
1284  + CphiMX_coeff[2]*tquad
1285  + CphiMX_coeff[3]*tinv
1286  + CphiMX_coeff[4]*tln;
1287  m_Theta_ij[counterIJ] = Theta_coeff[0]
1288  + Theta_coeff[1]*tlin
1289  + Theta_coeff[2]*tquad
1290  + Theta_coeff[3]*tinv
1291  + Theta_coeff[4]*tln;
1292  m_Beta0MX_ij_L[counterIJ] = beta0MX_coeff[1]
1293  + beta0MX_coeff[2]*twoT
1294  - beta0MX_coeff[3]*invT2
1295  + beta0MX_coeff[4]*invT;
1296  m_Beta1MX_ij_L[counterIJ] = beta1MX_coeff[1]
1297  + beta1MX_coeff[2]*twoT
1298  - beta1MX_coeff[3]*invT2
1299  + beta1MX_coeff[4]*invT;
1300  m_Beta2MX_ij_L[counterIJ] = beta2MX_coeff[1]
1301  + beta2MX_coeff[2]*twoT
1302  - beta2MX_coeff[3]*invT2
1303  + beta2MX_coeff[4]*invT;
1304  m_CphiMX_ij_L[counterIJ] = CphiMX_coeff[1]
1305  + CphiMX_coeff[2]*twoT
1306  - CphiMX_coeff[3]*invT2
1307  + CphiMX_coeff[4]*invT;
1308  m_Theta_ij_L[counterIJ] = Theta_coeff[1]
1309  + Theta_coeff[2]*twoT
1310  - Theta_coeff[3]*invT2
1311  + Theta_coeff[4]*invT;
1312  doDerivs = 2;
1313  if (doDerivs > 1) {
1314  m_Beta0MX_ij_LL[counterIJ] =
1315  + beta0MX_coeff[2]*2.0
1316  + beta0MX_coeff[3]*twoinvT3
1317  - beta0MX_coeff[4]*invT2;
1318  m_Beta1MX_ij_LL[counterIJ] =
1319  + beta1MX_coeff[2]*2.0
1320  + beta1MX_coeff[3]*twoinvT3
1321  - beta1MX_coeff[4]*invT2;
1322  m_Beta2MX_ij_LL[counterIJ] =
1323  + beta2MX_coeff[2]*2.0
1324  + beta2MX_coeff[3]*twoinvT3
1325  - beta2MX_coeff[4]*invT2;
1326  m_CphiMX_ij_LL[counterIJ] =
1327  + CphiMX_coeff[2]*2.0
1328  + CphiMX_coeff[3]*twoinvT3
1329  - CphiMX_coeff[4]*invT2;
1330  m_Theta_ij_LL[counterIJ] =
1331  + Theta_coeff[2]*2.0
1332  + Theta_coeff[3]*twoinvT3
1333  - Theta_coeff[4]*invT2;
1334  }
1335  break;
1336  }
1337  }
1338  }
1339 
1340  // Lambda interactions and Mu_nnn
1341  // i must be neutral for this term to be nonzero. We take advantage of this
1342  // here to lower the operation count.
1343  for (size_t i = 1; i < m_kk; i++) {
1344  if (charge(i) == 0.0) {
1345  for (size_t j = 1; j < m_kk; j++) {
1346  size_t n = i * m_kk + j;
1347  const double* Lambda_coeff = m_Lambda_nj_coeff.ptrColumn(n);
1348  switch (m_formPitzerTemp) {
1349  case PITZER_TEMP_CONSTANT:
1350  m_Lambda_nj(i,j) = Lambda_coeff[0];
1351  break;
1352  case PITZER_TEMP_LINEAR:
1353  m_Lambda_nj(i,j) = Lambda_coeff[0] + Lambda_coeff[1]*tlin;
1354  m_Lambda_nj_L(i,j) = Lambda_coeff[1];
1355  m_Lambda_nj_LL(i,j) = 0.0;
1356  break;
1357  case PITZER_TEMP_COMPLEX1:
1358  m_Lambda_nj(i,j) = Lambda_coeff[0]
1359  + Lambda_coeff[1]*tlin
1360  + Lambda_coeff[2]*tquad
1361  + Lambda_coeff[3]*tinv
1362  + Lambda_coeff[4]*tln;
1363 
1364  m_Lambda_nj_L(i,j) = Lambda_coeff[1]
1365  + Lambda_coeff[2]*twoT
1366  - Lambda_coeff[3]*invT2
1367  + Lambda_coeff[4]*invT;
1368 
1369  m_Lambda_nj_LL(i,j) =
1370  Lambda_coeff[2]*2.0
1371  + Lambda_coeff[3]*twoinvT3
1372  - Lambda_coeff[4]*invT2;
1373  }
1374 
1375  if (j == i) {
1376  const double* Mu_coeff = m_Mu_nnn_coeff.ptrColumn(i);
1377  switch (m_formPitzerTemp) {
1378  case PITZER_TEMP_CONSTANT:
1379  m_Mu_nnn[i] = Mu_coeff[0];
1380  break;
1381  case PITZER_TEMP_LINEAR:
1382  m_Mu_nnn[i] = Mu_coeff[0] + Mu_coeff[1]*tlin;
1383  m_Mu_nnn_L[i] = Mu_coeff[1];
1384  m_Mu_nnn_LL[i] = 0.0;
1385  break;
1386  case PITZER_TEMP_COMPLEX1:
1387  m_Mu_nnn[i] = Mu_coeff[0]
1388  + Mu_coeff[1]*tlin
1389  + Mu_coeff[2]*tquad
1390  + Mu_coeff[3]*tinv
1391  + Mu_coeff[4]*tln;
1392  m_Mu_nnn_L[i] = Mu_coeff[1]
1393  + Mu_coeff[2]*twoT
1394  - Mu_coeff[3]*invT2
1395  + Mu_coeff[4]*invT;
1396  m_Mu_nnn_LL[i] =
1397  Mu_coeff[2]*2.0
1398  + Mu_coeff[3]*twoinvT3
1399  - Mu_coeff[4]*invT2;
1400  }
1401  }
1402  }
1403  }
1404  }
1405 
1406  switch(m_formPitzerTemp) {
1407  case PITZER_TEMP_CONSTANT:
1408  for (size_t i = 1; i < m_kk; i++) {
1409  for (size_t j = 1; j < m_kk; j++) {
1410  for (size_t k = 1; k < m_kk; k++) {
1411  size_t n = i * m_kk *m_kk + j * m_kk + k;
1412  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1413  m_Psi_ijk[n] = Psi_coeff[0];
1414  }
1415  }
1416  }
1417  break;
1418  case PITZER_TEMP_LINEAR:
1419  for (size_t i = 1; i < m_kk; i++) {
1420  for (size_t j = 1; j < m_kk; j++) {
1421  for (size_t k = 1; k < m_kk; k++) {
1422  size_t n = i * m_kk *m_kk + j * m_kk + k;
1423  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1424  m_Psi_ijk[n] = Psi_coeff[0] + Psi_coeff[1]*tlin;
1425  m_Psi_ijk_L[n] = Psi_coeff[1];
1426  m_Psi_ijk_LL[n] = 0.0;
1427  }
1428  }
1429  }
1430  break;
1431  case PITZER_TEMP_COMPLEX1:
1432  for (size_t i = 1; i < m_kk; i++) {
1433  for (size_t j = 1; j < m_kk; j++) {
1434  for (size_t k = 1; k < m_kk; k++) {
1435  size_t n = i * m_kk *m_kk + j * m_kk + k;
1436  const double* Psi_coeff = m_Psi_ijk_coeff.ptrColumn(n);
1437  m_Psi_ijk[n] = Psi_coeff[0]
1438  + Psi_coeff[1]*tlin
1439  + Psi_coeff[2]*tquad
1440  + Psi_coeff[3]*tinv
1441  + Psi_coeff[4]*tln;
1442  m_Psi_ijk_L[n] = Psi_coeff[1]
1443  + Psi_coeff[2]*twoT
1444  - Psi_coeff[3]*invT2
1445  + Psi_coeff[4]*invT;
1446  m_Psi_ijk_LL[n] =
1447  Psi_coeff[2]*2.0
1448  + Psi_coeff[3]*twoinvT3
1449  - Psi_coeff[4]*invT2;
1450  }
1451  }
1452  }
1453  break;
1454  }
1455 }
1456 
1457 void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const
1458 {
1459  // HKM -> Assumption is made that the solvent is species 0.
1460  if (m_indexSolvent != 0) {
1461  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
1462  "Wrong index solvent value!");
1463  }
1464 
1465  // Use the CROPPED molality of the species in solution.
1466  const vector_fp& molality = m_molalitiesCropped;
1467 
1468  // These are data inputs about the Pitzer correlation. They come from the
1469  // input file for the Pitzer model.
1470  vector_fp& gamma_Unscaled = m_gamma_tmp;
1471 
1472  // Local variables defined by Coltrin
1473  double etheta[5][5], etheta_prime[5][5], sqrtIs;
1474 
1475  // Molality based ionic strength of the solution
1476  double Is = 0.0;
1477 
1478  // Molarcharge of the solution: In Pitzer's notation, this is his variable
1479  // called "Z".
1480  double molarcharge = 0.0;
1481 
1482  // molalitysum is the sum of the molalities over all solutes, even those
1483  // with zero charge.
1484  double molalitysumUncropped = 0.0;
1485  debuglog("\n Debugging information from hmw_act \n", m_debugCalc);
1486 
1487  // Make sure the counter variables are setup
1488  counterIJ_setup();
1489 
1490  // ---------- Calculate common sums over solutes ---------------------
1491  for (size_t n = 1; n < m_kk; n++) {
1492  // ionic strength
1493  Is += charge(n) * charge(n) * molality[n];
1494  // total molar charge
1495  molarcharge += fabs(charge(n)) * molality[n];
1496  molalitysumUncropped += m_molalities[n];
1497  }
1498  Is *= 0.5;
1499 
1500  // Store the ionic molality in the object for reference.
1501  m_IionicMolality = Is;
1502  sqrtIs = sqrt(Is);
1503  if (m_debugCalc) {
1504  writelog(" Step 1: \n");
1505  writelogf(" ionic strenth = %14.7le \n total molar "
1506  "charge = %14.7le \n", Is, molarcharge);
1507  }
1508 
1509  // The following call to calc_lambdas() calculates all 16 elements of the
1510  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
1511  calc_lambdas(Is);
1512 
1513  // Step 2: Find the coefficients E-theta and E-thetaprime for all
1514  // combinations of positive unlike charges up to 4
1515  debuglog(" Step 2: \n", m_debugCalc);
1516  for (int z1 = 1; z1 <=4; z1++) {
1517  for (int z2 =1; z2 <=4; z2++) {
1518  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
1519  if (m_debugCalc) {
1520  writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
1521  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
1522  }
1523  }
1524  }
1525 
1526  debuglog(" Step 3: \n"
1527  " Species Species g(x) hfunc(x)\n",
1528  m_debugCalc);
1529 
1530  // calculate g(x) and hfunc(x) for each cation-anion pair MX. In the
1531  // original literature, hfunc, was called gprime. However, it's not the
1532  // derivative of g(x), so I renamed it.
1533  for (size_t i = 1; i < (m_kk - 1); i++) {
1534  for (size_t j = (i+1); j < m_kk; j++) {
1535  // Find the counterIJ for the symmetric binary interaction
1536  size_t n = m_kk*i + j;
1537  size_t counterIJ = m_CounterIJ[n];
1538 
1539  // Only loop over oppositely charge species
1540  if (charge(i)*charge(j) < 0) {
1541  // x is a reduced function variable
1542  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
1543  if (x1 > 1.0E-100) {
1544  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
1545  m_hfunc_IJ[counterIJ] = -2.0 *
1546  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
1547  } else {
1548  m_gfunc_IJ[counterIJ] = 0.0;
1549  m_hfunc_IJ[counterIJ] = 0.0;
1550  }
1551 
1552  if (m_Beta2MX_ij[counterIJ] != 0.0) {
1553  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
1554  if (x2 > 1.0E-100) {
1555  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
1556  m_h2func_IJ[counterIJ] = -2.0 *
1557  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
1558  } else {
1559  m_g2func_IJ[counterIJ] = 0.0;
1560  m_h2func_IJ[counterIJ] = 0.0;
1561  }
1562  }
1563  } else {
1564  m_gfunc_IJ[counterIJ] = 0.0;
1565  m_hfunc_IJ[counterIJ] = 0.0;
1566  }
1567  if (m_debugCalc) {
1568  writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i),
1569  speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
1570  }
1571  }
1572  }
1573 
1574  // SUBSECTION TO CALCULATE BMX, BprimeMX, BphiMX
1575  // Agrees with Pitzer, Eq. (49), (51), (55)
1576  debuglog(" Step 4: \n"
1577  " Species Species BMX BprimeMX BphiMX\n",
1578  m_debugCalc);
1579 
1580  for (size_t i = 1; i < m_kk - 1; i++) {
1581  for (size_t j = i+1; j < m_kk; j++) {
1582  // Find the counterIJ for the symmetric binary interaction
1583  size_t n = m_kk*i + j;
1584  size_t counterIJ = m_CounterIJ[n];
1585 
1586  // both species have a non-zero charge, and one is positive and the
1587  // other is negative
1588  if (charge(i)*charge(j) < 0.0) {
1589  m_BMX_IJ[counterIJ] = m_Beta0MX_ij[counterIJ]
1590  + m_Beta1MX_ij[counterIJ] * m_gfunc_IJ[counterIJ]
1591  + m_Beta2MX_ij[counterIJ] * m_g2func_IJ[counterIJ];
1592 
1593  if (m_debugCalc) {
1594  writelogf("%d %g: %g %g %g %g\n",
1595  counterIJ, m_BMX_IJ[counterIJ], m_Beta0MX_ij[counterIJ],
1596  m_Beta1MX_ij[counterIJ], m_Beta2MX_ij[counterIJ], m_gfunc_IJ[counterIJ]);
1597  }
1598  if (Is > 1.0E-150) {
1599  m_BprimeMX_IJ[counterIJ] = (m_Beta1MX_ij[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
1600  m_Beta2MX_ij[counterIJ] * m_h2func_IJ[counterIJ]/Is);
1601  } else {
1602  m_BprimeMX_IJ[counterIJ] = 0.0;
1603  }
1604  m_BphiMX_IJ[counterIJ] = m_BMX_IJ[counterIJ] + Is*m_BprimeMX_IJ[counterIJ];
1605  } else {
1606  m_BMX_IJ[counterIJ] = 0.0;
1607  m_BprimeMX_IJ[counterIJ] = 0.0;
1608  m_BphiMX_IJ[counterIJ] = 0.0;
1609  }
1610  if (m_debugCalc) {
1611  writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
1612  speciesName(i), speciesName(j),
1613  m_BMX_IJ[counterIJ], m_BprimeMX_IJ[counterIJ], m_BphiMX_IJ[counterIJ]);
1614  }
1615  }
1616  }
1617 
1618  // SUBSECTION TO CALCULATE CMX
1619  // Agrees with Pitzer, Eq. (53).
1620  debuglog(" Step 5: \n"
1621  " Species Species CMX\n", m_debugCalc);
1622  for (size_t i = 1; i < m_kk-1; i++) {
1623  for (size_t j = i+1; j < m_kk; j++) {
1624  // Find the counterIJ for the symmetric binary interaction
1625  size_t n = m_kk*i + j;
1626  size_t counterIJ = m_CounterIJ[n];
1627 
1628  // both species have a non-zero charge, and one is positive
1629  // and the other is negative
1630  if (charge(i)*charge(j) < 0.0) {
1631  m_CMX_IJ[counterIJ] = m_CphiMX_ij[counterIJ]/
1632  (2.0* sqrt(fabs(charge(i)*charge(j))));
1633  } else {
1634  m_CMX_IJ[counterIJ] = 0.0;
1635  }
1636  if (m_debugCalc) {
1637  writelogf(" %-16s %-16s %11.7f \n",
1638  speciesName(i), speciesName(j), m_CMX_IJ[counterIJ]);
1639  }
1640  }
1641  }
1642 
1643  // SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi
1644  // Agrees with Pitzer, Eq. 72, 73, 74
1645  debuglog(" Step 6: \n"
1646  " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
1647  m_debugCalc);
1648  for (size_t i = 1; i < m_kk-1; i++) {
1649  for (size_t j = i+1; j < m_kk; j++) {
1650  // Find the counterIJ for the symmetric binary interaction
1651  size_t n = m_kk*i + j;
1652  size_t counterIJ = m_CounterIJ[n];
1653 
1654  // both species have a non-zero charge, and one is positive and the
1655  // other is negative
1656  if (charge(i)*charge(j) > 0) {
1657  int z1 = (int) fabs(charge(i));
1658  int z2 = (int) fabs(charge(j));
1659  m_Phi_IJ[counterIJ] = m_Theta_ij[counterIJ] + etheta[z1][z2];
1660  m_Phiprime_IJ[counterIJ] = etheta_prime[z1][z2];
1661  m_PhiPhi_IJ[counterIJ] = m_Phi_IJ[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
1662  } else {
1663  m_Phi_IJ[counterIJ] = 0.0;
1664  m_Phiprime_IJ[counterIJ] = 0.0;
1665  m_PhiPhi_IJ[counterIJ] = 0.0;
1666  }
1667  if (m_debugCalc) {
1668  writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
1669  speciesName(i), speciesName(j),
1670  m_Phi_IJ[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ[counterIJ]);
1671  }
1672  }
1673  }
1674 
1675  // SUBSECTION FOR CALCULATION OF F
1676  // Agrees with Pitzer Eqn. (65)
1677  debuglog(" Step 7: \n", m_debugCalc);
1678  double Aphi = A_Debye_TP() / 3.0;
1679  double F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
1680  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
1681  if (m_debugCalc) {
1682  writelogf(" initial value of F = %10.6f \n", F);
1683  }
1684  for (size_t i = 1; i < m_kk-1; i++) {
1685  for (size_t j = i+1; j < m_kk; j++) {
1686  // Find the counterIJ for the symmetric binary interaction
1687  size_t n = m_kk*i + j;
1688  size_t counterIJ = m_CounterIJ[n];
1689 
1690  // both species have a non-zero charge, and one is positive and the
1691  // other is negative
1692  if (charge(i)*charge(j) < 0) {
1693  F += molality[i]*molality[j] * m_BprimeMX_IJ[counterIJ];
1694  }
1695 
1696  // Both species have a non-zero charge, and they
1697  // have the same sign
1698  if (charge(i)*charge(j) > 0) {
1699  F += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
1700  }
1701  if (m_debugCalc) {
1702  writelogf(" F = %10.6f \n", F);
1703  }
1704  }
1705  }
1706  debuglog(" Step 8: Summing in All Contributions to Activity Coefficients \n",
1707  m_debugCalc);
1708 
1709  for (size_t i = 1; i < m_kk; i++) {
1710 
1711  // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR CATIONS
1712  // equations agree with my notes, Eqn. (118).
1713  // Equations agree with Pitzer, eqn.(63)
1714  if (charge(i) > 0.0) {
1715  if (m_debugCalc) {
1716  writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
1717  }
1718  // species i is the cation (positive) to calc the actcoeff
1719  double zsqF = charge(i)*charge(i)*F;
1720  if (m_debugCalc) {
1721  writelogf(" Unary term: z*z*F = %10.5f\n", zsqF);
1722  }
1723  double sum1 = 0.0;
1724  double sum2 = 0.0;
1725  double sum3 = 0.0;
1726  double sum4 = 0.0;
1727  double sum5 = 0.0;
1728  for (size_t j = 1; j < m_kk; j++) {
1729  // Find the counterIJ for the symmetric binary interaction
1730  size_t n = m_kk*i + j;
1731  size_t counterIJ = m_CounterIJ[n];
1732 
1733  if (charge(j) < 0.0) {
1734  // sum over all anions
1735  sum1 += molality[j] *
1736  (2.0*m_BMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
1737  if (m_debugCalc) {
1738  std::string snj = speciesName(j) + ":";
1739  writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj,
1740  molality[j]*2.0*m_BMX_IJ[counterIJ]);
1741  writelogf(" m_j Z CMX = %10.5f\n",
1742  molality[j]* molarcharge*m_CMX_IJ[counterIJ]);
1743  }
1744  if (j < m_kk-1) {
1745  // This term is the ternary interaction involving the
1746  // non-duplicate sum over double anions, j, k, with
1747  // respect to the cation, i.
1748  for (size_t k = j+1; k < m_kk; k++) {
1749  // an inner sum over all anions
1750  if (charge(k) < 0.0) {
1751  n = k + j * m_kk + i * m_kk * m_kk;
1752  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
1753  if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
1754  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1755  writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
1756  molality[j]*molality[k]*m_Psi_ijk[n]);
1757  }
1758  }
1759  }
1760  }
1761  }
1762 
1763  if (charge(j) > 0.0) {
1764  // sum over all cations
1765  if (j != i) {
1766  sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
1767  if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
1768  std::string snj = speciesName(j) + ":";
1769  writelogf(" Phi term with %-12s 2 m_j Phi_cc = %10.5f\n", snj,
1770  molality[j]*(2.0*m_Phi_IJ[counterIJ]));
1771  }
1772  }
1773  for (size_t k = 1; k < m_kk; k++) {
1774  if (charge(k) < 0.0) {
1775  // two inner sums over anions
1776  n = k + j * m_kk + i * m_kk * m_kk;
1777  sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
1778  if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
1779  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1780  writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
1781  molality[j]*molality[k]*m_Psi_ijk[n]);
1782  }
1783 
1784  // Find the counterIJ for the j,k interaction
1785  n = m_kk*j + k;
1786  size_t counterIJ2 = m_CounterIJ[n];
1787  sum4 += (fabs(charge(i))*
1788  molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
1789  if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
1790  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1791  writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj,
1792  fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
1793  }
1794  }
1795  }
1796  }
1797 
1798  // Handle neutral j species
1799  if (charge(j) == 0) {
1800  sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
1801  if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
1802  std::string snj = speciesName(j) + ":";
1803  writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj,
1804  molality[j]*2.0*m_Lambda_nj(j,i));
1805  }
1806 
1807  // Zeta interaction term
1808  for (size_t k = 1; k < m_kk; k++) {
1809  if (charge(k) < 0.0) {
1810  size_t izeta = j;
1811  size_t jzeta = i;
1812  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
1813  double zeta = m_Psi_ijk[n];
1814  if (zeta != 0.0) {
1815  sum5 += molality[j]*molality[k]*zeta;
1816  if (m_debugCalc) {
1817  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1818  writelogf(" Zeta term on %-16s m_n m_a zeta_nMa = %10.5f\n", snj,
1819  molality[j]*molality[k]*m_Psi_ijk[n]);
1820  }
1821  }
1822  }
1823  }
1824  }
1825  }
1826 
1827  // Add all of the contributions up to yield the log of the solute
1828  // activity coefficients (molality scale)
1829  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
1830  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
1831  if (m_debugCalc) {
1832  writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f \n",
1833  speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
1834  }
1835  }
1836 
1837  // SUBSECTION FOR CALCULATING THE ACTCOEFF FOR ANIONS
1838  // equations agree with my notes, Eqn. (119).
1839  // Equations agree with Pitzer, eqn.(64)
1840  if (charge(i) < 0) {
1841  if (m_debugCalc) {
1842  writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
1843  }
1844  // species i is an anion (negative)
1845  double zsqF = charge(i)*charge(i)*F;
1846  if (m_debugCalc) {
1847  writelogf(" Unary term: z*z*F = %10.5f\n", zsqF);
1848  }
1849  double sum1 = 0.0;
1850  double sum2 = 0.0;
1851  double sum3 = 0.0;
1852  double sum4 = 0.0;
1853  double sum5 = 0.0;
1854  for (size_t j = 1; j < m_kk; j++) {
1855  // Find the counterIJ for the symmetric binary interaction
1856  size_t n = m_kk*i + j;
1857  size_t counterIJ = m_CounterIJ[n];
1858 
1859  // For Anions, do the cation interactions.
1860  if (charge(j) > 0) {
1861  sum1 += molality[j]*
1862  (2.0*m_BMX_IJ[counterIJ]+molarcharge*m_CMX_IJ[counterIJ]);
1863  if (m_debugCalc) {
1864  std::string snj = speciesName(j) + ":";
1865  writelogf(" Bin term with %-13s 2 m_j BMX = %10.5f\n", snj,
1866  molality[j]*2.0*m_BMX_IJ[counterIJ]);
1867  writelogf(" m_j Z CMX = %10.5f\n",
1868  molality[j]* molarcharge*m_CMX_IJ[counterIJ]);
1869  }
1870  if (j < m_kk-1) {
1871  for (size_t k = j+1; k < m_kk; k++) {
1872  // an inner sum over all cations
1873  if (charge(k) > 0) {
1874  n = k + j * m_kk + i * m_kk * m_kk;
1875  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
1876  if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
1877  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1878  writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
1879  molality[j]*molality[k]*m_Psi_ijk[n]);
1880  }
1881  }
1882  }
1883  }
1884  }
1885 
1886  // For Anions, do the other anion interactions.
1887  if (charge(j) < 0.0) {
1888  // sum over all anions
1889  if (j != i) {
1890  sum2 += molality[j]*(2.0*m_Phi_IJ[counterIJ]);
1891  if (m_debugCalc && (molality[j] * m_Phi_IJ[counterIJ])!= 0.0) {
1892  std::string snj = speciesName(j) + ":";
1893  writelogf(" Phi term with %-12s 2 m_j Phi_aa = %10.5f\n", snj,
1894  molality[j]*(2.0*m_Phi_IJ[counterIJ]));
1895  }
1896  }
1897  for (size_t k = 1; k < m_kk; k++) {
1898  if (charge(k) > 0.0) {
1899  // two inner sums over cations
1900  n = k + j * m_kk + i * m_kk * m_kk;
1901  sum2 += molality[j]*molality[k]*m_Psi_ijk[n];
1902  if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
1903  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1904  writelogf(" Psi term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
1905  molality[j]*molality[k]*m_Psi_ijk[n]);
1906  }
1907  // Find the counterIJ for the symmetric binary interaction
1908  n = m_kk*j + k;
1909  size_t counterIJ2 = m_CounterIJ[n];
1910  sum4 += fabs(charge(i))*
1911  molality[j]*molality[k]*m_CMX_IJ[counterIJ2];
1912  if (m_debugCalc && (molality[j]*molality[k]*m_CMX_IJ[counterIJ2]) != 0.0) {
1913  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1914  writelogf(" Tern CMX term on %-16s abs(z_i) m_j m_k CMX = %10.5f\n", snj,
1915  fabs(charge(i))* molality[j]*molality[k]*m_CMX_IJ[counterIJ2]);
1916  }
1917  }
1918  }
1919  }
1920 
1921  // for Anions, do the neutral species interaction
1922  if (charge(j) == 0.0) {
1923  sum5 += molality[j]*2.0*m_Lambda_nj(j,i);
1924  if (m_debugCalc && (molality[j]*2.0*m_Lambda_nj(j,i)) != 0.0) {
1925  std::string snj = speciesName(j) + ":";
1926  writelogf(" Lambda term with %-12s 2 m_j lam_ji = %10.5f\n", snj,
1927  molality[j]*2.0*m_Lambda_nj(j,i));
1928  }
1929  // Zeta interaction term
1930  for (size_t k = 1; k < m_kk; k++) {
1931  if (charge(k) > 0.0) {
1932  size_t izeta = j;
1933  size_t jzeta = k;
1934  size_t kzeta = i;
1935  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
1936  double zeta = m_Psi_ijk[n];
1937  if (zeta != 0.0) {
1938  sum5 += molality[j]*molality[k]*zeta;
1939  if (m_debugCalc) {
1940  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1941  writelogf(" Zeta term on %-16s m_n m_c zeta_ncX = %10.5f\n", snj,
1942  molality[j]*molality[k]*m_Psi_ijk[n]);
1943  }
1944  }
1945  }
1946  }
1947  }
1948  }
1949  m_lnActCoeffMolal_Unscaled[i] = zsqF + sum1 + sum2 + sum3 + sum4 + sum5;
1950  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
1951  if (m_debugCalc) {
1952  writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
1953  speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
1954  }
1955  }
1956 
1957  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
1958  // equations agree with my notes,
1959  // Equations agree with Pitzer,
1960  if (charge(i) == 0.0) {
1961  if (m_debugCalc) {
1962  writelogf(" Contributions to ln(ActCoeff_%s):\n", speciesName(i));
1963  }
1964  double sum1 = 0.0;
1965  double sum3 = 0.0;
1966  for (size_t j = 1; j < m_kk; j++) {
1967  sum1 += molality[j]*2.0*m_Lambda_nj(i,j);
1968  if (m_debugCalc && m_Lambda_nj(i,j) != 0.0) {
1969  std::string snj = speciesName(j) + ":";
1970  writelogf(" Lambda_n term on %-16s 2 m_j lambda_n_j = %10.5f\n", snj,
1971  molality[j]*2.0*m_Lambda_nj(i,j));
1972  }
1973  // Zeta term -> we piggyback on the psi term
1974  if (charge(j) > 0.0) {
1975  for (size_t k = 1; k < m_kk; k++) {
1976  if (charge(k) < 0.0) {
1977  size_t n = k + j * m_kk + i * m_kk * m_kk;
1978  sum3 += molality[j]*molality[k]*m_Psi_ijk[n];
1979  if (m_debugCalc && m_Psi_ijk[n] != 0.0) {
1980  std::string snj = speciesName(j) + "," + speciesName(k) + ":";
1981  writelogf(" Zeta term on %-16s m_j m_k psi_ijk = %10.5f\n", snj,
1982  molality[j]*molality[k]*m_Psi_ijk[n]);
1983  }
1984  }
1985  }
1986  }
1987  }
1988  double sum2 = 3.0 * molality[i]* molality[i] * m_Mu_nnn[i];
1989  if (m_debugCalc && m_Mu_nnn[i] != 0.0) {
1990  writelogf(" Mu_nnn term 3 m_n m_n Mu_n_n = %10.5f\n",
1991  3.0 * molality[i]* molality[i] * m_Mu_nnn[i]);
1992  }
1993  m_lnActCoeffMolal_Unscaled[i] = sum1 + sum2 + sum3;
1994  gamma_Unscaled[i] = exp(m_lnActCoeffMolal_Unscaled[i]);
1995  if (m_debugCalc) {
1996  writelogf(" Net %-16s lngamma[i] = %9.5f gamma[i]=%10.6f\n",
1997  speciesName(i), m_lnActCoeffMolal_Unscaled[i], gamma_Unscaled[i]);
1998  }
1999  }
2000  }
2001  debuglog(" Step 9: \n", m_debugCalc);
2002 
2003  // SUBSECTION FOR CALCULATING THE OSMOTIC COEFF
2004  // equations agree with my notes, Eqn. (117).
2005  // Equations agree with Pitzer, eqn.(62)
2006  double sum1 = 0.0;
2007  double sum2 = 0.0;
2008  double sum3 = 0.0;
2009  double sum4 = 0.0;
2010  double sum5 = 0.0;
2011  double sum6 = 0.0;
2012  double sum7 = 0.0;
2013 
2014  // term1 is the DH term in the osmotic coefficient expression
2015  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer
2016  // implementations.
2017  // Is = Ionic strength on the molality scale (units of (gmol/kg))
2018  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2019  double term1 = -Aphi * pow(Is,1.5) / (1.0 + 1.2 * sqrt(Is));
2020 
2021  for (size_t j = 1; j < m_kk; j++) {
2022  // Loop Over Cations
2023  if (charge(j) > 0.0) {
2024  for (size_t k = 1; k < m_kk; k++) {
2025  if (charge(k) < 0.0) {
2026  // Find the counterIJ for the symmetric j,k binary interaction
2027  size_t n = m_kk*j + k;
2028  size_t counterIJ = m_CounterIJ[n];
2029 
2030  sum1 += molality[j]*molality[k]*
2031  (m_BphiMX_IJ[counterIJ] + molarcharge*m_CMX_IJ[counterIJ]);
2032  }
2033  }
2034 
2035  for (size_t k = j+1; k < m_kk; k++) {
2036  if (j == (m_kk-1)) {
2037  // we should never reach this step
2038  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2039  "logic error 1 in Step 9 of hmw_act");
2040  }
2041  if (charge(k) > 0.0) {
2042  // Find the counterIJ for the symmetric j,k binary interaction
2043  // between 2 cations.
2044  size_t n = m_kk*j + k;
2045  size_t counterIJ = m_CounterIJ[n];
2046  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2047  for (size_t m = 1; m < m_kk; m++) {
2048  if (charge(m) < 0.0) {
2049  // species m is an anion
2050  n = m + k * m_kk + j * m_kk * m_kk;
2051  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2052  }
2053  }
2054  }
2055  }
2056  }
2057 
2058  // Loop Over Anions
2059  if (charge(j) < 0) {
2060  for (size_t k = j+1; k < m_kk; k++) {
2061  if (j == m_kk-1) {
2062  // we should never reach this step
2063  throw CanteraError("HMWSoln::s_updatePitzer_lnMolalityActCoeff",
2064  "logic error 2 in Step 9 of hmw_act");
2065  }
2066  if (charge(k) < 0) {
2067  // Find the counterIJ for the symmetric j,k binary interaction
2068  // between two anions
2069  size_t n = m_kk*j + k;
2070  size_t counterIJ = m_CounterIJ[n];
2071  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ[counterIJ];
2072  for (size_t m = 1; m < m_kk; m++) {
2073  if (charge(m) > 0.0) {
2074  n = m + k * m_kk + j * m_kk * m_kk;
2075  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk[n];
2076  }
2077  }
2078  }
2079  }
2080  }
2081 
2082  // Loop Over Neutral Species
2083  if (charge(j) == 0) {
2084  for (size_t k = 1; k < m_kk; k++) {
2085  if (charge(k) < 0.0) {
2086  sum4 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2087  }
2088  if (charge(k) > 0.0) {
2089  sum5 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2090  }
2091  if (charge(k) == 0.0) {
2092  if (k > j) {
2093  sum6 += molality[j]*molality[k]*m_Lambda_nj(j,k);
2094  } else if (k == j) {
2095  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj(j,k);
2096  }
2097  }
2098  if (charge(k) < 0.0) {
2099  size_t izeta = j;
2100  for (size_t m = 1; m < m_kk; m++) {
2101  if (charge(m) > 0.0) {
2102  size_t jzeta = m;
2103  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2104  double zeta = m_Psi_ijk[n];
2105  if (zeta != 0.0) {
2106  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta;
2107  }
2108  }
2109  }
2110  }
2111  }
2112  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn[j];
2113  }
2114  }
2115  double sum_m_phi_minus_1 = 2.0 *
2116  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2117  // Calculate the osmotic coefficient from
2118  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2119  double osmotic_coef;
2120  if (molalitysumUncropped > 1.0E-150) {
2121  osmotic_coef = 1.0 + (sum_m_phi_minus_1 / molalitysumUncropped);
2122  } else {
2123  osmotic_coef = 1.0;
2124  }
2125  if (m_debugCalc) {
2126  writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
2127  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
2128  term1, sum1, sum2, sum3, sum4, sum5);
2129  writelogf(" sum_m_phi_minus_1=%10.6f osmotic_coef=%10.6f\n",
2130  sum_m_phi_minus_1, osmotic_coef);
2131  writelog(" Step 10: \n");
2132  }
2133  double lnwateract = -(m_weightSolvent/1000.0) * molalitysumUncropped * osmotic_coef;
2134 
2135  // In Cantera, we define the activity coefficient of the solvent as
2136  //
2137  // act_0 = actcoeff_0 * Xmol_0
2138  //
2139  // We have just computed act_0. However, this routine returns
2140  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2141  double xmolSolvent = moleFraction(m_indexSolvent);
2142  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
2143  m_lnActCoeffMolal_Unscaled[0] = lnwateract - log(xx);
2144  if (m_debugCalc) {
2145  double wateract = exp(lnwateract);
2146  writelogf(" Weight of Solvent = %16.7g\n", m_weightSolvent);
2147  writelogf(" molalitySumUncropped = %16.7g\n", molalitysumUncropped);
2148  writelogf(" ln_a_water=%10.6f a_water=%10.6f\n\n",
2149  lnwateract, wateract);
2150  }
2151 }
2152 
2153 void HMWSoln::s_update_dlnMolalityActCoeff_dT() const
2154 {
2155  static const int cacheId = m_cache.getId();
2156  CachedScalar cached = m_cache.getScalar(cacheId);
2157  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2158  return;
2159  }
2160 
2161  // Zero the unscaled 2nd derivatives
2162  m_dlnActCoeffMolaldT_Unscaled.assign(m_kk, 0.0);
2163 
2164  // Do the actual calculation of the unscaled temperature derivatives
2165  s_updatePitzer_dlnMolalityActCoeff_dT();
2166 
2167  for (size_t k = 1; k < m_kk; k++) {
2168  if (CROP_speciesCropped_[k] == 2) {
2169  m_dlnActCoeffMolaldT_Unscaled[k] = 0.0;
2170  }
2171  }
2172 
2173  if (CROP_speciesCropped_[0]) {
2174  m_dlnActCoeffMolaldT_Unscaled[0] = 0.0;
2175  }
2176 
2177  // Do the pH scaling to the derivatives
2178  s_updateScaling_pHScaling_dT();
2179 }
2180 
2181 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const
2182 {
2183  // It may be assumed that the Pitzer activity coefficient routine is called
2184  // immediately preceding the calling of this routine. Therefore, some
2185  // quantities do not need to be recalculated in this routine.
2186 
2187  // HKM -> Assumption is made that the solvent is species 0.
2188  if (m_indexSolvent != 0) {
2189  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2190  "Wrong index solvent value!");
2191  }
2192 
2193  const vector_fp& molality = m_molalitiesCropped;
2194  double* d_gamma_dT_Unscaled = m_gamma_tmp.data();
2195 
2196  // Local variables defined by Coltrin
2197  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2198 
2199  // Molality based ionic strength of the solution
2200  double Is = 0.0;
2201 
2202  // Molarcharge of the solution: In Pitzer's notation, this is his variable
2203  // called "Z".
2204  double molarcharge = 0.0;
2205 
2206  // molalitysum is the sum of the molalities over all solutes, even those
2207  // with zero charge.
2208  double molalitysum = 0.0;
2209 
2210  debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dT()\n",
2211  m_debugCalc);
2212 
2213  // Make sure the counter variables are setup
2214  counterIJ_setup();
2215 
2216  // ---------- Calculate common sums over solutes ---------------------
2217  for (size_t n = 1; n < m_kk; n++) {
2218  // ionic strength
2219  Is += charge(n) * charge(n) * molality[n];
2220  // total molar charge
2221  molarcharge += fabs(charge(n)) * molality[n];
2222  molalitysum += molality[n];
2223  }
2224  Is *= 0.5;
2225 
2226  // Store the ionic molality in the object for reference.
2227  m_IionicMolality = Is;
2228  sqrtIs = sqrt(Is);
2229  if (m_debugCalc) {
2230  writelog(" Step 1: \n");
2231  writelogf(" ionic strenth = %14.7le \n total molar "
2232  "charge = %14.7le \n", Is, molarcharge);
2233  }
2234 
2235  // The following call to calc_lambdas() calculates all 16 elements of the
2236  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2237  calc_lambdas(Is);
2238 
2239  // Step 2: Find the coefficients E-theta and E-thetaprime for all
2240  // combinations of positive unlike charges up to 4
2241  debuglog(" Step 2: \n", m_debugCalc);
2242  for (int z1 = 1; z1 <=4; z1++) {
2243  for (int z2 =1; z2 <=4; z2++) {
2244  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2245  if (m_debugCalc) {
2246  writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2247  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2248  }
2249  }
2250  }
2251 
2252  debuglog(" Step 3: \n"
2253  " Species Species g(x) hfunc(x) \n",
2254  m_debugCalc);
2255 
2256  // calculate g(x) and hfunc(x) for each cation-anion pair MX
2257  // In the original literature, hfunc, was called gprime. However,
2258  // it's not the derivative of g(x), so I renamed it.
2259  for (size_t i = 1; i < (m_kk - 1); i++) {
2260  for (size_t j = (i+1); j < m_kk; j++) {
2261  // Find the counterIJ for the symmetric binary interaction
2262  size_t n = m_kk*i + j;
2263  size_t counterIJ = m_CounterIJ[n];
2264 
2265  // Only loop over oppositely charge species
2266  if (charge(i)*charge(j) < 0) {
2267  // x is a reduced function variable
2268  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2269  if (x1 > 1.0E-100) {
2270  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
2271  m_hfunc_IJ[counterIJ] = -2.0 *
2272  (1.0-(1.0 + x1 + 0.5 * x1 *x1) * exp(-x1)) / (x1 * x1);
2273  } else {
2274  m_gfunc_IJ[counterIJ] = 0.0;
2275  m_hfunc_IJ[counterIJ] = 0.0;
2276  }
2277 
2278  if (m_Beta2MX_ij_L[counterIJ] != 0.0) {
2279  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2280  if (x2 > 1.0E-100) {
2281  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2282  m_h2func_IJ[counterIJ] = -2.0 *
2283  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2284  } else {
2285  m_g2func_IJ[counterIJ] = 0.0;
2286  m_h2func_IJ[counterIJ] = 0.0;
2287  }
2288  }
2289  } else {
2290  m_gfunc_IJ[counterIJ] = 0.0;
2291  m_hfunc_IJ[counterIJ] = 0.0;
2292  }
2293  if (m_debugCalc) {
2294  std::string sni = speciesName(i);
2295  std::string snj = speciesName(j);
2296  writelogf(" %-16s %-16s %9.5f %9.5f \n", sni.c_str(), snj.c_str(),
2297  m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
2298  }
2299  }
2300  }
2301 
2302  // SUBSECTION TO CALCULATE BMX_L, BprimeMX_L, BphiMX_L
2303  // These are now temperature derivatives of the previously calculated
2304  // quantities.
2305  debuglog(" Step 4: \n"
2306  " Species Species BMX BprimeMX BphiMX \n",
2307  m_debugCalc);
2308 
2309  for (size_t i = 1; i < m_kk - 1; i++) {
2310  for (size_t j = i+1; j < m_kk; j++) {
2311  // Find the counterIJ for the symmetric binary interaction
2312  size_t n = m_kk*i + j;
2313  size_t counterIJ = m_CounterIJ[n];
2314 
2315  // both species have a non-zero charge, and one is positive
2316  // and the other is negative
2317  if (charge(i)*charge(j) < 0.0) {
2318  m_BMX_IJ_L[counterIJ] = m_Beta0MX_ij_L[counterIJ]
2319  + m_Beta1MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ]
2320  + m_Beta2MX_ij_L[counterIJ] * m_gfunc_IJ[counterIJ];
2321  if (m_debugCalc) {
2322  writelogf("%d %g: %g %g %g %g\n",
2323  counterIJ, m_BMX_IJ_L[counterIJ], m_Beta0MX_ij_L[counterIJ],
2324  m_Beta1MX_ij_L[counterIJ], m_Beta2MX_ij_L[counterIJ], m_gfunc_IJ[counterIJ]);
2325  }
2326  if (Is > 1.0E-150) {
2327  m_BprimeMX_IJ_L[counterIJ] = (m_Beta1MX_ij_L[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2328  m_Beta2MX_ij_L[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2329  } else {
2330  m_BprimeMX_IJ_L[counterIJ] = 0.0;
2331  }
2332  m_BphiMX_IJ_L[counterIJ] = m_BMX_IJ_L[counterIJ] + Is*m_BprimeMX_IJ_L[counterIJ];
2333  } else {
2334  m_BMX_IJ_L[counterIJ] = 0.0;
2335  m_BprimeMX_IJ_L[counterIJ] = 0.0;
2336  m_BphiMX_IJ_L[counterIJ] = 0.0;
2337  }
2338  if (m_debugCalc) {
2339  writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
2340  speciesName(i), speciesName(j),
2341  m_BMX_IJ_L[counterIJ], m_BprimeMX_IJ_L[counterIJ], m_BphiMX_IJ_L[counterIJ]);
2342  }
2343  }
2344  }
2345 
2346  // --------- SUBSECTION TO CALCULATE CMX_L ----------
2347  debuglog(" Step 5: \n"
2348  " Species Species CMX \n", m_debugCalc);
2349  for (size_t i = 1; i < m_kk-1; i++) {
2350  for (size_t j = i+1; j < m_kk; j++) {
2351  // Find the counterIJ for the symmetric binary interaction
2352  size_t n = m_kk*i + j;
2353  size_t counterIJ = m_CounterIJ[n];
2354 
2355  // both species have a non-zero charge, and one is positive
2356  // and the other is negative
2357  if (charge(i)*charge(j) < 0.0) {
2358  m_CMX_IJ_L[counterIJ] = m_CphiMX_ij_L[counterIJ]/
2359  (2.0* sqrt(fabs(charge(i)*charge(j))));
2360  } else {
2361  m_CMX_IJ_L[counterIJ] = 0.0;
2362  }
2363  if (m_debugCalc) {
2364  writelogf(" %-16s %-16s %11.7f \n",
2365  speciesName(i), speciesName(j), m_CMX_IJ_L[counterIJ]);
2366  }
2367  }
2368  }
2369 
2370  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2371  debuglog(" Step 6: \n"
2372  " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
2373  m_debugCalc);
2374  for (size_t i = 1; i < m_kk-1; i++) {
2375  for (size_t j = i+1; j < m_kk; j++) {
2376  // Find the counterIJ for the symmetric binary interaction
2377  size_t n = m_kk*i + j;
2378  size_t counterIJ = m_CounterIJ[n];
2379 
2380  // both species have a non-zero charge, and one is positive
2381  // and the other is negative
2382  if (charge(i)*charge(j) > 0) {
2383  m_Phi_IJ_L[counterIJ] = m_Theta_ij_L[counterIJ];
2384  m_Phiprime_IJ[counterIJ] = 0.0;
2385  m_PhiPhi_IJ_L[counterIJ] = m_Phi_IJ_L[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
2386  } else {
2387  m_Phi_IJ_L[counterIJ] = 0.0;
2388  m_Phiprime_IJ[counterIJ] = 0.0;
2389  m_PhiPhi_IJ_L[counterIJ] = 0.0;
2390  }
2391  if (m_debugCalc) {
2392  writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
2393  speciesName(i), speciesName(j),
2394  m_Phi_IJ_L[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_L[counterIJ]);
2395  }
2396  }
2397  }
2398 
2399  // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
2400  debuglog(" Step 7: \n", m_debugCalc);
2401  double dA_DebyedT = dA_DebyedT_TP();
2402  double dAphidT = dA_DebyedT /3.0;
2403  double dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
2404  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
2405  if (m_debugCalc) {
2406  writelogf(" initial value of dFdT = %10.6f \n", dFdT);
2407  }
2408  for (size_t i = 1; i < m_kk-1; i++) {
2409  for (size_t j = i+1; j < m_kk; j++) {
2410  // Find the counterIJ for the symmetric binary interaction
2411  size_t n = m_kk*i + j;
2412  size_t counterIJ = m_CounterIJ[n];
2413 
2414  // both species have a non-zero charge, and one is positive
2415  // and the other is negative
2416  if (charge(i)*charge(j) < 0) {
2417  dFdT += molality[i]*molality[j] * m_BprimeMX_IJ_L[counterIJ];
2418  }
2419 
2420  // Both species have a non-zero charge, and they
2421  // have the same sign, e.g., both positive or both negative.
2422  if (charge(i)*charge(j) > 0) {
2423  dFdT += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
2424  }
2425  if (m_debugCalc) {
2426  writelogf(" dFdT = %10.6f \n", dFdT);
2427  }
2428  }
2429  }
2430  debuglog(" Step 8: \n", m_debugCalc);
2431 
2432  for (size_t i = 1; i < m_kk; i++) {
2433  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
2434  if (charge(i) > 0) {
2435  // species i is the cation (positive) to calc the actcoeff
2436  double zsqdFdT = charge(i)*charge(i)*dFdT;
2437  double sum1 = 0.0;
2438  double sum2 = 0.0;
2439  double sum3 = 0.0;
2440  double sum4 = 0.0;
2441  double sum5 = 0.0;
2442  for (size_t j = 1; j < m_kk; j++) {
2443  // Find the counterIJ for the symmetric binary interaction
2444  size_t n = m_kk*i + j;
2445  size_t counterIJ = m_CounterIJ[n];
2446 
2447  if (charge(j) < 0.0) {
2448  // sum over all anions
2449  sum1 += molality[j]*
2450  (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2451  if (j < m_kk-1) {
2452  // This term is the ternary interaction involving the
2453  // non-duplicate sum over double anions, j, k, with
2454  // respect to the cation, i.
2455  for (size_t k = j+1; k < m_kk; k++) {
2456  // an inner sum over all anions
2457  if (charge(k) < 0.0) {
2458  n = k + j * m_kk + i * m_kk * m_kk;
2459  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2460  }
2461  }
2462  }
2463  }
2464 
2465  if (charge(j) > 0.0) {
2466  // sum over all cations
2467  if (j != i) {
2468  sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2469  }
2470  for (size_t k = 1; k < m_kk; k++) {
2471  if (charge(k) < 0.0) {
2472  // two inner sums over anions
2473  n = k + j * m_kk + i * m_kk * m_kk;
2474  sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2475 
2476  // Find the counterIJ for the j,k interaction
2477  n = m_kk*j + k;
2478  size_t counterIJ2 = m_CounterIJ[n];
2479  sum4 += fabs(charge(i))*
2480  molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2481  }
2482  }
2483  }
2484 
2485  // Handle neutral j species
2486  if (charge(j) == 0) {
2487  sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2488  }
2489 
2490  // Zeta interaction term
2491  for (size_t k = 1; k < m_kk; k++) {
2492  if (charge(k) < 0.0) {
2493  size_t izeta = j;
2494  size_t jzeta = i;
2495  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
2496  double zeta_L = m_Psi_ijk_L[n];
2497  if (zeta_L != 0.0) {
2498  sum5 += molality[j]*molality[k]*zeta_L;
2499  }
2500  }
2501  }
2502  }
2503 
2504  // Add all of the contributions up to yield the log of the
2505  // solute activity coefficients (molality scale)
2506  m_dlnActCoeffMolaldT_Unscaled[i] =
2507  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2508  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2509  if (m_debugCalc) {
2510  writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
2511  speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
2512  writelogf(" %12g %12g %12g %12g %12g %12g\n",
2513  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
2514  }
2515  }
2516 
2517  // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR ANIONS ------
2518  if (charge(i) < 0) {
2519  // species i is an anion (negative)
2520  double zsqdFdT = charge(i)*charge(i)*dFdT;
2521  double sum1 = 0.0;
2522  double sum2 = 0.0;
2523  double sum3 = 0.0;
2524  double sum4 = 0.0;
2525  double sum5 = 0.0;
2526  for (size_t j = 1; j < m_kk; j++) {
2527  // Find the counterIJ for the symmetric binary interaction
2528  size_t n = m_kk*i + j;
2529  size_t counterIJ = m_CounterIJ[n];
2530 
2531  // For Anions, do the cation interactions.
2532  if (charge(j) > 0) {
2533  sum1 += molality[j]*
2534  (2.0*m_BMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2535  if (j < m_kk-1) {
2536  for (size_t k = j+1; k < m_kk; k++) {
2537  // an inner sum over all cations
2538  if (charge(k) > 0) {
2539  n = k + j * m_kk + i * m_kk * m_kk;
2540  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2541  }
2542  }
2543  }
2544  }
2545 
2546  // For Anions, do the other anion interactions.
2547  if (charge(j) < 0.0) {
2548  // sum over all anions
2549  if (j != i) {
2550  sum2 += molality[j]*(2.0*m_Phi_IJ_L[counterIJ]);
2551  }
2552  for (size_t k = 1; k < m_kk; k++) {
2553  if (charge(k) > 0.0) {
2554  // two inner sums over cations
2555  n = k + j * m_kk + i * m_kk * m_kk;
2556  sum2 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2557  // Find the counterIJ for the symmetric binary interaction
2558  n = m_kk*j + k;
2559  size_t counterIJ2 = m_CounterIJ[n];
2560  sum4 += fabs(charge(i)) *
2561  molality[j]*molality[k]*m_CMX_IJ_L[counterIJ2];
2562  }
2563  }
2564  }
2565 
2566  // for Anions, do the neutral species interaction
2567  if (charge(j) == 0.0) {
2568  sum5 += molality[j]*2.0*m_Lambda_nj_L(j,i);
2569  for (size_t k = 1; k < m_kk; k++) {
2570  if (charge(k) > 0.0) {
2571  size_t izeta = j;
2572  size_t jzeta = k;
2573  size_t kzeta = i;
2574  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
2575  double zeta_L = m_Psi_ijk_L[n];
2576  if (zeta_L != 0.0) {
2577  sum5 += molality[j]*molality[k]*zeta_L;
2578  }
2579  }
2580  }
2581  }
2582  }
2583  m_dlnActCoeffMolaldT_Unscaled[i] =
2584  zsqdFdT + sum1 + sum2 + sum3 + sum4 + sum5;
2585  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2586  if (m_debugCalc) {
2587  writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f\n",
2588  speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
2589  writelogf(" %12g %12g %12g %12g %12g %12g\n",
2590  zsqdFdT, sum1, sum2, sum3, sum4, sum5);
2591  }
2592  }
2593 
2594  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
2595  // equations agree with my notes,
2596  // Equations agree with Pitzer,
2597  if (charge(i) == 0.0) {
2598  double sum1 = 0.0;
2599  double sum3 = 0.0;
2600  for (size_t j = 1; j < m_kk; j++) {
2601  sum1 += molality[j]*2.0*m_Lambda_nj_L(i,j);
2602  // Zeta term -> we piggyback on the psi term
2603  if (charge(j) > 0.0) {
2604  for (size_t k = 1; k < m_kk; k++) {
2605  if (charge(k) < 0.0) {
2606  size_t n = k + j * m_kk + i * m_kk * m_kk;
2607  sum3 += molality[j]*molality[k]*m_Psi_ijk_L[n];
2608  }
2609  }
2610  }
2611  }
2612  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_L[i];
2613  m_dlnActCoeffMolaldT_Unscaled[i] = sum1 + sum2 + sum3;
2614  d_gamma_dT_Unscaled[i] = exp(m_dlnActCoeffMolaldT_Unscaled[i]);
2615  if (m_debugCalc) {
2616  writelogf(" %-16s lngamma[i]=%10.6f gamma[i]=%10.6f \n",
2617  speciesName(i), m_dlnActCoeffMolaldT_Unscaled[i], d_gamma_dT_Unscaled[i]);
2618  }
2619  }
2620  }
2621  debuglog(" Step 9: \n", m_debugCalc);
2622 
2623  // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dT ---------
2624  double sum1 = 0.0;
2625  double sum2 = 0.0;
2626  double sum3 = 0.0;
2627  double sum4 = 0.0;
2628  double sum5 = 0.0;
2629  double sum6 = 0.0;
2630  double sum7 = 0.0;
2631 
2632  // term1 is the temperature derivative of the DH term in the osmotic
2633  // coefficient expression
2634  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
2635  // Is = Ionic strength on the molality scale (units of (gmol/kg))
2636  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
2637  double term1 = -dAphidT * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
2638 
2639  for (size_t j = 1; j < m_kk; j++) {
2640  // Loop Over Cations
2641  if (charge(j) > 0.0) {
2642  for (size_t k = 1; k < m_kk; k++) {
2643  if (charge(k) < 0.0) {
2644  // Find the counterIJ for the symmetric j,k binary interaction
2645  size_t n = m_kk*j + k;
2646  size_t counterIJ = m_CounterIJ[n];
2647  sum1 += molality[j]*molality[k]*
2648  (m_BphiMX_IJ_L[counterIJ] + molarcharge*m_CMX_IJ_L[counterIJ]);
2649  }
2650  }
2651 
2652  for (size_t k = j+1; k < m_kk; k++) {
2653  if (j == (m_kk-1)) {
2654  // we should never reach this step
2655  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2656  "logic error 1 in Step 9 of hmw_act");
2657  }
2658  if (charge(k) > 0.0) {
2659  // Find the counterIJ for the symmetric j,k binary interaction
2660  // between 2 cations.
2661  size_t n = m_kk*j + k;
2662  size_t counterIJ = m_CounterIJ[n];
2663  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2664  for (size_t m = 1; m < m_kk; m++) {
2665  if (charge(m) < 0.0) {
2666  // species m is an anion
2667  n = m + k * m_kk + j * m_kk * m_kk;
2668  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2669  }
2670  }
2671  }
2672  }
2673  }
2674 
2675  // Loop Over Anions
2676  if (charge(j) < 0) {
2677  for (size_t k = j+1; k < m_kk; k++) {
2678  if (j == m_kk-1) {
2679  // we should never reach this step
2680  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT",
2681  "logic error 2 in Step 9 of hmw_act");
2682  }
2683  if (charge(k) < 0) {
2684  // Find the counterIJ for the symmetric j,k binary interaction
2685  // between two anions
2686  size_t n = m_kk*j + k;
2687  size_t counterIJ = m_CounterIJ[n];
2688  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_L[counterIJ];
2689  for (size_t m = 1; m < m_kk; m++) {
2690  if (charge(m) > 0.0) {
2691  n = m + k * m_kk + j * m_kk * m_kk;
2692  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_L[n];
2693  }
2694  }
2695  }
2696  }
2697  }
2698 
2699  // Loop Over Neutral Species
2700  if (charge(j) == 0) {
2701  for (size_t k = 1; k < m_kk; k++) {
2702  if (charge(k) < 0.0) {
2703  sum4 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2704  }
2705  if (charge(k) > 0.0) {
2706  sum5 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2707  }
2708  if (charge(k) == 0.0) {
2709  if (k > j) {
2710  sum6 += molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2711  } else if (k == j) {
2712  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_L(j,k);
2713  }
2714  }
2715  if (charge(k) < 0.0) {
2716  size_t izeta = j;
2717  for (size_t m = 1; m < m_kk; m++) {
2718  if (charge(m) > 0.0) {
2719  size_t jzeta = m;
2720  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
2721  double zeta_L = m_Psi_ijk_L[n];
2722  if (zeta_L != 0.0) {
2723  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_L;
2724  }
2725  }
2726  }
2727  }
2728  }
2729  sum7 += molality[j]*molality[j]*molality[j]*m_Mu_nnn_L[j];
2730  }
2731  }
2732  double sum_m_phi_minus_1 = 2.0 *
2733  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
2734  // Calculate the osmotic coefficient from
2735  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
2736  double d_osmotic_coef_dT;
2737  if (molalitysum > 1.0E-150) {
2738  d_osmotic_coef_dT = 0.0 + (sum_m_phi_minus_1 / molalitysum);
2739  } else {
2740  d_osmotic_coef_dT = 0.0;
2741  }
2742 
2743  if (m_debugCalc) {
2744  writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
2745  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
2746  term1, sum1, sum2, sum3, sum4, sum5);
2747  writelogf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dT =%10.6f\n",
2748  sum_m_phi_minus_1, d_osmotic_coef_dT);
2749  writelog(" Step 10: \n");
2750  }
2751  double d_lnwateract_dT = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dT;
2752 
2753  // In Cantera, we define the activity coefficient of the solvent as
2754  //
2755  // act_0 = actcoeff_0 * Xmol_0
2756  //
2757  // We have just computed act_0. However, this routine returns
2758  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
2759  m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT;
2760  if (m_debugCalc) {
2761  double d_wateract_dT = exp(d_lnwateract_dT);
2762  writelogf(" d_ln_a_water_dT = %10.6f d_a_water_dT=%10.6f\n\n",
2763  d_lnwateract_dT, d_wateract_dT);
2764  }
2765 }
2766 
2767 void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const
2768 {
2769  static const int cacheId = m_cache.getId();
2770  CachedScalar cached = m_cache.getScalar(cacheId);
2771  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
2772  return;
2773  }
2774 
2775  // Zero the unscaled 2nd derivatives
2776  m_d2lnActCoeffMolaldT2_Unscaled.assign(m_kk, 0.0);
2777 
2778  //! Calculate the unscaled 2nd derivatives
2779  s_updatePitzer_d2lnMolalityActCoeff_dT2();
2780 
2781  for (size_t k = 1; k < m_kk; k++) {
2782  if (CROP_speciesCropped_[k] == 2) {
2783  m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0;
2784  }
2785  }
2786 
2787  if (CROP_speciesCropped_[0]) {
2788  m_d2lnActCoeffMolaldT2_Unscaled[0] = 0.0;
2789  }
2790 
2791  // Scale the 2nd derivatives
2792  s_updateScaling_pHScaling_dT2();
2793 }
2794 
2795 void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const
2796 {
2797  // HKM -> Assumption is made that the solvent is species 0.
2798  if (m_indexSolvent != 0) {
2799  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
2800  "Wrong index solvent value!");
2801  }
2802 
2803  const double* molality = m_molalitiesCropped.data();
2804 
2805  // Local variables defined by Coltrin
2806  double etheta[5][5], etheta_prime[5][5], sqrtIs;
2807 
2808  // Molality based ionic strength of the solution
2809  double Is = 0.0;
2810 
2811  // Molarcharge of the solution: In Pitzer's notation, this is his variable
2812  // called "Z".
2813  double molarcharge = 0.0;
2814 
2815  // molalitysum is the sum of the molalities over all solutes, even those
2816  // with zero charge.
2817  double molalitysum = 0.0;
2818 
2819  debuglog("\n Debugging information from s_Pitzer_d2lnMolalityActCoeff_dT2()\n",
2820  m_debugCalc);
2821 
2822  // Make sure the counter variables are setup
2823  counterIJ_setup();
2824 
2825  // ---------- Calculate common sums over solutes ---------------------
2826  for (size_t n = 1; n < m_kk; n++) {
2827  // ionic strength
2828  Is += charge(n) * charge(n) * molality[n];
2829  // total molar charge
2830  molarcharge += fabs(charge(n)) * molality[n];
2831  molalitysum += molality[n];
2832  }
2833  Is *= 0.5;
2834 
2835  // Store the ionic molality in the object for reference.
2836  m_IionicMolality = Is;
2837  sqrtIs = sqrt(Is);
2838  if (m_debugCalc) {
2839  writelog(" Step 1: \n");
2840  writelogf(" ionic strenth = %14.7le \n total molar "
2841  "charge = %14.7le \n", Is, molarcharge);
2842  }
2843 
2844  // The following call to calc_lambdas() calculates all 16 elements of the
2845  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
2846  calc_lambdas(Is);
2847 
2848  // Step 2: Find the coefficients E-theta and E-thetaprime for all
2849  // combinations of positive unlike charges up to 4
2850  debuglog(" Step 2: \n", m_debugCalc);
2851  for (int z1 = 1; z1 <=4; z1++) {
2852  for (int z2 =1; z2 <=4; z2++) {
2853  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
2854  if (m_debugCalc) {
2855  writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
2856  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
2857  }
2858  }
2859  }
2860 
2861  debuglog(" Step 3: \n"
2862  " Species Species g(x) hfunc(x) \n",
2863  m_debugCalc);
2864 
2865  // calculate gfunc(x) and hfunc(x) for each cation-anion pair MX. In the
2866  // original literature, hfunc, was called gprime. However, it's not the
2867  // derivative of gfunc(x), so I renamed it.
2868  for (size_t i = 1; i < (m_kk - 1); i++) {
2869  for (size_t j = (i+1); j < m_kk; j++) {
2870  // Find the counterIJ for the symmetric binary interaction
2871  size_t n = m_kk*i + j;
2872  size_t counterIJ = m_CounterIJ[n];
2873 
2874  // Only loop over oppositely charge species
2875  if (charge(i)*charge(j) < 0) {
2876  // x is a reduced function variable
2877  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
2878  if (x1 > 1.0E-100) {
2879  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 *x1);
2880  m_hfunc_IJ[counterIJ] = -2.0*
2881  (1.0-(1.0 + x1 + 0.5*x1 * x1) * exp(-x1)) / (x1 * x1);
2882  } else {
2883  m_gfunc_IJ[counterIJ] = 0.0;
2884  m_hfunc_IJ[counterIJ] = 0.0;
2885  }
2886 
2887  if (m_Beta2MX_ij_LL[counterIJ] != 0.0) {
2888  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
2889  if (x2 > 1.0E-100) {
2890  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
2891  m_h2func_IJ[counterIJ] = -2.0 *
2892  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
2893  } else {
2894  m_g2func_IJ[counterIJ] = 0.0;
2895  m_h2func_IJ[counterIJ] = 0.0;
2896  }
2897  }
2898  } else {
2899  m_gfunc_IJ[counterIJ] = 0.0;
2900  m_hfunc_IJ[counterIJ] = 0.0;
2901  }
2902  if (m_debugCalc) {
2903  writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i), speciesName(j),
2904  m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
2905  }
2906  }
2907  }
2908 
2909  // SUBSECTION TO CALCULATE BMX_L, BprimeMX_LL, BphiMX_L
2910  // These are now temperature derivatives of the previously calculated
2911  // quantities.
2912  debuglog(" Step 4: \n"
2913  " Species Species BMX BprimeMX BphiMX \n",
2914  m_debugCalc);
2915 
2916  for (size_t i = 1; i < m_kk - 1; i++) {
2917  for (size_t j = i+1; j < m_kk; j++) {
2918  // Find the counterIJ for the symmetric binary interaction
2919  size_t n = m_kk*i + j;
2920  size_t counterIJ = m_CounterIJ[n];
2921 
2922  // both species have a non-zero charge, and one is positive
2923  // and the other is negative
2924  if (charge(i)*charge(j) < 0.0) {
2925  m_BMX_IJ_LL[counterIJ] = m_Beta0MX_ij_LL[counterIJ]
2926  + m_Beta1MX_ij_LL[counterIJ] * m_gfunc_IJ[counterIJ]
2927  + m_Beta2MX_ij_LL[counterIJ] * m_g2func_IJ[counterIJ];
2928  if (m_debugCalc) {
2929  writelogf("%d %g: %g %g %g %g\n",
2930  counterIJ, m_BMX_IJ_LL[counterIJ], m_Beta0MX_ij_LL[counterIJ],
2931  m_Beta1MX_ij_LL[counterIJ], m_Beta2MX_ij_LL[counterIJ], m_gfunc_IJ[counterIJ]);
2932  }
2933  if (Is > 1.0E-150) {
2934  m_BprimeMX_IJ_LL[counterIJ] = (m_Beta1MX_ij_LL[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
2935  m_Beta2MX_ij_LL[counterIJ] * m_h2func_IJ[counterIJ]/Is);
2936  } else {
2937  m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2938  }
2939  m_BphiMX_IJ_LL[counterIJ] = m_BMX_IJ_LL[counterIJ] + Is*m_BprimeMX_IJ_LL[counterIJ];
2940  } else {
2941  m_BMX_IJ_LL[counterIJ] = 0.0;
2942  m_BprimeMX_IJ_LL[counterIJ] = 0.0;
2943  m_BphiMX_IJ_LL[counterIJ] = 0.0;
2944  }
2945  if (m_debugCalc) {
2946  writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
2947  speciesName(i), speciesName(j),
2948  m_BMX_IJ_LL[counterIJ], m_BprimeMX_IJ_LL[counterIJ], m_BphiMX_IJ_LL[counterIJ]);
2949  }
2950  }
2951  }
2952 
2953  // --------- SUBSECTION TO CALCULATE CMX_LL ----------
2954  debuglog(" Step 5: \n"
2955  " Species Species CMX \n", m_debugCalc);
2956  for (size_t i = 1; i < m_kk-1; i++) {
2957  for (size_t j = i+1; j < m_kk; j++) {
2958  // Find the counterIJ for the symmetric binary interaction
2959  size_t n = m_kk*i + j;
2960  size_t counterIJ = m_CounterIJ[n];
2961 
2962  // both species have a non-zero charge, and one is positive
2963  // and the other is negative
2964  if (charge(i)*charge(j) < 0.0) {
2965  m_CMX_IJ_LL[counterIJ] = m_CphiMX_ij_LL[counterIJ]/
2966  (2.0* sqrt(fabs(charge(i)*charge(j))));
2967  } else {
2968  m_CMX_IJ_LL[counterIJ] = 0.0;
2969  }
2970  if (m_debugCalc) {
2971  writelogf(" %-16s %-16s %11.7f \n",
2972  speciesName(i), speciesName(j), m_CMX_IJ_LL[counterIJ]);
2973  }
2974  }
2975  }
2976 
2977  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
2978  debuglog(" Step 6: \n"
2979  " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
2980  m_debugCalc);
2981  for (size_t i = 1; i < m_kk-1; i++) {
2982  for (size_t j = i+1; j < m_kk; j++) {
2983  // Find the counterIJ for the symmetric binary interaction
2984  size_t n = m_kk*i + j;
2985  size_t counterIJ = m_CounterIJ[n];
2986 
2987  // both species have a non-zero charge, and one is positive
2988  // and the other is negative
2989  if (charge(i)*charge(j) > 0) {
2990  m_Phi_IJ_LL[counterIJ] = m_Theta_ij_LL[counterIJ];
2991  m_Phiprime_IJ[counterIJ] = 0.0;
2992  m_PhiPhi_IJ_LL[counterIJ] = m_Phi_IJ_LL[counterIJ];
2993  } else {
2994  m_Phi_IJ_LL[counterIJ] = 0.0;
2995  m_Phiprime_IJ[counterIJ] = 0.0;
2996  m_PhiPhi_IJ_LL[counterIJ] = 0.0;
2997  }
2998  if (m_debugCalc) {
2999  writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
3000  speciesName(i), speciesName(j),
3001  m_Phi_IJ_LL[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_LL[counterIJ]);
3002  }
3003  }
3004  }
3005 
3006  // ----------- SUBSECTION FOR CALCULATION OF d2FdT2 ---------------------
3007  debuglog(" Step 7: \n", m_debugCalc);
3008  double d2AphidT2 = d2A_DebyedT2_TP() / 3.0;
3009  double d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3010  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3011  if (m_debugCalc) {
3012  writelogf(" initial value of d2FdT2 = %10.6f \n", d2FdT2);
3013  }
3014  for (size_t i = 1; i < m_kk-1; i++) {
3015  for (size_t j = i+1; j < m_kk; j++) {
3016  // Find the counterIJ for the symmetric binary interaction
3017  size_t n = m_kk*i + j;
3018  size_t counterIJ = m_CounterIJ[n];
3019 
3020  // both species have a non-zero charge, and one is positive
3021  // and the other is negative
3022  if (charge(i)*charge(j) < 0) {
3023  d2FdT2 += molality[i]*molality[j] * m_BprimeMX_IJ_LL[counterIJ];
3024  }
3025 
3026  // Both species have a non-zero charge, and they
3027  // have the same sign, e.g., both positive or both negative.
3028  if (charge(i)*charge(j) > 0) {
3029  d2FdT2 += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3030  }
3031  if (m_debugCalc) {
3032  writelogf(" d2FdT2 = %10.6f \n", d2FdT2);
3033  }
3034  }
3035  }
3036  debuglog(" Step 8: \n", m_debugCalc);
3037 
3038  for (size_t i = 1; i < m_kk; i++) {
3039  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdT FOR CATIONS -----
3040  if (charge(i) > 0) {
3041  // species i is the cation (positive) to calc the actcoeff
3042  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3043  double sum1 = 0.0;
3044  double sum2 = 0.0;
3045  double sum3 = 0.0;
3046  double sum4 = 0.0;
3047  double sum5 = 0.0;
3048  for (size_t j = 1; j < m_kk; j++) {
3049  // Find the counterIJ for the symmetric binary interaction
3050  size_t n = m_kk*i + j;
3051  size_t counterIJ = m_CounterIJ[n];
3052 
3053  if (charge(j) < 0.0) {
3054  // sum over all anions
3055  sum1 += molality[j]*
3056  (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3057  if (j < m_kk-1) {
3058  // This term is the ternary interaction involving the
3059  // non-duplicate sum over double anions, j, k, with
3060  // respect to the cation, i.
3061  for (size_t k = j+1; k < m_kk; k++) {
3062  // an inner sum over all anions
3063  if (charge(k) < 0.0) {
3064  n = k + j * m_kk + i * m_kk * m_kk;
3065  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3066  }
3067  }
3068  }
3069  }
3070 
3071  if (charge(j) > 0.0) {
3072  // sum over all cations
3073  if (j != i) {
3074  sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3075  }
3076  for (size_t k = 1; k < m_kk; k++) {
3077  if (charge(k) < 0.0) {
3078  // two inner sums over anions
3079  n = k + j * m_kk + i * m_kk * m_kk;
3080  sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3081 
3082  // Find the counterIJ for the j,k interaction
3083  n = m_kk*j + k;
3084  size_t counterIJ2 = m_CounterIJ[n];
3085  sum4 += fabs(charge(i)) *
3086  molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3087  }
3088  }
3089  }
3090 
3091  // Handle neutral j species
3092  if (charge(j) == 0) {
3093  sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3094  // Zeta interaction term
3095  for (size_t k = 1; k < m_kk; k++) {
3096  if (charge(k) < 0.0) {
3097  size_t izeta = j;
3098  size_t jzeta = i;
3099  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3100  double zeta_LL = m_Psi_ijk_LL[n];
3101  if (zeta_LL != 0.0) {
3102  sum5 += molality[j]*molality[k]*zeta_LL;
3103  }
3104  }
3105  }
3106  }
3107  }
3108  // Add all of the contributions up to yield the log of the
3109  // solute activity coefficients (molality scale)
3110  m_d2lnActCoeffMolaldT2_Unscaled[i] =
3111  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3112  if (m_debugCalc) {
3113  writelogf(" %-16s d2lngammadT2[i]=%10.6f \n",
3114  speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3115  writelogf(" %12g %12g %12g %12g %12g %12g\n",
3116  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
3117  }
3118  }
3119 
3120  // ------ SUBSECTION FOR CALCULATING THE d2ACTCOEFFdT2 FOR ANIONS ------
3121  if (charge(i) < 0) {
3122  // species i is an anion (negative)
3123  double zsqd2FdT2 = charge(i)*charge(i)*d2FdT2;
3124  double sum1 = 0.0;
3125  double sum2 = 0.0;
3126  double sum3 = 0.0;
3127  double sum4 = 0.0;
3128  double sum5 = 0.0;
3129  for (size_t j = 1; j < m_kk; j++) {
3130  // Find the counterIJ for the symmetric binary interaction
3131  size_t n = m_kk*i + j;
3132  size_t counterIJ = m_CounterIJ[n];
3133 
3134  // For Anions, do the cation interactions.
3135  if (charge(j) > 0) {
3136  sum1 += molality[j]*
3137  (2.0*m_BMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3138  if (j < m_kk-1) {
3139  for (size_t k = j+1; k < m_kk; k++) {
3140  // an inner sum over all cations
3141  if (charge(k) > 0) {
3142  n = k + j * m_kk + i * m_kk * m_kk;
3143  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3144  }
3145  }
3146  }
3147  }
3148 
3149  // For Anions, do the other anion interactions.
3150  if (charge(j) < 0.0) {
3151  // sum over all anions
3152  if (j != i) {
3153  sum2 += molality[j]*(2.0*m_Phi_IJ_LL[counterIJ]);
3154  }
3155  for (size_t k = 1; k < m_kk; k++) {
3156  if (charge(k) > 0.0) {
3157  // two inner sums over cations
3158  n = k + j * m_kk + i * m_kk * m_kk;
3159  sum2 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3160  // Find the counterIJ for the symmetric binary interaction
3161  n = m_kk*j + k;
3162  size_t counterIJ2 = m_CounterIJ[n];
3163  sum4 += fabs(charge(i)) *
3164  molality[j]*molality[k]*m_CMX_IJ_LL[counterIJ2];
3165  }
3166  }
3167  }
3168 
3169  // for Anions, do the neutral species interaction
3170  if (charge(j) == 0.0) {
3171  sum5 += molality[j]*2.0*m_Lambda_nj_LL(j,i);
3172  // Zeta interaction term
3173  for (size_t k = 1; k < m_kk; k++) {
3174  if (charge(k) > 0.0) {
3175  size_t izeta = j;
3176  size_t jzeta = k;
3177  size_t kzeta = i;
3178  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3179  double zeta_LL = m_Psi_ijk_LL[n];
3180  if (zeta_LL != 0.0) {
3181  sum5 += molality[j]*molality[k]*zeta_LL;
3182  }
3183  }
3184  }
3185  }
3186  }
3187  m_d2lnActCoeffMolaldT2_Unscaled[i] =
3188  zsqd2FdT2 + sum1 + sum2 + sum3 + sum4 + sum5;
3189  if (m_debugCalc) {
3190  writelogf(" %-16s d2lngammadT2[i]=%10.6f\n",
3191  speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3192  writelogf(" %12g %12g %12g %12g %12g %12g\n",
3193  zsqd2FdT2, sum1, sum2, sum3, sum4, sum5);
3194  }
3195  }
3196 
3197  // SUBSECTION FOR CALCULATING NEUTRAL SOLUTE ACT COEFF
3198  // equations agree with my notes,
3199  // Equations agree with Pitzer,
3200  if (charge(i) == 0.0) {
3201  double sum1 = 0.0;
3202  double sum3 = 0.0;
3203  for (size_t j = 1; j < m_kk; j++) {
3204  sum1 += molality[j]*2.0*m_Lambda_nj_LL(i,j);
3205  // Zeta term -> we piggyback on the psi term
3206  if (charge(j) > 0.0) {
3207  for (size_t k = 1; k < m_kk; k++) {
3208  if (charge(k) < 0.0) {
3209  size_t n = k + j * m_kk + i * m_kk * m_kk;
3210  sum3 += molality[j]*molality[k]*m_Psi_ijk_LL[n];
3211  }
3212  }
3213  }
3214  }
3215  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_LL[i];
3216  m_d2lnActCoeffMolaldT2_Unscaled[i] = sum1 + sum2 + sum3;
3217  if (m_debugCalc) {
3218  writelog(" %-16s d2lngammadT2[i]=%10.6f \n",
3219  speciesName(i), m_d2lnActCoeffMolaldT2_Unscaled[i]);
3220  }
3221  }
3222 
3223  }
3224  debuglog(" Step 9: \n", m_debugCalc);
3225 
3226  // ------ SUBSECTION FOR CALCULATING THE d2 OSMOTIC COEFF dT2 ---------
3227  double sum1 = 0.0;
3228  double sum2 = 0.0;
3229  double sum3 = 0.0;
3230  double sum4 = 0.0;
3231  double sum5 = 0.0;
3232  double sum6 = 0.0;
3233  double sum7 = 0.0;
3234 
3235  // term1 is the temperature derivative of the DH term in the osmotic
3236  // coefficient expression
3237  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3238  // Is = Ionic strength on the molality scale (units of (gmol/kg))
3239  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3240  double term1 = -d2AphidT2 * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3241 
3242  for (size_t j = 1; j < m_kk; j++) {
3243  // Loop Over Cations
3244  if (charge(j) > 0.0) {
3245  for (size_t k = 1; k < m_kk; k++) {
3246  if (charge(k) < 0.0) {
3247  // Find the counterIJ for the symmetric j,k binary interaction
3248  size_t n = m_kk*j + k;
3249  size_t counterIJ = m_CounterIJ[n];
3250 
3251  sum1 += molality[j]*molality[k] *
3252  (m_BphiMX_IJ_LL[counterIJ] + molarcharge*m_CMX_IJ_LL[counterIJ]);
3253  }
3254  }
3255 
3256  for (size_t k = j+1; k < m_kk; k++) {
3257  if (j == (m_kk-1)) {
3258  // we should never reach this step
3259  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3260  "logic error 1 in Step 9 of hmw_act");
3261  }
3262  if (charge(k) > 0.0) {
3263  // Find the counterIJ for the symmetric j,k binary interaction
3264  // between 2 cations.
3265  size_t n = m_kk*j + k;
3266  size_t counterIJ = m_CounterIJ[n];
3267  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3268  for (size_t m = 1; m < m_kk; m++) {
3269  if (charge(m) < 0.0) {
3270  // species m is an anion
3271  n = m + k * m_kk + j * m_kk * m_kk;
3272  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3273  }
3274  }
3275  }
3276  }
3277  }
3278 
3279  // Loop Over Anions
3280  if (charge(j) < 0) {
3281  for (size_t k = j+1; k < m_kk; k++) {
3282  if (j == m_kk-1) {
3283  // we should never reach this step
3284  throw CanteraError("HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2",
3285  "logic error 2 in Step 9 of hmw_act");
3286  }
3287  if (charge(k) < 0) {
3288  // Find the counterIJ for the symmetric j,k binary interaction
3289  // between two anions
3290  size_t n = m_kk*j + k;
3291  size_t counterIJ = m_CounterIJ[n];
3292 
3293  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_LL[counterIJ];
3294  for (size_t m = 1; m < m_kk; m++) {
3295  if (charge(m) > 0.0) {
3296  n = m + k * m_kk + j * m_kk * m_kk;
3297  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_LL[n];
3298  }
3299  }
3300  }
3301  }
3302  }
3303 
3304  // Loop Over Neutral Species
3305  if (charge(j) == 0) {
3306  for (size_t k = 1; k < m_kk; k++) {
3307  if (charge(k) < 0.0) {
3308  sum4 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3309  }
3310  if (charge(k) > 0.0) {
3311  sum5 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3312  }
3313  if (charge(k) == 0.0) {
3314  if (k > j) {
3315  sum6 += molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3316  } else if (k == j) {
3317  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_LL(j,k);
3318  }
3319  }
3320  if (charge(k) < 0.0) {
3321  size_t izeta = j;
3322  for (size_t m = 1; m < m_kk; m++) {
3323  if (charge(m) > 0.0) {
3324  size_t jzeta = m;
3325  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3326  double zeta_LL = m_Psi_ijk_LL[n];
3327  if (zeta_LL != 0.0) {
3328  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_LL;
3329  }
3330  }
3331  }
3332  }
3333  }
3334 
3335  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_LL[j];
3336  }
3337  }
3338  double sum_m_phi_minus_1 = 2.0 *
3339  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3340  // Calculate the osmotic coefficient from
3341  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3342  double d2_osmotic_coef_dT2;
3343  if (molalitysum > 1.0E-150) {
3344  d2_osmotic_coef_dT2 = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3345  } else {
3346  d2_osmotic_coef_dT2 = 0.0;
3347  }
3348  if (m_debugCalc) {
3349  writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3350  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3351  term1, sum1, sum2, sum3, sum4, sum5);
3352  writelogf(" sum_m_phi_minus_1=%10.6f d2_osmotic_coef_dT2=%10.6f\n",
3353  sum_m_phi_minus_1, d2_osmotic_coef_dT2);
3354  writelog(" Step 10: \n");
3355  }
3356  double d2_lnwateract_dT2 = -(m_weightSolvent/1000.0) * molalitysum * d2_osmotic_coef_dT2;
3357 
3358  // In Cantera, we define the activity coefficient of the solvent as
3359  //
3360  // act_0 = actcoeff_0 * Xmol_0
3361  //
3362  // We have just computed act_0. However, this routine returns
3363  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3364  m_d2lnActCoeffMolaldT2_Unscaled[0] = d2_lnwateract_dT2;
3365 
3366  if (m_debugCalc) {
3367  double d2_wateract_dT2 = exp(d2_lnwateract_dT2);
3368  writelogf(" d2_ln_a_water_dT2 = %10.6f d2_a_water_dT2=%10.6f\n\n",
3369  d2_lnwateract_dT2, d2_wateract_dT2);
3370  }
3371 }
3372 
3373 void HMWSoln::s_update_dlnMolalityActCoeff_dP() const
3374 {
3375  static const int cacheId = m_cache.getId();
3376  CachedScalar cached = m_cache.getScalar(cacheId);
3377  if( cached.validate(temperature(), pressure(), stateMFNumber()) ) {
3378  return;
3379  }
3380 
3381  m_dlnActCoeffMolaldP_Unscaled.assign(m_kk, 0.0);
3382  s_updatePitzer_dlnMolalityActCoeff_dP();
3383 
3384  for (size_t k = 1; k < m_kk; k++) {
3385  if (CROP_speciesCropped_[k] == 2) {
3386  m_dlnActCoeffMolaldP_Unscaled[k] = 0.0;
3387  }
3388  }
3389 
3390  if (CROP_speciesCropped_[0]) {
3391  m_dlnActCoeffMolaldP_Unscaled[0] = 0.0;
3392  }
3393 
3394  s_updateScaling_pHScaling_dP();
3395 }
3396 
3397 void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const
3398 {
3399  // HKM -> Assumption is made that the solvent is species 0.
3400  if (m_indexSolvent != 0) {
3401  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3402  "Wrong index solvent value!");
3403  }
3404 
3405  const double* molality = m_molalitiesCropped.data();
3406 
3407  // Local variables defined by Coltrin
3408  double etheta[5][5], etheta_prime[5][5], sqrtIs;
3409 
3410  // Molality based ionic strength of the solution
3411  double Is = 0.0;
3412 
3413  // Molarcharge of the solution: In Pitzer's notation, this is his variable
3414  // called "Z".
3415  double molarcharge = 0.0;
3416 
3417  // molalitysum is the sum of the molalities over all solutes, even those
3418  // with zero charge.
3419  double molalitysum = 0.0;
3420  double currTemp = temperature();
3421  double currPres = pressure();
3422 
3423  debuglog("\n Debugging information from s_Pitzer_dlnMolalityActCoeff_dP()\n",
3424  m_debugCalc);
3425 
3426  // Make sure the counter variables are setup
3427  counterIJ_setup();
3428 
3429  // ---------- Calculate common sums over solutes ---------------------
3430  for (size_t n = 1; n < m_kk; n++) {
3431  // ionic strength
3432  Is += charge(n) * charge(n) * molality[n];
3433  // total molar charge
3434  molarcharge += fabs(charge(n)) * molality[n];
3435  molalitysum += molality[n];
3436  }
3437  Is *= 0.5;
3438 
3439  // Store the ionic molality in the object for reference.
3440  m_IionicMolality = Is;
3441  sqrtIs = sqrt(Is);
3442  if (m_debugCalc) {
3443  writelog(" Step 1: \n");
3444  writelogf(" ionic strenth = %14.7le \n total molar "
3445  "charge = %14.7le \n", Is, molarcharge);
3446  }
3447 
3448  // The following call to calc_lambdas() calculates all 16 elements of the
3449  // elambda and elambda1 arrays, given the value of the ionic strength (Is)
3450  calc_lambdas(Is);
3451 
3452 
3453  // Step 2: Find the coefficients E-theta and E-thetaprime for all
3454  // combinations of positive unlike charges up to 4
3455  debuglog(" Step 2: \n", m_debugCalc);
3456  for (int z1 = 1; z1 <=4; z1++) {
3457  for (int z2 =1; z2 <=4; z2++) {
3458  calc_thetas(z1, z2, &etheta[z1][z2], &etheta_prime[z1][z2]);
3459  if (m_debugCalc) {
3460  writelogf(" z1=%3d z2=%3d E-theta(I) = %f, E-thetaprime(I) = %f\n",
3461  z1, z2, etheta[z1][z2], etheta_prime[z1][z2]);
3462  }
3463  }
3464  }
3465 
3466  debuglog(" Step 3: \n"
3467  " Species Species g(x) hfunc(x)\n",
3468  m_debugCalc);
3469 
3470  // calculate g(x) and hfunc(x) for each cation-anion pair MX
3471  // In the original literature, hfunc, was called gprime. However,
3472  // it's not the derivative of g(x), so I renamed it.
3473  for (size_t i = 1; i < (m_kk - 1); i++) {
3474  for (size_t j = (i+1); j < m_kk; j++) {
3475  // Find the counterIJ for the symmetric binary interaction
3476  size_t n = m_kk*i + j;
3477  size_t counterIJ = m_CounterIJ[n];
3478 
3479  // Only loop over oppositely charge species
3480  if (charge(i)*charge(j) < 0) {
3481  // x is a reduced function variable
3482  double x1 = sqrtIs * m_Alpha1MX_ij[counterIJ];
3483  if (x1 > 1.0E-100) {
3484  m_gfunc_IJ[counterIJ] = 2.0*(1.0-(1.0 + x1) * exp(-x1)) / (x1 * x1);
3485  m_hfunc_IJ[counterIJ] = -2.0*
3486  (1.0-(1.0 + x1 + 0.5 * x1 * x1) * exp(-x1)) / (x1 * x1);
3487  } else {
3488  m_gfunc_IJ[counterIJ] = 0.0;
3489  m_hfunc_IJ[counterIJ] = 0.0;
3490  }
3491 
3492  if (m_Beta2MX_ij_P[counterIJ] != 0.0) {
3493  double x2 = sqrtIs * m_Alpha2MX_ij[counterIJ];
3494  if (x2 > 1.0E-100) {
3495  m_g2func_IJ[counterIJ] = 2.0*(1.0-(1.0 + x2) * exp(-x2)) / (x2 * x2);
3496  m_h2func_IJ[counterIJ] = -2.0 *
3497  (1.0-(1.0 + x2 + 0.5 * x2 * x2) * exp(-x2)) / (x2 * x2);
3498  } else {
3499  m_g2func_IJ[counterIJ] = 0.0;
3500  m_h2func_IJ[counterIJ] = 0.0;
3501  }
3502  }
3503  } else {
3504  m_gfunc_IJ[counterIJ] = 0.0;
3505  m_hfunc_IJ[counterIJ] = 0.0;
3506  }
3507  if (m_debugCalc) {
3508  writelogf(" %-16s %-16s %9.5f %9.5f \n", speciesName(i),
3509  speciesName(j), m_gfunc_IJ[counterIJ], m_hfunc_IJ[counterIJ]);
3510  }
3511  }
3512  }
3513 
3514  // SUBSECTION TO CALCULATE BMX_P, BprimeMX_P, BphiMX_P
3515  // These are now temperature derivatives of the previously calculated
3516  // quantities.
3517  debuglog(" Step 4: \n"
3518  " Species Species BMX BprimeMX BphiMX \n",
3519  m_debugCalc);
3520 
3521  for (size_t i = 1; i < m_kk - 1; i++) {
3522  for (size_t j = i+1; j < m_kk; j++) {
3523  // Find the counterIJ for the symmetric binary interaction
3524  size_t n = m_kk*i + j;
3525  size_t counterIJ = m_CounterIJ[n];
3526 
3527  // both species have a non-zero charge, and one is positive
3528  // and the other is negative
3529  if (charge(i)*charge(j) < 0.0) {
3530  m_BMX_IJ_P[counterIJ] = m_Beta0MX_ij_P[counterIJ]
3531  + m_Beta1MX_ij_P[counterIJ] * m_gfunc_IJ[counterIJ]
3532  + m_Beta2MX_ij_P[counterIJ] * m_g2func_IJ[counterIJ];
3533  if (m_debugCalc) {
3534  writelogf("%d %g: %g %g %g %g\n",
3535  counterIJ, m_BMX_IJ_P[counterIJ], m_Beta0MX_ij_P[counterIJ],
3536  m_Beta1MX_ij_P[counterIJ], m_Beta2MX_ij_P[counterIJ], m_gfunc_IJ[counterIJ]);
3537  }
3538  if (Is > 1.0E-150) {
3539  m_BprimeMX_IJ_P[counterIJ] = (m_Beta1MX_ij_P[counterIJ] * m_hfunc_IJ[counterIJ]/Is +
3540  m_Beta2MX_ij_P[counterIJ] * m_h2func_IJ[counterIJ]/Is);
3541  } else {
3542  m_BprimeMX_IJ_P[counterIJ] = 0.0;
3543  }
3544  m_BphiMX_IJ_P[counterIJ] = m_BMX_IJ_P[counterIJ] + Is*m_BprimeMX_IJ_P[counterIJ];
3545  } else {
3546  m_BMX_IJ_P[counterIJ] = 0.0;
3547  m_BprimeMX_IJ_P[counterIJ] = 0.0;
3548  m_BphiMX_IJ_P[counterIJ] = 0.0;
3549  }
3550  if (m_debugCalc) {
3551  writelogf(" %-16s %-16s %11.7f %11.7f %11.7f \n",
3552  speciesName(i), speciesName(j),
3553  m_BMX_IJ_P[counterIJ], m_BprimeMX_IJ_P[counterIJ], m_BphiMX_IJ_P[counterIJ]);
3554  }
3555  }
3556  }
3557 
3558  // --------- SUBSECTION TO CALCULATE CMX_P ----------
3559  debuglog(" Step 5: \n"
3560  " Species Species CMX \n", m_debugCalc);
3561  for (size_t i = 1; i < m_kk-1; i++) {
3562  for (size_t j = i+1; j < m_kk; j++) {
3563  // Find the counterIJ for the symmetric binary interaction
3564  size_t n = m_kk*i + j;
3565  size_t counterIJ = m_CounterIJ[n];
3566 
3567  // both species have a non-zero charge, and one is positive
3568  // and the other is negative
3569  if (charge(i)*charge(j) < 0.0) {
3570  m_CMX_IJ_P[counterIJ] = m_CphiMX_ij_P[counterIJ]/
3571  (2.0* sqrt(fabs(charge(i)*charge(j))));
3572  } else {
3573  m_CMX_IJ_P[counterIJ] = 0.0;
3574  }
3575  if (m_debugCalc) {
3576  writelogf(" %-16s %-16s %11.7f \n",
3577  speciesName(i), speciesName(j), m_CMX_IJ_P[counterIJ]);
3578  }
3579  }
3580  }
3581 
3582  // ------- SUBSECTION TO CALCULATE Phi, PhiPrime, and PhiPhi ----------
3583  debuglog(" Step 6: \n"
3584  " Species Species Phi_ij Phiprime_ij Phi^phi_ij \n",
3585  m_debugCalc);
3586  for (size_t i = 1; i < m_kk-1; i++) {
3587  for (size_t j = i+1; j < m_kk; j++) {
3588  // Find the counterIJ for the symmetric binary interaction
3589  size_t n = m_kk*i + j;
3590  size_t counterIJ = m_CounterIJ[n];
3591 
3592  // both species have a non-zero charge, and one is positive
3593  // and the other is negative
3594  if (charge(i)*charge(j) > 0) {
3595  m_Phi_IJ_P[counterIJ] = m_Theta_ij_P[counterIJ];
3596  m_Phiprime_IJ[counterIJ] = 0.0;
3597  m_PhiPhi_IJ_P[counterIJ] = m_Phi_IJ_P[counterIJ] + Is * m_Phiprime_IJ[counterIJ];
3598  } else {
3599  m_Phi_IJ_P[counterIJ] = 0.0;
3600  m_Phiprime_IJ[counterIJ] = 0.0;
3601  m_PhiPhi_IJ_P[counterIJ] = 0.0;
3602  }
3603  if (m_debugCalc) {
3604  writelogf(" %-16s %-16s %10.6f %10.6f %10.6f \n",
3605  speciesName(i), speciesName(j),
3606  m_Phi_IJ_P[counterIJ], m_Phiprime_IJ[counterIJ], m_PhiPhi_IJ_P[counterIJ]);
3607  }
3608  }
3609  }
3610 
3611  // ----------- SUBSECTION FOR CALCULATION OF dFdT ---------------------
3612  debuglog(" Step 7: \n", m_debugCalc);
3613  double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres);
3614  double dAphidP = dA_DebyedP /3.0;
3615  double dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is))
3616  + (2.0/1.2) * log(1.0+1.2*(sqrtIs)));
3617  if (m_debugCalc) {
3618  writelogf(" initial value of dFdP = %10.6f \n", dFdP);
3619  }
3620  for (size_t i = 1; i < m_kk-1; i++) {
3621  for (size_t j = i+1; j < m_kk; j++) {
3622  // Find the counterIJ for the symmetric binary interaction
3623  size_t n = m_kk*i + j;
3624  size_t counterIJ = m_CounterIJ[n];
3625 
3626  // both species have a non-zero charge, and one is positive
3627  // and the other is negative
3628  if (charge(i)*charge(j) < 0) {
3629  dFdP += molality[i]*molality[j] * m_BprimeMX_IJ_P[counterIJ];
3630  }
3631 
3632  // Both species have a non-zero charge, and they
3633  // have the same sign, e.g., both positive or both negative.
3634  if (charge(i)*charge(j) > 0) {
3635  dFdP += molality[i]*molality[j] * m_Phiprime_IJ[counterIJ];
3636  }
3637  if (m_debugCalc) {
3638  writelogf(" dFdP = %10.6f \n", dFdP);
3639  }
3640  }
3641  }
3642  debuglog(" Step 8: \n", m_debugCalc);
3643 
3644  for (size_t i = 1; i < m_kk; i++) {
3645  // -------- SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR CATIONS -----
3646  if (charge(i) > 0) {
3647  // species i is the cation (positive) to calc the actcoeff
3648  double zsqdFdP = charge(i)*charge(i)*dFdP;
3649  double sum1 = 0.0;
3650  double sum2 = 0.0;
3651  double sum3 = 0.0;
3652  double sum4 = 0.0;
3653  double sum5 = 0.0;
3654  for (size_t j = 1; j < m_kk; j++) {
3655  // Find the counterIJ for the symmetric binary interaction
3656  size_t n = m_kk*i + j;
3657  size_t counterIJ = m_CounterIJ[n];
3658 
3659  if (charge(j) < 0.0) {
3660  // sum over all anions
3661  sum1 += molality[j]*
3662  (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3663  if (j < m_kk-1) {
3664  // This term is the ternary interaction involving the
3665  // non-duplicate sum over double anions, j, k, with
3666  // respect to the cation, i.
3667  for (size_t k = j+1; k < m_kk; k++) {
3668  // an inner sum over all anions
3669  if (charge(k) < 0.0) {
3670  n = k + j * m_kk + i * m_kk * m_kk;
3671  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3672  }
3673  }
3674  }
3675  }
3676 
3677  if (charge(j) > 0.0) {
3678  // sum over all cations
3679  if (j != i) {
3680  sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3681  }
3682  for (size_t k = 1; k < m_kk; k++) {
3683  if (charge(k) < 0.0) {
3684  // two inner sums over anions
3685  n = k + j * m_kk + i * m_kk * m_kk;
3686  sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3687 
3688  // Find the counterIJ for the j,k interaction
3689  n = m_kk*j + k;
3690  size_t counterIJ2 = m_CounterIJ[n];
3691  sum4 += fabs(charge(i)) *
3692  molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3693  }
3694  }
3695  }
3696 
3697  // for Anions, do the neutral species interaction
3698  if (charge(j) == 0) {
3699  sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3700  // Zeta interaction term
3701  for (size_t k = 1; k < m_kk; k++) {
3702  if (charge(k) < 0.0) {
3703  size_t izeta = j;
3704  size_t jzeta = i;
3705  n = izeta * m_kk * m_kk + jzeta * m_kk + k;
3706  double zeta_P = m_Psi_ijk_P[n];
3707  if (zeta_P != 0.0) {
3708  sum5 += molality[j]*molality[k]*zeta_P;
3709  }
3710  }
3711  }
3712  }
3713  }
3714 
3715  // Add all of the contributions up to yield the log of the
3716  // solute activity coefficients (molality scale)
3717  m_dlnActCoeffMolaldP_Unscaled[i] =
3718  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3719 
3720  if (m_debugCalc) {
3721  writelogf(" %-16s lngamma[i]=%10.6f \n",
3722  speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
3723  writelogf(" %12g %12g %12g %12g %12g %12g\n",
3724  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
3725  }
3726  }
3727 
3728  // ------ SUBSECTION FOR CALCULATING THE dACTCOEFFdP FOR ANIONS ------
3729  if (charge(i) < 0) {
3730  // species i is an anion (negative)
3731  double zsqdFdP = charge(i)*charge(i)*dFdP;
3732  double sum1 = 0.0;
3733  double sum2 = 0.0;
3734  double sum3 = 0.0;
3735  double sum4 = 0.0;
3736  double sum5 = 0.0;
3737  for (size_t j = 1; j < m_kk; j++) {
3738  // Find the counterIJ for the symmetric binary interaction
3739  size_t n = m_kk*i + j;
3740  size_t counterIJ = m_CounterIJ[n];
3741 
3742  // For Anions, do the cation interactions.
3743  if (charge(j) > 0) {
3744  sum1 += molality[j] *
3745  (2.0*m_BMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3746  if (j < m_kk-1) {
3747  for (size_t k = j+1; k < m_kk; k++) {
3748  // an inner sum over all cations
3749  if (charge(k) > 0) {
3750  n = k + j * m_kk + i * m_kk * m_kk;
3751  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3752  }
3753  }
3754  }
3755  }
3756 
3757  // For Anions, do the other anion interactions.
3758  if (charge(j) < 0.0) {
3759  // sum over all anions
3760  if (j != i) {
3761  sum2 += molality[j]*(2.0*m_Phi_IJ_P[counterIJ]);
3762  }
3763  for (size_t k = 1; k < m_kk; k++) {
3764  if (charge(k) > 0.0) {
3765  // two inner sums over cations
3766  n = k + j * m_kk + i * m_kk * m_kk;
3767  sum2 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3768  // Find the counterIJ for the symmetric binary interaction
3769  n = m_kk*j + k;
3770  size_t counterIJ2 = m_CounterIJ[n];
3771  sum4 += fabs(charge(i))*
3772  molality[j]*molality[k]*m_CMX_IJ_P[counterIJ2];
3773  }
3774  }
3775  }
3776 
3777  // for Anions, do the neutral species interaction
3778  if (charge(j) == 0.0) {
3779  sum5 += molality[j]*2.0*m_Lambda_nj_P(j,i);
3780  // Zeta interaction term
3781  for (size_t k = 1; k < m_kk; k++) {
3782  if (charge(k) > 0.0) {
3783  size_t izeta = j;
3784  size_t jzeta = k;
3785  size_t kzeta = i;
3786  n = izeta * m_kk * m_kk + jzeta * m_kk + kzeta;
3787  double zeta_P = m_Psi_ijk_P[n];
3788  if (zeta_P != 0.0) {
3789  sum5 += molality[j]*molality[k]*zeta_P;
3790  }
3791  }
3792  }
3793  }
3794  }
3795  m_dlnActCoeffMolaldP_Unscaled[i] =
3796  zsqdFdP + sum1 + sum2 + sum3 + sum4 + sum5;
3797  if (m_debugCalc) {
3798  writelogf(" %-16s lndactcoeffmolaldP[i]=%10.6f \n",
3799  speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
3800  writelogf(" %12g %12g %12g %12g %12g %12g\n",
3801  zsqdFdP, sum1, sum2, sum3, sum4, sum5);
3802  }
3803  }
3804 
3805  // ------ SUBSECTION FOR CALCULATING d NEUTRAL SOLUTE ACT COEFF dP -----
3806  if (charge(i) == 0.0) {
3807  double sum1 = 0.0;
3808  double sum3 = 0.0;
3809  for (size_t j = 1; j < m_kk; j++) {
3810  sum1 += molality[j]*2.0*m_Lambda_nj_P(i,j);
3811  // Zeta term -> we piggyback on the psi term
3812  if (charge(j) > 0.0) {
3813  for (size_t k = 1; k < m_kk; k++) {
3814  if (charge(k) < 0.0) {
3815  size_t n = k + j * m_kk + i * m_kk * m_kk;
3816  sum3 += molality[j]*molality[k]*m_Psi_ijk_P[n];
3817  }
3818  }
3819  }
3820  }
3821  double sum2 = 3.0 * molality[i] * molality[i] * m_Mu_nnn_P[i];
3822  m_dlnActCoeffMolaldP_Unscaled[i] = sum1 + sum2 + sum3;
3823  if (m_debugCalc) {
3824  writelogf(" %-16s dlnActCoeffMolaldP[i]=%10.6f \n",
3825  speciesName(i), m_dlnActCoeffMolaldP_Unscaled[i]);
3826  }
3827  }
3828  }
3829  debuglog(" Step 9: \n", m_debugCalc);
3830 
3831  // ------ SUBSECTION FOR CALCULATING THE d OSMOTIC COEFF dP ---------
3832  double sum1 = 0.0;
3833  double sum2 = 0.0;
3834  double sum3 = 0.0;
3835  double sum4 = 0.0;
3836  double sum5 = 0.0;
3837  double sum6 = 0.0;
3838  double sum7 = 0.0;
3839 
3840  // term1 is the temperature derivative of the DH term in the osmotic
3841  // coefficient expression
3842  // b = 1.2 sqrt(kg/gmol) <- arbitrarily set in all Pitzer implementations.
3843  // Is = Ionic strength on the molality scale (units of (gmol/kg))
3844  // Aphi = A_Debye / 3 (units of sqrt(kg/gmol))
3845  double term1 = -dAphidP * Is * sqrt(Is) / (1.0 + 1.2 * sqrt(Is));
3846 
3847  for (size_t j = 1; j < m_kk; j++) {
3848  // Loop Over Cations
3849  if (charge(j) > 0.0) {
3850  for (size_t k = 1; k < m_kk; k++) {
3851  if (charge(k) < 0.0) {
3852  // Find the counterIJ for the symmetric j,k binary interaction
3853  size_t n = m_kk*j + k;
3854  size_t counterIJ = m_CounterIJ[n];
3855  sum1 += molality[j]*molality[k]*
3856  (m_BphiMX_IJ_P[counterIJ] + molarcharge*m_CMX_IJ_P[counterIJ]);
3857  }
3858  }
3859 
3860  for (size_t k = j+1; k < m_kk; k++) {
3861  if (j == (m_kk-1)) {
3862  // we should never reach this step
3863  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3864  "logic error 1 in Step 9 of hmw_act");
3865  }
3866  if (charge(k) > 0.0) {
3867  // Find the counterIJ for the symmetric j,k binary interaction
3868  // between 2 cations.
3869  size_t n = m_kk*j + k;
3870  size_t counterIJ = m_CounterIJ[n];
3871  sum2 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3872  for (size_t m = 1; m < m_kk; m++) {
3873  if (charge(m) < 0.0) {
3874  // species m is an anion
3875  n = m + k * m_kk + j * m_kk * m_kk;
3876  sum2 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3877  }
3878  }
3879  }
3880  }
3881  }
3882 
3883  // Loop Over Anions
3884  if (charge(j) < 0) {
3885  for (size_t k = j+1; k < m_kk; k++) {
3886  if (j == m_kk-1) {
3887  // we should never reach this step
3888  throw CanteraError("HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP",
3889  "logic error 2 in Step 9 of hmw_act");
3890  }
3891  if (charge(k) < 0) {
3892  // Find the counterIJ for the symmetric j,k binary interaction
3893  // between two anions
3894  size_t n = m_kk*j + k;
3895  size_t counterIJ = m_CounterIJ[n];
3896 
3897  sum3 += molality[j]*molality[k]*m_PhiPhi_IJ_P[counterIJ];
3898  for (size_t m = 1; m < m_kk; m++) {
3899  if (charge(m) > 0.0) {
3900  n = m + k * m_kk + j * m_kk * m_kk;
3901  sum3 += molality[j]*molality[k]*molality[m]*m_Psi_ijk_P[n];
3902  }
3903  }
3904  }
3905  }
3906  }
3907 
3908  // Loop Over Neutral Species
3909  if (charge(j) == 0) {
3910  for (size_t k = 1; k < m_kk; k++) {
3911  if (charge(k) < 0.0) {
3912  sum4 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3913  }
3914  if (charge(k) > 0.0) {
3915  sum5 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3916  }
3917  if (charge(k) == 0.0) {
3918  if (k > j) {
3919  sum6 += molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3920  } else if (k == j) {
3921  sum6 += 0.5 * molality[j]*molality[k]*m_Lambda_nj_P(j,k);
3922  }
3923  }
3924  if (charge(k) < 0.0) {
3925  size_t izeta = j;
3926  for (size_t m = 1; m < m_kk; m++) {
3927  if (charge(m) > 0.0) {
3928  size_t jzeta = m;
3929  size_t n = k + jzeta * m_kk + izeta * m_kk * m_kk;
3930  double zeta_P = m_Psi_ijk_P[n];
3931  if (zeta_P != 0.0) {
3932  sum7 += molality[izeta]*molality[jzeta]*molality[k]*zeta_P;
3933  }
3934  }
3935  }
3936  }
3937  }
3938 
3939  sum7 += molality[j] * molality[j] * molality[j] * m_Mu_nnn_P[j];
3940  }
3941  }
3942  double sum_m_phi_minus_1 = 2.0 *
3943  (term1 + sum1 + sum2 + sum3 + sum4 + sum5 + sum6 + sum7);
3944 
3945  // Calculate the osmotic coefficient from
3946  // osmotic_coeff = 1 + dGex/d(M0noRT) / sum(molality_i)
3947  double d_osmotic_coef_dP;
3948  if (molalitysum > 1.0E-150) {
3949  d_osmotic_coef_dP = 0.0 + (sum_m_phi_minus_1 / molalitysum);
3950  } else {
3951  d_osmotic_coef_dP = 0.0;
3952  }
3953  if (m_debugCalc) {
3954  writelogf(" term1=%10.6f sum1=%10.6f sum2=%10.6f "
3955  "sum3=%10.6f sum4=%10.6f sum5=%10.6f\n",
3956  term1, sum1, sum2, sum3, sum4, sum5);
3957  writelogf(" sum_m_phi_minus_1=%10.6f d_osmotic_coef_dP =%10.6f\n",
3958  sum_m_phi_minus_1, d_osmotic_coef_dP);
3959  writelog(" Step 10: \n");
3960  }
3961  double d_lnwateract_dP = -(m_weightSolvent/1000.0) * molalitysum * d_osmotic_coef_dP;
3962 
3963  // In Cantera, we define the activity coefficient of the solvent as
3964  //
3965  // act_0 = actcoeff_0 * Xmol_0
3966  //
3967  // We have just computed act_0. However, this routine returns
3968  // ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0).
3969  m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP;
3970  if (m_debugCalc) {
3971  writelogf(" d_ln_a_water_dP = %10.6f d_a_water_dP=%10.6f\n\n",
3972  d_lnwateract_dP, exp(d_lnwateract_dP));
3973  }
3974 }
3975 
3976 void HMWSoln::calc_lambdas(double is) const
3977 {
3978  if( m_last_is == is ) {
3979  return;
3980  }
3981  m_last_is = is;
3982 
3983  // Coefficients c1-c4 are used to approximate the integral function "J";
3984  // aphi is the Debye-Huckel constant at 25 C
3985  double c1 = 4.581, c2 = 0.7237, c3 = 0.0120, c4 = 0.528;
3986  double aphi = 0.392; /* Value at 25 C */
3987  if (m_debugCalc) {
3988  writelogf(" Is = %g\n", is);
3989  }
3990  if (is < 1.0E-150) {
3991  for (int i = 0; i < 17; i++) {
3992  elambda[i] = 0.0;
3993  elambda1[i] = 0.0;
3994  }
3995  return;
3996  }
3997 
3998  // Calculate E-lambda terms for charge combinations of like sign,
3999  // using method of Pitzer (1975). Charges up to 4 are calculated.
4000  for (int i=1; i<=4; i++) {
4001  for (int j=i; j<=4; j++) {
4002  int ij = i*j;
4003 
4004  // calculate the product of the charges
4005  double zprod = (double)ij;
4006 
4007  // calculate Xmn (A1) from Harvie, Weare (1980).
4008  double x = 6.0* zprod * aphi * sqrt(is); // eqn 23
4009 
4010  double jfunc = x / (4.0 + c1*pow(x,-c2)*exp(-c3*pow(x,c4))); // eqn 47
4011 
4012  double t = c3 * c4 * pow(x,c4);
4013  double dj = c1* pow(x,(-c2-1.0)) * (c2+t) * exp(-c3*pow(x,c4));
4014  double jprime = (jfunc/x)*(1.0 + jfunc*dj);
4015 
4016  elambda[ij] = zprod*jfunc / (4.0*is); // eqn 14
4017  elambda1[ij] = (3.0*zprod*zprod*aphi*jprime/(4.0*sqrt(is))
4018  - elambda[ij])/is;
4019  if (m_debugCalc) {
4020  writelogf(" ij = %d, elambda = %g, elambda1 = %g\n",
4021  ij, elambda[ij], elambda1[ij]);
4022  }
4023  }
4024  }
4025 }
4026 
4027 void HMWSoln::calc_thetas(int z1, int z2,
4028  double* etheta, double* etheta_prime) const
4029 {
4030  // Calculate E-theta(i) and E-theta'(I) using method of Pitzer (1987)
4031  int i = abs(z1);
4032  int j = abs(z2);
4033 
4034  AssertThrowMsg(i <= 4 && j <= 4, "HMWSoln::calc_thetas",
4035  "we shouldn't be here");
4036  AssertThrowMsg(i != 0 && j != 0, "HMWSoln::calc_thetas",
4037  "called with one species being neutral");
4038 
4039  // Check to see if the charges are of opposite sign. If they are of opposite
4040  // sign then their etheta interaction is zero.
4041  if (z1*z2 < 0) {
4042  *etheta = 0.0;
4043  *etheta_prime = 0.0;
4044  } else {
4045  // Actually calculate the interaction.
4046  double f1 = (double)i / (2.0 * j);
4047  double f2 = (double)j / (2.0 * i);
4048  *etheta = elambda[i*j] - f1*elambda[j*j] - f2*elambda[i*i];
4049  *etheta_prime = elambda1[i*j] - f1*elambda1[j*j] - f2*elambda1[i*i];
4050  }
4051 }
4052 
4053 void HMWSoln::s_updateIMS_lnMolalityActCoeff() const
4054 {
4055  // Calculate the molalities. Currently, the molalities may not be current
4056  // with respect to the contents of the State objects' data.
4057  calcMolalities();
4058  double xmolSolvent = moleFraction(m_indexSolvent);
4059  double xx = std::max(m_xmolSolventMIN, xmolSolvent);
4060  if (IMS_typeCutoff_ == 0) {
4061  for (size_t k = 1; k < m_kk; k++) {
4062  IMS_lnActCoeffMolal_[k]= 0.0;
4063  }
4064  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
4065  return;
4066  } else if (IMS_typeCutoff_ == 1) {
4067  if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
4068  for (size_t k = 1; k < m_kk; k++) {
4069  IMS_lnActCoeffMolal_[k]= 0.0;
4070  }
4071  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
4072  return;
4073  } else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
4074  double tmp = log(xx * IMS_gamma_k_min_);
4075  for (size_t k = 1; k < m_kk; k++) {
4076  IMS_lnActCoeffMolal_[k]= tmp;
4077  }
4078  IMS_lnActCoeffMolal_[m_indexSolvent] = log(IMS_gamma_o_min_);
4079  return;
4080  } else {
4081  // If we are in the middle region, calculate the connecting polynomials
4082  double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
4083  double xminus2 = xminus * xminus;
4084  double xminus3 = xminus2 * xminus;
4085  double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
4086  double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
4087 
4088  double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
4089  double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
4090 
4091  double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
4092  double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
4093 
4094  double h1_g = h1 / IMS_gamma_o_min_;
4095  double h1_g_prime = h1_prime / IMS_gamma_o_min_;
4096 
4097  double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
4098  double h1_f = h1 * alpha;
4099  double h1_f_prime = h1_prime * alpha;
4100 
4101  double f = h2 + h1_f;
4102  double f_prime = h2_prime + h1_f_prime;
4103 
4104  double g = h2 + h1_g;
4105  double g_prime = h2_prime + h1_g_prime;
4106 
4107  double tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
4108  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
4109  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
4110 
4111  tmp = log(xmolSolvent) + lngammak;
4112  for (size_t k = 1; k < m_kk; k++) {
4113  IMS_lnActCoeffMolal_[k]= tmp;
4114  }
4115  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
4116  }
4117  } else if (IMS_typeCutoff_ == 2) {
4118  // Exponentials - trial 2
4119  if (xmolSolvent > IMS_X_o_cutoff_) {
4120  for (size_t k = 1; k < m_kk; k++) {
4121  IMS_lnActCoeffMolal_[k]= 0.0;
4122  }
4123  IMS_lnActCoeffMolal_[m_indexSolvent] = - log(xx) + (xx - 1.0)/xx;
4124  return;
4125  } else {
4126  double xoverc = xmolSolvent/IMS_cCut_;
4127  double eterm = std::exp(-xoverc);
4128 
4129  double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
4130  + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
4131  double f_prime = 1.0 + eterm*fptmp;
4132  double f = xmolSolvent + IMS_efCut_
4133  + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
4134 
4135  double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
4136  + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
4137  double g_prime = 1.0 + eterm*gptmp;
4138  double g = xmolSolvent + IMS_egCut_
4139  + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
4140 
4141  double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
4142  double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
4143  double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
4144 
4145  tmp = log(xx) + lngammak;
4146  for (size_t k = 1; k < m_kk; k++) {
4147  IMS_lnActCoeffMolal_[k]= tmp;
4148  }
4149  IMS_lnActCoeffMolal_[m_indexSolvent] = lngammao;
4150  }
4151  }
4152  return;
4153 }
4154 
4155 void HMWSoln::printCoeffs() const
4156 {
4157  calcMolalities();
4158  vector_fp& moleF = m_tmpV;
4159 
4160  // Update the coefficients wrt Temperature. Calculate the derivatives as well
4161  s_updatePitzer_CoeffWRTemp(2);
4162  getMoleFractions(moleF.data());
4163 
4164  writelog("Index Name MoleF MolalityCropped Charge\n");
4165  for (size_t k = 0; k < m_kk; k++) {
4166  writelogf("%2d %-16s %14.7le %14.7le %5.1f \n",
4167  k, speciesName(k), moleF[k], m_molalitiesCropped[k], charge(k));
4168  }
4169 
4170  writelog("\n Species Species beta0MX "
4171  "beta1MX beta2MX CphiMX alphaMX thetaij\n");
4172  for (size_t i = 1; i < m_kk - 1; i++) {
4173  for (size_t j = i+1; j < m_kk; j++) {
4174  size_t n = i * m_kk + j;
4175  size_t ct = m_CounterIJ[n];
4176  writelogf(" %-16s %-16s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f \n",
4177  speciesName(i), speciesName(j),
4178  m_Beta0MX_ij[ct], m_Beta1MX_ij[ct],
4179  m_Beta2MX_ij[ct], m_CphiMX_ij[ct],
4180  m_Alpha1MX_ij[ct], m_Theta_ij[ct]);
4181  }
4182  }
4183 
4184  writelog("\n Species Species Species psi \n");
4185  for (size_t i = 1; i < m_kk; i++) {
4186  for (size_t j = 1; j < m_kk; j++) {
4187  for (size_t k = 1; k < m_kk; k++) {
4188  size_t n = k + j * m_kk + i * m_kk * m_kk;
4189  if (m_Psi_ijk[n] != 0.0) {
4190  writelogf(" %-16s %-16s %-16s %9.5f \n",
4191  speciesName(i), speciesName(j),
4192  speciesName(k), m_Psi_ijk[n]);
4193  }
4194  }
4195  }
4196  }
4197 }
4198 
4199 void HMWSoln::applyphScale(doublereal* acMolality) const
4200 {
4201  if (m_pHScalingType == PHSCALE_PITZER) {
4202  return;
4203  }
4204  AssertTrace(m_pHScalingType == PHSCALE_NBS);
4205  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4206  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4207  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4208  for (size_t k = 0; k < m_kk; k++) {
4209  acMolality[k] *= exp(charge(k) * afac);
4210  }
4211 }
4212 
4213 void HMWSoln::s_updateScaling_pHScaling() const
4214 {
4215  if (m_pHScalingType == PHSCALE_PITZER) {
4216  m_lnActCoeffMolal_Scaled = m_lnActCoeffMolal_Unscaled;
4217  return;
4218  }
4219  AssertTrace(m_pHScalingType == PHSCALE_NBS);
4220  doublereal lnGammaClMs2 = s_NBS_CLM_lnMolalityActCoeff();
4221  doublereal lnGammaCLMs1 = m_lnActCoeffMolal_Unscaled[m_indexCLM];
4222  doublereal afac = -1.0 *(lnGammaClMs2 - lnGammaCLMs1);
4223  for (size_t k = 0; k < m_kk; k++) {
4224  m_lnActCoeffMolal_Scaled[k] = m_lnActCoeffMolal_Unscaled[k] + charge(k) * afac;
4225  }
4226 }
4227 
4228 void HMWSoln::s_updateScaling_pHScaling_dT() const
4229 {
4230  if (m_pHScalingType == PHSCALE_PITZER) {
4231  m_dlnActCoeffMolaldT_Scaled = m_dlnActCoeffMolaldT_Unscaled;
4232  return;
4233  }
4234  AssertTrace(m_pHScalingType == PHSCALE_NBS);
4235  doublereal dlnGammaClM_dT_s2 = s_NBS_CLM_dlnMolalityActCoeff_dT();
4236  doublereal dlnGammaCLM_dT_s1 = m_dlnActCoeffMolaldT_Unscaled[m_indexCLM];
4237  doublereal afac = -1.0 *(dlnGammaClM_dT_s2 - dlnGammaCLM_dT_s1);
4238  for (size_t k = 0; k < m_kk; k++) {
4239  m_dlnActCoeffMolaldT_Scaled[k] = m_dlnActCoeffMolaldT_Unscaled[k] + charge(k) * afac;
4240  }
4241 }
4242 
4243 void HMWSoln::s_updateScaling_pHScaling_dT2() const
4244 {
4245  if (m_pHScalingType == PHSCALE_PITZER) {
4246  m_d2lnActCoeffMolaldT2_Scaled = m_d2lnActCoeffMolaldT2_Unscaled;
4247  return;
4248  }
4249  AssertTrace(m_pHScalingType == PHSCALE_NBS);
4250  doublereal d2lnGammaClM_dT2_s2 = s_NBS_CLM_d2lnMolalityActCoeff_dT2();
4251  doublereal d2lnGammaCLM_dT2_s1 = m_d2lnActCoeffMolaldT2_Unscaled[m_indexCLM];
4252  doublereal afac = -1.0 *(d2lnGammaClM_dT2_s2 - d2lnGammaCLM_dT2_s1);
4253  for (size_t k = 0; k < m_kk; k++) {
4254  m_d2lnActCoeffMolaldT2_Scaled[k] = m_d2lnActCoeffMolaldT2_Unscaled[k] + charge(k) * afac;
4255  }
4256 }
4257 
4258 void HMWSoln::s_updateScaling_pHScaling_dP() const
4259 {
4260  if (m_pHScalingType == PHSCALE_PITZER) {
4261  m_dlnActCoeffMolaldP_Scaled = m_dlnActCoeffMolaldP_Unscaled;
4262  return;
4263  }
4264  AssertTrace(m_pHScalingType == PHSCALE_NBS);
4265  doublereal dlnGammaClM_dP_s2 = s_NBS_CLM_dlnMolalityActCoeff_dP();
4266  doublereal dlnGammaCLM_dP_s1 = m_dlnActCoeffMolaldP_Unscaled[m_indexCLM];
4267  doublereal afac = -1.0 *(dlnGammaClM_dP_s2 - dlnGammaCLM_dP_s1);
4268  for (size_t k = 0; k < m_kk; k++) {
4269  m_dlnActCoeffMolaldP_Scaled[k] = m_dlnActCoeffMolaldP_Unscaled[k] + charge(k) * afac;
4270  }
4271 }
4272 
4273 doublereal HMWSoln::s_NBS_CLM_lnMolalityActCoeff() const
4274 {
4275  doublereal sqrtIs = sqrt(m_IionicMolality);
4276  doublereal A = A_Debye_TP();
4277  doublereal lnGammaClMs2 = - A * sqrtIs /(1.0 + 1.5 * sqrtIs);
4278  return lnGammaClMs2;
4279 }
4280 
4281 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dT() const
4282 {
4283  doublereal sqrtIs = sqrt(m_IionicMolality);
4284  doublereal dAdT = dA_DebyedT_TP();
4285  return - dAdT * sqrtIs /(1.0 + 1.5 * sqrtIs);
4286 }
4287 
4288 doublereal HMWSoln::s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
4289 {
4290  doublereal sqrtIs = sqrt(m_IionicMolality);
4291  doublereal d2AdT2 = d2A_DebyedT2_TP();
4292  return - d2AdT2 * sqrtIs /(1.0 + 1.5 * sqrtIs);
4293 }
4294 
4295 doublereal HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP() const
4296 {
4297  doublereal sqrtIs = sqrt(m_IionicMolality);
4298  doublereal dAdP = dA_DebyedP_TP();
4299  return - dAdP * sqrtIs /(1.0 + 1.5 * sqrtIs);
4300 }
4301 
4302 int HMWSoln::debugPrinting()
4303 {
4304  return m_debugCalc;
4305 }
4306 
4307 }
Cantera::HMWSoln::m_Theta_ij_P
vector_fp m_Theta_ij_P
Derivative of Theta_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2141
Cantera::HMWSoln::m_BMX_IJ_L
vector_fp m_BMX_IJ_L
Derivative of BMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2357
Cantera::HMWSoln::m_PhiPhi_IJ_P
vector_fp m_PhiPhi_IJ_P
Derivative of m_PhiPhi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2419
Cantera::VPStandardStateTP::updateStandardStateThermo
virtual void updateStandardStateThermo() const
Updates the standard state thermodynamic functions at the current T and P of the solution.
Definition: VPStandardStateTP.cpp:313
Cantera::HMWSoln::calcMolalitiesCropped
void calcMolalitiesCropped() const
Calculate the cropped molalities.
Definition: HMWSoln.cpp:1047
Cantera::HMWSoln::getChemPotentials
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: HMWSoln.cpp:574
Cantera::HMWSoln::s_updateScaling_pHScaling
void s_updateScaling_pHScaling() const
Apply the current phScale to a set of activity Coefficients.
Definition: HMWSoln.cpp:4213
Cantera::HMWSoln::m_CphiMX_ij_LL
vector_fp m_CphiMX_ij_LL
Derivative of Cphi_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2108
Cantera::Phase::molarVolume
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:676
Cantera::HMWSoln::s_update_dlnMolalityActCoeff_dT
void s_update_dlnMolalityActCoeff_dT() const
This function calculates the temperature derivative of the natural logarithm of the molality activity...
Definition: HMWSoln.cpp:2153
Cantera::HMWSoln::m_Theta_ij_L
vector_fp m_Theta_ij_L
Derivative of Theta_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2135
Cantera::HMWSoln::m_speciesCharge_Stoich
vector_fp m_speciesCharge_Stoich
Stoichiometric species charge -> This is for calculations of the ionic strength which ignore ion-ion ...
Definition: HMWSoln.h:2004
Cantera::HMWSoln::m_CMX_IJ_P
vector_fp m_CMX_IJ_P
Derivative of m_CMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2432
Cantera::HMWSoln::m_gfunc_IJ
vector_fp m_gfunc_IJ
Various temporary arrays used in the calculation of the Pitzer activity coefficients.
Definition: HMWSoln.h:2339
Cantera::HMWSoln::cp_mole
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
Definition: HMWSoln.cpp:465
Cantera::ValueCache::getId
int getId()
Get a unique id for a cached value.
Definition: ValueCache.cpp:21
Cantera::HMWSoln::m_Beta2MX_ij_L
vector_fp m_Beta2MX_ij_L
Derivative of Beta2_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2064
Cantera::HMWSoln::m_d2lnActCoeffMolaldT2_Unscaled
vector_fp m_d2lnActCoeffMolaldT2_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2299
Cantera::ThermoPhase::thermalExpansionCoeff
virtual doublereal thermalExpansionCoeff() const
Return the volumetric thermal expansion coefficient. Units: 1/K.
Definition: ThermoPhase.h:304
Cantera::HMWSoln::MC_X_o_cutoff_
doublereal MC_X_o_cutoff_
value of the solvent mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2486
Cantera::HMWSoln::applyphScale
virtual void applyphScale(doublereal *acMolality) const
Apply the current phScale to a set of activity Coefficients or activities.
Definition: HMWSoln.cpp:4199
Cantera::HMWSoln::IMS_lnActCoeffMolal_
vector_fp IMS_lnActCoeffMolal_
Logarithm of the molal activity coefficients.
Definition: HMWSoln.h:2440
Cantera::Phase::temperature
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
Cantera::Array2D::resize
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with &#39;v&#39;.
Definition: Array.h:112
Cantera::HMWSoln::m_Lambda_nj_coeff
Array2D m_Lambda_nj_coeff
Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2228
Cantera::HMWSoln::m_maxIionicStrength
double m_maxIionicStrength
Maximum value of the ionic strength allowed in the calculation of the activity coefficients.
Definition: HMWSoln.h:1914
Cantera::HMWSoln::entropy_mole
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Definition: HMWSoln.cpp:453
Cantera::HMWSoln::m_Beta1MX_ij_LL
vector_fp m_Beta1MX_ij_LL
Derivative of Beta1_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2043
Cantera::MolalityVPSSTP::m_Mnaught
doublereal m_Mnaught
This is the multiplication factor that goes inside log expressions involving the molalities of specie...
Definition: MolalityVPSSTP.h:606
Cantera::HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dP
doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const
Calculate the pressure derivative of the Chlorine activity coefficient.
Definition: HMWSoln.cpp:4295
Cantera::HMWSoln::m_Mu_nnn_L
vector_fp m_Mu_nnn_L
Mu coefficient temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2246
Cantera::HMWSoln::m_dlnActCoeffMolaldP_Unscaled
vector_fp m_dlnActCoeffMolaldP_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2307
Cantera::HMWSoln::m_Beta0MX_ij_coeff
Array2D m_Beta0MX_ij_coeff
Array of coefficients for Beta0, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:2030
Cantera::HMWSoln::m_Beta2MX_ij_P
vector_fp m_Beta2MX_ij_P
Derivative of Beta2_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2070
Cantera::HMWSoln::AionicRadius
double AionicRadius(int k=0) const
Reports the ionic radius of the kth species.
Definition: HMWSoln.cpp:844
Cantera::HMWSoln::m_Beta2MX_ij
vector_fp m_Beta2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2061
Cantera::Phase::moleFraction
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:547
Cantera::HMWSoln::m_Lambda_nj_P
Array2D m_Lambda_nj_P
Derivative of Lambda_nj[i][j] wrt P.
Definition: HMWSoln.h:2213
Cantera::npos
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
Cantera::HMWSoln::eosType
virtual int eosType() const
Equation of state type flag.
Definition: HMWSoln.cpp:372
Cantera::HMWSoln::m_PhiPhi_IJ
vector_fp m_PhiPhi_IJ
Intermediate variable called PhiPhi in Pitzer&#39;s paper.
Definition: HMWSoln.h:2410
Cantera::PDSS::density
virtual doublereal density() const
Return the standard state density at standard state.
Definition: PDSS.cpp:261
Cantera::HMWSoln::m_Phi_IJ_L
vector_fp m_Phi_IJ_L
Derivative of m_Phi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2396
electrolytes.h
Header file for a common definitions used in electrolytes thermodynamics.
Cantera::HMWSoln::IMS_typeCutoff_
int IMS_typeCutoff_
IMS Cutoff type.
Definition: HMWSoln.h:2443
ThermoFactory.h
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
Cantera::HMWSoln::getPartialMolarEntropies
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
Definition: HMWSoln.cpp:616
Cantera::HMWSoln::m_dlnActCoeffMolaldT_Scaled
vector_fp m_dlnActCoeffMolaldT_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2287
Cantera::PDSS::satPressure
virtual doublereal satPressure(doublereal T)
saturation pressure
Definition: PDSS.cpp:375
Cantera::HMWSoln::setMolarDensity
virtual void setMolarDensity(const doublereal conc)
Set the internally stored molar density (kmol/m^3) for the phase.
Definition: HMWSoln.cpp:511
Cantera::HMWSoln::counterIJ_setup
void counterIJ_setup() const
Set up a counter variable for keeping track of symmetric binary interactions amongst the solute speci...
Definition: HMWSoln.cpp:1188
Cantera::MolalityVPSSTP::m_indexCLM
size_t m_indexCLM
Index of the phScale species.
Definition: MolalityVPSSTP.h:589
Cantera::writelog
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
Cantera::HMWSoln::m_BMX_IJ_P
vector_fp m_BMX_IJ_P
Derivative of BMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2363
Cantera::HMWSoln::s_updateIMS_lnMolalityActCoeff
void s_updateIMS_lnMolalityActCoeff() const
This function will be called to update the internally stored natural logarithm of the molality activi...
Definition: HMWSoln.cpp:4053
Cantera::ThermoPhase::isothermalCompressibility
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
Definition: ThermoPhase.h:293
Cantera::HMWSoln::m_Beta0MX_ij_L
vector_fp m_Beta0MX_ij_L
Derivative of Beta0_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2016
Cantera::HMWSoln::elambda
double elambda[17]
This is elambda, MEC.
Definition: HMWSoln.h:2327
Cantera::XML_Node
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
PDSS_Water.h
Implementation of a pressure dependent standard state virtual function for a Pure Water Phase (see Sp...
Cantera::HMWSoln::m_dlnActCoeffMolaldP_Scaled
vector_fp m_dlnActCoeffMolaldP_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt P...
Definition: HMWSoln.h:2303
Cantera::HMWSoln::debugPrinting
int debugPrinting()
Return int specifying the amount of debug printing.
Definition: HMWSoln.cpp:4302
Cantera::HMWSoln::m_Mu_nnn_P
vector_fp m_Mu_nnn_P
Mu coefficient pressure derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2266
Cantera::HMWSoln::s_NBS_CLM_dlnMolalityActCoeff_dT
doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const
Calculate the temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:4281
Cantera::HMWSoln::m_Lambda_nj_L
Array2D m_Lambda_nj_L
Derivative of Lambda_nj[i][j] wrt T. see m_Lambda_ij.
Definition: HMWSoln.h:2207
Cantera::VPStandardStateTP::getCp_R
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the species standard states at the cu...
Definition: VPStandardStateTP.cpp:137
Cantera::HMWSoln::setDensity
virtual void setDensity(const doublereal rho)
Set the internally stored density (kg/m^3) of the phase.
Definition: HMWSoln.cpp:502
Cantera::HMWSoln::m_BprimeMX_IJ_P
vector_fp m_BprimeMX_IJ_P
Derivative of BprimeMX wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2376
Cantera::HMWSoln::m_CphiMX_ij_P
vector_fp m_CphiMX_ij_P
Derivative of Cphi_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2111
Cantera::warn_deprecated
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
Cantera::Phase::m_speciesSize
vector_fp m_speciesSize
Vector of species sizes.
Definition: Phase.h:798
Cantera::Phase::density
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
Cantera::HMWSoln::m_Psi_ijk_L
vector_fp m_Psi_ijk_L
Derivative of Psi_ijk[n] wrt T.
Definition: HMWSoln.h:2170
Cantera::HMWSoln::m_lnActCoeffMolal_Unscaled
vector_fp m_lnActCoeffMolal_Unscaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2283
Cantera::Array2D::ptrColumn
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:314
Cantera::HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2
void s_updatePitzer_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2795
Cantera::HMWSoln::m_Theta_ij
vector_fp m_Theta_ij
Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2132
Cantera::HMWSoln::m_hfunc_IJ
vector_fp m_hfunc_IJ
hfunc, was called gprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2346
Cantera::HMWSoln::elambda1
double elambda1[17]
This is elambda1, MEC.
Definition: HMWSoln.h:2330
Cantera::HMWSoln::m_Beta1MX_ij_P
vector_fp m_Beta1MX_ij_P
Derivative of Beta1_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2046
Cantera::HMWSoln::m_BprimeMX_IJ_LL
vector_fp m_BprimeMX_IJ_LL
Derivative of BprimeMX wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2373
Cantera::HMWSoln::ADebye_V
double ADebye_V(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_V.
Definition: HMWSoln.cpp:795
Cantera::HMWSoln::m_PhiPhi_IJ_LL
vector_fp m_PhiPhi_IJ_LL
Derivative of m_PhiPhi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2416
Cantera::HMWSoln::IMS_cCut_
doublereal IMS_cCut_
Parameter in the polyExp cutoff treatment having to do with rate of exp decay.
Definition: HMWSoln.h:2456
Cantera::HMWSoln::m_Psi_ijk_LL
vector_fp m_Psi_ijk_LL
Derivative of Psi_ijk[n] wrt TT.
Definition: HMWSoln.h:2174
Cantera::HMWSoln::m_IionicMolality
double m_IionicMolality
Current value of the ionic strength on the molality scale Associated Salts, if present in the mechani...
Definition: HMWSoln.h:1910
Cantera::HMWSoln::duplMyselfAsThermoPhase
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine for objects which inherit from ThermoPhase.
Definition: HMWSoln.cpp:367
Cantera::HMWSoln::m_Phi_IJ
vector_fp m_Phi_IJ
Intermediate variable called Phi in Pitzer&#39;s paper.
Definition: HMWSoln.h:2393
Cantera::HMWSoln::m_BphiMX_IJ
vector_fp m_BphiMX_IJ
Intermediate variable called BphiMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2380
Cantera::HMWSoln::m_CphiMX_ij_coeff
Array2D m_CphiMX_ij_coeff
Array of coefficients for CphiMX, a parameter in the activity coefficient formulation.
Definition: HMWSoln.h:2120
Cantera::HMWSoln::getPartialMolarVolumes
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
Definition: HMWSoln.cpp:653
Cantera::HMWSoln::s_updateScaling_pHScaling_dP
void s_updateScaling_pHScaling_dP() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt pressure...
Definition: HMWSoln.cpp:4258
Cantera::Phase::mean_X
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:690
Cantera::HMWSoln::m_Beta1MX_ij
vector_fp m_Beta1MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2037
Cantera::HMWSoln::m_tmpV
vector_fp m_tmpV
vector of size m_kk, used as a temporary holding area.
Definition: HMWSoln.h:1989
Cantera::ThermoPhase::RT
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:809
Cantera::HMWSoln::getPartialMolarCp
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
Definition: HMWSoln.cpp:666
Cantera::ThermoPhase
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
Cantera::HMWSoln::HMWSoln
HMWSoln()
Default Constructor.
Definition: HMWSoln.cpp:27
Cantera::HMWSoln::m_lnActCoeffMolal_Scaled
vector_fp m_lnActCoeffMolal_Scaled
Logarithm of the activity coefficients on the molality scale.
Definition: HMWSoln.h:2276
Cantera::HMWSoln::MC_X_o_min_
doublereal MC_X_o_min_
gamma_o value for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2489
Cantera::HMWSoln::m_Phi_IJ_P
vector_fp m_Phi_IJ_P
Derivative of m_Phi_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2402
Cantera::HMWSoln::m_waterSS
PDSS * m_waterSS
Water standard state calculator.
Definition: HMWSoln.h:1977
Cantera::HMWSoln::m_Phiprime_IJ
vector_fp m_Phiprime_IJ
Intermediate variable called Phiprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2406
Cantera::HMWSoln::m_h2func_IJ
vector_fp m_h2func_IJ
hfunc2, was called gprime in Pitzer&#39;s paper.
Definition: HMWSoln.h:2350
Cantera::HMWSoln::printCoeffs
void printCoeffs() const
Print out all of the input Pitzer coefficients.
Definition: HMWSoln.cpp:4155
Cantera::HMWSoln::m_Psi_ijk_P
vector_fp m_Psi_ijk_P
Derivative of Psi_ijk[n] wrt P.
Definition: HMWSoln.h:2178
Cantera::HMWSoln::m_CMX_IJ_L
vector_fp m_CMX_IJ_L
Derivative of m_CMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2426
Cantera::HMWSoln::m_molalitiesAreCropped
bool m_molalitiesAreCropped
Boolean indicating whether the molalities are cropped or are modified.
Definition: HMWSoln.h:2316
Cantera::HMWSoln::m_dlnActCoeffMolaldT_Unscaled
vector_fp m_dlnActCoeffMolaldT_Unscaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt T...
Definition: HMWSoln.h:2291
Cantera::HMWSoln::m_Psi_ijk_coeff
Array2D m_Psi_ijk_coeff
Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2194
Cantera::HMWSoln::s_NBS_CLM_lnMolalityActCoeff
doublereal s_NBS_CLM_lnMolalityActCoeff() const
Calculate the Chlorine activity coefficient on the NBS scale.
Definition: HMWSoln.cpp:4273
HMWSoln.h
Headers for the HMWSoln ThermoPhase object, which models concentrated electrolyte solutions (see Ther...
Cantera::HMWSoln::m_waterProps
std::unique_ptr< WaterProps > m_waterProps
Pointer to the water property calculator.
Definition: HMWSoln.h:1986
Cantera::HMWSoln::relative_enthalpy
virtual doublereal relative_enthalpy() const
Excess molar enthalpy of the solution from the mixing process.
Definition: HMWSoln.cpp:401
Cantera::HMWSoln::calc_thetas
void calc_thetas(int z1, int z2, double *etheta, double *etheta_prime) const
Calculate etheta and etheta_prime.
Definition: HMWSoln.cpp:4027
Cantera::HMWSoln::m_Mu_nnn
vector_fp m_Mu_nnn
Mu coefficient for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2236
Cantera::HMWSoln::m_molalitiesCropped
vector_fp m_molalitiesCropped
Cropped and modified values of the molalities used in activity coefficient calculations.
Definition: HMWSoln.h:2313
PITZERFORM_BASE
#define PITZERFORM_BASE
Major Parameters: The form of the Pitzer expression refers to the form of the Gibbs free energy expre...
Definition: HMWSoln.h:37
Cantera::HMWSoln::m_Beta2MX_ij_LL
vector_fp m_Beta2MX_ij_LL
Derivative of Beta2_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2067
Cantera::HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT
void s_updatePitzer_dlnMolalityActCoeff_dT() const
Calculates the temperature derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:2181
Cantera::ValueCache::getScalar
CachedScalar getScalar(int id)
Get a reference to a CachedValue object representing a scalar (doublereal) with the given id...
Definition: ValueCache.h:165
Cantera::HMWSoln::CROP_speciesCropped_
vector_int CROP_speciesCropped_
This is a boolean-type vector indicating whether a species&#39;s activity coefficient is in the cropped r...
Definition: HMWSoln.h:2519
Cantera::Phase::speciesName
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:267
Cantera::HMWSoln::m_CMX_IJ_LL
vector_fp m_CMX_IJ_LL
Derivative of m_CMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2429
Cantera::HMWSoln::m_Beta2MX_ij_coeff
Array2D m_Beta2MX_ij_coeff
Array of coefficients for Beta2, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:2080
Cantera::HMWSoln::m_Lambda_nj_LL
Array2D m_Lambda_nj_LL
Derivative of Lambda_nj[i][j] wrt TT.
Definition: HMWSoln.h:2210
Cantera::HMWSoln::getUnscaledMolalityActivityCoefficients
void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const
Get the array of unscaled non-dimensional molality based activity coefficients at the current solutio...
Definition: HMWSoln.cpp:561
Cantera::HMWSoln::IMS_slopegCut_
doublereal IMS_slopegCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:2470
Cantera::HMWSoln::IMS_slopefCut_
doublereal IMS_slopefCut_
Parameter in the polyExp cutoff treatment.
Definition: HMWSoln.h:2463
Cantera::HMWSoln::m_Beta0MX_ij_LL
vector_fp m_Beta0MX_ij_LL
Derivative of Beta0_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2019
Cantera::HMWSoln::s_updatePitzer_CoeffWRTemp
void s_updatePitzer_CoeffWRTemp(int doDerivs=2) const
Calculates the Pitzer coefficients&#39; dependence on the temperature.
Definition: HMWSoln.cpp:1211
Cantera::PHSCALE_NBS
const int PHSCALE_NBS
Scale to be used for evaluation of single-ion activity coefficients is that used by the NBS standard ...
Definition: MolalityVPSSTP.h:657
Cantera::HMWSoln::m_form_A_Debye
int m_form_A_Debye
Form of the constant outside the Debye-Huckel term called A.
Definition: HMWSoln.h:1939
Cantera::Phase::m_cache
ValueCache m_cache
Cached for saved calculations within each ThermoPhase.
Definition: Phase.h:761
Cantera::HMWSoln::s_NBS_CLM_d2lnMolalityActCoeff_dT2
doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const
Calculate the second temperature derivative of the Chlorine activity coefficient on the NBS scale...
Definition: HMWSoln.cpp:4288
Cantera::HMWSoln::enthalpy_mole
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
Definition: HMWSoln.cpp:395
Cantera::HMWSoln::m_CounterIJ
vector_int m_CounterIJ
a counter variable for keeping track of symmetric binary interactions amongst the solute species...
Definition: HMWSoln.h:2324
Cantera::HMWSoln::MC_slopepCut_
doublereal MC_slopepCut_
Parameter in the Molality Exp cutoff treatment.
Definition: HMWSoln.h:2496
Cantera::VPStandardStateTP::getStandardChemPotentials
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity for the species at their standard states at the...
Definition: VPStandardStateTP.cpp:99
Cantera::CanteraError
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
Cantera::HMWSoln::m_CMX_IJ
vector_fp m_CMX_IJ
Intermediate variable called CMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2423
Cantera::Phase::stateMFNumber
int stateMFNumber() const
Return the State Mole Fraction Number.
Definition: Phase.h:747
Cantera::HMWSoln::getPartialMolarEnthalpies
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
Definition: HMWSoln.cpp:597
Cantera::HMWSoln::m_TempPitzerRef
double m_TempPitzerRef
Reference Temperature for the Pitzer formulations.
Definition: HMWSoln.h:1917
Cantera::HMWSoln::m_CphiMX_ij_L
vector_fp m_CphiMX_ij_L
Derivative of Cphi_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2105
Cantera::PHSCALE_PITZER
const int PHSCALE_PITZER
Scale to be used for the output of single-ion activity coefficients is that used by Pitzer...
Definition: MolalityVPSSTP.h:632
Cantera::HMWSoln::m_BphiMX_IJ_L
vector_fp m_BphiMX_IJ_L
Derivative of BphiMX_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2383
Cantera::HMWSoln::m_g2func_IJ
vector_fp m_g2func_IJ
This is the value of g2(x2) in Pitzer&#39;s papers. Vector index is counterIJ.
Definition: HMWSoln.h:2342
Cantera::MolalityVPSSTP::m_pHScalingType
int m_pHScalingType
Scaling to be used for output of single-ion species activity coefficients.
Definition: MolalityVPSSTP.h:582
Cantera::HMWSoln::d2A_DebyedT2_TP
virtual double d2A_DebyedT2_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the 2nd derivative of the Debye Huckel constant with respect to temperature as a function of...
Definition: HMWSoln.cpp:818
Cantera::HMWSoln::m_densWaterSS
double m_densWaterSS
density of standard-state water
Definition: HMWSoln.h:1983
Cantera::Phase::getMoleFractions
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:542
Cantera::HMWSoln::m_Lambda_nj
Array2D m_Lambda_nj
Lambda coefficient for the ij interaction.
Definition: HMWSoln.h:2204
Cantera::HMWSoln::s_update_d2lnMolalityActCoeff_dT2
void s_update_d2lnMolalityActCoeff_dT2() const
This function calculates the temperature second derivative of the natural logarithm of the molality a...
Definition: HMWSoln.cpp:2767
Cantera::HMWSoln::cv_mole
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
Definition: HMWSoln.cpp:471
Cantera::importPhase
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
Definition: ThermoFactory.cpp:240
Cantera::HMWSoln::m_Phi_IJ_LL
vector_fp m_Phi_IJ_LL
Derivative of m_Phi_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2399
Cantera::HMWSoln::m_formPitzer
int m_formPitzer
This is the form of the Pitzer parameterization used in this model.
Definition: HMWSoln.h:1850
Cantera::VPStandardStateTP::pressure
virtual doublereal pressure() const
Returns the current pressure of the phase.
Definition: VPStandardStateTP.h:144
Cantera::HMWSoln::IMS_gamma_k_min_
doublereal IMS_gamma_k_min_
gamma_k minimum for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2453
Cantera::HMWSoln::getActivityConcentrations
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized activity concentrations.
Definition: HMWSoln.cpp:519
Cantera::HMWSoln::m_Alpha2MX_ij
vector_fp m_Alpha2MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2095
Cantera::HMWSoln::m_electrolyteSpeciesType
vector_int m_electrolyteSpeciesType
Vector containing the electrolyte species type.
Definition: HMWSoln.h:1902
Cantera::debuglog
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:161
Cantera::HMWSoln::m_Beta0MX_ij
vector_fp m_Beta0MX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2013
Cantera::HMWSoln::m_BprimeMX_IJ_L
vector_fp m_BprimeMX_IJ_L
Derivative of BprimeMX wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2370
Cantera::writelogf
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:200
Cantera::cHMWSoln0
const int cHMWSoln0
eosTypes returned for this ThermoPhase Object
Definition: electrolytes.h:41
Cantera::HMWSoln::IMS_gamma_o_min_
doublereal IMS_gamma_o_min_
gamma_o value for the cutoff process at the zero solvent point
Definition: HMWSoln.h:2450
AssertTrace
#define AssertTrace(expr)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:239
AssertThrowMsg
#define AssertThrowMsg(expr, procedure,...)
Assertion must be true or an error is thrown.
Definition: ctexceptions.h:270
Cantera::HMWSoln::m_BMX_IJ_LL
vector_fp m_BMX_IJ_LL
Derivative of BMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2360
Cantera::HMWSoln::m_formGC
int m_formGC
Format for the generalized concentration:
Definition: HMWSoln.h:1890
Cantera::MolalityVPSSTP::m_molalities
vector_fp m_molalities
Current value of the molalities of the species in the phase.
Definition: MolalityVPSSTP.h:610
Cantera::HMWSoln::m_A_Debye
double m_A_Debye
A_Debye: this expression appears on the top of the ln actCoeff term in the general Debye-Huckel expre...
Definition: HMWSoln.h:1971
Cantera::HMWSoln::standardConcentration
virtual doublereal standardConcentration(size_t k=0) const
Return the standard concentration for the kth species.
Definition: HMWSoln.cpp:532
Cantera::HMWSoln::initLengths
void initLengths()
Initialize all of the species-dependent lengths in the object.
Definition: HMWSoln.cpp:851
Cantera::SmallNumber
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
Cantera::vector_fp
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
Cantera::Phase::meanMolecularWeight
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:661
Cantera::HMWSoln::relative_molal_enthalpy
virtual doublereal relative_molal_enthalpy() const
Excess molar enthalpy of the solution from the mixing process on a molality basis.
Definition: HMWSoln.cpp:413
Cantera::CachedValue::validate
bool validate(double state1New)
Check whether the currently cached value is valid based on a single state variable.
Definition: ValueCache.h:44
Cantera::HMWSoln::calcDensity
void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
Definition: HMWSoln.cpp:483
Cantera::HMWSoln::m_gamma_tmp
vector_fp m_gamma_tmp
Intermediate storage of the activity coefficient itself.
Definition: HMWSoln.h:2436
Cantera::GasConstant
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Cantera::HMWSoln::s_update_dlnMolalityActCoeff_dP
void s_update_dlnMolalityActCoeff_dP() const
This function calculates the pressure derivative of the natural logarithm of the molality activity co...
Definition: HMWSoln.cpp:3373
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
Cantera::HMWSoln::getActivities
virtual void getActivities(doublereal *ac) const
Get the array of non-dimensional activities at the current solution temperature, pressure, and solution concentration.
Definition: HMWSoln.cpp:542
Cantera::MolalityVPSSTP::m_weightSolvent
doublereal m_weightSolvent
Molecular weight of the Solvent.
Definition: MolalityVPSSTP.h:592
Cantera::HMWSoln::IMS_X_o_cutoff_
doublereal IMS_X_o_cutoff_
value of the solute mole fraction that centers the cutoff polynomials for the cutoff =1 process; ...
Definition: HMWSoln.h:2447
Cantera::HMWSoln::m_Psi_ijk
vector_fp m_Psi_ijk
Array of 3D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2166
Cantera::HMWSoln
Class HMWSoln represents a dilute or concentrated liquid electrolyte phase which obeys the Pitzer for...
Definition: HMWSoln.h:1145
Cantera::HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP
void s_updatePitzer_dlnMolalityActCoeff_dP() const
Calculates the Pressure derivative of the natural logarithm of the molality activity coefficients...
Definition: HMWSoln.cpp:3397
Cantera::HMWSoln::m_BMX_IJ
vector_fp m_BMX_IJ
Intermediate variable called BMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2354
Cantera::HMWSoln::s_updateScaling_pHScaling_dT
void s_updateScaling_pHScaling_dT() const
Apply the current phScale to a set of derivatives of the activity Coefficients wrt temperature...
Definition: HMWSoln.cpp:4228
Cantera::Phase::m_kk
size_t m_kk
Number of species in the phase.
Definition: Phase.h:784
Cantera::VPStandardStateTP::getEntropy_R
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Entropy functions for the standard state species at the current T and...
Definition: VPStandardStateTP.cpp:113
Cantera::HMWSoln::dA_DebyedP_TP
virtual double dA_DebyedP_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to pressure, as a function of tempe...
Definition: HMWSoln.cpp:752
Cantera::HMWSoln::dA_DebyedT_TP
virtual double dA_DebyedT_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the derivative of the Debye Huckel constant with respect to temperature as a function of tem...
Definition: HMWSoln.cpp:728
Cantera::HMWSoln::m_IionicMolalityStoich
double m_IionicMolalityStoich
Stoichiometric ionic strength on the molality scale.
Definition: HMWSoln.h:1922
Cantera::HMWSoln::s_updateScaling_pHScaling_dT2
void s_updateScaling_pHScaling_dT2() const
Apply the current phScale to a set of 2nd derivatives of the activity Coefficients wrt temperature...
Definition: HMWSoln.cpp:4243
Cantera::HMWSoln::m_debugCalc
int m_debugCalc
Turn on copious debug printing when this is true.
Definition: HMWSoln.h:2774
Cantera::ThermoPhase::initThermoFile
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
Definition: ThermoPhase.cpp:655
Cantera::HMWSoln::m_BprimeMX_IJ
vector_fp m_BprimeMX_IJ
Intermediate variable called BprimeMX in Pitzer&#39;s paper.
Definition: HMWSoln.h:2367
Cantera::CachedValue::value
T value
The value of the cached property.
Definition: ValueCache.h:118
Cantera::HMWSoln::m_d2lnActCoeffMolaldT2_Scaled
vector_fp m_d2lnActCoeffMolaldT2_Scaled
Derivative of the Logarithm of the activity coefficients on the molality scale wrt TT...
Definition: HMWSoln.h:2295
Cantera::MolalityVPSSTP::calcMolalities
void calcMolalities() const
Calculates the molality of all species and stores the result internally.
Definition: MolalityVPSSTP.cpp:118
Cantera
Namespace for the Cantera kernel.
Definition: application.cpp:29
Cantera::HMWSoln::m_CphiMX_ij
vector_fp m_CphiMX_ij
Array of 2D data used in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2102
Cantera::MolalityVPSSTP::m_indexSolvent
size_t m_indexSolvent
Index of the solvent.
Definition: MolalityVPSSTP.h:573
Cantera::HMWSoln::m_formPitzerTemp
int m_formPitzerTemp
This is the form of the temperature dependence of Pitzer parameterization used in the model...
Definition: HMWSoln.h:1860
Cantera::HMWSoln::calc_lambdas
void calc_lambdas(double is) const
Calculate the lambda interactions.
Definition: HMWSoln.cpp:3976
Cantera::HMWSoln::m_Beta1MX_ij_L
vector_fp m_Beta1MX_ij_L
Derivative of Beta1_ij[i][j] wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2040
Cantera::HMWSoln::ADebye_J
double ADebye_J(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_J.
Definition: HMWSoln.cpp:806
Cantera::HMWSoln::m_Theta_ij_LL
vector_fp m_Theta_ij_LL
Derivative of Theta_ij[i][j] wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2138
Cantera::PDSS::setState_TP
virtual void setState_TP(doublereal temp, doublereal pres)
Set the internal temperature and pressure.
Definition: PDSS.cpp:365
Cantera::HMWSoln::m_Theta_ij_coeff
Array2D m_Theta_ij_coeff
Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
Definition: HMWSoln.h:2153
Cantera::HMWSoln::satPressure
virtual doublereal satPressure(doublereal T)
Get the saturation pressure for a given temperature.
Definition: HMWSoln.cpp:686
Cantera::HMWSoln::m_Mu_nnn_coeff
Array2D m_Mu_nnn_coeff
Array of coefficients form_Mu_nnn term.
Definition: HMWSoln.h:2269
Cantera::HMWSoln::m_Mu_nnn_LL
vector_fp m_Mu_nnn_LL
Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient.
Definition: HMWSoln.h:2256
Cantera::HMWSoln::s_updatePitzer_lnMolalityActCoeff
void s_updatePitzer_lnMolalityActCoeff() const
Calculate the Pitzer portion of the activity coefficients.
Definition: HMWSoln.cpp:1457
Cantera::HMWSoln::m_Aionic
vector_fp m_Aionic
a_k = Size of the ionic species in the DH formulation. units = meters
Definition: HMWSoln.h:1905
Cantera::VPStandardStateTP::getEnthalpy_RT
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
Definition: VPStandardStateTP.cpp:107
Cantera::Phase::charge
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:577
Cantera::HMWSoln::m_BphiMX_IJ_P
vector_fp m_BphiMX_IJ_P
Derivative of BphiMX_IJ wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2389
Cantera::HMWSoln::A_Debye_TP
virtual double A_Debye_TP(double temperature=-1.0, double pressure=-1.0) const
Value of the Debye Huckel constant as a function of temperature and pressure.
Definition: HMWSoln.cpp:696
Cantera::VPStandardStateTP::getStandardVolumes
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of the species standard states at the current T and P of the solution...
Definition: VPStandardStateTP.cpp:143
Cantera::HMWSoln::m_Beta1MX_ij_coeff
Array2D m_Beta1MX_ij_coeff
Array of coefficients for Beta1, a variable in Pitzer&#39;s papers.
Definition: HMWSoln.h:2054
Cantera::HMWSoln::m_PhiPhi_IJ_L
vector_fp m_PhiPhi_IJ_L
Derivative of m_PhiPhi_IJ wrt T. Vector index is counterIJ.
Definition: HMWSoln.h:2413
Cantera::Phase::setDensity
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
Definition: Phase.h:622
Cantera::HMWSoln::ADebye_L
double ADebye_L(double temperature=-1.0, double pressure=-1.0) const
Return Pitzer&#39;s definition of A_L.
Definition: HMWSoln.cpp:784
Cantera::HMWSoln::m_Beta0MX_ij_P
vector_fp m_Beta0MX_ij_P
Derivative of Beta0_ij[i][j] wrt P. Vector index is counterIJ.
Definition: HMWSoln.h:2022
Cantera::HMWSoln::m_BphiMX_IJ_LL
vector_fp m_BphiMX_IJ_LL
Derivative of BphiMX_IJ wrt TT. Vector index is counterIJ.
Definition: HMWSoln.h:2386
Cantera::HMWSoln::gibbs_mole
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: HMWSoln.cpp:459
Generated by   doxygen 1.8.14