Cantera  2.3.0
RedlichKwongMFTP.cpp
Go to the documentation of this file.
1 //! @file RedlichKwongMFTP.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at http://www.cantera.org/license.txt for license and copyright information.
5 
10 #include "cantera/base/ctml.h"
11 
12 #include <boost/math/tools/roots.hpp>
13 
14 using namespace std;
15 namespace bmt = boost::math::tools;
16 
17 namespace Cantera
18 {
19 
20 const doublereal RedlichKwongMFTP::omega_a = 4.27480233540E-01;
21 const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
22 const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
23 
24 RedlichKwongMFTP::RedlichKwongMFTP() :
25  m_standardMixingRules(0),
26  m_formTempParam(0),
27  m_b_current(0.0),
28  m_a_current(0.0),
29  NSolns_(0),
30  Vroot_{0.0, 0.0, 0.0},
31  dpdV_(0.0),
32  dpdT_(0.0)
33 {
34 }
35 
36 RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string& id_) :
37  m_standardMixingRules(0),
38  m_formTempParam(0),
39  m_b_current(0.0),
40  m_a_current(0.0),
41  NSolns_(0),
42  Vroot_{0.0, 0.0, 0.0},
43  dpdV_(0.0),
44  dpdT_(0.0)
45 {
46  initThermoFile(infile, id_);
47 }
48 
49 RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id_) :
50  m_standardMixingRules(0),
51  m_formTempParam(0),
52  m_b_current(0.0),
53  m_a_current(0.0),
54  NSolns_(0),
55  Vroot_{0.0, 0.0, 0.0},
56  dpdV_(0.0),
57  dpdT_(0.0)
58 {
59  importPhase(phaseRefRoot, this);
60 }
61 
62 RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP& b) :
63  m_standardMixingRules(0),
64  m_formTempParam(0),
65  m_b_current(0.0),
66  m_a_current(0.0),
67  NSolns_(0),
68  dpdV_(0.0),
69  dpdT_(0.0)
70 {
71  *this = b;
72 }
73 
74 RedlichKwongMFTP& RedlichKwongMFTP::operator=(const RedlichKwongMFTP& b)
75 {
76  if (&b != this) {
77  // Mostly, this is a passthrough to the underlying assignment operator
78  // for the ThermoPhae parent object.
79  MixtureFugacityTP::operator=(b);
80 
81  // However, we have to handle data that we own.
82  m_standardMixingRules = b.m_standardMixingRules;
83  m_formTempParam = b.m_formTempParam;
84  m_b_current = b.m_b_current;
85  m_a_current = b.m_a_current;
86  a_vec_Curr_ = b.a_vec_Curr_;
87  b_vec_Curr_ = b.b_vec_Curr_;
88  a_coeff_vec = b.a_coeff_vec;
89 
90  m_pc_Species = b.m_pc_Species;
91  m_tc_Species = b.m_tc_Species;
92  m_vc_Species = b.m_vc_Species;
93  NSolns_ = b.NSolns_;
94  Vroot_[0] = b.Vroot_[0];
95  Vroot_[1] = b.Vroot_[1];
96  Vroot_[2] = b.Vroot_[2];
97  m_pp = b.m_pp;
98  m_tmpV = b.m_tmpV;
99  m_partialMolarVolumes = b.m_partialMolarVolumes;
100  dpdV_ = b.dpdV_;
101  dpdT_ = b.dpdT_;
102  dpdni_ = b.dpdni_;
103  }
104  return *this;
105 }
106 
108 {
109  return new RedlichKwongMFTP(*this);
110 }
111 
113 {
114  warn_deprecated("RedlichKwongMFTP::eosType",
115  "To be removed after Cantera 2.3.");
116  return cRedlichKwongMFTP;
117 }
118 
119 // ------------Molar Thermodynamic Properties -------------------------
120 
122 {
124  doublereal h_ideal = RT() * mean_X(m_h0_RT);
125  doublereal h_nonideal = hresid();
126  return h_ideal + h_nonideal;
127 }
128 
130 {
132  doublereal sr_ideal = GasConstant * (mean_X(m_s0_R)
133  - sum_xlogx() - std::log(pressure()/m_spthermo->refPressure()));
134  doublereal sr_nonideal = sresid();
135  return sr_ideal + sr_nonideal;
136 }
137 
138 doublereal RedlichKwongMFTP::cp_mole() const
139 {
141  doublereal TKelvin = temperature();
142  doublereal sqt = sqrt(TKelvin);
143  doublereal mv = molarVolume();
144  doublereal vpb = mv + m_b_current;
146  doublereal cpref = GasConstant * mean_X(m_cp0_R);
147  doublereal dadt = da_dt();
148  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
149  doublereal dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac
150  +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
151  return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_;
152 }
153 
154 doublereal RedlichKwongMFTP::cv_mole() const
155 {
157  doublereal TKelvin = temperature();
158  doublereal sqt = sqrt(TKelvin);
159  doublereal mv = molarVolume();
160  doublereal vpb = mv + m_b_current;
161  doublereal cvref = GasConstant * (mean_X(m_cp0_R) - 1.0);
162  doublereal dadt = da_dt();
163  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
164  return (cvref - 1.0/(2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv)*fac
165  +1.0/(m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt));
166 }
167 
168 doublereal RedlichKwongMFTP::pressure() const
169 {
171 
172  // Get a copy of the private variables stored in the State object
173  doublereal T = temperature();
174  double molarV = meanMolecularWeight() / density();
175  double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current));
176  if (fabs(pp -m_Pcurrent) > 1.0E-5 * fabs(m_Pcurrent)) {
177  throw CanteraError(" RedlichKwongMFTP::pressure()", "setState broken down, maybe");
178  }
179 
180  return m_Pcurrent;
181 }
182 
184 {
185  // Calculate the molarVolume of the solution (m**3 kmol-1)
186  const doublereal* const dtmp = moleFractdivMMW();
188  double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp);
189 
190  // Set the density in the parent State object directly, by calling the
191  // Phase::setDensity() function.
192  Phase::setDensity(1.0/invDens);
193 }
194 
195 void RedlichKwongMFTP::setTemperature(const doublereal temp)
196 {
197  Phase::setTemperature(temp);
199  updateAB();
200 }
201 
203 {
205  updateAB();
206 }
207 
209 {
211  for (size_t k = 0; k < m_kk; k++) {
212  c[k] *= moleFraction(k)*pressure()/RT();
213  }
214 }
215 
216 doublereal RedlichKwongMFTP::standardConcentration(size_t k) const
217 {
218  getStandardVolumes(m_tmpV.data());
219  return 1.0 / m_tmpV[k];
220 }
221 
223 {
224  doublereal mv = molarVolume();
225  doublereal sqt = sqrt(temperature());
226  doublereal vpb = mv + m_b_current;
227  doublereal vmb = mv - m_b_current;
228 
229  for (size_t k = 0; k < m_kk; k++) {
230  m_pp[k] = 0.0;
231  for (size_t i = 0; i < m_kk; i++) {
232  size_t counter = k + m_kk*i;
233  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
234  }
235  }
236  doublereal pres = pressure();
237 
238  for (size_t k = 0; k < m_kk; k++) {
239  ac[k] = (- RT() * log(pres * mv / RT())
240  + RT() * log(mv / vmb)
241  + RT() * b_vec_Curr_[k] / vmb
242  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
243  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
244  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
245  );
246  }
247  for (size_t k = 0; k < m_kk; k++) {
248  ac[k] = exp(ac[k]/RT());
249  }
250 }
251 
252 // ---- Partial Molar Properties of the Solution -----------------
253 
254 void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
255 {
256  getChemPotentials(muRT);
257  for (size_t k = 0; k < m_kk; k++) {
258  muRT[k] *= 1.0 / RT();
259  }
260 }
261 
262 void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
263 {
264  getGibbs_ref(mu);
265  for (size_t k = 0; k < m_kk; k++) {
266  double xx = std::max(SmallNumber, moleFraction(k));
267  mu[k] += RT()*(log(xx));
268  }
269 
270  doublereal mv = molarVolume();
271  doublereal sqt = sqrt(temperature());
272  doublereal vpb = mv + m_b_current;
273  doublereal vmb = mv - m_b_current;
274 
275  for (size_t k = 0; k < m_kk; k++) {
276  m_pp[k] = 0.0;
277  for (size_t i = 0; i < m_kk; i++) {
278  size_t counter = k + m_kk*i;
279  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
280  }
281  }
282  doublereal pres = pressure();
283  doublereal refP = refPressure();
284 
285  for (size_t k = 0; k < m_kk; k++) {
286  mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
287  + RT() * log(mv / vmb)
288  + RT() * b_vec_Curr_[k] / vmb
289  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
290  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
291  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
292  );
293  }
294 }
295 
297 {
298  // First we get the reference state contributions
299  getEnthalpy_RT_ref(hbar);
300  scale(hbar, hbar+m_kk, hbar, RT());
301 
302  // We calculate dpdni_
303  doublereal TKelvin = temperature();
304  doublereal mv = molarVolume();
305  doublereal sqt = sqrt(TKelvin);
306  doublereal vpb = mv + m_b_current;
307  doublereal vmb = mv - m_b_current;
308  for (size_t k = 0; k < m_kk; k++) {
309  m_pp[k] = 0.0;
310  for (size_t i = 0; i < m_kk; i++) {
311  size_t counter = k + m_kk*i;
312  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
313  }
314  }
315  for (size_t k = 0; k < m_kk; k++) {
316  dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
317  + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
318  }
319  doublereal dadt = da_dt();
320  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
321 
322  for (size_t k = 0; k < m_kk; k++) {
323  m_tmpV[k] = 0.0;
324  for (size_t i = 0; i < m_kk; i++) {
325  size_t counter = k + m_kk*i;
326  m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
327  }
328  }
329 
331  doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_;
332  for (size_t k = 0; k < m_kk; k++) {
333  double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
334  + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
335  + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
336  hbar[k] = hbar[k] + hE_v;
337  hbar[k] -= fac2 * dpdni_[k];
338  }
339 }
340 
341 void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
342 {
343  getEntropy_R_ref(sbar);
344  scale(sbar, sbar+m_kk, sbar, GasConstant);
345  doublereal TKelvin = temperature();
346  doublereal sqt = sqrt(TKelvin);
347  doublereal mv = molarVolume();
348  doublereal refP = refPressure();
349 
350  for (size_t k = 0; k < m_kk; k++) {
351  doublereal xx = std::max(SmallNumber, moleFraction(k));
352  sbar[k] += GasConstant * (- log(xx));
353  }
354  for (size_t k = 0; k < m_kk; k++) {
355  m_pp[k] = 0.0;
356  for (size_t i = 0; i < m_kk; i++) {
357  size_t counter = k + m_kk*i;
358  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
359  }
360  }
361  for (size_t k = 0; k < m_kk; k++) {
362  m_tmpV[k] = 0.0;
363  for (size_t i = 0; i < m_kk; i++) {
364  size_t counter = k + m_kk*i;
365  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
366  }
367  }
368 
369  doublereal dadt = da_dt();
370  doublereal fac = dadt - m_a_current / (2.0 * TKelvin);
371  doublereal vmb = mv - m_b_current;
372  doublereal vpb = mv + m_b_current;
373  for (size_t k = 0; k < m_kk; k++) {
374  sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
375  + GasConstant
376  + GasConstant * log(mv/vmb)
377  + GasConstant * b_vec_Curr_[k]/vmb
378  + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
379  - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
380  + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
381  - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
382  );
383  }
384 
386  getPartialMolarVolumes(m_partialMolarVolumes.data());
387  for (size_t k = 0; k < m_kk; k++) {
388  sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
389  }
390 }
391 
393 {
394  getIntEnergy_RT(ubar);
395  scale(ubar, ubar+m_kk, ubar, RT());
396 }
397 
398 void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const
399 {
400  getCp_R(cpbar);
401  scale(cpbar, cpbar+m_kk, cpbar, GasConstant);
402 }
403 
404 void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
405 {
406  for (size_t k = 0; k < m_kk; k++) {
407  m_pp[k] = 0.0;
408  for (size_t i = 0; i < m_kk; i++) {
409  size_t counter = k + m_kk*i;
410  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
411  }
412  }
413  for (size_t k = 0; k < m_kk; k++) {
414  m_tmpV[k] = 0.0;
415  for (size_t i = 0; i < m_kk; i++) {
416  size_t counter = k + m_kk*i;
417  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
418  }
419  }
420 
421  doublereal sqt = sqrt(temperature());
422  doublereal mv = molarVolume();
423  doublereal vmb = mv - m_b_current;
424  doublereal vpb = mv + m_b_current;
425  for (size_t k = 0; k < m_kk; k++) {
426  doublereal num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
427  + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
428  - 2.0 * m_pp[k] / (sqt * vpb)
429  + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
430  );
431  doublereal denom = (m_Pcurrent + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
432  );
433  vbar[k] = num / denom;
434  }
435 }
436 
438 {
439  double pc, tc, vc;
440  double a0 = 0.0;
441  double aT = 0.0;
442  for (size_t i = 0; i < m_kk; i++) {
443  for (size_t j = 0; j <m_kk; j++) {
444  size_t counter = i + m_kk * j;
445  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
446  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
447  }
448  }
449  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
450  return tc;
451 }
452 
454 {
455  double pc, tc, vc;
456  double a0 = 0.0;
457  double aT = 0.0;
458  for (size_t i = 0; i < m_kk; i++) {
459  for (size_t j = 0; j <m_kk; j++) {
460  size_t counter = i + m_kk * j;
461  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
462  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
463  }
464  }
465  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
466  return pc;
467 }
468 
470 {
471  double pc, tc, vc;
472  double a0 = 0.0;
473  double aT = 0.0;
474  for (size_t i = 0; i < m_kk; i++) {
475  for (size_t j = 0; j <m_kk; j++) {
476  size_t counter = i + m_kk * j;
477  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
478  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
479  }
480  }
481  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
482  return vc;
483 }
484 
486 {
487  double pc, tc, vc;
488  double a0 = 0.0;
489  double aT = 0.0;
490  for (size_t i = 0; i < m_kk; i++) {
491  for (size_t j = 0; j <m_kk; j++) {
492  size_t counter = i + m_kk * j;
493  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
494  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
495  }
496  }
497  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
498  return pc*vc/tc/GasConstant;
499 }
500 
502 {
503  double pc, tc, vc;
504  double a0 = 0.0;
505  double aT = 0.0;
506  for (size_t i = 0; i < m_kk; i++) {
507  for (size_t j = 0; j <m_kk; j++) {
508  size_t counter = i + m_kk * j;
509  a0 += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(0, counter);
510  aT += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(1, counter);
511  }
512  }
513  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
514  double mmw = meanMolecularWeight();
515  return mmw / vc;
516 }
517 
518 void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
519 {
520  double tmp, tmp2;
522  getGibbs_RT_ref(m_tmpV.data());
523 
524  // Within the method, we protect against inf results if the exponent is too
525  // high.
526  //
527  // If it is too low, we set the partial pressure to zero. This capability is
528  // needed by the elemental potential method.
529  doublereal pres = 0.0;
530  double m_p0 = refPressure();
531  for (size_t k = 0; k < m_kk; k++) {
532  tmp = -m_tmpV[k] + mu_RT[k];
533  if (tmp < -600.) {
534  m_pp[k] = 0.0;
535  } else if (tmp > 500.0) {
536  tmp2 = tmp / 500.;
537  tmp2 *= tmp2;
538  m_pp[k] = m_p0 * exp(500.) * tmp2;
539  } else {
540  m_pp[k] = m_p0 * exp(tmp);
541  }
542  pres += m_pp[k];
543  }
544  // set state
545  setState_PX(pres, &m_pp[0]);
546 }
547 
548 bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
549 {
550  bool added = MixtureFugacityTP::addSpecies(spec);
551  if (added) {
552  a_vec_Curr_.resize(m_kk * m_kk, 0.0);
553  b_vec_Curr_.push_back(0.0);
554 
555  a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
556 
557  m_pc_Species.push_back(0.0);
558  m_tc_Species.push_back(0.0);
559  m_vc_Species.push_back(0.0);
560 
561  m_pp.push_back(0.0);
562  m_tmpV.push_back(0.0);
563  m_partialMolarVolumes.push_back(0.0);
564  dpdni_.push_back(0.0);
565  }
566  return added;
567 }
568 
569 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
570 {
571  // Check the model parameter for the Redlich-Kwong equation of state
572  // two are allowed
573  // RedlichKwong mixture of species, each of which are RK fluids
574  // RedlichKwongMFTP mixture of species with cross term coefficients
575  if (phaseNode.hasChild("thermo")) {
576  XML_Node& thermoNode = phaseNode.child("thermo");
577  std::string model = thermoNode["model"];
578  if (model == "RedlichKwong") {
580  } else if (model == "RedlichKwongMFTP") {
582  } else {
583  throw CanteraError("RedlichKwongMFTP::initThermoXML",
584  "Unknown thermo model : " + model);
585  }
586 
587  // Go get all of the coefficients and factors in the
588  // activityCoefficients XML block
589  XML_Node* acNodePtr = 0;
590  if (thermoNode.hasChild("activityCoefficients")) {
591  XML_Node& acNode = thermoNode.child("activityCoefficients");
592  acNodePtr = &acNode;
593  size_t nC = acNode.nChildren();
594 
595  // Loop through the children getting multiple instances of
596  // parameters
597  for (size_t i = 0; i < nC; i++) {
598  XML_Node& xmlACChild = acNodePtr->child(i);
599  if (ba::iequals(xmlACChild.name(), "purefluidparameters")) {
600  readXMLPureFluid(xmlACChild);
601  }
602  }
603  if (m_standardMixingRules == 1) {
605  }
606 
607  // Loop through the children getting multiple instances of
608  // parameters
609  for (size_t i = 0; i < nC; i++) {
610  XML_Node& xmlACChild = acNodePtr->child(i);
611  if (ba::iequals(xmlACChild.name(), "crossfluidparameters")) {
612  readXMLCrossFluid(xmlACChild);
613  }
614  }
615  }
616  }
617 
618  for (size_t i = 0; i < m_kk; i++) {
619  double a0coeff = a_coeff_vec(0, i*m_kk + i);
620  double aTcoeff = a_coeff_vec(1, i*m_kk + i);
621  double ai = a0coeff + aTcoeff * 500.;
622  double bi = b_vec_Curr_[i];
623  calcCriticalConditions(ai, bi, a0coeff, aTcoeff, m_pc_Species[i], m_tc_Species[i], m_vc_Species[i]);
624  }
625 
626  MixtureFugacityTP::initThermoXML(phaseNode, id);
627 }
628 
630 {
631  vector_fp vParams;
632  string xname = pureFluidParam.name();
633  if (xname != "pureFluidParameters") {
634  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
635  "Incorrect name for processing this routine: " + xname);
636  }
637 
638  // Read the species. Find the index of the species in the current phase.
639  // It's not an error to not find the species
640  string iName = pureFluidParam.attrib("species");
641  if (iName == "") {
642  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute");
643  }
644  size_t iSpecies = speciesIndex(iName);
645  if (iSpecies == npos) {
646  return;
647  }
648  size_t counter = iSpecies + m_kk * iSpecies;
649  size_t nParamsExpected, nParamsFound;
650  size_t num = pureFluidParam.nChildren();
651  for (size_t iChild = 0; iChild < num; iChild++) {
652  XML_Node& xmlChild = pureFluidParam.child(iChild);
653  string nodeName = ba::to_lower_copy(xmlChild.name());
654 
655  if (nodeName == "a_coeff") {
656  string iModel = ba::to_lower_copy(xmlChild.attrib("model"));
657  if (iModel == "constant") {
658  nParamsExpected = 1;
659  } else if (iModel == "linear_a") {
660  nParamsExpected = 2;
661  if (m_formTempParam == 0) {
662  m_formTempParam = 1;
663  }
664  } else {
665  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "unknown model");
666  }
667 
668  getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
669  nParamsFound = vParams.size();
670  if (nParamsFound != nParamsExpected) {
671  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for a_coeff" + iName + ")",
672  "wrong number of params found");
673  }
674 
675  for (size_t i = 0; i < nParamsFound; i++) {
676  a_coeff_vec(i, counter) = vParams[i];
677  }
678  } else if (nodeName == "b_coeff") {
679  getFloatArray(xmlChild, vParams, true, "m3/kmol", "b_coeff");
680  nParamsFound = vParams.size();
681  if (nParamsFound != 1) {
682  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for b_coeff" + iName + ")",
683  "wrong number of params found");
684  }
685  b_vec_Curr_[iSpecies] = vParams[0];
686  }
687  }
688 }
689 
691 {
692  int nParam = 2;
693  for (size_t i = 0; i < m_kk; i++) {
694  size_t icounter = i + m_kk * i;
695  for (size_t j = 0; j < m_kk; j++) {
696  if (i != j) {
697  size_t counter = i + m_kk * j;
698  size_t jcounter = j + m_kk * j;
699  for (int n = 0; n < nParam; n++) {
700  a_coeff_vec(n, counter) = sqrt(a_coeff_vec(n, icounter) * a_coeff_vec(n, jcounter));
701  }
702  }
703  }
704  }
705 }
706 
708 {
709  vector_fp vParams;
710  string xname = CrossFluidParam.name();
711  if (xname != "crossFluidParameters") {
712  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
713  "Incorrect name for processing this routine: " + xname);
714  }
715 
716  // Read the species. Find the index of the species in the current phase.
717  // It's not an error to not find the species
718  string iName = CrossFluidParam.attrib("species1");
719  if (iName == "") {
720  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute");
721  }
722  size_t iSpecies = speciesIndex(iName);
723  if (iSpecies == npos) {
724  return;
725  }
726  string jName = CrossFluidParam.attrib("species2");
727  if (iName == "") {
728  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species2 attribute");
729  }
730  size_t jSpecies = speciesIndex(jName);
731  if (jSpecies == npos) {
732  return;
733  }
734 
735  size_t counter = iSpecies + m_kk * jSpecies;
736  size_t counter0 = jSpecies + m_kk * iSpecies;
737  size_t nParamsExpected, nParamsFound;
738  size_t num = CrossFluidParam.nChildren();
739  for (size_t iChild = 0; iChild < num; iChild++) {
740  XML_Node& xmlChild = CrossFluidParam.child(iChild);
741  string nodeName = ba::to_lower_copy(xmlChild.name());
742 
743  if (nodeName == "a_coeff") {
744  string iModel = ba::to_lower_copy(xmlChild.attrib("model"));
745  if (iModel == "constant") {
746  nParamsExpected = 1;
747  } else if (iModel == "linear_a") {
748  nParamsExpected = 2;
749  if (m_formTempParam == 0) {
750  m_formTempParam = 1;
751  }
752  } else {
753  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "unknown model");
754  }
755 
756  getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
757  nParamsFound = vParams.size();
758  if (nParamsFound != nParamsExpected) {
759  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid(for a_coeff" + iName + ")",
760  "wrong number of params found");
761  }
762 
763  for (size_t i = 0; i < nParamsFound; i++) {
764  a_coeff_vec(i, counter) = vParams[i];
765  a_coeff_vec(i, counter0) = vParams[i];
766  }
767  }
768  }
769 }
770 
772 {
774  std::string model = thermoNode["model"];
775 }
776 
777 doublereal RedlichKwongMFTP::sresid() const
778 {
779  // note this agrees with tpx
780  doublereal rho = density();
781  doublereal mmw = meanMolecularWeight();
782  doublereal molarV = mmw / rho;
783  double hh = m_b_current / molarV;
784  doublereal zz = z();
785  doublereal dadt = da_dt();
786  doublereal T = temperature();
787  doublereal sqT = sqrt(T);
788  doublereal fac = dadt - m_a_current / (2.0 * T);
789  double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
790  return GasConstant * sresid_mol_R;
791 }
792 
793 doublereal RedlichKwongMFTP::hresid() const
794 {
795  // note this agrees with tpx
796  doublereal rho = density();
797  doublereal mmw = meanMolecularWeight();
798  doublereal molarV = mmw / rho;
799  double hh = m_b_current / molarV;
800  doublereal zz = z();
801  doublereal dadt = da_dt();
802  doublereal T = temperature();
803  doublereal sqT = sqrt(T);
804  doublereal fac = T * dadt - 3.0 *m_a_current / (2.0);
805  return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
806 }
807 
808 doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
809 {
810  double v = m_b_current * 1.1;
811  double atmp;
812  double btmp;
813  calculateAB(TKelvin, atmp, btmp);
814  doublereal pres = std::max(psatEst(TKelvin), presGuess);
815  double Vroot[3];
816  bool foundLiq = false;
817  int m = 0;
818  while (m < 100 && !foundLiq) {
819  int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot);
820  if (nsol == 1 || nsol == 2) {
821  double pc = critPressure();
822  if (pres > pc) {
823  foundLiq = true;
824  }
825  pres *= 1.04;
826  } else {
827  foundLiq = true;
828  }
829  }
830 
831  if (foundLiq) {
832  v = Vroot[0];
833  presGuess = pres;
834  } else {
835  v = -1.0;
836  }
837  return v;
838 }
839 
840 doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
841 {
842  // It's necessary to set the temperature so that m_a_current is set correctly.
843  setTemperature(TKelvin);
844  double tcrit = critTemperature();
845  doublereal mmw = meanMolecularWeight();
846  if (rhoguess == -1.0) {
847  if (phaseRequested != FLUID_GAS) {
848  if (TKelvin > tcrit) {
849  rhoguess = presPa * mmw / (GasConstant * TKelvin);
850  } else {
851  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
852  rhoguess = presPa * mmw / (GasConstant * TKelvin);
853  } else if (phaseRequested >= FLUID_LIQUID_0) {
854  double lqvol = liquidVolEst(TKelvin, presPa);
855  rhoguess = mmw / lqvol;
856  }
857  }
858  } else {
859  // Assume the Gas phase initial guess, if nothing is specified to
860  // the routine
861  rhoguess = presPa * mmw / (GasConstant * TKelvin);
862  }
863  }
864 
865  doublereal volguess = mmw / rhoguess;
866  NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
867 
868  doublereal molarVolLast = Vroot_[0];
869  if (NSolns_ >= 2) {
870  if (phaseRequested >= FLUID_LIQUID_0) {
871  molarVolLast = Vroot_[0];
872  } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
873  molarVolLast = Vroot_[2];
874  } else {
875  if (volguess > Vroot_[1]) {
876  molarVolLast = Vroot_[2];
877  } else {
878  molarVolLast = Vroot_[0];
879  }
880  }
881  } else if (NSolns_ == 1) {
882  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
883  molarVolLast = Vroot_[0];
884  } else {
885  return -2.0;
886  }
887  } else if (NSolns_ == -1) {
888  if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
889  molarVolLast = Vroot_[0];
890  } else if (TKelvin > tcrit) {
891  molarVolLast = Vroot_[0];
892  } else {
893  return -2.0;
894  }
895  } else {
896  molarVolLast = Vroot_[0];
897  return -1.0;
898  }
899  return mmw / molarVolLast;
900 }
901 
903 {
904  double Vroot[3];
905  double T = temperature();
906  int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
907  if (nsol != 3) {
908  return critDensity();
909  }
910 
911  auto resid = [this, T](double v) {
912  double pp;
913  return dpdVCalc(T, v, pp);
914  };
915 
916  boost::uintmax_t maxiter = 100;
917  std::pair<double, double> vv = bmt::toms748_solve(
918  resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
919 
920  doublereal mmw = meanMolecularWeight();
921  return mmw / (0.5 * (vv.first + vv.second));
922 }
923 
925 {
926  double Vroot[3];
927  double T = temperature();
928  int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
929  if (nsol != 3) {
930  return critDensity();
931  }
932 
933  auto resid = [this, T](double v) {
934  double pp;
935  return dpdVCalc(T, v, pp);
936  };
937 
938  boost::uintmax_t maxiter = 100;
939  std::pair<double, double> vv = bmt::toms748_solve(
940  resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
941 
942  doublereal mmw = meanMolecularWeight();
943  return mmw / (0.5 * (vv.first + vv.second));
944 }
945 
946 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
947 {
948  doublereal sqt = sqrt(TKelvin);
949  double pres = GasConstant * TKelvin / (molarVol - m_b_current)
950  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
951  return pres;
952 }
953 
954 doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
955 {
956  doublereal sqt = sqrt(TKelvin);
957  presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
958  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
959 
960  doublereal vpb = molarVol + m_b_current;
961  doublereal vmb = molarVol - m_b_current;
962  doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb)
963  + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
964  return dpdv;
965 }
966 
968 {
969  doublereal TKelvin = temperature();
970  doublereal mv = molarVolume();
971  doublereal pres;
972 
973  dpdV_ = dpdVCalc(TKelvin, mv, pres);
974  doublereal sqt = sqrt(TKelvin);
975  doublereal vpb = mv + m_b_current;
976  doublereal vmb = mv - m_b_current;
977  doublereal dadt = da_dt();
978  doublereal fac = dadt - m_a_current/(2.0 * TKelvin);
979  dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
980 }
981 
982 void RedlichKwongMFTP::updateMixingExpressions()
983 {
984  updateAB();
985 }
986 
988 {
989  double temp = temperature();
990  if (m_formTempParam == 1) {
991  for (size_t i = 0; i < m_kk; i++) {
992  for (size_t j = 0; j < m_kk; j++) {
993  size_t counter = i * m_kk + j;
994  a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
995  }
996  }
997  }
998 
999  m_b_current = 0.0;
1000  m_a_current = 0.0;
1001  for (size_t i = 0; i < m_kk; i++) {
1002  m_b_current += moleFractions_[i] * b_vec_Curr_[i];
1003  for (size_t j = 0; j < m_kk; j++) {
1004  m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
1005  }
1006  }
1007 }
1008 
1009 void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
1010 {
1011  bCalc = 0.0;
1012  aCalc = 0.0;
1013  if (m_formTempParam == 1) {
1014  for (size_t i = 0; i < m_kk; i++) {
1015  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1016  for (size_t j = 0; j < m_kk; j++) {
1017  size_t counter = i * m_kk + j;
1018  doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
1019  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1020  }
1021  }
1022  } else {
1023  for (size_t i = 0; i < m_kk; i++) {
1024  bCalc += moleFractions_[i] * b_vec_Curr_[i];
1025  for (size_t j = 0; j < m_kk; j++) {
1026  size_t counter = i * m_kk + j;
1027  doublereal a_vec_Curr = a_coeff_vec(0,counter);
1028  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
1029  }
1030  }
1031  }
1032 }
1033 
1034 doublereal RedlichKwongMFTP::da_dt() const
1035 {
1036  doublereal dadT = 0.0;
1037  if (m_formTempParam == 1) {
1038  for (size_t i = 0; i < m_kk; i++) {
1039  for (size_t j = 0; j < m_kk; j++) {
1040  size_t counter = i * m_kk + j;
1041  dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
1042  }
1043  }
1044  }
1045  return dadT;
1046 }
1047 
1048 void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
1049  doublereal& pc, doublereal& tc, doublereal& vc) const
1050 {
1051  if (m_formTempParam != 0) {
1052  a = a0_coeff;
1053  }
1054  if (b <= 0.0) {
1055  tc = 1000000.;
1056  pc = 1.0E13;
1057  vc = omega_vc * GasConstant * tc / pc;
1058  return;
1059  }
1060  if (a <= 0.0) {
1061  tc = 0.0;
1062  pc = 0.0;
1063  vc = 2.0 * b;
1064  return;
1065  }
1066  double tmp = a * omega_b / (b * omega_a * GasConstant);
1067  double pp = 2./3.;
1068  doublereal sqrttc, f, dfdt, deltatc;
1069 
1070  if (m_formTempParam == 0) {
1071  tc = pow(tmp, pp);
1072  } else {
1073  tc = pow(tmp, pp);
1074  for (int j = 0; j < 10; j++) {
1075  sqrttc = sqrt(tc);
1076  f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT_coeff * tc - a0_coeff;
1077  dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT_coeff;
1078  deltatc = - f / dfdt;
1079  tc += deltatc;
1080  }
1081  if (deltatc > 0.1) {
1082  throw CanteraError("RedlichKwongMFTP::calcCriticalConditions", "didn't converge");
1083  }
1084  }
1085 
1086  pc = omega_b * GasConstant * tc / b;
1087  vc = omega_vc * GasConstant * tc / pc;
1088 }
1089 
1090 int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
1091  doublereal Vroot[3]) const
1092 {
1093  Vroot[0] = 0.0;
1094  Vroot[1] = 0.0;
1095  Vroot[2] = 0.0;
1096  if (TKelvin <= 0.0) {
1097  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature");
1098  }
1099 
1100  // Derive the coefficients of the cubic polynomial to solve.
1101  doublereal an = 1.0;
1102  doublereal bn = - GasConstant * TKelvin / pres;
1103  doublereal sqt = sqrt(TKelvin);
1104  doublereal cn = - (GasConstant * TKelvin * b / pres - a/(pres * sqt) + b * b);
1105  doublereal dn = - (a * b / (pres * sqt));
1106 
1107  double tmp = a * omega_b / (b * omega_a * GasConstant);
1108  double pp = 2./3.;
1109  double tc = pow(tmp, pp);
1110  double pc = omega_b * GasConstant * tc / b;
1111  double vc = omega_vc * GasConstant * tc / pc;
1112  // Derive the center of the cubic, x_N
1113  doublereal xN = - bn /(3 * an);
1114 
1115  // Derive the value of delta**2. This is a key quantity that determines the
1116  // number of turning points
1117  doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
1118  doublereal delta = 0.0;
1119 
1120  // Calculate a couple of ratios
1121  doublereal ratio1 = 3.0 * an * cn / (bn * bn);
1122  doublereal ratio2 = pres * b / (GasConstant * TKelvin);
1123  if (fabs(ratio1) < 1.0E-7) {
1124  doublereal ratio3 = a / (GasConstant * sqt) * pres / (GasConstant * TKelvin);
1125  if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
1126  doublereal zz = 1.0;
1127  for (int i = 0; i < 10; i++) {
1128  doublereal znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
1129  doublereal deltaz = znew - zz;
1130  zz = znew;
1131  if (fabs(deltaz) < 1.0E-14) {
1132  break;
1133  }
1134  }
1135  doublereal v = zz * GasConstant * TKelvin / pres;
1136  Vroot[0] = v;
1137  return 1;
1138  }
1139  }
1140 
1141  int nSolnValues;
1142  double h2 = 4. * an * an * delta2 * delta2 * delta2;
1143  if (delta2 > 0.0) {
1144  delta = sqrt(delta2);
1145  }
1146 
1147  doublereal h = 2.0 * an * delta * delta2;
1148  doublereal yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn;
1149  doublereal desc = yN * yN - h2;
1150 
1151  if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
1152  if (desc != 0.0) {
1153  // this is for getting to other cases
1154  throw CanteraError("NicholsSolve()", "numerical issues");
1155  }
1156  desc = 0.0;
1157  }
1158 
1159  if (desc < 0.0) {
1160  nSolnValues = 3;
1161  } else if (desc == 0.0) {
1162  nSolnValues = 2;
1163  // We are here as p goes to zero.
1164  } else if (desc > 0.0) {
1165  nSolnValues = 1;
1166  }
1167 
1168  // One real root -> have to determine whether gas or liquid is the root
1169  if (desc > 0.0) {
1170  doublereal tmpD = sqrt(desc);
1171  doublereal tmp1 = (- yN + tmpD) / (2.0 * an);
1172  doublereal sgn1 = 1.0;
1173  if (tmp1 < 0.0) {
1174  sgn1 = -1.0;
1175  tmp1 = -tmp1;
1176  }
1177  doublereal tmp2 = (- yN - tmpD) / (2.0 * an);
1178  doublereal sgn2 = 1.0;
1179  if (tmp2 < 0.0) {
1180  sgn2 = -1.0;
1181  tmp2 = -tmp2;
1182  }
1183  doublereal p1 = pow(tmp1, 1./3.);
1184  doublereal p2 = pow(tmp2, 1./3.);
1185  doublereal alpha = xN + sgn1 * p1 + sgn2 * p2;
1186  Vroot[0] = alpha;
1187  Vroot[1] = 0.0;
1188  Vroot[2] = 0.0;
1189  tmp = an * Vroot[0] * Vroot[0] * Vroot[0] + bn * Vroot[0] * Vroot[0] + cn * Vroot[0] + dn;
1190  } else if (desc < 0.0) {
1191  doublereal tmp = - yN/h;
1192  doublereal val = acos(tmp);
1193  doublereal theta = val / 3.0;
1194  doublereal oo = 2. * Pi / 3.;
1195  doublereal alpha = xN + 2. * delta * cos(theta);
1196  doublereal beta = xN + 2. * delta * cos(theta + oo);
1197  doublereal gamma = xN + 2. * delta * cos(theta + 2.0 * oo);
1198  Vroot[0] = beta;
1199  Vroot[1] = gamma;
1200  Vroot[2] = alpha;
1201 
1202  for (int i = 0; i < 3; i++) {
1203  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1204  if (fabs(tmp) > 1.0E-4) {
1205  for (int j = 0; j < 3; j++) {
1206  if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
1207  writelog("RedlichKwongMFTP::NicholsSolve(T = {}, p = {}):"
1208  " WARNING roots have merged: {}, {}\n",
1209  TKelvin, pres, Vroot[i], Vroot[j]);
1210  }
1211  }
1212  }
1213  }
1214  } else if (desc == 0.0) {
1215  if (yN == 0.0 && h == 0.0) {
1216  Vroot[0] = xN;
1217  Vroot[1] = xN;
1218  Vroot[2] = xN;
1219  } else {
1220  // need to figure out whether delta is pos or neg
1221  if (yN > 0.0) {
1222  tmp = pow(yN/(2*an), 1./3.);
1223  if (fabs(tmp - delta) > 1.0E-9) {
1224  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1225  }
1226  Vroot[1] = xN + delta;
1227  Vroot[0] = xN - 2.0*delta; // liquid phase root
1228  } else {
1229  tmp = pow(yN/(2*an), 1./3.);
1230  if (fabs(tmp - delta) > 1.0E-9) {
1231  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1232  }
1233  delta = -delta;
1234  Vroot[0] = xN + delta;
1235  Vroot[1] = xN - 2.0*delta; // gas phase root
1236  }
1237  }
1238  for (int i = 0; i < 2; i++) {
1239  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1240  }
1241  }
1242 
1243  // Unfortunately, there is a heavy amount of roundoff error due to bad
1244  // conditioning in this
1245  double res, dresdV = 0.0;
1246  for (int i = 0; i < nSolnValues; i++) {
1247  for (int n = 0; n < 20; n++) {
1248  res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1249  if (fabs(res) < 1.0E-14) {
1250  break;
1251  }
1252  dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
1253  double del = - res / dresdV;
1254  Vroot[i] += del;
1255  if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
1256  break;
1257  }
1258  double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1259  if (fabs(res2) < fabs(res)) {
1260  continue;
1261  } else {
1262  Vroot[i] -= del;
1263  Vroot[i] += 0.1 * del;
1264  }
1265  }
1266  if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
1267  writelog("RedlichKwongMFTP::NicholsSolve(T = {}, p = {}): "
1268  "WARNING root didn't converge V = {}", TKelvin, pres, Vroot[i]);
1269  writelogendl();
1270  }
1271  }
1272 
1273  if (nSolnValues == 1) {
1274  if (TKelvin > tc) {
1275  if (Vroot[0] < vc) {
1276  nSolnValues = -1;
1277  }
1278  } else {
1279  if (Vroot[0] < xN) {
1280  nSolnValues = -1;
1281  }
1282  }
1283  } else {
1284  if (nSolnValues == 2 && delta > 0.0) {
1285  nSolnValues = -2;
1286  }
1287  }
1288  return nSolnValues;
1289 }
1290 
1291 }
size_t getFloatArray(const XML_Node &node, vector_fp &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:299
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:676
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 getEnthalpy_RT_ref(doublereal *hrt) const
doublereal m_b_current
Value of b in the equation of state.
virtual void getPartialMolarIntEnergies(doublereal *ubar) const
Return an array of partial molar internal energies for the species in the mixture.
doublereal dpdV_
The derivative of the pressure wrt the volume.
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
std::string name() const
Returns the name of the XML node.
Definition: xml.h:370
virtual void getGibbs_RT_ref(doublereal *grt) const
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with &#39;v&#39;.
Definition: Array.h:112
vector_fp m_tmpV
Temporary storage - length = m_kk.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:251
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
doublereal moleFraction(size_t k) const
Return the mole fraction of a single species.
Definition: Phase.cpp:547
void pressureDerivatives() const
Calculate dpdV and dpdT at the current conditions.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture...
RedlichKwongMFTP()
Base constructor.
doublereal z() const
Calculate the value of z.
doublereal sum_xlogx() const
Evaluate .
Definition: Phase.cpp:700
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
virtual doublereal critVolume() const
Critical volume (m3/kmol).
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.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine for objects which inherit from ThermoPhase.
virtual void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional activity coefficients at the current solution temperature, pressure, and solution concentration.
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:54
STL namespace.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
void readXMLPureFluid(XML_Node &pureFluidParam)
Read the pure species RedlichKwong input parameters.
virtual doublereal refPressure(size_t k=npos) const
The reference-state pressure for species k.
virtual doublereal cp_mole() const
Molar heat capacity at constant pressure. Units: J/kmol/K.
const doublereal Pi
Pi.
Definition: ct_defs.h:51
virtual void getPartialMolarVolumes(doublereal *vbar) const
Return an array of partial molar volumes for the species in the mixture.
int NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b, doublereal Vroot[3]) const
Solve the cubic equation of state.
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture...
doublereal mean_X(const doublereal *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
Definition: Phase.cpp:690
virtual doublereal critTemperature() const
Critical temperature (K).
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:809
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
virtual bool addSpecies(shared_ptr< Species > spec)
virtual 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 getGibbs_ref(doublereal *g) const
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
virtual doublereal critDensity() const
Critical density (kg/m3).
virtual doublereal cv_mole() const
Molar heat capacity at constant volume. Units: J/kmol/K.
doublereal m_Pcurrent
Current value of the pressure.
int m_standardMixingRules
boolean indicating whether standard mixing rules are applied
static const doublereal omega_b
Omega constant for b.
int m_formTempParam
Form of the temperature parameterization.
void calculateAB(doublereal temp, doublereal &aCalc, doublereal &bCalc) const
Calculate the a and the b parameters given the temperature.
doublereal dpdT_
The derivative of the pressure wrt the temperature.
doublereal m_a_current
Value of a in the equation of state.
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
vector_fp m_pp
Temporary storage - length = m_kk.
virtual void setParametersFromXML(const XML_Node &thermoNode)
Set equation of state parameter values from XML entries.
virtual doublereal critPressure() const
Critical pressure (Pa).
const doublereal * moleFractdivMMW() const
Returns a const pointer to the start of the moleFraction/MW array.
Definition: Phase.cpp:563
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.
virtual void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials.
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
void readXMLCrossFluid(XML_Node &pureFluidParam)
Read the cross species RedlichKwong input parameters.
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
virtual doublereal enthalpy_mole() const
Molar enthalpy. Units: J/kmol.
void importPhase(XML_Node &phase, ThermoPhase *th)
Import a phase information into an empty ThermoPhase object.
virtual doublereal critCompressibility() const
Critical compressibility (unitless).
doublereal dot(InputIter x_begin, InputIter x_end, InputIter2 y_begin)
Function that calculates a templated inner product.
Definition: utilities.h:107
static const doublereal omega_a
Omega constant for a -> value of a in terms of critical properties.
virtual doublereal standardConcentration(size_t k=0) const
Returns the standard concentration , which is used to normalize the generalized concentration.
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:536
virtual void getPartialMolarEntropies(doublereal *sbar) const
Returns an array of partial molar entropies of the species in the solution.
MultiSpeciesThermo * m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1693
XML_Node & child(const size_t n) const
Return a changeable reference to the n&#39;th child of the current node.
Definition: xml.cpp:546
void writelogendl()
Write an end of line character to the screen and flush output.
Definition: global.cpp:38
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const
Calculate the pressure given the temperature and the molar volume.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:126
std::string attrib(const std::string &attr) const
Function returns the value of an attribute.
Definition: xml.cpp:500
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
virtual void getEntropy_R_ref(doublereal *er) const
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
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:157
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:661
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:637
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:130
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:149
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Definition: ThermoPhase.h:1554
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
virtual void setState_PX(doublereal p, doublereal *x)
Set the pressure (Pa) and mole fractions.
virtual doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
Contains declarations for string manipulation functions within Cantera.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
static const doublereal omega_vc
Omega constant for the critical molar volume.
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
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...
size_t m_kk
Number of species in the phase.
Definition: Phase.h:784
void applyStandardMixingRules()
Apply mixing rules for a coefficients.
virtual void getPartialMolarEnthalpies(doublereal *hbar) const
Returns an array of partial molar enthalpies for the species in the mixture.
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 entropy_mole() const
Molar entropy. Units: J/kmol/K.
Namespace for the Cantera kernel.
Definition: application.cpp:29
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
virtual int eosType() const
Equation of state type flag.
vector_fp dpdni_
Vector of derivatives of pressure wrt mole number.
size_t nChildren(bool discardComments=false) const
Return the number of children.
Definition: xml.cpp:556
void updateAB()
Update the a and b parameters.
virtual void getIntEnergy_RT(doublereal *urt) const
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
virtual void setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase.
Definition: Phase.h:622
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.
virtual void getActivityConcentrations(doublereal *c) const
This method returns an array of generalized concentrations.