Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
10 
13 #include "cantera/base/ctml.h"
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 
20 MixtureFugacityTP::MixtureFugacityTP() :
21  iState_(FLUID_GAS),
22  forcedState_(FLUID_UNDEFINED),
23  m_Tlast_ref(-1.0)
24 {
25 }
26 
28 {
30 }
31 
33 {
34  forcedState_ = solnBranch;
35 }
36 
38 {
39  return forcedState_;
40 }
41 
43 {
44  return iState_;
45 }
46 
47 // ---- Partial Molar Properties of the Solution -----------------
48 
49 void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
50 {
51  getChemPotentials(muRT);
52  for (size_t k = 0; k < m_kk; k++) {
53  muRT[k] *= 1.0 / RT();
54  }
55 }
56 
57 // ----- Thermodynamic Values for the Species Standard States States ----
58 
60 {
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(doublereal* hrt) const
70 {
71  getEnthalpy_RT_ref(hrt);
72 }
73 
74 void MixtureFugacityTP::getEntropy_R(doublereal* sr) const
75 {
77  copy(m_s0_R.begin(), m_s0_R.end(), sr);
78  double tmp = log(pressure() / refPressure());
79  for (size_t k = 0; k < m_kk; k++) {
80  sr[k] -= tmp;
81  }
82 }
83 
84 void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const
85 {
87  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
88  double tmp = log(pressure() / refPressure());
89  for (size_t k = 0; k < m_kk; k++) {
90  grt[k] += tmp;
91  }
92 }
93 
94 void MixtureFugacityTP::getPureGibbs(doublereal* g) const
95 {
97  scale(m_g0_RT.begin(), m_g0_RT.end(), g, RT());
98  double tmp = log(pressure() / refPressure()) * RT();
99  for (size_t k = 0; k < m_kk; k++) {
100  g[k] += tmp;
101  }
102 }
103 
104 void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const
105 {
107  copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
108  for (size_t i = 0; i < m_kk; i++) {
109  urt[i] -= 1.0;
110  }
111 }
112 
113 void MixtureFugacityTP::getCp_R(doublereal* cpr) const
114 {
116  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
117 }
118 
119 void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const
120 {
122  for (size_t i = 0; i < m_kk; i++) {
123  vol[i] = RT() / pressure();
124  }
125 }
126 
127 // ----- Thermodynamic Values for the Species Reference States ----
128 
129 void MixtureFugacityTP::getEnthalpy_RT_ref(doublereal* hrt) const
130 {
132  copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
133 }
134 
135 void MixtureFugacityTP::getGibbs_RT_ref(doublereal* grt) const
136 {
138  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
139 }
140 
141 void MixtureFugacityTP::getGibbs_ref(doublereal* g) const
142 {
143  const vector_fp& gibbsrt = gibbs_RT_ref();
144  scale(gibbsrt.begin(), gibbsrt.end(), g, RT());
145 }
146 
148 {
150  return m_g0_RT;
151 }
152 
153 void MixtureFugacityTP::getEntropy_R_ref(doublereal* er) const
154 {
156  copy(m_s0_R.begin(), m_s0_R.end(), er);
157 }
158 
159 void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const
160 {
162  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
163 }
164 
165 void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const
166 {
168  for (size_t i = 0; i < m_kk; i++) {
169  vol[i]= RT() / refPressure();
170  }
171 }
172 
174 {
175  int doTP = 0;
176  string comp = getChildValue(state,"moleFractions");
177  if (comp != "") {
178  // not overloaded in current object -> phase state is not calculated.
180  doTP = 1;
181  } else {
182  comp = getChildValue(state,"massFractions");
183  if (comp != "") {
184  // not overloaded in current object -> phase state is not calculated.
186  doTP = 1;
187  }
188  }
189  double t = temperature();
190  if (state.hasChild("temperature")) {
191  t = getFloat(state, "temperature", "temperature");
192  doTP = 1;
193  }
194  if (state.hasChild("pressure")) {
195  double p = getFloat(state, "pressure", "pressure");
196  setState_TP(t, p);
197  } else if (state.hasChild("density")) {
198  double rho = getFloat(state, "density", "density");
199  setState_TR(t, rho);
200  } else if (doTP) {
201  double rho = Phase::density();
202  setState_TR(t, rho);
203  }
204 }
205 
206 bool MixtureFugacityTP::addSpecies(shared_ptr<Species> spec)
207 {
208  bool added = ThermoPhase::addSpecies(spec);
209  if (added) {
210  if (m_kk == 1) {
211  moleFractions_.push_back(1.0);
212  } else {
213  moleFractions_.push_back(0.0);
214  }
215  m_h0_RT.push_back(0.0);
216  m_cp0_R.push_back(0.0);
217  m_g0_RT.push_back(0.0);
218  m_s0_R.push_back(0.0);
219  }
220  return added;
221 }
222 
223 void MixtureFugacityTP::setTemperature(const doublereal temp)
224 {
226  setState_TR(temperature(), density());
227 }
228 
230 {
231  setState_TP(temperature(), p);
232  }
233 
235 {
238 }
239 
240 void MixtureFugacityTP::setMoleFractions_NoState(const doublereal* const x)
241 {
244  updateMixingExpressions();
245 }
246 
248 {
249  throw NotImplementedError("MixtureFugacityTP::calcDensity() "
250  "called, but EOS for phase is not known");
251 }
252 
253 void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres)
254 {
255  // A pretty tricky algorithm is needed here, due to problems involving
256  // standard states of real fluids. For those cases you need to combine the T
257  // and P specification for the standard state, or else you may venture into
258  // the forbidden zone, especially when nearing the triple point. Therefore,
259  // we need to do the standard state thermo calc with the (t, pres) combo.
261 
264  // Depends on the mole fractions and the temperature
265  updateMixingExpressions();
266 
267  if (forcedState_ == FLUID_UNDEFINED) {
268  double rhoNow = Phase::density();
269  double rho = densityCalc(t, pres, iState_, rhoNow);
270  if (rho > 0.0) {
271  Phase::setDensity(rho);
272  iState_ = phaseState(true);
273  } else {
274  if (rho < -1.5) {
275  rho = densityCalc(t, pres, FLUID_UNDEFINED , rhoNow);
276  if (rho > 0.0) {
277  Phase::setDensity(rho);
278  iState_ = phaseState(true);
279  } else {
280  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
281  }
282  } else {
283  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
284  }
285  }
286  } else if (forcedState_ == FLUID_GAS) {
287  // Normal density calculation
288  if (iState_ < FLUID_LIQUID_0) {
289  double rhoNow = Phase::density();
290  double rho = densityCalc(t, pres, iState_, rhoNow);
291  if (rho > 0.0) {
292  Phase::setDensity(rho);
293  iState_ = phaseState(true);
294  if (iState_ >= FLUID_LIQUID_0) {
295  throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state");
296  }
297  } else {
298  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
299  }
300  }
301  } else if (forcedState_ > FLUID_LIQUID_0) {
302  if (iState_ >= FLUID_LIQUID_0) {
303  double rhoNow = Phase::density();
304  double rho = densityCalc(t, pres, iState_, rhoNow);
305  if (rho > 0.0) {
306  Phase::setDensity(rho);
307  iState_ = phaseState(true);
308  if (iState_ == FLUID_GAS) {
309  throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state");
310  }
311  } else {
312  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
313  }
314  }
315  }
316 }
317 
318 void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho)
319 {
323  Phase::setDensity(rho);
324  // depends on mole fraction and temperature
325  updateMixingExpressions();
326  iState_ = phaseState(true);
327 }
328 
329 void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublereal* x)
330 {
331  setMoleFractions_NoState(x);
332  setState_TP(t,p);
333 }
334 
335 doublereal MixtureFugacityTP::z() const
336 {
337  return pressure() * meanMolecularWeight() / (density() * RT());
338 }
339 
340 doublereal MixtureFugacityTP::sresid() const
341 {
342  throw CanteraError("MixtureFugacityTP::sresid()", "Base Class: not implemented");
343 }
344 
345 doublereal MixtureFugacityTP::hresid() const
346 {
347  throw CanteraError("MixtureFugacityTP::hresid()", "Base Class: not implemented");
348 }
349 
350 doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
351 {
352  doublereal pcrit = critPressure();
353  doublereal tt = critTemperature() / TKelvin;
354  if (tt < 1.0) {
355  return pcrit;
356  }
357  doublereal lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
358  return pcrit*exp(lpr);
359 }
360 
361 doublereal MixtureFugacityTP::liquidVolEst(doublereal TKelvin, doublereal& pres) const
362 {
363  throw CanteraError("MixtureFugacityTP::liquidVolEst()", "unimplemented");
364 }
365 
366 doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
367  int phase, doublereal rhoguess)
368 {
369  doublereal tcrit = critTemperature();
370  doublereal mmw = meanMolecularWeight();
371  if (rhoguess == -1.0) {
372  if (phase != -1) {
373  if (TKelvin > tcrit) {
374  rhoguess = presPa * mmw / (GasConstant * TKelvin);
375  } else {
376  if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
377  rhoguess = presPa * mmw / (GasConstant * TKelvin);
378  } else if (phase >= FLUID_LIQUID_0) {
379  double lqvol = liquidVolEst(TKelvin, presPa);
380  rhoguess = mmw / lqvol;
381  }
382  }
383  } else {
384  // Assume the Gas phase initial guess, if nothing is specified to
385  // the routine
386  rhoguess = presPa * mmw / (GasConstant * TKelvin);
387  }
388  }
389 
390  double molarVolBase = mmw / rhoguess;
391  double molarVolLast = molarVolBase;
392  double vc = mmw / critDensity();
393 
394  // molar volume of the spinodal at the current temperature and mole
395  // fractions. this will be updated as we go.
396  double molarVolSpinodal = vc;
397  bool conv = false;
398 
399  // We start on one side of the vc and stick with that side
400  bool gasSide = molarVolBase > vc;
401  if (gasSide) {
402  molarVolLast = (GasConstant * TKelvin)/presPa;
403  } else {
404  molarVolLast = liquidVolEst(TKelvin, presPa);
405  }
406 
407  // OK, now we do a small solve to calculate the molar volume given the T,P
408  // value. The algorithm is taken from dfind()
409  for (int n = 0; n < 200; n++) {
410  // Calculate the predicted reduced pressure, pred0, based on the current
411  // tau and dd. Calculate the derivative of the predicted pressure wrt
412  // the molar volume. This routine also returns the pressure, presBase
413  double presBase;
414  double dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
415 
416  // If dpdV is positive, then we are in the middle of the 2 phase region
417  // and beyond the spinodal stability curve. We need to adjust the
418  // initial guess outwards and start a new iteration.
419  if (dpdVBase >= 0.0) {
420  if (TKelvin > tcrit) {
421  throw CanteraError("MixtureFugacityTP::densityCalc",
422  "T > tcrit unexpectedly");
423  }
424 
425  // TODO Spawn a calculation for the value of the spinodal point that
426  // is very accurate. Answer the question as to whether a
427  // solution is possible on the current side of the vapor dome.
428  if (gasSide) {
429  if (molarVolBase >= vc) {
430  molarVolSpinodal = molarVolBase;
431  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
432  } else {
433  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
434  }
435  } else {
436  if (molarVolBase <= vc) {
437  molarVolSpinodal = molarVolBase;
438  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
439  } else {
440  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
441  }
442  }
443  continue;
444  }
445 
446  // Check for convergence
447  if (fabs(presBase-presPa) < 1.0E-30 + 1.0E-8 * presPa) {
448  conv = true;
449  break;
450  }
451 
452  // Dampen and crop the update
453  doublereal dpdV = dpdVBase;
454  if (n < 10) {
455  dpdV = dpdVBase * 1.5;
456  }
457 
458  // Formulate the update to the molar volume by Newton's method. Then,
459  // crop it to a max value of 0.1 times the current volume
460  double delMV = - (presBase - presPa) / dpdV;
461  if ((!gasSide || delMV < 0.0) && fabs(delMV) > 0.2 * molarVolBase) {
462  delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
463  }
464  // Only go 1/10 the way towards the spinodal at any one time.
465  if (TKelvin < tcrit) {
466  if (gasSide) {
467  if (delMV < 0.0 && -delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
468  delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
469  }
470  } else {
471  if (delMV > 0.0 && delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
472  delMV = 0.5 * (molarVolSpinodal - molarVolBase);
473  }
474  }
475  }
476  // updated the molar volume value
477  molarVolLast = molarVolBase;
478  molarVolBase += delMV;
479 
480  if (fabs(delMV/molarVolBase) < 1.0E-14) {
481  conv = true;
482  break;
483  }
484 
485  // Check for negative molar volumes
486  if (molarVolBase <= 0.0) {
487  molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
488  }
489  }
490 
491  // Check for convergence, and return 0.0 if it wasn't achieved.
492  double densBase = 0.0;
493  if (! conv) {
494  molarVolBase = 0.0;
495  throw CanteraError("MixtureFugacityTP::densityCalc()", "Process did not converge");
496  } else {
497  densBase = mmw / molarVolBase;
498  }
499  return densBase;
500 }
501 
502 void MixtureFugacityTP::updateMixingExpressions()
503 {
504 }
505 
506 int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
507  doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
508 {
509  int retn = 0;
510  doublereal densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
511  if (densLiq <= 0.0) {
512  retn = -1;
513  } else {
514  densLiqGuess = densLiq;
515  setState_TR(TKelvin, densLiq);
516  liqGRT = gibbs_mole() / RT();
517  }
518 
519  doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
520  if (densGas <= 0.0) {
521  if (retn == -1) {
522  throw CanteraError("MixtureFugacityTP::corr0",
523  "Error occurred trying to find gas density at (T,P) = {} {}",
524  TKelvin, pres);
525  }
526  retn = -2;
527  } else {
528  densGasGuess = densGas;
529  setState_TR(TKelvin, densGas);
530  gasGRT = gibbs_mole() / RT();
531  }
532  return retn;
533 }
534 
535 int MixtureFugacityTP::phaseState(bool checkState) const
536 {
537  int state = iState_;
538  if (checkState) {
539  double t = temperature();
540  double tcrit = critTemperature();
541  double rhocrit = critDensity();
542  if (t >= tcrit) {
543  return FLUID_SUPERCRIT;
544  }
545  double tmid = tcrit - 100.;
546  if (tmid < 0.0) {
547  tmid = tcrit / 2.0;
548  }
549  double pp = psatEst(tmid);
550  double mmw = meanMolecularWeight();
551  double molVolLiqTmid = liquidVolEst(tmid, pp);
552  double molVolGasTmid = GasConstant * tmid / pp;
553  double densLiqTmid = mmw / molVolLiqTmid;
554  double densGasTmid = mmw / molVolGasTmid;
555  double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
556  doublereal rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
557 
558  double rho = density();
559  int iStateGuess = FLUID_LIQUID_0;
560  if (rho < rhoMid) {
561  iStateGuess = FLUID_GAS;
562  }
563  double molarVol = mmw / rho;
564  double presCalc;
565 
566  double dpdv = dpdVCalc(t, molarVol, presCalc);
567  if (dpdv < 0.0) {
568  state = iStateGuess;
569  } else {
570  state = FLUID_UNSTABLE;
571  }
572  }
573  return state;
574 }
575 
577 {
578  throw CanteraError("MixtureFugacityTP::densSpinodalLiquid", "unimplemented");
579 }
580 
582 {
583  throw CanteraError("MixtureFugacityTP::densSpinodalGas", "unimplemented");
584 }
585 
586 doublereal MixtureFugacityTP::satPressure(doublereal TKelvin)
587 {
588  doublereal molarVolGas;
589  doublereal molarVolLiquid;
590  return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
591 }
592 
593 doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& molarVolGas,
594  doublereal& molarVolLiquid)
595 {
596  // The algorithm for this routine has undergone quite a bit of work. It
597  // probably needs more work. However, it seems now to be fairly robust. The
598  // key requirement is to find an initial pressure where both the liquid and
599  // the gas exist. This is not as easy as it sounds, and it gets exceedingly
600  // hard as the critical temperature is approached from below. Once we have
601  // this initial state, then we seek to equilibrate the Gibbs free energies
602  // of the gas and liquid and use the formula
603  //
604  // dp = VdG
605  //
606  // to create an update condition for deltaP using
607  //
608  // - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
609  //
610  // @TODO Suggestions for the future would be to switch it to an algorithm
611  // that uses the gas molar volume and the liquid molar volumes as the
612  // fundamental unknowns.
613 
614  // we need this because this is a non-const routine that is public
615  setTemperature(TKelvin);
616  double densSave = density();
617  double tempSave = temperature();
618  double pres;
619  doublereal mw = meanMolecularWeight();
620  if (TKelvin < critTemperature()) {
621  pres = psatEst(TKelvin);
622  // trial value = Psat from correlation
623  doublereal volLiquid = liquidVolEst(TKelvin, pres);
624  double RhoLiquidGood = mw / volLiquid;
625  double RhoGasGood = pres * mw / (GasConstant * TKelvin);
626  doublereal delGRT = 1.0E6;
627  doublereal liqGRT, gasGRT;
628 
629  // First part of the calculation involves finding a pressure at which
630  // the gas and the liquid state coexists.
631  doublereal presLiquid = 0.;
632  doublereal presGas;
633  doublereal presBase = pres;
634  bool foundLiquid = false;
635  bool foundGas = false;
636 
637  doublereal densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
638  if (densLiquid > 0.0) {
639  foundLiquid = true;
640  presLiquid = pres;
641  RhoLiquidGood = densLiquid;
642  }
643  if (!foundLiquid) {
644  for (int i = 0; i < 50; i++) {
645  pres = 1.1 * pres;
646  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
647  if (densLiquid > 0.0) {
648  foundLiquid = true;
649  presLiquid = pres;
650  RhoLiquidGood = densLiquid;
651  break;
652  }
653  }
654  }
655 
656  pres = presBase;
657  doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
658  if (densGas <= 0.0) {
659  foundGas = false;
660  } else {
661  foundGas = true;
662  presGas = pres;
663  RhoGasGood = densGas;
664  }
665  if (!foundGas) {
666  for (int i = 0; i < 50; i++) {
667  pres = 0.9 * pres;
668  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
669  if (densGas > 0.0) {
670  foundGas = true;
671  presGas = pres;
672  RhoGasGood = densGas;
673  break;
674  }
675  }
676  }
677 
678  if (foundGas && foundLiquid && presGas != presLiquid) {
679  pres = 0.5 * (presLiquid + presGas);
680  bool goodLiq;
681  bool goodGas;
682  for (int i = 0; i < 50; i++) {
683  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
684  if (densLiquid <= 0.0) {
685  goodLiq = false;
686  } else {
687  goodLiq = true;
688  RhoLiquidGood = densLiquid;
689  presLiquid = pres;
690  }
691  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
692  if (densGas <= 0.0) {
693  goodGas = false;
694  } else {
695  goodGas = true;
696  RhoGasGood = densGas;
697  presGas = pres;
698  }
699  if (goodGas && goodLiq) {
700  break;
701  }
702  if (!goodLiq && !goodGas) {
703  pres = 0.5 * (pres + presLiquid);
704  }
705  if (goodLiq || goodGas) {
706  pres = 0.5 * (presLiquid + presGas);
707  }
708  }
709  }
710  if (!foundGas || !foundLiquid) {
711  writelog("error couldn't find a starting pressure\n");
712  return 0.0;
713  }
714  if (presGas != presLiquid) {
715  writelog("error couldn't find a starting pressure\n");
716  return 0.0;
717  }
718 
719  pres = presGas;
720  double presLast = pres;
721  double RhoGas = RhoGasGood;
722  double RhoLiquid = RhoLiquidGood;
723 
724  // Now that we have found a good pressure we can proceed with the algorithm.
725  for (int i = 0; i < 20; i++) {
726  int stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
727  if (stab == 0) {
728  presLast = pres;
729  delGRT = liqGRT - gasGRT;
730  doublereal delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
731  doublereal dp = - delGRT * GasConstant * TKelvin / delV;
732 
733  if (fabs(dp) > 0.1 * pres) {
734  if (dp > 0.0) {
735  dp = 0.1 * pres;
736  } else {
737  dp = -0.1 * pres;
738  }
739  }
740  pres += dp;
741  } else if (stab == -1) {
742  delGRT = 1.0E6;
743  if (presLast > pres) {
744  pres = 0.5 * (presLast + pres);
745  } else {
746  // we are stuck here - try this
747  pres = 1.1 * pres;
748  }
749  } else if (stab == -2) {
750  if (presLast < pres) {
751  pres = 0.5 * (presLast + pres);
752  } else {
753  // we are stuck here - try this
754  pres = 0.9 * pres;
755  }
756  }
757  molarVolGas = mw / RhoGas;
758  molarVolLiquid = mw / RhoLiquid;
759 
760  if (fabs(delGRT) < 1.0E-8) {
761  // converged
762  break;
763  }
764  }
765 
766  molarVolGas = mw / RhoGas;
767  molarVolLiquid = mw / RhoLiquid;
768  // Put the fluid in the desired end condition
769  setState_TR(tempSave, densSave);
770  return pres;
771  } else {
772  pres = critPressure();
773  setState_TP(TKelvin, pres);
774  molarVolGas = mw / density();
775  molarVolLiquid = molarVolGas;
776  setState_TR(tempSave, densSave);
777  }
778  return pres;
779 }
780 
781 doublereal MixtureFugacityTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
782 {
783  throw CanteraError("MixtureFugacityTP::pressureCalc", "unimplemented");
784 }
785 
786 doublereal MixtureFugacityTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
787 {
788  throw CanteraError("MixtureFugacityTP::dpdVCalc", "unimplemented");
789 }
790 
792 {
793  double Tnow = temperature();
794 
795  // If the temperature has changed since the last time these
796  // properties were computed, recompute them.
797  if (m_Tlast_ref != Tnow) {
798  m_spthermo.update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
799  m_Tlast_ref = Tnow;
800 
801  // update the species Gibbs functions
802  for (size_t k = 0; k < m_kk; k++) {
803  m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
804  }
805  doublereal pref = refPressure();
806  if (pref <= 0.0) {
807  throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo()", "neg ref pressure");
808  }
809  }
810 }
811 
813 {
815  m_Tlast_ref += 0.001234;
816 }
817 
818 }
virtual int forcedSolutionBranch() const
Report the solution branch which the solution is restricted to.
virtual doublereal satPressure(doublereal TKelvin)
Calculate the saturation pressure at the current mixture content for the given temperature.
virtual void calcDensity()
Calculate the density of the mixture using the partial molar volumes and mole fractions as input...
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
virtual bool addSpecies(shared_ptr< Species > spec)
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture...
virtual doublereal densityCalc(doublereal TKelvin, doublereal pressure, int phaseRequested, doublereal rhoguess)
Calculates the density given the temperature and the pressure and a guess at the density.
virtual int standardStateConvention() const
This method returns the convention used in specification of the standard state, of which there are cu...
CTML ("Cantera Markup Language") is the variant of XML that Cantera uses to store data...
MultiSpeciesThermo m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1607
virtual void getGibbs_RT_ref(doublereal *grt) const
Returns the vector of nondimensional Gibbs Free Energies of the reference state at the current temper...
std::string getChildValue(const XML_Node &parent, const std::string &nameString)
This function reads a child node with the name, nameString, and returns its XML value as the return s...
Definition: ctml.cpp:131
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture...
An error indicating that an unimplemented function has been called.
Definition: ctexceptions.h:187
virtual void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials These are partial molar Gibbs free energ...
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual void update(doublereal T, doublereal *cp_R, doublereal *h_RT, doublereal *s_R) const
Compute the reference-state properties for all species.
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 getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:461
doublereal z() const
Calculate the value of z.
int forcedState_
Force the system to be on a particular side of the spinodal curve.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
const int cSS_CONVENTION_TEMPERATURE
Standard state uses the molar convention.
Definition: ThermoPhase.h:36
vector_fp m_g0_RT
Temporary storage for dimensionless reference state Gibbs energies.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:97
STL namespace.
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:607
doublereal calculatePsat(doublereal TKelvin, doublereal &molarVolGas, doublereal &molarVolLiquid)
Calculate the saturation pressure at the current mixture content for the given temperature.
#define FLUID_UNSTABLE
Various states of the Fugacity object.
virtual void getPureGibbs(doublereal *gpure) const
Get the pure Gibbs free energies of each species.
doublereal RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:748
virtual int reportSolnBranchActual() const
Report the solution branch which the solution is actually on.
virtual bool addSpecies(shared_ptr< Species > spec)
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
virtual void invalidateCache()
Invalidate any cached values which are normally updated only when a change in state is detected...
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 critPressure() const
Critical pressure (Pa).
Definition: ThermoPhase.h:1275
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: ThermoPhase.h:220
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
virtual void getCp_R_ref(doublereal *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual void compositionChanged()
Apply changes to the state which are needed after the composition changes.
Definition: Phase.cpp:834
void setMoleFractionsByName(const compositionMap &xMap)
Set the species mole fractions by name.
Definition: Phase.cpp:292
void setMassFractionsByName(const compositionMap &yMap)
Set the species mass fractions by name.
Definition: Phase.cpp:336
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values.
Definition: Phase.cpp:251
virtual void setForcedSolutionBranch(int solnBranch)
Set the solution branch to force the ThermoPhase to exist on one branch or another.
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:466
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:245
virtual doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
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 doublereal critTemperature() const
Critical temperature (K).
Definition: ThermoPhase.h:1270
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const
Calculate the pressure given the temperature and the molar volume.
int iState_
Current state of the fluid.
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 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
doublereal m_Tlast_ref
The last temperature at which the reference state thermodynamic properties were calculated at...
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
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
Contains declarations for string manipulation functions within Cantera.
int corr0(doublereal TKelvin, doublereal pres, doublereal &densLiq, doublereal &densGas, doublereal &liqGRT, doublereal &gasGRT)
Utility routine in the calculation of the saturation pressure.
doublereal getFloat(const XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:164
const vector_fp & gibbs_RT_ref() const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
vector_fp moleFractions_
Storage for the current values of the mole fractions of the species.
virtual void getEnthalpy_RT(doublereal *hrt) const
Get the nondimensional Enthalpy functions for the species at their standard states at the current T a...
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:788
virtual void setState_TP(doublereal T, doublereal pres)
Set the temperature (K) and pressure (Pa)
virtual void getStandardVolumes_ref(doublereal *vol) const
Get the molar volumes of the species reference states at the current T and P_ref of the solution...
virtual void getGibbs_RT(doublereal *grt) const
Get the nondimensional Gibbs functions for the species at their standard states of solution at the cu...
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
virtual doublereal critDensity() const
Critical density (kg/m3).
Definition: ThermoPhase.h:1290
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.
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species at the current T an...
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
virtual void setStateFromXML(const XML_Node &state)
Set the initial state of the phase to the conditions specified in the state XML element.
virtual doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.