Cantera  3.1.0a1
MixtureFugacityTP.cpp
Go to the documentation of this file.
1 /**
2  * @file MixtureFugacityTP.cpp
3  * Methods file for a derived class of ThermoPhase that handles
4  * non-ideal mixtures based on the fugacity models (see @ref thermoprops and
5  * class @link Cantera::MixtureFugacityTP MixtureFugacityTP@endlink).
6  */
7 
8 // This file is part of Cantera. See License.txt in the top-level directory or
9 // at https://cantera.org/license.txt for license and copyright information.
10 
13 #include "cantera/base/utilities.h"
14 #include "cantera/base/global.h"
15 
16 using namespace std;
17 
18 namespace Cantera
19 {
20 
21 int MixtureFugacityTP::standardStateConvention() const
22 {
24 }
25 
26 void MixtureFugacityTP::setForcedSolutionBranch(int solnBranch)
27 {
28  forcedState_ = solnBranch;
29 }
30 
31 int MixtureFugacityTP::forcedSolutionBranch() const
32 {
33  return forcedState_;
34 }
35 
36 int MixtureFugacityTP::reportSolnBranchActual() const
37 {
38  return iState_;
39 }
40 
41 // ---- Molar Thermodynamic Properties ---------------------------
42 double MixtureFugacityTP::enthalpy_mole() const
43 {
44  double h_ideal = RT() * mean_X(m_h0_RT);
45  double h_nonideal = hresid();
46  return h_ideal + h_nonideal;
47 }
48 
49 
50 double MixtureFugacityTP::entropy_mole() const
51 {
52  double s_ideal = GasConstant * (mean_X(m_s0_R) - sum_xlogx()
53  - std::log(pressure()/refPressure()));
54  double s_nonideal = sresid();
55  return s_ideal + s_nonideal;
56 }
57 
58 // ----- Thermodynamic Values for the Species Standard States States ----
59 
60 void MixtureFugacityTP::getStandardChemPotentials(double* g) const
61 {
62  copy(m_g0_RT.begin(), m_g0_RT.end(), g);
63  double tmp = log(pressure() / refPressure());
64  for (size_t k = 0; k < m_kk; k++) {
65  g[k] = RT() * (g[k] + tmp);
66  }
67 }
68 
69 void MixtureFugacityTP::getEnthalpy_RT(double* hrt) const
70 {
71  getEnthalpy_RT_ref(hrt);
72 }
73 
74 void MixtureFugacityTP::getEntropy_R(double* sr) const
75 {
76  copy(m_s0_R.begin(), m_s0_R.end(), sr);
77  double tmp = log(pressure() / refPressure());
78  for (size_t k = 0; k < m_kk; k++) {
79  sr[k] -= tmp;
80  }
81 }
82 
83 void MixtureFugacityTP::getGibbs_RT(double* grt) const
84 {
85  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
86  double tmp = log(pressure() / refPressure());
87  for (size_t k = 0; k < m_kk; k++) {
88  grt[k] += tmp;
89  }
90 }
91 
92 void MixtureFugacityTP::getPureGibbs(double* g) const
93 {
94  scale(m_g0_RT.begin(), m_g0_RT.end(), g, RT());
95  double tmp = log(pressure() / refPressure()) * RT();
96  for (size_t k = 0; k < m_kk; k++) {
97  g[k] += tmp;
98  }
99 }
100 
101 void MixtureFugacityTP::getIntEnergy_RT(double* urt) const
102 {
103  copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
104  for (size_t i = 0; i < m_kk; i++) {
105  urt[i] -= 1.0;
106  }
107 }
108 
109 void MixtureFugacityTP::getCp_R(double* cpr) const
110 {
111  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
112 }
113 
114 void MixtureFugacityTP::getStandardVolumes(double* vol) const
115 {
116  for (size_t i = 0; i < m_kk; i++) {
117  vol[i] = RT() / pressure();
118  }
119 }
120 
121 // ----- Thermodynamic Values for the Species Reference States ----
122 
123 void MixtureFugacityTP::getEnthalpy_RT_ref(double* hrt) const
124 {
125  copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
126 }
127 
128 void MixtureFugacityTP::getGibbs_RT_ref(double* grt) const
129 {
130  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
131 }
132 
133 void MixtureFugacityTP::getGibbs_ref(double* g) const
134 {
135  const vector<double>& gibbsrt = gibbs_RT_ref();
136  scale(gibbsrt.begin(), gibbsrt.end(), g, RT());
137 }
138 
139 const vector<double>& MixtureFugacityTP::gibbs_RT_ref() const
140 {
141  return m_g0_RT;
142 }
143 
144 void MixtureFugacityTP::getEntropy_R_ref(double* er) const
145 {
146  copy(m_s0_R.begin(), m_s0_R.end(), er);
147 }
148 
149 void MixtureFugacityTP::getCp_R_ref(double* cpr) const
150 {
151  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
152 }
153 
154 void MixtureFugacityTP::getStandardVolumes_ref(double* vol) const
155 {
156  for (size_t i = 0; i < m_kk; i++) {
157  vol[i]= RT() / refPressure();
158  }
159 }
160 
161 bool MixtureFugacityTP::addSpecies(shared_ptr<Species> spec)
162 {
163  bool added = ThermoPhase::addSpecies(spec);
164  if (added) {
165  if (m_kk == 1) {
166  moleFractions_.push_back(1.0);
167  } else {
168  moleFractions_.push_back(0.0);
169  }
170  m_h0_RT.push_back(0.0);
171  m_cp0_R.push_back(0.0);
172  m_g0_RT.push_back(0.0);
173  m_s0_R.push_back(0.0);
174  m_tmpV.push_back(0.0);
175  }
176  return added;
177 }
178 
179 void MixtureFugacityTP::setTemperature(const double temp)
180 {
181  Phase::setTemperature(temp);
182  _updateReferenceStateThermo();
183  // depends on mole fraction and temperature
184  updateMixingExpressions();
185  iState_ = phaseState(true);
186 }
187 
188 void MixtureFugacityTP::setPressure(double p)
189 {
190  // A pretty tricky algorithm is needed here, due to problems involving
191  // standard states of real fluids. For those cases you need to combine the T
192  // and P specification for the standard state, or else you may venture into
193  // the forbidden zone, especially when nearing the triple point. Therefore,
194  // we need to do the standard state thermo calc with the (t, pres) combo.
195 
196  double t = temperature();
197  double rhoNow = density();
198  if (forcedState_ == FLUID_UNDEFINED) {
199  double rho = densityCalc(t, p, iState_, rhoNow);
200  if (rho > 0.0) {
201  setDensity(rho);
202  iState_ = phaseState(true);
203  } else {
204  if (rho < -1.5) {
205  rho = densityCalc(t, p, FLUID_UNDEFINED , rhoNow);
206  if (rho > 0.0) {
207  setDensity(rho);
208  iState_ = phaseState(true);
209  } else {
210  throw CanteraError("MixtureFugacityTP::setPressure",
211  "neg rho");
212  }
213  } else {
214  throw CanteraError("MixtureFugacityTP::setPressure",
215  "neg rho");
216  }
217  }
218  } else if (forcedState_ == FLUID_GAS) {
219  // Normal density calculation
220  if (iState_ < FLUID_LIQUID_0) {
221  double rho = densityCalc(t, p, iState_, rhoNow);
222  if (rho > 0.0) {
223  setDensity(rho);
224  iState_ = phaseState(true);
225  if (iState_ >= FLUID_LIQUID_0) {
226  throw CanteraError("MixtureFugacityTP::setPressure",
227  "wrong state");
228  }
229  } else {
230  throw CanteraError("MixtureFugacityTP::setPressure",
231  "neg rho");
232  }
233  }
234  } else if (forcedState_ > FLUID_LIQUID_0) {
235  if (iState_ >= FLUID_LIQUID_0) {
236  double rho = densityCalc(t, p, iState_, rhoNow);
237  if (rho > 0.0) {
238  setDensity(rho);
239  iState_ = phaseState(true);
240  if (iState_ == FLUID_GAS) {
241  throw CanteraError("MixtureFugacityTP::setPressure",
242  "wrong state");
243  }
244  } else {
245  throw CanteraError("MixtureFugacityTP::setPressure",
246  "neg rho");
247  }
248  }
249  }
250 }
251 
252 void MixtureFugacityTP::compositionChanged()
253 {
254  Phase::compositionChanged();
255  getMoleFractions(moleFractions_.data());
256  updateMixingExpressions();
257 }
258 
259 void MixtureFugacityTP::getActivityConcentrations(double* c) const
260 {
261  getActivityCoefficients(c);
262  double p_RT = pressure() / RT();
263  for (size_t k = 0; k < m_kk; k++) {
264  c[k] *= moleFraction(k)*p_RT;
265  }
266 }
267 
268 double MixtureFugacityTP::z() const
269 {
270  return pressure() * meanMolecularWeight() / (density() * RT());
271 }
272 
273 double MixtureFugacityTP::sresid() const
274 {
275  throw NotImplementedError("MixtureFugacityTP::sresid");
276 }
277 
278 double MixtureFugacityTP::hresid() const
279 {
280  throw NotImplementedError("MixtureFugacityTP::hresid");
281 }
282 
283 double MixtureFugacityTP::psatEst(double TKelvin) const
284 {
285  double pcrit = critPressure();
286  double tt = critTemperature() / TKelvin;
287  if (tt < 1.0) {
288  return pcrit;
289  }
290  double lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
291  return pcrit*exp(lpr);
292 }
293 
294 double MixtureFugacityTP::liquidVolEst(double TKelvin, double& pres) const
295 {
296  throw NotImplementedError("MixtureFugacityTP::liquidVolEst");
297 }
298 
299 double MixtureFugacityTP::densityCalc(double TKelvin, double presPa,
300  int phase, double rhoguess)
301 {
302  double tcrit = critTemperature();
303  double mmw = meanMolecularWeight();
304  if (rhoguess == -1.0) {
305  if (phase != -1) {
306  if (TKelvin > tcrit) {
307  rhoguess = presPa * mmw / (GasConstant * TKelvin);
308  } else {
309  if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
310  rhoguess = presPa * mmw / (GasConstant * TKelvin);
311  } else if (phase >= FLUID_LIQUID_0) {
312  double lqvol = liquidVolEst(TKelvin, presPa);
313  rhoguess = mmw / lqvol;
314  }
315  }
316  } else {
317  // Assume the Gas phase initial guess, if nothing is specified to
318  // the routine
319  rhoguess = presPa * mmw / (GasConstant * TKelvin);
320  }
321  }
322 
323  double molarVolBase = mmw / rhoguess;
324  double molarVolLast = molarVolBase;
325  double vc = mmw / critDensity();
326 
327  // molar volume of the spinodal at the current temperature and mole
328  // fractions. this will be updated as we go.
329  double molarVolSpinodal = vc;
330  bool conv = false;
331 
332  // We start on one side of the vc and stick with that side
333  bool gasSide = molarVolBase > vc;
334  if (gasSide) {
335  molarVolLast = (GasConstant * TKelvin)/presPa;
336  } else {
337  molarVolLast = liquidVolEst(TKelvin, presPa);
338  }
339 
340  // OK, now we do a small solve to calculate the molar volume given the T,P
341  // value. The algorithm is taken from dfind()
342  for (int n = 0; n < 200; n++) {
343  // Calculate the predicted reduced pressure, pred0, based on the current
344  // tau and dd. Calculate the derivative of the predicted pressure wrt
345  // the molar volume. This routine also returns the pressure, presBase
346  double presBase;
347  double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
348 
349  // If dpdV is positive, then we are in the middle of the 2 phase region
350  // and beyond the spinodal stability curve. We need to adjust the
351  // initial guess outwards and start a new iteration.
352  if (dpdVBase >= 0.0) {
353  if (TKelvin > tcrit) {
354  throw CanteraError("MixtureFugacityTP::densityCalc",
355  "T > tcrit unexpectedly");
356  }
357 
358  // TODO Spawn a calculation for the value of the spinodal point that
359  // is very accurate. Answer the question as to whether a
360  // solution is possible on the current side of the vapor dome.
361  if (gasSide) {
362  if (molarVolBase >= vc) {
363  molarVolSpinodal = molarVolBase;
364  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
365  } else {
366  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
367  }
368  } else {
369  if (molarVolBase <= vc) {
370  molarVolSpinodal = molarVolBase;
371  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
372  } else {
373  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
374  }
375  }
376  continue;
377  }
378 
379  // Check for convergence
380  if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
381  conv = true;
382  break;
383  }
384 
385  // Dampen and crop the update
386  double dpdV = dpdVBase;
387  if (n < 10) {
388  dpdV = dpdVBase * 1.5;
389  }
390 
391  // Formulate the update to the molar volume by Newton's method. Then,
392  // crop it to a max value of 0.1 times the current volume
393  double delMV = - (presBase - presPa) / dpdV;
394  if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
395  delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
396  }
397  // Only go 1/10 the way towards the spinodal at any one time.
398  if (TKelvin < tcrit) {
399  if (gasSide) {
400  if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
401  delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
402  }
403  } else {
404  if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
405  delMV = 0.5 * (molarVolSpinodal - molarVolBase);
406  }
407  }
408  }
409  // updated the molar volume value
410  molarVolLast = molarVolBase;
411  molarVolBase += delMV;
412 
413  if (fabs(delMV/molarVolBase) < 1.0E-14) {
414  conv = true;
415  break;
416  }
417 
418  // Check for negative molar volumes
419  if (molarVolBase <= 0.0) {
420  molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
421  }
422  }
423 
424  // Check for convergence, and return 0.0 if it wasn't achieved.
425  double densBase = 0.0;
426  if (! conv) {
427  molarVolBase = 0.0;
428  throw CanteraError("MixtureFugacityTP::densityCalc",
429  "Process did not converge");
430  } else {
431  densBase = mmw / molarVolBase;
432  }
433  return densBase;
434 }
435 
436 void MixtureFugacityTP::updateMixingExpressions()
437 {
438 }
439 
440 int MixtureFugacityTP::corr0(double TKelvin, double pres, double& densLiqGuess,
441  double& densGasGuess, double& liqGRT, double& gasGRT)
442 {
443  int retn = 0;
444  double densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
445  if (densLiq <= 0.0) {
446  retn = -1;
447  } else {
448  densLiqGuess = densLiq;
449  setState_TD(TKelvin, densLiq);
450  liqGRT = gibbs_mole() / RT();
451  }
452 
453  double densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
454  if (densGas <= 0.0) {
455  if (retn == -1) {
456  throw CanteraError("MixtureFugacityTP::corr0",
457  "Error occurred trying to find gas density at (T,P) = {} {}",
458  TKelvin, pres);
459  }
460  retn = -2;
461  } else {
462  densGasGuess = densGas;
463  setState_TD(TKelvin, densGas);
464  gasGRT = gibbs_mole() / RT();
465  }
466  return retn;
467 }
468 
469 int MixtureFugacityTP::phaseState(bool checkState) const
470 {
471  int state = iState_;
472  if (checkState) {
473  double t = temperature();
474  double tcrit = critTemperature();
475  double rhocrit = critDensity();
476  if (t >= tcrit) {
477  return FLUID_SUPERCRIT;
478  }
479  double tmid = tcrit - 100.;
480  if (tmid < 0.0) {
481  tmid = tcrit / 2.0;
482  }
483  double pp = psatEst(tmid);
484  double mmw = meanMolecularWeight();
485  double molVolLiqTmid = liquidVolEst(tmid, pp);
486  double molVolGasTmid = GasConstant * tmid / pp;
487  double densLiqTmid = mmw / molVolLiqTmid;
488  double densGasTmid = mmw / molVolGasTmid;
489  double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
490  double rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
491 
492  double rho = density();
493  int iStateGuess = FLUID_LIQUID_0;
494  if (rho < rhoMid) {
495  iStateGuess = FLUID_GAS;
496  }
497  double molarVol = mmw / rho;
498  double presCalc;
499 
500  double dpdv = dpdVCalc(t, molarVol, presCalc);
501  if (dpdv < 0.0) {
502  state = iStateGuess;
503  } else {
504  state = FLUID_UNSTABLE;
505  }
506  }
507  return state;
508 }
509 
510 double MixtureFugacityTP::densSpinodalLiquid() const
511 {
512  throw NotImplementedError("MixtureFugacityTP::densSpinodalLiquid");
513 }
514 
515 double MixtureFugacityTP::densSpinodalGas() const
516 {
517  throw NotImplementedError("MixtureFugacityTP::densSpinodalGas");
518 }
519 
520 double MixtureFugacityTP::satPressure(double TKelvin)
521 {
522  double molarVolGas;
523  double molarVolLiquid;
524  return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
525 }
526 
527 double MixtureFugacityTP::calculatePsat(double TKelvin, double& molarVolGas,
528  double& molarVolLiquid)
529 {
530  // The algorithm for this routine has undergone quite a bit of work. It
531  // probably needs more work. However, it seems now to be fairly robust. The
532  // key requirement is to find an initial pressure where both the liquid and
533  // the gas exist. This is not as easy as it sounds, and it gets exceedingly
534  // hard as the critical temperature is approached from below. Once we have
535  // this initial state, then we seek to equilibrate the Gibbs free energies
536  // of the gas and liquid and use the formula
537  //
538  // dp = VdG
539  //
540  // to create an update condition for deltaP using
541  //
542  // - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
543  //
544  // @todo Suggestions for the future would be to switch it to an algorithm
545  // that uses the gas molar volume and the liquid molar volumes as the
546  // fundamental unknowns.
547 
548  // we need this because this is a non-const routine that is public
549  setTemperature(TKelvin);
550  double densSave = density();
551  double tempSave = temperature();
552  double pres;
553  double mw = meanMolecularWeight();
554  if (TKelvin < critTemperature()) {
555  pres = psatEst(TKelvin);
556  // trial value = Psat from correlation
557  double volLiquid = liquidVolEst(TKelvin, pres);
558  double RhoLiquidGood = mw / volLiquid;
559  double RhoGasGood = pres * mw / (GasConstant * TKelvin);
560  double delGRT = 1.0E6;
561  double liqGRT, gasGRT;
562 
563  // First part of the calculation involves finding a pressure at which
564  // the gas and the liquid state coexists.
565  double presLiquid = 0.;
566  double presGas;
567  double presBase = pres;
568  bool foundLiquid = false;
569  bool foundGas = false;
570 
571  double densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
572  if (densLiquid > 0.0) {
573  foundLiquid = true;
574  presLiquid = pres;
575  RhoLiquidGood = densLiquid;
576  }
577  if (!foundLiquid) {
578  for (int i = 0; i < 50; i++) {
579  pres = 1.1 * pres;
580  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
581  if (densLiquid > 0.0) {
582  foundLiquid = true;
583  presLiquid = pres;
584  RhoLiquidGood = densLiquid;
585  break;
586  }
587  }
588  }
589 
590  pres = presBase;
591  double densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
592  if (densGas <= 0.0) {
593  foundGas = false;
594  } else {
595  foundGas = true;
596  presGas = pres;
597  RhoGasGood = densGas;
598  }
599  if (!foundGas) {
600  for (int i = 0; i < 50; i++) {
601  pres = 0.9 * pres;
602  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
603  if (densGas > 0.0) {
604  foundGas = true;
605  presGas = pres;
606  RhoGasGood = densGas;
607  break;
608  }
609  }
610  }
611 
612  if (foundGas && foundLiquid && presGas != presLiquid) {
613  pres = 0.5 * (presLiquid + presGas);
614  bool goodLiq;
615  bool goodGas;
616  for (int i = 0; i < 50; i++) {
617  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
618  if (densLiquid <= 0.0) {
619  goodLiq = false;
620  } else {
621  goodLiq = true;
622  RhoLiquidGood = densLiquid;
623  presLiquid = pres;
624  }
625  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
626  if (densGas <= 0.0) {
627  goodGas = false;
628  } else {
629  goodGas = true;
630  RhoGasGood = densGas;
631  presGas = pres;
632  }
633  if (goodGas && goodLiq) {
634  break;
635  }
636  if (!goodLiq && !goodGas) {
637  pres = 0.5 * (pres + presLiquid);
638  }
639  if (goodLiq || goodGas) {
640  pres = 0.5 * (presLiquid + presGas);
641  }
642  }
643  }
644  if (!foundGas || !foundLiquid) {
645  warn_user("MixtureFugacityTP::calculatePsat",
646  "could not find a starting pressure; exiting.");
647  return 0.0;
648  }
649  if (presGas != presLiquid) {
650  warn_user("MixtureFugacityTP::calculatePsat",
651  "could not find a starting pressure; exiting");
652  return 0.0;
653  }
654 
655  pres = presGas;
656  double presLast = pres;
657  double RhoGas = RhoGasGood;
658  double RhoLiquid = RhoLiquidGood;
659 
660  // Now that we have found a good pressure we can proceed with the algorithm.
661  for (int i = 0; i < 20; i++) {
662  int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
663  if (stab == 0) {
664  presLast = pres;
665  delGRT = liqGRT - gasGRT;
666  double delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
667  double dp = - delGRT * GasConstant * TKelvin / delV;
668 
669  if (fabs(dp) > 0.1 * pres) {
670  if (dp > 0.0) {
671  dp = 0.1 * pres;
672  } else {
673  dp = -0.1 * pres;
674  }
675  }
676  pres += dp;
677  } else if (stab == -1) {
678  delGRT = 1.0E6;
679  if (presLast > pres) {
680  pres = 0.5 * (presLast + pres);
681  } else {
682  // we are stuck here - try this
683  pres = 1.1 * pres;
684  }
685  } else if (stab == -2) {
686  if (presLast < pres) {
687  pres = 0.5 * (presLast + pres);
688  } else {
689  // we are stuck here - try this
690  pres = 0.9 * pres;
691  }
692  }
693  molarVolGas = mw / RhoGas;
694  molarVolLiquid = mw / RhoLiquid;
695 
696  if (fabs(delGRT) < 1.0E-8) {
697  // converged
698  break;
699  }
700  }
701 
702  molarVolGas = mw / RhoGas;
703  molarVolLiquid = mw / RhoLiquid;
704  // Put the fluid in the desired end condition
705  setState_TD(tempSave, densSave);
706  return pres;
707  } else {
708  pres = critPressure();
709  setState_TP(TKelvin, pres);
710  molarVolGas = mw / density();
711  molarVolLiquid = molarVolGas;
712  setState_TD(tempSave, densSave);
713  }
714  return pres;
715 }
716 
717 double MixtureFugacityTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const
718 {
719  throw NotImplementedError("MixtureFugacityTP::dpdVCalc");
720 }
721 
722 void MixtureFugacityTP::_updateReferenceStateThermo() const
723 {
724  double Tnow = temperature();
725 
726  // If the temperature has changed since the last time these
727  // properties were computed, recompute them.
728  if (m_tlast != Tnow) {
729  m_spthermo.update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
730  m_tlast = Tnow;
731 
732  // update the species Gibbs functions
733  for (size_t k = 0; k < m_kk; k++) {
734  m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
735  }
736  double pref = refPressure();
737  if (pref <= 0.0) {
738  throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo",
739  "negative reference pressure");
740  }
741  }
742 }
743 
744 double MixtureFugacityTP::critTemperature() const
745 {
746  double pc, tc, vc;
747  calcCriticalConditions(pc, tc, vc);
748  return tc;
749 }
750 
751 double MixtureFugacityTP::critPressure() const
752 {
753  double pc, tc, vc;
754  calcCriticalConditions(pc, tc, vc);
755  return pc;
756 }
757 
758 double MixtureFugacityTP::critVolume() const
759 {
760  double pc, tc, vc;
761  calcCriticalConditions(pc, tc, vc);
762  return vc;
763 }
764 
765 double MixtureFugacityTP::critCompressibility() const
766 {
767  double pc, tc, vc;
768  calcCriticalConditions(pc, tc, vc);
769  return pc*vc/tc/GasConstant;
770 }
771 
772 double MixtureFugacityTP::critDensity() const
773 {
774  double pc, tc, vc;
775  calcCriticalConditions(pc, tc, vc);
776  double mmw = meanMolecularWeight();
777  return mmw / vc;
778 }
779 
780 void MixtureFugacityTP::calcCriticalConditions(double& pc, double& tc, double& vc) const
781 {
782  throw NotImplementedError("MixtureFugacityTP::calcCriticalConditions");
783 }
784 
785 int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b,
786  double aAlpha, double Vroot[3], double an,
787  double bn, double cn, double dn, double tc, double vc) const
788 {
789  fill_n(Vroot, 3, 0.0);
790  if (T <= 0.0) {
791  throw CanteraError("MixtureFugacityTP::solveCubic",
792  "negative temperature T = {}", T);
793  }
794 
795  // Derive the center of the cubic, x_N
796  double xN = - bn /(3 * an);
797 
798  // Derive the value of delta**2. This is a key quantity that determines the number of turning points
799  double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an);
800  double delta = 0.0;
801 
802  // Calculate a couple of ratios
803  // Cubic equation in z : z^3 - (1-B) z^2 + (A -2B -3B^2)z - (AB- B^2- B^3) = 0
804  double ratio1 = 3.0 * an * cn / (bn * bn);
805  double ratio2 = pres * b / (GasConstant * T); // B
806  if (fabs(ratio1) < 1.0E-7) {
807  double ratio3 = aAlpha / (GasConstant * T) * pres / (GasConstant * T); // A
808  if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) {
809  // A and B terms in cubic equation for z are almost zero, then z is near to 1
810  double zz = 1.0;
811  for (int i = 0; i < 10; i++) {
812  double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1);
813  double deltaz = znew - zz;
814  zz = znew;
815  if (fabs(deltaz) < 1.0E-14) {
816  break;
817  }
818  }
819  double v = zz * GasConstant * T / pres;
820  Vroot[0] = v;
821  return 1;
822  }
823  }
824 
825  int nSolnValues = -1; // Represents number of solutions to the cubic equation
826  double h2 = 4. * an * an * delta2 * delta2 * delta2; // h^2
827  if (delta2 > 0.0) {
828  delta = sqrt(delta2);
829  }
830 
831  double h = 2.0 * an * delta * delta2;
832  double yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn; // y_N term
833  double disc = yN * yN - h2; // discriminant
834 
835  //check if y = h
836  if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) {
837  if (disc > 1e-10) {
838  throw CanteraError("MixtureFugacityTP::solveCubic",
839  "value of yN and h are too high, unrealistic roots may be obtained");
840  }
841  disc = 0.0;
842  }
843 
844  if (disc < -1e-14) {
845  // disc<0 then we have three distinct roots.
846  nSolnValues = 3;
847  } else if (fabs(disc) < 1e-14) {
848  // disc=0 then we have two distinct roots (third one is repeated root)
849  nSolnValues = 2;
850  // We are here as p goes to zero.
851  } else if (disc > 1e-14) {
852  // disc> 0 then we have one real root.
853  nSolnValues = 1;
854  }
855 
856  double tmp;
857  // One real root -> have to determine whether gas or liquid is the root
858  if (disc > 0.0) {
859  double tmpD = sqrt(disc);
860  double tmp1 = (- yN + tmpD) / (2.0 * an);
861  double sgn1 = 1.0;
862  if (tmp1 < 0.0) {
863  sgn1 = -1.0;
864  tmp1 = -tmp1;
865  }
866  double tmp2 = (- yN - tmpD) / (2.0 * an);
867  double sgn2 = 1.0;
868  if (tmp2 < 0.0) {
869  sgn2 = -1.0;
870  tmp2 = -tmp2;
871  }
872  double p1 = pow(tmp1, 1./3.);
873  double p2 = pow(tmp2, 1./3.);
874  double alpha = xN + sgn1 * p1 + sgn2 * p2;
875  Vroot[0] = alpha;
876  Vroot[1] = 0.0;
877  Vroot[2] = 0.0;
878  } else if (disc < 0.0) {
879  // Three real roots alpha, beta, gamma are obtained.
880  double val = acos(-yN / h);
881  double theta = val / 3.0;
882  double twoThirdPi = 2. * Pi / 3.;
883  double alpha = xN + 2. * delta * cos(theta);
884  double beta = xN + 2. * delta * cos(theta + twoThirdPi);
885  double gamma = xN + 2. * delta * cos(theta + 2.0 * twoThirdPi);
886  Vroot[0] = beta;
887  Vroot[1] = gamma;
888  Vroot[2] = alpha;
889 
890  for (int i = 0; i < 3; i++) {
891  tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
892  if (fabs(tmp) > 1.0E-4) {
893  for (int j = 0; j < 3; j++) {
894  if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) {
895  writelog("MixtureFugacityTP::solveCubic(T ={}, p ={}):"
896  " WARNING roots have merged: {}, {}\n",
897  T, pres, Vroot[i], Vroot[j]);
898  }
899  }
900  }
901  }
902  } else if (disc == 0.0) {
903  //Three equal roots are obtained, that is, alpha = beta = gamma
904  if (yN < 1e-18 && h < 1e-18) {
905  // yN = 0.0 and h = 0 (that is, disc = 0)
906  Vroot[0] = xN;
907  Vroot[1] = xN;
908  Vroot[2] = xN;
909  } else {
910  // h and yN need to figure out whether delta^3 is positive or negative
911  if (yN > 0.0) {
912  tmp = pow(yN/(2*an), 1./3.);
913  // In this case, tmp and delta must be equal.
914  if (fabs(tmp - delta) > 1.0E-9) {
915  throw CanteraError("MixtureFugacityTP::solveCubic",
916  "Inconsistency in solver: solver is ill-conditioned.");
917  }
918  Vroot[1] = xN + delta;
919  Vroot[0] = xN - 2.0*delta; // liquid phase root
920  } else {
921  tmp = pow(yN/(2*an), 1./3.);
922  // In this case, tmp and delta must be equal.
923  if (fabs(tmp - delta) > 1.0E-9) {
924  throw CanteraError("MixtureFugacityTP::solveCubic",
925  "Inconsistency in solver: solver is ill-conditioned.");
926  }
927  delta = -delta;
928  Vroot[0] = xN + delta;
929  Vroot[1] = xN - 2.0*delta; // gas phase root
930  }
931  }
932  }
933 
934  // Find an accurate root, since there might be a heavy amount of roundoff error due to bad conditioning in this solver.
935  double res, dresdV = 0.0;
936  for (int i = 0; i < nSolnValues; i++) {
937  for (int n = 0; n < 20; n++) {
938  res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
939  if (fabs(res) < 1.0E-14) { // accurate root is obtained
940  break;
941  }
942  dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn; // derivative of the residual
943  double del = - res / dresdV;
944  Vroot[i] += del;
945  if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) {
946  break;
947  }
948  double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn;
949  if (fabs(res2) < fabs(res)) {
950  continue;
951  } else {
952  Vroot[i] -= del; // Go back to previous value of Vroot.
953  Vroot[i] += 0.1 * del; // under-relax by 0.1
954  }
955  }
956  if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) {
957  writelog("MixtureFugacityTP::solveCubic(T = {}, p = {}): "
958  "WARNING root didn't converge V = {}", T, pres, Vroot[i]);
959  writelogendl();
960  }
961  }
962 
963  if (nSolnValues == 1) {
964  // Determine the phase of the single root.
965  // nSolnValues = 1 represents the gas phase by default.
966  if (T > tc) {
967  if (Vroot[0] < vc) {
968  // Supercritical phase
969  nSolnValues = -1;
970  }
971  } else {
972  if (Vroot[0] < xN) {
973  //Liquid phase
974  nSolnValues = -1;
975  }
976  }
977  } else {
978  // Determine if we have two distinct roots or three equal roots
979  // nSolnValues = 2 represents 2 equal roots by default.
980  if (nSolnValues == 2 && delta > 1e-14) {
981  //If delta > 0, we have two distinct roots (and one repeated root)
982  nSolnValues = -2;
983  }
984  }
985  return nSolnValues;
986 }
987 
988 }
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
#define FLUID_UNSTABLE
Various states of the Fugacity object.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:195
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:175
void writelogendl()
Write an end of line character to the screen and flush output.
Definition: global.cpp:41
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:104
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
const double Pi
Pi.
Definition: ct_defs.h:68
void warn_user(const string &method, const string &msg, const Args &... args)
Print a user warning raised from method as CanteraWarning.
Definition: global.h:267
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const int cSS_CONVENTION_TEMPERATURE
Standard state uses the molar convention.
Definition: ThermoPhase.h:144
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...