Cantera  2.0
RedlichKwongMFTP.cpp
Go to the documentation of this file.
1 /**
2  * @file RedlichKwongMFTP.cpp
3  * Definition file for a derived class of ThermoPhase that assumes either
4  * an ideal gas or ideal solution approximation and handles
5  * variable pressure standard state methods for calculating
6  * thermodynamic properties (see \ref thermoprops and
7  * class \link Cantera::RedlichKwongMFTP RedlichKwongMFTP\endlink).
8  */
9 /*
10  * Copyright (2005) Sandia Corporation. Under the terms of
11  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
12  * U.S. Government retains certain rights in this software.
13  */
14 
16 
17 #include "cantera/thermo/mix_defs.h"
21 
22 using namespace std;
23 
24 namespace Cantera
25 {
26 
27 const doublereal RedlichKwongMFTP::omega_a = 4.27480233540E-01;
28 const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
29 const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
30 
31 //====================================================================================================================
32 /*
33  * Default constructor
34  */
35 RedlichKwongMFTP::RedlichKwongMFTP() :
37  m_standardMixingRules(0),
38  m_formTempParam(0),
39  m_b_current(0.0),
40  m_a_current(0.0),
41  a_vec_Curr_(0),
42  b_vec_Curr_(0),
43  a_coeff_vec(0,0),
44  m_pc_Species(0),
45  m_tc_Species(0),
46  m_vc_Species(0),
47  NSolns_(0),
48  m_pp(0),
49  m_tmpV(0),
50  m_partialMolarVolumes(0),
51  dpdV_(0.0),
52  dpdT_(0.0),
53  dpdni_(0)
54 {
55  Vroot_[0] = 0.0;
56  Vroot_[1] = 0.0;
57  Vroot_[2] = 0.0;
58 }
59 //====================================================================================================================
60 RedlichKwongMFTP::RedlichKwongMFTP(std::string infile, std::string id) :
62  m_standardMixingRules(0),
63  m_formTempParam(0),
64  m_b_current(0.0),
65  m_a_current(0.0),
66  a_vec_Curr_(0),
67  b_vec_Curr_(0),
68  a_coeff_vec(0,0),
69  m_pc_Species(0),
70  m_tc_Species(0),
71  m_vc_Species(0),
72  NSolns_(0),
73  m_pp(0),
74  m_tmpV(0),
75  m_partialMolarVolumes(0),
76  dpdV_(0.0),
77  dpdT_(0.0),
78  dpdni_(0)
79 {
80  Vroot_[0] = 0.0;
81  Vroot_[1] = 0.0;
82  Vroot_[2] = 0.0;
83  XML_Node* root = get_XML_File(infile);
84  if (id == "-") {
85  id = "";
86  }
87  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root);
88  if (!xphase) {
89  throw CanteraError("newPhase",
90  "Couldn't find phase named \"" + id + "\" in file, " + infile);
91  }
92  importPhase(*xphase, this);
93 }
94 //====================================================================================================================
95 RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, std::string id) :
97  m_standardMixingRules(0),
98  m_formTempParam(0),
99  m_b_current(0.0),
100  m_a_current(0.0),
101  a_vec_Curr_(0),
102  b_vec_Curr_(0),
103  a_coeff_vec(0,0),
104  m_pc_Species(0),
105  m_tc_Species(0),
106  m_vc_Species(0),
107  NSolns_(0),
108  m_pp(0),
109  m_tmpV(0),
110  m_partialMolarVolumes(0),
111  dpdV_(0.0),
112  dpdT_(0.0),
113  dpdni_(0)
114 {
115  Vroot_[0] = 0.0;
116  Vroot_[1] = 0.0;
117  Vroot_[2] = 0.0;
118  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, &phaseRefRoot);
119  if (!xphase) {
120  throw CanteraError("RedlichKwongMFTP::RedlichKwongMFTP()","Couldn't find phase named \"" + id + "\" in XML node");
121  }
122  importPhase(*xphase, this);
123 }
124 
125 //====================================================================================================================
128  m_standardMixingRules(0),
129  m_formTempParam(0),
130  m_b_current(0.0),
131  m_a_current(0.0),
132  a_vec_Curr_(0),
133  b_vec_Curr_(0),
134  a_coeff_vec(0,0),
135  m_pc_Species(0),
136  m_tc_Species(0),
137  m_vc_Species(0),
138  NSolns_(0),
139  m_pp(0),
140  m_tmpV(0),
141  m_partialMolarVolumes(0),
142  dpdV_(0.0),
143  dpdT_(0.0),
144  dpdni_(0)
145 {
146  std::string infile = "co2_redlichkwong.xml";
147  std::string id;
148  if (testProb == 1) {
149  infile = "co2_redlichkwong.xml";
150  id = "carbondioxide";
151  } else {
152  throw CanteraError("", "test prob = 1 only");
153  }
154  XML_Node* root = get_XML_File(infile);
155  if (id == "-") {
156  id = "";
157  }
158  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id, root);
159  if (!xphase) {
160  throw CanteraError("newPhase", "Couldn't find phase named \"" + id + "\" in file, " + infile);
161  }
162  importPhase(*xphase, this);
163 }
164 //====================================================================================================================
165 /*
166  * Copy Constructor:
167  *
168  * Note this stuff will not work until the underlying phase
169  * has a working copy constructor.
170  *
171  * The copy constructor just calls the assignment operator
172  * to do the heavy lifting.
173  */
176  m_standardMixingRules(0),
177  m_formTempParam(0),
178  m_b_current(0.0),
179  m_a_current(0.0),
180  a_vec_Curr_(0),
181  b_vec_Curr_(0),
182  a_coeff_vec(0,0),
183  m_pc_Species(0),
184  m_tc_Species(0),
185  m_vc_Species(0),
186  NSolns_(0),
187  m_pp(0),
188  m_tmpV(0),
189  m_partialMolarVolumes(0),
190  dpdV_(0.0),
191  dpdT_(0.0),
192  dpdni_(0)
193 {
194  *this = b;
195 }
196 
197 //====================================================================================================================
198 /*
199  * operator=()
200  *
201  * Note this stuff will not work until the underlying phase
202  * has a working assignment operator
203  */
206 {
207  if (&b != this) {
208  /*
209  * Mostly, this is a passthrough to the underlying
210  * assignment operator for the ThermoPhae parent object.
211  */
213  /*
214  * However, we have to handle data that we own.
215  */
220  a_vec_Curr_ = b.a_vec_Curr_;
221  b_vec_Curr_ = b.b_vec_Curr_;
222  a_coeff_vec = b.a_coeff_vec;
223 
224  m_pc_Species = b.m_pc_Species;
225  m_tc_Species = b.m_tc_Species;
226  m_vc_Species = b.m_vc_Species;
227  NSolns_ = b.NSolns_;
228  Vroot_[0] = b.Vroot_[0];
229  Vroot_[1] = b.Vroot_[1];
230  Vroot_[2] = b.Vroot_[2];
231  m_pp = b.m_pp;
232  m_tmpV = b.m_tmpV;
233  m_partialMolarVolumes = b.m_partialMolarVolumes;
234  dpdV_ = b.dpdV_;
235  dpdT_ = b.dpdT_;
236  dpdni_ = b.dpdni_;
237  }
238  return *this;
239 }
240 //====================================================================================================================
241 /*
242  * ~RedlichKwongMFTP(): (virtual)
243  *
244  */
246 {
247 }
248 //====================================================================================================================
249 /*
250  * Duplication function.
251  * This calls the copy constructor for this object.
252  */
254 {
255  RedlichKwongMFTP* vptp = new RedlichKwongMFTP(*this);
256  return (ThermoPhase*) vptp;
257 }
258 //====================================================================================================================
260 {
261  return cRedlichKwongMFTP;
262 }
263 
264 //====================================================================================================================
265 /*
266  * ------------Molar Thermodynamic Properties -------------------------
267  */
268 //====================================================================================================================
269 // Molar enthalpy. Units: J/kmol.
271 {
273  doublereal rt = _RT();
274  doublereal h_ideal = rt * mean_X(DATA_PTR(m_h0_RT));
275  doublereal h_nonideal = hresid();
276  return (h_ideal + h_nonideal);
277 }
278 //====================================================================================================================
279 // Molar internal energy. Units: J/kmol.
281 {
282  doublereal p0 = pressure();
283  doublereal md = molarDensity();
284  return (enthalpy_mole() - p0 / md);
285 }
286 //====================================================================================================================
287 // Molar entropy. Units: J/kmol/K.
289 {
291  doublereal sr_ideal = GasConstant * (mean_X(DATA_PTR(m_s0_R))
292  - sum_xlogx() - std::log(pressure()/m_spthermo->refPressure()));
293  doublereal sr_nonideal = sresid();
294  return (sr_ideal + sr_nonideal);
295 }
296 //====================================================================================================================
297 // Molar Gibbs function. Units: J/kmol.
299 {
300  return enthalpy_mole() - temperature() * entropy_mole();
301 }
302 //====================================================================================================================
303 /// Molar heat capacity at constant pressure. Units: J/kmol/K.
304 doublereal RedlichKwongMFTP::cp_mole() const
305 {
307  doublereal TKelvin = temperature();
308  doublereal sqt = sqrt(TKelvin);
309  doublereal mv = molarVolume();
310  doublereal vpb = mv + m_b_current;
312  doublereal cpref = GasConstant * mean_X(DATA_PTR(m_cp0_R));
313  doublereal dadt = da_dt();
314  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
315  doublereal dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
316  +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
317  double cp = dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
318  return cp;
319 }
320 //====================================================================================================================
321 /// Molar heat capacity at constant volume. Units: J/kmol/K.
322 doublereal RedlichKwongMFTP::cv_mole() const
323 {
324  throw CanteraError("", "unimplemented");
325  return cp_mole() - GasConstant;
326 }
327 //====================================================================================================================
328 // Return the thermodynamic pressure (Pa).
329 /*
330  * Since the mass density, temperature, and mass fractions are stored,
331  * this method uses these values to implement the
332  * mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$.
333  *
334  * \f[
335  * P = \frac{RT}{v-b_{mix}} - \frac{a_{mix}}{T^{0.5} v \left( v + b_{mix} \right) }
336  * \f]
337  *
338  */
339 doublereal RedlichKwongMFTP::pressure() const
340 {
341 
342 
343 #ifdef DEBUG_MODE
345 
346  // Get a copy of the private variables stored in the State object
347  double rho = density();
348  doublereal T = temperature();
349  doublereal mmw = meanMolecularWeight();
350  double molarV = mmw / rho;
351 
352  double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
353 
354  if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) {
355  throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe");
356  }
357 #endif
358 
359  return m_Pcurrent;
360 }
361 //====================================================================================================================
363 {
364  /*
365  * Calculate the molarVolume of the solution (m**3 kmol-1)
366  */
367 
368  const doublereal* const dtmp = moleFractdivMMW();
370  double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp);
371  /*
372  * Set the density in the parent State object directly,
373  * by calling the Phase::setDensity() function.
374  */
375  double dens = 1.0/invDens;
376  Phase::setDensity(dens);
377 
378 }
379 
380 //====================================================================================================================
381 void RedlichKwongMFTP::setTemperature(const doublereal temp)
382 {
383  Phase::setTemperature(temp);
385  updateAB();
386 }
387 //====================================================================================================================
388 void RedlichKwongMFTP::setMassFractions(const doublereal* const x)
389 {
391  updateAB();
392 }
393 //====================================================================================================================
394 void RedlichKwongMFTP::setMassFractions_NoNorm(const doublereal* const x)
395 {
397  updateAB();
398 }
399 //====================================================================================================================
400 void RedlichKwongMFTP::setMoleFractions(const doublereal* const x)
401 {
403  updateAB();
404 }
405 //====================================================================================================================
406 void RedlichKwongMFTP::setMoleFractions_NoNorm(const doublereal* const x)
407 {
409  updateAB();
410 }
411 //====================================================================================================================
412 void RedlichKwongMFTP::setConcentrations(const doublereal* const c)
413 {
415  updateAB();
416 }
417 
418 //====================================================================================================================
420 {
421 
422 
423  throw CanteraError("RedlichKwongMFTP::isothermalCompressibility() ",
424  "not implemented");
425 
426  return 0.0;
427 }
428 //====================================================================================================================
430 {
431  getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes));
432  for (size_t k = 0; k < m_kk; k++) {
433  c[k] = moleFraction(k) / m_partialMolarVolumes[k];
434  }
435 }
436 //====================================================================================================================
437 /*
438  * Returns the standard concentration \f$ C^0_k \f$, which is used to normalize
439  * the generalized concentration.
440  */
441 doublereal RedlichKwongMFTP::standardConcentration(size_t k) const
442 {
443 
445 
446  return 1.0 / m_tmpV[k];
447 
448 
449 
450 }
451 //====================================================================================================================
452 /*
453  * Returns the natural logarithm of the standard
454  * concentration of the kth species
455  */
456 doublereal RedlichKwongMFTP::logStandardConc(size_t k) const
457 {
458  double c = standardConcentration(k);
459  double lc = std::log(c);
460  return lc;
461 }
462 //====================================================================================================================
463 /*
464  *
465  * getUnitsStandardConcentration()
466  *
467  * Returns the units of the standard and general concentrations
468  * Note they have the same units, as their divisor is
469  * defined to be equal to the activity of the kth species
470  * in the solution, which is unitless.
471  *
472  * This routine is used in print out applications where the
473  * units are needed. Usually, MKS units are assumed throughout
474  * the program and in the XML input files.
475  *
476  * uA[0] = kmol units - default = 1
477  * uA[1] = m units - default = -nDim(), the number of spatial
478  * dimensions in the Phase class.
479  * uA[2] = kg units - default = 0;
480  * uA[3] = Pa(pressure) units - default = 0;
481  * uA[4] = Temperature units - default = 0;
482  * uA[5] = time units - default = 0
483  *
484  * For EOS types other than cIdealSolidSolnPhase1, the default
485  * kmol/m3 holds for standard concentration units. For
486  * cIdealSolidSolnPhase0 type, the standard concentration is
487  * unitless.
488  */
489 void RedlichKwongMFTP::getUnitsStandardConc(double* uA, int, int sizeUA) const
490 {
491  //int eos = eosType();
492 
493  for (int i = 0; i < sizeUA; i++) {
494  if (i == 0) {
495  uA[0] = 1.0;
496  }
497  if (i == 1) {
498  uA[1] = -static_cast<int>(nDim());
499  }
500  if (i == 2) {
501  uA[2] = 0.0;
502  }
503  if (i == 3) {
504  uA[3] = 0.0;
505  }
506  if (i == 4) {
507  uA[4] = 0.0;
508  }
509  if (i == 5) {
510  uA[5] = 0.0;
511  }
512  }
513 
514 }
515 
516 //====================================================================================================================
517 //! Get the array of non-dimensional activity coefficients at
518 //! the current solution temperature, pressure, and solution concentration.
519 /*!
520  * For ideal gases, the activity coefficients are all equal to one.
521  *
522  * @param ac Output vector of activity coefficients. Length: m_kk.
523  */
525 {
526  doublereal TKelvin = temperature();
527  doublereal rt = TKelvin * GasConstant;
528  doublereal mv = molarVolume();
529  doublereal sqt = sqrt(TKelvin);
530  doublereal vpb = mv + m_b_current;
531  doublereal vmb = mv - m_b_current;
532 
533  for (size_t k = 0; k < m_kk; k++) {
534  m_pp[k] = 0.0;
535  for (size_t i = 0; i < m_kk; i++) {
536  size_t counter = k + m_kk*i;
537  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
538  }
539  }
540  doublereal pres = pressure();
541 
542  for (size_t k = 0; k < m_kk; k++) {
543  ac[k] = (- rt * log(pres * mv / rt)
544  + rt * log(mv / vmb)
545  + rt * b_vec_Curr_[k] / vmb
546  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
547  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
548  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
549  );
550  }
551  for (size_t k = 0; k < m_kk; k++) {
552  ac[k] = exp(ac[k]/rt);
553  }
554 }
555 //====================================================================================================================
556 /*
557  * ---- Partial Molar Properties of the Solution -----------------
558  */
559 //====================================================================================================================
560 /*
561  * Get the array of non-dimensional species chemical potentials
562  * These are partial molar Gibbs free energies.
563  * \f$ \mu_k / \hat R T \f$.
564  * Units: unitless
565  *
566  * We close the loop on this function, here, calling
567  * getChemPotentials() and then dividing by RT.
568  */
569 void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
570 {
571  getChemPotentials(muRT);
572  doublereal invRT = 1.0 / _RT();
573  for (size_t k = 0; k < m_kk; k++) {
574  muRT[k] *= invRT;
575  }
576 }
577 //====================================================================================================================
578 void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
579 {
580  getGibbs_ref(mu);
581  doublereal xx;
582  doublereal rt = temperature() * GasConstant;
583  for (size_t k = 0; k < m_kk; k++) {
585  mu[k] += rt*(log(xx));
586  }
587 
588  doublereal TKelvin = temperature();
589  doublereal mv = molarVolume();
590  doublereal sqt = sqrt(TKelvin);
591  doublereal vpb = mv + m_b_current;
592  doublereal vmb = mv - m_b_current;
593 
594  for (size_t k = 0; k < m_kk; k++) {
595  m_pp[k] = 0.0;
596  for (size_t i = 0; i < m_kk; i++) {
597  size_t counter = k + m_kk*i;
598  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
599  }
600  }
601  doublereal pres = pressure();
602  doublereal refP = refPressure();
603 
604  for (size_t k = 0; k < m_kk; k++) {
605  mu[k] += (rt * log(pres/refP) - rt * log(pres * mv / rt)
606  + rt * log(mv / vmb)
607  + rt * b_vec_Curr_[k] / vmb
608  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
609  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
610  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
611  );
612  }
613 }
614 //====================================================================================================================
616 {
617  /*
618  * First we get the reference state contributions
619  */
620  getEnthalpy_RT_ref(hbar);
621  doublereal rt = GasConstant * temperature();
622  scale(hbar, hbar+m_kk, hbar, rt);
623 
624  /*
625  * We calculate dpdni_
626  */
627  doublereal TKelvin = temperature();
628  doublereal mv = molarVolume();
629  doublereal sqt = sqrt(TKelvin);
630 
631  doublereal vpb = mv + m_b_current;
632  doublereal vmb = mv - m_b_current;
633 
634  for (size_t k = 0; k < m_kk; k++) {
635  m_pp[k] = 0.0;
636  for (size_t i = 0; i < m_kk; i++) {
637  size_t counter = k + m_kk*i;
638  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
639  }
640  }
641 
642 
643 
644  for (size_t k = 0; k < m_kk; k++) {
645  dpdni_[k] = rt/vmb + rt * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
646  + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
647  }
648  doublereal dadt = da_dt();
649  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
650 
651  for (size_t k = 0; k < m_kk; k++) {
652  m_tmpV[k] = 0.0;
653  for (size_t i = 0; i < m_kk; i++) {
654  size_t counter = k + m_kk*i;
655  m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
656  }
657  }
658 
660  doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_;
661 
662  for (size_t k = 0; k < m_kk; k++) {
663  double hE_v = (mv * dpdni_[k] - rt - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
664  + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
665  + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
666  hbar[k] = hbar[k] + hE_v;
667 
668 
669  hbar[k] -= fac2 * dpdni_[k];
670  }
671 
672 }
673 //====================================================================================================================
674 void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
675 {
676  getEntropy_R_ref(sbar);
677  doublereal r = GasConstant;
678  scale(sbar, sbar+m_kk, sbar, r);
679  doublereal TKelvin = temperature();
680  doublereal sqt = sqrt(TKelvin);
681  doublereal mv = molarVolume();
682  doublereal refP = refPressure();
683 
684  for (size_t k = 0; k < m_kk; k++) {
685  doublereal xx = std::max(SmallNumber, moleFraction(k));
686  sbar[k] += r * (- log(xx));
687  }
688 
689  for (size_t k = 0; k < m_kk; k++) {
690  m_pp[k] = 0.0;
691  for (size_t i = 0; i < m_kk; i++) {
692  size_t counter = k + m_kk*i;
693  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
694  }
695  }
696 
697  for (size_t k = 0; k < m_kk; k++) {
698  m_tmpV[k] = 0.0;
699  for (size_t i = 0; i < m_kk; i++) {
700  size_t counter = k + m_kk*i;
701  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
702  }
703  }
704 
705 
706  doublereal dadt = da_dt();
707  doublereal fac = dadt - m_a_current / (2.0 * TKelvin);
708  doublereal vmb = mv - m_b_current;
709  doublereal vpb = mv + m_b_current;
710 
711 
712  for (size_t k = 0; k < m_kk; k++) {
713  sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
714  + GasConstant
715  + GasConstant * log(mv/vmb)
716  + GasConstant * b_vec_Curr_[k]/vmb
717  + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
718  - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
719  + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
720  - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
721  ) ;
722  }
723 
725  getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes));
726  for (size_t k = 0; k < m_kk; k++) {
727  sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
728  }
729 }
730 //====================================================================================================================
732 {
733  getIntEnergy_RT(ubar);
734  doublereal rt = GasConstant * temperature();
735  scale(ubar, ubar+m_kk, ubar, rt);
736 }
737 //====================================================================================================================
738 void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const
739 {
740  getCp_R(cpbar);
741  doublereal r = GasConstant;
742  scale(cpbar, cpbar+m_kk, cpbar, r);
743 }
744 //====================================================================================================================
745 void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
746 {
747  // getStandardVolumes(vbar);
748 
749 
750  for (size_t k = 0; k < m_kk; k++) {
751  m_pp[k] = 0.0;
752  for (size_t i = 0; i < m_kk; i++) {
753  size_t counter = k + m_kk*i;
754  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
755  }
756  }
757 
758  for (size_t k = 0; k < m_kk; k++) {
759  m_tmpV[k] = 0.0;
760  for (size_t i = 0; i < m_kk; i++) {
761  size_t counter = k + m_kk*i;
762  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
763  }
764  }
765 
766  doublereal TKelvin = temperature();
767  doublereal sqt = sqrt(TKelvin);
768  doublereal mv = molarVolume();
769 
770  doublereal rt = GasConstant * TKelvin;
771 
772  doublereal vmb = mv - m_b_current;
773  doublereal vpb = mv + m_b_current;
774 
775  for (size_t k = 0; k < m_kk; k++) {
776 
777  doublereal num = (rt + rt * m_b_current/ vmb + rt * b_vec_Curr_[k] / vmb
778  + rt * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
779  - 2.0 * m_pp[k] / (sqt * vpb)
780  + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
781  );
782 
783  doublereal denom = (m_Pcurrent + rt * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
784  );
785 
786  vbar[k] = num / denom;
787  }
788 
789 }
790 //====================================================================================================================
792 {
793  double pc, tc, vc;
794  double a0 = 0.0;
795  double aT = 0.0;
796  for (size_t i = 0; i < m_kk; i++) {
797  for (size_t j = 0; j <m_kk; j++) {
798  size_t counter = i + m_kk * j;
799  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
800  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
801  }
802  }
803  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
804  return tc;
805 }
806 //====================================================================================================================
808 {
809  double pc, tc, vc;
810  double a0 = 0.0;
811  double aT = 0.0;
812  for (size_t i = 0; i < m_kk; i++) {
813  for (size_t j = 0; j <m_kk; j++) {
814  size_t counter = i + m_kk * j;
815  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
816  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
817  }
818  }
819  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
820 
821  return pc;
822 }
823 //====================================================================================================================
825 {
826  double pc, tc, vc;
827  double a0 = 0.0;
828  double aT = 0.0;
829  for (size_t i = 0; i < m_kk; i++) {
830  for (size_t j = 0; j <m_kk; j++) {
831  size_t counter = i + m_kk * j;
832  a0 += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(0, counter);
833  aT += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(1, counter);
834  }
835  }
836  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
837 
838  double mmw = meanMolecularWeight();
839  return mmw / vc;
840 }
841 //====================================================================================================================
842 
843 /*
844  * ----- Thermodynamic Values for the Species Reference States ----
845  */
846 
847 
848 //====================================================================================================================
849 
850 /*
851  * Perform initializations after all species have been
852  * added.
853  */
855 {
856  initLengths();
858 }
859 
860 //====================================================================================================================
861 void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
862 {
863  double tmp, tmp2;
865 
867 
868 
869  /*
870  * Within the method, we protect against inf results if the
871  * exponent is too high.
872  *
873  * If it is too low, we set
874  * the partial pressure to zero. This capability is needed
875  * by the elemental potential method.
876  */
877  doublereal pres = 0.0;
878  double m_p0 = refPressure();
879  for (size_t k = 0; k < m_kk; k++) {
880  tmp = -m_tmpV[k] + mu_RT[k];
881  if (tmp < -600.) {
882  m_pp[k] = 0.0;
883  } else if (tmp > 500.0) {
884  tmp2 = tmp / 500.;
885  tmp2 *= tmp2;
886  m_pp[k] = m_p0 * exp(500.) * tmp2;
887  } else {
888  m_pp[k] = m_p0 * exp(tmp);
889  }
890  pres += m_pp[k];
891  }
892  // set state
893  setState_PX(pres, &m_pp[0]);
894 }
895 //====================================================================================================================
896 /*
897  * Initialize the internal lengths.
898  * (this is not a virtual function)
899  */
901 {
902 
903 
904  a_vec_Curr_.resize(m_kk * m_kk, 0.0);
905  b_vec_Curr_.resize(m_kk, 0.0);
906 
907  a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
908 
909 
910  m_pc_Species.resize(m_kk, 0.0);
911  m_tc_Species.resize(m_kk, 0.0);
912  m_vc_Species.resize(m_kk, 0.0);
913 
914 
915  m_pp.resize(m_kk, 0.0);
916  m_tmpV.resize(m_kk, 0.0);
917  m_partialMolarVolumes.resize(m_kk, 0.0);
918  dpdni_.resize(m_kk, 0.0);
919 }
920 //====================================================================================================================
921 /*
922  * Import and initialize a ThermoPhase object
923  *
924  * param phaseNode This object must be the phase node of a
925  * complete XML tree
926  * description of the phase, including all of the
927  * species data. In other words while "phase" must
928  * point to an XML phase object, it must have
929  * sibling nodes "speciesData" that describe
930  * the species in the phase.
931  * param id ID of the phase. If nonnull, a check is done
932  * to see if phaseNode is pointing to the phase
933  * with the correct id.
934  *
935  * This routine initializes the lengths in the current object and
936  * then calls the parent routine.
937  */
938 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, std::string id)
939 {
941 
942  /*
943  * Check the model parameter for the Redlich-Kwong equation of state
944  * two are allowed
945  * RedlichKwong mixture of species, each of which are RK fluids
946  * RedlichKwongMFTP mixture of species with cross term coefficients
947  */
948  if (phaseNode.hasChild("thermo")) {
949  XML_Node& thermoNode = phaseNode.child("thermo");
950  std::string model = thermoNode["model"];
951  if (model == "RedlichKwong") {
953  } else if (model == "RedlichKwongMFTP") {
955  } else {
956  throw CanteraError("RedlichKwongMFTP::initThermoXML",
957  "Unknown thermo model : " + model);
958  }
959 
960 
961  /*
962  * Go get all of the coefficients and factors in the
963  * activityCoefficients XML block
964  */
965  XML_Node* acNodePtr = 0;
966  if (thermoNode.hasChild("activityCoefficients")) {
967  XML_Node& acNode = thermoNode.child("activityCoefficients");
968  acNodePtr = &acNode;
969  size_t nC = acNode.nChildren();
970 
971  /*
972  * Loop through the children getting multiple instances of
973  * parameters
974  */
975  for (size_t i = 0; i < nC; i++) {
976  XML_Node& xmlACChild = acNodePtr->child(i);
977  string stemp = xmlACChild.name();
978  string nodeName = lowercase(stemp);
979  /*
980  * Process a binary salt field, or any of the other XML fields
981  * that make up the Pitzer Database. Entries will be ignored
982  * if any of the species in the entry isn't in the solution.
983  */
984  if (nodeName == "purefluidparameters") {
985  readXMLPureFluid(xmlACChild);
986  }
987  }
988  if (m_standardMixingRules == 1) {
990  }
991  /*
992  * Loop through the children getting multiple instances of
993  * parameters
994  */
995  for (size_t i = 0; i < nC; i++) {
996  XML_Node& xmlACChild = acNodePtr->child(i);
997  string stemp = xmlACChild.name();
998  string nodeName = lowercase(stemp);
999  /*
1000  * Process a binary salt field, or any of the other XML fields
1001  * that make up the Pitzer Database. Entries will be ignored
1002  * if any of the species in the entry isn't in the solution.
1003  */
1004  if (nodeName == "crossfluidparameters") {
1005  readXMLCrossFluid(xmlACChild);
1006  }
1007  }
1008 
1009  }
1010  }
1011 
1012  for (size_t i = 0; i < m_kk; i++) {
1013  double a0coeff = a_coeff_vec(0, i*m_kk + i);
1014  double aTcoeff = a_coeff_vec(1, i*m_kk + i);
1015  double ai = a0coeff + aTcoeff * 500.;
1016  double bi = b_vec_Curr_[i];
1017  calcCriticalConditions(ai, bi, a0coeff, aTcoeff, m_pc_Species[i], m_tc_Species[i], m_vc_Species[i]);
1018  }
1019 
1020  MixtureFugacityTP::initThermoXML(phaseNode, id);
1021 }
1022 //====================================================================================================================
1023 
1025 {
1026  vector_fp vParams;
1027  string xname = pureFluidParam.name();
1028  if (xname != "pureFluidParameters") {
1029  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
1030  "Incorrect name for processing this routine: " + xname);
1031  }
1032 
1033  /*
1034  * Read the species
1035  * Find the index of the species in the current phase. It's not an error to not find the species
1036  */
1037  string iName = pureFluidParam.attrib("species");
1038  if (iName == "") {
1039  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute");
1040  }
1041  size_t iSpecies = speciesIndex(iName);
1042  if (iSpecies == npos) {
1043  return;
1044  }
1045  size_t counter = iSpecies + m_kk * iSpecies;
1046  size_t nParamsExpected, nParamsFound;
1047  size_t num = pureFluidParam.nChildren();
1048  for (size_t iChild = 0; iChild < num; iChild++) {
1049  XML_Node& xmlChild = pureFluidParam.child(iChild);
1050  string stemp = xmlChild.name();
1051  string nodeName = lowercase(stemp);
1052 
1053  if (nodeName == "a_coeff") {
1054  string iModel = lowercase(xmlChild.attrib("model"));
1055  if (iModel == "constant") {
1056  nParamsExpected = 1;
1057  } else if (iModel == "linear_a") {
1058  nParamsExpected = 2;
1059  if (m_formTempParam == 0) {
1060  m_formTempParam = 1;
1061  }
1062  } else {
1063  throw CanteraError("", "unknown model");
1064  }
1065 
1066  ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
1067  nParamsFound = vParams.size();
1068  if (nParamsFound != nParamsExpected) {
1069  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for a_coeff" + iName + ")",
1070  "wrong number of params found");
1071  }
1072 
1073  for (size_t i = 0; i < nParamsFound; i++) {
1074  a_coeff_vec(i, counter) = vParams[i];
1075  }
1076  } else if (nodeName == "b_coeff") {
1077  ctml::getFloatArray(xmlChild, vParams, true, "m3/kmol", "b_coeff");
1078  nParamsFound = vParams.size();
1079  if (nParamsFound != 1) {
1080  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for b_coeff" + iName + ")",
1081  "wrong number of params found");
1082  }
1083  b_vec_Curr_[iSpecies] = vParams[0];
1084  }
1085  }
1086 }
1087 //====================================================================================================================
1089 {
1090  int nParam = 2;
1091  for (size_t i = 0; i < m_kk; i++) {
1092  size_t icounter = i + m_kk * i;
1093  for (size_t j = 0; j < m_kk; j++) {
1094  if (i != j) {
1095  size_t counter = i + m_kk * j;
1096  size_t jcounter = j + m_kk * j;
1097  for (int n = 0; n < nParam; n++) {
1098  a_coeff_vec(n, counter) = sqrt(a_coeff_vec(n, icounter) * a_coeff_vec(n, jcounter));
1099  }
1100  }
1101  }
1102  }
1103 }
1104 //====================================================================================================================
1105 
1107 {
1108  vector_fp vParams;
1109  string xname = CrossFluidParam.name();
1110  if (xname != "crossFluidParameters") {
1111  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
1112  "Incorrect name for processing this routine: " + xname);
1113  }
1114 
1115  /*
1116  * Read the species
1117  * Find the index of the species in the current phase. It's not an error to not find the species
1118  */
1119  string iName = CrossFluidParam.attrib("species1");
1120  if (iName == "") {
1121  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute");
1122  }
1123  size_t iSpecies = speciesIndex(iName);
1124  if (iSpecies == npos) {
1125  return;
1126  }
1127  string jName = CrossFluidParam.attrib("species2");
1128  if (iName == "") {
1129  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species2 attribute");
1130  }
1131  size_t jSpecies = speciesIndex(jName);
1132  if (jSpecies == npos) {
1133  return;
1134  }
1135 
1136  size_t counter = iSpecies + m_kk * jSpecies;
1137  size_t counter0 = jSpecies + m_kk * iSpecies;
1138  size_t nParamsExpected, nParamsFound;
1139  size_t num = CrossFluidParam.nChildren();
1140  for (size_t iChild = 0; iChild < num; iChild++) {
1141  XML_Node& xmlChild = CrossFluidParam.child(iChild);
1142  string stemp = xmlChild.name();
1143  string nodeName = lowercase(stemp);
1144 
1145  if (nodeName == "a_coeff") {
1146  string iModel = lowercase(xmlChild.attrib("model"));
1147  if (iModel == "constant") {
1148  nParamsExpected = 1;
1149  } else if (iModel == "linear_a") {
1150  nParamsExpected = 2;
1151  if (m_formTempParam == 0) {
1152  m_formTempParam = 1;
1153  }
1154  } else {
1155  throw CanteraError("", "unknown model");
1156  }
1157 
1158  ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
1159  nParamsFound = vParams.size();
1160  if (nParamsFound != nParamsExpected) {
1161  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid(for a_coeff" + iName + ")",
1162  "wrong number of params found");
1163  }
1164 
1165  for (size_t i = 0; i < nParamsFound; i++) {
1166  a_coeff_vec(i, counter) = vParams[i];
1167  a_coeff_vec(i, counter0) = vParams[i];
1168  }
1169  }
1170  }
1171 }
1172 //====================================================================================================================
1174 {
1176  std::string model = thermoNode["model"];
1177 
1178 
1179 }
1180 //====================================================================================================================
1181 // Calculate the deviation terms for the total entropy of the mixture from the
1182 // ideal gas mixture
1183 /*
1184  * Here we use the current state conditions
1185  *
1186  * @return Returns the change in entropy in units of J kmol-1 K-1.
1187  */
1188 doublereal RedlichKwongMFTP::sresid() const
1189 {
1190  // note this agrees with tpx
1191  doublereal rho = density();
1192  doublereal mmw = meanMolecularWeight();
1193  doublereal molarV = mmw / rho;
1194  double hh = m_b_current / molarV;
1195  doublereal zz = z();
1196  doublereal dadt = da_dt();
1197  doublereal T = temperature();
1198  doublereal sqT = sqrt(T);
1199  doublereal fac = dadt - m_a_current / (2.0 * T);
1200  double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
1201  double sp = GasConstant * sresid_mol_R;
1202  return sp;
1203 }
1204 //====================================================================================================================
1205 // Calculate the deviation terms for the total enthalpy of the mixture from the
1206 // ideal gas mixture
1207 /*
1208  * Here we use the current state conditions
1209  *
1210  * @return Returns the change in entropy in units of J kmol-1.
1211  */
1212 doublereal RedlichKwongMFTP::hresid() const
1213 {
1214  // note this agrees with tpx
1215  doublereal rho = density();
1216  doublereal mmw = meanMolecularWeight();
1217  doublereal molarV = mmw / rho;
1218  double hh = m_b_current / molarV;
1219  doublereal zz = z();
1220  doublereal dadt = da_dt();
1221  doublereal T = temperature();
1222  doublereal sqT = sqrt(T);
1223  doublereal fac = T * dadt - 3.0 *m_a_current / (2.0);
1224  double hresid_mol = GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
1225  return hresid_mol;
1226 }
1227 //====================================================================================================================
1228 // Estimate for the molar volume of the liquid
1229 /*
1230  * Note: this is only used as a starting guess for later routines that actually calculate an
1231  * accurate value for the liquid molar volume.
1232  * This routine doesn't change the state of the system.
1233  *
1234  * @param TKelvin temperature in kelvin
1235  * @param pres Pressure in Pa. This is used as an initial guess. If the routine
1236  * needs to change the pressure to find a stable liquid state, the
1237  * new pressure is returned in this variable.
1238  *
1239  * @return Returns the estimate of the liquid volume. If the liquid can't be found, this
1240  * routine returns -1.
1241  */
1242 doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
1243 {
1244  double v = m_b_current * 1.1;
1245  double atmp;
1246  double btmp;
1247  calculateAB(TKelvin, atmp, btmp);
1248 
1249  doublereal pres = presGuess;
1250  double pp = psatEst(TKelvin);
1251  if (pres < pp) {
1252  pres = pp;
1253  }
1254  double Vroot[3];
1255 
1256  bool foundLiq = false;
1257  int m = 0;
1258  do {
1259 
1260  int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot);
1261 
1262  // printf("nsol = %d\n", nsol);
1263  // printf("liquidVolEst start: T = %g , p = %g, a = %g, b = %g\n", TKelvin, pres, m_a_current, m_b_current);
1264 
1265  if (nsol == 1 || nsol == 2) {
1266  double pc = critPressure();
1267  if (pres > pc) {
1268  foundLiq = true;
1269  }
1270  pres *= 1.04;
1271 
1272  } else {
1273  foundLiq = true;
1274  }
1275  } while ((m < 100) && (!foundLiq));
1276 
1277  if (foundLiq) {
1278  v = Vroot[0];
1279  presGuess = pres;
1280  } else {
1281  v = -1.0;
1282  }
1283  //printf (" RedlichKwongMFTP::liquidVolEst %g %g converged in %d its\n", TKelvin, pres, i);
1284  return v;
1285 }
1286 //====================================================================================================================
1287 // Calculates the density given the temperature and the pressure and a guess at the density.
1288 /*
1289  * Note, below T_c, this is a multivalued function. We do not cross the vapor dome in this.
1290  * This is protected because it is called during setState_TP() routines. Infinite loops would result
1291  * if it were not protected.
1292  *
1293  * -> why is this not const?
1294  *
1295  * parameters:
1296  * @param TKelvin Temperature in Kelvin
1297  * @param pressure Pressure in Pascals (Newton/m**2)
1298  * @param phaseReqested int representing the phase whose density we are requesting. If we put
1299  * a gas or liquid phase here, we will attempt to find a volume in that
1300  * part of the volume space, only, in this routine. A value of FLUID_UNDEFINED
1301  * means that we will accept anything.
1302  *
1303  * @param rhoguess Guessed density of the fluid. A value of -1.0 indicates that there
1304  * is no guessed density
1305  *
1306  *
1307  * @return We return the density of the fluid at the requested phase. If we have not found any
1308  * acceptable density we return a -1. If we have found an accectable density at a
1309  * different phase, we return a -2.
1310  */
1311 doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
1312 {
1313 
1314  /*
1315  * It's necessary to set the temperature so that m_a_current is set correctly.
1316  */
1317  setTemperature(TKelvin);
1318  double tcrit = critTemperature();
1319  doublereal mmw = meanMolecularWeight();
1320  double densBase = 0.0;
1321  if (rhoguess == -1.0) {
1322  if (phaseRequested != FLUID_GAS) {
1323  if (TKelvin > tcrit) {
1324  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1325  } else {
1326  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
1327  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1328  } else if (phaseRequested >= FLUID_LIQUID_0) {
1329  double lqvol = liquidVolEst(TKelvin, presPa);
1330  rhoguess = mmw / lqvol;
1331  }
1332  }
1333  } else {
1334  /*
1335  * Assume the Gas phase initial guess, if nothing is
1336  * specified to the routine
1337  */
1338  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1339  }
1340 
1341  }
1342 
1343 
1344  doublereal volguess = mmw / rhoguess;
1345  NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
1346 
1347  doublereal molarVolLast = Vroot_[0];
1348  if (NSolns_ >= 2) {
1349  if (phaseRequested >= FLUID_LIQUID_0) {
1350  molarVolLast = Vroot_[0];
1351  } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
1352  molarVolLast = Vroot_[2];
1353  } else {
1354  if (volguess > Vroot_[1]) {
1355  molarVolLast = Vroot_[2];
1356  } else {
1357  molarVolLast = Vroot_[0];
1358  }
1359  }
1360  } else if (NSolns_ == 1) {
1361  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
1362  molarVolLast = Vroot_[0];
1363  } else {
1364  //molarVolLast = Vroot_[0];
1365  //printf("DensityCalc(): Possible problem encountered\n");
1366  return -2.0;
1367  }
1368  } else if (NSolns_ == -1) {
1369  if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
1370  molarVolLast = Vroot_[0];
1371  } else if (TKelvin > tcrit) {
1372  molarVolLast = Vroot_[0];
1373  } else {
1374  // molarVolLast = Vroot_[0];
1375  // printf("DensityCalc(): Possible problem encountered\n");
1376  return -2.0;
1377  }
1378  } else {
1379  molarVolLast = Vroot_[0];
1380  //printf("DensityCalc(): Possible problem encountered\n");
1381  return -1.0;
1382  }
1383  densBase = mmw / molarVolLast;
1384  return densBase;
1385 }
1386 //====================================================================================================================
1387 // Return the value of the density at the liquid spinodal point (on the liquid side)
1388 // for the current temperature.
1389 /*
1390  * @return returns the density with units of kg m-3
1391  */
1393 {
1394  if (NSolns_ != 3) {
1395  double dens = critDensity();
1396  return dens;
1397  }
1398  double vmax = Vroot_[1];
1399  double vmin = Vroot_[0];
1400  RootFind rf(fdpdv_);
1401  rf.setPrintLvl(10);
1402  rf.setTol(1.0E-5, 1.0E-10);
1404 
1405  double vbest = 0.5 * (Vroot_[0]+Vroot_[1]);
1406  double funcNeeded = 0.0;
1407 
1408  int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
1409  if (status != ROOTFIND_SUCCESS) {
1410  throw CanteraError(" RedlichKwongMFTP::densSpinodalLiquid() ", "didn't converge");
1411  }
1412  doublereal mmw = meanMolecularWeight();
1413  doublereal rho = mmw / vbest;
1414  return rho;
1415 }
1416 //====================================================================================================================
1417 // Return the value of the density at the gas spinodal point (on the gas side)
1418 // for the current temperature.
1419 /*
1420  * @return returns the density with units of kg m-3
1421  */
1423 {
1424  if (NSolns_ != 3) {
1425  double dens = critDensity();
1426  return dens;
1427  }
1428  double vmax = Vroot_[2];
1429  double vmin = Vroot_[1];
1430  RootFind rf(fdpdv_);
1431  rf.setPrintLvl(10);
1432  rf.setTol(1.0E-5, 1.0E-10);
1434 
1435  double vbest = 0.5 * (Vroot_[1]+Vroot_[2]);
1436  double funcNeeded = 0.0;
1437 
1438  int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
1439  if (status != ROOTFIND_SUCCESS) {
1440  throw CanteraError(" RedlichKwongMFTP::densSpinodalGas() ", "didn't converge");
1441  }
1442  doublereal mmw = meanMolecularWeight();
1443  doublereal rho = mmw / vbest;
1444  return rho;
1445 }
1446 //====================================================================================================================
1447 // Calculate the pressure given the temperature and the molar volume
1448 /*
1449  * Calculate the pressure given the temperature and the molar volume
1450  *
1451  * @param TKelvin temperature in kelvin
1452  * @param molarVol molar volume ( m3/kmol)
1453  *
1454  * @return Returns the pressure.
1455  */
1456 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
1457 {
1458  doublereal sqt = sqrt(TKelvin);
1459  double pres = GasConstant * TKelvin / (molarVol - m_b_current)
1460  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
1461  return pres;
1462 }
1463 //====================================================================================================================
1464 // Calculate the pressure and the pressure derivative given the temperature and the molar volume
1465 /*
1466  * Temperature and mole number are held constant
1467  *
1468  * @param TKelvin temperature in kelvin
1469  * @param molarVol molar volume ( m3/kmol)
1470  *
1471  * @param presCalc Returns the pressure.
1472  *
1473  * @return Returns the derivative of the pressure wrt the molar volume
1474  */
1475 doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
1476 {
1477  doublereal sqt = sqrt(TKelvin);
1478  presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
1479  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
1480 
1481  doublereal vpb = molarVol + m_b_current;
1482  doublereal vmb = molarVol - m_b_current;
1483  doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb)
1484  + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
1485  return dpdv;
1486 }
1487 //====================================================================================================================
1488 
1490 {
1491  doublereal TKelvin = temperature();
1492  doublereal mv = molarVolume();
1493  doublereal pres;
1494 
1495  dpdV_ = dpdVCalc(TKelvin, mv, pres);
1496 
1497  doublereal sqt = sqrt(TKelvin);
1498  doublereal vpb = mv + m_b_current;
1499  doublereal vmb = mv - m_b_current;
1500  doublereal dadt = da_dt();
1501  doublereal fac = dadt - m_a_current/(2.0 * TKelvin);
1502 
1503  dpdT_ = (GasConstant / (vmb) - fac / (sqt * mv * vpb));
1504 }
1505 //====================================================================================================================
1506 void RedlichKwongMFTP::updateMixingExpressions()
1507 {
1508  updateAB();
1509 }
1510 //====================================================================================================================
1512 {
1513  double temp = temperature();
1514  if (m_formTempParam == 1) {
1515  for (size_t i = 0; i < m_kk; i++) {
1516  for (size_t j = 0; j < m_kk; j++) {
1517  size_t counter = i * m_kk + j;
1518  a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
1519  }
1520  }
1521  }
1522 
1523  m_b_current = 0.0;
1524  m_a_current = 0.0;
1525  for (size_t i = 0; i < m_kk; i++) {
1526  m_b_current += moleFractions_[i] * b_vec_Curr_[i];
1527  for (size_t j = 0; j < m_kk; j++) {
1528  m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
1529  }
1530  }
1531 }
1532 //====================================================================================================================
1533 void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
1534 {
1535  bCalc = 0.0;
1536  aCalc = 0.0;
1537  if (m_formTempParam == 1) {
1538  for (size_t i = 0; i < m_kk; i++) {
1539  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1540  for (size_t j = 0; j < m_kk; j++) {
1541  size_t counter = i * m_kk + j;
1542  doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
1543  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1544  }
1545  }
1546  } else {
1547  for (size_t i = 0; i < m_kk; i++) {
1548  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1549  for (size_t j = 0; j < m_kk; j++) {
1550  size_t counter = i * m_kk + j;
1551  doublereal a_vec_Curr = a_coeff_vec(0,counter);
1552  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1553  }
1554  }
1555  }
1556 }
1557 //====================================================================================================================
1558 doublereal RedlichKwongMFTP::da_dt() const
1559 {
1560 
1561  doublereal dadT = 0.0;
1562  if (m_formTempParam == 1) {
1563  for (size_t i = 0; i < m_kk; i++) {
1564  for (size_t j = 0; j < m_kk; j++) {
1565  size_t counter = i * m_kk + j;
1566  dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
1567  }
1568  }
1569  }
1570  return dadT;
1571 }
1572 //====================================================================================================================
1573 void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
1574  doublereal& pc, doublereal& tc, doublereal& vc) const
1575 {
1576  if (m_formTempParam != 0) {
1577  a = a0_coeff;
1578  }
1579  if (b <= 0.0) {
1580  tc = 1000000.;
1581  pc = 1.0E13;
1582  vc = omega_vc * GasConstant * tc / pc;
1583  return;
1584  }
1585  if (a <= 0.0) {
1586  tc = 0.0;
1587  pc = 0.0;
1588  vc = 2.0 * b;
1589  return;
1590  }
1591  double tmp = a * omega_b / (b * omega_a * GasConstant);
1592  double pp = 2./3.;
1593  doublereal sqrttc, f, dfdt, deltatc;
1594 
1595  if (m_formTempParam == 0) {
1596 
1597  tc = pow(tmp, pp);
1598  } else {
1599  tc = pow(tmp, pp);
1600  for (int j = 0; j < 10; j++) {
1601  sqrttc = sqrt(tc);
1602  f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT_coeff * tc - a0_coeff;
1603  dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT_coeff;
1604  deltatc = - f / dfdt;
1605  tc += deltatc;
1606  }
1607  if (deltatc > 0.1) {
1608  throw CanteraError("RedlichKwongMFTP::calcCriticalConditions", "didn't converge");
1609  }
1610  }
1611 
1612  pc = omega_b * GasConstant * tc / b;
1613  vc = omega_vc * GasConstant * tc / pc;
1614 }
1615 
1616 //====================================================================================================================
1617 // Solve the cubic equation of state
1618 /*
1619  * The R-K equation of state may be solved via the following formula
1620  *
1621  * V**3 - V**2(RT/P) - V(RTb/P - a/(P T**.5) + b*b) - (a b / (P T**.5)) = 0
1622  *
1623 
1624  * Returns the number of solutions found. If it only finds the liquid branch solution, it will return a -1 or a -2
1625  * instead of 1 or 2. If it returns 0, then there is an error.
1626  *
1627  */
1628 int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
1629  doublereal Vroot[3]) const
1630 {
1631  Vroot[0] = 0.0;
1632  Vroot[1] = 0.0;
1633  Vroot[2] = 0.0;
1634  int nTurningPoints;
1635  bool lotsOfNumError = false;
1636  doublereal Vturn[2];
1637  if (TKelvin <= 0.0) {
1638  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature");
1639  }
1640  /*
1641  * Derive the coefficients of the cubic polynomial to solve.
1642  */
1643  doublereal an = 1.0;
1644  doublereal bn = - GasConstant * TKelvin / pres;
1645  doublereal sqt = sqrt(TKelvin);
1646  doublereal cn = - (GasConstant * TKelvin * b / pres - a/(pres * sqt) + b * b);
1647  doublereal dn = - (a * b / (pres * sqt));
1648 
1649  double tmp = a * omega_b / (b * omega_a * GasConstant);
1650  double pp = 2./3.;
1651  double tc = pow(tmp, pp);
1652  double pc = omega_b * GasConstant * tc / b;
1653  double vc = omega_vc * GasConstant * tc / pc;
1654  // Derive the center of the cubic, x_N
1655  doublereal xN = - bn /(3 * an);
1656 
1657 
1658  // Derive the value of delta**2. This is a key quantity that determines the number of turning points
1659  doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
1660  doublereal delta = 0.0;
1661 
1662  // Calculate a couple of ratios
1663  doublereal ratio1 = 3.0 * an * cn / (bn * bn);
1664  doublereal ratio2 = pres * b / (GasConstant * TKelvin);
1665  if (fabs(ratio1) < 1.0E-7) {
1666  //printf("NicholsSolve(): Alternative solution (p = %g T = %g)\n", pres, TKelvin);
1667  doublereal ratio3 = a / (GasConstant * sqt) * pres / (GasConstant * TKelvin);
1668  if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
1669  doublereal z = 1.0;
1670  for (int i = 0; i < 10; i++) {
1671  doublereal znew = z / (z - ratio2) - ratio3 / (z + ratio1);
1672  doublereal deltaz = znew - z;
1673  z = znew;
1674  if (fabs(deltaz) < 1.0E-14) {
1675  break;
1676  }
1677  }
1678  doublereal v = z * GasConstant * TKelvin / pres;
1679  Vroot[0] = v;
1680  return 1;
1681  }
1682  }
1683 
1684 
1685  int nSolnValues;
1686  nTurningPoints = 2;
1687 
1688 #ifdef PRINTPV
1689  double V[100];
1690  int n = 0;
1691  for (int i = 0; i < 90; i++) {
1692  V[n++] = 0.030 + 0.005 * i;
1693  }
1694  double p1, presCalc;
1695  for (int i = 0; i < n; i++) {
1696  p1 = dpdVCalc(TKelvin, V[i], presCalc);
1697  printf(" %13.5g %13.5g %13.5g \n", V[i], presCalc , p1);
1698  }
1699 #endif
1700 
1701  double h2 = 4. * an * an * delta2 * delta2 * delta2;
1702  if (delta2 == 0.0) {
1703  nTurningPoints = 1;
1704  Vturn[0] = xN;
1705  Vturn[1] = xN;
1706  } else if (delta2 < 0.0) {
1707  nTurningPoints = 0;
1708  Vturn[0] = xN;
1709  Vturn[1] = xN;
1710  } else {
1711  delta = sqrt(delta2);
1712  Vturn[0] = xN - delta;
1713  Vturn[1] = xN + delta;
1714 #ifdef PRINTPV
1715  double presCalc;
1716  double p1 = dpdVCalc(TKelvin, Vturn[0], presCalc);
1717 
1718  double p2 = dpdVCalc(TKelvin, Vturn[1], presCalc);
1719 
1720  printf("p1 = %g p2 = %g \n", p1, p2);
1721  p1 = dpdVCalc(TKelvin, 0.9*Vturn[0], presCalc);
1722  printf("0.9 p1 = %g \n", p1);
1723 #endif
1724  }
1725 
1726  doublereal h = 2.0 * an * delta * delta2;
1727 
1728  doublereal yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn;
1729 
1730  doublereal desc = yN * yN - h2;
1731 
1732  if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
1733  if (desc != 0.0) {
1734  // this is for getting to other cases
1735  printf("NicholsSolve(): numerical issues\n");
1736  throw CanteraError("NicholsSolve()", "numerical issues");
1737  }
1738  desc = 0.0;
1739  }
1740 
1741  if (desc < 0.0) {
1742  nSolnValues = 3;
1743  } else if (desc == 0.0) {
1744  nSolnValues = 2;
1745  // We are here as p goes to zero.
1746  // double hleft = 3.0 * an * cn / (bn * bn);
1747  //double ynleft = 9.0 * an * cn / (2.0 * bn * bn) - 27.0 * an * an * dn / (2.0 * bn * bn * bn);
1748  //printf("hleft = %g , ynleft = %g\n", -3. / 2. * hleft, -ynleft);
1749  //double h2left = - 3 * hleft + 3 * hleft * hleft - hleft * hleft * hleft;
1750  //double y2left = - 2.0 * ynleft + ynleft * ynleft;
1751  //printf("h2left = %g , yn2left = %g\n", h2left, y2left);
1752 
1753  } else if (desc > 0.0) {
1754  nSolnValues = 1;
1755  }
1756 
1757  /*
1758  * One real root -> have to determine whether gas or liquid is the root
1759  */
1760  if (desc > 0.0) {
1761  doublereal tmpD = sqrt(desc);
1762  doublereal tmp1 = (- yN + tmpD) / (2.0 * an);
1763  doublereal sgn1 = 1.0;
1764  if (tmp1 < 0.0) {
1765  sgn1 = -1.0;
1766  tmp1 = -tmp1;
1767  }
1768  doublereal tmp2 = (- yN - tmpD) / (2.0 * an);
1769  doublereal sgn2 = 1.0;
1770  if (tmp2 < 0.0) {
1771  sgn2 = -1.0;
1772  tmp2 = -tmp2;
1773  }
1774  doublereal p1 = pow(tmp1, 1./3.);
1775  doublereal p2 = pow(tmp2, 1./3.);
1776 
1777  doublereal alpha = xN + sgn1 * p1 + sgn2 * p2;
1778  Vroot[0] = alpha;
1779  Vroot[1] = 0.0;
1780  Vroot[2] = 0.0;
1781 
1782  double tmp = an * Vroot[0] * Vroot[0] * Vroot[0] + bn * Vroot[0] * Vroot[0] + cn * Vroot[0] + dn;
1783  if (fabs(tmp) > 1.0E-4) {
1784  lotsOfNumError = true;
1785  }
1786 
1787  } else if (desc < 0.0) {
1788  doublereal tmp = - yN/h;
1789 
1790  doublereal val = acos(tmp);
1791  doublereal theta = val / 3.0;
1792 
1793  doublereal oo = 2. * Cantera::Pi / 3.;
1794  doublereal alpha = xN + 2. * delta * cos(theta);
1795 
1796  doublereal beta = xN + 2. * delta * cos(theta + oo);
1797 
1798  doublereal gamma = xN + 2. * delta * cos(theta + 2.0 * oo);
1799 
1800 
1801  Vroot[0] = beta;
1802  Vroot[1] = gamma;
1803  Vroot[2] = alpha;
1804 
1805  for (int i = 0; i < 3; i++) {
1806  double tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1807  if (fabs(tmp) > 1.0E-4) {
1808  lotsOfNumError = true;
1809  for (int j = 0; j < 3; j++) {
1810  if (j != i) {
1811  if (fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
1812  writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " +
1813  fp2str(pres) + "): WARNING roots have merged: " +
1814  fp2str(Vroot[i]) + ", " + fp2str(Vroot[j]));
1815  writelogendl();
1816  }
1817  }
1818  }
1819  }
1820  }
1821  } else if (desc == 0.0) {
1822  if (yN == 0.0 && h == 0.0) {
1823  Vroot[0] = xN;
1824  Vroot[1] = xN;
1825  Vroot[2] = xN;
1826  } else {
1827  // need to figure out whether delta is pos or neg
1828  if (yN > 0.0) {
1829  double tmp = pow(yN/(2*an), 1./3.);
1830  if (fabs(tmp - delta) > 1.0E-9) {
1831  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1832  }
1833  Vroot[1] = xN + delta;
1834  Vroot[0] = xN - 2.0*delta; // liquid phase root
1835  } else {
1836  double tmp = pow(yN/(2*an), 1./3.);
1837  if (fabs(tmp - delta) > 1.0E-9) {
1838  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1839  }
1840  delta = -delta;
1841  Vroot[0] = xN + delta;
1842  Vroot[1] = xN - 2.0*delta; // gas phase root
1843  }
1844  }
1845  for (int i = 0; i < 2; i++) {
1846  double tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1847  if (fabs(tmp) > 1.0E-4) {
1848  lotsOfNumError = true;
1849  }
1850  }
1851  }
1852 
1853  /*
1854  * Unfortunately, there is a heavy amount of roundoff error due to bad conditioning in this
1855  */
1856  double res, dresdV;
1857  for (int i = 0; i < nSolnValues; i++) {
1858  for (int n = 0; n < 20; n++) {
1859  res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1860  if (fabs(res) < 1.0E-14) {
1861  break;
1862  }
1863  dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
1864  double del = - res / dresdV;
1865 
1866  Vroot[i] += del;
1867  if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
1868  break;
1869  }
1870  double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1871  if (fabs(res2) < fabs(res)) {
1872  continue;
1873  } else {
1874  Vroot[i] -= del;
1875  Vroot[i] += 0.1 * del;
1876  }
1877  }
1878  if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
1879  writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " +
1880  fp2str(pres) + "): WARNING root didn't converge V = " + fp2str(Vroot[i]));
1881  writelogendl();
1882  }
1883  }
1884 
1885  if (nSolnValues == 1) {
1886  if (TKelvin > tc) {
1887  if (Vroot[0] < vc) {
1888  nSolnValues = -1;
1889  }
1890  } else {
1891  if (Vroot[0] < xN) {
1892  nSolnValues = -1;
1893  }
1894  }
1895 
1896  } else {
1897  if (nSolnValues == 2) {
1898  if (delta > 0.0) {
1899  nSolnValues = -2;
1900  }
1901  }
1902  }
1903  // writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + fp2str(pres) + "): finished");
1904  // writelogendl();
1905  return nSolnValues;
1906 }
1907 
1908 }