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