Cantera  2.0
MixedSolventElectrolyte.cpp
1 /**
2  * @file MargulesVPSSTP.cpp
3  * Definitions for ThermoPhase object for phases which
4  * employ excess gibbs free energy formulations related to Margules
5  * expansions (see \ref thermoprops
6  * and class \link Cantera::MargulesVPSSTP MargulesVPSSTP\endlink).
7  *
8  */
9 /*
10  * Copyright (2009) Sandia Corporation. Under the terms of
11  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
12  * U.S. Government retains certain rights in this software.
13  */
14 
15 #include "cantera/thermo/MixedSolventElectrolyte.h"
18 
19 #include <iomanip>
20 #include <fstream>
21 
22 using namespace std;
23 
24 namespace Cantera
25 {
26 
27 static const double xxSmall = 1.0E-150;
28 /*
29  * Default constructor.
30  *
31  */
32 MixedSolventElectrolyte::MixedSolventElectrolyte() :
34  numBinaryInteractions_(0),
35  formMargules_(0),
36  formTempModel_(0)
37 {
38 }
39 
40 /*
41  * Working constructors
42  *
43  * The two constructors below are the normal way
44  * the phase initializes itself. They are shells that call
45  * the routine initThermo(), with a reference to the
46  * XML database to get the info for the phase.
47 
48  */
49 MixedSolventElectrolyte::MixedSolventElectrolyte(std::string inputFile, std::string id) :
51  numBinaryInteractions_(0),
52  formMargules_(0),
53  formTempModel_(0)
54 {
55  constructPhaseFile(inputFile, id);
56 }
57 
60  numBinaryInteractions_(0),
61  formMargules_(0),
62  formTempModel_(0)
63 {
64  constructPhaseXML(phaseRoot, id);
65 }
66 
67 
68 /*
69  * Copy Constructor:
70  *
71  * Note this stuff will not work until the underlying phase
72  * has a working copy constructor
73  */
76 {
78 }
79 
80 /*
81  * operator=()
82  *
83  * Note this stuff will not work until the underlying phase
84  * has a working assignment operator
85  */
88 {
89  if (&b == this) {
90  return *this;
91  }
92 
94 
96  m_HE_b_ij = b.m_HE_b_ij;
97  m_HE_c_ij = b.m_HE_c_ij;
98  m_HE_d_ij = b.m_HE_d_ij;
99  m_SE_b_ij = b.m_SE_b_ij;
100  m_SE_c_ij = b.m_SE_c_ij;
101  m_SE_d_ij = b.m_SE_d_ij;
112 
113  return *this;
114 }
115 
116 /**
117  *
118  * ~MixedSolventElectrolyte(): (virtual)
119  *
120  * Destructor: does nothing:
121  *
122  */
124 {
125 }
126 
127 /*
128  * This routine duplicates the current object and returns
129  * a pointer to ThermoPhase.
130  */
133 {
135  return (ThermoPhase*) mtp;
136 }
137 
138 // Special constructor for a hard-coded problem
139 /*
140  *
141  * LiKCl treating the PseudoBinary layer as passthrough.
142  * -> test to predict the eutectic and liquidus correctly.
143  *
144  */
147  numBinaryInteractions_(0),
148  formMargules_(0),
149  formTempModel_(0)
150 {
151 
152 
153  constructPhaseFile("LiKCl_liquid.xml", "");
154 
155 
157 
158  m_HE_b_ij.resize(1);
159  m_HE_c_ij.resize(1);
160  m_HE_d_ij.resize(1);
161 
162  m_SE_b_ij.resize(1);
163  m_SE_c_ij.resize(1);
164  m_SE_d_ij.resize(1);
165 
166  m_VHE_b_ij.resize(1);
167  m_VHE_c_ij.resize(1);
168  m_VHE_d_ij.resize(1);
169 
170  m_VSE_b_ij.resize(1);
171  m_VSE_c_ij.resize(1);
172  m_VSE_d_ij.resize(1);
173 
174  m_pSpecies_A_ij.resize(1);
175  m_pSpecies_B_ij.resize(1);
176 
177 
178 
179  m_HE_b_ij[0] = -17570E3;
180  m_HE_c_ij[0] = -377.0E3;
181  m_HE_d_ij[0] = 0.0;
182 
183  m_SE_b_ij[0] = -7.627E3;
184  m_SE_c_ij[0] = 4.958E3;
185  m_SE_d_ij[0] = 0.0;
186 
187 
188  size_t iLiCl = speciesIndex("LiCl(L)");
189  if (iLiCl == npos) {
190  throw CanteraError("MixedSolventElectrolyte test1 constructor",
191  "Unable to find LiCl(L)");
192  }
193  m_pSpecies_B_ij[0] = iLiCl;
194 
195 
196  size_t iKCl = speciesIndex("KCl(L)");
197  if (iKCl == npos) {
198  throw CanteraError("MixedSolventElectrolyte test1 constructor",
199  "Unable to find KCl(L)");
200  }
201  m_pSpecies_A_ij[0] = iKCl;
202 }
203 
204 
205 /*
206  * -------------- Utilities -------------------------------
207  */
208 
209 
210 // Equation of state type flag.
211 /*
212  * The ThermoPhase base class returns
213  * zero. Subclasses should define this to return a unique
214  * non-zero value. Known constants defined for this purpose are
215  * listed in mix_defs.h. The MixedSolventElectrolyte class also returns
216  * zero, as it is a non-complete class.
217  */
219 {
220  return 0;
221 }
222 
223 /*
224  * Import, construct, and initialize a phase
225  * specification from an XML tree into the current object.
226  *
227  * This routine is a precursor to constructPhaseXML(XML_Node*)
228  * routine, which does most of the work.
229  *
230  * @param infile XML file containing the description of the
231  * phase
232  *
233  * @param id Optional parameter identifying the name of the
234  * phase. If none is given, the first XML
235  * phase element will be used.
236  */
237 void MixedSolventElectrolyte::constructPhaseFile(std::string inputFile, std::string id)
238 {
239 
240  if ((int) inputFile.size() == 0) {
241  throw CanteraError("MixedSolventElectrolyte:constructPhaseFile",
242  "input file is null");
243  }
244  string path = findInputFile(inputFile);
245  std::ifstream fin(path.c_str());
246  if (!fin) {
247  throw CanteraError("MixedSolventElectrolyte:constructPhaseFile","could not open "
248  +path+" for reading.");
249  }
250  /*
251  * The phase object automatically constructs an XML object.
252  * Use this object to store information.
253  */
254  XML_Node& phaseNode_XML = xml();
255  XML_Node* fxml = new XML_Node();
256  fxml->build(fin);
257  XML_Node* fxml_phase = findXMLPhase(fxml, id);
258  if (!fxml_phase) {
259  throw CanteraError("MixedSolventElectrolyte:constructPhaseFile",
260  "ERROR: Can not find phase named " +
261  id + " in file named " + inputFile);
262  }
263  fxml_phase->copy(&phaseNode_XML);
264  constructPhaseXML(*fxml_phase, id);
265  delete fxml;
266 }
267 
268 /*
269  * Import, construct, and initialize a HMWSoln phase
270  * specification from an XML tree into the current object.
271  *
272  * Most of the work is carried out by the cantera base
273  * routine, importPhase(). That routine imports all of the
274  * species and element data, including the standard states
275  * of the species.
276  *
277  * Then, In this routine, we read the information
278  * particular to the specification of the activity
279  * coefficient model for the Pitzer parameterization.
280  *
281  * We also read information about the molar volumes of the
282  * standard states if present in the XML file.
283  *
284  * @param phaseNode This object must be the phase node of a
285  * complete XML tree
286  * description of the phase, including all of the
287  * species data. In other words while "phase" must
288  * point to an XML phase object, it must have
289  * sibling nodes "speciesData" that describe
290  * the species in the phase.
291  * @param id ID of the phase. If nonnull, a check is done
292  * to see if phaseNode is pointing to the phase
293  * with the correct id.
294  */
295 void MixedSolventElectrolyte::constructPhaseXML(XML_Node& phaseNode, std::string id)
296 {
297  string stemp;
298  if ((int) id.size() > 0) {
299  string idp = phaseNode.id();
300  if (idp != id) {
301  throw CanteraError("MixedSolventElectrolyte::constructPhaseXML",
302  "phasenode and Id are incompatible");
303  }
304  }
305 
306  /*
307  * Find the Thermo XML node
308  */
309  if (!phaseNode.hasChild("thermo")) {
310  throw CanteraError("MixedSolventElectrolyte::constructPhaseXML",
311  "no thermo XML node");
312  }
313  XML_Node& thermoNode = phaseNode.child("thermo");
314 
315  /*
316  * Make sure that the thermo model is Margules
317  */
318  stemp = thermoNode.attrib("model");
319  string formString = lowercase(stemp);
320  if (formString != "margules") {
321  throw CanteraError("MixedSolventElectrolyte::constructPhaseXML",
322  "model name isn't Margules: " + formString);
323 
324  }
325 
326  /*
327  * Call the Cantera importPhase() function. This will import
328  * all of the species into the phase. This will also handle
329  * all of the solvent and solute standard states
330  */
331  bool m_ok = importPhase(phaseNode, this);
332  if (!m_ok) {
333  throw CanteraError("MixedSolventElectrolyte::constructPhaseXML","importPhase failed ");
334  }
335 
336 }
337 //====================================================================================================================
338 /*
339  * ------------ Molar Thermodynamic Properties ----------------------
340  */
341 //====================================================================================================================
342 /*
343  * - Activities, Standard States, Activity Concentrations -----------
344  */
345 //====================================================================================================================
346 // Get the array of non-dimensional molar-based activity coefficients at
347 // the current solution temperature, pressure, and solution concentration.
348 /*
349  * @param ac Output vector of activity coefficients. Length: m_kk.
350  */
352 {
353  /*
354  * Update the activity coefficients
355  */
357 
358  /*
359  * take the exp of the internally stored coefficients.
360  */
361  for (size_t k = 0; k < m_kk; k++) {
362  ac[k] = exp(lnActCoeff_Scaled_[k]);
363  }
364 }
365 //====================================================================================================================
366 /*
367  * ------------ Partial Molar Properties of the Solution ------------
368  */
369 
370 
371 
373 {
374  getChemPotentials(mu);
375  double ve = Faraday * electricPotential();
376  for (size_t k = 0; k < m_kk; k++) {
377  mu[k] += ve*charge(k);
378  }
379 }
380 
381 
383 {
384  doublereal xx;
385  /*
386  * First get the standard chemical potentials in
387  * molar form.
388  * -> this requires updates of standard state as a function
389  * of T and P
390  */
392  /*
393  * Update the activity coefficients
394  */
396  /*
397  *
398  */
399  doublereal RT = GasConstant * temperature();
400  for (size_t k = 0; k < m_kk; k++) {
401  xx = std::max(moleFractions_[k], xxSmall);
402  mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
403  }
404 }
405 
406 /// Molar enthalpy. Units: J/kmol.
408 {
409  size_t kk = nSpecies();
410  double h = 0;
411  vector_fp hbar(kk);
412  getPartialMolarEnthalpies(&hbar[0]);
413  for (size_t i = 0; i < kk; i++) {
414  h += moleFractions_[i]*hbar[i];
415  }
416  return h;
417 }
418 
419 /// Molar entropy. Units: J/kmol.
421 {
422  size_t kk = nSpecies();
423  double s = 0;
424  vector_fp sbar(kk);
425  getPartialMolarEntropies(&sbar[0]);
426  for (size_t i = 0; i < kk; i++) {
427  s += moleFractions_[i]*sbar[i];
428  }
429  return s;
430 }
431 
432 /// Molar heat capacity at constant pressure. Units: J/kmol/K.
434 {
435  size_t kk = nSpecies();
436  double cp = 0;
437  vector_fp cpbar(kk);
438  getPartialMolarCp(&cpbar[0]);
439  for (size_t i = 0; i < kk; i++) {
440  cp += moleFractions_[i]*cpbar[i];
441  }
442  return cp;
443 }
444 
445 /// Molar heat capacity at constant volume. Units: J/kmol/K.
447 {
448  return cp_mole() - GasConstant;
449 }
450 
451 // Returns an array of partial molar enthalpies for the species
452 // in the mixture.
453 /*
454  * Units (J/kmol)
455  *
456  * For this phase, the partial molar enthalpies are equal to the
457  * standard state enthalpies modified by the derivative of the
458  * molality-based activity coefficient wrt temperature
459  *
460  * \f[
461  * \bar h_k(T,P) = h^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
462  * \f]
463  *
464  */
466 {
467  /*
468  * Get the nondimensional standard state enthalpies
469  */
470  getEnthalpy_RT(hbar);
471  /*
472  * dimensionalize it.
473  */
474  double T = temperature();
475  double RT = GasConstant * T;
476  for (size_t k = 0; k < m_kk; k++) {
477  hbar[k] *= RT;
478  }
479  /*
480  * Update the activity coefficients, This also update the
481  * internally stored molalities.
482  */
485  double RTT = RT * T;
486  for (size_t k = 0; k < m_kk; k++) {
487  hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
488  }
489 }
490 
491 // Returns an array of partial molar heat capacities for the species
492 // in the mixture.
493 /*
494  * Units (J/kmol)
495  *
496  * For this phase, the partial molar enthalpies are equal to the
497  * standard state enthalpies modified by the derivative of the
498  * activity coefficient wrt temperature
499  *
500  * \f[
501  * ??????????? \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
502  * \f]
503  *
504  */
505 void MixedSolventElectrolyte::getPartialMolarCp(doublereal* cpbar) const
506 {
507  /*
508  * Get the nondimensional standard state entropies
509  */
510  getCp_R(cpbar);
511  double T = temperature();
512  /*
513  * Update the activity coefficients, This also update the
514  * internally stored molalities.
515  */
518 
519  for (size_t k = 0; k < m_kk; k++) {
520  cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
521  }
522  /*
523  * dimensionalize it.
524  */
525  for (size_t k = 0; k < m_kk; k++) {
526  cpbar[k] *= GasConstant;
527  }
528 }
529 
530 // Returns an array of partial molar entropies for the species
531 // in the mixture.
532 /*
533  * Units (J/kmol)
534  *
535  * For this phase, the partial molar enthalpies are equal to the
536  * standard state enthalpies modified by the derivative of the
537  * activity coefficient wrt temperature
538  *
539  * \f[
540  * \bar s_k(T,P) = s^o_k(T,P) - R T^2 \frac{d \ln(\gamma_k)}{dT}
541  * \f]
542  *
543  */
545 {
546  double xx;
547  /*
548  * Get the nondimensional standard state entropies
549  */
550  getEntropy_R(sbar);
551  double T = temperature();
552  /*
553  * Update the activity coefficients, This also update the
554  * internally stored molalities.
555  */
558 
559  for (size_t k = 0; k < m_kk; k++) {
560  xx = std::max(moleFractions_[k], xxSmall);
561  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
562  }
563  /*
564  * dimensionalize it.
565  */
566  for (size_t k = 0; k < m_kk; k++) {
567  sbar[k] *= GasConstant;
568  }
569 }
570 
571 /*
572  * ------------ Partial Molar Properties of the Solution ------------
573  */
574 
575 // Return an array of partial molar volumes for the
576 // species in the mixture. Units: m^3/kmol.
577 /*
578  * Frequently, for this class of thermodynamics representations,
579  * the excess Volume due to mixing is zero. Here, we set it as
580  * a default. It may be overridden in derived classes.
581  *
582  * @param vbar Output vector of species partial molar volumes.
583  * Length = m_kk. units are m^3/kmol.
584  */
586 {
587  int delAK, delBK;
588  double XA, XB, g0 , g1;
589  double T = temperature();
590 
591  /*
592  * Get the standard state values in m^3 kmol-1
593  */
594  getStandardVolumes(vbar);
595 
596  for (size_t iK = 0; iK < m_kk; iK++) {
597  delAK = 0;
598  delBK = 0;
599  for (size_t i = 0; i < numBinaryInteractions_; i++) {
600  size_t iA = m_pSpecies_A_ij[i];
601  size_t iB = m_pSpecies_B_ij[i];
602 
603  if (iA==iK) {
604  delAK = 1;
605  } else if (iB==iK) {
606  delBK = 1;
607  }
608 
609  XA = moleFractions_[iA];
610  XB = moleFractions_[iB];
611 
612  g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
613  g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
614 
615  vbar[iK] += XA*XB*(g0+g1*XB)+((delAK-XA)*XB+XA*(delBK-XB))*(g0+g1*XB)+XA*XB*(delBK-XB)*g1;
616  }
617  }
618 }
619 
620 doublereal MixedSolventElectrolyte::err(std::string msg) const
621 {
622  throw CanteraError("MixedSolventElectrolyte","Base class method "
623  +msg+" called. Equation of state type: "+int2str(eosType()));
624  return 0;
625 }
626 
627 
628 /*
629  * @internal Initialize. This method is provided to allow
630  * subclasses to perform any initialization required after all
631  * species have been added. For example, it might be used to
632  * resize internal work arrays that must have an entry for
633  * each species. The base class implementation does nothing,
634  * and subclasses that do not require initialization do not
635  * need to overload this method. When importing a CTML phase
636  * description, this method is called just prior to returning
637  * from function importPhase.
638  *
639  * @see importCTML.cpp
640  */
642 {
643  initLengths();
645 }
646 
647 
648 // Initialize lengths of local variables after all species have
649 // been identified.
651 {
652  m_kk = nSpecies();
654 }
655 
656 /*
657  * initThermoXML() (virtual from ThermoPhase)
658  * Import and initialize a ThermoPhase object
659  *
660  * @param phaseNode This object must be the phase node of a
661  * complete XML tree
662  * description of the phase, including all of the
663  * species data. In other words while "phase" must
664  * point to an XML phase object, it must have
665  * sibling nodes "speciesData" that describe
666  * the species in the phase.
667  * @param id ID of the phase. If nonnull, a check is done
668  * to see if phaseNode is pointing to the phase
669  * with the correct id.
670  */
671 void MixedSolventElectrolyte::initThermoXML(XML_Node& phaseNode, std::string id)
672 {
673  string subname = "MixedSolventElectrolyte::initThermoXML";
674  string stemp;
675 
676  /*
677  * Check on the thermo field. Must have:
678  * <thermo model="IdealSolidSolution" />
679  */
680 
681  XML_Node& thermoNode = phaseNode.child("thermo");
682  string mStringa = thermoNode.attrib("model");
683  string mString = lowercase(mStringa);
684  if (mString != "margules") {
685  throw CanteraError(subname.c_str(),
686  "Unknown thermo model: " + mStringa);
687  }
688 
689 
690  /*
691  * Go get all of the coefficients and factors in the
692  * activityCoefficients XML block
693  */
694  /*
695  * Go get all of the coefficients and factors in the
696  * activityCoefficients XML block
697  */
698  XML_Node* acNodePtr = 0;
699  if (thermoNode.hasChild("activityCoefficients")) {
700  XML_Node& acNode = thermoNode.child("activityCoefficients");
701  acNodePtr = &acNode;
702  string mStringa = acNode.attrib("model");
703  string mString = lowercase(mStringa);
704  if (mString != "margules") {
705  throw CanteraError(subname.c_str(),
706  "Unknown activity coefficient model: " + mStringa);
707  }
708  size_t n = acNodePtr->nChildren();
709  for (size_t i = 0; i < n; i++) {
710  XML_Node& xmlACChild = acNodePtr->child(i);
711  stemp = xmlACChild.name();
712  string nodeName = lowercase(stemp);
713  /*
714  * Process a binary salt field, or any of the other XML fields
715  * that make up the Pitzer Database. Entries will be ignored
716  * if any of the species in the entry isn't in the solution.
717  */
718  if (nodeName == "binaryneutralspeciesparameters") {
719  readXMLBinarySpecies(xmlACChild);
720 
721  }
722  }
723  }
724 
725  /*
726  * Go down the chain
727  */
728  MolarityIonicVPSSTP::initThermoXML(phaseNode, id);
729 
730 
731 }
732 //===================================================================================================================
733 
734 // Update the activity coefficients
735 /*
736  * This function will be called to update the internally stored
737  * natural logarithm of the activity coefficients
738  *
739  * he = X_A X_B(B + C X_B)
740  */
742 {
743  int delAK, delBK;
744  double XA, XB, g0, g1;
745  double T = temperature();
746  double RT = GasConstant*T;
747  lnActCoeff_Scaled_.assign(m_kk, 0.0);
748  for (size_t iK = 0; iK < m_kk; iK++) {
749  for (size_t i = 0; i < numBinaryInteractions_; i++) {
750  size_t iA = m_pSpecies_A_ij[i];
751  size_t iB = m_pSpecies_B_ij[i];
752  delAK = 0;
753  delBK = 0;
754  if (iA==iK) {
755  delAK = 1;
756  } else if (iB==iK) {
757  delBK = 1;
758  }
759  XA = moleFractions_[iA];
760  XB = moleFractions_[iB];
761  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
762  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
763  lnActCoeff_Scaled_[iK] += (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
764  }
765  }
766 }
767 //===================================================================================================================
768 // Update the derivative of the log of the activity coefficients wrt T
769 /*
770  * This function will be called to update the internally stored
771  * natural logarithm of the activity coefficients
772  *
773  * he = X_A X_B(B + C X_B)
774  */
776 {
777  int delAK, delBK;
778  doublereal XA, XB, g0, g1;
779  doublereal T = temperature();
780  doublereal RTT = GasConstant*T*T;
781  dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
782  d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0);
783  for (size_t iK = 0; iK < m_kk; iK++) {
784  for (size_t i = 0; i < numBinaryInteractions_; i++) {
785  size_t iA = m_pSpecies_A_ij[i];
786  size_t iB = m_pSpecies_B_ij[i];
787  delAK = 0;
788  delBK = 0;
789  if (iA==iK) {
790  delAK = 1;
791  } else if (iB==iK) {
792  delBK = 1;
793  }
794  XA = moleFractions_[iA];
795  XB = moleFractions_[iB];
796  g0 = -m_HE_b_ij[i] / RTT;
797  g1 = -m_HE_c_ij[i] / RTT;
798  double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
799  dlnActCoeffdT_Scaled_[iK] += temp;
800  d2lnActCoeffdT2_Scaled_[iK] -= 2.0 * temp / T;
801  }
802  }
803 }
804 //====================================================================================================================
805 void MixedSolventElectrolyte::getdlnActCoeffdT(doublereal* dlnActCoeffdT) const
806 {
808  for (size_t k = 0; k < m_kk; k++) {
809  dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
810  }
811 }
812 //====================================================================================================================
813 void MixedSolventElectrolyte::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const
814 {
816  for (size_t k = 0; k < m_kk; k++) {
817  d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
818  }
819 }
820 //====================================================================================================================
821 
822 // Get the change in activity coefficients w.r.t. change in state (temp, mole fraction, etc.) along
823 // a line in parameter space or along a line in physical space
824 /*
825  *
826  * @param dTds Input of temperature change along the path
827  * @param dXds Input vector of changes in mole fraction along the path. length = m_kk
828  * Along the path length it must be the case that the mole fractions sum to one.
829  * @param dlnActCoeffds Output vector of the directional derivatives of the
830  * log Activity Coefficients along the path. length = m_kk
831  * units are 1/units(s). if s is a physical coordinate then the units are 1/m.
832  */
833 void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
834  doublereal* dlnActCoeffds) const
835 {
836  int delAK, delBK;
837  double XA, XB, g0, g1, dXA, dXB;
838  double T = temperature();
839  double RT = GasConstant*T;
840 
841  //fvo_zero_dbl_1(dlnActCoeff, m_kk);
843 
844  for (size_t iK = 0; iK < m_kk; iK++) {
845  dlnActCoeffds[iK] = 0.0;
846  for (size_t i = 0; i < numBinaryInteractions_; i++) {
847 
848  size_t iA = m_pSpecies_A_ij[i];
849  size_t iB = m_pSpecies_B_ij[i];
850 
851  delAK = 0;
852  delBK = 0;
853 
854  if (iA==iK) {
855  delAK = 1;
856  } else if (iB==iK) {
857  delBK = 1;
858  }
859 
860  XA = moleFractions_[iA];
861  XB = moleFractions_[iB];
862 
863  dXA = dXds[iA];
864  dXB = dXds[iB];
865 
866  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
867  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
868 
869  dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB
870  + dlnActCoeffdT_Scaled_[iK]*dTds;
871  }
872  }
873 }
874 //====================================================================================================================
875 // Update the derivative of the log of the activity coefficients wrt dlnN
876 /*
877  * This function will be called to update the internally stored gradients of the
878  * logarithm of the activity coefficients. These are used in the determination
879  * of the diffusion coefficients.
880  *
881  * he = X_A X_B(B + C X_B)
882  */
884 {
885  int delAK, delBK;
886  double XA, XB, XK, g0, g1;
887  double T = temperature();
888  double RT = GasConstant*T;
889 
890  dlnActCoeffdlnN_diag_.assign(m_kk, 0);
891 
892  for (size_t iK = 0; iK < m_kk; iK++) {
893 
894  XK = moleFractions_[iK];
895 
896  for (size_t i = 0; i < numBinaryInteractions_; i++) {
897 
898  size_t iA = m_pSpecies_A_ij[i];
899  size_t iB = m_pSpecies_B_ij[i];
900 
901  delAK = 0;
902  delBK = 0;
903 
904  if (iA==iK) {
905  delAK = 1;
906  } else if (iB==iK) {
907  delBK = 1;
908  }
909 
910  XA = moleFractions_[iA];
911  XB = moleFractions_[iB];
912 
913  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
914  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
915 
916  dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB)));
917 
918  // double gfac = g0 + g1 * XB;
919  // double gggg = (delBK - XB) * g1;
920 
921 
922  // dlnActCoeffdlnN_diag_[iK] += gfac * delAK * ( - XB + delBK);
923 
924  // dlnActCoeffdlnN_diag_[iK] += gfac * delBK * ( - XA + delAK);
925 
926  // dlnActCoeffdlnN_diag_[iK] += gfac * (2.0 * XA * XB - delAK * XB - XA * delBK);
927 
928  // dlnActCoeffdlnN_diag_[iK] += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBK);
929 
930  // dlnActCoeffdlnN_diag_[iK] += gggg * ( - 2.0 * XA * XB + delAK * XB + XA * delBK);
931 
932  // dlnActCoeffdlnN_diag_[iK] += - g1 * XA * XB * (- XB + delBK);
933  }
934  dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK];//-XK;
935  }
936 }
937 
938 //====================================================================================================================
939 // Update the derivative of the log of the activity coefficients wrt dlnN
940 /*
941  * This function will be called to update the internally stored gradients of the
942  * logarithm of the activity coefficients. These are used in the determination
943  * of the diffusion coefficients.
944  *
945  */
947 {
948  doublereal delAK, delBK;
949  double XA, XB, g0, g1,XM;
950  double T = temperature();
951  double RT = GasConstant*T;
952 
953  doublereal delAM, delBM;
955 
956  /*
957  * Loop over the activity coefficient gamma_k
958  */
959  for (size_t iK = 0; iK < m_kk; iK++) {
960  for (size_t iM = 0; iM < m_kk; iM++) {
961  XM = moleFractions_[iM];
962  for (size_t i = 0; i < numBinaryInteractions_; i++) {
963 
964  size_t iA = m_pSpecies_A_ij[i];
965  size_t iB = m_pSpecies_B_ij[i];
966 
967  delAK = 0.0;
968  delBK = 0.0;
969  delAM = 0.0;
970  delBM = 0.0;
971  if (iA==iK) {
972  delAK = 1.0;
973  } else if (iB==iK) {
974  delBK = 1.0;
975  }
976  if (iA==iM) {
977  delAM = 1.0;
978  } else if (iB==iM) {
979  delBM = 1.0;
980  }
981 
982  XA = moleFractions_[iA];
983  XB = moleFractions_[iB];
984 
985  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
986  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
987 
988  dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB));
989  dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA);
990 
991  // double gfac = g0 + g1 * XB;
992  // double gggg = (delBK - XB) * g1;
993 
994 
995  // dlnActCoeffdlnN_(iK, iM) += gfac * delAK * ( - XB + delBM);
996 
997  // dlnActCoeffdlnN_(iK, iM) += gfac * delBK * ( - XA + delAM);
998 
999  // dlnActCoeffdlnN_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM);
1000 
1001  // dlnActCoeffdlnN_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM);
1002 
1003  // dlnActCoeffdlnN_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM);
1004 
1005  // dlnActCoeffdlnN_(iK, iM) += - g1 * XA * XB * (- XB + delBM);
1006  }
1007  dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM);
1008  }
1009  }
1010 }
1011 //====================================================================================================================
1013 {
1014  doublereal XA, XB, g0 , g1;
1015  doublereal T = temperature();
1016 
1017  dlnActCoeffdlnX_diag_.assign(m_kk, 0);
1018  doublereal RT = GasConstant * T;
1019 
1020  for (size_t i = 0; i < numBinaryInteractions_; i++) {
1021 
1022  size_t iA = m_pSpecies_A_ij[i];
1023  size_t iB = m_pSpecies_B_ij[i];
1024 
1025  XA = moleFractions_[iA];
1026  XB = moleFractions_[iB];
1027 
1028  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
1029  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
1030 
1031  dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB);
1032  dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB);
1033  }
1034 }
1035 
1036 //====================================================================================================================
1037 void MixedSolventElectrolyte::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
1038 {
1040  for (size_t k = 0; k < m_kk; k++) {
1041  dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
1042  }
1043 }
1044 //====================================================================================================================
1045 void MixedSolventElectrolyte::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
1046 {
1048  for (size_t k = 0; k < m_kk; k++) {
1049  dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
1050  }
1051 }
1052 //====================================================================================================================
1053 void MixedSolventElectrolyte::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
1054 {
1056  double* data = & dlnActCoeffdlnN_(0,0);
1057  for (size_t k = 0; k < m_kk; k++) {
1058  for (size_t m = 0; m < m_kk; m++) {
1059  dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
1060  }
1061  }
1062 }
1063 //====================================================================================================================
1065 {
1066  numBinaryInteractions_ = num;
1067  m_HE_b_ij.resize(num, 0.0);
1068  m_HE_c_ij.resize(num, 0.0);
1069  m_HE_d_ij.resize(num, 0.0);
1070  m_SE_b_ij.resize(num, 0.0);
1071  m_SE_c_ij.resize(num, 0.0);
1072  m_SE_d_ij.resize(num, 0.0);
1073  m_VHE_b_ij.resize(num, 0.0);
1074  m_VHE_c_ij.resize(num, 0.0);
1075  m_VHE_d_ij.resize(num, 0.0);
1076  m_VSE_b_ij.resize(num, 0.0);
1077  m_VSE_c_ij.resize(num, 0.0);
1078  m_VSE_d_ij.resize(num, 0.0);
1079 
1080  m_pSpecies_A_ij.resize(num, npos);
1081  m_pSpecies_B_ij.resize(num, npos);
1082 
1083 }
1084 //====================================================================================================================
1085 
1086 /*
1087  * Process an XML node called "binaryNeutralSpeciesParameters"
1088  * This node contains all of the parameters necessary to describe
1089  * the Margules Interaction for a single binary interaction
1090  * This function reads the XML file and writes the coefficients
1091  * it finds to an internal data structures.
1092  */
1094 {
1095  string xname = xmLBinarySpecies.name();
1096  if (xname != "binaryNeutralSpeciesParameters") {
1097  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies",
1098  "Incorrect name for processing this routine: " + xname);
1099  }
1100  double* charge = DATA_PTR(m_speciesCharge);
1101  string stemp;
1102  size_t nParamsFound;
1103  vector_fp vParams;
1104  string iName = xmLBinarySpecies.attrib("speciesA");
1105  if (iName == "") {
1106  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "no speciesA attrib");
1107  }
1108  string jName = xmLBinarySpecies.attrib("speciesB");
1109  if (jName == "") {
1110  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "no speciesB attrib");
1111  }
1112  /*
1113  * Find the index of the species in the current phase. It's not
1114  * an error to not find the species
1115  */
1116  size_t iSpecies = speciesIndex(iName);
1117  if (iSpecies == npos) {
1118  return;
1119  }
1120  string ispName = speciesName(iSpecies);
1121  if (charge[iSpecies] != 0) {
1122  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "speciesA charge problem");
1123  }
1124  size_t jSpecies = speciesIndex(jName);
1125  if (jSpecies == npos) {
1126  return;
1127  }
1128  string jspName = speciesName(jSpecies);
1129  if (charge[jSpecies] != 0) {
1130  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies", "speciesB charge problem");
1131  }
1132 
1134  size_t iSpot = numBinaryInteractions_ - 1;
1135  m_pSpecies_A_ij[iSpot] = iSpecies;
1136  m_pSpecies_B_ij[iSpot] = jSpecies;
1137 
1138  size_t num = xmLBinarySpecies.nChildren();
1139  for (size_t iChild = 0; iChild < num; iChild++) {
1140  XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
1141  stemp = xmlChild.name();
1142  string nodeName = lowercase(stemp);
1143  /*
1144  * Process the binary species interaction child elements
1145  */
1146  if (nodeName == "excessenthalpy") {
1147  /*
1148  * Get the string containing all of the values
1149  */
1150  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
1151  nParamsFound = vParams.size();
1152 
1153  if (nParamsFound != 2) {
1154  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEnthalpy for " + ispName
1155  + "::" + jspName,
1156  "wrong number of params found");
1157  }
1158  m_HE_b_ij[iSpot] = vParams[0];
1159  m_HE_c_ij[iSpot] = vParams[1];
1160  }
1161 
1162  if (nodeName == "excessentropy") {
1163  /*
1164  * Get the string containing all of the values
1165  */
1166  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
1167  nParamsFound = vParams.size();
1168 
1169  if (nParamsFound != 2) {
1170  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessEntropy for " + ispName
1171  + "::" + jspName,
1172  "wrong number of params found");
1173  }
1174  m_SE_b_ij[iSpot] = vParams[0];
1175  m_SE_c_ij[iSpot] = vParams[1];
1176  }
1177 
1178  if (nodeName == "excessvolume_enthalpy") {
1179  /*
1180  * Get the string containing all of the values
1181  */
1182  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
1183  nParamsFound = vParams.size();
1184 
1185  if (nParamsFound != 2) {
1186  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
1187  + "::" + jspName,
1188  "wrong number of params found");
1189  }
1190  m_VHE_b_ij[iSpot] = vParams[0];
1191  m_VHE_c_ij[iSpot] = vParams[1];
1192  }
1193 
1194  if (nodeName == "excessvolume_entropy") {
1195  /*
1196  * Get the string containing all of the values
1197  */
1198  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
1199  nParamsFound = vParams.size();
1200 
1201  if (nParamsFound != 2) {
1202  throw CanteraError("MixedSolventElectrolyte::readXMLBinarySpecies::excessVolume_Entropy for " + ispName
1203  + "::" + jspName,
1204  "wrong number of params found");
1205  }
1206  m_VSE_b_ij[iSpot] = vParams[0];
1207  m_VSE_c_ij[iSpot] = vParams[1];
1208  }
1209 
1210 
1211  }
1212 
1213 }
1214 
1215 }
1216