18string propertySymbols[] = {
"H",
"S",
"U",
"V",
"P",
"T"};
24void Substance::setStdState(
double h0,
double s0,
double t0,
double p0)
26 Set(PropertyPair::TP, t0, p0);
29 double hoff = h0 - hh;
30 double soff = s0 - ss;
31 m_entropy_offset += soff;
32 m_energy_offset += hoff;
40const double DeltaT = 0.000001;
48 return std::numeric_limits<double>::quiet_NaN();
51 double Tsave = T, dt = 1.e-4 * T;
53 double T1 = std::max(
Tmin(), Tsave - dt);
54 double T2 = std::min(
Tmax(), Tsave + dt);
58 if ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
68 if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
77 return T*(s2 - s1)/(T2-T1);
84 return std::numeric_limits<double>::infinity();
87 double Tsave = T, dt = 1.e-4 * T;
89 double T1, T2, s1, s2;
94 T1 = std::max(
Tmin(), Tsave - dt);
95 Set(PropertyPair::TP, T1, p0);
98 T2 = std::min(
Tsat(p0), Tsave + dt);
103 if (T2 < Tsave + dt) {
104 Set(PropertyPair::TX, T2, 0.);
106 Set(PropertyPair::TP, T2, p0);
112 T1 = std::max(
Tsat(p0), Tsave - dt);
117 if (T1 > Tsave - dt) {
118 Set(PropertyPair::TX, T1, 1.);
120 Set(PropertyPair::TP, T1, p0);
123 T2 = std::min(
Tmax(), Tsave + dt);
124 Set(PropertyPair::TP, T2, p0);
128 Set(PropertyPair::TV, Tsave, 1.0 / RhoSave);
129 return T * (s2 - s1) / (T2 - T1);
132double Substance::thermalExpansionCoeff()
137 return std::numeric_limits<double>::infinity();
140 double Tsave = T, dt = 1.e-4 * T;
141 double RhoSave = Rho;
142 double T1, T2, v1, v2;
147 T1 = std::max(
Tmin(), Tsave - dt);
148 Set(PropertyPair::TP, T1, p0);
151 T2 = std::min(
Tsat(p0), Tsave + dt);
156 if (T2 < Tsave + dt) {
157 Set(PropertyPair::TX, T2, 0.);
159 Set(PropertyPair::TP, T2, p0);
165 T1 = std::max(
Tsat(p0), Tsave - dt);
170 if (T1 > Tsave - dt) {
171 Set(PropertyPair::TX, T1, 1.);
173 Set(PropertyPair::TP, T1, p0);
176 T2 = std::min(
Tmax(), Tsave + dt);
177 Set(PropertyPair::TP, T2, p0);
181 Set(PropertyPair::TV, Tsave, 1.0 / RhoSave);
182 return 2.0 * (v2 - v1) / ((v2 + v1) * (T2 - T1));
185double Substance::isothermalCompressibility()
189 return std::numeric_limits<double>::infinity();
192 double Psave =
P(), dp = 1.e-4 * Psave;
193 double RhoSave = Rho;
194 double P1, P2, v1, v2;
201 P1 = std::max(Ps(), P1);
203 if (P1 > Psave - dp) {
204 Set(PropertyPair::PX, P1, 0.);
206 Set(PropertyPair::TP, T, P1);
210 Set(PropertyPair::TP, T, P2);
214 P1 = std::max(1.e-7, Psave - dp);
215 Set(PropertyPair::TP, T, P1);
219 P2 = std::min(Ps(), P2);
221 if (P2 < Psave + dp) {
222 Set(PropertyPair::PX, P2, 1.);
224 Set(PropertyPair::TP, T, P2);
229 Set(PropertyPair::TV, T, 1.0 / RhoSave);
230 return -(v2 - v1) / (v0 * (P2 - P1));
238 if (T + DeltaT <
Tcrit()) {
240 dpdt = (Ps() - ps1)/DeltaT;
244 dpdt = (ps1 - Ps())/DeltaT;
257 return ((Rho > Rhv) && (Rho < Rhf) ? 1 : 0);
259 return ((Rho >= Rhv) && (Rho <= Rhf) ? 1 : 0);
265 return (1.0/Rho <
Vcrit() ? 0.0 : 1.0);
270 }
else if (Rho >= Rhf) {
275 return (1.0/Rho - vl)/(vv - vl);
283 if (p <= 0. || p >
Pcrit()) {
291 "Illegal pressure value: {}", p);
297 double tol = 1.e-6*p;
302 while (fabs(dp) > tol) {
303 T = std::min(T,
Tcrit());
304 T = std::max(T,
Tmin());
306 double dt = dp/
dPsdT();
307 double dta = fabs(dt);
314 if (LoopCount > 100) {
316 throw CanteraError(
"Substance::Tsat",
"No convergence: p = {}", p);
325static const double TolAbsH = 1.e-4;
326static const double TolAbsU = 1.e-4;
327static const double TolAbsS = 1.e-7;
328static const double TolAbsP = 0.0;
329static const double TolAbsV = 1.e-8;
330static const double TolAbsT = 1.e-3;
331static const double TolRel = 1.e-8;
340 XY =
static_cast<PropertyPair::type
>(-XY);
344 case PropertyPair::TV:
348 case PropertyPair::HP:
349 if (
Lever(Pgiven, y0, x0, propertyFlag::H)) {
352 set_xy(propertyFlag::H, propertyFlag::P,
353 x0, y0, TolAbsH, TolAbsP, TolRel, TolRel);
355 case PropertyPair::SP:
356 if (
Lever(Pgiven, y0, x0, propertyFlag::S)) {
359 set_xy(propertyFlag::S, propertyFlag::P,
360 x0, y0, TolAbsS, TolAbsP, TolRel, TolRel);
362 case PropertyPair::PV:
363 if (
Lever(Pgiven, x0, y0, propertyFlag::V)) {
366 set_xy(propertyFlag::P, propertyFlag::V,
367 x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
369 case PropertyPair::TP:
372 if (fabs(y0 - Ps()) / y0 < TolRel) {
374 "Saturated mixture detected: use vapor "
375 "fraction to specify state instead");
376 }
else if (y0 < Ps()) {
377 Set(PropertyPair::TX, x0, 1.0);
379 Set(PropertyPair::TX, x0, 0.0);
382 set_xy(propertyFlag::T, propertyFlag::P,
383 x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);
385 case PropertyPair::UV:
386 set_xy(propertyFlag::U, propertyFlag::V,
387 x0, y0, TolAbsU, TolAbsV, TolRel, TolRel);
389 case PropertyPair::ST:
390 if (
Lever(Tgiven, y0, x0, propertyFlag::S)) {
393 set_xy(propertyFlag::S, propertyFlag::T,
394 x0, y0, TolAbsS, TolAbsT, TolRel, TolRel);
396 case PropertyPair::SV:
397 set_xy(propertyFlag::S, propertyFlag::V,
398 x0, y0, TolAbsS, TolAbsV, TolRel, TolRel);
400 case PropertyPair::UP:
401 if (
Lever(Pgiven, y0, x0, propertyFlag::U)) {
404 set_xy(propertyFlag::U, propertyFlag::P,
405 x0, y0, TolAbsU, TolAbsP, TolRel, TolRel);
407 case PropertyPair::VH:
408 set_xy(propertyFlag::V, propertyFlag::H,
409 x0, y0, TolAbsV, TolAbsH, TolRel, TolRel);
411 case PropertyPair::TH:
412 set_xy(propertyFlag::T, propertyFlag::H,
413 x0, y0, TolAbsT, TolAbsH, TolRel, TolRel);
415 case PropertyPair::SH:
416 set_xy(propertyFlag::S, propertyFlag::H,
417 x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
419 case PropertyPair::PX:
422 if (y0 > 1.0 || y0 < 0.0) {
424 "Invalid vapor fraction, {}", y0);
426 if (x0 < Ps() || x0 >
Pcrit()) {
428 "Illegal pressure value: {} (supercritical "
429 "or below triple point)", x0);
434 Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
436 case PropertyPair::TX:
437 if (y0 > 1.0 || y0 < 0.0) {
439 "Invalid vapor fraction, {}", y0);
443 "Illegal temperature value: {} "
444 "(supercritical or below triple point)", x0);
448 Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
457void Substance::set_Rho(
double r0)
462 throw CanteraError(
"Substance::set_Rho",
"Invalid density: {}", r0);
466void Substance::set_T(
double t0)
468 if ((t0 >=
Tmin()) && (t0 <=
Tmax())) {
471 throw CanteraError(
"Substance::set_T",
"Illegal temperature: {}", t0);
475void Substance::set_v(
double v0)
481 "Negative specific volume: {}", v0);
485double Substance::Ps()
489 "Illegal temperature value: {}", T);
498 double Rho_save = Rho;
500 double lps = log(pp);
503 for (i = 0; i<20; i++) {
512 double gf =
hp() - T*
sp();
521 double gv =
hp() - T*
sp();
528 if (fabs(dg) < 0.001) {
531 double dp = dg/(1.0/Rhv - 1.0/Rhf);
534 lps -= dg/(pp*(1.0/Rhv - 1.0/Rhf));
541 pp = psold + 0.5*(
Pcrit() - psold);
543 }
else if (pp < 0.0) {
550 throw CanteraError(
"Substance::update_sat",
"No convergence");
559double Substance::vprop(propertyFlag::type ijob)
562 case propertyFlag::H:
564 case propertyFlag::S:
566 case propertyFlag::U:
568 case propertyFlag::V:
570 case propertyFlag::P:
573 throw CanteraError(
"Substance::vprop",
"Invalid job index");
581 double Rhosave = Rho;
583 if (sat >=
Tcrit()) {
588 }
else if (itp == Pgiven) {
589 if (sat >=
Pcrit()) {
602 throw CanteraError(
"Substance::Lever",
"General error");
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;
620void Substance::set_xy(propertyFlag::type ifx, propertyFlag::type ify,
622 double atx,
double aty,
623 double rtx,
double rty)
625 double v_here, t_here;
626 double dvs1 = 2.0*
Vcrit();
627 double dvs2 = 0.7*
Vcrit();
630 double v_save = 1.0/Rho;
633 if ((T == Undef) && (Rho == Undef)) {
638 }
else if (Rho == Undef) {
640 Set(PropertyPair::TV,T,
Vcrit()*1.1);
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);
656 if ((err_x < atx + rtx*Xa) && (err_y < aty + rty*Ya)) {
661 double dt = 0.001*t_here;
662 if (t_here + dt >
Tmax()) {
667 double dv = 0.001*v_here;
668 if (v_here <=
Vcrit()) {
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;
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;
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;
686 double dvm = 0.2*v_here;
693 double dtm = 0.1*t_here;
694 double dva = fabs(dv);
695 double dta = fabs(dt);
708 Set(PropertyPair::TV, t_here, v_here);
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());
723double Substance::prop(propertyFlag::type ijob)
725 if (ijob == propertyFlag::P) {
728 if (ijob == propertyFlag::T) {
732 if ((xx > 0.0) && (xx < 1.0)) {
733 double Rho_save = Rho;
735 double vp = vprop(ijob);
737 double lp = vprop(ijob);
738 double pp = (1.0 - xx)*lp + xx*vp;
746static const double ErrP = 1.e-7;
747static const double Big = 1.e30;
749void Substance::BracketSlope(
double Pressure)
752 dv = (v_here <
Vcrit() ? -0.05*v_here : 0.2*v_here);
760 double dpdv = (Pmax - Pmin)/(Vmax - Vmin);
763 dv = dvbf*(Pressure - P_here)/dpdv;
776 double dvs1 = 2.0*
Vcrit();
777 double dvs2 = 0.7*
Vcrit();
780 double v_save = 1.0/Rho;
785 while (P_here = Pp(),
786 fabs(Pressure - P_here) >= ErrP* Pressure || LoopCount == 0) {
788 BracketSlope(Pressure);
791 if (v_here <=
Vcrit()) {
794 Set(PropertyPair::TV,
Temp, v_here+dv);
795 double dpdv = (Pp() - P_here)/dv;
797 BracketSlope(Pressure);
799 if ((P_here > Pressure) && (v_here > Vmin)) {
801 }
else if ((P_here < Pressure) && (v_here < Vmax)) {
804 if (v_here == Vmin) {
807 if (v_here == Vmax) {
811 throw CanteraError(
"Substance::set_TPp",
"Vmin >= Vmax");
812 }
else if ((Vmin > 0.0) && (Vmax < Big)) {
818 BracketSlope(Pressure);
820 dv = (Pressure - P_here)/dpdv;
824 double dvm = 0.2*v_here;
832 double vt = v_here + dv;
833 if ((vt < Vmin) || (vt > Vmax)) {
834 dv = Vmin + (Pressure - Pmin)*(Vmax - Vmin)/(Pmax - Pmin) - v_here;
837 double dva = fabs(dv);
843 throw CanteraError(
"Substance::set_TPp",
"dv = 0 and no convergence");
845 Set(PropertyPair::TV,
Temp, v_here);
847 if (LoopCount > 200) {
848 Set(PropertyPair::TV,
Temp, v_save);
850 "No convergence for P = {}, T = {}\n"
851 "(P* = {}, V* = {})", Pressure,
Temp,
855 Set(PropertyPair::TV,
Temp,v_here);
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.
int TwoPhase(bool strict=false)
Returns 1 if the current state is a liquid/vapor mixture, 0 otherwise.
virtual double cp()
Specific heat at constant pressure [J/kg/K].
double x()
Vapor mass fraction.
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].
virtual double cv()
Specific heat at constant volume [J/kg/K].
double Temp()
Temperature [K].
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.
virtual double sp()=0
Entropy of a single-phase state.
virtual double Psat()=0
Saturation pressure, Pa.
double v()
Specific volume [m^3/kg].
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...
virtual double Pcrit()=0
Critical pressure [Pa].
virtual double dPsdT()
The derivative of the saturation pressure with respect to temperature.
double Tsat(double p)
Saturation temperature at pressure p.
void update_sat()
Update saturated liquid and vapor densities and saturation pressure.
void set_TPp(double t0, double p0)
set T and P
double s()
Entropy [J/kg/K].
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.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
Contains declarations for string manipulation functions within Cantera.