Use actual enums for the pure fluid constants
This commit is contained in:
Ray Speth
parent
b7ee30aa04
commit
3e71897a9c
4 changed files with 180 additions and 216 deletions
|
|
@ -15,11 +15,7 @@
|
|||
|
||||
#include "ThermoPhase.h"
|
||||
#include "mix_defs.h"
|
||||
|
||||
namespace tpx
|
||||
{
|
||||
class Substance;
|
||||
}
|
||||
#include "cantera/tpx/Sub.h"
|
||||
|
||||
namespace Cantera
|
||||
{
|
||||
|
|
@ -548,7 +544,7 @@ protected:
|
|||
* @param x Value of the first component
|
||||
* @param y Value of the second component
|
||||
*/
|
||||
void Set(int n, double x, double y) const;
|
||||
void Set(tpx::PropertyPair::type n, double x, double y) const;
|
||||
|
||||
//! Sets the state using a TPX::TV call
|
||||
void setTPXState() const;
|
||||
|
|
|
|||
|
|
@ -33,40 +33,22 @@ public:
|
|||
TPX_Error(std::string p, std::string e) : CanteraError(p, e) { }
|
||||
};
|
||||
|
||||
|
||||
const double OneAtm = 1.01325e5;
|
||||
const double Liquid = 0.0;
|
||||
const double Vapor = 1.0;
|
||||
|
||||
const int TV = 12, HP = 34, SP = 54, PV = 42, TP = 14, UV = 62,
|
||||
ST = 51, SV = 52, UP = 64, VH = 23, TH = 13, SH = 53,
|
||||
PX = 47, TX = 17;
|
||||
|
||||
const int VT = -12, PH = -34, PS = -54, VP = -42, PT = -14, VU = -62,
|
||||
TS = -51, VS = -52, PU = -64, HV = -23, HT = -13, HS = -53,
|
||||
XP = -47, XT = -17;
|
||||
namespace PropertyPair
|
||||
{
|
||||
enum type {
|
||||
TV = 12, HP = 34, SP = 54, PV = 42, TP = 14, UV = 62, ST = 51,
|
||||
SV = 52, UP = 64, VH = 23, TH = 13, SH = 53, PX = 47, TX = 17,
|
||||
VT = -12, PH = -34, PS = -54, VP = -42, PT = -14, VU = -62, TS = -51,
|
||||
VS = -52, PU = -64, HV = -23, HT = -13, HS = -53, XP = -47, XT = -17
|
||||
};
|
||||
}
|
||||
|
||||
const int Pgiven = 0, Tgiven = 1;
|
||||
|
||||
const int EvalH = 3;
|
||||
const int EvalS = 5;
|
||||
const int EvalU = 6;
|
||||
const int EvalV = 2;
|
||||
const int EvalP = 4;
|
||||
const int EvalT = 1;
|
||||
const int EvalX = 7;
|
||||
const int EvalF = 8;
|
||||
const int EvalG = 9;
|
||||
const int EvalC = 10;
|
||||
const int EvalM = 11;
|
||||
const int EvalN = 12;
|
||||
const int EvalMW = 13;
|
||||
const int EvalEA = 14;
|
||||
const int EvalCdot = 15;
|
||||
const int EvalDdot = 16;
|
||||
const int EvalWdot = 17;
|
||||
const int EvalTchem = 18;
|
||||
const int EvalRgas = 19;
|
||||
namespace propertyFlag
|
||||
{
|
||||
enum type { H, S, U, V, P, T };
|
||||
}
|
||||
|
||||
const double Undef = 999.1234;
|
||||
|
||||
|
|
@ -79,7 +61,7 @@ public:
|
|||
|
||||
void setStdState(double h0 = 0.0, double s0 = 0.0,
|
||||
double t0 = 298.15, double p0 = 1.01325e5) {
|
||||
Set(TP, t0, p0);
|
||||
Set(PropertyPair::TP, t0, p0);
|
||||
double hh = h();
|
||||
double ss = s();
|
||||
double hoff = h0 - hh;
|
||||
|
|
@ -106,16 +88,16 @@ public:
|
|||
return T; // temperature, K
|
||||
}
|
||||
double v() { // specific vol, m^3/kg
|
||||
return prop(EvalV);
|
||||
return prop(propertyFlag::V);
|
||||
}
|
||||
double u() { // int. energy, J/kg
|
||||
return prop(EvalU);
|
||||
return prop(propertyFlag::U);
|
||||
}
|
||||
double h() { // enthalpy, J/kg
|
||||
return prop(EvalH);
|
||||
return prop(propertyFlag::H);
|
||||
}
|
||||
double s() { // entropy, J/kg/K
|
||||
return prop(EvalS);
|
||||
return prop(propertyFlag::S);
|
||||
}
|
||||
double f() { // Helmholtz function, J/kg
|
||||
return u() - T*s();
|
||||
|
|
@ -137,36 +119,35 @@ public:
|
|||
virtual double cp() {
|
||||
double Tsave = T, dt = 1.e-4*T;
|
||||
double p0 = P();
|
||||
Set(TP, Tsave - dt, p0);
|
||||
Set(PropertyPair::TP, Tsave - dt, p0);
|
||||
double s1 = s();
|
||||
Set(TP, Tsave + dt, p0);
|
||||
Set(PropertyPair::TP, Tsave + dt, p0);
|
||||
double s2 = s();
|
||||
Set(TP, Tsave, p0);
|
||||
Set(PropertyPair::TP, Tsave, p0);
|
||||
return T*(s2 - s1)/(2.0*dt);
|
||||
}
|
||||
|
||||
virtual double thermalExpansionCoeff() {
|
||||
double Tsave = T, dt = 1.e-4*T;
|
||||
double p0 = P();
|
||||
Set(TP, Tsave - dt, p0);
|
||||
Set(PropertyPair::TP, Tsave - dt, p0);
|
||||
double v1 = v();
|
||||
Set(TP, Tsave + dt, p0);
|
||||
Set(PropertyPair::TP, Tsave + dt, p0);
|
||||
double v2 = v();
|
||||
Set(TP, Tsave, p0);
|
||||
Set(PropertyPair::TP, Tsave, p0);
|
||||
return (v2 - v1)/((v2 + v1)*dt);
|
||||
}
|
||||
|
||||
virtual double isothermalCompressibility() {
|
||||
double Psave = P(), dp = 1.e-4*Psave;
|
||||
Set(TP, T, Psave - dp);
|
||||
Set(PropertyPair::TP, T, Psave - dp);
|
||||
double v1 = v();
|
||||
Set(TP, T, Psave + dp);
|
||||
Set(PropertyPair::TP, T, Psave + dp);
|
||||
double v2 = v();
|
||||
Set(TP, T, Psave);
|
||||
Set(PropertyPair::TP, T, Psave);
|
||||
return -(v2 - v1)/((v2 + v1)*dp);
|
||||
}
|
||||
|
||||
|
||||
// saturation properties
|
||||
|
||||
double Ps();
|
||||
|
|
@ -182,14 +163,12 @@ public:
|
|||
return hp() - T*sp();
|
||||
}
|
||||
|
||||
double prop(int ijob);
|
||||
double prop(propertyFlag::type ijob);
|
||||
void set_TPp(double t0, double p0); // set T and P
|
||||
|
||||
|
||||
// functions to set or change state:
|
||||
|
||||
void Set(int XY, double x0, double y0);
|
||||
void Set_meta(double phase, double pp);
|
||||
void Set(PropertyPair::type XY, double x0, double y0);
|
||||
|
||||
protected:
|
||||
|
||||
|
|
@ -214,7 +193,7 @@ protected:
|
|||
double vp() {
|
||||
return 1.0/Rho;
|
||||
}
|
||||
int Lever(int itp, double sat, double val, int ifunc);
|
||||
int Lever(int itp, double sat, double val, propertyFlag::type ifunc);
|
||||
void update_sat();
|
||||
|
||||
private:
|
||||
|
|
@ -222,9 +201,9 @@ private:
|
|||
void set_T(double t0);
|
||||
void set_v(double v0);
|
||||
void BracketSlope(double p);
|
||||
double lprop(int ijob);
|
||||
double vprop(int ijob);
|
||||
void set_xy(int if1, int if2, double X, double Y,
|
||||
double vprop(propertyFlag::type ijob);
|
||||
void set_xy(propertyFlag::type if1, propertyFlag::type if2,
|
||||
double X, double Y,
|
||||
double atx, double aty, double rtx, double rty);
|
||||
|
||||
int kbr;
|
||||
|
|
|
|||
|
|
@ -102,13 +102,13 @@ initThermo()
|
|||
double cp0_R, h0_RT, s0_R, T0, p;
|
||||
T0 = 298.15;
|
||||
if (T0 < m_sub->Tcrit()) {
|
||||
m_sub->Set(tpx::TX, T0, 1.0);
|
||||
m_sub->Set(tpx::PropertyPair::TX, T0, 1.0);
|
||||
p = 0.01*m_sub->P();
|
||||
} else {
|
||||
p = 0.001*m_sub->Pcrit();
|
||||
}
|
||||
p = 0.001 * p;
|
||||
m_sub->Set(tpx::TP, T0, p);
|
||||
m_sub->Set(tpx::PropertyPair::TP, T0, p);
|
||||
|
||||
m_spthermo->update_one(0, T0, &cp0_R, &h0_RT, &s0_R);
|
||||
double s_R = s0_R - log(p/refPressure());
|
||||
|
|
@ -182,18 +182,18 @@ pressure() const
|
|||
void PureFluidPhase::
|
||||
setPressure(doublereal p)
|
||||
{
|
||||
Set(tpx::TP, temperature(), p);
|
||||
Set(tpx::PropertyPair::TP, temperature(), p);
|
||||
setDensity(1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::Set(int n, double x, double y) const
|
||||
void PureFluidPhase::Set(tpx::PropertyPair::type n, double x, double y) const
|
||||
{
|
||||
m_sub->Set(n, x, y);
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::setTPXState() const
|
||||
{
|
||||
Set(tpx::TV, temperature(), 1.0/density());
|
||||
Set(tpx::PropertyPair::TV, temperature(), 1.0/density());
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
||||
|
|
@ -356,9 +356,9 @@ void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const
|
|||
double t = temperature();
|
||||
//double pref = m_spthermo->refPressure();
|
||||
double plow = 1.0E-8;
|
||||
Set(tpx::TP, t, plow);
|
||||
Set(tpx::PropertyPair::TP, t, plow);
|
||||
getEnthalpy_RT(hrt);
|
||||
Set(tpx::TP, t, psave);
|
||||
Set(tpx::PropertyPair::TP, t, psave);
|
||||
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
|
@ -374,10 +374,10 @@ void PureFluidPhase::getGibbs_RT_ref(doublereal* grt) const
|
|||
double t = temperature();
|
||||
double pref = m_spthermo->refPressure();
|
||||
double plow = 1.0E-8;
|
||||
Set(tpx::TP, t, plow);
|
||||
Set(tpx::PropertyPair::TP, t, plow);
|
||||
getGibbs_RT(grt);
|
||||
grt[0] += log(pref/plow);
|
||||
Set(tpx::TP, t, psave);
|
||||
Set(tpx::PropertyPair::TP, t, psave);
|
||||
}
|
||||
//====================================================================================================================
|
||||
// Returns the vector of the gibbs function of the reference state at the current temperature
|
||||
|
|
@ -406,10 +406,10 @@ void PureFluidPhase::getEntropy_R_ref(doublereal* er) const
|
|||
double t = temperature();
|
||||
double pref = m_spthermo->refPressure();
|
||||
double plow = 1.0E-8;
|
||||
Set(tpx::TP, t, plow);
|
||||
Set(tpx::PropertyPair::TP, t, plow);
|
||||
getEntropy_R(er);
|
||||
er[0] -= log(pref/plow);
|
||||
Set(tpx::TP, t, psave);
|
||||
Set(tpx::PropertyPair::TP, t, psave);
|
||||
}
|
||||
//====================================================================================================================
|
||||
// critical temperature
|
||||
|
|
@ -441,28 +441,28 @@ doublereal PureFluidPhase::satTemperature(doublereal p) const
|
|||
void PureFluidPhase::setState_HP(doublereal h, doublereal p,
|
||||
doublereal tol)
|
||||
{
|
||||
Set(tpx::HP, h, p);
|
||||
Set(tpx::PropertyPair::HP, h, p);
|
||||
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::setState_UV(doublereal u, doublereal v,
|
||||
doublereal tol)
|
||||
{
|
||||
Set(tpx::UV, u, v);
|
||||
Set(tpx::PropertyPair::UV, u, v);
|
||||
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::setState_SV(doublereal s, doublereal v,
|
||||
doublereal tol)
|
||||
{
|
||||
Set(tpx::SV, s, v);
|
||||
Set(tpx::PropertyPair::SV, s, v);
|
||||
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::setState_SP(doublereal s, doublereal p,
|
||||
doublereal tol)
|
||||
{
|
||||
Set(tpx::SP, s, p);
|
||||
Set(tpx::PropertyPair::SP, s, p);
|
||||
setState_TR(m_sub->Temp(), 1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
|
|
@ -470,7 +470,7 @@ void PureFluidPhase::setState_SP(doublereal s, doublereal p,
|
|||
doublereal PureFluidPhase::satPressure(doublereal t) const
|
||||
{
|
||||
doublereal vsv = m_sub->v();
|
||||
Set(tpx::TV,t,vsv);
|
||||
Set(tpx::PropertyPair::TV,t,vsv);
|
||||
doublereal ps = m_sub->Ps();
|
||||
return ps;
|
||||
}
|
||||
|
|
@ -485,14 +485,14 @@ void PureFluidPhase::setState_Tsat(doublereal t, doublereal x)
|
|||
{
|
||||
setTemperature(t);
|
||||
setTPXState();
|
||||
Set(tpx::TX, t, x);
|
||||
Set(tpx::PropertyPair::TX, t, x);
|
||||
setDensity(1.0/m_sub->v());
|
||||
}
|
||||
//====================================================================================================================
|
||||
void PureFluidPhase::setState_Psat(doublereal p, doublereal x)
|
||||
{
|
||||
setTPXState();
|
||||
Set(tpx::PX, p, x);
|
||||
Set(tpx::PropertyPair::PX, p, x);
|
||||
setTemperature(m_sub->Temp());
|
||||
setDensity(1.0/m_sub->v());
|
||||
}
|
||||
|
|
|
|||
265
src/tpx/Sub.cpp
265
src/tpx/Sub.cpp
|
|
@ -42,11 +42,10 @@ const double DeltaT = 0.000001;
|
|||
/// with respect to temperature.
|
||||
double Substance::dPsdT()
|
||||
{
|
||||
double ps1, tsave, dpdt;
|
||||
tsave = T;
|
||||
ps1 = Ps();
|
||||
double tsave = T;
|
||||
double ps1 = Ps();
|
||||
set_T(T + DeltaT);
|
||||
dpdt = (Ps() - ps1)/DeltaT;
|
||||
double dpdt = (Ps() - ps1)/DeltaT;
|
||||
set_T(tsave);
|
||||
return dpdt;
|
||||
}
|
||||
|
|
@ -66,7 +65,6 @@ int Substance::TwoPhase()
|
|||
/// returned if v > Vcrit.
|
||||
double Substance::x()
|
||||
{
|
||||
double vv, vl;
|
||||
if (T >= Tcrit()) {
|
||||
return (1.0/Rho < Vcrit() ? 0.0 : 1.0);
|
||||
} else {
|
||||
|
|
@ -76,8 +74,8 @@ double Substance::x()
|
|||
} else if (Rho >= Rhf) {
|
||||
return 0.0;
|
||||
} else {
|
||||
vv = 1.0/Rhv;
|
||||
vl = 1.0/Rhf;
|
||||
double vv = 1.0/Rhv;
|
||||
double vl = 1.0/Rhf;
|
||||
return (1.0/Rho - vl)/(vv - vl);
|
||||
}
|
||||
}
|
||||
|
|
@ -86,20 +84,19 @@ double Substance::x()
|
|||
/// Saturation temperature at pressure p.
|
||||
double Substance::Tsat(double p)
|
||||
{
|
||||
double Tsave, p_here, dp, dt, dpdt, dta,
|
||||
dtm, tsat;
|
||||
if (p <= 0.0 || p > Pcrit()) {
|
||||
throw TPX_Error("Substance::Tsat", "illegal pressure value");
|
||||
}
|
||||
int LoopCount = 0;
|
||||
double tol = 1.e-6*p;
|
||||
Tsave = T;
|
||||
double Tsave = T;
|
||||
if (T < Tmin()) {
|
||||
T = 0.5*(Tcrit() - Tmin());
|
||||
}
|
||||
if (T >= Tcrit()) {
|
||||
T = 0.5*(Tcrit() - Tmin());
|
||||
}
|
||||
double dp;
|
||||
do {
|
||||
if (T > Tcrit()) {
|
||||
T = Tcrit() - 0.001;
|
||||
|
|
@ -107,12 +104,10 @@ double Substance::Tsat(double p)
|
|||
if (T < Tmin()) {
|
||||
T = Tmin() + 0.001;
|
||||
}
|
||||
p_here = Ps();
|
||||
dpdt = dPsdT();
|
||||
dp = p - p_here;
|
||||
dt = dp/dpdt;
|
||||
dta = fabs(dt);
|
||||
dtm = 0.1*T;
|
||||
dp = p - Ps();
|
||||
double dt = dp/dPsdT();
|
||||
double dta = fabs(dt);
|
||||
double dtm = 0.1*T;
|
||||
if (dta > dtm) {
|
||||
dt = dt*dtm/dta;
|
||||
}
|
||||
|
|
@ -123,7 +118,7 @@ double Substance::Tsat(double p)
|
|||
throw TPX_Error("Substance::Tsat", "No convergence");
|
||||
}
|
||||
} while (fabs(dp) > tol);
|
||||
tsat = T;
|
||||
double tsat = T;
|
||||
T = Tsave;
|
||||
return tsat;
|
||||
}
|
||||
|
|
@ -139,94 +134,103 @@ static const double TolAbsV = 1.e-8;
|
|||
static const double TolAbsT = 1.e-3;
|
||||
static const double TolRel = 3.e-8;
|
||||
|
||||
void Substance::Set(int XY, double x0, double y0)
|
||||
void Substance::Set(PropertyPair::type XY, double x0, double y0)
|
||||
{
|
||||
double temp;
|
||||
|
||||
/* if inverted (PT) switch order and change sign of XY (TP = -PT) */
|
||||
if (XY < 0) {
|
||||
double tmp = x0;
|
||||
x0 = y0;
|
||||
y0 = tmp;
|
||||
XY *= -1;
|
||||
std::swap(x0, y0);
|
||||
XY = static_cast<PropertyPair::type>(-XY);
|
||||
}
|
||||
|
||||
switch (XY) {
|
||||
case TV:
|
||||
case PropertyPair::TV:
|
||||
set_T(x0);
|
||||
set_v(y0);
|
||||
break;
|
||||
|
||||
case HP:
|
||||
if (Lever(Pgiven, y0, x0, EvalH)) {
|
||||
case PropertyPair::HP:
|
||||
if (Lever(Pgiven, y0, x0, propertyFlag::H)) {
|
||||
return;
|
||||
}
|
||||
set_xy(EvalH, EvalP, x0, y0, TolAbsH, TolAbsP, TolRel, TolRel);
|
||||
set_xy(propertyFlag::H, propertyFlag::P,
|
||||
x0, y0, TolAbsH, TolAbsP, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case SP:
|
||||
if (Lever(Pgiven, y0, x0, EvalS)) {
|
||||
case PropertyPair::SP:
|
||||
if (Lever(Pgiven, y0, x0, propertyFlag::S)) {
|
||||
return;
|
||||
}
|
||||
set_xy(EvalS, EvalP, x0, y0, TolAbsS, TolAbsP, TolRel, TolRel);
|
||||
set_xy(propertyFlag::S, propertyFlag::P,
|
||||
x0, y0, TolAbsS, TolAbsP, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case PV:
|
||||
if (Lever(Pgiven, x0, y0, EvalV)) {
|
||||
case PropertyPair::PV:
|
||||
if (Lever(Pgiven, x0, y0, propertyFlag::V)) {
|
||||
return;
|
||||
}
|
||||
set_xy(EvalP, EvalV, x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
|
||||
set_xy(propertyFlag::P, propertyFlag::V,
|
||||
x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case TP:
|
||||
case PropertyPair::TP:
|
||||
if (x0 < Tcrit()) {
|
||||
set_T(x0);
|
||||
if (y0 < Ps()) {
|
||||
Set(TX, x0, Vapor);
|
||||
Set(PropertyPair::TX, x0, 1.0);
|
||||
} else {
|
||||
Set(TX, x0, Liquid);
|
||||
Set(PropertyPair::TX, x0, 0.0);
|
||||
}
|
||||
} else {
|
||||
set_T(x0);
|
||||
}
|
||||
set_xy(EvalT, EvalP, x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);
|
||||
set_xy(propertyFlag::T, propertyFlag::P,
|
||||
x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case UV:
|
||||
set_xy(EvalU, EvalV, x0, y0, TolAbsU, TolAbsV, TolRel, TolRel);
|
||||
case PropertyPair::UV:
|
||||
set_xy(propertyFlag::U, propertyFlag::V,
|
||||
x0, y0, TolAbsU, TolAbsV, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case ST:
|
||||
if (Lever(Tgiven, y0, x0, EvalS)) {
|
||||
case PropertyPair::ST:
|
||||
if (Lever(Tgiven, y0, x0, propertyFlag::S)) {
|
||||
return;
|
||||
}
|
||||
set_xy(EvalS, EvalT, x0, y0, TolAbsS, TolAbsT, TolRel, TolRel);
|
||||
set_xy(propertyFlag::S, propertyFlag::T,
|
||||
x0, y0, TolAbsS, TolAbsT, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case SV:
|
||||
set_xy(EvalS, EvalV, x0, y0, TolAbsS, TolAbsV, TolRel, TolRel);
|
||||
case PropertyPair::SV:
|
||||
set_xy(propertyFlag::S, propertyFlag::V,
|
||||
x0, y0, TolAbsS, TolAbsV, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case UP:
|
||||
if (Lever(Pgiven, y0, x0, EvalU)) {
|
||||
case PropertyPair::UP:
|
||||
if (Lever(Pgiven, y0, x0, propertyFlag::U)) {
|
||||
return;
|
||||
}
|
||||
set_xy(EvalU, EvalP, x0, y0, TolAbsU, TolAbsP, TolRel, TolRel);
|
||||
set_xy(propertyFlag::U, propertyFlag::P,
|
||||
x0, y0, TolAbsU, TolAbsP, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case VH:
|
||||
set_xy(EvalV, EvalH, x0, y0, TolAbsV, TolAbsH, TolRel, TolRel);
|
||||
case PropertyPair::VH:
|
||||
set_xy(propertyFlag::V, propertyFlag::H,
|
||||
x0, y0, TolAbsV, TolAbsH, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case TH:
|
||||
set_xy(EvalT, EvalH, x0, y0, TolAbsT, TolAbsH, TolRel, TolRel);
|
||||
case PropertyPair::TH:
|
||||
set_xy(propertyFlag::T, propertyFlag::H,
|
||||
x0, y0, TolAbsT, TolAbsH, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case SH:
|
||||
set_xy(EvalS, EvalH, x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
|
||||
case PropertyPair::SH:
|
||||
set_xy(propertyFlag::S, propertyFlag::H,
|
||||
x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
|
||||
break;
|
||||
|
||||
case PX:
|
||||
case PropertyPair::PX:
|
||||
temp = Tsat(x0);
|
||||
if (y0 > 1.0 || y0 < 0.0) {
|
||||
throw TPX_Error("Substance::Set",
|
||||
|
|
@ -241,7 +245,7 @@ void Substance::Set(int XY, double x0, double y0)
|
|||
}
|
||||
break;
|
||||
|
||||
case TX:
|
||||
case PropertyPair::TX:
|
||||
if (y0 > 1.0 || y0 < 0.0) {
|
||||
throw TPX_Error("Substance::Set",
|
||||
"Invalid vapor fraction, " + fp2str(y0));
|
||||
|
|
@ -301,22 +305,10 @@ double Substance::Ps()
|
|||
}
|
||||
|
||||
// update saturated liquid and vapor densities and saturation pressure
|
||||
|
||||
void Substance::Set_meta(double phase, double pp)
|
||||
{
|
||||
if (phase == Liquid) {
|
||||
Rho = ldens(); // trial value = liquid dens
|
||||
} else {
|
||||
Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas
|
||||
}
|
||||
set_TPp(T, pp);
|
||||
}
|
||||
|
||||
void Substance::update_sat()
|
||||
{
|
||||
if ((T != Tslast) && (T < Tcrit())) {
|
||||
double Rho_save = Rho;
|
||||
double gf, gv, dg, dp, dlp, psold;
|
||||
|
||||
double pp = Psat();
|
||||
double lps = log(pp);
|
||||
|
|
@ -325,24 +317,24 @@ void Substance::update_sat()
|
|||
|
||||
for (i = 0; i<20; i++) {
|
||||
if (i==0) {
|
||||
Rho = ldens(); // trial value = liquid density
|
||||
Rho = ldens(); // trial value = liquid density
|
||||
} else {
|
||||
Rho = Rhf;
|
||||
}
|
||||
set_TPp(T,pp);
|
||||
Rhf = Rho; // sat liquid density
|
||||
Rhf = Rho; // sat liquid density
|
||||
|
||||
gf = hp() - T*sp();
|
||||
double gf = hp() - T*sp();
|
||||
if (i==0) {
|
||||
Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas
|
||||
Rho = pp*MolWt()/(8314.0*T); // trial value = ideal gas
|
||||
} else {
|
||||
Rho = Rhv;
|
||||
}
|
||||
set_TPp(T,pp);
|
||||
|
||||
Rhv = Rho; // sat vapor density
|
||||
gv = hp() - T*sp();
|
||||
dg = gv - gf;
|
||||
Rhv = Rho; // sat vapor density
|
||||
double gv = hp() - T*sp();
|
||||
double dg = gv - gf;
|
||||
|
||||
if (Rhv > Rhf) {
|
||||
std::swap(Rhv, Rhf);
|
||||
|
|
@ -352,12 +344,11 @@ void Substance::update_sat()
|
|||
if (fabs(dg) < 0.001 && Rhf > Rhv) {
|
||||
break;
|
||||
}
|
||||
dp = dg/(1.0/Rhv - 1.0/Rhf);
|
||||
psold = pp;
|
||||
double dp = dg/(1.0/Rhv - 1.0/Rhf);
|
||||
double psold = pp;
|
||||
|
||||
if (fabs(dp) > pp) {
|
||||
dlp = dg/(pp*(1.0/Rhv - 1.0/Rhf));
|
||||
lps -= dlp;
|
||||
lps -= dg/(pp*(1.0/Rhv - 1.0/Rhf));
|
||||
pp = exp(lps);
|
||||
} else {
|
||||
pp -= dp;
|
||||
|
|
@ -386,33 +377,33 @@ void Substance::update_sat()
|
|||
}
|
||||
}
|
||||
|
||||
double Substance::vprop(int ijob)
|
||||
double Substance::vprop(propertyFlag::type ijob)
|
||||
{
|
||||
switch (ijob) {
|
||||
case EvalH:
|
||||
case propertyFlag::H:
|
||||
return hp();
|
||||
case EvalS:
|
||||
case propertyFlag::S:
|
||||
return sp();
|
||||
case EvalU:
|
||||
case propertyFlag::U:
|
||||
return up();
|
||||
case EvalV:
|
||||
case propertyFlag::V:
|
||||
return vp();
|
||||
case EvalP:
|
||||
case propertyFlag::P:
|
||||
return Pp();
|
||||
default:
|
||||
throw TPX_Error("Substance::vprop", "invalid job index");
|
||||
}
|
||||
}
|
||||
|
||||
int Substance::Lever(int itp, double sat, double val, int ifunc)
|
||||
int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc)
|
||||
{
|
||||
/*
|
||||
* uses lever rule to set state in the dome. Returns 1 if in dome,
|
||||
* 0 if not, in which case state not set.
|
||||
*/
|
||||
double Valf, Valg, Tsave, Rhosave, xx, vv, psat;
|
||||
Tsave = T;
|
||||
Rhosave = Rho;
|
||||
double psat;
|
||||
double Tsave = T;
|
||||
double Rhosave = Rho;
|
||||
if (itp == Tgiven) {
|
||||
if (sat >= Tcrit()) {
|
||||
return 0;
|
||||
|
|
@ -435,13 +426,13 @@ int Substance::Lever(int itp, double sat, double val, int ifunc)
|
|||
} else {
|
||||
throw TPX_Error("Substance::Lever","general error");
|
||||
}
|
||||
Set(TX, T, Vapor);
|
||||
Valg = vprop(ifunc);
|
||||
Set(TX, T, Liquid);
|
||||
Valf = vprop(ifunc);
|
||||
Set(PropertyPair::TX, T, 1.0);
|
||||
double Valg = vprop(ifunc);
|
||||
Set(PropertyPair::TX, T, 0.0);
|
||||
double Valf = vprop(ifunc);
|
||||
if (val >= Valf && val <= Valg) {
|
||||
xx = (val - Valf)/(Valg - Valf);
|
||||
vv = (1.0 - xx)/Rhf + xx/Rhv;
|
||||
double xx = (val - Valf)/(Valg - Valf);
|
||||
double vv = (1.0 - xx)/Rhf + xx/Rhv;
|
||||
set_v(vv);
|
||||
return 1;
|
||||
} else {
|
||||
|
|
@ -452,14 +443,12 @@ int Substance::Lever(int itp, double sat, double val, int ifunc)
|
|||
}
|
||||
|
||||
|
||||
void Substance::set_xy(int ifx, int ify, double X, double Y,
|
||||
void Substance::set_xy(propertyFlag::type ifx, propertyFlag::type ify,
|
||||
double X, double Y,
|
||||
double atx, double aty,
|
||||
double rtx, double rty)
|
||||
{
|
||||
|
||||
double v_here, t_here, dv, dt, dxdt, dydt, dxdv, dydv,
|
||||
det, x_here, y_here, dvm, dtm, dva, dta;
|
||||
double Xa, Ya, err_x, err_y;
|
||||
double v_here, t_here;
|
||||
double dvs1 = 2.0*Vcrit();
|
||||
double dvs2 = 0.7*Vcrit();
|
||||
int LoopCount = 0;
|
||||
|
|
@ -467,12 +456,14 @@ void Substance::set_xy(int ifx, int ify, double X, double Y,
|
|||
double v_save = 1.0/Rho;
|
||||
double t_save = T;
|
||||
|
||||
if ((T == Undef) && (Rho == Undef)) { // new object, try to pick
|
||||
Set(TV,Tcrit()*1.1,Vcrit()*1.1); // "reasonable" starting point
|
||||
if ((T == Undef) && (Rho == Undef)) {
|
||||
// new object, try to pick a "reasonable" starting point
|
||||
Set(PropertyPair::TV,Tcrit()*1.1,Vcrit()*1.1);
|
||||
t_here = T;
|
||||
v_here = 1.0/Rho;
|
||||
} else if (Rho == Undef) { // new object, try to pick
|
||||
Set(TV,T,Vcrit()*1.1); // "reasonable" starting point
|
||||
} else if (Rho == Undef) {
|
||||
// new object, try to pick a "reasonable" starting point
|
||||
Set(PropertyPair::TV,T,Vcrit()*1.1);
|
||||
t_here = T;
|
||||
v_here = 1.0/Rho;
|
||||
} else {
|
||||
|
|
@ -480,56 +471,55 @@ void Substance::set_xy(int ifx, int ify, double X, double Y,
|
|||
t_here = t_save;
|
||||
}
|
||||
|
||||
Xa = fabs(X);
|
||||
Ya = fabs(Y);
|
||||
double Xa = fabs(X);
|
||||
double Ya = fabs(Y);
|
||||
|
||||
// loop
|
||||
do {
|
||||
x_here = prop(ifx);
|
||||
y_here = prop(ify);
|
||||
err_x = fabs(X - x_here);
|
||||
err_y = fabs(Y - y_here);
|
||||
while (true) {
|
||||
double x_here = prop(ifx);
|
||||
double y_here = prop(ify);
|
||||
double err_x = fabs(X - x_here);
|
||||
double err_y = fabs(Y - y_here);
|
||||
|
||||
if ((err_x < atx + rtx*Xa) && (err_y < aty + rty*Ya)) {
|
||||
break;
|
||||
}
|
||||
|
||||
/* perturb t */
|
||||
dt = 0.001*t_here;
|
||||
double dt = 0.001*t_here;
|
||||
if (t_here + dt > Tmax()) {
|
||||
dt *= -1.0;
|
||||
}
|
||||
|
||||
/* perturb v */
|
||||
dv = 0.001*v_here;
|
||||
double dv = 0.001*v_here;
|
||||
if (v_here <= Vcrit()) {
|
||||
dv *= -1.0;
|
||||
}
|
||||
|
||||
/* derivatives with respect to T */
|
||||
Set(TV, t_here + dt, v_here);
|
||||
dxdt = (prop(ifx) - x_here)/dt;
|
||||
dydt = (prop(ify) - y_here)/dt;
|
||||
Set(PropertyPair::TV, t_here + dt, v_here);
|
||||
double dxdt = (prop(ifx) - x_here)/dt;
|
||||
double dydt = (prop(ify) - y_here)/dt;
|
||||
|
||||
/* derivatives with respect to v */
|
||||
Set(TV, t_here, v_here + dv);
|
||||
dxdv = (prop(ifx) - x_here)/dv;
|
||||
dydv = (prop(ify) - y_here)/dv;
|
||||
Set(PropertyPair::TV, t_here, v_here + dv);
|
||||
double dxdv = (prop(ifx) - x_here)/dv;
|
||||
double dydv = (prop(ify) - y_here)/dv;
|
||||
|
||||
det = dxdt*dydv - dydt*dxdv;
|
||||
double det = dxdt*dydv - dydt*dxdv;
|
||||
dt = ((X - x_here)*dydv - (Y - y_here)*dxdv)/det;
|
||||
dv = ((Y - y_here)*dxdt - (X - x_here)*dydt)/det;
|
||||
|
||||
dvm = 0.2*v_here;
|
||||
double dvm = 0.2*v_here;
|
||||
if (v_here < dvs1) {
|
||||
dvm *= 0.5;
|
||||
}
|
||||
if (v_here < dvs2) {
|
||||
dvm *= 0.5;
|
||||
}
|
||||
dtm = 0.1*t_here;
|
||||
dva = fabs(dv);
|
||||
dta = fabs(dt);
|
||||
double dtm = 0.1*t_here;
|
||||
double dva = fabs(dv);
|
||||
double dta = fabs(dt);
|
||||
if (dva > dvm) {
|
||||
dv *= dvm/dva;
|
||||
}
|
||||
|
|
@ -546,32 +536,31 @@ void Substance::set_xy(int ifx, int ify, double X, double Y,
|
|||
if (v_here <= 0.0) {
|
||||
v_here = 0.0001;
|
||||
}
|
||||
Set(TV, t_here, v_here);
|
||||
Set(PropertyPair::TV, t_here, v_here);
|
||||
LoopCount++;
|
||||
if (LoopCount > 200) {
|
||||
throw TPX_Error("Substance::set_xy","no convergence");
|
||||
}
|
||||
} while (1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double Substance::prop(int ijob)
|
||||
double Substance::prop(propertyFlag::type ijob)
|
||||
{
|
||||
double xx, pp, lp, vp, Rho_save;
|
||||
if (ijob == EvalP) {
|
||||
if (ijob == propertyFlag::P) {
|
||||
return P();
|
||||
}
|
||||
if (ijob == EvalT) {
|
||||
if (ijob == propertyFlag::T) {
|
||||
return T;
|
||||
}
|
||||
xx = x();
|
||||
double xx = x();
|
||||
if ((xx > 0.0) && (xx < 1.0)) {
|
||||
Rho_save = Rho;
|
||||
double Rho_save = Rho;
|
||||
Rho = Rhv;
|
||||
vp = vprop(ijob);
|
||||
double vp = vprop(ijob);
|
||||
Rho = Rhf;
|
||||
lp = vprop(ijob);
|
||||
pp = (1.0 - xx)*lp + xx*vp;
|
||||
double lp = vprop(ijob);
|
||||
double pp = (1.0 - xx)*lp + xx*vp;
|
||||
Rho = Rho_save;
|
||||
return pp;
|
||||
} else {
|
||||
|
|
@ -627,7 +616,7 @@ void Substance::set_TPp(double Temp, double Pressure)
|
|||
if (v_here <= Vcrit()) {
|
||||
dv *= -1.0;
|
||||
}
|
||||
Set(TV, Temp, v_here+dv);
|
||||
Set(PropertyPair::TV, Temp, v_here+dv);
|
||||
double dpdv = (Pp() - P_here)/dv;
|
||||
if (dpdv > 0.0) {
|
||||
BracketSlope(Pressure);
|
||||
|
|
@ -678,15 +667,15 @@ void Substance::set_TPp(double Temp, double Pressure)
|
|||
if (dv == 0.0) {
|
||||
throw TPX_Error("Substance::set_TPp","dv = 0 and no convergence");
|
||||
}
|
||||
Set(TV, Temp, v_here);
|
||||
Set(PropertyPair::TV, Temp, v_here);
|
||||
LoopCount++;
|
||||
if (LoopCount > 100) {
|
||||
Set(TV, Temp, v_save);
|
||||
Set(PropertyPair::TV, Temp, v_save);
|
||||
throw TPX_Error("Substance::set_TPp",string("no convergence for ")
|
||||
+"P* = "+fp2str(Pressure/Pcrit())+". V* = "
|
||||
+fp2str(v_save/Vcrit()));
|
||||
}
|
||||
}
|
||||
Set(TV, Temp,v_here);
|
||||
Set(PropertyPair::TV, Temp,v_here);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue