Cantera  2.1.2
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 
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 RedlichKwongMFTP::RedlichKwongMFTP() :
33  m_standardMixingRules(0),
34  m_formTempParam(0),
35  m_b_current(0.0),
36  m_a_current(0.0),
37  a_vec_Curr_(0),
38  b_vec_Curr_(0),
39  a_coeff_vec(0,0),
40  m_pc_Species(0),
41  m_tc_Species(0),
42  m_vc_Species(0),
43  NSolns_(0),
44  m_pp(0),
45  m_tmpV(0),
46  m_partialMolarVolumes(0),
47  dpdV_(0.0),
48  dpdT_(0.0),
49  dpdni_(0)
50 {
51  Vroot_[0] = 0.0;
52  Vroot_[1] = 0.0;
53  Vroot_[2] = 0.0;
54 }
55 
56 RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, std::string id_) :
58  m_standardMixingRules(0),
59  m_formTempParam(0),
60  m_b_current(0.0),
61  m_a_current(0.0),
62  a_vec_Curr_(0),
63  b_vec_Curr_(0),
64  a_coeff_vec(0,0),
65  m_pc_Species(0),
66  m_tc_Species(0),
67  m_vc_Species(0),
68  NSolns_(0),
69  m_pp(0),
70  m_tmpV(0),
71  m_partialMolarVolumes(0),
72  dpdV_(0.0),
73  dpdT_(0.0),
74  dpdni_(0)
75 {
76  Vroot_[0] = 0.0;
77  Vroot_[1] = 0.0;
78  Vroot_[2] = 0.0;
79  XML_Node* root = get_XML_File(infile);
80  if (id_ == "-") {
81  id_ = "";
82  }
83  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id_, root);
84  if (!xphase) {
85  throw CanteraError("newPhase",
86  "Couldn't find phase named \"" + id_ + "\" in file, " + infile);
87  }
88  importPhase(*xphase, this);
89 }
90 
91 RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id_) :
93  m_standardMixingRules(0),
94  m_formTempParam(0),
95  m_b_current(0.0),
96  m_a_current(0.0),
97  a_vec_Curr_(0),
98  b_vec_Curr_(0),
99  a_coeff_vec(0,0),
100  m_pc_Species(0),
101  m_tc_Species(0),
102  m_vc_Species(0),
103  NSolns_(0),
104  m_pp(0),
105  m_tmpV(0),
106  m_partialMolarVolumes(0),
107  dpdV_(0.0),
108  dpdT_(0.0),
109  dpdni_(0)
110 {
111  Vroot_[0] = 0.0;
112  Vroot_[1] = 0.0;
113  Vroot_[2] = 0.0;
114  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id_, &phaseRefRoot);
115  if (!xphase) {
116  throw CanteraError("RedlichKwongMFTP::RedlichKwongMFTP()","Couldn't find phase named \"" + id_ + "\" in XML node");
117  }
118  importPhase(*xphase, this);
119 }
120 
123  m_standardMixingRules(0),
124  m_formTempParam(0),
125  m_b_current(0.0),
126  m_a_current(0.0),
127  a_vec_Curr_(0),
128  b_vec_Curr_(0),
129  a_coeff_vec(0,0),
130  m_pc_Species(0),
131  m_tc_Species(0),
132  m_vc_Species(0),
133  NSolns_(0),
134  m_pp(0),
135  m_tmpV(0),
136  m_partialMolarVolumes(0),
137  dpdV_(0.0),
138  dpdT_(0.0),
139  dpdni_(0)
140 {
141  std::string infile = "co2_redlichkwong.xml";
142  std::string id_;
143  if (testProb == 1) {
144  infile = "co2_redlichkwong.xml";
145  id_ = "carbondioxide";
146  } else {
147  throw CanteraError("", "test prob = 1 only");
148  }
149  XML_Node* root = get_XML_File(infile);
150  if (id_ == "-") {
151  id_ = "";
152  }
153  XML_Node* xphase = get_XML_NameID("phase", std::string("#")+id_, root);
154  if (!xphase) {
155  throw CanteraError("newPhase", "Couldn't find phase named \"" + id_ + "\" in file, " + infile);
156  }
157  importPhase(*xphase, this);
158 }
159 
162  m_standardMixingRules(0),
163  m_formTempParam(0),
164  m_b_current(0.0),
165  m_a_current(0.0),
166  a_vec_Curr_(0),
167  b_vec_Curr_(0),
168  a_coeff_vec(0,0),
169  m_pc_Species(0),
170  m_tc_Species(0),
171  m_vc_Species(0),
172  NSolns_(0),
173  m_pp(0),
174  m_tmpV(0),
175  m_partialMolarVolumes(0),
176  dpdV_(0.0),
177  dpdT_(0.0),
178  dpdni_(0)
179 {
180  *this = b;
181 }
182 
185 {
186  if (&b != this) {
187  /*
188  * Mostly, this is a passthrough to the underlying
189  * assignment operator for the ThermoPhae parent object.
190  */
192  /*
193  * However, we have to handle data that we own.
194  */
199  a_vec_Curr_ = b.a_vec_Curr_;
200  b_vec_Curr_ = b.b_vec_Curr_;
201  a_coeff_vec = b.a_coeff_vec;
202 
203  m_pc_Species = b.m_pc_Species;
204  m_tc_Species = b.m_tc_Species;
205  m_vc_Species = b.m_vc_Species;
206  NSolns_ = b.NSolns_;
207  Vroot_[0] = b.Vroot_[0];
208  Vroot_[1] = b.Vroot_[1];
209  Vroot_[2] = b.Vroot_[2];
210  m_pp = b.m_pp;
211  m_tmpV = b.m_tmpV;
212  m_partialMolarVolumes = b.m_partialMolarVolumes;
213  dpdV_ = b.dpdV_;
214  dpdT_ = b.dpdT_;
215  dpdni_ = b.dpdni_;
216  }
217  return *this;
218 }
219 
221 {
222  return new RedlichKwongMFTP(*this);
223 }
224 
226 {
227  return cRedlichKwongMFTP;
228 }
229 
230 /*
231  * ------------Molar Thermodynamic Properties -------------------------
232  */
233 
235 {
237  doublereal rt = _RT();
238  doublereal h_ideal = rt * mean_X(DATA_PTR(m_h0_RT));
239  doublereal h_nonideal = hresid();
240  return h_ideal + h_nonideal;
241 }
242 
244 {
245  doublereal p0 = pressure();
246  doublereal md = molarDensity();
247  return enthalpy_mole() - p0 / md;
248 }
249 
251 {
253  doublereal sr_ideal = GasConstant * (mean_X(DATA_PTR(m_s0_R))
254  - sum_xlogx() - std::log(pressure()/m_spthermo->refPressure()));
255  doublereal sr_nonideal = sresid();
256  return sr_ideal + sr_nonideal;
257 }
258 
260 {
261  return enthalpy_mole() - temperature() * entropy_mole();
262 }
263 
264 doublereal RedlichKwongMFTP::cp_mole() const
265 {
267  doublereal TKelvin = temperature();
268  doublereal sqt = sqrt(TKelvin);
269  doublereal mv = molarVolume();
270  doublereal vpb = mv + m_b_current;
272  doublereal cpref = GasConstant * mean_X(DATA_PTR(m_cp0_R));
273  doublereal dadt = da_dt();
274  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
275  doublereal dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
276  +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
277  return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
278 }
279 
280 doublereal RedlichKwongMFTP::cv_mole() const
281 {
282  throw CanteraError("", "unimplemented");
283  return cp_mole() - GasConstant;
284 }
285 
286 doublereal RedlichKwongMFTP::pressure() const
287 {
288 #ifdef DEBUG_MODE
290 
291  // Get a copy of the private variables stored in the State object
292  double rho = density();
293  doublereal T = temperature();
294  doublereal mmw = meanMolecularWeight();
295  double molarV = mmw / rho;
296 
297  double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
298 
299  if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) {
300  throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe");
301  }
302 #endif
303  return m_Pcurrent;
304 }
305 
307 {
308  /*
309  * Calculate the molarVolume of the solution (m**3 kmol-1)
310  */
311 
312  const doublereal* const dtmp = moleFractdivMMW();
314  double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp);
315  /*
316  * Set the density in the parent State object directly,
317  * by calling the Phase::setDensity() function.
318  */
319  double dens = 1.0/invDens;
320  Phase::setDensity(dens);
321 }
322 
323 void RedlichKwongMFTP::setTemperature(const doublereal temp)
324 {
325  Phase::setTemperature(temp);
327  updateAB();
328 }
329 
330 void RedlichKwongMFTP::setMassFractions(const doublereal* const x)
331 {
333  updateAB();
334 }
335 
336 void RedlichKwongMFTP::setMassFractions_NoNorm(const doublereal* const x)
337 {
339  updateAB();
340 }
341 
342 void RedlichKwongMFTP::setMoleFractions(const doublereal* const x)
343 {
345  updateAB();
346 }
347 
348 void RedlichKwongMFTP::setMoleFractions_NoNorm(const doublereal* const x)
349 {
351  updateAB();
352 }
353 
354 void RedlichKwongMFTP::setConcentrations(const doublereal* const c)
355 {
357  updateAB();
358 }
359 
361 {
362  throw CanteraError("RedlichKwongMFTP::isothermalCompressibility() ",
363  "not implemented");
364  return 0.0;
365 }
366 
368 {
369  getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes));
370  for (size_t k = 0; k < m_kk; k++) {
371  c[k] = moleFraction(k) / m_partialMolarVolumes[k];
372  }
373 }
374 
375 doublereal RedlichKwongMFTP::standardConcentration(size_t k) const
376 {
378  return 1.0 / m_tmpV[k];
379 }
380 
381 doublereal RedlichKwongMFTP::logStandardConc(size_t k) const
382 {
383  double c = standardConcentration(k);
384  return std::log(c);
385 }
386 
387 void RedlichKwongMFTP::getUnitsStandardConc(double* uA, int, int sizeUA) const
388 {
389  //int eos = eosType();
390 
391  for (int i = 0; i < sizeUA; i++) {
392  if (i == 0) {
393  uA[0] = 1.0;
394  }
395  if (i == 1) {
396  uA[1] = -static_cast<int>(nDim());
397  }
398  if (i == 2) {
399  uA[2] = 0.0;
400  }
401  if (i == 3) {
402  uA[3] = 0.0;
403  }
404  if (i == 4) {
405  uA[4] = 0.0;
406  }
407  if (i == 5) {
408  uA[5] = 0.0;
409  }
410  }
411 
412 }
413 
415 {
416  doublereal TKelvin = temperature();
417  doublereal rt = TKelvin * GasConstant;
418  doublereal mv = molarVolume();
419  doublereal sqt = sqrt(TKelvin);
420  doublereal vpb = mv + m_b_current;
421  doublereal vmb = mv - m_b_current;
422 
423  for (size_t k = 0; k < m_kk; k++) {
424  m_pp[k] = 0.0;
425  for (size_t i = 0; i < m_kk; i++) {
426  size_t counter = k + m_kk*i;
427  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
428  }
429  }
430  doublereal pres = pressure();
431 
432  for (size_t k = 0; k < m_kk; k++) {
433  ac[k] = (- rt * log(pres * mv / rt)
434  + rt * log(mv / vmb)
435  + rt * b_vec_Curr_[k] / vmb
436  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
437  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
438  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
439  );
440  }
441  for (size_t k = 0; k < m_kk; k++) {
442  ac[k] = exp(ac[k]/rt);
443  }
444 }
445 
446 /*
447  * ---- Partial Molar Properties of the Solution -----------------
448  */
449 
450 void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
451 {
452  getChemPotentials(muRT);
453  doublereal invRT = 1.0 / _RT();
454  for (size_t k = 0; k < m_kk; k++) {
455  muRT[k] *= invRT;
456  }
457 }
458 
459 void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
460 {
461  getGibbs_ref(mu);
462  doublereal xx;
463  doublereal rt = temperature() * GasConstant;
464  for (size_t k = 0; k < m_kk; k++) {
465  xx = std::max(SmallNumber, moleFraction(k));
466  mu[k] += rt*(log(xx));
467  }
468 
469  doublereal TKelvin = temperature();
470  doublereal mv = molarVolume();
471  doublereal sqt = sqrt(TKelvin);
472  doublereal vpb = mv + m_b_current;
473  doublereal vmb = mv - m_b_current;
474 
475  for (size_t k = 0; k < m_kk; k++) {
476  m_pp[k] = 0.0;
477  for (size_t i = 0; i < m_kk; i++) {
478  size_t counter = k + m_kk*i;
479  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
480  }
481  }
482  doublereal pres = pressure();
483  doublereal refP = refPressure();
484 
485  for (size_t k = 0; k < m_kk; k++) {
486  mu[k] += (rt * log(pres/refP) - rt * log(pres * mv / rt)
487  + rt * log(mv / vmb)
488  + rt * b_vec_Curr_[k] / vmb
489  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
490  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
491  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
492  );
493  }
494 }
495 
497 {
498  /*
499  * First we get the reference state contributions
500  */
501  getEnthalpy_RT_ref(hbar);
502  doublereal rt = GasConstant * temperature();
503  scale(hbar, hbar+m_kk, hbar, rt);
504 
505  /*
506  * We calculate dpdni_
507  */
508  doublereal TKelvin = temperature();
509  doublereal mv = molarVolume();
510  doublereal sqt = sqrt(TKelvin);
511 
512  doublereal vpb = mv + m_b_current;
513  doublereal vmb = mv - m_b_current;
514 
515  for (size_t k = 0; k < m_kk; k++) {
516  m_pp[k] = 0.0;
517  for (size_t i = 0; i < m_kk; i++) {
518  size_t counter = k + m_kk*i;
519  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
520  }
521  }
522 
523 
524 
525  for (size_t k = 0; k < m_kk; k++) {
526  dpdni_[k] = rt/vmb + rt * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
527  + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
528  }
529  doublereal dadt = da_dt();
530  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
531 
532  for (size_t k = 0; k < m_kk; k++) {
533  m_tmpV[k] = 0.0;
534  for (size_t i = 0; i < m_kk; i++) {
535  size_t counter = k + m_kk*i;
536  m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
537  }
538  }
539 
541  doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_;
542 
543  for (size_t k = 0; k < m_kk; k++) {
544  double hE_v = (mv * dpdni_[k] - rt - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
545  + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
546  + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
547  hbar[k] = hbar[k] + hE_v;
548 
549 
550  hbar[k] -= fac2 * dpdni_[k];
551  }
552 
553 }
554 
555 void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
556 {
557  getEntropy_R_ref(sbar);
558  doublereal r = GasConstant;
559  scale(sbar, sbar+m_kk, sbar, r);
560  doublereal TKelvin = temperature();
561  doublereal sqt = sqrt(TKelvin);
562  doublereal mv = molarVolume();
563  doublereal refP = refPressure();
564 
565  for (size_t k = 0; k < m_kk; k++) {
566  doublereal xx = std::max(SmallNumber, moleFraction(k));
567  sbar[k] += r * (- log(xx));
568  }
569 
570  for (size_t k = 0; k < m_kk; k++) {
571  m_pp[k] = 0.0;
572  for (size_t i = 0; i < m_kk; i++) {
573  size_t counter = k + m_kk*i;
574  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
575  }
576  }
577 
578  for (size_t k = 0; k < m_kk; k++) {
579  m_tmpV[k] = 0.0;
580  for (size_t i = 0; i < m_kk; i++) {
581  size_t counter = k + m_kk*i;
582  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
583  }
584  }
585 
586 
587  doublereal dadt = da_dt();
588  doublereal fac = dadt - m_a_current / (2.0 * TKelvin);
589  doublereal vmb = mv - m_b_current;
590  doublereal vpb = mv + m_b_current;
591 
592 
593  for (size_t k = 0; k < m_kk; k++) {
594  sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
595  + GasConstant
596  + GasConstant * log(mv/vmb)
597  + GasConstant * b_vec_Curr_[k]/vmb
598  + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
599  - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
600  + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
601  - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
602  ) ;
603  }
604 
606  getPartialMolarVolumes(DATA_PTR(m_partialMolarVolumes));
607  for (size_t k = 0; k < m_kk; k++) {
608  sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
609  }
610 }
611 
613 {
614  getIntEnergy_RT(ubar);
615  doublereal rt = GasConstant * temperature();
616  scale(ubar, ubar+m_kk, ubar, rt);
617 }
618 
619 void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const
620 {
621  getCp_R(cpbar);
622  doublereal r = GasConstant;
623  scale(cpbar, cpbar+m_kk, cpbar, r);
624 }
625 
626 void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
627 {
628  // getStandardVolumes(vbar);
629 
630 
631  for (size_t k = 0; k < m_kk; k++) {
632  m_pp[k] = 0.0;
633  for (size_t i = 0; i < m_kk; i++) {
634  size_t counter = k + m_kk*i;
635  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
636  }
637  }
638 
639  for (size_t k = 0; k < m_kk; k++) {
640  m_tmpV[k] = 0.0;
641  for (size_t i = 0; i < m_kk; i++) {
642  size_t counter = k + m_kk*i;
643  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
644  }
645  }
646 
647  doublereal TKelvin = temperature();
648  doublereal sqt = sqrt(TKelvin);
649  doublereal mv = molarVolume();
650 
651  doublereal rt = GasConstant * TKelvin;
652 
653  doublereal vmb = mv - m_b_current;
654  doublereal vpb = mv + m_b_current;
655 
656  for (size_t k = 0; k < m_kk; k++) {
657 
658  doublereal num = (rt + rt * m_b_current/ vmb + rt * b_vec_Curr_[k] / vmb
659  + rt * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
660  - 2.0 * m_pp[k] / (sqt * vpb)
661  + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
662  );
663 
664  doublereal denom = (m_Pcurrent + rt * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
665  );
666 
667  vbar[k] = num / denom;
668  }
669 
670 }
671 
673 {
674  double pc, tc, vc;
675  double a0 = 0.0;
676  double aT = 0.0;
677  for (size_t i = 0; i < m_kk; i++) {
678  for (size_t j = 0; j <m_kk; j++) {
679  size_t counter = i + m_kk * j;
680  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
681  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
682  }
683  }
684  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
685  return tc;
686 }
687 
689 {
690  double pc, tc, vc;
691  double a0 = 0.0;
692  double aT = 0.0;
693  for (size_t i = 0; i < m_kk; i++) {
694  for (size_t j = 0; j <m_kk; j++) {
695  size_t counter = i + m_kk * j;
696  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
697  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
698  }
699  }
700  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
701 
702  return pc;
703 }
704 
706 {
707  double pc, tc, vc;
708  double a0 = 0.0;
709  double aT = 0.0;
710  for (size_t i = 0; i < m_kk; i++) {
711  for (size_t j = 0; j <m_kk; j++) {
712  size_t counter = i + m_kk * j;
713  a0 += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(0, counter);
714  aT += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(1, counter);
715  }
716  }
717  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
718 
719  double mmw = meanMolecularWeight();
720  return mmw / vc;
721 }
722 
724 {
725  initLengths();
727 }
728 
729 void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
730 {
731  double tmp, tmp2;
733 
735 
736 
737  /*
738  * Within the method, we protect against inf results if the
739  * exponent is too high.
740  *
741  * If it is too low, we set
742  * the partial pressure to zero. This capability is needed
743  * by the elemental potential method.
744  */
745  doublereal pres = 0.0;
746  double m_p0 = refPressure();
747  for (size_t k = 0; k < m_kk; k++) {
748  tmp = -m_tmpV[k] + mu_RT[k];
749  if (tmp < -600.) {
750  m_pp[k] = 0.0;
751  } else if (tmp > 500.0) {
752  tmp2 = tmp / 500.;
753  tmp2 *= tmp2;
754  m_pp[k] = m_p0 * exp(500.) * tmp2;
755  } else {
756  m_pp[k] = m_p0 * exp(tmp);
757  }
758  pres += m_pp[k];
759  }
760  // set state
761  setState_PX(pres, &m_pp[0]);
762 }
763 
765 {
766 
767 
768  a_vec_Curr_.resize(m_kk * m_kk, 0.0);
769  b_vec_Curr_.resize(m_kk, 0.0);
770 
771  a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
772 
773 
774  m_pc_Species.resize(m_kk, 0.0);
775  m_tc_Species.resize(m_kk, 0.0);
776  m_vc_Species.resize(m_kk, 0.0);
777 
778 
779  m_pp.resize(m_kk, 0.0);
780  m_tmpV.resize(m_kk, 0.0);
781  m_partialMolarVolumes.resize(m_kk, 0.0);
782  dpdni_.resize(m_kk, 0.0);
783 }
784 
785 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
786 {
788 
789  /*
790  * Check the model parameter for the Redlich-Kwong equation of state
791  * two are allowed
792  * RedlichKwong mixture of species, each of which are RK fluids
793  * RedlichKwongMFTP mixture of species with cross term coefficients
794  */
795  if (phaseNode.hasChild("thermo")) {
796  XML_Node& thermoNode = phaseNode.child("thermo");
797  std::string model = thermoNode["model"];
798  if (model == "RedlichKwong") {
800  } else if (model == "RedlichKwongMFTP") {
802  } else {
803  throw CanteraError("RedlichKwongMFTP::initThermoXML",
804  "Unknown thermo model : " + model);
805  }
806 
807 
808  /*
809  * Go get all of the coefficients and factors in the
810  * activityCoefficients XML block
811  */
812  XML_Node* acNodePtr = 0;
813  if (thermoNode.hasChild("activityCoefficients")) {
814  XML_Node& acNode = thermoNode.child("activityCoefficients");
815  acNodePtr = &acNode;
816  size_t nC = acNode.nChildren();
817 
818  /*
819  * Loop through the children getting multiple instances of
820  * parameters
821  */
822  for (size_t i = 0; i < nC; i++) {
823  XML_Node& xmlACChild = acNodePtr->child(i);
824  string stemp = xmlACChild.name();
825  string nodeName = lowercase(stemp);
826  /*
827  * Process a binary salt field, or any of the other XML fields
828  * that make up the Pitzer Database. Entries will be ignored
829  * if any of the species in the entry isn't in the solution.
830  */
831  if (nodeName == "purefluidparameters") {
832  readXMLPureFluid(xmlACChild);
833  }
834  }
835  if (m_standardMixingRules == 1) {
837  }
838  /*
839  * Loop through the children getting multiple instances of
840  * parameters
841  */
842  for (size_t i = 0; i < nC; i++) {
843  XML_Node& xmlACChild = acNodePtr->child(i);
844  string stemp = xmlACChild.name();
845  string nodeName = lowercase(stemp);
846  /*
847  * Process a binary salt field, or any of the other XML fields
848  * that make up the Pitzer Database. Entries will be ignored
849  * if any of the species in the entry isn't in the solution.
850  */
851  if (nodeName == "crossfluidparameters") {
852  readXMLCrossFluid(xmlACChild);
853  }
854  }
855 
856  }
857  }
858 
859  for (size_t i = 0; i < m_kk; i++) {
860  double a0coeff = a_coeff_vec(0, i*m_kk + i);
861  double aTcoeff = a_coeff_vec(1, i*m_kk + i);
862  double ai = a0coeff + aTcoeff * 500.;
863  double bi = b_vec_Curr_[i];
864  calcCriticalConditions(ai, bi, a0coeff, aTcoeff, m_pc_Species[i], m_tc_Species[i], m_vc_Species[i]);
865  }
866 
867  MixtureFugacityTP::initThermoXML(phaseNode, id);
868 }
869 
871 {
872  vector_fp vParams;
873  string xname = pureFluidParam.name();
874  if (xname != "pureFluidParameters") {
875  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
876  "Incorrect name for processing this routine: " + xname);
877  }
878 
879  /*
880  * Read the species
881  * Find the index of the species in the current phase. It's not an error to not find the species
882  */
883  string iName = pureFluidParam.attrib("species");
884  if (iName == "") {
885  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute");
886  }
887  size_t iSpecies = speciesIndex(iName);
888  if (iSpecies == npos) {
889  return;
890  }
891  size_t counter = iSpecies + m_kk * iSpecies;
892  size_t nParamsExpected, nParamsFound;
893  size_t num = pureFluidParam.nChildren();
894  for (size_t iChild = 0; iChild < num; iChild++) {
895  XML_Node& xmlChild = pureFluidParam.child(iChild);
896  string stemp = xmlChild.name();
897  string nodeName = lowercase(stemp);
898 
899  if (nodeName == "a_coeff") {
900  string iModel = lowercase(xmlChild.attrib("model"));
901  if (iModel == "constant") {
902  nParamsExpected = 1;
903  } else if (iModel == "linear_a") {
904  nParamsExpected = 2;
905  if (m_formTempParam == 0) {
906  m_formTempParam = 1;
907  }
908  } else {
909  throw CanteraError("", "unknown model");
910  }
911 
912  ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
913  nParamsFound = vParams.size();
914  if (nParamsFound != nParamsExpected) {
915  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for a_coeff" + iName + ")",
916  "wrong number of params found");
917  }
918 
919  for (size_t i = 0; i < nParamsFound; i++) {
920  a_coeff_vec(i, counter) = vParams[i];
921  }
922  } else if (nodeName == "b_coeff") {
923  ctml::getFloatArray(xmlChild, vParams, true, "m3/kmol", "b_coeff");
924  nParamsFound = vParams.size();
925  if (nParamsFound != 1) {
926  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for b_coeff" + iName + ")",
927  "wrong number of params found");
928  }
929  b_vec_Curr_[iSpecies] = vParams[0];
930  }
931  }
932 }
933 
935 {
936  int nParam = 2;
937  for (size_t i = 0; i < m_kk; i++) {
938  size_t icounter = i + m_kk * i;
939  for (size_t j = 0; j < m_kk; j++) {
940  if (i != j) {
941  size_t counter = i + m_kk * j;
942  size_t jcounter = j + m_kk * j;
943  for (int n = 0; n < nParam; n++) {
944  a_coeff_vec(n, counter) = sqrt(a_coeff_vec(n, icounter) * a_coeff_vec(n, jcounter));
945  }
946  }
947  }
948  }
949 }
950 
952 {
953  vector_fp vParams;
954  string xname = CrossFluidParam.name();
955  if (xname != "crossFluidParameters") {
956  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
957  "Incorrect name for processing this routine: " + xname);
958  }
959 
960  /*
961  * Read the species
962  * Find the index of the species in the current phase. It's not an error to not find the species
963  */
964  string iName = CrossFluidParam.attrib("species1");
965  if (iName == "") {
966  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute");
967  }
968  size_t iSpecies = speciesIndex(iName);
969  if (iSpecies == npos) {
970  return;
971  }
972  string jName = CrossFluidParam.attrib("species2");
973  if (iName == "") {
974  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species2 attribute");
975  }
976  size_t jSpecies = speciesIndex(jName);
977  if (jSpecies == npos) {
978  return;
979  }
980 
981  size_t counter = iSpecies + m_kk * jSpecies;
982  size_t counter0 = jSpecies + m_kk * iSpecies;
983  size_t nParamsExpected, nParamsFound;
984  size_t num = CrossFluidParam.nChildren();
985  for (size_t iChild = 0; iChild < num; iChild++) {
986  XML_Node& xmlChild = CrossFluidParam.child(iChild);
987  string stemp = xmlChild.name();
988  string nodeName = lowercase(stemp);
989 
990  if (nodeName == "a_coeff") {
991  string iModel = lowercase(xmlChild.attrib("model"));
992  if (iModel == "constant") {
993  nParamsExpected = 1;
994  } else if (iModel == "linear_a") {
995  nParamsExpected = 2;
996  if (m_formTempParam == 0) {
997  m_formTempParam = 1;
998  }
999  } else {
1000  throw CanteraError("", "unknown model");
1001  }
1002 
1003  ctml::getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
1004  nParamsFound = vParams.size();
1005  if (nParamsFound != nParamsExpected) {
1006  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid(for a_coeff" + iName + ")",
1007  "wrong number of params found");
1008  }
1009 
1010  for (size_t i = 0; i < nParamsFound; i++) {
1011  a_coeff_vec(i, counter) = vParams[i];
1012  a_coeff_vec(i, counter0) = vParams[i];
1013  }
1014  }
1015  }
1016 }
1017 
1019 {
1021  std::string model = thermoNode["model"];
1022 }
1023 
1024 doublereal RedlichKwongMFTP::sresid() const
1025 {
1026  // note this agrees with tpx
1027  doublereal rho = density();
1028  doublereal mmw = meanMolecularWeight();
1029  doublereal molarV = mmw / rho;
1030  double hh = m_b_current / molarV;
1031  doublereal zz = z();
1032  doublereal dadt = da_dt();
1033  doublereal T = temperature();
1034  doublereal sqT = sqrt(T);
1035  doublereal fac = dadt - m_a_current / (2.0 * T);
1036  double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
1037  return GasConstant * sresid_mol_R;
1038 }
1039 
1040 doublereal RedlichKwongMFTP::hresid() const
1041 {
1042  // note this agrees with tpx
1043  doublereal rho = density();
1044  doublereal mmw = meanMolecularWeight();
1045  doublereal molarV = mmw / rho;
1046  double hh = m_b_current / molarV;
1047  doublereal zz = z();
1048  doublereal dadt = da_dt();
1049  doublereal T = temperature();
1050  doublereal sqT = sqrt(T);
1051  doublereal fac = T * dadt - 3.0 *m_a_current / (2.0);
1052  return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
1053 }
1054 
1055 doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
1056 {
1057  double v = m_b_current * 1.1;
1058  double atmp;
1059  double btmp;
1060  calculateAB(TKelvin, atmp, btmp);
1061 
1062  doublereal pres = presGuess;
1063  double pp = psatEst(TKelvin);
1064  if (pres < pp) {
1065  pres = pp;
1066  }
1067  double Vroot[3];
1068 
1069  bool foundLiq = false;
1070  int m = 0;
1071  do {
1072 
1073  int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot);
1074 
1075  // printf("nsol = %d\n", nsol);
1076  // printf("liquidVolEst start: T = %g , p = %g, a = %g, b = %g\n", TKelvin, pres, m_a_current, m_b_current);
1077 
1078  if (nsol == 1 || nsol == 2) {
1079  double pc = critPressure();
1080  if (pres > pc) {
1081  foundLiq = true;
1082  }
1083  pres *= 1.04;
1084 
1085  } else {
1086  foundLiq = true;
1087  }
1088  } while ((m < 100) && (!foundLiq));
1089 
1090  if (foundLiq) {
1091  v = Vroot[0];
1092  presGuess = pres;
1093  } else {
1094  v = -1.0;
1095  }
1096  //printf (" RedlichKwongMFTP::liquidVolEst %g %g converged in %d its\n", TKelvin, pres, i);
1097  return v;
1098 }
1099 
1100 doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
1101 {
1102 
1103  /*
1104  * It's necessary to set the temperature so that m_a_current is set correctly.
1105  */
1106  setTemperature(TKelvin);
1107  double tcrit = critTemperature();
1108  doublereal mmw = meanMolecularWeight();
1109  if (rhoguess == -1.0) {
1110  if (phaseRequested != FLUID_GAS) {
1111  if (TKelvin > tcrit) {
1112  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1113  } else {
1114  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
1115  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1116  } else if (phaseRequested >= FLUID_LIQUID_0) {
1117  double lqvol = liquidVolEst(TKelvin, presPa);
1118  rhoguess = mmw / lqvol;
1119  }
1120  }
1121  } else {
1122  /*
1123  * Assume the Gas phase initial guess, if nothing is
1124  * specified to the routine
1125  */
1126  rhoguess = presPa * mmw / (GasConstant * TKelvin);
1127  }
1128 
1129  }
1130 
1131 
1132  doublereal volguess = mmw / rhoguess;
1133  NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
1134 
1135  doublereal molarVolLast = Vroot_[0];
1136  if (NSolns_ >= 2) {
1137  if (phaseRequested >= FLUID_LIQUID_0) {
1138  molarVolLast = Vroot_[0];
1139  } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
1140  molarVolLast = Vroot_[2];
1141  } else {
1142  if (volguess > Vroot_[1]) {
1143  molarVolLast = Vroot_[2];
1144  } else {
1145  molarVolLast = Vroot_[0];
1146  }
1147  }
1148  } else if (NSolns_ == 1) {
1149  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
1150  molarVolLast = Vroot_[0];
1151  } else {
1152  //molarVolLast = Vroot_[0];
1153  //printf("DensityCalc(): Possible problem encountered\n");
1154  return -2.0;
1155  }
1156  } else if (NSolns_ == -1) {
1157  if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
1158  molarVolLast = Vroot_[0];
1159  } else if (TKelvin > tcrit) {
1160  molarVolLast = Vroot_[0];
1161  } else {
1162  // molarVolLast = Vroot_[0];
1163  // printf("DensityCalc(): Possible problem encountered\n");
1164  return -2.0;
1165  }
1166  } else {
1167  molarVolLast = Vroot_[0];
1168  //printf("DensityCalc(): Possible problem encountered\n");
1169  return -1.0;
1170  }
1171  return mmw / molarVolLast;
1172 }
1173 
1175 {
1176  if (NSolns_ != 3) {
1177  return critDensity();
1178  }
1179  double vmax = Vroot_[1];
1180  double vmin = Vroot_[0];
1181  RootFind rf(fdpdv_);
1182  rf.setPrintLvl(10);
1183  rf.setTol(1.0E-5, 1.0E-10);
1185 
1186  double vbest = 0.5 * (Vroot_[0]+Vroot_[1]);
1187  double funcNeeded = 0.0;
1188 
1189  int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
1190  if (status != ROOTFIND_SUCCESS) {
1191  throw CanteraError(" RedlichKwongMFTP::densSpinodalLiquid() ", "didn't converge");
1192  }
1193  doublereal mmw = meanMolecularWeight();
1194  return mmw / vbest;
1195 }
1196 
1198 {
1199  if (NSolns_ != 3) {
1200  return critDensity();
1201  }
1202  double vmax = Vroot_[2];
1203  double vmin = Vroot_[1];
1204  RootFind rf(fdpdv_);
1205  rf.setPrintLvl(10);
1206  rf.setTol(1.0E-5, 1.0E-10);
1208 
1209  double vbest = 0.5 * (Vroot_[1]+Vroot_[2]);
1210  double funcNeeded = 0.0;
1211 
1212  int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
1213  if (status != ROOTFIND_SUCCESS) {
1214  throw CanteraError(" RedlichKwongMFTP::densSpinodalGas() ", "didn't converge");
1215  }
1216  doublereal mmw = meanMolecularWeight();
1217  return mmw / vbest;
1218 }
1219 
1220 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
1221 {
1222  doublereal sqt = sqrt(TKelvin);
1223  double pres = GasConstant * TKelvin / (molarVol - m_b_current)
1224  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
1225  return pres;
1226 }
1227 
1228 doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
1229 {
1230  doublereal sqt = sqrt(TKelvin);
1231  presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
1232  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
1233 
1234  doublereal vpb = molarVol + m_b_current;
1235  doublereal vmb = molarVol - m_b_current;
1236  doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb)
1237  + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
1238  return dpdv;
1239 }
1240 
1242 {
1243  doublereal TKelvin = temperature();
1244  doublereal mv = molarVolume();
1245  doublereal pres;
1246 
1247  dpdV_ = dpdVCalc(TKelvin, mv, pres);
1248 
1249  doublereal sqt = sqrt(TKelvin);
1250  doublereal vpb = mv + m_b_current;
1251  doublereal vmb = mv - m_b_current;
1252  doublereal dadt = da_dt();
1253  doublereal fac = dadt - m_a_current/(2.0 * TKelvin);
1254 
1255  dpdT_ = (GasConstant / (vmb) - fac / (sqt * mv * vpb));
1256 }
1257 
1258 void RedlichKwongMFTP::updateMixingExpressions()
1259 {
1260  updateAB();
1261 }
1262 
1264 {
1265  double temp = temperature();
1266  if (m_formTempParam == 1) {
1267  for (size_t i = 0; i < m_kk; i++) {
1268  for (size_t j = 0; j < m_kk; j++) {
1269  size_t counter = i * m_kk + j;
1270  a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
1271  }
1272  }
1273  }
1274 
1275  m_b_current = 0.0;
1276  m_a_current = 0.0;
1277  for (size_t i = 0; i < m_kk; i++) {
1278  m_b_current += moleFractions_[i] * b_vec_Curr_[i];
1279  for (size_t j = 0; j < m_kk; j++) {
1280  m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
1281  }
1282  }
1283 }
1284 
1285 void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
1286 {
1287  bCalc = 0.0;
1288  aCalc = 0.0;
1289  if (m_formTempParam == 1) {
1290  for (size_t i = 0; i < m_kk; i++) {
1291  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1292  for (size_t j = 0; j < m_kk; j++) {
1293  size_t counter = i * m_kk + j;
1294  doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
1295  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1296  }
1297  }
1298  } else {
1299  for (size_t i = 0; i < m_kk; i++) {
1300  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1301  for (size_t j = 0; j < m_kk; j++) {
1302  size_t counter = i * m_kk + j;
1303  doublereal a_vec_Curr = a_coeff_vec(0,counter);
1304  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1305  }
1306  }
1307  }
1308 }
1309 
1310 doublereal RedlichKwongMFTP::da_dt() const
1311 {
1312 
1313  doublereal dadT = 0.0;
1314  if (m_formTempParam == 1) {
1315  for (size_t i = 0; i < m_kk; i++) {
1316  for (size_t j = 0; j < m_kk; j++) {
1317  size_t counter = i * m_kk + j;
1318  dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
1319  }
1320  }
1321  }
1322  return dadT;
1323 }
1324 
1325 void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
1326  doublereal& pc, doublereal& tc, doublereal& vc) const
1327 {
1328  if (m_formTempParam != 0) {
1329  a = a0_coeff;
1330  }
1331  if (b <= 0.0) {
1332  tc = 1000000.;
1333  pc = 1.0E13;
1334  vc = omega_vc * GasConstant * tc / pc;
1335  return;
1336  }
1337  if (a <= 0.0) {
1338  tc = 0.0;
1339  pc = 0.0;
1340  vc = 2.0 * b;
1341  return;
1342  }
1343  double tmp = a * omega_b / (b * omega_a * GasConstant);
1344  double pp = 2./3.;
1345  doublereal sqrttc, f, dfdt, deltatc;
1346 
1347  if (m_formTempParam == 0) {
1348 
1349  tc = pow(tmp, pp);
1350  } else {
1351  tc = pow(tmp, pp);
1352  for (int j = 0; j < 10; j++) {
1353  sqrttc = sqrt(tc);
1354  f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT_coeff * tc - a0_coeff;
1355  dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT_coeff;
1356  deltatc = - f / dfdt;
1357  tc += deltatc;
1358  }
1359  if (deltatc > 0.1) {
1360  throw CanteraError("RedlichKwongMFTP::calcCriticalConditions", "didn't converge");
1361  }
1362  }
1363 
1364  pc = omega_b * GasConstant * tc / b;
1365  vc = omega_vc * GasConstant * tc / pc;
1366 }
1367 
1368 int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
1369  doublereal Vroot[3]) const
1370 {
1371  Vroot[0] = 0.0;
1372  Vroot[1] = 0.0;
1373  Vroot[2] = 0.0;
1374  int nTurningPoints;
1375  bool lotsOfNumError = false;
1376  doublereal Vturn[2];
1377  if (TKelvin <= 0.0) {
1378  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature");
1379  }
1380  /*
1381  * Derive the coefficients of the cubic polynomial to solve.
1382  */
1383  doublereal an = 1.0;
1384  doublereal bn = - GasConstant * TKelvin / pres;
1385  doublereal sqt = sqrt(TKelvin);
1386  doublereal cn = - (GasConstant * TKelvin * b / pres - a/(pres * sqt) + b * b);
1387  doublereal dn = - (a * b / (pres * sqt));
1388 
1389  double tmp = a * omega_b / (b * omega_a * GasConstant);
1390  double pp = 2./3.;
1391  double tc = pow(tmp, pp);
1392  double pc = omega_b * GasConstant * tc / b;
1393  double vc = omega_vc * GasConstant * tc / pc;
1394  // Derive the center of the cubic, x_N
1395  doublereal xN = - bn /(3 * an);
1396 
1397 
1398  // Derive the value of delta**2. This is a key quantity that determines the number of turning points
1399  doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
1400  doublereal delta = 0.0;
1401 
1402  // Calculate a couple of ratios
1403  doublereal ratio1 = 3.0 * an * cn / (bn * bn);
1404  doublereal ratio2 = pres * b / (GasConstant * TKelvin);
1405  if (fabs(ratio1) < 1.0E-7) {
1406  //printf("NicholsSolve(): Alternative solution (p = %g T = %g)\n", pres, TKelvin);
1407  doublereal ratio3 = a / (GasConstant * sqt) * pres / (GasConstant * TKelvin);
1408  if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
1409  doublereal zz = 1.0;
1410  for (int i = 0; i < 10; i++) {
1411  doublereal znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
1412  doublereal deltaz = znew - zz;
1413  zz = znew;
1414  if (fabs(deltaz) < 1.0E-14) {
1415  break;
1416  }
1417  }
1418  doublereal v = zz * GasConstant * TKelvin / pres;
1419  Vroot[0] = v;
1420  return 1;
1421  }
1422  }
1423 
1424 
1425  int nSolnValues;
1426  nTurningPoints = 2;
1427 
1428 #ifdef PRINTPV
1429  double V[100];
1430  int n = 0;
1431  for (int i = 0; i < 90; i++) {
1432  V[n++] = 0.030 + 0.005 * i;
1433  }
1434  double p1, presCalc;
1435  for (int i = 0; i < n; i++) {
1436  p1 = dpdVCalc(TKelvin, V[i], presCalc);
1437  printf(" %13.5g %13.5g %13.5g \n", V[i], presCalc , p1);
1438  }
1439 #endif
1440 
1441  double h2 = 4. * an * an * delta2 * delta2 * delta2;
1442  if (delta2 == 0.0) {
1443  nTurningPoints = 1;
1444  Vturn[0] = xN;
1445  Vturn[1] = xN;
1446  } else if (delta2 < 0.0) {
1447  nTurningPoints = 0;
1448  Vturn[0] = xN;
1449  Vturn[1] = xN;
1450  } else {
1451  delta = sqrt(delta2);
1452  Vturn[0] = xN - delta;
1453  Vturn[1] = xN + delta;
1454 #ifdef PRINTPV
1455  double presCalc;
1456  double p1 = dpdVCalc(TKelvin, Vturn[0], presCalc);
1457 
1458  double p2 = dpdVCalc(TKelvin, Vturn[1], presCalc);
1459 
1460  printf("p1 = %g p2 = %g \n", p1, p2);
1461  p1 = dpdVCalc(TKelvin, 0.9*Vturn[0], presCalc);
1462  printf("0.9 p1 = %g \n", p1);
1463 #endif
1464  }
1465 
1466  doublereal h = 2.0 * an * delta * delta2;
1467 
1468  doublereal yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn;
1469 
1470  doublereal desc = yN * yN - h2;
1471 
1472  if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
1473  if (desc != 0.0) {
1474  // this is for getting to other cases
1475  printf("NicholsSolve(): numerical issues\n");
1476  throw CanteraError("NicholsSolve()", "numerical issues");
1477  }
1478  desc = 0.0;
1479  }
1480 
1481  if (desc < 0.0) {
1482  nSolnValues = 3;
1483  } else if (desc == 0.0) {
1484  nSolnValues = 2;
1485  // We are here as p goes to zero.
1486  // double hleft = 3.0 * an * cn / (bn * bn);
1487  //double ynleft = 9.0 * an * cn / (2.0 * bn * bn) - 27.0 * an * an * dn / (2.0 * bn * bn * bn);
1488  //printf("hleft = %g , ynleft = %g\n", -3. / 2. * hleft, -ynleft);
1489  //double h2left = - 3 * hleft + 3 * hleft * hleft - hleft * hleft * hleft;
1490  //double y2left = - 2.0 * ynleft + ynleft * ynleft;
1491  //printf("h2left = %g , yn2left = %g\n", h2left, y2left);
1492 
1493  } else if (desc > 0.0) {
1494  nSolnValues = 1;
1495  }
1496 
1497  /*
1498  * One real root -> have to determine whether gas or liquid is the root
1499  */
1500  if (desc > 0.0) {
1501  doublereal tmpD = sqrt(desc);
1502  doublereal tmp1 = (- yN + tmpD) / (2.0 * an);
1503  doublereal sgn1 = 1.0;
1504  if (tmp1 < 0.0) {
1505  sgn1 = -1.0;
1506  tmp1 = -tmp1;
1507  }
1508  doublereal tmp2 = (- yN - tmpD) / (2.0 * an);
1509  doublereal sgn2 = 1.0;
1510  if (tmp2 < 0.0) {
1511  sgn2 = -1.0;
1512  tmp2 = -tmp2;
1513  }
1514  doublereal p1 = pow(tmp1, 1./3.);
1515  doublereal p2 = pow(tmp2, 1./3.);
1516 
1517  doublereal alpha = xN + sgn1 * p1 + sgn2 * p2;
1518  Vroot[0] = alpha;
1519  Vroot[1] = 0.0;
1520  Vroot[2] = 0.0;
1521 
1522  tmp = an * Vroot[0] * Vroot[0] * Vroot[0] + bn * Vroot[0] * Vroot[0] + cn * Vroot[0] + dn;
1523  if (fabs(tmp) > 1.0E-4) {
1524  lotsOfNumError = true;
1525  }
1526 
1527  } else if (desc < 0.0) {
1528  doublereal tmp = - yN/h;
1529 
1530  doublereal val = acos(tmp);
1531  doublereal theta = val / 3.0;
1532 
1533  doublereal oo = 2. * Cantera::Pi / 3.;
1534  doublereal alpha = xN + 2. * delta * cos(theta);
1535 
1536  doublereal beta = xN + 2. * delta * cos(theta + oo);
1537 
1538  doublereal gamma = xN + 2. * delta * cos(theta + 2.0 * oo);
1539 
1540 
1541  Vroot[0] = beta;
1542  Vroot[1] = gamma;
1543  Vroot[2] = alpha;
1544 
1545  for (int i = 0; i < 3; i++) {
1546  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1547  if (fabs(tmp) > 1.0E-4) {
1548  lotsOfNumError = true;
1549  for (int j = 0; j < 3; j++) {
1550  if (j != i) {
1551  if (fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
1552  writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " +
1553  fp2str(pres) + "): WARNING roots have merged: " +
1554  fp2str(Vroot[i]) + ", " + fp2str(Vroot[j]));
1555  writelogendl();
1556  }
1557  }
1558  }
1559  }
1560  }
1561  } else if (desc == 0.0) {
1562  if (yN == 0.0 && h == 0.0) {
1563  Vroot[0] = xN;
1564  Vroot[1] = xN;
1565  Vroot[2] = xN;
1566  } else {
1567  // need to figure out whether delta is pos or neg
1568  if (yN > 0.0) {
1569  tmp = pow(yN/(2*an), 1./3.);
1570  if (fabs(tmp - delta) > 1.0E-9) {
1571  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1572  }
1573  Vroot[1] = xN + delta;
1574  Vroot[0] = xN - 2.0*delta; // liquid phase root
1575  } else {
1576  tmp = pow(yN/(2*an), 1./3.);
1577  if (fabs(tmp - delta) > 1.0E-9) {
1578  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1579  }
1580  delta = -delta;
1581  Vroot[0] = xN + delta;
1582  Vroot[1] = xN - 2.0*delta; // gas phase root
1583  }
1584  }
1585  for (int i = 0; i < 2; i++) {
1586  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1587  if (fabs(tmp) > 1.0E-4) {
1588  lotsOfNumError = true;
1589  }
1590  }
1591  }
1592 
1593  /*
1594  * Unfortunately, there is a heavy amount of roundoff error due to bad conditioning in this
1595  */
1596  double res, dresdV = 0.0;
1597  for (int i = 0; i < nSolnValues; i++) {
1598  for (int n = 0; n < 20; n++) {
1599  res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1600  if (fabs(res) < 1.0E-14) {
1601  break;
1602  }
1603  dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
1604  double del = - res / dresdV;
1605 
1606  Vroot[i] += del;
1607  if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
1608  break;
1609  }
1610  double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1611  if (fabs(res2) < fabs(res)) {
1612  continue;
1613  } else {
1614  Vroot[i] -= del;
1615  Vroot[i] += 0.1 * del;
1616  }
1617  }
1618  if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
1619  writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " +
1620  fp2str(pres) + "): WARNING root didn't converge V = " + fp2str(Vroot[i]));
1621  writelogendl();
1622  }
1623  }
1624 
1625  if (nSolnValues == 1) {
1626  if (TKelvin > tc) {
1627  if (Vroot[0] < vc) {
1628  nSolnValues = -1;
1629  }
1630  } else {
1631  if (Vroot[0] < xN) {
1632  nSolnValues = -1;
1633  }
1634  }
1635 
1636  } else {
1637  if (nSolnValues == 2) {
1638  if (delta > 0.0) {
1639  nSolnValues = -2;
1640  }
1641  }
1642  }
1643  // writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + fp2str(pres) + "): finished");
1644  // writelogendl();
1645  return nSolnValues;
1646 }
1647 
1648 }
int NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b, doublereal Vroot[3]) const
Solve the cubic equation of state.
RedlichKwongMFTP & operator=(const RedlichKwongMFTP &right)
Assignment operator.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:534
This class can handle either an ideal solution or an ideal gas approximation of a phase...
doublereal m_b_current
Value of b in the equation of state.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize a ThermoPhase object, potentially reading activity coefficient information from an XML dat...
doublereal dpdV_
The derivative of the pressure wrt the volume.
XML_Node * get_XML_File(const std::string &file, int debug)
Return a pointer to the XML tree for a Cantera input file.
Definition: global.cpp:144
virtual doublereal isothermalCompressibility() const
Returns the isothermal compressibility. Units: 1/Pa.
virtual void getPartialMolarIntEnergies(doublereal *ubar) const
Get the species partial molar enthalpies. Units: J/kmol.
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:534
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 pressure() const
Return the thermodynamic pressure (Pa).
vector_fp m_tmpV
Temporary storage - length = m_kk.
virtual void getGibbs_ref(doublereal *g) const
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
doublereal _RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:977
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
RedlichKwongMFTP()
Base constructor.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
doublereal z() const
Calculate the value of z.
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phase, doublereal rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values, and then normalize them so that they sum to 1...
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
This is a filter class for ThermoPhase that implements some preparatory steps for efficiently handlin...
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplicator from the ThermoPhase parent class.
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
MixtureFugacityTP & operator=(const MixtureFugacityTP &b)
Assignment operator.
Header file for implicit nonlinear solver of a one dimensional function (see Numerical Utilities with...
void readXMLPureFluid(XML_Node &pureFluidParam)
Read the pure species RedlichKwong input parameters.
virtual void getPartialMolarCp(doublereal *cpbar) const
Get the partial molar heat capacities Units: J/kmol/K.
virtual int eosType() const
Equation of state type flag.
void setFuncIsGenerallyDecreasing(bool value)
Set the function behavior flag.
Definition: RootFind.cpp:1169
virtual void getGibbs_RT_ref(doublereal *grt) const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
const doublereal Pi
Pi.
Definition: ct_defs.h:51
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:597
std::string lowercase(const std::string &s)
Cast a copy of a string to lower case.
Definition: stringUtils.cpp:58
virtual void setMassFractions(const doublereal *const y)
Set the mass fractions to the specified values, and then normalize them so that they sum to 1...
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
virtual void getPartialMolarVolumes(doublereal *vbar) const
Get the species partial molar volumes. Units: m^3/kmol.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize a ThermoPhase object, potentially reading activity coefficient information from an XML dat...
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
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
void setToEquilState(const doublereal *lambda_RT)
This method is used by the ChemEquil equilibrium solver.
vector_fp m_s0_R
Temporary storage for dimensionless reference state entropies.
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...
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:623
doublereal m_Pcurrent
Current value of the pressures.
int m_standardMixingRules
boolean indicating whether standard mixing rules are applied
size_t speciesIndex(const std::string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:229
static const doublereal omega_b
Omega constant for b.
bool importPhase(XML_Node &phase, ThermoPhase *th, SpeciesThermoFactory *spfactory)
Import a phase information into an empty thermophase object.
virtual void setMoleFractions_NoNorm(const doublereal *const x)
Set the mole fractions to the specified values without normalizing.
int m_formTempParam
Form of the temperature parameterization.
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
void setTol(doublereal rtolf, doublereal atolf, doublereal rtolx=0.0, doublereal atolx=0.0)
Set the tolerance parameters for the rootfinder.
Definition: RootFind.cpp:1140
virtual void setMassFractions(const doublereal *const y)
Set the mass fractions to the specified values, and then normalize them so that they sum to 1...
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
Root finder for 1D problems.
Definition: RootFind.h:131
void setFuncIsGenerallyIncreasing(bool value)
Set the function behavior flag.
Definition: RootFind.cpp:1161
doublereal dpdT_
The derivative of the pressure wrt the temperature.
doublereal m_a_current
Value of a in the equation of state.
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture...
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:607
virtual void setConcentrations(const doublereal *const c)
Set the concentrations to the specified values within the phase.
doublereal sum_xlogx() const
Evaluate .
Definition: Phase.cpp:633
void pressureDerivatives() const
Calculate dpdV and dpdT at the current conditions.
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:29
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values, and then normalize them so that they sum to 1...
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
virtual doublereal logStandardConc(size_t k=0) const
Returns the natural logarithm of the standard concentration of the kth species.
vector_fp m_pp
Temporary storage - length = m_kk.
std::string name() const
Returns the name of the XML node.
Definition: xml.h:390
virtual void setParametersFromXML(const XML_Node &thermoNode)
Set equation of state parameter values from XML entries.
Definition file for a derived class of ThermoPhase that assumes either an ideal gas or ideal solution...
virtual doublereal dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal &presCalc) const
Calculate the pressure and the pressure derivative given the temperature and the molar volume...
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
virtual doublereal standardConcentration(size_t k=0) const
Returns the standard concentration , which is used to normalize the generalized concentration.
std::vector< doublereal > moleFractions_
Storage for the current values of the mole fractions of the species.
virtual doublereal refPressure(size_t k=npos) const =0
The reference-state pressure for species k.
virtual void setConcentrations(const doublereal *const c)
Set the concentrations to the specified values within the phase.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
virtual void setTemperature(const doublereal temp)
Set the temperature (K)
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture...
virtual doublereal critDensity() const
Critical density (kg/m3).
void readXMLCrossFluid(XML_Node &pureFluidParam)
Read the cross species RedlichKwong input parameters.
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:574
doublereal dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition: utilities.h:131
static const doublereal omega_a
Omega constant for a -> value of a in terms of critical properties.
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.
virtual void getIntEnergy_RT(doublereal *urt) const
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
virtual doublereal intEnergy_mole() const
Molar internal energy. Units: J/kmol.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:159
int solve(doublereal xmin, doublereal xmax, int itmax, doublereal &funcTargetValue, doublereal *xbest)
Using a line search method, find the root of a 1D function.
Definition: RootFind.cpp:184
virtual doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
void writelogendl()
Write an end of line character to the screen and flush output.
Definition: global.cpp:66
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:524
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
size_t nDim() const
Returns the number of spatial dimensions (1, 2, or 3)
Definition: Phase.h:512
virtual doublereal critPressure() const
Critical pressure (Pa).
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:139
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
void calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const
Calculate the a and the b parameters given the temperature.
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:562
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:154
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Definition: ThermoPhase.h:1489
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:588
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
virtual void setState_PX(doublereal p, doublereal *x)
Set the pressure (Pa) and mole fractions.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Get the species partial molar enthalpies. Units: J/kmol.
static const doublereal omega_vc
Omega constant for the critical molar volume.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
void writelog(const std::string &msg)
Write a message to the screen.
Definition: global.cpp:43
virtual void getUnitsStandardConc(double *uA, int k=0, int sizeUA=6) const
Returns the units of the standard and generalized concentrations.
size_t m_kk
Number of species in the phase.
Definition: Phase.h:716
void applyStandardMixingRules()
Apply mixing rules for a coefficients.
virtual void getEntropy_R_ref(doublereal *er) const
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const
Calculate the pressure given the temperature and the molar volume.
void setPrintLvl(int printLvl)
Set the print level from the rootfinder.
Definition: RootFind.cpp:1156
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials.
virtual doublereal critTemperature() const
Critical temperature (K).
virtual void getPartialMolarEntropies(doublereal *sbar) const
Get the species partial molar entropies. Units: J/kmol/K.
#define ROOTFIND_SUCCESS
This means that the root solver was a success.
Definition: RootFind.h:29
vector_fp dpdni_
Vector of derivatives of pressure wrt mole number.
SpeciesThermo * m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1625
void updateAB()
Update the a and b parameters.
const doublereal * moleFractdivMMW() const
Returns a const pointer to the start of the moleFraction/MW array.
Definition: Phase.cpp:540
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
XML_Node * get_XML_NameID(const std::string &nameTarget, const std::string &file_ID, XML_Node *root)
This routine will locate an XML node in either the input XML tree or in another input file specified ...
Definition: global.cpp:271
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:594
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional activity coefficients at the current solution temperature, pressure, and solution concentration.
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase Note the density of a phase is an independent...
Definition: Phase.h:549
vector_fp m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.