Cantera  2.1.2
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  * Copyright (2009) Sandia Corporation. Under the terms of
10  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
11  * U.S. Government retains certain rights in this software.
12  */
16 
17 #include <iomanip>
18 #include <fstream>
19 
20 using namespace std;
21 
22 namespace Cantera
23 {
24 MargulesVPSSTP::MargulesVPSSTP() :
26  numBinaryInteractions_(0),
27  formMargules_(0),
28  formTempModel_(0)
29 {
30 }
31 
32 MargulesVPSSTP::MargulesVPSSTP(const std::string& inputFile, const std::string& id_) :
34  numBinaryInteractions_(0),
35  formMargules_(0),
36  formTempModel_(0)
37 {
38  initThermoFile(inputFile, id_);
39 }
40 
41 MargulesVPSSTP::MargulesVPSSTP(XML_Node& phaseRoot, const std::string& id_) :
43  numBinaryInteractions_(0),
44  formMargules_(0),
45  formTempModel_(0)
46 {
47  importPhase(*findXMLPhase(&phaseRoot, id_), this);
48 }
49 
52 {
54 }
55 
58 {
59  if (&b == this) {
60  return *this;
61  }
62 
64 
66  m_HE_b_ij = b.m_HE_b_ij;
67  m_HE_c_ij = b.m_HE_c_ij;
68  m_HE_d_ij = b.m_HE_d_ij;
69  m_SE_b_ij = b.m_SE_b_ij;
70  m_SE_c_ij = b.m_SE_c_ij;
71  m_SE_d_ij = b.m_SE_d_ij;
82 
83  return *this;
84 }
85 
88 {
89  return new MargulesVPSSTP(*this);
90 }
91 
94  numBinaryInteractions_(0),
95  formMargules_(0),
96  formTempModel_(0)
97 {
98 
99  initThermoFile("LiKCl_liquid.xml", "");
100 
101 
103 
104  m_HE_b_ij.resize(1);
105  m_HE_c_ij.resize(1);
106  m_HE_d_ij.resize(1);
107 
108  m_SE_b_ij.resize(1);
109  m_SE_c_ij.resize(1);
110  m_SE_d_ij.resize(1);
111 
112  m_VHE_b_ij.resize(1);
113  m_VHE_c_ij.resize(1);
114  m_VHE_d_ij.resize(1);
115 
116  m_VSE_b_ij.resize(1);
117  m_VSE_c_ij.resize(1);
118  m_VSE_d_ij.resize(1);
119 
120  m_pSpecies_A_ij.resize(1);
121  m_pSpecies_B_ij.resize(1);
122 
123 
124 
125  m_HE_b_ij[0] = -17570E3;
126  m_HE_c_ij[0] = -377.0E3;
127  m_HE_d_ij[0] = 0.0;
128 
129  m_SE_b_ij[0] = -7.627E3;
130  m_SE_c_ij[0] = 4.958E3;
131  m_SE_d_ij[0] = 0.0;
132 
133 
134  size_t iLiCl = speciesIndex("LiCl(L)");
135  if (iLiCl == npos) {
136  throw CanteraError("MargulesVPSSTP test1 constructor",
137  "Unable to find LiCl(L)");
138  }
139  m_pSpecies_B_ij[0] = iLiCl;
140 
141 
142  size_t iKCl = speciesIndex("KCl(L)");
143  if (iKCl == npos) {
144  throw CanteraError("MargulesVPSSTP test1 constructor",
145  "Unable to find KCl(L)");
146  }
147  m_pSpecies_A_ij[0] = iKCl;
148 }
149 
150 /*
151  * -------------- Utilities -------------------------------
152  */
153 
155 {
156  return 0;
157 }
158 
159 /*
160  * - Activities, Standard States, Activity Concentrations -----------
161  */
162 
163 void MargulesVPSSTP::getLnActivityCoefficients(doublereal* lnac) const
164 {
165  /*
166  * Update the activity coefficients
167  */
169 
170  /*
171  * take the exp of the internally stored coefficients.
172  */
173  for (size_t k = 0; k < m_kk; k++) {
174  lnac[k] = lnActCoeff_Scaled_[k];
175  }
176 }
177 
178 /*
179  * ------------ Partial Molar Properties of the Solution ------------
180  */
181 
183 {
184  getChemPotentials(mu);
185  double ve = Faraday * electricPotential();
186  for (size_t k = 0; k < m_kk; k++) {
187  mu[k] += ve*charge(k);
188  }
189 }
190 
191 void MargulesVPSSTP::getChemPotentials(doublereal* mu) const
192 {
193  doublereal xx;
194  /*
195  * First get the standard chemical potentials in
196  * molar form.
197  * -> this requires updates of standard state as a function
198  * of T and P
199  */
201  /*
202  * Update the activity coefficients
203  */
205  doublereal RT = GasConstant * temperature();
206  for (size_t k = 0; k < m_kk; k++) {
207  xx = std::max(moleFractions_[k], SmallNumber);
208  mu[k] += RT * (log(xx) + lnActCoeff_Scaled_[k]);
209  }
210 }
211 
213 {
214  size_t kk = nSpecies();
215  double h = 0;
216  vector_fp hbar(kk);
217  getPartialMolarEnthalpies(&hbar[0]);
218  for (size_t i = 0; i < kk; i++) {
219  h += moleFractions_[i]*hbar[i];
220  }
221  return h;
222 }
223 
225 {
226  size_t kk = nSpecies();
227  double s = 0;
228  vector_fp sbar(kk);
229  getPartialMolarEntropies(&sbar[0]);
230  for (size_t i = 0; i < kk; i++) {
231  s += moleFractions_[i]*sbar[i];
232  }
233  return s;
234 }
235 
236 doublereal MargulesVPSSTP::cp_mole() const
237 {
238  size_t kk = nSpecies();
239  double cp = 0;
240  vector_fp cpbar(kk);
241  getPartialMolarCp(&cpbar[0]);
242  for (size_t i = 0; i < kk; i++) {
243  cp += moleFractions_[i]*cpbar[i];
244  }
245  return cp;
246 }
247 
248 doublereal MargulesVPSSTP::cv_mole() const
249 {
250  return cp_mole() - GasConstant;
251 }
252 
253 void MargulesVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
254 {
255  /*
256  * Get the nondimensional standard state enthalpies
257  */
258  getEnthalpy_RT(hbar);
259  /*
260  * dimensionalize it.
261  */
262  double T = temperature();
263  double RT = GasConstant * T;
264  for (size_t k = 0; k < m_kk; k++) {
265  hbar[k] *= RT;
266  }
267  /*
268  * Update the activity coefficients, This also update the
269  * internally stored molalities.
270  */
273  double RTT = RT * T;
274  for (size_t k = 0; k < m_kk; k++) {
275  hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
276  }
277 }
278 
279 void MargulesVPSSTP::getPartialMolarCp(doublereal* cpbar) const
280 {
281  /*
282  * Get the nondimensional standard state entropies
283  */
284  getCp_R(cpbar);
285  double T = temperature();
286  /*
287  * Update the activity coefficients, This also update the
288  * internally stored molalities.
289  */
292 
293  for (size_t k = 0; k < m_kk; k++) {
294  cpbar[k] -= 2 * T * dlnActCoeffdT_Scaled_[k] + T * T * d2lnActCoeffdT2_Scaled_[k];
295  }
296  /*
297  * dimensionalize it.
298  */
299  for (size_t k = 0; k < m_kk; k++) {
300  cpbar[k] *= GasConstant;
301  }
302 }
303 
304 void MargulesVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
305 {
306  double xx;
307  /*
308  * Get the nondimensional standard state entropies
309  */
310  getEntropy_R(sbar);
311  double T = temperature();
312  /*
313  * Update the activity coefficients, This also update the
314  * internally stored molalities.
315  */
318 
319  for (size_t k = 0; k < m_kk; k++) {
320  xx = std::max(moleFractions_[k], SmallNumber);
321  sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
322  }
323  /*
324  * dimensionalize it.
325  */
326  for (size_t k = 0; k < m_kk; k++) {
327  sbar[k] *= GasConstant;
328  }
329 }
330 
331 void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const
332 {
333 
334  size_t iA, iB, delAK, delBK;
335  double XA, XB, g0 , g1;
336  double T = temperature();
337 
338  /*
339  * Get the standard state values in m^3 kmol-1
340  */
341  getStandardVolumes(vbar);
342 
343  for (size_t iK = 0; iK < m_kk; iK++) {
344  delAK = 0;
345  delBK = 0;
346  }
347  for (size_t i = 0; i < numBinaryInteractions_; i++) {
348  iA = m_pSpecies_A_ij[i];
349  iB = m_pSpecies_B_ij[i];
350  XA = moleFractions_[iA];
351  XB = moleFractions_[iB];
352  g0 = (m_VHE_b_ij[i] - T * m_VSE_b_ij[i]);
353  g1 = (m_VHE_c_ij[i] - T * m_VSE_c_ij[i]);
354  const doublereal temp1 = g0 + g1 * XB;
355  const doublereal all = -1.0*XA*XB*temp1 - XA*XB*XB*g1;
356 
357  for (size_t iK = 0; iK < m_kk; iK++) {
358  vbar[iK] += all;
359  // vbar[iK] += XA*XB*temp1+((delAK-XA)*XB+XA*(delBK-XB))*temp1+XB*XA*(delBK-XB)*g1;
360  }
361  vbar[iA] += XB * temp1;
362  vbar[iB] += XA * temp1 + XA*XB*g1;
363  }
364 }
365 
366 doublereal MargulesVPSSTP::err(const std::string& msg) const
367 {
368  throw CanteraError("MargulesVPSSTP","Base class method "
369  +msg+" called. Equation of state type: "+int2str(eosType()));
370  return 0;
371 }
372 
374 {
375  initLengths();
377 }
378 
380 {
381  m_kk = nSpecies();
383 }
384 
385 void MargulesVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
386 {
387  string stemp;
388  string subname = "MargulesVPSSTP::initThermoXML";
389  if ((int) id_.size() > 0) {
390  string idp = phaseNode.id();
391  if (idp != id_) {
392  throw CanteraError(subname, "phasenode and Id are incompatible");
393  }
394  }
395 
396  /*
397  * Find the Thermo XML node
398  */
399  if (!phaseNode.hasChild("thermo")) {
400  throw CanteraError(subname,
401  "no thermo XML node");
402  }
403  XML_Node& thermoNode = phaseNode.child("thermo");
404 
405  /*
406  * Make sure that the thermo model is Margules
407  */
408  stemp = thermoNode.attrib("model");
409  string formString = lowercase(stemp);
410  if (formString != "margules") {
411  throw CanteraError(subname,
412  "model name isn't Margules: " + formString);
413 
414  }
415 
416  /*
417  * Go get all of the coefficients and factors in the
418  * activityCoefficients XML block
419  */
420  XML_Node* acNodePtr = 0;
421  if (thermoNode.hasChild("activityCoefficients")) {
422  XML_Node& acNode = thermoNode.child("activityCoefficients");
423  acNodePtr = &acNode;
424  string mStringa = acNode.attrib("model");
425  string mString = lowercase(mStringa);
426  if (mString != "margules") {
427  throw CanteraError(subname.c_str(),
428  "Unknown activity coefficient model: " + mStringa);
429  }
430  for (size_t i = 0; i < acNodePtr->nChildren(); i++) {
431  XML_Node& xmlACChild = acNodePtr->child(i);
432  stemp = xmlACChild.name();
433  string nodeName = lowercase(stemp);
434  /*
435  * Process a binary salt field, or any of the other XML fields
436  * that make up the Pitzer Database. Entries will be ignored
437  * if any of the species in the entry isn't in the solution.
438  */
439  if (nodeName == "binaryneutralspeciesparameters") {
440  readXMLBinarySpecies(xmlACChild);
441 
442  }
443  }
444  }
445 
446  /*
447  * Go down the chain
448  */
449  GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
450 
451 
452 }
453 
455 {
456  size_t iA, iB, iK;
457  double XA, XB, g0 , g1;
458  double T = temperature();
459  double invRT = 1.0 / (GasConstant*T);
460  lnActCoeff_Scaled_.resize(m_kk);
461  for (iK = 0; iK < m_kk; iK++) {
462  lnActCoeff_Scaled_[iK] = 0.0;
463  }
464  for (size_t i = 0; i < numBinaryInteractions_; i++) {
465  iA = m_pSpecies_A_ij[i];
466  iB = m_pSpecies_B_ij[i];
467  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) * invRT;
468  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) * invRT;
469  XA = moleFractions_[iA];
470  XB = moleFractions_[iB];
471  const doublereal XAXB = XA * XB;
472  const doublereal g0g1XB = (g0 + g1 * XB);
473  const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
474  for (iK = 0; iK < m_kk; iK++) {
475  lnActCoeff_Scaled_[iK] += all;
476  }
477  lnActCoeff_Scaled_[iA] += XB * g0g1XB;
478  lnActCoeff_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
479  }
480 }
481 
483 {
484  size_t iA, iB, iK;
485  doublereal XA, XB, g0, g1;
486  doublereal invT = 1.0 / temperature();
487  doublereal invRTT = 1.0 / (GasConstant)*invT*invT;
488  dlnActCoeffdT_Scaled_.resize(m_kk);
490  for (iK = 0; iK < m_kk; iK++) {
491  dlnActCoeffdT_Scaled_[iK] = 0.0;
492  d2lnActCoeffdT2_Scaled_[iK] = 0.0;
493  }
494  for (size_t i = 0; i < numBinaryInteractions_; i++) {
495  iA = m_pSpecies_A_ij[i];
496  iB = m_pSpecies_B_ij[i];
497  XA = moleFractions_[iA];
498  XB = moleFractions_[iB];
499  g0 = -m_HE_b_ij[i] * invRTT;
500  g1 = -m_HE_c_ij[i] * invRTT;
501  const doublereal XAXB = XA * XB;
502  const doublereal g0g1XB = (g0 + g1 * XB);
503  const doublereal all = -1.0 * XAXB * g0g1XB - XAXB * XB * g1;
504  const doublereal mult = 2.0 * invT;
505  const doublereal dT2all = mult * all;
506  for (iK = 0; iK < m_kk; iK++) {
507  // double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1;
508  dlnActCoeffdT_Scaled_[iK] += all;
509  d2lnActCoeffdT2_Scaled_[iK] -= dT2all;
510  }
511  dlnActCoeffdT_Scaled_[iA] += XB * g0g1XB;
512  dlnActCoeffdT_Scaled_[iB] += XA * g0g1XB + XAXB * g1;
513  d2lnActCoeffdT2_Scaled_[iA] -= mult * XB * g0g1XB;
514  d2lnActCoeffdT2_Scaled_[iB] -= mult * (XA * g0g1XB + XAXB * g1);
515  }
516 }
517 
518 void MargulesVPSSTP::getdlnActCoeffdT(doublereal* dlnActCoeffdT) const
519 {
521  for (size_t k = 0; k < m_kk; k++) {
522  dlnActCoeffdT[k] = dlnActCoeffdT_Scaled_[k];
523  }
524 }
525 
526 void MargulesVPSSTP::getd2lnActCoeffdT2(doublereal* d2lnActCoeffdT2) const
527 {
529  for (size_t k = 0; k < m_kk; k++) {
530  d2lnActCoeffdT2[k] = d2lnActCoeffdT2_Scaled_[k];
531  }
532 }
533 
534 void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
535  doublereal* dlnActCoeffds) const
536 {
537  size_t iA, iB, iK;
538  double XA, XB, g0 , g1, dXA, dXB;
539  double T = temperature();
540  double RT = GasConstant*T;
541 
542  //fvo_zero_dbl_1(dlnActCoeff, m_kk);
544 
545  for (iK = 0; iK < m_kk; iK++) {
546  dlnActCoeffds[iK] = 0.0;
547  }
548 
549  for (size_t i = 0; i < numBinaryInteractions_; i++) {
550  iA = m_pSpecies_A_ij[i];
551  iB = m_pSpecies_B_ij[i];
552  XA = moleFractions_[iA];
553  XB = moleFractions_[iB];
554  dXA = dXds[iA];
555  dXB = dXds[iB];
556  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
557  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
558  const doublereal g02g1XB = g0 + 2*g1*XB;
559  const doublereal g2XAdXB = 2*g1*XA*dXB;
560  const doublereal all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB;
561  for (iK = 0; iK < m_kk; iK++) {
562  // dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB
563  // + dlnActCoeffdT_Scaled_[iK]*dTds;
564  dlnActCoeffds[iK] += all + dlnActCoeffdT_Scaled_[iK]*dTds;
565  }
566  dlnActCoeffds[iA] += dXB * g02g1XB;
567  dlnActCoeffds[iB] += dXA * g02g1XB + g2XAdXB;
568  }
569 }
570 
572 {
573  size_t iA, iB, iK, delAK, delBK;
574  double XA, XB, XK, g0 , g1;
575  double T = temperature();
576  double RT = GasConstant*T;
577 
578  dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
579 
580  for (iK = 0; iK < m_kk; iK++) {
581 
582  XK = moleFractions_[iK];
583 
584  for (size_t i = 0; i < numBinaryInteractions_; i++) {
585 
586  iA = m_pSpecies_A_ij[i];
587  iB = m_pSpecies_B_ij[i];
588 
589  delAK = 0;
590  delBK = 0;
591 
592  if (iA==iK) {
593  delAK = 1;
594  } else if (iB==iK) {
595  delBK = 1;
596  }
597 
598  XA = moleFractions_[iA];
599  XB = moleFractions_[iB];
600 
601  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
602  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
603 
604  dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB)));
605 
606  // double gfac = g0 + g1 * XB;
607  // double gggg = (delBK - XB) * g1;
608 
609 
610  // dlnActCoeffdlnN_diag_[iK] += gfac * delAK * ( - XB + delBK);
611 
612  // dlnActCoeffdlnN_diag_[iK] += gfac * delBK * ( - XA + delAK);
613 
614  // dlnActCoeffdlnN_diag_[iK] += gfac * (2.0 * XA * XB - delAK * XB - XA * delBK);
615 
616  // dlnActCoeffdlnN_diag_[iK] += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBK);
617 
618  // dlnActCoeffdlnN_diag_[iK] += gggg * ( - 2.0 * XA * XB + delAK * XB + XA * delBK);
619 
620  // dlnActCoeffdlnN_diag_[iK] += - g1 * XA * XB * (- XB + delBK);
621  }
622  dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK];//-XK;
623  }
624 }
625 
627 {
628  size_t iA, iB;
629  doublereal delAK, delBK;
630  double XA, XB, g0, g1,XM;
631  double T = temperature();
632  double RT = GasConstant*T;
633 
634  doublereal delAM, delBM;
635 
637 
638  /*
639  * Loop over the activity coefficient gamma_k
640  */
641  for (size_t iK = 0; iK < m_kk; iK++) {
642  for (size_t iM = 0; iM < m_kk; iM++) {
643  XM = moleFractions_[iM];
644  for (size_t i = 0; i < numBinaryInteractions_; i++) {
645 
646  iA = m_pSpecies_A_ij[i];
647  iB = m_pSpecies_B_ij[i];
648 
649  delAK = 0.0;
650  delBK = 0.0;
651  delAM = 0.0;
652  delBM = 0.0;
653  if (iA==iK) {
654  delAK = 1.0;
655  } else if (iB==iK) {
656  delBK = 1.0;
657  }
658  if (iA==iM) {
659  delAM = 1.0;
660  } else if (iB==iM) {
661  delBM = 1.0;
662  }
663 
664  XA = moleFractions_[iA];
665  XB = moleFractions_[iB];
666 
667  g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
668  g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
669 
670  dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB));
671  dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA);
672 
673  // double gfac = g0 + g1 * XB;
674  // double gggg = (delBK - XB) * g1;
675 
676 
677  // dlnActCoeffdlnN_(iK, iM) += gfac * delAK * ( - XB + delBM);
678 
679  // dlnActCoeffdlnN_(iK, iM) += gfac * delBK * ( - XA + delAM);
680 
681  // dlnActCoeffdlnN_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM);
682 
683  // dlnActCoeffdlnN_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM);
684 
685  // dlnActCoeffdlnN_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM);
686 
687  // dlnActCoeffdlnN_(iK, iM) += - g1 * XA * XB * (- XB + delBM);
688  }
689  dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM);
690  }
691  }
692 }
693 
695 {
696  doublereal T = temperature();
697  dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
698  doublereal RT = GasConstant * T;
699 
700  for (size_t i = 0; i < numBinaryInteractions_; i++) {
701  size_t iA = m_pSpecies_A_ij[i];
702  size_t iB = m_pSpecies_B_ij[i];
703 
704  doublereal XA = moleFractions_[iA];
705  doublereal XB = moleFractions_[iB];
706 
707  doublereal g0 = (m_HE_b_ij[i] - T * m_SE_b_ij[i]) / RT;
708  doublereal g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT;
709 
710  dlnActCoeffdlnX_diag_[iA] += XA*XB*(2*g1*-2*g0-6*g1*XB);
711  dlnActCoeffdlnX_diag_[iB] += XA*XB*(2*g1*-2*g0-6*g1*XB);
712  }
713 }
714 
715 void MargulesVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
716 {
718  for (size_t k = 0; k < m_kk; k++) {
719  dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
720  }
721 }
722 
723 void MargulesVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
724 {
726  for (size_t k = 0; k < m_kk; k++) {
727  dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
728  }
729 }
730 
731 void MargulesVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
732 {
734  double* data = & dlnActCoeffdlnN_(0,0);
735  for (size_t k = 0; k < m_kk; k++) {
736  for (size_t m = 0; m < m_kk; m++) {
737  dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
738  }
739  }
740 }
741 
743 {
745  m_HE_b_ij.resize(num, 0.0);
746  m_HE_c_ij.resize(num, 0.0);
747  m_HE_d_ij.resize(num, 0.0);
748  m_SE_b_ij.resize(num, 0.0);
749  m_SE_c_ij.resize(num, 0.0);
750  m_SE_d_ij.resize(num, 0.0);
751  m_VHE_b_ij.resize(num, 0.0);
752  m_VHE_c_ij.resize(num, 0.0);
753  m_VHE_d_ij.resize(num, 0.0);
754  m_VSE_b_ij.resize(num, 0.0);
755  m_VSE_c_ij.resize(num, 0.0);
756  m_VSE_d_ij.resize(num, 0.0);
757 
758  m_pSpecies_A_ij.resize(num, npos);
759  m_pSpecies_B_ij.resize(num, npos);
760 }
761 
763 {
764  string xname = xmLBinarySpecies.name();
765  if (xname != "binaryNeutralSpeciesParameters") {
766  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies",
767  "Incorrect name for processing this routine: " + xname);
768  }
769  string stemp;
770  size_t nParamsFound;
771  vector_fp vParams;
772  string iName = xmLBinarySpecies.attrib("speciesA");
773  if (iName == "") {
774  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesA attrib");
775  }
776  string jName = xmLBinarySpecies.attrib("speciesB");
777  if (jName == "") {
778  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "no speciesB attrib");
779  }
780  /*
781  * Find the index of the species in the current phase. It's not
782  * an error to not find the species
783  */
784  size_t iSpecies = speciesIndex(iName);
785  if (iSpecies == npos) {
786  return;
787  }
788  string ispName = speciesName(iSpecies);
789  if (charge(iSpecies) != 0) {
790  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesA charge problem");
791  }
792  size_t jSpecies = speciesIndex(jName);
793  if (jSpecies == npos) {
794  return;
795  }
796  string jspName = speciesName(jSpecies);
797  if (charge(jSpecies) != 0) {
798  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies", "speciesB charge problem");
799  }
800 
802  size_t iSpot = numBinaryInteractions_ - 1;
803  m_pSpecies_A_ij[iSpot] = iSpecies;
804  m_pSpecies_B_ij[iSpot] = jSpecies;
805 
806  for (size_t iChild = 0; iChild < xmLBinarySpecies.nChildren(); iChild++) {
807  XML_Node& xmlChild = xmLBinarySpecies.child(iChild);
808  stemp = xmlChild.name();
809  string nodeName = lowercase(stemp);
810  /*
811  * Process the binary species interaction child elements
812  */
813  if (nodeName == "excessenthalpy") {
814  /*
815  * Get the string containing all of the values
816  */
817  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEnthalpy");
818  nParamsFound = vParams.size();
819 
820  if (nParamsFound != 2) {
821  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEnthalpy for " + ispName
822  + "::" + jspName,
823  "wrong number of params found. Need 2");
824  }
825  m_HE_b_ij[iSpot] = vParams[0];
826  m_HE_c_ij[iSpot] = vParams[1];
827  }
828 
829  if (nodeName == "excessentropy") {
830  /*
831  * Get the string containing all of the values
832  */
833  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessEntropy");
834  nParamsFound = vParams.size();
835 
836  if (nParamsFound != 2) {
837  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessEntropy for " + ispName
838  + "::" + jspName,
839  "wrong number of params found. Need 2");
840  }
841  m_SE_b_ij[iSpot] = vParams[0];
842  m_SE_c_ij[iSpot] = vParams[1];
843  }
844 
845  if (nodeName == "excessvolume_enthalpy") {
846  /*
847  * Get the string containing all of the values
848  */
849  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Enthalpy");
850  nParamsFound = vParams.size();
851 
852  if (nParamsFound != 2) {
853  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Enthalpy for " + ispName
854  + "::" + jspName,
855  "wrong number of params found. Need 2");
856  }
857  m_VHE_b_ij[iSpot] = vParams[0];
858  m_VHE_c_ij[iSpot] = vParams[1];
859  }
860 
861  if (nodeName == "excessvolume_entropy") {
862  /*
863  * Get the string containing all of the values
864  */
865  ctml::getFloatArray(xmlChild, vParams, true, "toSI", "excessVolume_Entropy");
866  nParamsFound = vParams.size();
867 
868  if (nParamsFound != 2) {
869  throw CanteraError("MargulesVPSSTP::readXMLBinarySpecies::excessVolume_Entropy for " + ispName
870  + "::" + jspName,
871  "wrong number of params found. Need 2");
872  }
873  m_VSE_b_ij[iSpot] = vParams[0];
874  m_VSE_c_ij[iSpot] = vParams[1];
875  }
876  }
877 }
878 
879 }
int formMargules_
form of the Margules interaction expression
vector_fp m_VSE_c_ij
Entropy term for the ternary mole fraction interaction of the excess gibbs free energy expression...
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
void readXMLBinarySpecies(XML_Node &xmlBinarySpecies)
Process an XML node called "binaryNeutralSpeciesParameters".
doublereal electricPotential() const
Returns the electric potential of this phase (V).
Definition: ThermoPhase.h:391
vector_fp m_VSE_d_ij
Entropy term for the quaternary mole fraction interaction of the excess gibbs free energy expression...
vector_fp m_HE_b_ij
Enthalpy term for the binary mole fraction interaction of the excess gibbs free energy expression...
XML_Node * findXMLPhase(XML_Node *root, const std::string &idtarget)
Search an XML_Node tree for a named phase XML_Node.
Definition: xml.cpp:1104
vector_fp m_VHE_d_ij
Enthalpy term for the quaternary mole fraction interaction of the excess gibbs free energy expression...
Header for intermediate ThermoPhase object for phases which employ gibbs excess free energy based for...
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:534
GibbsExcessVPSSTP & operator=(const GibbsExcessVPSSTP &b)
Assignment operator.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.h:128
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
std::vector< size_t > m_pSpecies_A_ij
vector of species indices representing species A in the interaction
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
vector_fp m_HE_d_ij
Enthalpy term for the quaternary mole fraction interaction of the excess gibbs free energy expression...
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
virtual void getd2lnActCoeffdT2(doublereal *d2lnActCoeffdT2) const
Get the array of temperature second derivatives of the log activity coefficients. ...
size_t numBinaryInteractions_
number of binary interaction expressions
std::string lowercase(const std::string &s)
Cast a copy of a string to lower case.
Definition: stringUtils.cpp:58
virtual void getdlnActCoeffdlnN(const size_t ld, doublereal *const dlnActCoeffdlnN)
Get the array of derivatives of the ln activity coefficients with respect to the ln species mole numb...
std::vector< doublereal > d2lnActCoeffdT2_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
virtual void getdlnActCoeffdT(doublereal *dlnActCoeffdT) const
Get the array of temperature derivatives of the log activity coefficients.
std::vector< doublereal > dlnActCoeffdlnN_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
void resizeNumInteractions(const size_t num)
Resize internal arrays within the object that depend upon the number of binary Margules interaction t...
void initLengths()
Initialize lengths of local variables after all species have been identified.
virtual void getPartialMolarCp(doublereal *cpbar) const
Returns an array of partial molar entropies for the species in the mixture.
XML_Node & child(const size_t n) const
Return a changeable reference to the n'th child of the current node.
Definition: xml.cpp:584
virtual void getdlnActCoeffds(const doublereal dTds, const doublereal *const dXds, doublereal *dlnActCoeffds) const
Get the change in activity coefficients w.r.t.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
vector_fp m_SE_d_ij
Entropy term for the quaternary mole fraction interaction of the excess gibbs free energy expression...
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:229
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty thermophase object.
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine for objects which inherit from ThermoPhase.
vector_fp m_SE_b_ij
Entropy term for the binary mole fraction interaction of the excess gibbs free energy expression...
void s_update_dlnActCoeff_dlnN() const
Update the derivative of the log of the activity coefficients wrt log(moles_m)
void getElectrochemPotentials(doublereal *mu) const
Get the species electrochemical potentials.
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
std::string name() const
Returns the name of the XML node.
Definition: xml.h:390
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies for the species in the mixture.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:574
vector_fp m_VHE_b_ij
Enthalpy term for the binary mole fraction interaction of the excess gibbs free energy expression...
void s_update_dlnActCoeff_dlnN_diag() const
Update the derivative of the log of the activity coefficients wrt log(moles) - diagonal only...
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual void getLnActivityCoefficients(doublereal *lnac) const
Get the array of non-dimensional molar-based ln activity coefficients at the current solution tempera...
vector_fp m_VSE_b_ij
Entropy term for the binary mole fraction interaction of the excess gibbs free energy expression...
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
std::vector< doublereal > dlnActCoeffdlnX_diag_
Storage for the current derivative values of the gradients with respect to logarithm of the mole frac...
std::vector< size_t > m_pSpecies_B_ij
vector of species indices representing species B in the interaction
virtual void getdlnActCoeffdlnX_diag(doublereal *dlnActCoeffdlnX_diag) const
Get the array of log concentration-like derivatives of the log activity coefficients - diagonal compo...
MargulesVPSSTP & operator=(const MargulesVPSSTP &b)
Assignment operator.
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
vector_fp m_VHE_c_ij
Enthalpy term for the ternary mole fraction interaction of the excess gibbs free energy expression...
std::string id() const
Return the id attribute, if present.
Definition: xml.cpp:467
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:139
std::vector< doublereal > lnActCoeff_Scaled_
Storage for the current values of the activity coefficients of the species.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
MargulesVPSSTP()
Constructor.
std::vector< doublereal > moleFractions_
Storage for the current values of the mole fractions of the species.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
doublereal err(const std::string &msg) const
Error function.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize a ThermoPhase object, potentially reading activity coefficient information from an XML dat...
Contains declarations for string manipulation functions within Cantera.
void s_update_lnActCoeff() const
Update the activity coefficients.
vector_fp m_SE_c_ij
Entropy term for the ternary mole fraction interaction of the excess gibbs free energy expression...
virtual int eosType() const
Equation of state type flag.
MargulesVPSSTP is a derived class of GibbsExcessVPSSTP that employs the Margules approximation for th...
vector_fp m_HE_c_ij
Enthalpy term for the ternary mole fraction interaction of the excess gibbs free energy expression...
int formTempModel_
form of the temperature dependence of the Margules interaction expression
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
size_t m_kk
Number of species in the phase.
Definition: Phase.h:716
Array2D dlnActCoeffdlnN_
Storage for the current derivative values of the gradients with respect to logarithm of the species m...
void zero()
Set all of the entries to zero.
Definition: Array.h:256
virtual void initThermoFile(const std::string &inputFile, const std::string &id)
void s_update_dlnActCoeff_dT() const
Update the derivative of the log of the activity coefficients wrt T.
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
void s_update_dlnActCoeff_dlnX_diag() const
Update the derivative of the log of the activity coefficients wrt log(mole fraction) ...
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:246
size_t getFloatArray(const Cantera::XML_Node &node, std::vector< doublereal > &v, const bool convert, const std::string &unitsString, const std::string &nodeName)
This function reads the current node or a child node of the current node with the default name...
Definition: ctml.cpp:419
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:594
std::vector< doublereal > dlnActCoeffdT_Scaled_
Storage for the current derivative values of the gradients with respect to temperature of the log of ...
virtual void getdlnActCoeffdlnN_diag(doublereal *dlnActCoeffdlnN_diag) const
Get the array of derivatives of the log activity coefficients wrt mole numbers - diagonal only...
doublereal charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:504