Cantera  2.4.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 
9 #include "cantera/base/ctml.h"
10 
11 #include <boost/math/tools/roots.hpp>
12 
13 #include <algorithm>
14 
15 using namespace std;
16 namespace bmt = boost::math::tools;
17 
18 namespace Cantera
19 {
20 
21 const doublereal RedlichKwongMFTP::omega_a = 4.27480233540E-01;
22 const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
23 const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
24 
25 RedlichKwongMFTP::RedlichKwongMFTP() :
26  m_formTempParam(0),
27  m_b_current(0.0),
28  m_a_current(0.0),
29  NSolns_(0),
30  dpdV_(0.0),
31  dpdT_(0.0)
32 {
33  fill_n(Vroot_, 3, 0.0);
34 }
35 
36 RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string& id_) :
37  m_formTempParam(0),
38  m_b_current(0.0),
39  m_a_current(0.0),
40  NSolns_(0),
41  dpdV_(0.0),
42  dpdT_(0.0)
43 {
44  fill_n(Vroot_, 3, 0.0);
45  initThermoFile(infile, id_);
46 }
47 
48 RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id_) :
49  m_formTempParam(0),
50  m_b_current(0.0),
51  m_a_current(0.0),
52  NSolns_(0),
53  dpdV_(0.0),
54  dpdT_(0.0)
55 {
56  fill_n(Vroot_, 3, 0.0);
57  importPhase(phaseRefRoot, this);
58 }
59 
60 void RedlichKwongMFTP::setSpeciesCoeffs(const std::string& species,
61  double a0, double a1, double b)
62 {
63  size_t k = speciesIndex(species);
64  if (k == npos) {
65  throw CanteraError("RedlichKwongMFTP::setSpeciesCoeffs",
66  "Unknown species '{}'.", species);
67  }
68 
69  if (a1 != 0.0) {
70  m_formTempParam = 1; // expression is temperature-dependent
71  }
72 
73  size_t counter = k + m_kk * k;
74  a_coeff_vec(0, counter) = a0;
75  a_coeff_vec(1, counter) = a1;
76 
77  // standard mixing rule for cross-species interaction term
78  for (size_t j = 0; j < m_kk; j++) {
79  if (k == j) {
80  continue;
81  }
82  double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0);
83  double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1);
84  if (a_coeff_vec(0, j + m_kk * k) == 0) {
85  a_coeff_vec(0, j + m_kk * k) = a0kj;
86  a_coeff_vec(1, j + m_kk * k) = a1kj;
87  a_coeff_vec(0, k + m_kk * j) = a0kj;
88  a_coeff_vec(1, k + m_kk * j) = a1kj;
89  }
90  }
91  a_coeff_vec.getRow(0, a_vec_Curr_.data());
92  b_vec_Curr_[k] = b;
93 }
94 
95 void RedlichKwongMFTP::setBinaryCoeffs(const std::string& species_i,
96  const std::string& species_j, double a0, double a1)
97 {
98  size_t ki = speciesIndex(species_i);
99  if (ki == npos) {
100  throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
101  "Unknown species '{}'.", species_i);
102  }
103  size_t kj = speciesIndex(species_j);
104  if (kj == npos) {
105  throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
106  "Unknown species '{}'.", species_j);
107  }
108 
109  if (a1 != 0.0) {
110  m_formTempParam = 1; // expression is temperature-dependent
111  }
112  size_t counter1 = ki + m_kk * kj;
113  size_t counter2 = kj + m_kk * ki;
114  a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
115  a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
116  a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0;
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()/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  return pp;
177 }
178 
180 {
181  // Calculate the molarVolume of the solution (m**3 kmol-1)
182  const doublereal* const dtmp = moleFractdivMMW();
184  double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp);
185 
186  // Set the density in the parent State object directly, by calling the
187  // Phase::setDensity() function.
188  Phase::setDensity(1.0/invDens);
189 }
190 
191 void RedlichKwongMFTP::setTemperature(const doublereal temp)
192 {
193  Phase::setTemperature(temp);
195  updateAB();
196 }
197 
199 {
201  updateAB();
202 }
203 
205 {
207  for (size_t k = 0; k < m_kk; k++) {
208  c[k] *= moleFraction(k)*pressure()/RT();
209  }
210 }
211 
212 doublereal RedlichKwongMFTP::standardConcentration(size_t k) const
213 {
214  getStandardVolumes(m_tmpV.data());
215  return 1.0 / m_tmpV[k];
216 }
217 
219 {
220  doublereal mv = molarVolume();
221  doublereal sqt = sqrt(temperature());
222  doublereal vpb = mv + m_b_current;
223  doublereal vmb = mv - m_b_current;
224 
225  for (size_t k = 0; k < m_kk; k++) {
226  m_pp[k] = 0.0;
227  for (size_t i = 0; i < m_kk; i++) {
228  size_t counter = k + m_kk*i;
229  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
230  }
231  }
232  doublereal pres = pressure();
233 
234  for (size_t k = 0; k < m_kk; k++) {
235  ac[k] = (- RT() * log(pres * mv / RT())
236  + RT() * log(mv / vmb)
237  + RT() * b_vec_Curr_[k] / vmb
238  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
239  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
240  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
241  );
242  }
243  for (size_t k = 0; k < m_kk; k++) {
244  ac[k] = exp(ac[k]/RT());
245  }
246 }
247 
248 // ---- Partial Molar Properties of the Solution -----------------
249 
250 void RedlichKwongMFTP::getChemPotentials_RT(doublereal* muRT) const
251 {
252  getChemPotentials(muRT);
253  for (size_t k = 0; k < m_kk; k++) {
254  muRT[k] *= 1.0 / RT();
255  }
256 }
257 
258 void RedlichKwongMFTP::getChemPotentials(doublereal* mu) const
259 {
260  getGibbs_ref(mu);
261  for (size_t k = 0; k < m_kk; k++) {
262  double xx = std::max(SmallNumber, moleFraction(k));
263  mu[k] += RT()*(log(xx));
264  }
265 
266  doublereal mv = molarVolume();
267  doublereal sqt = sqrt(temperature());
268  doublereal vpb = mv + m_b_current;
269  doublereal vmb = mv - m_b_current;
270 
271  for (size_t k = 0; k < m_kk; k++) {
272  m_pp[k] = 0.0;
273  for (size_t i = 0; i < m_kk; i++) {
274  size_t counter = k + m_kk*i;
275  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
276  }
277  }
278  doublereal pres = pressure();
279  doublereal refP = refPressure();
280 
281  for (size_t k = 0; k < m_kk; k++) {
282  mu[k] += (RT() * log(pres/refP) - RT() * log(pres * mv / RT())
283  + RT() * log(mv / vmb)
284  + RT() * b_vec_Curr_[k] / vmb
285  - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv)
286  + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv)
287  - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb)
288  );
289  }
290 }
291 
293 {
294  // First we get the reference state contributions
295  getEnthalpy_RT_ref(hbar);
296  scale(hbar, hbar+m_kk, hbar, RT());
297 
298  // We calculate dpdni_
299  doublereal TKelvin = temperature();
300  doublereal mv = molarVolume();
301  doublereal sqt = sqrt(TKelvin);
302  doublereal vpb = mv + m_b_current;
303  doublereal vmb = mv - m_b_current;
304  for (size_t k = 0; k < m_kk; k++) {
305  m_pp[k] = 0.0;
306  for (size_t i = 0; i < m_kk; i++) {
307  size_t counter = k + m_kk*i;
308  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
309  }
310  }
311  for (size_t k = 0; k < m_kk; k++) {
312  dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb)
313  + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
314  }
315  doublereal dadt = da_dt();
316  doublereal fac = TKelvin * dadt - 3.0 * m_a_current / 2.0;
317 
318  for (size_t k = 0; k < m_kk; k++) {
319  m_tmpV[k] = 0.0;
320  for (size_t i = 0; i < m_kk; i++) {
321  size_t counter = k + m_kk*i;
322  m_tmpV[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter];
323  }
324  }
325 
327  doublereal fac2 = mv + TKelvin * dpdT_ / dpdV_;
328  for (size_t k = 0; k < m_kk; k++) {
329  double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac
330  + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_tmpV[k]
331  + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac);
332  hbar[k] = hbar[k] + hE_v;
333  hbar[k] -= fac2 * dpdni_[k];
334  }
335 }
336 
337 void RedlichKwongMFTP::getPartialMolarEntropies(doublereal* sbar) const
338 {
339  getEntropy_R_ref(sbar);
340  scale(sbar, sbar+m_kk, sbar, GasConstant);
341  doublereal TKelvin = temperature();
342  doublereal sqt = sqrt(TKelvin);
343  doublereal mv = molarVolume();
344  doublereal refP = refPressure();
345 
346  for (size_t k = 0; k < m_kk; k++) {
347  doublereal xx = std::max(SmallNumber, moleFraction(k));
348  sbar[k] += GasConstant * (- log(xx));
349  }
350  for (size_t k = 0; k < m_kk; k++) {
351  m_pp[k] = 0.0;
352  for (size_t i = 0; i < m_kk; i++) {
353  size_t counter = k + m_kk*i;
354  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
355  }
356  }
357  for (size_t k = 0; k < m_kk; k++) {
358  m_tmpV[k] = 0.0;
359  for (size_t i = 0; i < m_kk; i++) {
360  size_t counter = k + m_kk*i;
361  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
362  }
363  }
364 
365  doublereal dadt = da_dt();
366  doublereal fac = dadt - m_a_current / (2.0 * TKelvin);
367  doublereal vmb = mv - m_b_current;
368  doublereal vpb = mv + m_b_current;
369  for (size_t k = 0; k < m_kk; k++) {
370  sbar[k] -=(GasConstant * log(GasConstant * TKelvin / (refP * mv))
371  + GasConstant
372  + GasConstant * log(mv/vmb)
373  + GasConstant * b_vec_Curr_[k]/vmb
374  + m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv)
375  - 2.0 * m_tmpV[k]/(m_b_current * sqt) * log(vpb/mv)
376  + b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac
377  - 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
378  );
379  }
380 
382  getPartialMolarVolumes(m_partialMolarVolumes.data());
383  for (size_t k = 0; k < m_kk; k++) {
384  sbar[k] -= -m_partialMolarVolumes[k] * dpdT_;
385  }
386 }
387 
389 {
390  getIntEnergy_RT(ubar);
391  scale(ubar, ubar+m_kk, ubar, RT());
392 }
393 
394 void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const
395 {
396  getCp_R(cpbar);
397  scale(cpbar, cpbar+m_kk, cpbar, GasConstant);
398 }
399 
400 void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const
401 {
402  for (size_t k = 0; k < m_kk; k++) {
403  m_pp[k] = 0.0;
404  for (size_t i = 0; i < m_kk; i++) {
405  size_t counter = k + m_kk*i;
406  m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter];
407  }
408  }
409  for (size_t k = 0; k < m_kk; k++) {
410  m_tmpV[k] = 0.0;
411  for (size_t i = 0; i < m_kk; i++) {
412  size_t counter = k + m_kk*i;
413  m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1,counter);
414  }
415  }
416 
417  doublereal sqt = sqrt(temperature());
418  doublereal mv = molarVolume();
419  doublereal vmb = mv - m_b_current;
420  doublereal vpb = mv + m_b_current;
421  for (size_t k = 0; k < m_kk; k++) {
422  doublereal num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb
423  + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb)
424  - 2.0 * m_pp[k] / (sqt * vpb)
425  + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb)
426  );
427  doublereal denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb)
428  );
429  vbar[k] = num / denom;
430  }
431 }
432 
434 {
435  double pc, tc, vc;
436  double a0 = 0.0;
437  double aT = 0.0;
438  for (size_t i = 0; i < m_kk; i++) {
439  for (size_t j = 0; j <m_kk; j++) {
440  size_t counter = i + m_kk * j;
441  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
442  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
443  }
444  }
445  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
446  return tc;
447 }
448 
450 {
451  double pc, tc, vc;
452  double a0 = 0.0;
453  double aT = 0.0;
454  for (size_t i = 0; i < m_kk; i++) {
455  for (size_t j = 0; j <m_kk; j++) {
456  size_t counter = i + m_kk * j;
457  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
458  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
459  }
460  }
461  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
462  return pc;
463 }
464 
466 {
467  double pc, tc, vc;
468  double a0 = 0.0;
469  double aT = 0.0;
470  for (size_t i = 0; i < m_kk; i++) {
471  for (size_t j = 0; j <m_kk; j++) {
472  size_t counter = i + m_kk * j;
473  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
474  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
475  }
476  }
477  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
478  return vc;
479 }
480 
482 {
483  double pc, tc, vc;
484  double a0 = 0.0;
485  double aT = 0.0;
486  for (size_t i = 0; i < m_kk; i++) {
487  for (size_t j = 0; j <m_kk; j++) {
488  size_t counter = i + m_kk * j;
489  a0 += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(0, counter);
490  aT += moleFractions_[i] * moleFractions_[j] * a_coeff_vec(1, counter);
491  }
492  }
493  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
494  return pc*vc/tc/GasConstant;
495 }
496 
498 {
499  double pc, tc, vc;
500  double a0 = 0.0;
501  double aT = 0.0;
502  for (size_t i = 0; i < m_kk; i++) {
503  for (size_t j = 0; j <m_kk; j++) {
504  size_t counter = i + m_kk * j;
505  a0 += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(0, counter);
506  aT += moleFractions_[i] * moleFractions_[j] *a_coeff_vec(1, counter);
507  }
508  }
509  calcCriticalConditions(m_a_current, m_b_current, a0, aT, pc, tc, vc);
510  double mmw = meanMolecularWeight();
511  return mmw / vc;
512 }
513 
514 void RedlichKwongMFTP::setToEquilState(const doublereal* mu_RT)
515 {
516  double tmp, tmp2;
518  getGibbs_RT_ref(m_tmpV.data());
519 
520  // Within the method, we protect against inf results if the exponent is too
521  // high.
522  //
523  // If it is too low, we set the partial pressure to zero. This capability is
524  // needed by the elemental potential method.
525  doublereal pres = 0.0;
526  double m_p0 = refPressure();
527  for (size_t k = 0; k < m_kk; k++) {
528  tmp = -m_tmpV[k] + mu_RT[k];
529  if (tmp < -600.) {
530  m_pp[k] = 0.0;
531  } else if (tmp > 500.0) {
532  tmp2 = tmp / 500.;
533  tmp2 *= tmp2;
534  m_pp[k] = m_p0 * exp(500.) * tmp2;
535  } else {
536  m_pp[k] = m_p0 * exp(tmp);
537  }
538  pres += m_pp[k];
539  }
540  // set state
541  setState_PX(pres, &m_pp[0]);
542 }
543 
544 bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
545 {
546  bool added = MixtureFugacityTP::addSpecies(spec);
547  if (added) {
548  a_vec_Curr_.resize(m_kk * m_kk, 0.0);
549  b_vec_Curr_.push_back(0.0);
550 
551  a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
552 
553  m_pp.push_back(0.0);
554  m_tmpV.push_back(0.0);
555  m_partialMolarVolumes.push_back(0.0);
556  dpdni_.push_back(0.0);
557  }
558  return added;
559 }
560 
561 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
562 {
563  if (phaseNode.hasChild("thermo")) {
564  XML_Node& thermoNode = phaseNode.child("thermo");
565  std::string model = thermoNode["model"];
566  if (model != "RedlichKwong" && model != "RedlichKwongMFTP") {
567  throw CanteraError("RedlichKwongMFTP::initThermoXML",
568  "Unknown thermo model : " + model);
569  }
570 
571  // Go get all of the coefficients and factors in the
572  // activityCoefficients XML block
573  if (thermoNode.hasChild("activityCoefficients")) {
574  XML_Node& acNode = thermoNode.child("activityCoefficients");
575 
576  // Loop through the children getting multiple instances of
577  // parameters
578  for (size_t i = 0; i < acNode.nChildren(); i++) {
579  XML_Node& xmlACChild = acNode.child(i);
580  if (caseInsensitiveEquals(xmlACChild.name(), "purefluidparameters")) {
581  readXMLPureFluid(xmlACChild);
582  } else if (caseInsensitiveEquals(xmlACChild.name(), "crossfluidparameters")) {
583  readXMLCrossFluid(xmlACChild);
584  }
585  }
586  }
587  }
588 
589  MixtureFugacityTP::initThermoXML(phaseNode, id);
590 }
591 
593 {
594  string xname = pureFluidParam.name();
595  if (xname != "pureFluidParameters") {
596  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
597  "Incorrect name for processing this routine: " + xname);
598  }
599 
600  double a0 = 0.0;
601  double a1 = 0.0;
602  double b = 0.0;
603  for (size_t iChild = 0; iChild < pureFluidParam.nChildren(); iChild++) {
604  XML_Node& xmlChild = pureFluidParam.child(iChild);
605  string nodeName = toLowerCopy(xmlChild.name());
606 
607  if (nodeName == "a_coeff") {
608  vector_fp vParams;
609  string iModel = toLowerCopy(xmlChild.attrib("model"));
610  getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
611 
612  if (iModel == "constant" && vParams.size() == 1) {
613  a0 = vParams[0];
614  a1 = 0;
615  } else if (iModel == "linear_a" && vParams.size() == 2) {
616  a0 = vParams[0];
617  a1 = vParams[1];
618  } else {
619  throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
620  "unknown model or incorrect number of parameters");
621  }
622 
623  } else if (nodeName == "b_coeff") {
624  b = getFloatCurrent(xmlChild, "toSI");
625  }
626  }
627  setSpeciesCoeffs(pureFluidParam.attrib("species"), a0, a1, b);
628 }
629 
631 {
632  string xname = CrossFluidParam.name();
633  if (xname != "crossFluidParameters") {
634  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
635  "Incorrect name for processing this routine: " + xname);
636  }
637 
638  string iName = CrossFluidParam.attrib("species1");
639  string jName = CrossFluidParam.attrib("species2");
640 
641  size_t num = CrossFluidParam.nChildren();
642  for (size_t iChild = 0; iChild < num; iChild++) {
643  XML_Node& xmlChild = CrossFluidParam.child(iChild);
644  string nodeName = toLowerCopy(xmlChild.name());
645 
646  if (nodeName == "a_coeff") {
647  vector_fp vParams;
648  getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
649  string iModel = toLowerCopy(xmlChild.attrib("model"));
650  if (iModel == "constant" && vParams.size() == 1) {
651  setBinaryCoeffs(iName, jName, vParams[0], 0.0);
652  } else if (iModel == "linear_a") {
653  setBinaryCoeffs(iName, jName, vParams[0], vParams[1]);
654  } else {
655  throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
656  "unknown model ({}) or wrong number of parameters ({})",
657  iModel, vParams.size());
658  }
659  }
660  }
661 }
662 
664 {
666  std::string model = thermoNode["model"];
667 }
668 
669 doublereal RedlichKwongMFTP::sresid() const
670 {
671  // note this agrees with tpx
672  doublereal rho = density();
673  doublereal mmw = meanMolecularWeight();
674  doublereal molarV = mmw / rho;
675  double hh = m_b_current / molarV;
676  doublereal zz = z();
677  doublereal dadt = da_dt();
678  doublereal T = temperature();
679  doublereal sqT = sqrt(T);
680  doublereal fac = dadt - m_a_current / (2.0 * T);
681  double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current);
682  return GasConstant * sresid_mol_R;
683 }
684 
685 doublereal RedlichKwongMFTP::hresid() const
686 {
687  // note this agrees with tpx
688  doublereal rho = density();
689  doublereal mmw = meanMolecularWeight();
690  doublereal molarV = mmw / rho;
691  double hh = m_b_current / molarV;
692  doublereal zz = z();
693  doublereal dadt = da_dt();
694  doublereal T = temperature();
695  doublereal sqT = sqrt(T);
696  doublereal fac = T * dadt - 3.0 *m_a_current / (2.0);
697  return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current);
698 }
699 
700 doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGuess) const
701 {
702  double v = m_b_current * 1.1;
703  double atmp;
704  double btmp;
705  calculateAB(TKelvin, atmp, btmp);
706  doublereal pres = std::max(psatEst(TKelvin), presGuess);
707  double Vroot[3];
708  bool foundLiq = false;
709  int m = 0;
710  while (m < 100 && !foundLiq) {
711  int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot);
712  if (nsol == 1 || nsol == 2) {
713  double pc = critPressure();
714  if (pres > pc) {
715  foundLiq = true;
716  }
717  pres *= 1.04;
718  } else {
719  foundLiq = true;
720  }
721  }
722 
723  if (foundLiq) {
724  v = Vroot[0];
725  presGuess = pres;
726  } else {
727  v = -1.0;
728  }
729  return v;
730 }
731 
732 doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, int phaseRequested, doublereal rhoguess)
733 {
734  // It's necessary to set the temperature so that m_a_current is set correctly.
735  setTemperature(TKelvin);
736  double tcrit = critTemperature();
737  doublereal mmw = meanMolecularWeight();
738  if (rhoguess == -1.0) {
739  if (phaseRequested != FLUID_GAS) {
740  if (TKelvin > tcrit) {
741  rhoguess = presPa * mmw / (GasConstant * TKelvin);
742  } else {
743  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
744  rhoguess = presPa * mmw / (GasConstant * TKelvin);
745  } else if (phaseRequested >= FLUID_LIQUID_0) {
746  double lqvol = liquidVolEst(TKelvin, presPa);
747  rhoguess = mmw / lqvol;
748  }
749  }
750  } else {
751  // Assume the Gas phase initial guess, if nothing is specified to
752  // the routine
753  rhoguess = presPa * mmw / (GasConstant * TKelvin);
754  }
755  }
756 
757  doublereal volguess = mmw / rhoguess;
758  NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, Vroot_);
759 
760  doublereal molarVolLast = Vroot_[0];
761  if (NSolns_ >= 2) {
762  if (phaseRequested >= FLUID_LIQUID_0) {
763  molarVolLast = Vroot_[0];
764  } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
765  molarVolLast = Vroot_[2];
766  } else {
767  if (volguess > Vroot_[1]) {
768  molarVolLast = Vroot_[2];
769  } else {
770  molarVolLast = Vroot_[0];
771  }
772  }
773  } else if (NSolns_ == 1) {
774  if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
775  molarVolLast = Vroot_[0];
776  } else {
777  return -2.0;
778  }
779  } else if (NSolns_ == -1) {
780  if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
781  molarVolLast = Vroot_[0];
782  } else if (TKelvin > tcrit) {
783  molarVolLast = Vroot_[0];
784  } else {
785  return -2.0;
786  }
787  } else {
788  molarVolLast = Vroot_[0];
789  return -1.0;
790  }
791  return mmw / molarVolLast;
792 }
793 
795 {
796  double Vroot[3];
797  double T = temperature();
798  int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
799  if (nsol != 3) {
800  return critDensity();
801  }
802 
803  auto resid = [this, T](double v) {
804  double pp;
805  return dpdVCalc(T, v, pp);
806  };
807 
808  boost::uintmax_t maxiter = 100;
809  std::pair<double, double> vv = bmt::toms748_solve(
810  resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
811 
812  doublereal mmw = meanMolecularWeight();
813  return mmw / (0.5 * (vv.first + vv.second));
814 }
815 
817 {
818  double Vroot[3];
819  double T = temperature();
820  int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
821  if (nsol != 3) {
822  return critDensity();
823  }
824 
825  auto resid = [this, T](double v) {
826  double pp;
827  return dpdVCalc(T, v, pp);
828  };
829 
830  boost::uintmax_t maxiter = 100;
831  std::pair<double, double> vv = bmt::toms748_solve(
832  resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
833 
834  doublereal mmw = meanMolecularWeight();
835  return mmw / (0.5 * (vv.first + vv.second));
836 }
837 
838 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
839 {
840  doublereal sqt = sqrt(TKelvin);
841  double pres = GasConstant * TKelvin / (molarVol - m_b_current)
842  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
843  return pres;
844 }
845 
846 doublereal RedlichKwongMFTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
847 {
848  doublereal sqt = sqrt(TKelvin);
849  presCalc = GasConstant * TKelvin / (molarVol - m_b_current)
850  - m_a_current / (sqt * molarVol * (molarVol + m_b_current));
851 
852  doublereal vpb = molarVol + m_b_current;
853  doublereal vmb = molarVol - m_b_current;
854  doublereal dpdv = (- GasConstant * TKelvin / (vmb * vmb)
855  + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb));
856  return dpdv;
857 }
858 
860 {
861  doublereal TKelvin = temperature();
862  doublereal mv = molarVolume();
863  doublereal pres;
864 
865  dpdV_ = dpdVCalc(TKelvin, mv, pres);
866  doublereal sqt = sqrt(TKelvin);
867  doublereal vpb = mv + m_b_current;
868  doublereal vmb = mv - m_b_current;
869  doublereal dadt = da_dt();
870  doublereal fac = dadt - m_a_current/(2.0 * TKelvin);
871  dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb));
872 }
873 
874 void RedlichKwongMFTP::updateMixingExpressions()
875 {
876  updateAB();
877 }
878 
880 {
881  double temp = temperature();
882  if (m_formTempParam == 1) {
883  for (size_t i = 0; i < m_kk; i++) {
884  for (size_t j = 0; j < m_kk; j++) {
885  size_t counter = i * m_kk + j;
886  a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
887  }
888  }
889  }
890 
891  m_b_current = 0.0;
892  m_a_current = 0.0;
893  for (size_t i = 0; i < m_kk; i++) {
894  m_b_current += moleFractions_[i] * b_vec_Curr_[i];
895  for (size_t j = 0; j < m_kk; j++) {
896  m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j];
897  }
898  }
899 }
900 
901 void RedlichKwongMFTP::calculateAB(doublereal temp, doublereal& aCalc, doublereal& bCalc) const
902 {
903  bCalc = 0.0;
904  aCalc = 0.0;
905  if (m_formTempParam == 1) {
906  for (size_t i = 0; i < m_kk; i++) {
907  bCalc += moleFractions_[i] * b_vec_Curr_[i];
908  for (size_t j = 0; j < m_kk; j++) {
909  size_t counter = i * m_kk + j;
910  doublereal a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
911  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
912  }
913  }
914  } else {
915  for (size_t i = 0; i < m_kk; i++) {
916  bCalc += moleFractions_[i] * b_vec_Curr_[i];
917  for (size_t j = 0; j < m_kk; j++) {
918  size_t counter = i * m_kk + j;
919  doublereal a_vec_Curr = a_coeff_vec(0,counter);
920  aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j];
921  }
922  }
923  }
924 }
925 
926 doublereal RedlichKwongMFTP::da_dt() const
927 {
928  doublereal dadT = 0.0;
929  if (m_formTempParam == 1) {
930  for (size_t i = 0; i < m_kk; i++) {
931  for (size_t j = 0; j < m_kk; j++) {
932  size_t counter = i * m_kk + j;
933  dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j];
934  }
935  }
936  }
937  return dadT;
938 }
939 
940 void RedlichKwongMFTP::calcCriticalConditions(doublereal a, doublereal b, doublereal a0_coeff, doublereal aT_coeff,
941  doublereal& pc, doublereal& tc, doublereal& vc) const
942 {
943  if (m_formTempParam != 0) {
944  a = a0_coeff;
945  }
946  if (b <= 0.0) {
947  tc = 1000000.;
948  pc = 1.0E13;
949  vc = omega_vc * GasConstant * tc / pc;
950  return;
951  }
952  if (a <= 0.0) {
953  tc = 0.0;
954  pc = 0.0;
955  vc = 2.0 * b;
956  return;
957  }
958  double tmp = a * omega_b / (b * omega_a * GasConstant);
959  double pp = 2./3.;
960  doublereal sqrttc, f, dfdt, deltatc;
961 
962  if (m_formTempParam == 0) {
963  tc = pow(tmp, pp);
964  } else {
965  tc = pow(tmp, pp);
966  for (int j = 0; j < 10; j++) {
967  sqrttc = sqrt(tc);
968  f = omega_a * b * GasConstant * tc * sqrttc / omega_b - aT_coeff * tc - a0_coeff;
969  dfdt = 1.5 * omega_a * b * GasConstant * sqrttc / omega_b - aT_coeff;
970  deltatc = - f / dfdt;
971  tc += deltatc;
972  }
973  if (deltatc > 0.1) {
974  throw CanteraError("RedlichKwongMFTP::calcCriticalConditions", "didn't converge");
975  }
976  }
977 
978  pc = omega_b * GasConstant * tc / b;
979  vc = omega_vc * GasConstant * tc / pc;
980 }
981 
982 int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, doublereal b,
983  doublereal Vroot[3]) const
984 {
985  Vroot[0] = 0.0;
986  Vroot[1] = 0.0;
987  Vroot[2] = 0.0;
988  if (TKelvin <= 0.0) {
989  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "neg temperature");
990  }
991 
992  // Derive the coefficients of the cubic polynomial to solve.
993  doublereal an = 1.0;
994  doublereal bn = - GasConstant * TKelvin / pres;
995  doublereal sqt = sqrt(TKelvin);
996  doublereal cn = - (GasConstant * TKelvin * b / pres - a/(pres * sqt) + b * b);
997  doublereal dn = - (a * b / (pres * sqt));
998 
999  double tmp = a * omega_b / (b * omega_a * GasConstant);
1000  double pp = 2./3.;
1001  double tc = pow(tmp, pp);
1002  double pc = omega_b * GasConstant * tc / b;
1003  double vc = omega_vc * GasConstant * tc / pc;
1004  // Derive the center of the cubic, x_N
1005  doublereal xN = - bn /(3 * an);
1006 
1007  // Derive the value of delta**2. This is a key quantity that determines the
1008  // number of turning points
1009  doublereal delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
1010  doublereal delta = 0.0;
1011 
1012  // Calculate a couple of ratios
1013  doublereal ratio1 = 3.0 * an * cn / (bn * bn);
1014  doublereal ratio2 = pres * b / (GasConstant * TKelvin);
1015  if (fabs(ratio1) < 1.0E-7) {
1016  doublereal ratio3 = a / (GasConstant * sqt) * pres / (GasConstant * TKelvin);
1017  if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
1018  doublereal zz = 1.0;
1019  for (int i = 0; i < 10; i++) {
1020  doublereal znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
1021  doublereal deltaz = znew - zz;
1022  zz = znew;
1023  if (fabs(deltaz) < 1.0E-14) {
1024  break;
1025  }
1026  }
1027  doublereal v = zz * GasConstant * TKelvin / pres;
1028  Vroot[0] = v;
1029  return 1;
1030  }
1031  }
1032 
1033  int nSolnValues;
1034  double h2 = 4. * an * an * delta2 * delta2 * delta2;
1035  if (delta2 > 0.0) {
1036  delta = sqrt(delta2);
1037  }
1038 
1039  doublereal h = 2.0 * an * delta * delta2;
1040  doublereal yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn;
1041  doublereal desc = yN * yN - h2;
1042 
1043  if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
1044  if (desc != 0.0) {
1045  // this is for getting to other cases
1046  throw CanteraError("NicholsSolve()", "numerical issues");
1047  }
1048  desc = 0.0;
1049  }
1050 
1051  if (desc < 0.0) {
1052  nSolnValues = 3;
1053  } else if (desc == 0.0) {
1054  nSolnValues = 2;
1055  // We are here as p goes to zero.
1056  } else if (desc > 0.0) {
1057  nSolnValues = 1;
1058  }
1059 
1060  // One real root -> have to determine whether gas or liquid is the root
1061  if (desc > 0.0) {
1062  doublereal tmpD = sqrt(desc);
1063  doublereal tmp1 = (- yN + tmpD) / (2.0 * an);
1064  doublereal sgn1 = 1.0;
1065  if (tmp1 < 0.0) {
1066  sgn1 = -1.0;
1067  tmp1 = -tmp1;
1068  }
1069  doublereal tmp2 = (- yN - tmpD) / (2.0 * an);
1070  doublereal sgn2 = 1.0;
1071  if (tmp2 < 0.0) {
1072  sgn2 = -1.0;
1073  tmp2 = -tmp2;
1074  }
1075  doublereal p1 = pow(tmp1, 1./3.);
1076  doublereal p2 = pow(tmp2, 1./3.);
1077  doublereal alpha = xN + sgn1 * p1 + sgn2 * p2;
1078  Vroot[0] = alpha;
1079  Vroot[1] = 0.0;
1080  Vroot[2] = 0.0;
1081  tmp = an * Vroot[0] * Vroot[0] * Vroot[0] + bn * Vroot[0] * Vroot[0] + cn * Vroot[0] + dn;
1082  } else if (desc < 0.0) {
1083  doublereal tmp = - yN/h;
1084  doublereal val = acos(tmp);
1085  doublereal theta = val / 3.0;
1086  doublereal oo = 2. * Pi / 3.;
1087  doublereal alpha = xN + 2. * delta * cos(theta);
1088  doublereal beta = xN + 2. * delta * cos(theta + oo);
1089  doublereal gamma = xN + 2. * delta * cos(theta + 2.0 * oo);
1090  Vroot[0] = beta;
1091  Vroot[1] = gamma;
1092  Vroot[2] = alpha;
1093 
1094  for (int i = 0; i < 3; i++) {
1095  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1096  if (fabs(tmp) > 1.0E-4) {
1097  for (int j = 0; j < 3; j++) {
1098  if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
1099  writelog("RedlichKwongMFTP::NicholsSolve(T = {}, p = {}):"
1100  " WARNING roots have merged: {}, {}\n",
1101  TKelvin, pres, Vroot[i], Vroot[j]);
1102  }
1103  }
1104  }
1105  }
1106  } else if (desc == 0.0) {
1107  if (yN == 0.0 && h == 0.0) {
1108  Vroot[0] = xN;
1109  Vroot[1] = xN;
1110  Vroot[2] = xN;
1111  } else {
1112  // need to figure out whether delta is pos or neg
1113  if (yN > 0.0) {
1114  tmp = pow(yN/(2*an), 1./3.);
1115  if (fabs(tmp - delta) > 1.0E-9) {
1116  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1117  }
1118  Vroot[1] = xN + delta;
1119  Vroot[0] = xN - 2.0*delta; // liquid phase root
1120  } else {
1121  tmp = pow(yN/(2*an), 1./3.);
1122  if (fabs(tmp - delta) > 1.0E-9) {
1123  throw CanteraError("RedlichKwongMFTP::NicholsSolve()", "unexpected");
1124  }
1125  delta = -delta;
1126  Vroot[0] = xN + delta;
1127  Vroot[1] = xN - 2.0*delta; // gas phase root
1128  }
1129  }
1130  for (int i = 0; i < 2; i++) {
1131  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1132  }
1133  }
1134 
1135  // Unfortunately, there is a heavy amount of roundoff error due to bad
1136  // conditioning in this
1137  double res, dresdV = 0.0;
1138  for (int i = 0; i < nSolnValues; i++) {
1139  for (int n = 0; n < 20; n++) {
1140  res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1141  if (fabs(res) < 1.0E-14) {
1142  break;
1143  }
1144  dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn;
1145  double del = - res / dresdV;
1146  Vroot[i] += del;
1147  if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
1148  break;
1149  }
1150  double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
1151  if (fabs(res2) < fabs(res)) {
1152  continue;
1153  } else {
1154  Vroot[i] -= del;
1155  Vroot[i] += 0.1 * del;
1156  }
1157  }
1158  if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
1159  writelog("RedlichKwongMFTP::NicholsSolve(T = {}, p = {}): "
1160  "WARNING root didn't converge V = {}", TKelvin, pres, Vroot[i]);
1161  writelogendl();
1162  }
1163  }
1164 
1165  if (nSolnValues == 1) {
1166  if (TKelvin > tc) {
1167  if (Vroot[0] < vc) {
1168  nSolnValues = -1;
1169  }
1170  } else {
1171  if (Vroot[0] < xN) {
1172  nSolnValues = -1;
1173  }
1174  }
1175  } else {
1176  if (nSolnValues == 2 && delta > 0.0) {
1177  nSolnValues = -2;
1178  }
1179  }
1180  return nSolnValues;
1181 }
1182 
1183 }
void setSpeciesCoeffs(const std::string &species, double a0, double a1, double b)
Set the pure fluid interaction parameters for a species.
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:600
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 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:256
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:175
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:471
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:624
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:153
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 void getActivityCoefficients(doublereal *ac) const
Get the array of non-dimensional activity coefficients at the current solution temperature, pressure, and solution concentration.
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 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:614
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:748
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.
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.
shared_ptr< Species > species(const std::string &name) const
Return the Species object for the named species.
Definition: Phase.cpp:803
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
doublereal getFloatCurrent(const XML_Node &node, const std::string &type)
Get a floating-point value from the current XML element.
Definition: ctml.cpp:177
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:487
bool caseInsensitiveEquals(const std::string &input, const std::string &test)
Case insensitive equality predicate.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
void setBinaryCoeffs(const std::string &species_i, const std::string &species_j, double a0, double a1)
Set values for the interaction parameter between two species.
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.
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:116
virtual void setParametersFromXML(const XML_Node &eosdata)
Set equation of state parameter values from XML entries.
Definition: ThermoPhase.h:1468
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...
void getRow(size_t n, doublereal *const rw)
Get the nth row and return it in a vector.
Definition: Array.h:165
size_t m_kk
Number of species in the phase.
Definition: Phase.h:788
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 void initThermoFile(const std::string &inputFile, const std::string &id)
std::string toLowerCopy(const std::string &input)
Convert to lower case.
virtual doublereal entropy_mole() const
Molar entropy. Units: J/kmol/K.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
virtual void getPartialMolarCp(doublereal *cpbar) const
Return an array of partial molar heat capacities for the species in the mixture.
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.