Cantera  2.4.0
Sub.cpp
Go to the documentation of this file.
1 //! @file Sub.cpp
2 /*
3  * The Substance class
4  * D. Goodwin, Caltech Nov. 1996
5  */
6 
7 // This file is part of Cantera. See License.txt in the top-level directory or
8 // at http://www.cantera.org/license.txt for license and copyright information.
9 
10 #include "cantera/tpx/Sub.h"
12 #include "cantera/base/global.h"
13 
14 using std::string;
15 using namespace Cantera;
16 
17 namespace {
18 // these correspond to ordering withing propertyFlag::type
19 std::string propertySymbols[] = {"H", "S", "U", "V", "P", "T"};
20 }
21 
22 namespace tpx
23 {
24 Substance::Substance() :
25  T(Undef),
26  Rho(Undef),
27  Tslast(Undef),
28  Rhf(Undef),
29  Rhv(Undef),
30  Pst(Undef),
31  m_energy_offset(0.0),
32  m_entropy_offset(0.0),
33  kbr(0)
34 {
35 }
36 
37 void Substance::setStdState(double h0, double s0, double t0, double p0)
38 {
39  Set(PropertyPair::TP, t0, p0);
40  double hh = h();
41  double ss = s();
42  double hoff = h0 - hh;
43  double soff = s0 - ss;
44  m_entropy_offset += soff;
45  m_energy_offset += hoff;
46 }
47 
48 double Substance::P()
49 {
50  return TwoPhase() ? Ps() : Pp();
51 }
52 
53 const double DeltaT = 0.000001;
54 
55 double Substance::cv()
56 {
57  double Tsave = T, dt = 1.e-4*T;
58  double x0 = x();
59  double T1 = std::max(Tmin(), Tsave - dt);
60  double T2 = std::min(Tmax(), Tsave + dt);
61 
62  set_T(T1);
63  double x1 = x();
64  if ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
65  // If the initial state was pure liquid or pure vapor, and the state at
66  // T-dT is not, just take a one-sided difference
67  T1 = Tsave;
68  set_T(T1);
69  }
70  double s1 = s();
71 
72  set_T(T2);
73  double x2 = x();
74  if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
75  // If the initial state was pure liquid or pure vapor, and the state at
76  // T+dT is not, just take a one-sided difference
77  T2 = Tsave;
78  set_T(T2);
79  }
80  double s2 = s();
81 
82  set_T(Tsave);
83  return T*(s2 - s1)/(T2-T1);
84 }
85 
86 double Substance::cp()
87 {
88  double Tsave = T, dt = 1.e-4*T;
89  double T1 = std::max(Tmin(), Tsave - dt);
90  double T2 = std::min(Tmax(), Tsave + dt);
91  double p0 = P();
92  double x0 = x();
93  if (TwoPhase()) {
94  // In the two-phase region, cp is infinite
95  return std::numeric_limits<double>::infinity();
96  }
97 
98  Set(PropertyPair::TP, T1, p0);
99  double x1 = x();
100  if (x1 != x0) {
101  // If the initial state was pure liquid or pure vapor, and the state at
102  // T-dT is not, just take a one-sided difference
103  T1 = Tsave;
104  Set(PropertyPair::TP, T1, p0);
105  }
106  double s1 = s();
107 
108  Set(PropertyPair::TP, T2, p0);
109  double x2 = x();
110  if (x2 != x0) {
111  // If the initial state was pure liquid or pure vapor, and the state at
112  // T+dT is not, just take a one-sided difference
113  T2 = Tsave;
114  Set(PropertyPair::TP, T2, p0);
115  }
116  double s2 = s();
117 
118  Set(PropertyPair::TP, Tsave, p0);
119  return T*(s2 - s1)/(T2-T1);
120 }
121 
122 double Substance::thermalExpansionCoeff()
123 {
124  double Tsave = T, dt = 1.e-4*T;
125  double T1 = std::max(Tmin(), Tsave - dt);
126  double T2 = std::min(Tmax(), Tsave + dt);
127  double p0 = P();
128  double x0 = x();
129 
130  if (TwoPhase()) {
131  // In the two-phase region, the thermal expansion coefficient is
132  // infinite
133  return std::numeric_limits<double>::infinity();
134  }
135 
136  Set(PropertyPair::TP, T1, p0);
137  double x1 = x();
138  if (x1 != x0) {
139  // If the initial state was pure liquid or pure vapor, and the state at
140  // T-dT is not, just take a one-sided difference
141  T1 = Tsave;
142  Set(PropertyPair::TP, T1, p0);
143  }
144  double v1 = v();
145 
146  Set(PropertyPair::TP, T2, p0);
147  double x2 = x();
148  if (x2 != x0) {
149  // If the initial state was pure liquid or pure vapor, and the state at
150  // T+dT is not, just take a one-sided difference
151  T2 = Tsave;
152  Set(PropertyPair::TP, T2, p0);
153  }
154  double v2 = v();
155 
156  Set(PropertyPair::TP, Tsave, p0);
157  return 2.0*(v2 - v1)/((v2 + v1)*(T2-T1));
158 }
159 
160 double Substance::isothermalCompressibility()
161 {
162  double Psave = P(), dp = 1.e-4*Psave;
163  double x0 = x();
164 
165  if (TwoPhase()) {
166  // In the two-phase region, the isothermal compressibility is infinite
167  return std::numeric_limits<double>::infinity();
168  }
169 
170  double v0 = v();
171  double P1 = Psave - dp;
172  double P2 = Psave + dp;
173 
174  Set(PropertyPair::TP, T, P1);
175  double x1 = x();
176  if (x1 != x0) {
177  // If the initial state was pure liquid or pure vapor, and the state at
178  // P-dP is not, just take a one-sided difference
179  P1 = Psave;
180  Set(PropertyPair::TP, T, P1);
181  }
182  double v1 = v();
183 
184  Set(PropertyPair::TP, T, P2);
185  double x2 = x();
186  if (x2 != x0) {
187  // If the initial state was pure liquid or pure vapor, and the state at
188  // P+dP is not, just take a one-sided difference
189  P2 = Psave;
190  Set(PropertyPair::TP, T, P2);
191  }
192  double v2 = v();
193 
194  Set(PropertyPair::TP, T, Psave);
195  return -(v2 - v1)/(v0*(P2-P1));
196 }
197 
198 double Substance::dPsdT()
199 {
200  double tsave = T;
201  double ps1 = Ps();
202  T += DeltaT;
203  double dpdt = (Ps() - ps1)/DeltaT;
204  T = tsave;
205  return dpdt;
206 }
207 
208 int Substance::TwoPhase()
209 {
210  if (T >= Tcrit()) {
211  return 0;
212  }
213  update_sat();
214  return ((Rho < Rhf) && (Rho > Rhv) ? 1 : 0);
215 }
216 
217 double Substance::x()
218 {
219  if (T >= Tcrit()) {
220  return (1.0/Rho < Vcrit() ? 0.0 : 1.0);
221  } else {
222  update_sat();
223  if (Rho <= Rhv) {
224  return 1.0;
225  } else if (Rho >= Rhf) {
226  return 0.0;
227  } else {
228  double vv = 1.0/Rhv;
229  double vl = 1.0/Rhf;
230  return (1.0/Rho - vl)/(vv - vl);
231  }
232  }
233 }
234 
235 double Substance::Tsat(double p)
236 {
237  if (p <= 0.0 || p > Pcrit()) {
238  throw CanteraError("Substance::Tsat", "illegal pressure value");
239  }
240  int LoopCount = 0;
241  double tol = 1.e-6*p;
242  double Tsave = T;
243  if (T < Tmin()) {
244  T = 0.5*(Tcrit() - Tmin());
245  }
246  if (T >= Tcrit()) {
247  T = 0.5*(Tcrit() - Tmin());
248  }
249  double dp = 10*tol;
250  while (fabs(dp) > tol) {
251  if (T > Tcrit()) {
252  T = Tcrit() - 0.001;
253  }
254  if (T < Tmin()) {
255  T = Tmin() + 0.001;
256  }
257  dp = p - Ps();
258  double dt = dp/dPsdT();
259  double dta = fabs(dt);
260  double dtm = 0.1*T;
261  if (dta > dtm) {
262  dt = dt*dtm/dta;
263  }
264  T += dt;
265  LoopCount++;
266  if (LoopCount > 100) {
267  T = Tsave;
268  throw CanteraError("Substance::Tsat", "No convergence");
269  }
270  }
271  double tsat = T;
272  T = Tsave;
273  return tsat;
274 }
275 
276 // property tolerances
277 static const double TolAbsH = 1.e-4; // J/kg
278 static const double TolAbsU = 1.e-4;
279 static const double TolAbsS = 1.e-7;
280 static const double TolAbsP = 0.0; // Pa, this is supposed to be zero
281 static const double TolAbsV = 1.e-8;
282 static const double TolAbsT = 1.e-3;
283 static const double TolRel = 1.e-8;
284 
285 void Substance::Set(PropertyPair::type XY, double x0, double y0)
286 {
287  double temp;
288 
289  /* if inverted (PT) switch order and change sign of XY (TP = -PT) */
290  if (XY < 0) {
291  std::swap(x0, y0);
292  XY = static_cast<PropertyPair::type>(-XY);
293  }
294 
295  switch (XY) {
296  case PropertyPair::TV:
297  set_T(x0);
298  set_v(y0);
299  break;
300  case PropertyPair::HP:
301  if (Lever(Pgiven, y0, x0, propertyFlag::H)) {
302  return;
303  }
304  set_xy(propertyFlag::H, propertyFlag::P,
305  x0, y0, TolAbsH, TolAbsP, TolRel, TolRel);
306  break;
307  case PropertyPair::SP:
308  if (Lever(Pgiven, y0, x0, propertyFlag::S)) {
309  return;
310  }
311  set_xy(propertyFlag::S, propertyFlag::P,
312  x0, y0, TolAbsS, TolAbsP, TolRel, TolRel);
313  break;
314  case PropertyPair::PV:
315  if (Lever(Pgiven, x0, y0, propertyFlag::V)) {
316  return;
317  }
318  set_xy(propertyFlag::P, propertyFlag::V,
319  x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
320  break;
321  case PropertyPair::TP:
322  if (x0 < Tcrit()) {
323  set_T(x0);
324  if (y0 < Ps()) {
325  Set(PropertyPair::TX, x0, 1.0);
326  } else {
327  Set(PropertyPair::TX, x0, 0.0);
328  }
329  } else {
330  set_T(x0);
331  }
332  set_xy(propertyFlag::T, propertyFlag::P,
333  x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);
334  break;
335  case PropertyPair::UV:
336  set_xy(propertyFlag::U, propertyFlag::V,
337  x0, y0, TolAbsU, TolAbsV, TolRel, TolRel);
338  break;
339  case PropertyPair::ST:
340  if (Lever(Tgiven, y0, x0, propertyFlag::S)) {
341  return;
342  }
343  set_xy(propertyFlag::S, propertyFlag::T,
344  x0, y0, TolAbsS, TolAbsT, TolRel, TolRel);
345  break;
346  case PropertyPair::SV:
347  set_xy(propertyFlag::S, propertyFlag::V,
348  x0, y0, TolAbsS, TolAbsV, TolRel, TolRel);
349  break;
350  case PropertyPair::UP:
351  if (Lever(Pgiven, y0, x0, propertyFlag::U)) {
352  return;
353  }
354  set_xy(propertyFlag::U, propertyFlag::P,
355  x0, y0, TolAbsU, TolAbsP, TolRel, TolRel);
356  break;
357  case PropertyPair::VH:
358  set_xy(propertyFlag::V, propertyFlag::H,
359  x0, y0, TolAbsV, TolAbsH, TolRel, TolRel);
360  break;
361  case PropertyPair::TH:
362  set_xy(propertyFlag::T, propertyFlag::H,
363  x0, y0, TolAbsT, TolAbsH, TolRel, TolRel);
364  break;
365  case PropertyPair::SH:
366  set_xy(propertyFlag::S, propertyFlag::H,
367  x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
368  break;
369  case PropertyPair::PX:
370  temp = Tsat(x0);
371  if (y0 > 1.0 || y0 < 0.0) {
372  throw CanteraError("Substance::Set",
373  "Invalid vapor fraction, {}", y0);
374  } else if (temp >= Tcrit()) {
375  throw CanteraError("Substance::Set",
376  "Can't set vapor fraction above the critical point");
377  } else {
378  set_T(temp);
379  update_sat();
380  Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
381  }
382  break;
383  case PropertyPair::TX:
384  if (y0 > 1.0 || y0 < 0.0) {
385  throw CanteraError("Substance::Set",
386  "Invalid vapor fraction, {}", y0);
387  } else if (x0 >= Tcrit()) {
388  throw CanteraError("Substance::Set",
389  "Can't set vapor fraction above the critical point");
390  } else {
391  set_T(x0);
392  update_sat();
393  Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
394  }
395  break;
396  default:
397  throw CanteraError("Substance::Set", "Invalid input.");
398  }
399 }
400 
401 //------------------ Protected and Private Functions -------------------
402 
403 void Substance::set_Rho(double r0)
404 {
405  if (r0 > 0.0) {
406  Rho = r0;
407  } else {
408  throw CanteraError("Substance::set_Rho", "Invalid density: {}", r0);
409  }
410 }
411 
412 void Substance::set_T(double t0)
413 {
414  if ((t0 >= Tmin()) && (t0 <= Tmax())) {
415  T = t0;
416  } else {
417  throw CanteraError("Substance::set_T", "illegal temperature: {}", t0);
418  }
419 }
420 
421 void Substance::set_v(double v0)
422 {
423  if (v0 > 0) {
424  Rho = 1.0/v0;
425  } else {
426  throw CanteraError("Substance::set_v",
427  "negative specific volume: {}", v0);
428  }
429 }
430 
431 double Substance::Ps()
432 {
433  if (T < Tmin() || T > Tcrit()) {
434  throw CanteraError("Substance::Ps",
435  "illegal temperature value {}", T);
436  }
437  update_sat();
438  return Pst;
439 }
440 
441 void Substance::update_sat()
442 {
443  if ((T != Tslast) && (T < Tcrit())) {
444  double Rho_save = Rho;
445  double pp = Psat();
446  double lps = log(pp);
447  // trial value = Psat from correlation
448  int i;
449  for (i = 0; i<20; i++) {
450  if (i==0) {
451  Rho = ldens(); // trial value = liquid density
452  } else {
453  Rho = Rhf;
454  }
455  set_TPp(T,pp);
456  Rhf = Rho; // sat liquid density
457 
458  double gf = hp() - T*sp();
459  if (i==0) {
460  Rho = pp*MolWt()/(GasConstant*T); // trial value = ideal gas
461  } else {
462  Rho = Rhv;
463  }
464  set_TPp(T,pp);
465 
466  Rhv = Rho; // sat vapor density
467  double gv = hp() - T*sp();
468  double dg = gv - gf;
469  if (Rhv > Rhf) {
470  std::swap(Rhv, Rhf);
471  dg = - dg;
472  }
473 
474  if (fabs(dg) < 0.001 && Rhf > Rhv) {
475  break;
476  }
477  double dp = dg/(1.0/Rhv - 1.0/Rhf);
478  double psold = pp;
479  if (fabs(dp) > pp) {
480  lps -= dg/(pp*(1.0/Rhv - 1.0/Rhf));
481  pp = exp(lps);
482  } else {
483  pp -= dp;
484  lps = log(pp);
485  }
486  if (pp > Pcrit()) {
487  pp = psold + 0.5*(Pcrit() - psold);
488  lps = log(pp);
489  } else if (pp < 0.0) {
490  pp = psold/2.0;
491  lps = log(pp);
492  }
493  }
494  if (Rhf <= Rhv) {
495  throw CanteraError("Substance::update_sat",
496  "wrong root found for sat. liquid or vapor at P = {}", pp);
497  }
498 
499  if (i >= 20) {
500  throw CanteraError("substance::update_sat","no convergence");
501  } else {
502  Pst = pp;
503  Tslast = T;
504  }
505  Rho = Rho_save;
506  }
507 }
508 
509 double Substance::vprop(propertyFlag::type ijob)
510 {
511  switch (ijob) {
512  case propertyFlag::H:
513  return hp();
514  case propertyFlag::S:
515  return sp();
516  case propertyFlag::U:
517  return up();
518  case propertyFlag::V:
519  return vp();
520  case propertyFlag::P:
521  return Pp();
522  default:
523  throw CanteraError("Substance::vprop", "invalid job index");
524  }
525 }
526 
527 int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc)
528 {
529  double psat;
530  double Tsave = T;
531  double Rhosave = Rho;
532  if (itp == Tgiven) {
533  if (sat >= Tcrit()) {
534  return 0;
535  }
536  T = sat;
537  psat = Ps();
538  } else if (itp == Pgiven) {
539  if (sat >= Pcrit()) {
540  return 0;
541  }
542  psat = sat;
543  try {
544  T = Tsat(psat);
545  } catch (CanteraError&) {
546  // Failure to converge here is not an error
547  T = Tsave;
548  Rho = Rhosave;
549  return 0;
550  }
551  } else {
552  throw CanteraError("Substance::Lever","general error");
553  }
554  Set(PropertyPair::TX, T, 1.0);
555  double Valg = vprop(ifunc);
556  Set(PropertyPair::TX, T, 0.0);
557  double Valf = vprop(ifunc);
558  if (val >= Valf && val <= Valg) {
559  double xx = (val - Valf)/(Valg - Valf);
560  double vv = (1.0 - xx)/Rhf + xx/Rhv;
561  set_v(vv);
562  return 1;
563  } else {
564  T = Tsave;
565  Rho = Rhosave;
566  return 0;
567  }
568 }
569 
570 void Substance::set_xy(propertyFlag::type ifx, propertyFlag::type ify,
571  double X, double Y,
572  double atx, double aty,
573  double rtx, double rty)
574 {
575  double v_here, t_here;
576  double dvs1 = 2.0*Vcrit();
577  double dvs2 = 0.7*Vcrit();
578  int LoopCount = 0;
579 
580  double v_save = 1.0/Rho;
581  double t_save = T;
582 
583  if ((T == Undef) && (Rho == Undef)) {
584  // new object, try to pick a "reasonable" starting point
585  Set(PropertyPair::TV,Tcrit()*1.1,Vcrit()*1.1);
586  t_here = T;
587  v_here = 1.0/Rho;
588  } else if (Rho == Undef) {
589  // new object, try to pick a "reasonable" starting point
590  Set(PropertyPair::TV,T,Vcrit()*1.1);
591  t_here = T;
592  v_here = 1.0/Rho;
593  } else {
594  v_here = v_save;
595  t_here = t_save;
596  }
597 
598  double Xa = fabs(X);
599  double Ya = fabs(Y);
600  while (true) {
601  double x_here = prop(ifx);
602  double y_here = prop(ify);
603  double err_x = fabs(X - x_here);
604  double err_y = fabs(Y - y_here);
605 
606  if ((err_x < atx + rtx*Xa) && (err_y < aty + rty*Ya)) {
607  break;
608  }
609 
610  /* perturb t */
611  double dt = 0.001*t_here;
612  if (t_here + dt > Tmax()) {
613  dt *= -1.0;
614  }
615 
616  /* perturb v */
617  double dv = 0.001*v_here;
618  if (v_here <= Vcrit()) {
619  dv *= -1.0;
620  }
621 
622  /* derivatives with respect to T */
623  Set(PropertyPair::TV, t_here + dt, v_here);
624  double dxdt = (prop(ifx) - x_here)/dt;
625  double dydt = (prop(ify) - y_here)/dt;
626 
627  /* derivatives with respect to v */
628  Set(PropertyPair::TV, t_here, v_here + dv);
629  double dxdv = (prop(ifx) - x_here)/dv;
630  double dydv = (prop(ify) - y_here)/dv;
631 
632  double det = dxdt*dydv - dydt*dxdv;
633  dt = ((X - x_here)*dydv - (Y - y_here)*dxdv)/det;
634  dv = ((Y - y_here)*dxdt - (X - x_here)*dydt)/det;
635 
636  double dvm = 0.2*v_here;
637  if (v_here < dvs1) {
638  dvm *= 0.5;
639  }
640  if (v_here < dvs2) {
641  dvm *= 0.5;
642  }
643  double dtm = 0.1*t_here;
644  double dva = fabs(dv);
645  double dta = fabs(dt);
646  if (dva > dvm) {
647  dv *= dvm/dva;
648  }
649  if (dta > dtm) {
650  dt *= dtm/dta;
651  }
652  v_here += dv;
653  t_here += dt;
654  t_here = clip(t_here, Tmin(), Tmax());
655  if (v_here <= 0.0) {
656  v_here = 0.0001;
657  }
658  Set(PropertyPair::TV, t_here, v_here);
659  LoopCount++;
660  if (LoopCount > 200) {
661  std::string msg = fmt::format("No convergence. {} = {}, {} = {}",
662  propertySymbols[ifx], X, propertySymbols[ify], Y);
663  if (t_here == Tmin()) {
664  msg += fmt::format("\nAt temperature limit (Tmin = {})", Tmin());
665  } else if (t_here == Tmax()) {
666  msg += fmt::format("\nAt temperature limit (Tmax = {})", Tmax());
667  }
668  throw CanteraError("Substance::set_xy", msg);
669  }
670  }
671 }
672 
673 double Substance::prop(propertyFlag::type ijob)
674 {
675  if (ijob == propertyFlag::P) {
676  return P();
677  }
678  if (ijob == propertyFlag::T) {
679  return T;
680  }
681  double xx = x();
682  if ((xx > 0.0) && (xx < 1.0)) {
683  double Rho_save = Rho;
684  Rho = Rhv;
685  double vp = vprop(ijob);
686  Rho = Rhf;
687  double lp = vprop(ijob);
688  double pp = (1.0 - xx)*lp + xx*vp;
689  Rho = Rho_save;
690  return pp;
691  } else {
692  return vprop(ijob);
693  }
694 }
695 
696 static const double ErrP = 1.e-7;
697 static const double Big = 1.e30;
698 
699 void Substance::BracketSlope(double Pressure)
700 {
701  if (kbr == 0) {
702  dv = (v_here < Vcrit() ? -0.05*v_here : 0.2*v_here);
703  if (Vmin > 0.0) {
704  dv = 0.2*v_here;
705  }
706  if (Vmax < Big) {
707  dv = -0.05*v_here;
708  }
709  } else {
710  double dpdv = (Pmax - Pmin)/(Vmax - Vmin);
711  v_here = Vmax;
712  P_here = Pmax;
713  dv = dvbf*(Pressure - P_here)/dpdv;
714  dvbf = 0.5*dvbf;
715  }
716 }
717 
718 void Substance::set_TPp(double Temp, double Pressure)
719 {
720  kbr = 0;
721  dvbf = 1.0;
722  Vmin = 0.0;
723  Vmax = Big;
724  Pmin = Big;
725  Pmax = 0.0;
726  double dvs1 = 2.0*Vcrit();
727  double dvs2 = 0.7*Vcrit();
728  int LoopCount = 0;
729 
730  double v_save = 1.0/Rho;
731  T = Temp;
732  v_here = vp();
733 
734  // loop
735  while (P_here = Pp(),
736  fabs(Pressure - P_here) >= ErrP* Pressure || LoopCount == 0) {
737  if (P_here < 0.0) {
738  BracketSlope(Pressure);
739  } else {
740  dv = 0.001*v_here;
741  if (v_here <= Vcrit()) {
742  dv *= -1.0;
743  }
744  Set(PropertyPair::TV, Temp, v_here+dv);
745  double dpdv = (Pp() - P_here)/dv;
746  if (dpdv > 0.0) {
747  BracketSlope(Pressure);
748  } else {
749  if ((P_here > Pressure) && (v_here > Vmin)) {
750  Vmin = v_here;
751  } else if ((P_here < Pressure) && (v_here < Vmax)) {
752  Vmax = v_here;
753  }
754  if (v_here == Vmin) {
755  Pmin = P_here;
756  }
757  if (v_here == Vmax) {
758  Pmax = P_here;
759  }
760  if (Vmin >= Vmax) {
761  throw CanteraError("Substance::set_TPp","Vmin >= Vmax");
762  } else if ((Vmin > 0.0) && (Vmax < Big)) {
763  kbr = 1;
764  }
765  dvbf = 1.0;
766  if (dpdv == 0.0) {
767  dvbf = 0.5;
768  BracketSlope(Pressure);
769  } else {
770  dv = (Pressure - P_here)/dpdv;
771  }
772  }
773  }
774  double dvm = 0.2*v_here;
775  if (v_here < dvs1) {
776  dvm *= 0.5;
777  }
778  if (v_here < dvs2) {
779  dvm *= 0.5;
780  }
781  if (kbr != 0) {
782  double vt = v_here + dv;
783  if ((vt < Vmin) || (vt > Vmax)) {
784  dv = Vmin + (Pressure - Pmin)*(Vmax - Vmin)/(Pmax - Pmin) - v_here;
785  }
786  }
787  double dva = fabs(dv);
788  if (dva > dvm) {
789  dv *= dvm/dva;
790  }
791  v_here += dv;
792  if (dv == 0.0) {
793  throw CanteraError("Substance::set_TPp","dv = 0 and no convergence");
794  }
795  Set(PropertyPair::TV, Temp, v_here);
796  LoopCount++;
797  if (LoopCount > 100) {
798  Set(PropertyPair::TV, Temp, v_save);
799  throw CanteraError("Substance::set_TPp",
800  "no convergence for P* = {}, V* = {}",
801  Pressure/Pcrit(), v_save/Vcrit());
802  }
803  }
804  Set(PropertyPair::TV, Temp,v_here);
805 }
806 }
T clip(const T &value, const T &lower, const T &upper)
Clip value such that lower <= value <= upper.
Definition: global.h:269
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
const doublereal Undef
Fairly random number to be used to initialize variables against to see if they are subsequently defin...
Definition: ct_defs.h:134
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Contains declarations for string manipulation functions within Cantera.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8