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