Cantera  2.1.2
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  * Copyright (2005) Sandia Corporation. Under the terms of
9  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
10  * U.S. Government retains certain rights in this software.
11  */
12 
14 #include "cantera/thermo/VPSSMgr.h"
15 #include "cantera/thermo/PDSS.h"
17 #include "cantera/base/xml.h"
18 
19 using namespace std;
20 
21 namespace Cantera
22 {
23 
24 MixtureFugacityTP::MixtureFugacityTP() :
25  ThermoPhase(),
26  m_Pcurrent(-1.0),
27  moleFractions_(0),
28  iState_(FLUID_GAS),
29  forcedState_(FLUID_UNDEFINED),
30  m_Tlast_ref(-1.0),
31  m_logc0(0.0),
32  m_h0_RT(0),
33  m_cp0_R(0),
34  m_g0_RT(0),
35  m_s0_R(0)
36 {
37 }
38 
40  ThermoPhase(),
41  m_Pcurrent(-1.0),
42  moleFractions_(0),
43  iState_(FLUID_GAS),
44  forcedState_(FLUID_UNDEFINED),
45  m_Tlast_ref(-1.0),
46  m_logc0(0.0),
47  m_h0_RT(0),
48  m_cp0_R(0),
49  m_g0_RT(0),
50  m_s0_R(0)
51 {
53 }
54 
57 {
58  if (&b != this) {
59  /*
60  * Mostly, this is a passthrough to the underlying
61  * assignment operator for the ThermoPhase parent object.
62  */
64  /*
65  * However, we have to handle data that we own.
66  */
69  iState_ = b.iState_;
72  m_logc0 = b.m_logc0;
73  m_h0_RT = b.m_h0_RT;
74  m_cp0_R = b.m_cp0_R;
75  m_g0_RT = b.m_g0_RT;
76  m_s0_R = b.m_s0_R;
77  }
78  return *this;
79 }
80 
82 {
83  return new MixtureFugacityTP(*this);
84 }
85 
87 {
89 }
90 
92 {
93  forcedState_ = solnBranch;
94 }
95 
97 {
98  return forcedState_;
99 }
100 
102 {
103  return iState_;
104 }
105 
106 /*
107  * ------------Molar Thermodynamic Properties -------------------------
108  */
109 
110 doublereal MixtureFugacityTP::err(const std::string& msg) const
111 {
112  throw CanteraError("MixtureFugacityTP","Base class method "
113  +msg+" called. Equation of state type: "+int2str(eosType()));
114  return 0;
115 }
116 
117 /*
118  * ---- Partial Molar Properties of the Solution -----------------
119  */
120 
121 void MixtureFugacityTP::getChemPotentials_RT(doublereal* muRT) const
122 {
123  getChemPotentials(muRT);
124  doublereal invRT = 1.0 / _RT();
125  for (size_t k = 0; k < m_kk; k++) {
126  muRT[k] *= invRT;
127  }
128 }
129 
130 /*
131  * ----- Thermodynamic Values for the Species Standard States States ----
132  */
133 
135 {
137  copy(m_g0_RT.begin(), m_g0_RT.end(), g);
138  doublereal RT = _RT();
139  double tmp = log(pressure() /m_spthermo->refPressure());
140  for (size_t k = 0; k < m_kk; k++) {
141  g[k] = RT * (g[k] + tmp);
142  }
143 }
144 
145 void MixtureFugacityTP::getEnthalpy_RT(doublereal* hrt) const
146 {
147  getEnthalpy_RT_ref(hrt);
148 }
149 
150 #ifdef H298MODIFY_CAPABILITY
151 void MixtureFugacityTP::modifyOneHf298SS(const int k, const doublereal Hf298New)
152 {
153  m_spthermo->modifyOneHf298(k, Hf298New);
154  m_Tlast_ref += 0.0001234;
155 }
156 #endif
157 
158 void MixtureFugacityTP::getEntropy_R(doublereal* sr) const
159 {
161  copy(m_s0_R.begin(), m_s0_R.end(), sr);
162  double tmp = log(pressure() /m_spthermo->refPressure());
163  for (size_t k = 0; k < m_kk; k++) {
164  sr[k] -= tmp;
165  }
166 }
167 
168 void MixtureFugacityTP::getGibbs_RT(doublereal* grt) const
169 {
171  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
172  double tmp = log(pressure() /m_spthermo->refPressure());
173  for (size_t k = 0; k < m_kk; k++) {
174  grt[k] += tmp;
175  }
176 }
177 
178 void MixtureFugacityTP::getPureGibbs(doublereal* g) const
179 {
181  scale(m_g0_RT.begin(), m_g0_RT.end(), g, _RT());
182  double tmp = log(pressure() /m_spthermo->refPressure());
183  tmp *= _RT();
184  for (size_t k = 0; k < m_kk; k++) {
185  g[k] += tmp;
186  }
187 }
188 
189 void MixtureFugacityTP::getIntEnergy_RT(doublereal* urt) const
190 {
192  copy(m_h0_RT.begin(), m_h0_RT.end(), urt);
193  doublereal p = pressure();
194  doublereal tmp = p / _RT();
195  doublereal v0 = _RT() / p;
196  for (size_t i = 0; i < m_kk; i++) {
197  urt[i] -= tmp * v0;
198  }
199 }
200 
201 void MixtureFugacityTP::getCp_R(doublereal* cpr) const
202 {
204  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
205 }
206 
207 void MixtureFugacityTP::getStandardVolumes(doublereal* vol) const
208 {
210  doublereal v0 = _RT() / pressure();
211  for (size_t i = 0; i < m_kk; i++) {
212  vol[i]= v0;
213  }
214 }
215 
216 /*
217  * ----- Thermodynamic Values for the Species Reference States ----
218  */
219 
220 
221 void MixtureFugacityTP::getEnthalpy_RT_ref(doublereal* hrt) const
222 {
224  copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
225 }
226 
227 void MixtureFugacityTP::getGibbs_RT_ref(doublereal* grt) const
228 {
230  copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
231 }
232 
233 void MixtureFugacityTP::getGibbs_ref(doublereal* g) const
234 {
235  const vector_fp& gibbsrt = gibbs_RT_ref();
236  scale(gibbsrt.begin(), gibbsrt.end(), g, _RT());
237 }
238 
240 {
242  return m_g0_RT;
243 }
244 
245 void MixtureFugacityTP::getEntropy_R_ref(doublereal* er) const
246 {
248  copy(m_s0_R.begin(), m_s0_R.end(), er);
249  return;
250 }
251 
252 void MixtureFugacityTP::getCp_R_ref(doublereal* cpr) const
253 {
255  copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
256 }
257 
258 void MixtureFugacityTP::getStandardVolumes_ref(doublereal* vol) const
259 {
261  double pp = refPressure();
262  doublereal v0 = _RT() / pp;
263  for (size_t i = 0; i < m_kk; i++) {
264  vol[i]= v0;
265  }
266 }
267 
269 {
270  int doTP = 0;
271  string comp = ctml::getChildValue(state,"moleFractions");
272  if (comp != "") {
273  // not overloaded in current object -> phase state is not calculated.
275  doTP = 1;
276  } else {
277  comp = ctml::getChildValue(state,"massFractions");
278  if (comp != "") {
279  // not overloaded in current object -> phase state is not calculated.
281  doTP = 1;
282  }
283  }
284  double t = temperature();
285  if (state.hasChild("temperature")) {
286  t = ctml::getFloat(state, "temperature", "temperature");
287  doTP = 1;
288  }
289  if (state.hasChild("pressure")) {
290  double p = ctml::getFloat(state, "pressure", "pressure");
291  setState_TP(t, p);
292  } else if (state.hasChild("density")) {
293  double rho = ctml::getFloat(state, "density", "density");
294  setState_TR(t, rho);
295  } else if (doTP) {
296  double rho = Phase::density();
297  setState_TR(t, rho);
298  }
299 }
300 
302 {
303  initLengths();
305 }
306 
308 {
309  m_kk = nSpecies();
310  moleFractions_.resize(m_kk, 0.0);
311  moleFractions_[0] = 1.0;
312  m_h0_RT.resize(m_kk, 0.0);
313  m_cp0_R.resize(m_kk, 0.0);
314  m_g0_RT.resize(m_kk, 0.0);
315  m_s0_R.resize(m_kk, 0.0);
316 }
317 
318 void MixtureFugacityTP::setTemperature(const doublereal temp)
319 {
322 }
323 
325 {
326  setState_TP(temperature(), p);
327  // double chemPot[5], mf[5];
328  // getMoleFractions(mf);
329  // getChemPotentials(chemPot);
330  // for (int i = 0; i < m_kk; i++) {
331  // printf(" MixFug:setPres: mu(%d = %g) = %18.8g\n", i, mf[i], chemPot[i]);
332  // }
333 }
334 
335 void MixtureFugacityTP::setMassFractions(const doublereal* const y)
336 {
339 }
340 
341 void MixtureFugacityTP::setMassFractions_NoNorm(const doublereal* const y)
342 {
345 }
346 
347 void MixtureFugacityTP::setMoleFractions(const doublereal* const x)
348 {
351 }
352 
353 void MixtureFugacityTP::setMoleFractions_NoNorm(const doublereal* const x)
354 {
357 }
358 
359 void MixtureFugacityTP::setConcentrations(const doublereal* const c)
360 {
363 }
364 
365 void MixtureFugacityTP::setMoleFractions_NoState(const doublereal* const x)
366 {
369  updateMixingExpressions();
370 }
371 
373 {
374  err("MixtureFugacityTP::calcDensity() called, but EOS for phase is not known");
375 }
376 
377 void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres)
378 {
379  /*
380  * A pretty tricky algorithm is needed here, due to problems involving
381  * standard states of real fluids. For those cases you need
382  * to combine the T and P specification for the standard state, or else
383  * you may venture into the forbidden zone, especially when nearing the
384  * triple point.
385  * Therefore, we need to do the standard state thermo calc with the
386  * (t, pres) combo.
387  */
389 
390 
393  // Depends on the mole fractions and the temperature
394  updateMixingExpressions();
395  // setPressure(pres);
396  m_Pcurrent = pres;
397  // double mmw = meanMolecularWeight();
398 
399  if (forcedState_ == FLUID_UNDEFINED) {
400  double rhoNow = Phase::density();
401  double rho = densityCalc(t, pres, iState_, rhoNow);
402  if (rho > 0.0) {
403  Phase::setDensity(rho);
404  m_Pcurrent = pres;
405  iState_ = phaseState(true);
406  } else {
407  if (rho < -1.5) {
408  rho = densityCalc(t, pres, FLUID_UNDEFINED , rhoNow);
409  if (rho > 0.0) {
410  Phase::setDensity(rho);
411  m_Pcurrent = pres;
412  iState_ = phaseState(true);
413  } else {
414  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
415  }
416  } else {
417  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
418  }
419  }
420 
421 
422 
423  } else if (forcedState_ == FLUID_GAS) {
424  // Normal density calculation
425  if (iState_ < FLUID_LIQUID_0) {
426  double rhoNow = Phase::density();
427  double rho = densityCalc(t, pres, iState_, rhoNow);
428  if (rho > 0.0) {
429  Phase::setDensity(rho);
430  m_Pcurrent = pres;
431  iState_ = phaseState(true);
432  if (iState_ >= FLUID_LIQUID_0) {
433  throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state");
434  }
435  } else {
436  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
437  }
438 
439  }
440 
441 
442  } else if (forcedState_ > FLUID_LIQUID_0) {
443  if (iState_ >= FLUID_LIQUID_0) {
444  double rhoNow = Phase::density();
445  double rho = densityCalc(t, pres, iState_, rhoNow);
446  if (rho > 0.0) {
447  Phase::setDensity(rho);
448  m_Pcurrent = pres;
449  iState_ = phaseState(true);
450  if (iState_ == FLUID_GAS) {
451  throw CanteraError("MixtureFugacityTP::setState_TP()", "wrong state");
452  }
453  } else {
454  throw CanteraError("MixtureFugacityTP::setState_TP()", "neg rho");
455  }
456 
457  }
458  }
459 
460 
461 
462  //setTemperature(t);
463  //setPressure(pres);
464  //calcDensity();
465 }
466 
467 void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho)
468 {
472  Phase::setDensity(rho);
473  doublereal mv = molarVolume();
474  // depends on mole fraction and temperature
475  updateMixingExpressions();
476 
477  m_Pcurrent = pressureCalc(T, mv);
478  iState_ = phaseState(true);
479 
480  // printf("setState_TR: state at T = %g, rho = %g, mv = %g, P = %20.13g, iState = %d\n", T, rho, mv, m_Pcurrent, iState_);
481 }
482 
483 void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublereal* x)
484 {
485  setMoleFractions_NoState(x);
486  setState_TP(t,p);
487 }
488 
489 void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
490 {
492 
493  //m_VPSS_ptr->initThermo();
494 
495  // m_VPSS_ptr->initThermoXML(phaseNode, id);
496  ThermoPhase::initThermoXML(phaseNode, id_);
497 }
498 
499 doublereal MixtureFugacityTP::z() const
500 {
501  doublereal p = pressure();
502  doublereal rho = density();
503  doublereal mmw = meanMolecularWeight();
504  doublereal molarV = mmw / rho;
505  doublereal rt = _RT();
506  return p * molarV / rt;
507 }
508 
509 doublereal MixtureFugacityTP::sresid() const
510 {
511  throw CanteraError("MixtureFugacityTP::sresid()", "Base Class: not implemented");
512  return 0.0;
513 }
514 
515 doublereal MixtureFugacityTP::hresid() const
516 {
517  throw CanteraError("MixtureFugacityTP::hresid()", "Base Class: not implemented");
518  return 0.0;
519 }
520 
521 doublereal MixtureFugacityTP::psatEst(doublereal TKelvin) const
522 {
523  doublereal tcrit = critTemperature();
524  doublereal pcrit = critPressure();
525  doublereal tt = tcrit/TKelvin;
526  if (tt < 1.0) {
527  return pcrit;
528  }
529  doublereal lpr = -0.8734*tt*tt - 3.4522*tt + 4.2918;
530  return pcrit*exp(lpr);
531 }
532 
533 doublereal MixtureFugacityTP::liquidVolEst(doublereal TKelvin, doublereal& pres) const
534 {
535  throw CanteraError("MixtureFugacityTP::liquidVolEst()", "unimplemented");
536  return 0.0;
537 }
538 
539 doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa,
540  int phase, doublereal rhoguess)
541 {
542  double tcrit = critTemperature();
543  doublereal mmw = meanMolecularWeight();
544  // double pcrit = critPressure();
545  // doublereal deltaGuess = 0.0;
546  if (rhoguess == -1.0) {
547  if (phase != -1) {
548  if (TKelvin > tcrit) {
549  rhoguess = presPa * mmw / (GasConstant * TKelvin);
550  } else {
551  if (phase == FLUID_GAS || phase == FLUID_SUPERCRIT) {
552  rhoguess = presPa * mmw / (GasConstant * TKelvin);
553  } else if (phase >= FLUID_LIQUID_0) {
554  double lqvol = liquidVolEst(TKelvin, presPa);
555  rhoguess = mmw / lqvol;
556  }
557  }
558  } else {
559  /*
560  * Assume the Gas phase initial guess, if nothing is
561  * specified to the routine
562  */
563  rhoguess = presPa * mmw / (GasConstant * TKelvin);
564  }
565 
566  }
567 
568  double molarVolBase = mmw / rhoguess;
569  double molarVolLast = molarVolBase;
570  double vc = mmw / critDensity();
571  /*
572  * molar volume of the spinodal at the current temperature and mole fractions. this will
573  * be updated as we go.
574  */
575  double molarVolSpinodal = vc;
576  doublereal pcheck = 1.0E-30 + 1.0E-8 * presPa;
577  doublereal presBase, dpdVBase, delMV;
578  bool conv = false;
579  /*
580  * We start on one side of the vc and stick with that side
581  */
582  bool gasSide = molarVolBase > vc;
583  if (gasSide) {
584  molarVolLast = (GasConstant * TKelvin)/presPa;
585  } else {
586  molarVolLast = liquidVolEst(TKelvin, presPa);
587  }
588 
589  /*
590  * OK, now we do a small solve to calculate the molar volume given the T,P value.
591  * The algorithm is taken from dfind()
592  */
593  for (int n = 0; n < 200; n++) {
594 
595  /*
596  * Calculate the predicted reduced pressure, pred0, based on the
597  * current tau and dd.
598  */
599 
600  /*
601  * Calculate the derivative of the predicted pressure
602  * wrt the molar volume.
603  * This routine also returns the pressure, presBase
604  */
605  dpdVBase = dpdVCalc(TKelvin, molarVolBase, presBase);
606 
607  /*
608  * If dpdV is positve, then we are in the middle of the
609  * 2 phase region and beyond the spinodal stability curve. We need to adjust
610  * the initial guess outwards and start a new iteration.
611  */
612  if (dpdVBase >= 0.0) {
613  if (TKelvin > tcrit) {
614  throw CanteraError("", "confused");
615  }
616  /*
617  * TODO Spawn a calculation for the value of the spinodal point that is
618  * very accurate. Answer the question as to wethera solution is
619  * possible on the current side of the vapor dome.
620  */
621  if (gasSide) {
622  if (molarVolBase >= vc) {
623  molarVolSpinodal = molarVolBase;
624  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
625  } else {
626  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
627  }
628  } else {
629  if (molarVolBase <= vc) {
630  molarVolSpinodal = molarVolBase;
631  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
632  } else {
633  molarVolBase = 0.5 * (molarVolLast + molarVolSpinodal);
634  }
635  }
636  continue;
637  }
638 
639  /*
640  * Check for convergence
641  */
642  if (fabs(presBase-presPa) < pcheck) {
643  conv = true;
644  break;
645  }
646 
647  /*
648  * Dampen and crop the update
649  */
650  doublereal dpdV = dpdVBase;
651  if (n < 10) {
652  dpdV = dpdVBase * 1.5;
653  }
654  // if (dpdV > -0.001) dpdV = -0.001;
655 
656  /*
657  * Formulate the update to the molar volume by
658  * Newton's method. Then, crop it to a max value
659  * of 0.1 times the current volume
660  */
661  delMV = - (presBase - presPa) / dpdV;
662  if (!gasSide || delMV < 0.0) {
663  if (fabs(delMV) > 0.2 * molarVolBase) {
664  delMV = delMV / fabs(delMV) * 0.2 * molarVolBase;
665  }
666  }
667  /*
668  * Only go 1/10 the way towards the spinodal at any one time.
669  */
670  if (TKelvin < tcrit) {
671  if (gasSide) {
672  if (delMV < 0.0) {
673  if (-delMV > 0.5 * (molarVolBase - molarVolSpinodal)) {
674  delMV = - 0.5 * (molarVolBase - molarVolSpinodal);
675  }
676  }
677  } else {
678  if (delMV > 0.0) {
679  if (delMV > 0.5 * (molarVolSpinodal - molarVolBase)) {
680  delMV = 0.5 * (molarVolSpinodal - molarVolBase);
681  }
682  }
683  }
684  }
685  /*
686  * updated the molar volume value
687  */
688  molarVolLast = molarVolBase;
689  molarVolBase += delMV;
690 
691 
692  if (fabs(delMV/molarVolBase) < 1.0E-14) {
693  conv = true;
694  break;
695  }
696 
697  /*
698  * Check for negative molar volumes
699  */
700  if (molarVolBase <= 0.0) {
701  molarVolBase = std::min(1.0E-30, fabs(delMV*1.0E-4));
702  }
703 
704  }
705 
706 
707  /*
708  * Check for convergence, and return 0.0 if it wasn't achieved.
709  */
710  double densBase = 0.0;
711  if (! conv) {
712  molarVolBase = 0.0;
713  throw CanteraError("MixtureFugacityTP::densityCalc()", "Process didnot converge");
714  } else {
715  densBase = mmw / molarVolBase;
716  }
717  return densBase;
718 }
719 
720 void MixtureFugacityTP::updateMixingExpressions()
721 {
722 }
723 
724 MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
725  ResidEval(),
726  m_tp(tp)
727 {
728 }
729 
730 int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal* const y,
731  doublereal* const r)
732 {
733  int status = 0;
734  doublereal molarVol = y[0];
735  doublereal tt = m_tp->temperature();
736  doublereal pp;
737  doublereal val = m_tp->dpdVCalc(tt, molarVol, pp);
738  r[0] = val;
739  return status;
740 }
741 
742 int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
743  doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
744 {
745 
746  int retn = 0;
747  doublereal densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess);
748  if (densLiq <= 0.0) {
749  // throw Cantera::CanteraError("MixtureFugacityTP::corr0",
750  // "Error occurred trying to find liquid density at (T,P) = "
751  // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres));
752  retn = -1;
753  } else {
754  densLiqGuess = densLiq;
755  setState_TR(TKelvin, densLiq);
756  liqGRT = gibbs_mole() / _RT();
757  }
758 
759  doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess);
760  if (densGas <= 0.0) {
761  //throw Cantera::CanteraError("MixtureFugacityTP::corr0",
762  // "Error occurred trying to find gas density at (T,P) = "
763  // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres));
764  if (retn == -1) {
765  throw Cantera::CanteraError("MixtureFugacityTP::corr0",
766  "Error occurred trying to find gas density at (T,P) = "
767  + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres));
768  }
769  retn = -2;
770  } else {
771  densGasGuess = densGas;
772  setState_TR(TKelvin, densGas);
773  gasGRT = gibbs_mole() / _RT();
774  }
775  // delGRT = gibbsLiqRT - gibbsGasRT;
776  return retn;
777 }
778 
779 int MixtureFugacityTP::phaseState(bool checkState) const
780 {
781  int state = iState_;
782  if (checkState) {
783  double t = temperature();
784  double tcrit = critTemperature();
785  double rhocrit = critDensity();
786  if (t >= tcrit) {
787  return FLUID_SUPERCRIT;
788  }
789  double tmid = tcrit - 100.;
790  if (tmid < 0.0) {
791  tmid = tcrit / 2.0;
792  }
793  double pp = psatEst(tmid);
794  double mmw = meanMolecularWeight();
795  double molVolLiqTmid = liquidVolEst(tmid, pp);
796  double molVolGasTmid = GasConstant * tmid / (pp);
797  double densLiqTmid = mmw / molVolLiqTmid;
798  double densGasTmid = mmw / molVolGasTmid;
799  double densMidTmid = 0.5 * (densLiqTmid + densGasTmid);
800  doublereal rhoMid = rhocrit + (t - tcrit) * (rhocrit - densMidTmid) / (tcrit - tmid);
801 
802  double rho = density();
803  int iStateGuess = FLUID_LIQUID_0;
804  if (rho < rhoMid) {
805  iStateGuess = FLUID_GAS;
806  }
807  double molarVol = mmw / rho;
808  double presCalc;
809 
810  double dpdv = dpdVCalc(t, molarVol, presCalc);
811  if (dpdv < 0.0) {
812  state = iStateGuess;
813  } else {
814  state = FLUID_UNSTABLE;
815  }
816 
817  }
818  return state;
819 }
820 
822 {
823  throw CanteraError("", "unimplmented");
824  return 0.0;
825 }
826 
828 {
829  throw CanteraError("", "unimplmented");
830  return 0.0;
831 }
832 
833 doublereal MixtureFugacityTP::satPressure(doublereal TKelvin)
834 {
835  doublereal molarVolGas;
836  doublereal molarVolLiquid;
837  return calculatePsat(TKelvin, molarVolGas, molarVolLiquid);
838 }
839 
840 
841 doublereal MixtureFugacityTP::calculatePsat(doublereal TKelvin, doublereal& molarVolGas,
842  doublereal& molarVolLiquid)
843 {
844  /*
845  * The algorithm for this routine has undergone quite a bit of work. It probably needs more work.
846  * However, it seems now to be fairly robust.
847  * The key requirement is to find an initial pressure where both the liquid and the gas exist. This
848  * is not as easy as it sounds, and it gets exceedingly hard as the critical temperature is approached
849  * from below.
850  * Once we have this initial state, then we seek to equilibrate the gibbs free energies of the
851  * gas and liquid and use the formula
852  *
853  * dp = VdG
854  *
855  * to create an update condition for deltaP using
856  *
857  * - (Gliq - Ggas) = (Vliq - Vgas) (deltaP)
858  *
859  * @TODO Suggestions for the future would be to switch it to an algorithm that uses the gas molar volume
860  * and the liquid molar volumes as the fundamental unknowns.
861  */
862 
863  // we need this because this is a non-const routine that is public
864  setTemperature(TKelvin);
865  double tcrit = critTemperature();
866  double RhoLiquid, RhoGas;
867  double RhoLiquidGood, RhoGasGood;
868  double densSave = density();
869  double tempSave = temperature();
870  double pres;
871  doublereal mw = meanMolecularWeight();
872  if (TKelvin < tcrit) {
873 
874  pres = psatEst(TKelvin);
875  // trial value = Psat from correlation
876  doublereal volLiquid = liquidVolEst(TKelvin, pres);
877  RhoLiquidGood = mw / volLiquid;
878  RhoGasGood = pres * mw / (GasConstant * TKelvin);
879  doublereal delGRT = 1.0E6;
880  doublereal liqGRT, gasGRT;
881  int stab;
882  doublereal presLast = pres;
883 
884  /*
885  * First part of the calculation involves finding a pressure at which the
886  * gas and the liquid state coexists.
887  */
888  doublereal presLiquid = 0.;
889  doublereal presGas;
890  doublereal presBase = pres;
891  bool foundLiquid = false;
892  bool foundGas = false;
893 
894  doublereal densLiquid = densityCalc(TKelvin, presBase, FLUID_LIQUID_0, RhoLiquidGood);
895  if (densLiquid > 0.0) {
896  foundLiquid = true;
897  presLiquid = pres;
898  RhoLiquidGood = densLiquid;
899  }
900  if (!foundLiquid) {
901  for (int i = 0; i < 50; i++) {
902  pres = 1.1 * pres;
903  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
904  if (densLiquid > 0.0) {
905  foundLiquid = true;
906  presLiquid = pres;
907  RhoLiquidGood = densLiquid;
908  break;
909  }
910  }
911  }
912 
913  pres = presBase;
914  doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
915  if (densGas <= 0.0) {
916  foundGas = false;
917  } else {
918  foundGas = true;
919  presGas = pres;
920  RhoGasGood = densGas;
921  }
922  if (!foundGas) {
923  for (int i = 0; i < 50; i++) {
924  pres = 0.9 * pres;
925  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
926  if (densGas > 0.0) {
927  foundGas = true;
928  presGas = pres;
929  RhoGasGood = densGas;
930  break;
931  }
932  }
933  }
934 
935  if (foundGas && foundLiquid) {
936  if (presGas != presLiquid) {
937  pres = 0.5 * (presLiquid + presGas);
938  bool goodLiq;
939  bool goodGas;
940  for (int i = 0; i < 50; i++) {
941  densLiquid = densityCalc(TKelvin, pres, FLUID_LIQUID_0, RhoLiquidGood);
942  if (densLiquid <= 0.0) {
943  goodLiq = false;
944  } else {
945  goodLiq = true;
946  RhoLiquidGood = densLiquid;
947  presLiquid = pres;
948  }
949  densGas = densityCalc(TKelvin, pres, FLUID_GAS, RhoGasGood);
950  if (densGas <= 0.0) {
951  goodGas = false;
952  } else {
953  goodGas = true;
954  RhoGasGood = densGas;
955  presGas = pres;
956  }
957  if (goodGas && goodLiq) {
958  break;
959  }
960  if (!goodLiq && !goodGas) {
961  pres = 0.5 * (pres + presLiquid);
962  }
963  if (goodLiq || goodGas) {
964  pres = 0.5 * (presLiquid + presGas);
965  }
966  }
967  }
968  }
969  if (!foundGas || !foundLiquid) {
970  printf("error coundn't find a starting pressure\n");
971  return 0.0;
972  }
973  if (presGas != presLiquid) {
974  printf("error coundn't find a starting pressure\n");
975  return 0.0;
976  }
977 
978  pres = presGas;
979  presLast = pres;
980  RhoGas = RhoGasGood;
981  RhoLiquid = RhoLiquidGood;
982 
983 
984  /*
985  * Now that we have found a good pressure we can proceed with the algorithm.
986  */
987 
988  for (int i = 0; i < 20; i++) {
989 
990  stab = corr0(TKelvin, pres, RhoLiquid, RhoGas, liqGRT, gasGRT);
991  if (stab == 0) {
992  presLast = pres;
993  delGRT = liqGRT - gasGRT;
994  doublereal delV = mw * (1.0/RhoLiquid - 1.0/RhoGas);
995  doublereal dp = - delGRT * GasConstant * TKelvin / delV;
996 
997  if (fabs(dp) > 0.1 * pres) {
998  if (dp > 0.0) {
999  dp = 0.1 * pres;
1000  } else {
1001  dp = -0.1 * pres;
1002  }
1003  }
1004  pres += dp;
1005 
1006  } else if (stab == -1) {
1007  delGRT = 1.0E6;
1008  if (presLast > pres) {
1009  pres = 0.5 * (presLast + pres);
1010  } else {
1011  // we are stuck here - try this
1012  pres = 1.1 * pres;
1013  }
1014  } else if (stab == -2) {
1015  if (presLast < pres) {
1016  pres = 0.5 * (presLast + pres);
1017  } else {
1018  // we are stuck here - try this
1019  pres = 0.9 * pres;
1020  }
1021  }
1022  molarVolGas = mw / RhoGas;
1023  molarVolLiquid = mw / RhoLiquid;
1024 
1025 
1026  if (fabs(delGRT) < 1.0E-8) {
1027  // converged
1028  break;
1029  }
1030  }
1031 
1032  molarVolGas = mw / RhoGas;
1033  molarVolLiquid = mw / RhoLiquid;
1034  // Put the fluid in the desired end condition
1035  setState_TR(tempSave, densSave);
1036 
1037  return pres;
1038 
1039 
1040  } else {
1041  pres = critPressure();
1042  setState_TP(TKelvin, pres);
1043  RhoGas = density();
1044  molarVolGas = mw / RhoGas;
1045  molarVolLiquid = molarVolGas;
1046  setState_TR(tempSave, densSave);
1047  }
1048  return pres;
1049 }
1050 
1051 doublereal MixtureFugacityTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const
1052 {
1053  throw CanteraError("MixtureFugacityTP::pressureCalc", "unimplemented");
1054  return 0.0;
1055 }
1056 
1057 doublereal MixtureFugacityTP::dpdVCalc(doublereal TKelvin, doublereal molarVol, doublereal& presCalc) const
1058 {
1059  throw CanteraError("MixtureFugacityTP::dpdVCalc", "unimplemented");
1060  return 0.0;
1061 }
1062 
1064 {
1065  double Tnow = temperature();
1066 
1067  // If the temperature has changed since the last time these
1068  // properties were computed, recompute them.
1069  if (m_Tlast_ref != Tnow) {
1070  m_spthermo->update(Tnow, &m_cp0_R[0], &m_h0_RT[0], &m_s0_R[0]);
1071  m_Tlast_ref = Tnow;
1072 
1073  // update the species Gibbs functions
1074  for (size_t k = 0; k < m_kk; k++) {
1075  m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
1076  }
1077  doublereal pref = refPressure();
1078  if (pref <= 0.0) {
1079  throw CanteraError("MixtureFugacityTP::_updateReferenceStateThermo()", "neg ref pressure");
1080  }
1081  m_logc0 = log(pref/(GasConstant * Tnow));
1082  }
1083 }
1084 
1085 }
virtual void setMoleFractions_NoNorm(const doublereal *const x)
Set the mole fractions to the specified values without normalizing.
Definition: Phase.cpp:344
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 doublereal densSpinodalGas() const
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
virtual doublereal density() const
Density (kg/m^3).
Definition: Phase.h:534
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Initialize a ThermoPhase object, potentially reading activity coefficient information from an XML dat...
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 doublereal liquidVolEst(doublereal TKelvin, doublereal &pres) const
Estimate for the molar volume of the liquid.
virtual doublereal gibbs_mole() const
Molar Gibbs function. Units: J/kmol.
Definition: ThermoPhase.h:279
virtual int eosType() const
Equation of state type flag.
virtual void getGibbs_ref(doublereal *g) const
Declaration file for a virtual base class that manages the calculation of standard state properties f...
doublereal _RT() const
Return the Gas Constant multiplied by the current temperature.
Definition: ThermoPhase.h:977
virtual void getStandardChemPotentials(doublereal *mu) const
Get the array of chemical potentials at unit activity.
virtual void initThermo()
Initialize the ThermoPhase object after all species have been set up.
ThermoPhase & operator=(const ThermoPhase &right)
Assignment operator.
Definition: ThermoPhase.cpp:57
int forcedState_
Force the system to be on a particular side of the spinodal curve.
void getPureGibbs(doublereal *gpure) const
Get the pure Gibbs free energies of each species.
virtual void modifyOneHf298SS(const int k, const doublereal Hf298New)
Modify the value of the 298 K Heat of Formation of one species in the phase (J kmol-1) ...
Definition: ThermoPhase.h:227
doublereal z() const
Calculate the value of z.
virtual void getCp_R(doublereal *cpr) const
Get the nondimensional Heat Capacities at constant pressure for the standard state of the species at ...
const int cSS_CONVENTION_TEMPERATURE
Standard state uses the molar convention.
Definition: ThermoPhase.h:34
vector_fp m_g0_RT
Temporary storage for dimensionless reference state gibbs energies.
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values, and then normalize them so that they sum to 1...
doublereal err(const std::string &msg) const
MixtureFugacityTP has its own err routine.
Class XML_Node is a tree-based representation of the contents of an XML file.
Definition: xml.h:100
This is a filter class for ThermoPhase that implements some preparatory steps for efficiently handlin...
MixtureFugacityTP & operator=(const MixtureFugacityTP &b)
Assignment operator.
doublereal getFloat(const Cantera::XML_Node &parent, const std::string &name, const std::string &type)
Get a floating-point value from a child element.
Definition: ctml.cpp:267
virtual int standardStateConvention() const
This method returns the convention used in specification of the standard state, of which there are cu...
doublereal calculatePsat(doublereal TKelvin, doublereal &molarVolGas, doublereal &molarVolLiquid)
Calculate the saturation pressure at the current mixture content for the given temperature.
virtual doublereal critTemperature() const
Critical temperature (K).
Definition: ThermoPhase.h:1241
virtual void getGibbs_RT_ref(doublereal *grt) const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
doublereal pressure() const
Returns the current pressure of the phase.
virtual void getStandardVolumes_ref(doublereal *vol) const
Get the molar volumes of the species reference states at the current T and reference pressure of the ...
#define FLUID_UNSTABLE
Various states of the Fugacity object.
virtual void getEnthalpy_RT_ref(doublereal *hrt) const
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
virtual void getGibbs_RT(doublereal *grt) const
Get the nondimensional Gibbs functions for the species at their standard states of solution at the cu...
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
virtual void setState_TR(doublereal T, doublereal rho)
Set the internally stored temperature (K) and density (kg/m^3)
vector_fp m_s0_R
Temporary storage for dimensionless reference state entropies.
virtual void getStandardVolumes(doublereal *vol) const
Get the molar volumes of each species in their standard states at the current T and P of the solution...
void setMassFractionsByName(compositionMap &yMap)
Set the species mass fractions by name.
Definition: Phase.cpp:399
doublereal m_Pcurrent
Current value of the pressures.
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
virtual void setMassFractions(const doublereal *const y)
Set the mass fractions to the specified values, and then normalize them so that they sum to 1...
virtual void setTemperature(const doublereal temp)
Set the temperature of the phase.
Declarations for the virtual base class PDSS (pressure dependent standard state) which handles calcul...
virtual doublereal psatEst(doublereal TKelvin) const
Estimate for the saturation pressure.
virtual doublereal densSpinodalLiquid() const
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
virtual doublereal critDensity() const
Critical density (kg/m3).
Definition: ThermoPhase.h:1253
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
doublereal molarVolume() const
Molar volume (m^3/kmol).
Definition: Phase.cpp:607
virtual void setConcentrations(const doublereal *const c)
Set the concentrations to the specified values within the phase.
std::string fp2str(const double x, const std::string &fmt)
Convert a double into a c++ string.
Definition: stringUtils.cpp:29
doublereal m_logc0
Temporary storage for log of p/rt.
virtual void getChemPotentials(doublereal *mu) const
Get the species chemical potentials. Units: J/kmol.
Definition: ThermoPhase.h:595
Classes providing support for XML data files.
const vector_fp & gibbs_RT_ref() const
Returns the vector of nondimensional Gibbs free energies of the reference state at the current temper...
std::vector< doublereal > moleFractions_
Storage for the current values of the mole fractions of the species.
virtual doublereal refPressure(size_t k=npos) const =0
The reference-state pressure for species k.
virtual doublereal pressureCalc(doublereal TKelvin, doublereal molarVol) const
Calculate the pressure given the temperature and the molar volume.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
void setMoleFractionsByName(compositionMap &xMap)
Set the species mole fractions by name.
Definition: Phase.cpp:354
virtual void setConcentrations(const doublereal *const conc)
Set the concentrations to the specified values within the phase.
Definition: Phase.cpp:577
virtual void setMoleFractions(const doublereal *const x)
Set the mole fractions to the specified values There is no restriction on the sum of the mole fractio...
Definition: Phase.cpp:306
void getChemPotentials_RT(doublereal *mu) const
Get the array of non-dimensional species chemical potentials These are partial molar Gibbs free energ...
bool hasChild(const std::string &ch) const
Tests whether the current node has a child node with a particular name.
Definition: xml.cpp:574
virtual doublereal hresid() const
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture...
virtual void setForcedSolutionBranch(int solnBranch)
Set the solution branch to force the ThermoPhase to exist on one branch or another.
int phaseState(bool checkState=false) const
Returns the Phase State flag for the current state of the object.
virtual int forcedSolutionBranch() const
Report the solution branch which the solution is restricted to.
virtual void setPressure(doublereal p)
Set the internally stored pressure (Pa) at constant temperature and composition.
virtual void getIntEnergy_RT(doublereal *urt) const
Returns the vector of nondimensional internal Energies of the standard state at the current temperatu...
virtual void setMoleFractions_NoNorm(const doublereal *const x)
Set the mole fractions to the specified values without normalizing.
virtual doublereal refPressure() const
Returns the reference pressure in Pa.
Definition: ThermoPhase.h:159
virtual void setMassFractions_NoNorm(const doublereal *const y)
Set the mass fractions to the specified values without normalizing.
Definition: Phase.cpp:388
virtual void getEntropy_R(doublereal *sr) const
Get the array of nondimensional Enthalpy functions for the standard state species.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
virtual doublereal critPressure() const
Critical pressure (Pa).
Definition: ThermoPhase.h:1247
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
virtual ThermoPhase * duplMyselfAsThermoPhase() const
Duplication routine.
virtual void initThermoXML(XML_Node &phaseNode, const std::string &id)
Import and initialize a ThermoPhase object using an XML tree.
int iState_
Current state of the fluid.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
virtual void setTemperature(const doublereal temp)
Set the internally stored temperature of the phase (K).
Definition: Phase.h:562
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:154
doublereal m_Tlast_ref
The last temperature at which the reference state thermodynamic properties were calculated at...
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:588
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
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.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
virtual void setMassFractions(const doublereal *const y)
Set the mass fractions to the specified values and normalize them.
Definition: Phase.cpp:374
size_t m_kk
Number of species in the phase.
Definition: Phase.h:716
virtual void setState_TP(doublereal T, doublereal pres)
Set the temperature and pressure at the same time.
virtual void getEntropy_R_ref(doublereal *er) const
virtual void update(doublereal T, doublereal *cp_R, doublereal *h_RT, doublereal *s_R) const =0
Compute the reference-state properties for all species.
virtual doublereal sresid() const
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture...
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
std::string getChildValue(const Cantera::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:164
virtual int reportSolnBranchActual() const
Report the solution branch which the solution is actually on.
SpeciesThermo * m_spthermo
Pointer to the calculation manager for species reference-state thermodynamic properties.
Definition: ThermoPhase.h:1625
virtual void getCp_R_ref(doublereal *cprt) const
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 setDensity(const doublereal density_)
Set the internally stored density (kg/m^3) of the phase Note the density of a phase is an independent...
Definition: Phase.h:549
vector_fp m_h0_RT
Temporary storage for dimensionless reference state enthalpies.
vector_fp m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
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 setStateFromXML(const XML_Node &state)
Set the initial state of the phase to the conditions specified in the state XML element.