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