github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-13 07:15:14 -05:00
Code Issues Packages Projects Releases Wiki Activity

Replaced all tabs with spaces (finally) in C++ files

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-11-19 12:16:14 -05:00
parent 4249dc8f4a
commit 8327d54ea2
59 changed files with 8789 additions and 8789 deletions
Download Patch File Download Diff File Expand all files Collapse all files

238
src/Backends/Helmholtz/ExcessHEFunction.h
View File

@@ -18,17 +18,17 @@ typedef std::vector<std::vector<long double> > STLMatrix;
class DepartureFunction
{
public:
DepartureFunction(){};
virtual ~DepartureFunction(){};
DepartureFunction(){};
virtual ~DepartureFunction(){};
/// The excess Helmholtz energy of the binary pair
/// Pure-virtual function (must be implemented in derived class
virtual double alphar(double tau, double delta) = 0;
virtual double dalphar_dDelta(double tau, double delta) = 0;
virtual double d2alphar_dDelta2(double tau, double delta) = 0;
virtual double d2alphar_dDelta_dTau(double tau, double delta) = 0;
virtual double dalphar_dTau(double tau, double delta) = 0;
virtual double d2alphar_dTau2(double tau, double delta) = 0;
/// The excess Helmholtz energy of the binary pair
/// Pure-virtual function (must be implemented in derived class
virtual double alphar(double tau, double delta) = 0;
virtual double dalphar_dDelta(double tau, double delta) = 0;
virtual double d2alphar_dDelta2(double tau, double delta) = 0;
virtual double d2alphar_dDelta_dTau(double tau, double delta) = 0;
virtual double dalphar_dTau(double tau, double delta) = 0;
virtual double d2alphar_dTau2(double tau, double delta) = 0;
};
/** \brief The departure function used by the GERG-2008 formulation
@@ -42,10 +42,10 @@ public:
class GERG2008DepartureFunction : public DepartureFunction
{
protected:
bool using_gaussian;
ResidualHelmholtzGeneralizedExponential phi;
bool using_gaussian;
ResidualHelmholtzGeneralizedExponential phi;
public:
GERG2008DepartureFunction(){};
GERG2008DepartureFunction(){};
GERG2008DepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,
const std::vector<double> &eta,const std::vector<double> &epsilon,const std::vector<double> &beta,
const std::vector<double> &gamma, std::size_t Npower)
@@ -75,14 +75,14 @@ public:
phi.add_GERG2008Gaussian(_n, _d, _t, _eta, _epsilon, _beta, _gamma);
}
};
~GERG2008DepartureFunction(){};
~GERG2008DepartureFunction(){};
double alphar(double tau, double delta){return phi.base(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
};
/** \brief A polynomial/exponential departure function
@@ -96,9 +96,9 @@ public:
class ExponentialDepartureFunction : public DepartureFunction
{
protected:
ResidualHelmholtzGeneralizedExponential phi;
ResidualHelmholtzGeneralizedExponential phi;
public:
ExponentialDepartureFunction(){};
ExponentialDepartureFunction(){};
ExponentialDepartureFunction(const std::vector<double> &n, const std::vector<double> &d,
const std::vector<double> &t, const std::vector<double> &l)
{
@@ -108,14 +108,14 @@ public:
std::vector<long double> _l(l.begin(), l.begin()+l.size());
phi.add_Power(_n, _d, _t, _l);
};
~ExponentialDepartureFunction(){};
~ExponentialDepartureFunction(){};
double alphar(double tau, double delta){return phi.base(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
};
typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
@@ -123,12 +123,12 @@ typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
class ExcessTerm
{
public:
std::size_t N;
std::vector<std::vector<DepartureFunctionPointer> > DepartureFunctionMatrix;
std::vector<std::vector<long double> > F;
std::size_t N;
std::vector<std::vector<DepartureFunctionPointer> > DepartureFunctionMatrix;
std::vector<std::vector<long double> > F;
ExcessTerm(){};
~ExcessTerm(){};
~ExcessTerm(){};
/// Resize the parts of this term
void resize(std::size_t N){
@@ -137,127 +137,127 @@ public:
DepartureFunctionMatrix.resize(N);
for (std::size_t i = 0; i < N; ++i){
DepartureFunctionMatrix[i].resize(N);
}
}
};
double alphar(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
}
}
return summer;
}
double dalphar_dDelta(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta(tau,delta);
}
}
return summer;
}
double d2alphar_dDelta2(double tau, double delta, const std::vector<long double> &x)
double alphar(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2(tau,delta);
}
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
}
}
return summer;
};
double d2alphar_dDelta_dTau(double tau, double delta, const std::vector<long double> &x)
}
double dalphar_dDelta(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau(tau,delta);
}
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta(tau,delta);
}
}
return summer;
}
double dalphar_dTau(double tau, double delta, const std::vector<long double> &x)
double d2alphar_dDelta2(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau(tau,delta);
}
}
return summer;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2(tau,delta);
}
}
return summer;
};
double d2alphar_dTau2(double tau, double delta, const std::vector<long double> &x)
double d2alphar_dDelta_dTau(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dTau2(tau,delta);
}
}
return summer;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau(tau,delta);
}
}
return summer;
}
double dalphar_dTau(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau(tau,delta);
}
}
return summer;
};
double dalphar_dxi(double tau, double delta, const std::vector<long double> &x, std::size_t i)
double d2alphar_dTau2(double tau, double delta, const std::vector<long double> &x)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->alphar(tau,delta);
}
}
return summer;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dTau2(tau,delta);
}
}
return summer;
};
double dalphar_dxi(double tau, double delta, const std::vector<long double> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->alphar(tau,delta);
}
}
return summer;
};
double d2alphardxidxj(double tau, double delta, const std::vector<long double> &x, std::size_t i, std::size_t j)
{
if (i != j)
{
return F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
}
else
{
return 0;
}
{
return F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
}
else
{
return 0;
}
};
double d2alphar_dxi_dTau(double tau, double delta, const std::vector<long double> &x, std::size_t i)
double d2alphar_dxi_dTau(double tau, double delta, const std::vector<long double> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau(tau,delta);
}
}
return summer;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau(tau,delta);
}
}
return summer;
};
double d2alphar_dxi_dDelta(double tau, double delta, const std::vector<long double> &x, std::size_t i)
double d2alphar_dxi_dDelta(double tau, double delta, const std::vector<long double> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta(tau,delta);
}
}
return summer;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta(tau,delta);
}
}
return summer;
};
};

130
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -8,25 +8,25 @@ namespace CoolProp{
template<class T> T g_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// g function from Rachford-Rice
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += z[i]*(Ki-1)/(1-beta+beta*Ki);
}
return summer;
// g function from Rachford-Rice
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += z[i]*(Ki-1)/(1-beta+beta*Ki);
}
return summer;
}
template<class T> T dgdbeta_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// derivative of g function from Rachford-Rice with respect to beta
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += -z[i]*pow((Ki-1)/(1-beta+beta*Ki),2);
}
return summer;
// derivative of g function from Rachford-Rice with respect to beta
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += -z[i]*pow((Ki-1)/(1-beta+beta*Ki),2);
}
return summer;
}
void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
@@ -104,8 +104,8 @@ void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
}
void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function)
{
if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function)
{
if (HEOS.is_pure_or_pseudopure)
{
// At very low temperature (near the triple point temp), the isotherms are VERY steep
@@ -127,14 +127,14 @@ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS)
throw ValueError("twophase not implemented yet");
}
}
else{
else{
PT_flash_mixtures(HEOS);
}
}
// Find density
HEOS._rhomolar = HEOS.solver_rho_Tp(HEOS._T, HEOS._p);
HEOS._Q = -1;
}
// Find density
HEOS._rhomolar = HEOS.solver_rho_Tp(HEOS._T, HEOS._p);
HEOS._Q = -1;
}
void FlashRoutines::HQ_flash(HelmholtzEOSMixtureBackend &HEOS, long double Tguess)
{
@@ -212,15 +212,15 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
long double T = HEOS._T;
if (HEOS.is_pure_or_pseudopure)
{
// The maximum possible saturation temperature
// Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
long double Tmax_sat = HEOS.calc_Tmax_sat() + 1e-13;
// Check what the minimum limits for the equation of state are
long double Tmin_satL, Tmin_satV, Tmin_sat;
HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
// The maximum possible saturation temperature
// Critical point for pure fluids, slightly different for pseudo-pure, very different for mixtures
long double Tmax_sat = HEOS.calc_Tmax_sat() + 1e-13;
// Check what the minimum limits for the equation of state are
long double Tmin_satL, Tmin_satV, Tmin_sat;
HEOS.calc_Tmin_sat(Tmin_satL, Tmin_satV);
Tmin_sat = std::max(Tmin_satL, Tmin_satV) - 1e-13;
// Get a reference to keep the code a bit cleaner
CriticalRegionSplines &splines = HEOS.components[0]->pEOS->critical_region_splines;
@@ -232,8 +232,8 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
HEOS._p = HEOS.SatL->p();
}
else if (!is_in_closed_range(Tmin_sat, Tmax_sat, T)){
throw ValueError(format("Temperature to QT_flash [%0.8Lg K] must be in range [%0.8Lg K, %0.8Lg K]", T, Tmin_sat, Tmax_sat));
}
throw ValueError(format("Temperature to QT_flash [%0.8Lg K] must be in range [%0.8Lg K, %0.8Lg K]", T, Tmin_sat, Tmax_sat));
}
else if (get_config_bool(CRITICAL_SPLINES_ENABLED) && splines.enabled && HEOS._T > splines.T_min){
double rhoL = _HUGE, rhoV = _HUGE;
// Use critical region spline if it has it and temperature is in its range
@@ -248,7 +248,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
// Set some imput options
SaturationSolvers::saturation_T_pure_options options;
options.use_guesses = false;
double increment = 0.2;
double increment = 0.2;
try{
for (double omega = 1.0; omega > 0; omega -= increment){
try{
@@ -401,8 +401,8 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
options.specified_variable = SaturationSolvers::saturation_PHSU_pure_options::IMPOSED_PL;
// Use logarithm of delta as independent variables
options.use_logdelta = false;
double increment = 0.4;
double increment = 0.4;
try{
for (double omega = 1.0; omega > 0; omega -= increment){
@@ -929,7 +929,7 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE
{
if (!ValidNumber(HEOS._p)){throw ValueError("value for p in HSU_P_flash_singlephase_Brent is invalid");};
if (!ValidNumber(value)){throw ValueError("value for other in HSU_P_flash_singlephase_Brent is invalid");};
class solver_resid : public FuncWrapper1D
class solver_resid : public FuncWrapper1D
{
public:
@@ -949,20 +949,20 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE
};
double call(double T){
this->T = T;
this->T = T;
// Run the solver with T,P as inputs;
HEOS->update(PT_INPUTS, p, T);
// Run the solver with T,P as inputs;
HEOS->update(PT_INPUTS, p, T);
rhomolar = HEOS->rhomolar();
HEOS->update(DmolarT_INPUTS, rhomolar, T);
// Get the value of the desired variable
eos = HEOS->keyed_output(other);
// Get the value of the desired variable
eos = HEOS->keyed_output(other);
pp = HEOS->p();
// Difference between the two is to be driven to zero
// Difference between the two is to be driven to zero
r = eos - value;
// Store values for later use if there are errors
if (iter == 0){
r0 = r; T0 = T; eos0 = eos;
@@ -979,9 +979,9 @@ void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HE
return r;
};
};
solver_resid resid(&HEOS, HEOS._p, value, other);
std::string errstr;
solver_resid resid(&HEOS, HEOS._p, value, other);
std::string errstr;
try{
Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-12, 100, errstr);
// Un-specify the phase of the fluid
@@ -1036,8 +1036,8 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
// Find the phase, while updating all internal variables possible
HEOS.p_phase_determination_pure_or_pseudopure(other, value, saturation_called);
if (HEOS.isHomogeneousPhase())
{
if (HEOS.isHomogeneousPhase())
{
// Now we use the single-phase solver to find T,rho given P,Y using a
// bounded 1D solver by adjusting T and using given value of p
long double Tmin, Tmax;
@@ -1052,14 +1052,14 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
case iphase_liquid:
{
if (saturation_called){ Tmax = HEOS.SatL->T();}else{Tmax = HEOS._TLanc.pt();}
// Sometimes the minimum pressure for the melting line is a bit above the triple point pressure
if (HEOS.has_melting_line() && HEOS._p > HEOS.calc_melting_line(iP_min, -1, -1)){
if (HEOS.has_melting_line() && HEOS._p > HEOS.calc_melting_line(iP_min, -1, -1)){
Tmin = HEOS.calc_melting_line(iT, iP, HEOS._p)-1e-3;
}
else{
Tmin = HEOS.Tmin()-1e-3;
}
}
else{
Tmin = HEOS.Tmin()-1e-3;
}
break;
}
case iphase_supercritical_liquid:
@@ -1069,11 +1069,11 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
Tmax = 1.5*HEOS.Tmax();
// Sometimes the minimum pressure for the melting line is a bit above the triple point pressure
if (HEOS.has_melting_line() && HEOS._p > HEOS.calc_melting_line(iP_min, -1, -1)){
Tmin = HEOS.calc_melting_line(iT, iP, HEOS._p)-1e-3;
}
else{
Tmin = HEOS.Tmin()-1e-3;
}
Tmin = HEOS.calc_melting_line(iT, iP, HEOS._p)-1e-3;
}
else{
Tmin = HEOS.Tmin()-1e-3;
}
break;
}
default:
@@ -1085,8 +1085,8 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, parameters oth
catch(std::exception &e){
throw ValueError(format("unable to solve 1phase PY flash with Tmin=%Lg, Tmax=%Lg due to error: %s",Tmin, Tmax, e.what()));
}
HEOS._Q = -1;
}
HEOS._Q = -1;
}
}
else
{

4
src/Backends/Helmholtz/FlashRoutines.h
View File

@@ -75,8 +75,8 @@ public:
/// @param T0 The initial guess value for the temperature [K]
/// @param rhomolar0 The initial guess value for the density [mol/m^3]
static void HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend &HEOS, parameters other, long double T0, long double rhomolar0);
/// The single-phase flash routine for the pairs (P,H), (P,S), and (P,U). Similar analysis is needed
/// The single-phase flash routine for the pairs (P,H), (P,S), and (P,U). Similar analysis is needed
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
/// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
/// @param value The value of the other input

364
src/Backends/Helmholtz/Fluids/Ancillaries.cpp
View File

@@ -13,12 +13,12 @@ namespace CoolProp{
SaturationAncillaryFunction::SaturationAncillaryFunction(rapidjson::Value &json_code)
{
std::string type = cpjson::get_string(json_code,"type");
if (!type.compare("rational_polynomial"))
{
num_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["A"]));
den_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["B"]));
max_abs_error = cpjson::get_double(json_code,"max_abs_error");
std::string type = cpjson::get_string(json_code,"type");
if (!type.compare("rational_polynomial"))
{
num_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["A"]));
den_coeffs = vec_to_eigen(cpjson::get_double_array(json_code["B"]));
max_abs_error = cpjson::get_double(json_code,"max_abs_error");
try{
Tmin = cpjson::get_double(json_code,"Tmin");
Tmax = cpjson::get_double(json_code,"Tmax");
@@ -27,103 +27,103 @@ SaturationAncillaryFunction::SaturationAncillaryFunction(rapidjson::Value &json_
Tmin = _HUGE;
Tmax = _HUGE;
}
}
else
{
n = cpjson::get_double_array(json_code["n"]);
t = cpjson::get_double_array(json_code["t"]);
Tmin = cpjson::get_double(json_code,"Tmin");
Tmax = cpjson::get_double(json_code,"Tmax");
reducing_value = cpjson::get_double(json_code,"reducing_value");
using_tau_r = cpjson::get_bool(json_code,"using_tau_r");
T_r = cpjson::get_double(json_code,"T_r");
}
if (!type.compare("rational_polynomial"))
this->type = TYPE_RATIONAL_POLYNOMIAL;
else if (!type.compare("rhoLnoexp"))
this->type = TYPE_NOT_EXPONENTIAL;
else
this->type = TYPE_EXPONENTIAL;
this->N = n.size();
s = n;
}
else
{
n = cpjson::get_double_array(json_code["n"]);
t = cpjson::get_double_array(json_code["t"]);
Tmin = cpjson::get_double(json_code,"Tmin");
Tmax = cpjson::get_double(json_code,"Tmax");
reducing_value = cpjson::get_double(json_code,"reducing_value");
using_tau_r = cpjson::get_bool(json_code,"using_tau_r");
T_r = cpjson::get_double(json_code,"T_r");
}
if (!type.compare("rational_polynomial"))
this->type = TYPE_RATIONAL_POLYNOMIAL;
else if (!type.compare("rhoLnoexp"))
this->type = TYPE_NOT_EXPONENTIAL;
else
this->type = TYPE_EXPONENTIAL;
this->N = n.size();
s = n;
};
double SaturationAncillaryFunction::evaluate(double T)
{
if (type == TYPE_NOT_SET)
{
throw ValueError(format("type not set"));
}
else if (type == TYPE_RATIONAL_POLYNOMIAL)
{
Polynomial2D poly;
return poly.evaluate(num_coeffs, T)/poly.evaluate(den_coeffs, T);
}
else
{
double THETA = 1-T/T_r;
if (type == TYPE_NOT_SET)
{
throw ValueError(format("type not set"));
}
else if (type == TYPE_RATIONAL_POLYNOMIAL)
{
Polynomial2D poly;
return poly.evaluate(num_coeffs, T)/poly.evaluate(den_coeffs, T);
}
else
{
double THETA = 1-T/T_r;
for (std::size_t i = 0; i < N; ++i)
{
s[i] = n[i]*pow(THETA, t[i]);
}
double summer = std::accumulate(s.begin(), s.end(), 0.0);
for (std::size_t i = 0; i < N; ++i)
{
s[i] = n[i]*pow(THETA, t[i]);
}
double summer = std::accumulate(s.begin(), s.end(), 0.0);
if (type == TYPE_NOT_EXPONENTIAL)
{
return reducing_value*(1+summer);
}
else
{
double tau_r_value;
if (using_tau_r)
tau_r_value = T_r/T;
else
tau_r_value = 1.0;
return reducing_value*exp(tau_r_value*summer);
}
}
if (type == TYPE_NOT_EXPONENTIAL)
{
return reducing_value*(1+summer);
}
else
{
double tau_r_value;
if (using_tau_r)
tau_r_value = T_r/T;
else
tau_r_value = 1.0;
return reducing_value*exp(tau_r_value*summer);
}
}
}
double SaturationAncillaryFunction::invert(double value, double min_bound, double max_bound)
{
// Invert the ancillary curve to get the temperature as a function of the output variable
// Define the residual to be driven to zero
class solver_resid : public FuncWrapper1D
{
public:
int other;
SaturationAncillaryFunction *anc;
long double T, value, r, current_value;
// Invert the ancillary curve to get the temperature as a function of the output variable
// Define the residual to be driven to zero
class solver_resid : public FuncWrapper1D
{
public:
int other;
SaturationAncillaryFunction *anc;
long double T, value, r, current_value;
solver_resid(SaturationAncillaryFunction *anc, long double value) : anc(anc), value(value){};
solver_resid(SaturationAncillaryFunction *anc, long double value) : anc(anc), value(value){};
double call(double T){
this->T = T;
current_value = anc->evaluate(T);
r = current_value - value;
return r;
};
};
solver_resid resid(this, value);
std::string errstring;
double call(double T){
this->T = T;
current_value = anc->evaluate(T);
r = current_value - value;
return r;
};
};
solver_resid resid(this, value);
std::string errstring;
if (min_bound < 0){ min_bound = Tmin-0.01;}
if (max_bound < 0){ max_bound = Tmax;}
try{
try{
// Safe to expand the domain a little bit to lower temperature, absolutely cannot exceed Tmax
// because then you get (negative number)^(double) which is undefined.
return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-10,100,errstring);
}
catch(std::exception &e){
return Secant(resid,max_bound, -0.01, 1e-12, 100, errstring);
}
return Brent(resid,min_bound,max_bound,DBL_EPSILON,1e-10,100,errstring);
}
catch(std::exception &e){
return Secant(resid,max_bound, -0.01, 1e-12, 100, errstring);
}
}
void MeltingLineVariables::set_limits(void)
{
if (type == MELTING_LINE_SIMON_TYPE){
if (type == MELTING_LINE_SIMON_TYPE){
// Fill in the min and max pressures for each part
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
@@ -131,11 +131,11 @@ void MeltingLineVariables::set_limits(void)
part.p_max = part.p_0 + part.a*(pow(part.T_max/part.T_0,part.c)-1);
}
pmin = simon.parts.front().p_min;
pmax = simon.parts.back().p_max;
pmax = simon.parts.back().p_max;
Tmin = simon.parts.front().T_min;
Tmax = simon.parts.back().T_max;
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){
Tmax = simon.parts.back().T_max;
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){
// Fill in the min and max pressures for each part
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
@@ -144,10 +144,10 @@ void MeltingLineVariables::set_limits(void)
}
Tmin = polynomial_in_Tr.parts.front().T_min;
pmin = polynomial_in_Tr.parts.front().p_min;
Tmax = polynomial_in_Tr.parts.back().T_max;
Tmax = polynomial_in_Tr.parts.back().T_max;
pmax = polynomial_in_Tr.parts.back().p_max;
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){
// Fill in the min and max pressures for each part
for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i){
MeltingLinePiecewisePolynomialInThetaSegment &part = polynomial_in_Theta.parts[i];
@@ -156,90 +156,90 @@ void MeltingLineVariables::set_limits(void)
}
Tmin = polynomial_in_Theta.parts.front().T_min;
pmin = polynomial_in_Theta.parts.front().p_min;
Tmax = polynomial_in_Theta.parts.back().T_max;
Tmax = polynomial_in_Theta.parts.back().T_max;
pmax = polynomial_in_Theta.parts.back().p_max;
}
else{
throw ValueError("only Simon supported now");
}
}
else{
throw ValueError("only Simon supported now");
}
}
long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
{
if (type == MELTING_LINE_NOT_SET){throw ValueError("Melting line curve not set");}
if (type == MELTING_LINE_NOT_SET){throw ValueError("Melting line curve not set");}
if (OF == iP_max){ return pmax;}
else if (OF == iP_min){ return pmin;}
else if (OF == iT_max){ return Tmax;}
else if (OF == iT_min){ return Tmin;}
else if (OF == iP && GIVEN == iT){
long double T = value;
if (type == MELTING_LINE_SIMON_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.p_0 + part.a*(pow(T/part.T_0,part.c)-1);
}
}
throw ValueError("unable to calculate melting line (p,T) for Simon curve");
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.evaluate(T);
}
}
throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Tr curve");
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i){
MeltingLinePiecewisePolynomialInThetaSegment &part = polynomial_in_Theta.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.evaluate(T);
}
}
throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Theta curve");
}
else{
throw ValueError(format("Invalid melting line type [%d]",type));
}
}
else{
if (type == MELTING_LINE_SIMON_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
// p = part.p_0 + part.a*(pow(T/part.T_0,part.c)-1);
long double T = pow((value-part.p_0)/part.a+1,1/part.c)*part.T_0;
if (T >= part.T_0 && T <= part.T_max){
return T;
}
}
throw ValueError(format("unable to calculate melting line T(p) for Simon curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE)
{
else if (OF == iP && GIVEN == iT){
long double T = value;
if (type == MELTING_LINE_SIMON_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.p_0 + part.a*(pow(T/part.T_0,part.c)-1);
}
}
throw ValueError("unable to calculate melting line (p,T) for Simon curve");
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
MeltingLinePiecewisePolynomialInTrSegment &part = polynomial_in_Tr.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.evaluate(T);
}
}
throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Tr curve");
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i){
MeltingLinePiecewisePolynomialInThetaSegment &part = polynomial_in_Theta.parts[i];
if (is_in_closed_range(part.T_min, part.T_max, T)){
return part.evaluate(T);
}
}
throw ValueError("unable to calculate melting line (p,T) for polynomial_in_Theta curve");
}
else{
throw ValueError(format("Invalid melting line type [%d]",type));
}
}
else{
if (type == MELTING_LINE_SIMON_TYPE){
// Need to find the right segment
for (std::size_t i = 0; i < simon.parts.size(); ++i){
MeltingLinePiecewiseSimonSegment &part = simon.parts[i];
// p = part.p_0 + part.a*(pow(T/part.T_0,part.c)-1);
long double T = pow((value-part.p_0)/part.a+1,1/part.c)*part.T_0;
if (T >= part.T_0 && T <= part.T_max){
return T;
}
}
throw ValueError(format("unable to calculate melting line T(p) for Simon curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
}
else if (type == MELTING_LINE_POLYNOMIAL_IN_TR_TYPE)
{
class solver_resid : public FuncWrapper1D
{
public:
MeltingLinePiecewisePolynomialInTrSegment *part;
long double r, given_p, calc_p, T;
solver_resid(MeltingLinePiecewisePolynomialInTrSegment *part, long double p) : part(part), given_p(p){};
double call(double T){
{
public:
MeltingLinePiecewisePolynomialInTrSegment *part;
long double r, given_p, calc_p, T;
solver_resid(MeltingLinePiecewisePolynomialInTrSegment *part, long double p) : part(part), given_p(p){};
double call(double T){
this->T = T;
this->T = T;
calc_p = part->evaluate(T);
calc_p = part->evaluate(T);
// Difference between the two is to be driven to zero
r = given_p - calc_p;
// Difference between the two is to be driven to zero
r = given_p - calc_p;
return r;
};
};
return r;
};
};
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Tr.parts.size(); ++i){
@@ -257,23 +257,23 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
{
class solver_resid : public FuncWrapper1D
{
public:
MeltingLinePiecewisePolynomialInThetaSegment *part;
long double r, given_p, calc_p, T;
solver_resid(MeltingLinePiecewisePolynomialInThetaSegment *part, long double p) : part(part), given_p(p){};
double call(double T){
{
public:
MeltingLinePiecewisePolynomialInThetaSegment *part;
long double r, given_p, calc_p, T;
solver_resid(MeltingLinePiecewisePolynomialInThetaSegment *part, long double p) : part(part), given_p(p){};
double call(double T){
this->T = T;
this->T = T;
calc_p = part->evaluate(T);
calc_p = part->evaluate(T);
// Difference between the two is to be driven to zero
r = given_p - calc_p;
// Difference between the two is to be driven to zero
r = given_p - calc_p;
return r;
};
};
return r;
};
};
// Need to find the right segment
for (std::size_t i = 0; i < polynomial_in_Theta.parts.size(); ++i){
@@ -287,11 +287,11 @@ long double MeltingLineVariables::evaluate(int OF, int GIVEN, long double value)
}
throw ValueError(format("unable to calculate melting line T(p) for polynomial_in_Theta curve for p=%Lg; bounds are %Lg,%Lg Pa", value, pmin, pmax));
}
else{
throw ValueError(format("Invalid melting line type T(p) [%d]",type));
}
}
}
else{
throw ValueError(format("Invalid melting line type T(p) [%d]",type));
}
}
}
}; /* namespace CoolProp */
@@ -303,7 +303,7 @@ TEST_CASE("Water melting line", "[melting]")
int iT = CoolProp::iT, iP = CoolProp::iP;
SECTION("Ice Ih-liquid")
{
double actual = AS->melting_line(iT, iP, 138.268e6);
double actual = AS->melting_line(iT, iP, 138.268e6);
double expected = 260.0;
CAPTURE(actual);
CAPTURE(expected);
@@ -311,7 +311,7 @@ TEST_CASE("Water melting line", "[melting]")
}
SECTION("Ice III-liquid")
{
double actual = AS->melting_line(iT, iP, 268.685e6);
double actual = AS->melting_line(iT, iP, 268.685e6);
double expected = 254;
CAPTURE(actual);
CAPTURE(expected);
@@ -319,7 +319,7 @@ TEST_CASE("Water melting line", "[melting]")
}
SECTION("Ice V-liquid")
{
double actual = AS->melting_line(iT, iP, 479.640e6);
double actual = AS->melting_line(iT, iP, 479.640e6);
double expected = 265;
CAPTURE(actual);
CAPTURE(expected);
@@ -327,7 +327,7 @@ TEST_CASE("Water melting line", "[melting]")
}
SECTION("Ice VI-liquid")
{
double actual = AS->melting_line(iT, iP, 1356.76e6);
double actual = AS->melting_line(iT, iP, 1356.76e6);
double expected = 320;
CAPTURE(actual);
CAPTURE(expected);

8
src/Backends/Helmholtz/Fluids/FluidLibrary.cpp
View File

@@ -8,10 +8,10 @@ static JSONFluidLibrary library;
void load()
{
rapidjson::Document dd;
rapidjson::Document dd;
// This json formatted string comes from the all_fluids_JSON.h header which is a C++-escaped version of the JSON file
dd.Parse<0>(all_fluids_JSON.c_str());
if (dd.HasParseError()){
if (dd.HasParseError()){
throw ValueError("Unable to load all_fluids.json");
} else{
try{library.add_many(dd);}catch(std::exception &e){std::cout << e.what() << std::endl;}
@@ -19,8 +19,8 @@ void load()
}
JSONFluidLibrary & get_library(void){
if (library.is_empty()){ load(); }
return library;
if (library.is_empty()){ load(); }
return library;
}
CoolPropFluid& get_fluid(std::string fluid_string){

16
src/Backends/Helmholtz/Fluids/FluidLibrary.h
View File

@@ -886,7 +886,7 @@ protected:
if (!type.compare("Simon"))
{
rapidjson::Value &parts = melting_line["parts"];
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_SIMON_TYPE;
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_SIMON_TYPE;
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewiseSimonSegment data;
@@ -902,7 +902,7 @@ protected:
else if (!type.compare("polynomial_in_Tr"))
{
rapidjson::Value &parts = melting_line["parts"];
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_TR_TYPE;
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_TR_TYPE;
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewisePolynomialInTrSegment data;
@@ -918,7 +918,7 @@ protected:
else if (!type.compare("polynomial_in_Theta"))
{
rapidjson::Value &parts = melting_line["parts"];
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE;
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE;
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewisePolynomialInThetaSegment data;
@@ -934,8 +934,8 @@ protected:
else{
throw ValueError(format("melting line type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
// Set the limits for the melting line curve
fluid.ancillaries.melting_line.set_limits();
// Set the limits for the melting line curve
fluid.ancillaries.melting_line.set_limits();
}
else{
throw ValueError(format("melting line does not have \"type\" for fluid %s", fluid.name.c_str()));
@@ -1151,9 +1151,9 @@ public:
for (std::size_t i = 0; i < fluid.aliases.size(); ++i)
{
string_to_index_map[fluid.aliases[i]] = index;
// Add uppercase alias for EES compatibility
string_to_index_map[upper(fluid.aliases[i])] = index;
// Add uppercase alias for EES compatibility
string_to_index_map[upper(fluid.aliases[i])] = index;
}
if (get_debug_level() > 5){ std::cout << format("Loaded.\n"); }

150
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -555,65 +555,65 @@ long double HelmholtzEOSMixtureBackend::calc_rhomolar_critical(void)
}
long double HelmholtzEOSMixtureBackend::calc_pmax_sat(void)
{
if (is_pure_or_pseudopure)
{
if (components[0]->pEOS->pseudo_pure)
{
return components[0]->pEOS->max_sat_p.p;
}
else{
return p_critical();
}
}
else{
throw ValueError("calc_pmax_sat not yet defined for mixtures");
}
if (is_pure_or_pseudopure)
{
if (components[0]->pEOS->pseudo_pure)
{
return components[0]->pEOS->max_sat_p.p;
}
else{
return p_critical();
}
}
else{
throw ValueError("calc_pmax_sat not yet defined for mixtures");
}
}
long double HelmholtzEOSMixtureBackend::calc_Tmax_sat(void)
{
if (is_pure_or_pseudopure)
{
if (components[0]->pEOS->pseudo_pure)
{
return components[0]->pEOS->max_sat_T.T;
}
else{
return T_critical();
}
}
else{
throw ValueError("calc_Tmax_sat not yet defined for mixtures");
}
if (is_pure_or_pseudopure)
{
if (components[0]->pEOS->pseudo_pure)
{
return components[0]->pEOS->max_sat_T.T;
}
else{
return T_critical();
}
}
else{
throw ValueError("calc_Tmax_sat not yet defined for mixtures");
}
}
void HelmholtzEOSMixtureBackend::calc_Tmin_sat(long double &Tmin_satL, long double &Tmin_satV)
{
if (is_pure_or_pseudopure)
{
Tmin_satL = components[0]->pEOS->sat_min_liquid.T;
Tmin_satV = components[0]->pEOS->sat_min_vapor.T;
return;
}
else{
throw ValueError("calc_Tmin_sat not yet defined for mixtures");
}
if (is_pure_or_pseudopure)
{
Tmin_satL = components[0]->pEOS->sat_min_liquid.T;
Tmin_satV = components[0]->pEOS->sat_min_vapor.T;
return;
}
else{
throw ValueError("calc_Tmin_sat not yet defined for mixtures");
}
}
void HelmholtzEOSMixtureBackend::calc_pmin_sat(long double &pmin_satL, long double &pmin_satV)
{
if (is_pure_or_pseudopure)
{
pmin_satL = components[0]->pEOS->sat_min_liquid.p;
pmin_satV = components[0]->pEOS->sat_min_vapor.p;
return;
}
else{
throw ValueError("calc_pmin_sat not yet defined for mixtures");
}
if (is_pure_or_pseudopure)
{
pmin_satL = components[0]->pEOS->sat_min_liquid.p;
pmin_satV = components[0]->pEOS->sat_min_vapor.p;
return;
}
else{
throw ValueError("calc_pmin_sat not yet defined for mixtures");
}
}
// Minimum allowed saturation temperature the maximum of the saturation temperatures of liquid and vapor
// For pure fluids, both values are the same, for pseudo-pure they are probably the same, for mixtures they are definitely not the same
// For pure fluids, both values are the same, for pseudo-pure they are probably the same, for mixtures they are definitely not the same
long double HelmholtzEOSMixtureBackend::calc_Tmax(void)
{
@@ -991,26 +991,26 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
this->_phase = iphase_liquid; _Q = 1000; return;
}
else if (value > s_liq + s_liq_error_band && value < s_vap - s_vap_error_band){ definitely_two_phase = true;}
break;
break;
}
case iUmolar:
{
if (!component.ancillaries.hL.enabled()){break;}
// u = h-p/rho
// u = h-p/rho
// Ancillaries are h-h_anchor, so add back h_anchor
// Ancillaries are h-h_anchor, so add back h_anchor
long double h_liq = component.ancillaries.hL.evaluate(_TLanc) + component.EOSVector[0].hs_anchor.hmolar;
long double h_liq_error_band = component.ancillaries.hL.get_max_abs_error();
long double h_vap = h_liq + component.ancillaries.hLV.evaluate(_TLanc);
long double h_vap_error_band = h_liq_error_band + component.ancillaries.hLV.get_max_abs_error();
long double rho_vap = component.ancillaries.rhoV.evaluate(_TVanc);
long double rho_liq = component.ancillaries.rhoL.evaluate(_TLanc);
long double u_liq = h_liq-_p/rho_liq;
long double u_vap = h_vap-_p/rho_vap;
long double u_liq_error_band = 1.5*h_liq_error_band; // Most of error is in enthalpy
long double u_vap_error_band = 1.5*h_vap_error_band; // Most of error is in enthalpy
// Check if in range given the accuracy of the fit
long double rho_vap = component.ancillaries.rhoV.evaluate(_TVanc);
long double rho_liq = component.ancillaries.rhoL.evaluate(_TLanc);
long double u_liq = h_liq-_p/rho_liq;
long double u_vap = h_vap-_p/rho_vap;
long double u_liq_error_band = 1.5*h_liq_error_band; // Most of error is in enthalpy
long double u_vap_error_band = 1.5*h_vap_error_band; // Most of error is in enthalpy
// Check if in range given the accuracy of the fit
if (value > u_vap + u_vap_error_band){
this->_phase = iphase_gas; _Q = -1000; return;
}
@@ -1018,7 +1018,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
this->_phase = iphase_liquid; _Q = 1000; return;
}
else if (value > u_liq + u_liq_error_band && value < u_vap - u_vap_error_band){ definitely_two_phase = true;}
break;
break;
}
default:
@@ -1907,16 +1907,16 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
}
catch(std::exception &)
{
try{
try{
// Next we try with a Brent method bounded solver since the function is 1-1
double rhomolar = Brent(resid, 0.1*rhomolar_guess, 2*rhomolar_guess,DBL_EPSILON,1e-8,100,errstring);
if (!ValidNumber(rhomolar)){throw ValueError();}
return rhomolar;
}
catch(...){
throw ValueError(format("solver_rho_Tp was unable to find a solution for T=%10Lg, p=%10Lg, with guess value %10Lg",T,p,rhomolar_guess));
}
}
catch(...){
throw ValueError(format("solver_rho_Tp was unable to find a solution for T=%10Lg, p=%10Lg, with guess value %10Lg",T,p,rhomolar_guess));
}
return _HUGE;
}
}
@@ -2113,27 +2113,27 @@ long double HelmholtzEOSMixtureBackend::calc_umolar_nocache(long double T, long
}
long double HelmholtzEOSMixtureBackend::calc_umolar(void)
{
if (isTwoPhase())
if (isTwoPhase())
{
_umolar = _Q*SatV->umolar() + (1 - _Q)*SatL->umolar();
return static_cast<long double>(_umolar);
}
else if (isHomogeneousPhase())
{
// Calculate the reducing parameters
_delta = _rhomolar/_reducing.rhomolar;
_tau = _reducing.T/_T;
// Calculate the reducing parameters
_delta = _rhomolar/_reducing.rhomolar;
_tau = _reducing.T/_T;
// Calculate derivatives if needed, or just use cached values
long double da0_dTau = dalpha0_dTau();
long double dar_dTau = dalphar_dTau();
long double R_u = gas_constant();
// Calculate derivatives if needed, or just use cached values
long double da0_dTau = dalpha0_dTau();
long double dar_dTau = dalphar_dTau();
long double R_u = gas_constant();
// Get molar internal energy
_umolar = R_u*_T*_tau.pt()*(da0_dTau+dar_dTau);
// Get molar internal energy
_umolar = R_u*_T*_tau.pt()*(da0_dTau+dar_dTau);
return static_cast<long double>(_umolar);
}
return static_cast<long double>(_umolar);
}
else{
throw ValueError(format("phase is invalid in calc_umolar"));
}

24
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
View File

@@ -54,10 +54,10 @@ public:
bool using_mole_fractions(){return true;}
bool using_mass_fractions(){return false;}
bool using_volu_fractions(){return false;}
bool has_melting_line(){ return is_pure_or_pseudopure && components[0]->ancillaries.melting_line.enabled();};
long double calc_melting_line(int param, int given, long double value);
phases calc_phase(void){return _phase;};
bool has_melting_line(){ return is_pure_or_pseudopure && components[0]->ancillaries.melting_line.enabled();};
long double calc_melting_line(int param, int given, long double value);
phases calc_phase(void){return _phase;};
void calc_specify_phase(phases phase){ specify_phase(phase); }
void calc_unspecify_phase(){ unspecify_phase(); }
long double calc_saturation_ancillary(parameters param, int Q, parameters given, double value);
@@ -145,11 +145,11 @@ public:
long double calc_cpmolar_idealgas(void);
long double calc_pressure_nocache(long double T, long double rhomolar);
long double calc_smolar(void);
long double calc_smolar_nocache(long double T, long double rhomolar);
long double calc_smolar_nocache(long double T, long double rhomolar);
long double calc_hmolar(void);
long double calc_hmolar_nocache(long double T, long double rhomolar);
long double calc_hmolar(void);
long double calc_hmolar_nocache(long double T, long double rhomolar);
long double calc_umolar_nocache(long double T, long double rhomolar);
long double calc_umolar(void);
long double calc_speed_sound(void);
@@ -191,13 +191,13 @@ public:
long double calc_conductivity_background(void);
long double calc_Tmin(void);
long double calc_Tmax(void);
long double calc_Tmax(void);
long double calc_pmax(void);
long double calc_Ttriple(void);
long double calc_p_triple(void);
long double calc_pmax_sat(void);
long double calc_Tmax_sat(void);
void calc_Tmin_sat(long double &Tmin_satL, long double &Tmin_satV);
long double calc_pmax_sat(void);
long double calc_Tmax_sat(void);
void calc_Tmin_sat(long double &Tmin_satL, long double &Tmin_satV);
void calc_pmin_sat(long double &pmin_satL, long double &pmin_satV);
long double calc_T_critical(void);

300
src/Backends/Helmholtz/MixtureDerivatives.h
View File

@@ -32,22 +32,22 @@ class MixtureDerivatives{
/** \brief GERG 2004 Monograph equation 7.62
*
* The derivative term
* \f[
* n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta})
* \f]
* \f[
* n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta})
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS);
*/
static long double ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS);
/** \brief GERG 2004 Monograph equation 7.61
*
* The derivative term
* \f[
* \left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau})
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS);
* \f[
* \left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau})
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
*/
static long double dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS);
static long double dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS);
@@ -56,31 +56,31 @@ class MixtureDerivatives{
static long double d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.63
/** \brief GERG 2004 Monograph equation 7.63
*
* The derivative term
* \f[
* n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* \f[
* n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 monograph Eqn. 7.32
*
* The partial molar volume
* \f[
* \hat v_i = \left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \frac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Fugacity of the i-th component
*
@@ -96,10 +96,10 @@ class MixtureDerivatives{
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
*/
static long double ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Derivative of the natural logarithm of the fugacity with respect to T
/** \brief Derivative of the natural logarithm of the fugacity with respect to T
*
* From Witzke, Eqn. 3.14
* \f[
@@ -108,10 +108,10 @@ class MixtureDerivatives{
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double dln_fugacity_i_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double dln_fugacity_i_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Derivative of the natural logarithm of the fugacity with respect to T
/** \brief Derivative of the natural logarithm of the fugacity with respect to T
*
* From Witzke, Eqn. 3.15
* \f[
@@ -120,21 +120,21 @@ class MixtureDerivatives{
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double dln_fugacity_i_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double dln_fugacity_i_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Eqn. 7.29
*
* The derivative term
* The derivative term
*
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n}
* \left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
static long double dln_fugacity_i_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
@@ -148,217 +148,217 @@ class MixtureDerivatives{
* The derivative term
* \f[
* n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j}
* \f]
* which is equal to
* \f{eqnarray*}{
* n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k}
* \f}
* \f]
* which is equal to
* \f{eqnarray*}{
* n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/// GERG Equation 7.42
static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/// GERG Equation 7.42
static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Eqn. 7.30
*
* The derivative term
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p}
* \f]
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.31
/** \brief GERG 2004 Monograph Equation 7.31
*
* The derivative term
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* which is also equal to
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right]
* \f]
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* which is also equal to
* \f[
* n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param j The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
*/
static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.115
/** \brief Gernert Equation 3.115
*
* The derivative term
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* \f[
* \left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The first index of interest
* @param j The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*
*/
static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
*/
static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.130
/** \brief Gernert Equation 3.130
*
* The derivative term
* \f[
* \f[
* \left(\frac{\partial p}{\partial x_j} \right)_{T,V,x_{k\neq j}} = \rho RT\left(-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_{k\neq j}} \delta\alpha_{\delta}^r + \delta\left(\frac{\partial}{\partial x_j}\left(\left( \frac{\partial \alpha^r}{\partial \delta}\right)_{\tau,\bar x}\right)\right)_{T,V,x_{k\neq j}}\right)
* \f]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param j The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Gernert Equation 3.117
/** \brief Gernert Equation 3.117
*
* The derivative term
* \f[
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial x_j\partial n_i} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}}
* \f]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param j The first index of interest
* @param i The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
*/
static long double d2nalphar_dxj_dni__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, std::size_t i, x_N_dependency_flag xN_flag){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j, xN_flag) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j, xN_flag);};
/// Gernert Equation 3.119
/// Catch test provided
static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.119
/// Catch test provided
static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.118
/// Catch test provided
static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.118
/// Catch test provided
static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.134
/// Catch test provided
static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.134
/// Catch test provided
static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.121
/// Catch test provided
static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.121
/// Catch test provided
static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.122
/// Catch test provided
static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/// Gernert Equation 3.122
/// Catch test provided
static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph, equations 7.44 and 7.51
/** \brief GERG 2004 Monograph, equations 7.44 and 7.51
*
* The derivative term
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n}
* \f]
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right]
* \f]
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n}
* \f]
* \f[
* \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right]
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012
/** \brief GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012
*
* The derivative term
* \f{eqnarray*}{
* \frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau}
* \f}
* \f{eqnarray*}{
* \frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\
* && +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\
* && +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012
/** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012
*
* The derivative term
* \f{eqnarray*}{
* \left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\
* &+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k}
* \f}
* \f{eqnarray*}{
* \left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\
* &+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k}
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.41
/** \brief GERG 2004 Monograph equation 7.41
*
* The derivative term
* \f[
* n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i}
* \f]
* and
* GERG 2004 Monograph equation 7.46:
* \f[
* n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i}
* \f]
* GERG 2004 Monograph equation 7.47:
* \f{eqnarray*}{
* n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\
* \f}
* \f[
* n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i}
* \f]
* and
* GERG 2004 Monograph equation 7.46:
* \f[
* n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i}
* \f]
* GERG 2004 Monograph equation 7.47:
* \f{eqnarray*}{
* n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\
* &+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param j The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
*/
static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.48
*
* The derivative term
* \f[
* n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}
* \f]
* \f[
* n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.49
*
* The derivative term
* \f[
* n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}
* \f]
* \f[
* n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}
* \f]
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
*/
static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.52
*
* The derivative term
* \f{eqnarray*}{
* \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r
* \f}
* The derivative term
* \f{eqnarray*}{
* \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\
* &+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\
* &+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r
* \f}
* @param HEOS The HelmholtzEOSMixtureBackend to be used
* @param i The index of interest
* @param j The second index of interest
* @param xN_flag A flag specifying whether the all mole fractions are independent or only the first N-1
*/
static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
*/
static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
}; /* class MixtureDerivatives */

6
src/Backends/Helmholtz/MixtureParameters.cpp
View File

@@ -278,9 +278,9 @@ void MixtureParameters::set_mixture_parameters(HelmholtzEOSMixtureBackend &HEOS)
STLMatrix beta_v, gamma_v, beta_T, gamma_T;
beta_v.resize(N, std::vector<long double>(N, 0));
gamma_v.resize(N, std::vector<long double>(N, 0));
beta_T.resize(N, std::vector<long double>(N, 0));
gamma_T.resize(N, std::vector<long double>(N, 0));
gamma_v.resize(N, std::vector<long double>(N, 0));
beta_T.resize(N, std::vector<long double>(N, 0));
gamma_T.resize(N, std::vector<long double>(N, 0));
HEOS.Excess.resize(N);

70
src/Backends/Helmholtz/ReducingFunctions.cpp
View File

@@ -29,12 +29,12 @@ long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<long doub
}
long double ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
long double s = 0;
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2rhormolardxidxj(x,j,k, xN_flag);
}
return d2rhormolardxidxj(x,j,i, xN_flag)-drhormolardxi__constxj(x,j, xN_flag)-s;
long double s = 0;
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2rhormolardxidxj(x,j,k, xN_flag);
}
return d2rhormolardxidxj(x,j,i, xN_flag)-drhormolardxi__constxj(x,j, xN_flag)-s;
}
long double ReducingFunction::ndrhorbardni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
@@ -88,23 +88,23 @@ long double GERG2008ReducingFunction::Tr(const std::vector<long double> &x)
}
long double GERG2008ReducingFunction::dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return dYrdxi__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
return dYrdxi__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return d2Yrdxi2__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
return d2Yrdxi2__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
return d2Yrdxidxj(x, i, j, beta_T, gamma_T, T_c, Yc_T, xN_flag);
return d2Yrdxidxj(x, i, j, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::rhormolar(const std::vector<long double> &x)
{
return 1/Yr(x, beta_v, gamma_v, v_c, Yc_v);
return 1/Yr(x, beta_v, gamma_v, v_c, Yc_v);
}
long double GERG2008ReducingFunction::drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return -pow(rhormolar(x),2)*dvrmolardxi__constxj(x, i, xN_flag);
return -pow(rhormolar(x),2)*dvrmolardxi__constxj(x, i, xN_flag);
}
long double GERG2008ReducingFunction::dvrmolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
@@ -120,21 +120,21 @@ long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<long dou
}
long double GERG2008ReducingFunction::d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
long double rhor = this->rhormolar(x);
long double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
return 2*pow(rhor,3)*pow(dvrbardxi,2)-pow(rhor,2)*this->d2vrmolardxi2__constxj(x,i, xN_flag);
long double rhor = this->rhormolar(x);
long double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
return 2*pow(rhor,3)*pow(dvrbardxi,2)-pow(rhor,2)*this->d2vrmolardxi2__constxj(x,i, xN_flag);
}
long double GERG2008ReducingFunction::d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
double rhor = this->rhormolar(x);
double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
double dvrbardxj = this->dvrmolardxi__constxj(x,j, xN_flag);
return 2*pow(rhor,3)*dvrbardxi*dvrbardxj-pow(rhor,2)*this->d2vrmolardxidxj(x,i,j, xN_flag);
double rhor = this->rhormolar(x);
double dvrbardxi = this->dvrmolardxi__constxj(x,i, xN_flag);
double dvrbardxj = this->dvrmolardxi__constxj(x,j, xN_flag);
return 2*pow(rhor,3)*dvrbardxi*dvrbardxj-pow(rhor,2)*this->d2vrmolardxidxj(x,i,j, xN_flag);
}
long double GERG2008ReducingFunction::Yr(const std::vector<long double> &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc)
{
long double Yr = 0;
long double Yr = 0;
for (std::size_t i = 0; i < N; i++)
{
double xi = x[i];
@@ -148,7 +148,7 @@ long double GERG2008ReducingFunction::Yr(const std::vector<long double> &x, cons
Yr += c_Y_ij(i, j, beta, gamma, Y_c_ij)*f_Y_ij(x, i, j, beta);
}
}
return Yr;
return Yr;
}
long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag)
{
@@ -275,39 +275,39 @@ long double GERG2008ReducingFunction::d2Yrdxidxj(const std::vector<long double>
long double GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta)
{
double xk = x[k], xi = x[i], beta_Y = beta[k][i];
return xk*(xk+xi)/(beta_Y*beta_Y*xk+xi)+xk*xi/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi));
double xk = x[k], xi = x[i], beta_Y = beta[k][i];
return xk*(xk+xi)/(beta_Y*beta_Y*xk+xi)+xk*xi/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi));
}
long double GERG2008ReducingFunction::dfYikdxi__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta)
{
double xk = x[k], xi = x[i], beta_Y = beta[i][k];
return xk*(xi+xk)/(beta_Y*beta_Y*xi+xk)+xi*xk/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk));
double xk = x[k], xi = x[i], beta_Y = beta[i][k];
return xk*(xi+xk)/(beta_Y*beta_Y*xi+xk)+xi*xk/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk));
}
long double GERG2008ReducingFunction::c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c)
{
return 2*beta[i][j]*gamma[i][j]*Y_c[i][j];
return 2*beta[i][j]*gamma[i][j]*Y_c[i][j];
}
long double GERG2008ReducingFunction::f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta)
{
double xi = x[i], xj = x[j], beta_Y = beta[i][j];
return xi*xj*(xi+xj)/(beta_Y*beta_Y*xi+xj);
double xi = x[i], xj = x[j], beta_Y = beta[i][j];
return xi*xj*(xi+xj)/(beta_Y*beta_Y*xi+xj);
}
long double GERG2008ReducingFunction::d2fYikdxi2__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta)
{
double xi = x[i], xk = x[k], beta_Y = beta[i][k];
return 1/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk))*(2*xk-xi*xk*2*beta_Y*beta_Y/(beta_Y*beta_Y*xi+xk));
double xi = x[i], xk = x[k], beta_Y = beta[i][k];
return 1/(beta_Y*beta_Y*xi+xk)*(1-beta_Y*beta_Y*(xi+xk)/(beta_Y*beta_Y*xi+xk))*(2*xk-xi*xk*2*beta_Y*beta_Y/(beta_Y*beta_Y*xi+xk));
}
long double GERG2008ReducingFunction::d2fYkidxi2__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta)
{
double xi = x[i], xk = x[k], beta_Y = beta[k][i];
return 1/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi))*(2*xk-xk*xi*2/(beta_Y*beta_Y*xk+xi));
double xi = x[i], xk = x[k], beta_Y = beta[k][i];
return 1/(beta_Y*beta_Y*xk+xi)*(1-(xk+xi)/(beta_Y*beta_Y*xk+xi))*(2*xk-xk*xi*2/(beta_Y*beta_Y*xk+xi));
}
long double GERG2008ReducingFunction::d2fYijdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta)
{
double xi = x[i], xj = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y*beta_Y;
return (xi+xj)/(beta_Y2*xi+xj) + xj/(beta_Y2*xi+xj)*(1-(xi+xj)/(beta_Y2*xi+xj))
+xi/(beta_Y2*xi+xj)*(1-beta_Y2*(xi+xj)/(beta_Y2*xi+xj))
-xi*xj/pow(beta_Y2*xi+xj,2)*(1+beta_Y2-2*beta_Y2*(xi+xj)/(beta_Y2*xi+xj));
double xi = x[i], xj = x[j], beta_Y = beta[i][j], beta_Y2 = beta_Y*beta_Y;
return (xi+xj)/(beta_Y2*xi+xj) + xj/(beta_Y2*xi+xj)*(1-(xi+xj)/(beta_Y2*xi+xj))
+xi/(beta_Y2*xi+xj)*(1-beta_Y2*(xi+xj)/(beta_Y2*xi+xj))
-xi*xj/pow(beta_Y2*xi+xj,2)*(1+beta_Y2-2*beta_Y2*(xi+xj)/(beta_Y2*xi+xj));
}
} /* namespace CoolProp */

178
src/Backends/Helmholtz/ReducingFunctions.h
View File

@@ -29,57 +29,57 @@ std::string get_reducing_function_name(std::string CAS1, std::string CAS2);
class ReducingFunction
{
protected:
std::size_t N;
std::size_t N;
public:
ReducingFunction(){};
virtual ~ReducingFunction(){};
ReducingFunction(){};
virtual ~ReducingFunction(){};
/// A factory function to generate the requiredreducing function
static shared_ptr<ReducingFunction> factory(const std::vector<CoolPropFluid*> &components, std::vector< std::vector< long double> > &F);
/// The reduced temperature
virtual long double Tr(const std::vector<long double> &x) = 0;
/// The derivative of reduced temperature with respect to component i mole fraction
virtual long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
/// The molar reducing density
virtual long double rhormolar(const std::vector<long double> &x) = 0;
///Derivative of the molar reducing density with respect to component i mole fraction
virtual long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
/// The reduced temperature
virtual long double Tr(const std::vector<long double> &x) = 0;
/// The derivative of reduced temperature with respect to component i mole fraction
virtual long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
/// The molar reducing density
virtual long double rhormolar(const std::vector<long double> &x) = 0;
///Derivative of the molar reducing density with respect to component i mole fraction
virtual long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag) = 0;
virtual long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) = 0;
/** \brief GERG 2004 Monograph equation 7.56:
/** \brief GERG 2004 Monograph equation 7.56:
*
* If the \f$x_i\f$ are all independent
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
* \f]
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
* \f]
* If \f$x_N = 1-\sum x_i\f$:
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
* \f]
*/
long double d_ndTrdni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2T_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2T_r}{\partial x_j \partial x_k}\right)
* \f]
*/
long double d_ndTrdni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief
*
* GERG 2004 Monograph equation 7.55:
* If the \f$x_i\f$ are all independent
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
* \f]
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-1}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
* \f]
* Gernert, JPCRD, 2014, A28
* If \f$x_N = 1-\sum x_i\f$:
* \f[
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-2}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
* \f]
*/
long double d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
* \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j}\right)\right)_{x_i} = \left(\frac{\partial^2\rho_r}{\partial x_j \partial x_i}\right)-\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}-\sum_{k=0}^{N-2}x_k\left(\frac{\partial^2\rho_r}{\partial x_j \partial x_k}\right)
* \f]
*/
long double d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double ndrhorbardni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double ndTrdni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double ndrhorbardni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
long double ndTrdni__constnj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
};
/** \brief The reducing function model of GERG-2008
@@ -92,53 +92,53 @@ class GERG2008ReducingFunction : public ReducingFunction
private:
GERG2008ReducingFunction(const GERG2008ReducingFunction& that); // No copying
protected:
STLMatrix v_c; ///< \f$ v_{c,ij} = \frac{1}{8}\left(v_{c,i}^{1/3}+v_{c,j}^{1/3}\right)^{3}\f$ from GERG-2008
STLMatrix T_c; ///< \f$ T_{c,ij} = \sqrt{T_{c,i}T_{c,j}} \f$ from GERG=2008
STLMatrix beta_v; ///< \f$ \beta_{v,ij} \f$ from GERG-2008
STLMatrix gamma_v; ///< \f$ \gamma_{v,ij} \f$ from GERG-2008
STLMatrix beta_T; ///< \f$ \beta_{T,ij} \f$ from GERG-2008
STLMatrix gamma_T; ///< \f$ \gamma_{T,ij} \f$ from GERG-2008
STLMatrix v_c; ///< \f$ v_{c,ij} = \frac{1}{8}\left(v_{c,i}^{1/3}+v_{c,j}^{1/3}\right)^{3}\f$ from GERG-2008
STLMatrix T_c; ///< \f$ T_{c,ij} = \sqrt{T_{c,i}T_{c,j}} \f$ from GERG=2008
STLMatrix beta_v; ///< \f$ \beta_{v,ij} \f$ from GERG-2008
STLMatrix gamma_v; ///< \f$ \gamma_{v,ij} \f$ from GERG-2008
STLMatrix beta_T; ///< \f$ \beta_{T,ij} \f$ from GERG-2008
STLMatrix gamma_T; ///< \f$ \gamma_{T,ij} \f$ from GERG-2008
std::vector<long double> Yc_T; ///< Vector of critical temperatures for all components
std::vector<long double> Yc_v; ///< Vector of critical molar volumes for all components
std::vector<CoolPropFluid *> pFluids; ///< List of pointer to fluids
std::vector<CoolPropFluid *> pFluids; ///< List of pointer to fluids
public:
GERG2008ReducingFunction(std::vector<CoolPropFluid *> pFluids, STLMatrix beta_v, STLMatrix gamma_v, STLMatrix beta_T, STLMatrix gamma_T)
{
this->pFluids = pFluids;
this->beta_v = beta_v;
this->gamma_v = gamma_v;
this->beta_T = beta_T;
this->gamma_T = gamma_T;
this->N = pFluids.size();
T_c.resize(N,std::vector<long double>(N,0));
v_c.resize(N,std::vector<long double>(N,0));
GERG2008ReducingFunction(std::vector<CoolPropFluid *> pFluids, STLMatrix beta_v, STLMatrix gamma_v, STLMatrix beta_T, STLMatrix gamma_T)
{
this->pFluids = pFluids;
this->beta_v = beta_v;
this->gamma_v = gamma_v;
this->beta_T = beta_T;
this->gamma_T = gamma_T;
this->N = pFluids.size();
T_c.resize(N,std::vector<long double>(N,0));
v_c.resize(N,std::vector<long double>(N,0));
Yc_T.resize(N);
Yc_v.resize(N);
for (std::size_t i = 0; i < N; ++i)
{
for (std::size_t j = 0; j < N; j++)
{
T_c[i][j] = sqrt(pFluids[i]->pEOS->reduce.T*pFluids[j]->pEOS->reduce.T);
v_c[i][j] = 1.0/8.0*pow(pow(pFluids[i]->pEOS->reduce.rhomolar, -1.0/3.0)+pow(pFluids[j]->pEOS->reduce.rhomolar, -1.0/3.0),3);
}
for (std::size_t i = 0; i < N; ++i)
{
for (std::size_t j = 0; j < N; j++)
{
T_c[i][j] = sqrt(pFluids[i]->pEOS->reduce.T*pFluids[j]->pEOS->reduce.T);
v_c[i][j] = 1.0/8.0*pow(pow(pFluids[i]->pEOS->reduce.rhomolar, -1.0/3.0)+pow(pFluids[j]->pEOS->reduce.rhomolar, -1.0/3.0),3);
}
Yc_T[i] = pFluids[i]->pEOS->reduce.T;
Yc_v[i] = 1/pFluids[i]->pEOS->reduce.rhomolar;
}
};
}
};
/// Default destructor
~GERG2008ReducingFunction(){};
/** \brief The reducing temperature
/// Default destructor
~GERG2008ReducingFunction(){};
/** \brief The reducing temperature
* Calculated from \ref Yr with \f$T = Y\f$
*/
long double Tr(const std::vector<long double> &x);
/** \brief The derivative of reducing temperature with respect to component i mole fraction
long double Tr(const std::vector<long double> &x);
/** \brief The derivative of reducing temperature with respect to component i mole fraction
*
* Calculated from \ref dYrdxi__constxj with \f$T = Y\f$
*/
long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief The second derivative of reducing temperature with respect to component i mole fraction
long double dTrdxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief The second derivative of reducing temperature with respect to component i mole fraction
*
* Calculated from \ref d2Yrdxi2__constxj with \f$T = Y\f$
*/
@@ -147,7 +147,7 @@ public:
*
* Calculated from \ref d2Yrdxidxj with \f$T = Y\f$
*/
long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief The derivative of reducing molar volume with respect to component i mole fraction
*
@@ -168,7 +168,7 @@ public:
*
* Given by \f$ \rho_r = 1/v_r \f$
*/
long double rhormolar(const std::vector<long double> &x);
long double rhormolar(const std::vector<long double> &x);
/** \brief Derivative of the molar reducing density with respect to component i mole fraction
*
* See also GERG 2004, Eqn. 7.57
@@ -176,8 +176,8 @@ public:
* \left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_{i\neq j}} = -\rho_r^2\left(\frac{\partial v_r}{\partial x_i}\right)_{x_{i\neq j}}
* \f]
*/
long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Derivative of the molar reducing density with respect to component i mole fraction
long double drhormolardxi__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag);
/** \brief Derivative of the molar reducing density with respect to component i mole fraction
*
* See also GERG 2004, Eqn. 7.58
* \f[
@@ -192,8 +192,8 @@ public:
* \left(\frac{\partial^2 \rho_r}{\partial x_i\partial x_j}\right) = 2\rho_r^3\left(\left(\frac{\partial v_r}{\partial x_i}\right)_{x_{i\neq j}}\right)\left(\left(\frac{\partial v_r}{\partial x_j}\right)_{x_{i\neq j}}\right)-\rho_r^2\left(\left(\frac{\partial v_r}{\partial x_i\partial x_j}\right)\right)
* \f]
*/
long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
long double d2rhormolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief Generalized reducing term \f$Y_r\f$
*
* \f[
@@ -253,28 +253,28 @@ public:
* c_{Y,ij} = 2\beta_{Y,ij}\gamma_{Y,ij}Y_{c,ij}
* \f]
*/
long double c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
long double c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
/** \brief The function \f$ f_{Y,ij}(x_i,x_j) \f$
*
* \f[ f_{Y,ij}(x_i,x_j) = x_ix_j\frac{x_i+x_j}{\beta_{Y,ij}^2x_i+x_j} \f]
*/
long double f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta);
long double f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta);
/**
*
* \f[
* \left(\frac{\partial f_{Y,ki}(x_k, x_i)}{\partial x_i}\right)_{x_{k\neq i}} = x_k\frac{x_k+x_i}{\beta_{Y,ki}^2x_k+x_i} + \frac{x_kx_i}{\beta_{Y,ki}^2x_k+x_i}\left(1-\frac{x_k+x_i}{\beta_{Y,ki}^2x_k+x_i}\right)
* \f]
*/
long double dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta);
long double dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta);
/**
*
* \f[
* \left(\frac{\partial f_{Y,ik}(x_i, x_k)}{\partial x_i}\right)_{x_k} = x_k\frac{x_i+x_k}{\beta_{Y,ik}^2x_i+x_k} + \frac{x_ix_k}{\beta_{Y,ik}^2x_i+x_k}\left(1-\beta_{Y,ik}^2\frac{x_i+x_k}{\beta_{Y,ik}^2x_i+x_k}\right)
* \f]
*/
long double dfYikdxi__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
/**
long double dfYikdxi__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
/**
* \f[
* \left(\frac{\partial^2 f_{Y,ki}(x_k, x_i)}{\partial x_i^2}\right)_{x_{k\neq i}} = \frac{1}{\beta_{Y,ki}^2x_k+x_i}\left(1-\frac{x_k+x_i}{\beta_{Y,ki}^2x_k+x_i}\right)\left(2x_k-\frac{2x_kx_i}{\beta_{Y,ki}^2x_k+x_i}\right)
* \f]
@@ -285,8 +285,8 @@ public:
* \left(\frac{\partial^2 f_{Y,ik}(x_i, x_k)}{\partial x_i^2}\right)_{x_{k}} = \frac{1}{\beta_{Y,ik}^2x_i+x_k}\left(1-\beta_{Y,ik}^2\frac{x_i+x_k}{\beta_{Y,ik}^2x_i+x_k}\right)\left(2x_k-\frac{2x_ix_k\beta_{Y,ik}^2}{\beta_{Y,ik}^2x_i+x_k}\right)
* \f]
*/
long double d2fYikdxi2__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
/**
long double d2fYikdxi2__constxk(const std::vector<long double> &x, std::size_t i, std::size_t k, const STLMatrix &beta);
/**
* \f{eqnarray*}{
* \left(\frac{\partial^2 f_{Y,ki}(x_k, x_i)}{\partial x_i\partial x_j}\right)_{x_{k\neq j\neq i}} &=& \frac{x_i+x_j}{\beta_{Y,ij}^2x_i+x_j} + \frac{x_j}{\beta_{Y,ij}^2x_i+x_j}\left(1-\frac{x_i+x_j}{\beta_{Y,ij}^2x_i+x_j}\right) \\
* &+& \frac{x_i}{\beta_{Y,ij}^2x_i+x_j}\left(1-\beta_{Y,ij}^2\frac{x_i+x_j}{\beta_{Y,ij}^2x_i+x_j}\right) - \frac{x_ix_j}{(\beta_{Y,ij}^2x_i+x_j)^2}\left(1+\beta_{Y,ij}^2-2\beta_{Y,ij}^2\frac{x_i+x_j}{\beta_{Y,ij}^2x_i+x_j}\right)
@@ -323,18 +323,18 @@ class LemmonAirHFCReducingFunction
protected:
LemmonAirHFCReducingFunction(const LemmonAirHFCReducingFunction &);
public:
/// Set the coefficients based on reducing parameters loaded from JSON
static void convert_to_GERG(const std::vector<CoolPropFluid*> &pFluids, std::size_t i, std::size_t j, Dictionary d, long double &beta_T, long double &beta_v, long double &gamma_T, long double &gamma_v)
/// Set the coefficients based on reducing parameters loaded from JSON
static void convert_to_GERG(const std::vector<CoolPropFluid*> &pFluids, std::size_t i, std::size_t j, Dictionary d, long double &beta_T, long double &beta_v, long double &gamma_T, long double &gamma_v)
{
long double xi_ij = d.get_number("xi");
long double zeta_ij = d.get_number("zeta");
beta_T = 1;
beta_v = 1;
gamma_T = (pFluids[i]->pEOS->reduce.T + pFluids[j]->pEOS->reduce.T + xi_ij)/(2*sqrt(pFluids[i]->pEOS->reduce.T*pFluids[j]->pEOS->reduce.T));
long double v_i = 1/pFluids[i]->pEOS->reduce.rhomolar;
long double v_j = 1/pFluids[j]->pEOS->reduce.rhomolar;
long double one_third = 1.0/3.0;
gamma_v = (v_i + v_j + zeta_ij)/(0.25*pow(pow(v_i, one_third)+pow(v_j, one_third),3));
beta_T = 1;
beta_v = 1;
gamma_T = (pFluids[i]->pEOS->reduce.T + pFluids[j]->pEOS->reduce.T + xi_ij)/(2*sqrt(pFluids[i]->pEOS->reduce.T*pFluids[j]->pEOS->reduce.T));
long double v_i = 1/pFluids[i]->pEOS->reduce.rhomolar;
long double v_j = 1/pFluids[j]->pEOS->reduce.rhomolar;
long double one_third = 1.0/3.0;
gamma_v = (v_i + v_j + zeta_ij)/(0.25*pow(pow(v_i, one_third)+pow(v_j, one_third),3));
};
};

622
src/Backends/Helmholtz/TransportRoutines.cpp
View File

@@ -187,113 +187,113 @@ long double TransportRoutines::viscosity_initial_density_dependence_empirical(He
static void visc_Helper(double Tbar, double rhobar, double *mubar_0, double *mubar_1)
{
std::vector<std::vector<long double> > H(6,std::vector<long double>(7,0));
std::vector<std::vector<long double> > H(6,std::vector<long double>(7,0));
double sum;
int i,j;
int i,j;
// Dilute-gas component
*mubar_0=100.0*sqrt(Tbar)/(1.67752+2.20462/Tbar+0.6366564/powInt(Tbar,2)-0.241605/powInt(Tbar,3));
// Dilute-gas component
*mubar_0=100.0*sqrt(Tbar)/(1.67752+2.20462/Tbar+0.6366564/powInt(Tbar,2)-0.241605/powInt(Tbar,3));
//Fill in zeros in H
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
H[i][j]=0;
}
}
//Fill in zeros in H
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
H[i][j]=0;
}
}
//Set non-zero parameters of H
H[0][0]=5.20094e-1;
H[1][0]=8.50895e-2;
H[2][0]=-1.08374;
H[3][0]=-2.89555e-1;
//Set non-zero parameters of H
H[0][0]=5.20094e-1;
H[1][0]=8.50895e-2;
H[2][0]=-1.08374;
H[3][0]=-2.89555e-1;
H[0][1]=2.22531e-1;
H[1][1]=9.99115e-1;
H[2][1]=1.88797;
H[3][1]=1.26613;
H[5][1]=1.20573e-1;
H[0][1]=2.22531e-1;
H[1][1]=9.99115e-1;
H[2][1]=1.88797;
H[3][1]=1.26613;
H[5][1]=1.20573e-1;
H[0][2]=-2.81378e-1;
H[1][2]=-9.06851e-1;
H[2][2]=-7.72479e-1;
H[3][2]=-4.89837e-1;
H[4][2]=-2.57040e-1;
H[0][2]=-2.81378e-1;
H[1][2]=-9.06851e-1;
H[2][2]=-7.72479e-1;
H[3][2]=-4.89837e-1;
H[4][2]=-2.57040e-1;
H[0][3]=1.61913e-1;
H[1][3]=2.57399e-1;
H[0][3]=1.61913e-1;
H[1][3]=2.57399e-1;
H[0][4]=-3.25372e-2;
H[3][4]=6.98452e-2;
H[0][4]=-3.25372e-2;
H[3][4]=6.98452e-2;
H[4][5]=8.72102e-3;
H[4][5]=8.72102e-3;
H[3][6]=-4.35673e-3;
H[5][6]=-5.93264e-4;
H[3][6]=-4.35673e-3;
H[5][6]=-5.93264e-4;
// Finite density component
sum=0;
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
sum+=powInt(1/Tbar-1,i)*(H[i][j]*powInt(rhobar-1,j));
}
}
*mubar_1=exp(rhobar*sum);
// Finite density component
sum=0;
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
sum+=powInt(1/Tbar-1,i)*(H[i][j]*powInt(rhobar-1,j));
}
}
*mubar_1=exp(rhobar*sum);
}
long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double x_mu=0.068,qc=1/1.9,qd=1/1.1,nu=0.630,gamma=1.239,zeta_0=0.13,LAMBDA_0=0.06,Tbar_R=1.5, pstar, Tstar, rhostar;
double delta,tau,mubar_0,mubar_1,mubar_2,drhodp,drhodp_R,DeltaChibar,zeta,w,L,Y,psi_D,Tbar,rhobar;
double drhobar_dpbar,drhobar_dpbar_R,R_Water;
double x_mu=0.068,qc=1/1.9,qd=1/1.1,nu=0.630,gamma=1.239,zeta_0=0.13,LAMBDA_0=0.06,Tbar_R=1.5, pstar, Tstar, rhostar;
double delta,tau,mubar_0,mubar_1,mubar_2,drhodp,drhodp_R,DeltaChibar,zeta,w,L,Y,psi_D,Tbar,rhobar;
double drhobar_dpbar,drhobar_dpbar_R,R_Water;
pstar = 22.064e6; // [Pa]
Tstar = 647.096; // [K]
rhostar = 322; // [kg/m^3]
Tbar = HEOS.T()/Tstar;
Tbar = HEOS.T()/Tstar;
rhobar = HEOS.rhomass()/rhostar;
R_Water = HEOS.gas_constant()/HEOS.molar_mass(); // [J/kg/K]
// Dilute and finite gas portions
visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
// Dilute and finite gas portions
visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
// **********************************************************************
// ************************ Critical Enhancement ************************
// ************************ Critical Enhancement ************************
// **********************************************************************
delta=rhobar;
// "Normal" calculation
tau=1/Tbar;
drhodp=1/(R_Water*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2()));
drhobar_dpbar = pstar/rhostar*drhodp;
// "Reducing" calculation
tau=1/Tbar_R;
drhodp_R=1/(R_Water*Tbar_R*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta)));
drhobar_dpbar_R = pstar/rhostar*drhodp_R;
delta=rhobar;
// "Normal" calculation
tau=1/Tbar;
drhodp=1/(R_Water*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2()));
drhobar_dpbar = pstar/rhostar*drhodp;
// "Reducing" calculation
tau=1/Tbar_R;
drhodp_R=1/(R_Water*Tbar_R*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta)));
drhobar_dpbar_R = pstar/rhostar*drhodp_R;
DeltaChibar=rhobar*(drhobar_dpbar-drhobar_dpbar_R*Tbar_R/Tbar);
if (DeltaChibar<0)
DeltaChibar=0;
zeta=zeta_0*pow(DeltaChibar/LAMBDA_0,nu/gamma);
if (zeta<0.3817016416){
Y=1.0/5.0*qc*zeta*powInt(qd*zeta,5)*(1-qc*zeta+powInt(qc*zeta,2)-765.0/504.0*powInt(qd*zeta,2));
}
else
{
psi_D=acos(pow(1+powInt(qd*zeta,2),-1.0/2.0));
w=sqrt(std::abs((qc*zeta-1)/(qc*zeta+1)))*tan(psi_D/2.0);
if (qc*zeta>1){
L=log((1+w)/(1-w));
}
else{
L=2*atan(std::abs(w));
}
Y=1.0/12.0*sin(3*psi_D)-1/(4*qc*zeta)*sin(2*psi_D)+1.0/powInt(qc*zeta,2)*(1-5.0/4.0*powInt(qc*zeta,2))*sin(psi_D)-1.0/powInt(qc*zeta,3)*((1-3.0/2.0*powInt(qc*zeta,2))*psi_D-pow(std::abs(powInt(qc*zeta,2)-1),3.0/2.0)*L);
}
mubar_2=exp(x_mu*Y);
DeltaChibar=rhobar*(drhobar_dpbar-drhobar_dpbar_R*Tbar_R/Tbar);
if (DeltaChibar<0)
DeltaChibar=0;
zeta=zeta_0*pow(DeltaChibar/LAMBDA_0,nu/gamma);
if (zeta<0.3817016416){
Y=1.0/5.0*qc*zeta*powInt(qd*zeta,5)*(1-qc*zeta+powInt(qc*zeta,2)-765.0/504.0*powInt(qd*zeta,2));
}
else
{
psi_D=acos(pow(1+powInt(qd*zeta,2),-1.0/2.0));
w=sqrt(std::abs((qc*zeta-1)/(qc*zeta+1)))*tan(psi_D/2.0);
if (qc*zeta>1){
L=log((1+w)/(1-w));
}
else{
L=2*atan(std::abs(w));
}
Y=1.0/12.0*sin(3*psi_D)-1/(4*qc*zeta)*sin(2*psi_D)+1.0/powInt(qc*zeta,2)*(1-5.0/4.0*powInt(qc*zeta,2))*sin(psi_D)-1.0/powInt(qc*zeta,3)*((1-3.0/2.0*powInt(qc*zeta,2))*psi_D-pow(std::abs(powInt(qc*zeta,2)-1),3.0/2.0)*L);
}
mubar_2=exp(x_mu*Y);
return (mubar_0*mubar_1*mubar_2)/1e6;
return (mubar_0*mubar_1*mubar_2)/1e6;
}
long double TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
@@ -314,7 +314,7 @@ long double TransportRoutines::viscosity_hexane_higher_order_hardcoded(Helmholtz
// Output is in Pa-s
double c[] = {2.53402335/1e6, -9.724061002/1e6, 0.469437316, 158.5571631, 72.42916856/1e6, 10.60751253, 8.628373915, -6.61346441, -2.212724566};
return pow(rhor,static_cast<long double>(2.0/3.0))*sqrt(Tr)*(c[0]/Tr+c[1]/(c[2]+Tr+c[3]*rhor*rhor)+c[4]*(1+rhor)/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
return pow(rhor,static_cast<long double>(2.0/3.0))*sqrt(Tr)*(c[0]/Tr+c[1]/(c[2]+Tr+c[3]*rhor*rhor)+c[4]*(1+rhor)/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
}
long double TransportRoutines::viscosity_heptane_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
@@ -324,7 +324,7 @@ long double TransportRoutines::viscosity_heptane_higher_order_hardcoded(Helmholt
// Output is in Pa-s
double c[] = {0, 22.15000/1e6, -15.00870/1e6, 3.71791/1e6, 77.72818/1e6, 9.73449, 9.51900, -6.34076, -2.51909};
return pow(rhor,2.0L/3.0L)*sqrt(Tr)*(c[1]*rhor+c[2]*pow(rhor,2)+c[3]*pow(rhor,3)+c[4]*rhor/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
return pow(rhor,2.0L/3.0L)*sqrt(Tr)*(c[1]*rhor+c[2]*pow(rhor,2)+c[3]*pow(rhor,3)+c[4]*rhor/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
}
long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzEOSMixtureBackend &HEOS)
@@ -336,26 +336,26 @@ long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzE
long double tau = F.T_reduce/HEOS.T(), kii, krrr, kaaa, krr, kdrdr;
double psi1 = exp(tau)-F.c1;
double psi2 = exp(pow(tau,2))-F.c2;
double psi2 = exp(pow(tau,2))-F.c2;
double ki = (F.Ai[0] + F.Ai[1]*psi1 + F.Ai[2]*psi2)*tau;
double ki = (F.Ai[0] + F.Ai[1]*psi1 + F.Ai[2]*psi2)*tau;
double ka = (F.Aa[0] + F.Aa[1]*psi1 + F.Aa[2]*psi2)*pow(tau, F.Na);
double kr = (F.Ar[0] + F.Ar[1]*psi1 + F.Ar[2]*psi2)*pow(tau, F.Nr);
double kaa = (F.Aaa[0] + F.Aaa[1]*psi1 + F.Aaa[2]*psi2)*pow(tau, F.Naa);
double kr = (F.Ar[0] + F.Ar[1]*psi1 + F.Ar[2]*psi2)*pow(tau, F.Nr);
double kaa = (F.Aaa[0] + F.Aaa[1]*psi1 + F.Aaa[2]*psi2)*pow(tau, F.Naa);
if (F.Arr.empty()){
krr = 0;
kdrdr = (F.Adrdr[0] + F.Adrdr[1]*psi1 + F.Adrdr[2]*psi2)*pow(tau, F.Nrr);
}
else{
krr = (F.Arr[0] + F.Arr[1]*psi1 + F.Arr[2]*psi2)*pow(tau, F.Nrr);
krr = (F.Arr[0] + F.Arr[1]*psi1 + F.Arr[2]*psi2)*pow(tau, F.Nrr);
kdrdr = 0;
}
if (!F.Aii.empty() && !F.Arrr.empty() && !F.Aaaa.empty()){
kii = (F.Aii[0] + F.Aii[1]*psi1 + F.Aii[2]*psi2)*pow(tau, F.Nii);
krrr = (F.Arrr[0] + F.Arrr[1]*psi1 + F.Arrr[2]*psi2)*pow(tau, F.Nrrr);
kaaa = (F.Aaaa[0] + F.Aaaa[1]*psi1 + F.Aaaa[2]*psi2)*pow(tau, F.Naaa);
kii = (F.Aii[0] + F.Aii[1]*psi1 + F.Aii[2]*psi2)*pow(tau, F.Nii);
krrr = (F.Arrr[0] + F.Arrr[1]*psi1 + F.Arrr[2]*psi2)*pow(tau, F.Nrrr);
kaaa = (F.Aaaa[0] + F.Aaaa[1]*psi1 + F.Aaaa[2]*psi2)*pow(tau, F.Naaa);
}
else{
kii = 0;
@@ -364,12 +364,12 @@ long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzE
}
double p = HEOS.p()/1e5; // [bar]; 1e5 for conversion from Pa -> bar
double pr = HEOS.T()*HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar)/1e5; // [bar/K]; 1e5 for conversion from Pa -> bar
double pa = p - pr; //[bar]
double pr = HEOS.T()*HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar)/1e5; // [bar/K]; 1e5 for conversion from Pa -> bar
double pa = p - pr; //[bar]
double pid = HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
double deltapr = pr - pid;
double deltapr = pr - pid;
double eta_f = ka*pa + kr*deltapr + ki*pid + kaa*pa*pa + kdrdr*deltapr*deltapr + krr*pr*pr + kii*pid*pid + krrr*pr*pr*pr + kaaa*pa*pa*pa;
double eta_f = ka*pa + kr*deltapr + ki*pid + kaa*pa*pa + kdrdr*deltapr*deltapr + krr*pr*pr + kii*pid*pid + krrr*pr*pr*pr + kaaa*pa*pa*pa;
return eta_f; //[Pa-s]
}
@@ -381,43 +381,43 @@ long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzE
long double TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double eta_0,eta_0_slash, eta_E_slash, B,C,D,ln_eta,x;
//
// Arp, V.D., McCarty, R.D., and Friend, D.G.,
// "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
// NIST Technical Note 1334 (revised), 1998.
//
// Using Arp NIST report
// Report is not clear on viscosity, referring to REFPROP source code for clarity
//
// Arp, V.D., McCarty, R.D., and Friend, D.G.,
// "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
// NIST Technical Note 1334 (revised), 1998.
//
// Using Arp NIST report
// Report is not clear on viscosity, referring to REFPROP source code for clarity
// Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
// Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
long double rho = HEOS.keyed_output(CoolProp::iDmass)/1000.0, T = HEOS.T();
if (T <= 300){
x = log(T);
}
else{
x = log(300.0);
}
// Evaluate the terms B,C,D
B = -47.5295259/x+87.6799309-42.0741589*x+8.33128289*x*x-0.589252385*x*x*x;
C = 547.309267/x-904.870586+431.404928*x-81.4504854*x*x+5.37008433*x*x*x;
D = -1684.39324/x+3331.08630-1632.19172*x+308.804413*x*x-20.2936367*x*x*x;
eta_0_slash = -0.135311743/x+1.00347841+1.20654649*x-0.149564551*x*x+0.012520841*x*x*x;
eta_E_slash = rho*B+rho*rho*C+rho*rho*rho*D;
x = log(T);
}
else{
x = log(300.0);
}
// Evaluate the terms B,C,D
B = -47.5295259/x+87.6799309-42.0741589*x+8.33128289*x*x-0.589252385*x*x*x;
C = 547.309267/x-904.870586+431.404928*x-81.4504854*x*x+5.37008433*x*x*x;
D = -1684.39324/x+3331.08630-1632.19172*x+308.804413*x*x-20.2936367*x*x*x;
eta_0_slash = -0.135311743/x+1.00347841+1.20654649*x-0.149564551*x*x+0.012520841*x*x*x;
eta_E_slash = rho*B+rho*rho*C+rho*rho*rho*D;
if (T<=100)
{
ln_eta = eta_0_slash + eta_E_slash;
{
ln_eta = eta_0_slash + eta_E_slash;
// Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
return exp(ln_eta)/10.0/1e6;
}
else
{
ln_eta = eta_0_slash + eta_E_slash;
eta_0 = 196*pow(T,static_cast<long double>(0.71938))*exp(12.451/T-295.67/T/T-4.1249);
return exp(ln_eta)/10.0/1e6;
}
else
{
ln_eta = eta_0_slash + eta_E_slash;
eta_0 = 196*pow(T,static_cast<long double>(0.71938))*exp(12.451/T-295.67/T/T-4.1249);
// Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
return (exp(ln_eta)+eta_0-exp(eta_0_slash))/10.0/1e6;
}
return (exp(ln_eta)+eta_0-exp(eta_0_slash))/10.0/1e6;
}
}
long double TransportRoutines::viscosity_methanol_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
@@ -484,68 +484,68 @@ long double TransportRoutines::viscosity_methanol_hardcoded(HelmholtzEOSMixtureB
long double TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double C1 = 1.3163, //
C2 = 0.1832,
DeltaGstar = 771.23,
rhoL = 32.174,
rhocbar = 7.5114,
C2 = 0.1832,
DeltaGstar = 771.23,
rhoL = 32.174,
rhocbar = 7.5114,
Tc = 299.2793,
DELTAeta_max = 3.967,
Ru = 8.31451,
DELTAeta_max = 3.967,
Ru = 8.31451,
molar_mass = 70.014;
double a[] = {0.4425728, -0.5138403, 0.1547566, -0.02821844, 0.001578286};
double e_k = 243.91, sigma = 0.4278;
double Tstar = HEOS.T()/e_k;
double logTstar = log(Tstar);
double Omega = exp(a[0]+a[1]*logTstar+a[2]*pow(logTstar,2)+a[3]*pow(logTstar,3)+a[4]*pow(logTstar,4));
double a[] = {0.4425728, -0.5138403, 0.1547566, -0.02821844, 0.001578286};
double e_k = 243.91, sigma = 0.4278;
double Tstar = HEOS.T()/e_k;
double logTstar = log(Tstar);
double Omega = exp(a[0]+a[1]*logTstar+a[2]*pow(logTstar,2)+a[3]*pow(logTstar,3)+a[4]*pow(logTstar,4));
double eta_DG = 1.25*0.021357*sqrt(molar_mass*HEOS.T())/(sigma*sigma*Omega); // uPa-s
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double eta_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(HEOS.T())*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*HEOS.T()));
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double eta_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(HEOS.T())*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*HEOS.T()));
double chi = rhobar - rhocbar;
double tau = HEOS.T() - Tc;
double chi = rhobar - rhocbar;
double tau = HEOS.T() - Tc;
double DELTAeta_c = 4*DELTAeta_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
double DELTAeta_c = 4*DELTAeta_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
return (pow((rhoL-rhobar)/rhoL,C1)*eta_DG+pow(rhobar/rhoL,C1)*eta_L+DELTAeta_c)/1e6;
return (pow((rhoL-rhobar)/rhoL,C1)*eta_DG+pow(rhobar/rhoL,C1)*eta_L+DELTAeta_c)/1e6;
}
long double TransportRoutines::viscosity_dilute_ethane(HelmholtzEOSMixtureBackend &HEOS)
{
double C[] = {0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
double OMEGA_2_2 = 0, e_k = 245, Tstar;
double OMEGA_2_2 = 0, e_k = 245, Tstar;
Tstar = HEOS.T()/e_k;
for (int i = 1; i<= 9; i++)
{
OMEGA_2_2 += C[i]*pow(Tstar,(i-1)/3.0-1);
}
Tstar = HEOS.T()/e_k;
for (int i = 1; i<= 9; i++)
{
OMEGA_2_2 += C[i]*pow(Tstar,(i-1)/3.0-1);
}
return 12.0085*sqrt(Tstar)*OMEGA_2_2/1e6; //[Pa-s]
return 12.0085*sqrt(Tstar)*OMEGA_2_2/1e6; //[Pa-s]
}
long double TransportRoutines::viscosity_dilute_cyclohexane(HelmholtzEOSMixtureBackend &HEOS)
{
// From Tariq, JPCRD, 2014
long double T = HEOS.T();
long double S_eta = exp(-1.5093 + 364.87/T - 39537/pow(T, 2)); //[nm^2]
return 0.19592*sqrt(T)/S_eta/1e6; //[Pa-s]
return 0.19592*sqrt(T)/S_eta/1e6; //[Pa-s]
}
long double TransportRoutines::viscosity_ethane_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double r[] = {0,1,1,2,2,2,3,3,4,4,1,1};
double s[] = {0,0,1,0,1,1.5,0,2,0,1,0,1};
double g[] = {0, 0.47177003, -0.23950311, 0.39808301, -0.27343335, 0.35192260, -0.21101308, -0.00478579, 0.07378129, -0.030435255, -0.30435286, 0.001215675};
double r[] = {0,1,1,2,2,2,3,3,4,4,1,1};
double s[] = {0,0,1,0,1,1.5,0,2,0,1,0,1};
double g[] = {0, 0.47177003, -0.23950311, 0.39808301, -0.27343335, 0.35192260, -0.21101308, -0.00478579, 0.07378129, -0.030435255, -0.30435286, 0.001215675};
double sum1 = 0, sum2 = 0, tau = 305.33/HEOS.T(), delta = HEOS.rhomolar()/6870;
double sum1 = 0, sum2 = 0, tau = 305.33/HEOS.T(), delta = HEOS.rhomolar()/6870;
for (int i = 1; i<= 9; ++i){
sum1 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
for (int i = 10; i<= 11; ++i){
sum2 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
return 15.977*sum1/(1+sum2)/1e6;
for (int i = 1; i<= 9; ++i){
sum1 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
for (int i = 10; i<= 11; ++i){
sum2 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
return 15.977*sum1/(1+sum2)/1e6;
}
long double TransportRoutines::conductivity_dilute_ratio_polynomials(HelmholtzEOSMixtureBackend &HEOS){
@@ -615,9 +615,9 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
// Retrieve values from the state class
CoolProp::ConductivityCriticalSimplifiedOlchowySengersData &data = HEOS.components[0]->transport.conductivity_critical.Olchowy_Sengers;
double k = data.k,
R0 = data.R0,
nu = data.nu,
double k = data.k,
R0 = data.R0,
nu = data.nu,
gamma = data.gamma,
GAMMA = data.GAMMA,
zeta0 = data.zeta0,
@@ -625,8 +625,8 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
Tc = HEOS.get_reducing_state().T, // [K]
rhoc = HEOS.get_reducing_state().rhomolar, // [mol/m^3]
Pcrit = HEOS.get_reducing_state().p, // [Pa]
Tref, // [K]
cp,cv,delta,num,zeta,mu,pi=M_PI,
Tref, // [K]
cp,cv,delta,num,zeta,mu,pi=M_PI,
OMEGA_tilde,OMEGA_tilde0;
if (ValidNumber(data.T_ref))
@@ -634,32 +634,32 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
else
Tref = 1.5*Tc;
delta = HEOS.delta();
delta = HEOS.delta();
double dp_drho = HEOS.gas_constant()*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2());
double X = Pcrit/pow(rhoc,2)*HEOS.rhomolar()/dp_drho;
double X = Pcrit/pow(rhoc,2)*HEOS.rhomolar()/dp_drho;
double tau_ref = Tc/Tref;
double tau_ref = Tc/Tref;
double dp_drho_ref = HEOS.gas_constant()*Tref*(1+2*delta*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau_ref,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau_ref,delta));
double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T();
num = X - Xref;
double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T();
num = X - Xref;
// No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
// No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
// See also Lemmon, IJT, 2004, page 27
if (num < DBL_EPSILON*10)
return 0.0;
else
zeta = zeta0*pow(num/GAMMA, nu/gamma); //[m]
if (num < DBL_EPSILON*10)
return 0.0;
else
zeta = zeta0*pow(num/GAMMA, nu/gamma); //[m]
cp = HEOS.cpmolar(); //[J/mol/K]
cv = HEOS.cvmolar(); //[J/mol/K]
mu = HEOS.viscosity(); //[Pa-s]
cv = HEOS.cvmolar(); //[J/mol/K]
mu = HEOS.viscosity(); //[Pa-s]
OMEGA_tilde = 2.0/pi*((cp-cv)/cp*atan(zeta*qD)+cv/cp*(zeta*qD)); //[-]
OMEGA_tilde0 = 2.0/pi*(1.0-exp(-1.0/(1.0/(qD*zeta)+1.0/3.0*(zeta*qD)*(zeta*qD)/delta/delta))); //[-]
OMEGA_tilde = 2.0/pi*((cp-cv)/cp*atan(zeta*qD)+cv/cp*(zeta*qD)); //[-]
OMEGA_tilde0 = 2.0/pi*(1.0-exp(-1.0/(1.0/(qD*zeta)+1.0/3.0*(zeta*qD)*(zeta*qD)/delta/delta))); //[-]
double lambda = HEOS.rhomolar()*cp*R0*k*HEOS.T()/(6*pi*mu*zeta)*(OMEGA_tilde - OMEGA_tilde0); //[W/m/K]
return lambda; //[W/m/K]
double lambda = HEOS.rhomolar()*cp*R0*k*HEOS.T()/(6*pi*mu*zeta)*(OMEGA_tilde - OMEGA_tilde0); //[W/m/K]
return lambda; //[W/m/K]
}
else{
throw NotImplementedError("TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers is only for pure and pseudo-pure");
@@ -688,31 +688,31 @@ long double TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCR
long double TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMixtureBackend &HEOS){
double e_k = 251.196, Tstar;
double b[] = {0.4226159, 0.6280115, -0.5387661, 0.6735941, 0, 0, -0.4362677, 0.2255388};
double c[] = {0, 2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
double b[] = {0.4226159, 0.6280115, -0.5387661, 0.6735941, 0, 0, -0.4362677, 0.2255388};
double c[] = {0, 2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
//Vesovic Eq. 31 [no units]
double summer = 0;
for (int i=1; i<=5; i++)
summer += c[i]*pow(HEOS.T()/100.0, 2-i);
double cint_k = 1.0 + exp(-183.5/HEOS.T())*summer;
//Vesovic Eq. 31 [no units]
double summer = 0;
for (int i=1; i<=5; i++)
summer += c[i]*pow(HEOS.T()/100.0, 2-i);
double cint_k = 1.0 + exp(-183.5/HEOS.T())*summer;
//Vesovic Eq. 12 [no units]
double r = sqrt(2.0/5.0*cint_k);
//Vesovic Eq. 12 [no units]
double r = sqrt(2.0/5.0*cint_k);
// According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
// between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
// does not yield exactly the correct values
Tstar = HEOS.T()/e_k;
//Vesovic Eq. 30 [no units]
summer = 0;
for (int i=0; i<=7; i++)
summer += b[i]/pow(Tstar, i);
double Gstar_lambda = summer;
Tstar = HEOS.T()/e_k;
//Vesovic Eq. 30 [no units]
summer = 0;
for (int i=0; i<=7; i++)
summer += b[i]/pow(Tstar, i);
double Gstar_lambda = summer;
//Vesovic Eq. 29 [W/m/K]
double lambda_0 = 0.475598e-3*sqrt(HEOS.T())*(1+r*r)/Gstar_lambda;
//Vesovic Eq. 29 [W/m/K]
double lambda_0 = 0.475598e-3*sqrt(HEOS.T())*(1+r*r)/Gstar_lambda;
return lambda_0;
}
@@ -720,8 +720,8 @@ long double TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMix
long double TransportRoutines::conductivity_dilute_hardcoded_ethane(HelmholtzEOSMixtureBackend &HEOS){
double e_k = 245.0;
double tau = 305.33/HEOS.T(), Tstar = HEOS.T()/e_k;
double fint = 1.7104147-0.6936482/Tstar;
double tau = 305.33/HEOS.T(), Tstar = HEOS.T()/e_k;
double fint = 1.7104147-0.6936482/Tstar;
double lambda_0 = 0.276505e-3*(HEOS.calc_viscosity_dilute()*1e6)*(3.75-fint*(tau*tau*HEOS.d2alpha0_dTau2()+1.5)); //[W/m/K]
return lambda_0;
@@ -735,9 +735,9 @@ long double TransportRoutines::conductivity_dilute_eta0_and_poly(HelmholtzEOSMix
double eta0_uPas = HEOS.calc_viscosity_dilute()*1e6; // [uPa-s]
double summer = E.A[0]*eta0_uPas;
for (std::size_t i=1; i < E.A.size(); ++i)
summer += E.A[i]*pow(static_cast<long double>(HEOS.tau()), E.t[i]);
return summer;
for (std::size_t i=1; i < E.A.size(); ++i)
summer += E.A[i]*pow(static_cast<long double>(HEOS.tau()), E.t[i]);
return summer;
}
else{
throw NotImplementedError("TransportRoutines::conductivity_dilute_eta0_and_poly is only for pure and pseudo-pure");
@@ -747,146 +747,146 @@ long double TransportRoutines::conductivity_dilute_eta0_and_poly(HelmholtzEOSMix
long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureBackend &HEOS){
double L[5][6] = {{1.60397357,-0.646013523,0.111443906,0.102997357,-0.0504123634,0.00609859258},
{2.33771842,-2.78843778,1.53616167,-0.463045512,0.0832827019,-0.00719201245},
{2.19650529,-4.54580785,3.55777244,-1.40944978,0.275418278,-0.0205938816},
{-1.21051378,1.60812989,-0.621178141,0.0716373224,0,0},
{-2.7203370,4.57586331,-3.18369245,1.1168348,-0.19268305,0.012913842}};
{2.33771842,-2.78843778,1.53616167,-0.463045512,0.0832827019,-0.00719201245},
{2.19650529,-4.54580785,3.55777244,-1.40944978,0.275418278,-0.0205938816},
{-1.21051378,1.60812989,-0.621178141,0.0716373224,0,0},
{-2.7203370,4.57586331,-3.18369245,1.1168348,-0.19268305,0.012913842}};
double lambdabar_0,lambdabar_1,lambdabar_2,rhobar,Tbar,sum;
double Tstar=647.096,rhostar=322,pstar=22064000,lambdastar=1e-3,mustar=1e-6;
double xi;
int i,j;
double lambdabar_0,lambdabar_1,lambdabar_2,rhobar,Tbar,sum;
double Tstar=647.096,rhostar=322,pstar=22064000,lambdastar=1e-3,mustar=1e-6;
double xi;
int i,j;
double R = 461.51805; //[J/kg/K]
Tbar = HEOS.T()/Tstar;
rhobar = HEOS.keyed_output(CoolProp::iDmass)/rhostar;
Tbar = HEOS.T()/Tstar;
rhobar = HEOS.keyed_output(CoolProp::iDmass)/rhostar;
double rho = HEOS.keyed_output(CoolProp::iDmass);
// Dilute gas contribution
lambdabar_0 = sqrt(Tbar)/(2.443221e-3+1.323095e-2/Tbar+6.770357e-3/pow(Tbar,2)-3.454586e-3/pow(Tbar,3)+4.096266e-4/pow(Tbar,4));
// Dilute gas contribution
lambdabar_0 = sqrt(Tbar)/(2.443221e-3+1.323095e-2/Tbar+6.770357e-3/pow(Tbar,2)-3.454586e-3/pow(Tbar,3)+4.096266e-4/pow(Tbar,4));
sum=0;
for (i=0;i<=4;i++){
for (j=0;j<=5;j++){
sum+=L[i][j]*powInt(1.0/Tbar-1.0,i)*powInt(rhobar-1,j);
}
}
// Finite density contribution
lambdabar_1=exp(rhobar*sum);
sum=0;
for (i=0;i<=4;i++){
for (j=0;j<=5;j++){
sum+=L[i][j]*powInt(1.0/Tbar-1.0,i)*powInt(rhobar-1,j);
}
}
// Finite density contribution
lambdabar_1=exp(rhobar*sum);
double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5,
double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5,
qd_bar=1/0.4,pi=3.141592654, delta = HEOS.delta();
double drhodp = 1/(R*HEOS.T()*(1+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2()));
double drhobar_dpbar = pstar/rhostar*drhodp;
double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta)));
double drhobar_dpbar_Trbar = pstar/rhostar*drhodp_Trbar;
double drhodp = 1/(R*HEOS.T()*(1+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2()));
double drhobar_dpbar = pstar/rhostar*drhodp;
double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta)));
double drhobar_dpbar_Trbar = pstar/rhostar*drhodp_Trbar;
double cp = HEOS.cpmass(); // [J/kg/K]
double cv = HEOS.cvmass(); // [J/kg/K]
double cpbar = cp/R; //[-]
double mubar = HEOS.viscosity()/mustar;
double DELTAchibar_T = rhobar*(drhobar_dpbar-drhobar_dpbar_Trbar*Tr_bar/Tbar);
if (DELTAchibar_T<0)
xi = 0;
else
xi = xi_0*pow(DELTAchibar_T/Lambda_0,nu/gamma);
double y = qd_bar*xi;
double cv = HEOS.cvmass(); // [J/kg/K]
double cpbar = cp/R; //[-]
double mubar = HEOS.viscosity()/mustar;
double DELTAchibar_T = rhobar*(drhobar_dpbar-drhobar_dpbar_Trbar*Tr_bar/Tbar);
if (DELTAchibar_T<0)
xi = 0;
else
xi = xi_0*pow(DELTAchibar_T/Lambda_0,nu/gamma);
double y = qd_bar*xi;
double Z;
double kappa = cp/cv;
if (y < 1.2e-7)
Z = 0;
else
Z = 2/(pi*y)*(((1-1/kappa)*atan(y)+y/kappa)-(1-exp(-1/(1/y+y*y/3/rhobar/rhobar))));
double Z;
double kappa = cp/cv;
if (y < 1.2e-7)
Z = 0;
else
Z = 2/(pi*y)*(((1-1/kappa)*atan(y)+y/kappa)-(1-exp(-1/(1/y+y*y/3/rhobar/rhobar))));
lambdabar_2 = GAMMA*rhobar*cpbar*Tbar/mubar*Z;
lambdabar_2 = GAMMA*rhobar*cpbar*Tbar/mubar*Z;
return (lambdabar_0*lambdabar_1+lambdabar_2)*lambdastar;
return (lambdabar_0*lambdabar_1+lambdabar_2)*lambdastar;
}
long double TransportRoutines::conductivity_hardcoded_R23(HelmholtzEOSMixtureBackend &HEOS){
double B1 = -2.5370, // [mW/m/K]
B2 = 0.05366, // [mW/m/K^2]
C1 = 0.94215, // [-]
C2 = 0.14914, // [mW/m/K^2]
DeltaGstar = 2508.58, //[J/mol]
rhoL = 68.345, // [mol/dm^3] = [mol/L]
rhocbar = 7.5114, // [mol/dm^3]
DELTAlambda_max = 25, //[mW/m/K]
Ru = 8.31451, // [J/mol/K]
B2 = 0.05366, // [mW/m/K^2]
C1 = 0.94215, // [-]
C2 = 0.14914, // [mW/m/K^2]
DeltaGstar = 2508.58, //[J/mol]
rhoL = 68.345, // [mol/dm^3] = [mol/L]
rhocbar = 7.5114, // [mol/dm^3]
DELTAlambda_max = 25, //[mW/m/K]
Ru = 8.31451, // [J/mol/K]
Tc = 299.2793, //[K]
T = HEOS.T(); //[K]
double lambda_DG = B1 + B2*T;
double lambda_DG = B1 + B2*T;
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double lambda_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(T)*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*T));
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double lambda_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(T)*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*T));
double chi = rhobar - rhocbar;
double tau = T - Tc;
double chi = rhobar - rhocbar;
double tau = T - Tc;
double DELTAlambda_c = 4*DELTAlambda_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
double DELTAlambda_c = 4*DELTAlambda_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
return (pow((rhoL-rhobar)/rhoL,C1)*lambda_DG+pow(rhobar/rhoL,C1)*lambda_L+DELTAlambda_c)/1e3;
return (pow((rhoL-rhobar)/rhoL,C1)*lambda_DG+pow(rhobar/rhoL,C1)*lambda_L+DELTAlambda_c)/1e3;
}
long double TransportRoutines::conductivity_critical_hardcoded_ammonia(HelmholtzEOSMixtureBackend &HEOS){
/*
From "Thermal Conductivity of Ammonia in a Large
Temperature and Pressure Range Including the Critical Region"
by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
*/
From "Thermal Conductivity of Ammonia in a Large
Temperature and Pressure Range Including the Critical Region"
by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
*/
double T = HEOS.T(), Tc = 405.4, rhoc = 235, rho;
double LAMBDA=1.2, nu=0.63, gamma =1.24, DELTA=0.50,t,zeta_0_plus=1.34e-10,a_zeta=1,GAMMA_0_plus=0.423e-8;
double pi=3.141592654,a_chi,k_B=1.3806504e-23,X_T,DELTA_lambda,dPdT,eta_B,DELTA_lambda_id,DELTA_lambda_i;
double LAMBDA=1.2, nu=0.63, gamma =1.24, DELTA=0.50,t,zeta_0_plus=1.34e-10,a_zeta=1,GAMMA_0_plus=0.423e-8;
double pi=3.141592654,a_chi,k_B=1.3806504e-23,X_T,DELTA_lambda,dPdT,eta_B,DELTA_lambda_id,DELTA_lambda_i;
rho = HEOS.keyed_output(CoolProp::iDmass);
t = std::abs((T-Tc)/Tc);
a_chi = a_zeta/0.7;
eta_B = (2.60+1.6*t)*1e-5;
dPdT = (2.18-0.12/exp(17.8*t))*1e5; // [Pa-K]
X_T = 0.61*rhoc+16.5*log(t);
// Along the critical isochore (only a function of temperature) (Eq. 9)
DELTA_lambda_i = LAMBDA*(k_B*T*T)/(6*pi*eta_B*(zeta_0_plus*pow(t,-nu)*(1+a_zeta*pow(t,DELTA))))*dPdT*dPdT*GAMMA_0_plus*pow(t,-gamma)*(1+a_chi*pow(t,DELTA));
DELTA_lambda_id = DELTA_lambda_i*exp(-36*t*t);
if (rho < 0.6*rhoc)
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(0.6*rhoc-0.96*rhoc,2))*powInt(rho,2)/powInt(0.6*rhoc,2);
}
else
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(rho-0.96*rhoc,2));
}
t = std::abs((T-Tc)/Tc);
a_chi = a_zeta/0.7;
eta_B = (2.60+1.6*t)*1e-5;
dPdT = (2.18-0.12/exp(17.8*t))*1e5; // [Pa-K]
X_T = 0.61*rhoc+16.5*log(t);
// Along the critical isochore (only a function of temperature) (Eq. 9)
DELTA_lambda_i = LAMBDA*(k_B*T*T)/(6*pi*eta_B*(zeta_0_plus*pow(t,-nu)*(1+a_zeta*pow(t,DELTA))))*dPdT*dPdT*GAMMA_0_plus*pow(t,-gamma)*(1+a_chi*pow(t,DELTA));
DELTA_lambda_id = DELTA_lambda_i*exp(-36*t*t);
if (rho < 0.6*rhoc)
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(0.6*rhoc-0.96*rhoc,2))*powInt(rho,2)/powInt(0.6*rhoc,2);
}
else
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(rho-0.96*rhoc,2));
}
return DELTA_lambda;
return DELTA_lambda;
}
long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixtureBackend &HEOS){
/*
What an incredibly annoying formulation! Implied coefficients?? Not cool.
*/
What an incredibly annoying formulation! Implied coefficients?? Not cool.
*/
double rhoc = 68.0, lambda_e, lambda_c, T = HEOS.T(), rho = HEOS.rhomass();
double summer = 3.739232544/T-2.620316969e1/T/T+5.982252246e1/T/T/T-4.926397634e1/T/T/T/T;
double lambda_0 = 2.7870034e-3*pow(T, 7.034007057e-1)*exp(summer);
double c[]={ 1.862970530e-4,
-7.275964435e-7,
-1.427549651e-4,
3.290833592e-5,
-5.213335363e-8,
4.492659933e-8,
-5.924416513e-9,
7.087321137e-6,
-6.013335678e-6,
8.067145814e-7,
3.995125013e-7};
// Equation 17
lambda_e = (c[0]+c[1]*T+c[2]*pow(T,1/3.0)+c[3]*pow(T,2.0/3.0))*rho
+(c[4]+c[5]*pow(T,1.0/3.0)+c[6]*pow(T,2.0/3.0))*rho*rho*rho
+(c[7]+c[8]*pow(T,1.0/3.0)+c[9]*pow(T,2.0/3.0)+c[10]/T)*rho*rho*log(rho/rhoc);
double summer = 3.739232544/T-2.620316969e1/T/T+5.982252246e1/T/T/T-4.926397634e1/T/T/T/T;
double lambda_0 = 2.7870034e-3*pow(T, 7.034007057e-1)*exp(summer);
double c[]={ 1.862970530e-4,
-7.275964435e-7,
-1.427549651e-4,
3.290833592e-5,
-5.213335363e-8,
4.492659933e-8,
-5.924416513e-9,
7.087321137e-6,
-6.013335678e-6,
8.067145814e-7,
3.995125013e-7};
// Equation 17
lambda_e = (c[0]+c[1]*T+c[2]*pow(T,1/3.0)+c[3]*pow(T,2.0/3.0))*rho
+(c[4]+c[5]*pow(T,1.0/3.0)+c[6]*pow(T,2.0/3.0))*rho*rho*rho
+(c[7]+c[8]*pow(T,1.0/3.0)+c[9]*pow(T,2.0/3.0)+c[10]/T)*rho*rho*log(rho/rhoc);
// Critical component
lambda_c = 0.0;
@@ -928,9 +928,9 @@ long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixture
// 3.4685233d-17 and 3.726229668d0 are "magical" coefficients that are present in the REFPROP source to yield the right values. Not clear why these values are needed.
// Also, the form of the critical term in REFPROP does not agree with Hands paper. EL and MH from NIST are not sure where these coefficients come from.
lambda_c = 3.4685233e-17*3.726229668*sqrt(K_Tbar)*pow(T,2)/rho/eta*pow(dpdT,2)*exp(-18.66*pow(DeltaT,2)-4.25*pow(DeltaRho,4));
lambda_c = 3.4685233e-17*3.726229668*sqrt(K_Tbar)*pow(T,2)/rho/eta*pow(dpdT,2)*exp(-18.66*pow(DeltaT,2)-4.25*pow(DeltaRho,4));
}
return lambda_0+lambda_e+lambda_c;
return lambda_0+lambda_e+lambda_c;
}
long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference)
@@ -962,7 +962,7 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H
// The equivalent substance reducing ratios
long double f = Tc/Tc0*theta;
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
// To be solved for
long double T0 = HEOS.T()/f;
@@ -975,10 +975,10 @@ long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, H
long double eta_resid = HEOS_Reference.calc_viscosity_background();
// The F factor
long double F_eta = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M/M0);
long double F_eta = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M/M0);
// The total viscosity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
long double eta = eta_dilute + eta_resid*F_eta;
long double eta = eta_dilute + eta_resid*F_eta;
return eta;
}
@@ -1031,7 +1031,7 @@ long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS
// The equivalent substance reducing ratios
long double f = Tc/Tc0*theta;
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
// To be solved for
long double T0 = HEOS.T()/f;
@@ -1044,13 +1044,13 @@ long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS
long double lambda_resid = HEOS_Reference.calc_conductivity_background();
// The F factor
long double F_lambda = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M0/M);
long double F_lambda = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M0/M);
// The critical contribution from the fluid of interest [W/m/K]
long double lambda_critical = conductivity_critical_simplified_Olchowy_Sengers(HEOS);
// The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
long double lambda = lambda_int + lambda_dilute + lambda_resid*F_lambda + lambda_critical;
long double lambda = lambda_int + lambda_dilute + lambda_resid*F_lambda + lambda_critical;
return lambda;
}

102
src/Backends/Helmholtz/VLERoutines.cpp
View File

@@ -193,22 +193,22 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
{
// Invert liquid density ancillary to get temperature
// TODO: fit inverse ancillaries too
try{
T = HEOS.get_components()[0]->ancillaries.pL.invert(specified_value);
}
catch(std::exception &e)
{
throw ValueError("Unable to invert ancillary equation");
}
try{
T = HEOS.get_components()[0]->ancillaries.pL.invert(specified_value);
}
catch(std::exception &e)
{
throw ValueError("Unable to invert ancillary equation");
}
}
else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HL)
{
CoolProp::SimpleState hs_anchor = HEOS.get_state("hs_anchor");
// Ancillary is deltah = h - hs_anchor.h
try{ T = HEOS.get_components()[0]->ancillaries.hL.invert(specified_value - hs_anchor.hmolar); }
catch(std::exception &e){
throw ValueError("Unable to invert ancillary equation for hL");
}
catch(std::exception &e){
throw ValueError("Unable to invert ancillary equation for hL");
}
}
else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_HV)
{
@@ -235,12 +235,12 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
double Tmin = Tmin_satL;
double Tmax = HEOS.calc_Tmax_sat();
try{ T = Brent(resid, Tmin-3, Tmax + 1, DBL_EPSILON, 1e-10, 50, errstr); }
catch(std::exception &e){
catch(std::exception &e){
shared_ptr<HelmholtzEOSMixtureBackend> HEOS_copy(new HelmholtzEOSMixtureBackend(HEOS.get_components()));
HEOS_copy->update(QT_INPUTS, 1, Tmin); double hTmin = HEOS_copy->hmolar();
HEOS_copy->update(QT_INPUTS, 1, Tmax); double hTmax = HEOS_copy->hmolar();
T = (Tmax-Tmin)/(hTmax-hTmin)*(HEOS.hmolar()-hTmin) + Tmin;
}
}
}
else if (options.specified_variable == saturation_PHSU_pure_options::IMPOSED_SL)
{
@@ -457,7 +457,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
}
// dr_3/ddelta'' (liquid pressure not a function of vapor density)
J[2][2] = 0;
specified_parameter = CoolProp::iP;
specified_parameter = CoolProp::iP;
break;
case saturation_PHSU_pure_options::IMPOSED_PV:
// dr_3/dtau
@@ -472,7 +472,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
// dr_3/ddelta''
J[2][2] = SatV->first_partial_deriv(iP,iDelta,iTau)/specified_value;
}
specified_parameter = CoolProp::iP;
specified_parameter = CoolProp::iP;
break;
case saturation_PHSU_pure_options::IMPOSED_HL:
// dr_3/dtau
@@ -482,7 +482,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
if (options.use_logdelta){ J[2][1]*=deltaL;}
// dr_3/ddelta''
J[2][2] = 0; //(liquid enthalpy not a function of vapor density)
specified_parameter = CoolProp::iHmolar;
specified_parameter = CoolProp::iHmolar;
break;
case saturation_PHSU_pure_options::IMPOSED_HV:
// dr_3/dtau
@@ -492,7 +492,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
// dr_3/ddelta''
J[2][2] = SatV->first_partial_deriv(iHmolar,iDelta,iTau);
if (options.use_logdelta){ J[2][2]*=deltaV;}
specified_parameter = CoolProp::iHmolar;
specified_parameter = CoolProp::iHmolar;
break;
case saturation_PHSU_pure_options::IMPOSED_SL:
// dr_3/dtau
@@ -502,7 +502,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
if (options.use_logdelta){ J[2][1] *= deltaL; }
// dr_3/ddelta''
J[2][2] = 0; //(liquid entropy not a function of vapor density)
specified_parameter = CoolProp::iSmolar;
specified_parameter = CoolProp::iSmolar;
break;
case saturation_PHSU_pure_options::IMPOSED_SV:
// dr_3/dtau
@@ -512,7 +512,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
// dr_3/ddelta''
J[2][2] = SatV->first_partial_deriv(iSmolar,iDelta,iTau);
if (options.use_logdelta){ J[2][2]*=deltaV;}
specified_parameter = CoolProp::iSmolar;
specified_parameter = CoolProp::iSmolar;
break;
default:
throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable));
@@ -557,7 +557,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, l
// Set values back into the options structure for use in next solver
options.rhoL = rhoL; options.rhoV = rhoV; options.T = T;
// Error out
std::string info = get_parameter_information(specified_parameter, "short");
std::string info = get_parameter_information(specified_parameter, "short");
throw SolutionError(format("saturation_PHSU_pure solver did not converge after 50 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error));
}
}
@@ -704,8 +704,8 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long
rhoL = deltaL*reduce.rhomolar;
rhoV = deltaV*reduce.rhomolar;
T = reduce.T/tau;
p_error = (pL-pV)/pL;
p_error = (pL-pV)/pL;
error = sqrt(pow(r[0], 2)+pow(r[1], 2));
iter++;
@@ -718,10 +718,10 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long
}
}
while (error > 1e-9);
long double p_error_limit = 1e-3;
if (std::abs(p_error) > p_error_limit){
throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
}
long double p_error_limit = 1e-3;
if (std::abs(p_error) > p_error_limit){
throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit));
}
}
void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options)
{
@@ -783,8 +783,8 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE
}
else
{
// Use the density ancillary function as the starting point for the solver
// Use the density ancillary function as the starting point for the solver
// If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
if (T > 0.99*HEOS.get_reducing_state().T){
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T-0.1);
@@ -793,13 +793,13 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE
else{
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
// Update the guess for liquid density using density solver with vapor pressure
// and liquid density guess from ancillaries
@@ -827,8 +827,8 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE
//tau = reduce.T/T;
}
//if (get_debug_level()>5){
// std::cout << format("%s:%d: Akasaka guess values deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
// }
// std::cout << format("%s:%d: Akasaka guess values deltaL = %g deltaV = %g tau = %g\n",__FILE__,__LINE__,deltaL, deltaV, tau).c_str();
// }
do{
/*if (get_debug_level()>8){
@@ -888,16 +888,16 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HE
}
}
while (error > 1e-10 && std::abs(stepL) > 10*DBL_EPSILON*std::abs(stepL) && std::abs(stepV) > 10*DBL_EPSILON*std::abs(stepV));
long double p_error_limit = 1e-3;
long double p_error = (PL - PV)/PL;
if (std::abs(p_error) > p_error_limit){
long double p_error_limit = 1e-3;
long double p_error = (PL - PV)/PL;
if (std::abs(p_error) > p_error_limit){
options.pL = PL;
options.pV = PV;
options.rhoL = rhoL;
options.rhoV = rhoV;
throw SolutionError(format("saturation_T_pure_Akasaka solver abs error on p [%g] > limit [%g]", std::abs(p_error), p_error_limit));
}
throw SolutionError(format("saturation_T_pure_Akasaka solver abs error on p [%g] > limit [%g]", std::abs(p_error), p_error_limit));
}
}
long double sign(long double x)
@@ -936,8 +936,8 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
}
else
{
// Use the density ancillary function as the starting point for the solver
// Use the density ancillary function as the starting point for the solver
// If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature
if (T > 0.9999*HEOS.get_reducing_state().T){
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T-0.1);
@@ -946,13 +946,13 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
else{
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
// Apply a single step of Newton's method to improve guess value for liquid
// based on the error between the gas pressure (which is usually very close already)
// and the liquid pressure, which can sometimes (especially at low pressure),
// be way off, and often times negative
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
// Update the guess for liquid density using density solver with vapor pressure
// and liquid density guess from ancillaries, but only if the pressures are not

182
src/Backends/Helmholtz/VLERoutines.h
View File

@@ -37,15 +37,15 @@ namespace SaturationSolvers
};
/*! Returns the natural logarithm of K for component i using the method from Wilson as in
\f[
\ln K_i = \ln\left(\frac{p_{c,i}}{p}\right)+5.373(1+\omega_i)\left(1-\frac{T_{c,i}}{T}\right)
\f]
\f[
\ln K_i = \ln\left(\frac{p_{c,i}}{p}\right)+5.373(1+\omega_i)\left(1-\frac{T_{c,i}}{T}\right)
\f]
@param HEOS The Helmholtz EOS mixture backend
@param T Temperature [K]
@param p Pressure [Pa]
@param i Index of component [-]
*/
static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend &HEOS, long double T, long double p, int i){
@param T Temperature [K]
@param p Pressure [Pa]
@param i Index of component [-]
*/
static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend &HEOS, long double T, long double p, int i){
EquationOfState *EOS = (HEOS.get_components())[i]->pEOS;
return log(EOS->reduce.p/p)+5.373*(1 + EOS->accentric)*(1-EOS->reduce.T/T);
};
@@ -130,12 +130,12 @@ namespace SaturationSolvers
{
public:
int input_type;
double T, p, beta;
const std::vector<long double> *z;
double T, p, beta;
const std::vector<long double> *z;
std::vector<long double> *K;
HelmholtzEOSMixtureBackend *HEOS;
HelmholtzEOSMixtureBackend *HEOS;
WilsonK_resid(HelmholtzEOSMixtureBackend &HEOS, double beta, double imposed_value, int input_type, const std::vector<long double> &z, std::vector<long double> &K){
WilsonK_resid(HelmholtzEOSMixtureBackend &HEOS, double beta, double imposed_value, int input_type, const std::vector<long double> &z, std::vector<long double> &K){
this->z = &z; this->K = &K; this->HEOS = &HEOS; this->beta = beta; this->input_type = input_type;
if (input_type == imposed_T){
this->T = imposed_value;
@@ -143,35 +143,35 @@ namespace SaturationSolvers
else{
this->p = imposed_value;
}
};
double call(double input_value){
double summer = 0;
};
double call(double input_value){
double summer = 0;
if (input_type == imposed_T){
p = input_value; // Iterate on pressure
}
else{
T = input_value; // Iterate on temperature, pressure imposed
}
for (unsigned int i = 0; i< (*z).size(); i++) {
(*K)[i] = exp(Wilson_lnK_factor(*HEOS,T,p,i));
summer += (*z)[i]*((*K)[i]-1)/(1-beta+beta*(*K)[i]);
}
return summer;
};
for (unsigned int i = 0; i< (*z).size(); i++) {
(*K)[i] = exp(Wilson_lnK_factor(*HEOS,T,p,i));
summer += (*z)[i]*((*K)[i]-1)/(1-beta+beta*(*K)[i]);
}
return summer;
};
};
inline double saturation_preconditioner(HelmholtzEOSMixtureBackend &HEOS, double input_value, int input_type, const std::vector<long double> &z)
{
double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;
for (unsigned int i = 0; i < z.size(); i++)
{
double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0;
for (unsigned int i = 0; i < z.size(); i++)
{
EquationOfState *EOS = (HEOS.get_components())[i]->pEOS;
ptriple += EOS->sat_min_liquid.p*z[i];
ptriple += EOS->sat_min_liquid.p*z[i];
pcrit += EOS->reduce.p*z[i];
Ttriple += EOS->sat_min_liquid.T*z[i];
Tcrit += EOS->reduce.T*z[i];
}
Ttriple += EOS->sat_min_liquid.T*z[i];
Tcrit += EOS->reduce.T*z[i];
}
if (input_type == imposed_T)
{
@@ -185,16 +185,16 @@ namespace SaturationSolvers
}
inline double saturation_Wilson(HelmholtzEOSMixtureBackend &HEOS, double beta, double input_value, int input_type, const std::vector<long double> &z, double guess)
{
double T;
double T;
std::string errstr;
std::string errstr;
// Find first guess for T using Wilson K-factors
// Find first guess for T using Wilson K-factors
WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K());
T = Secant(Resid, guess, 0.001, 1e-10, 100, errstr);
if (!ValidNumber(T)){throw ValueError("saturation_p_Wilson failed to get good T");}
return T;
T = Secant(Resid, guess, 0.001, 1e-10, 100, errstr);
if (!ValidNumber(T)){throw ValueError("saturation_p_Wilson failed to get good T");}
return T;
}
struct SuccessiveSubstitutionStep
{
@@ -245,41 +245,41 @@ namespace SaturationSolvers
{
public:
newton_raphson_twophase_options::imposed_variable_options imposed_variable;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta;
unsigned int N;
bool logging;
int Nsteps;
STLMatrix J;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change, beta;
unsigned int N;
bool logging;
int Nsteps;
STLMatrix J;
HelmholtzEOSMixtureBackend *HEOS;
std::vector<long double> K, x, y, z, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
std::vector<long double> K, x, y, z, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
newton_raphson_twophase(){};
newton_raphson_twophase(){};
void resize(unsigned int N);
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear(); step_logger.clear(); error_rms = 1e99; Nsteps = 0;
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
void resize(unsigned int N);
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear(); step_logger.clear(); error_rms = 1e99; Nsteps = 0;
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
/** \brief Call the Newton-Raphson VLE Solver
/** \brief Call the Newton-Raphson VLE Solver
*
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
*
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param IO The input/output data structure
*/
void call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO);
*/
void call(HelmholtzEOSMixtureBackend &HEOS, newton_raphson_twophase_options &IO);
/* \brief Build the arrays for the Newton-Raphson solve
/* \brief Build the arrays for the Newton-Raphson solve
*
*/
void build_arrays();
*/
void build_arrays();
};
@@ -325,48 +325,48 @@ namespace SaturationSolvers
{
public:
newton_raphson_saturation_options::imposed_variable_options imposed_variable;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change;
unsigned int N;
bool logging;
long double error_rms, rhomolar_liq, rhomolar_vap, T, p, min_rel_change;
unsigned int N;
bool logging;
bool bubble_point;
int Nsteps;
STLMatrix J;
int Nsteps;
STLMatrix J;
HelmholtzEOSMixtureBackend *HEOS;
long double dTsat_dPsat, dPsat_dTsat;
std::vector<long double> K, x, y, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
std::vector<long double> K, x, y, r, negative_r, err_rel;
std::vector<SuccessiveSubstitutionStep> step_logger;
newton_raphson_saturation(){};
newton_raphson_saturation(){};
void resize(unsigned int N);
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear();
void resize(unsigned int N);
// Reset the state of all the internal variables
void pre_call()
{
K.clear(); x.clear(); y.clear();
step_logger.clear(); error_rms = 1e99; Nsteps = 0;
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
/** \brief Call the Newton-Raphson VLE Solver
/** \brief Call the Newton-Raphson VLE Solver
*
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
*
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param z Bulk mole fractions [-]
* @param z_incipient Initial guesses for the mole fractions of the incipient phase [-]
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param z Bulk mole fractions [-]
* @param z_incipient Initial guesses for the mole fractions of the incipient phase [-]
* @param IO The input/output data structure
*/
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO);
*/
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO);
/** \brief Build the arrays for the Newton-Raphson solve
/** \brief Build the arrays for the Newton-Raphson solve
*
* This method builds the Jacobian matrix, the sensitivity matrix, etc.
* This method builds the Jacobian matrix, the sensitivity matrix, etc.
*
*/
void build_arrays();
*/
void build_arrays();
/** \brief Check the derivatives in the Jacobian using numerical derivatives.
*/

988
src/Backends/Incompressible/IncompressibleBackend.cpp
View File

File diff suppressed because it is too large Load Diff

8
src/Backends/Incompressible/IncompressibleBackend.h
View File

@@ -36,8 +36,8 @@ public:
/// @param fluid_name the string with the fluid name
IncompressibleBackend(const std::string &fluid_name);
/// The instantiator
/// @param component_names The vector of strings of the fluid components, without file ending
IncompressibleBackend(const std::vector<std::string> &component_names);
/// @param component_names The vector of strings of the fluid components, without file ending
IncompressibleBackend(const std::vector<std::string> &component_names);
// Incompressible backend uses different compositions
bool using_mole_fractions(void){return this->fluid->getxid()==IFRAC_MOLE;};
@@ -117,8 +117,8 @@ public:
long double calc_umass(void){return fluid->u(_T, _p, _fractions[0]);};
long double calc_cpmass(void){return fluid->cp(_T, _p, _fractions[0]);};
long double calc_cvmass(void){return fluid->cv(_T, _p, _fractions[0]);};
long double calc_Tmax(void){return fluid->getTmax();};
long double calc_Tmin(void){return fluid->getTmin();};
long double calc_Tmax(void){return fluid->getTmax();};
long double calc_Tmin(void){return fluid->getTmin();};
};
} /* namespace CoolProp */

1416
src/Backends/Incompressible/IncompressibleFluid.cpp
View File

File diff suppressed because it is too large Load Diff

648
src/Backends/Incompressible/IncompressibleLibrary.cpp
View File

@@ -35,217 +35,217 @@ double const LiBrSolution::cpt_H2O = 76.0226; /* J/(mol.K), molar isobaric heat
double LiBrSolution::ps_mix(double T, double x)
/* Equation (1) */
{
static double m[8] = { 3.0, 4.0, 4.0, 8.0, 1.0, 1.0, 4.0, 6.0 };
static double n[8] = { 0.0, 5.0, 6.0, 3.0, 0.0, 2.0, 6.0, 0.0 };
static double t[8] = { 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 };
static double a[8] = { -2.41303e2, 1.91750e7, -1.75521e8, 3.25430e7,
3.92571e2, -2.12626e3, 1.85127e8, 1.91216e3 };
double tau, suma = 0.0;
int i;
static double m[8] = { 3.0, 4.0, 4.0, 8.0, 1.0, 1.0, 4.0, 6.0 };
static double n[8] = { 0.0, 5.0, 6.0, 3.0, 0.0, 2.0, 6.0, 0.0 };
static double t[8] = { 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 };
static double a[8] = { -2.41303e2, 1.91750e7, -1.75521e8, 3.25430e7,
3.92571e2, -2.12626e3, 1.85127e8, 1.91216e3 };
double tau, suma = 0.0;
int i;
tau = T / Tc_H2O;
for (i = 0; i <= 7; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
return (ps_H2O(T - suma));
tau = T / Tc_H2O;
for (i = 0; i <= 7; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
return (ps_H2O(T - suma));
} /* end function ps_mix */
double LiBrSolution::rho_mix(double T, double x)
/* Equation (2) */
{
static double m[2] = { 1.0, 1.0 };
static double n[2] = { 0.0, 6.0 };
static double a[2] = { 1.746, 4.709 };
static double m[2] = { 1.0, 1.0 };
static double n[2] = { 0.0, 6.0 };
static double a[2] = { 1.746, 4.709 };
double tau, suma = 0.0;
int i;
double tau, suma = 0.0;
int i;
tau = T / Tc_H2O;
for (i = 0; i <= 1; i++)
suma += a[i] * pow(x, m[i]) * pow(tau, n[i]);
tau = T / Tc_H2O;
for (i = 0; i <= 1; i++)
suma += a[i] * pow(x, m[i]) * pow(tau, n[i]);
return ((1.0 - x) * rho_H2O(T) + rhoc_H2O * suma);
return ((1.0 - x) * rho_H2O(T) + rhoc_H2O * suma);
} /* end function rho_mix */
double LiBrSolution::cp_mix(double T, double x)
/* Equation (3) */
{
static double m[8] = { 2.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0 };
static double n[8] = { 0.0, 0.0, 1.0, 2.0, 3.0, 0.0, 3.0, 2.0 };
static double t[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0 };
static double a[8] = { -1.42094e1, 4.04943e1, 1.11135e2, 2.29980e2,
1.34526e3, -1.41010e-2, 1.24977e-2, -6.83209e-4 };
static double m[8] = { 2.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0 };
static double n[8] = { 0.0, 0.0, 1.0, 2.0, 3.0, 0.0, 3.0, 2.0 };
static double t[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 3.0, 4.0 };
static double a[8] = { -1.42094e1, 4.04943e1, 1.11135e2, 2.29980e2,
1.34526e3, -1.41010e-2, 1.24977e-2, -6.83209e-4 };
double tau, suma = 0.0;
int i;
double tau, suma = 0.0;
int i;
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 7; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 7; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
return ((1.0 - x) * cp_H2O(T) + cpt_H2O * suma);
return ((1.0 - x) * cp_H2O(T) + cpt_H2O * suma);
} /* end function cp_mix */
double LiBrSolution::h_mix(double T, double x)
/* Equation (4) */
{
static double m[30] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 4.0,
5.0, 5.0, 6.0, 6.0, 1.0, 2.0, 2.0, 2.0, 5.0, 6.0, 7.0, 1.0, 1.0,
2.0, 2.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 };
static double m[30] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 4.0,
5.0, 5.0, 6.0, 6.0, 1.0, 2.0, 2.0, 2.0, 5.0, 6.0, 7.0, 1.0, 1.0,
2.0, 2.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 };
static double n[30] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0,
4.0, 5.0, 5.0, 6.0, 0.0, 3.0, 5.0, 7.0, 0.0, 3.0, 1.0, 0.0, 4.0,
2.0, 6.0, 7.0, 0.0, 0.0, 1.0, 2.0, 3.0 };
static double n[30] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0,
4.0, 5.0, 5.0, 6.0, 0.0, 3.0, 5.0, 7.0, 0.0, 3.0, 1.0, 0.0, 4.0,
2.0, 6.0, 7.0, 0.0, 0.0, 1.0, 2.0, 3.0 };
static double t[30] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0,
2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0,
4.0, 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 };
static double t[30] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0,
2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0,
4.0, 4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 };
static double a[30] = { 2.27431, -7.99511, 3.85239e2, -1.63940e4,
-4.22562e2, 1.13314e-1, -8.33474, -1.73833e4, 6.49763,
3.24552e3, -1.34643e4, 3.99322e4, -2.58877e5, -1.93046e-3,
2.80616, -4.04479e1, 1.45342e2, -2.74873, -4.49743e2,
-1.21794e1, -5.83739e-3, 2.33910e-1, 3.41888e-1, 8.85259,
-1.78731e1, 7.35179e-2, -1.79430e-4, 1.84261e-3, -6.24282e-3,
6.84765e-3 };
static double a[30] = { 2.27431, -7.99511, 3.85239e2, -1.63940e4,
-4.22562e2, 1.13314e-1, -8.33474, -1.73833e4, 6.49763,
3.24552e3, -1.34643e4, 3.99322e4, -2.58877e5, -1.93046e-3,
2.80616, -4.04479e1, 1.45342e2, -2.74873, -4.49743e2,
-1.21794e1, -5.83739e-3, 2.33910e-1, 3.41888e-1, 8.85259,
-1.78731e1, 7.35179e-2, -1.79430e-4, 1.84261e-3, -6.24282e-3,
6.84765e-3 };
double tau, suma = 0.0;
int i;
double tau, suma = 0.0;
int i;
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 29; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 29; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
return ((1.0 - x) * h_H2O(T) + hc_H2O * suma);
return ((1.0 - x) * h_H2O(T) + hc_H2O * suma);
} /* end function h_mix */
double LiBrSolution::s_mix(double T, double x)
/* Equation (5) */
{
static double m[29] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 1.0,
2.0, 2.0, 4.0, 5.0, 5.0, 6.0, 6.0, 1.0, 3.0, 5.0, 7.0, 1.0, 1.0,
1.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 };
static double m[29] = { 1.0, 1.0, 2.0, 3.0, 6.0, 1.0, 3.0, 5.0, 1.0,
2.0, 2.0, 4.0, 5.0, 5.0, 6.0, 6.0, 1.0, 3.0, 5.0, 7.0, 1.0, 1.0,
1.0, 2.0, 3.0, 1.0, 1.0, 1.0, 1.0 };
static double n[29] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0,
0.0, 4.0, 0.0, 4.0, 5.0, 2.0, 5.0, 0.0, 4.0, 0.0, 1.0, 0.0, 2.0,
4.0, 7.0, 1.0, 0.0, 1.0, 2.0, 3.0 };
static double n[29] = { 0.0, 1.0, 6.0, 6.0, 2.0, 0.0, 0.0, 4.0, 0.0,
0.0, 4.0, 0.0, 4.0, 5.0, 2.0, 5.0, 0.0, 4.0, 0.0, 1.0, 0.0, 2.0,
4.0, 7.0, 1.0, 0.0, 1.0, 2.0, 3.0 };
static double t[29] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0,
2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0,
4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 };
static double t[29] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0,
2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0,
4.0, 4.0, 4.0, 5.0, 5.0, 5.0, 5.0 };
static double a[29] = { 1.53091, -4.52564, 6.98302e+2, -2.16664e+4,
-1.47533e+3, 8.47012e-2, -6.59523, -2.95331e+4, 9.56314e-3,
-1.88679e-1, 9.31752, 5.78104, 1.38931e+4, -1.71762e+4,
4.15108e+2, -5.55647e+4, -4.23409e-3, 3.05242e+1, -1.67620,
1.48283e+1, 3.03055e-3, -4.01810e-2, 1.49252e-1, 2.59240,
-1.77421e-1, -6.99650e-5, 6.05007e-4, -1.65228e-3, 1.22966e-3 };
static double a[29] = { 1.53091, -4.52564, 6.98302e+2, -2.16664e+4,
-1.47533e+3, 8.47012e-2, -6.59523, -2.95331e+4, 9.56314e-3,
-1.88679e-1, 9.31752, 5.78104, 1.38931e+4, -1.71762e+4,
4.15108e+2, -5.55647e+4, -4.23409e-3, 3.05242e+1, -1.67620,
1.48283e+1, 3.03055e-3, -4.01810e-2, 1.49252e-1, 2.59240,
-1.77421e-1, -6.99650e-5, 6.05007e-4, -1.65228e-3, 1.22966e-3 };
double tau, suma = 0.0;
int i;
double tau, suma = 0.0;
int i;
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 28; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
tau = Tc_H2O / (T - T0);
for (i = 0; i <= 28; i++)
suma += a[i] * pow(x, m[i]) * pow(0.4 - x, n[i]) * pow(tau, t[i]);
return ((1.0 - x) * s_H2O(T) + sc_H2O * suma);
return ((1.0 - x) * s_H2O(T) + sc_H2O * suma);
} /* end function s_mix */
double LiBrSolution::ps_H2O(double T)
/* Equation (28) */
{
static double a[7] = { 0.0, -7.85951783, 1.84408259, -11.7866497,
22.6807411, -15.9618719, 1.80122502 };
static double a[7] = { 0.0, -7.85951783, 1.84408259, -11.7866497,
22.6807411, -15.9618719, 1.80122502 };
double tau, ps;
double tau, ps;
tau = 1 - T / Tc_H2O;
tau = 1 - T / Tc_H2O;
ps = pc_H2O
* exp(
Tc_H2O / T
* (a[1] * tau + a[2] * pow(tau, 1.5)
+ a[3] * pow(tau, 3.0)
+ a[4] * pow(tau, 3.5)
+ a[5] * pow(tau, 4.0)
+ a[6] * pow(tau, 7.5)));
ps = pc_H2O
* exp(
Tc_H2O / T
* (a[1] * tau + a[2] * pow(tau, 1.5)
+ a[3] * pow(tau, 3.0)
+ a[4] * pow(tau, 3.5)
+ a[5] * pow(tau, 4.0)
+ a[6] * pow(tau, 7.5)));
return (ps * 1.0e6);
return (ps * 1.0e6);
} /* end function ps_H2O */
double LiBrSolution::rho_H2O(double T)
/* Equation (29) */
{
static double b[7] = { 0.0, 1.99274064, 1.09965342, -0.510839303,
-1.75493479, -45.5170352, -6.7469445e5 };
double theta, rho;
static double b[7] = { 0.0, 1.99274064, 1.09965342, -0.510839303,
-1.75493479, -45.5170352, -6.7469445e5 };
double theta, rho;
theta = 1.0 - T / Tc_H2O;
theta = 1.0 - T / Tc_H2O;
rho = rhoc_H2O
* (1.0 + b[1] * pow(theta, 1.0 / 3.0)
+ b[2] * pow(theta, 2.0 / 3.0)
+ b[3] * pow(theta, 5.0 / 3.0)
+ b[4] * pow(theta, 16.0 / 3.0)
+ b[5] * pow(theta, 43.0 / 3.0)
+ b[6] * pow(theta, 110.0 / 3.0));
rho = rhoc_H2O
* (1.0 + b[1] * pow(theta, 1.0 / 3.0)
+ b[2] * pow(theta, 2.0 / 3.0)
+ b[3] * pow(theta, 5.0 / 3.0)
+ b[4] * pow(theta, 16.0 / 3.0)
+ b[5] * pow(theta, 43.0 / 3.0)
+ b[6] * pow(theta, 110.0 / 3.0));
return (rho);
return (rho);
} /* end function rho_H2O */
double LiBrSolution::cp_H2O(double T)
/* Equation (30) */
{
static double a[5] =
{ 1.38801, -2.95318, 3.18721, -0.645473, 9.18946e5 };
static double b[5] = { 0.0, 2.0, 3.0, 6.0, 34.0 };
static double c[5] = { 0.0, 2.0, 3.0, 5.0, 0.0 };
static double a[5] =
{ 1.38801, -2.95318, 3.18721, -0.645473, 9.18946e5 };
static double b[5] = { 0.0, 2.0, 3.0, 6.0, 34.0 };
static double c[5] = { 0.0, 2.0, 3.0, 5.0, 0.0 };
double suma = 0.0;
int i;
double suma = 0.0;
int i;
for (i = 0; i <= 4; i++)
suma += a[i] * exp(b[i] * log(1.0 - T / Tc_H2O))
* exp(c[i] * log(T / Tt_H2O));
for (i = 0; i <= 4; i++)
suma += a[i] * exp(b[i] * log(1.0 - T / Tc_H2O))
* exp(c[i] * log(T / Tt_H2O));
return (cpt_H2O * suma);
return (cpt_H2O * suma);
} /* end function cp_H2O */
double LiBrSolution::h_H2O(double T)
/* Equation (31) */
{
static double a[4] = { -4.37196e-1, 3.03440e-1, -1.29582, -1.76410e-1 };
static double alpha[4] = { 1.0 / 3.0, 2.0 / 3.0, 5.0 / 6.0, 21.0 / 6.0 };
static double a[4] = { -4.37196e-1, 3.03440e-1, -1.29582, -1.76410e-1 };
static double alpha[4] = { 1.0 / 3.0, 2.0 / 3.0, 5.0 / 6.0, 21.0 / 6.0 };
double suma = 0.0;
int i;
double suma = 0.0;
int i;
for (i = 0; i <= 3; i++)
suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O));
for (i = 0; i <= 3; i++)
suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O));
return (hc_H2O * (1.0 + suma));
return (hc_H2O * (1.0 + suma));
} /* end function h_H2O */
double LiBrSolution::s_H2O(double T)
/* Equation (32) */
{
static double a[4] = { -3.34112e-1, -8.47987e-1, -9.11980e-1, -1.64046 };
static double alpha[4] = { 1.0 / 3.0, 3.0 / 3.0, 8.0 / 3.0, 24.0 / 3.0 };
static double a[4] = { -3.34112e-1, -8.47987e-1, -9.11980e-1, -1.64046 };
static double alpha[4] = { 1.0 / 3.0, 3.0 / 3.0, 8.0 / 3.0, 24.0 / 3.0 };
double suma = 0.0;
int i;
double suma = 0.0;
int i;
for (i = 0; i <= 3; i++)
suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O));
for (i = 0; i <= 3; i++)
suma += a[i] * exp(alpha[i] * log(1.0 - T / Tc_H2O));
return (sc_H2O * (1.0 + suma));
return (sc_H2O * (1.0 + suma));
} /* end function s_H2O */
@@ -255,16 +255,16 @@ double LiBrSolution::s_H2O(double T)
double LiBrSolution::massToMole(double w)
/* Equation (7) */
{
return (w/M_LiBr)/(w/M_LiBr+(1.-w)/M_H2O);
//return (w*M_LiBr)/(w*M_LiBr+(1.-w)*M_H2O);
return (w/M_LiBr)/(w/M_LiBr+(1.-w)/M_H2O);
//return (w*M_LiBr)/(w*M_LiBr+(1.-w)*M_H2O);
}
double LiBrSolution::molarToSpecific(double w, double value)
/* Equation (7,8) */
{
double x = massToMole(w);
//return w/(x*M_LiBr) * value;
return 1. / ( x*M_LiBr + (1.-x)*M_H2O ) * value;
double x = massToMole(w);
//return w/(x*M_LiBr) * value;
return 1. / ( x*M_LiBr + (1.-x)*M_H2O ) * value;
}
bool const LiBrSolution::debug = false;
@@ -272,68 +272,68 @@ bool const LiBrSolution::debug = false;
LiBrSolution::LiBrSolution():IncompressibleFluid(){
name = std::string("LiBr");
description = std::string("Lithium-Bromide solution from Patek2006");
reference = std::string("Patek2006");
name = std::string("LiBr");
description = std::string("Lithium-Bromide solution from Patek2006");
reference = std::string("Patek2006");
Tmin = 273.00;
Tmax = 500.00;
TminPsat = Tmin;
Tmin = 273.00;
Tmax = 500.00;
TminPsat = Tmin;
xmin = 0.0;
xmax = 1.0;
xmin = 0.0;
xmax = 1.0;
xref = 0.0;
Tref = 0.0;
pref = 0.0;
href = 0.0;
sref = 0.0;
uref = 0.0;
rhoref = 0.0;
xbase = 0.0;
Tbase = 0.0;
xref = 0.0;
Tref = 0.0;
pref = 0.0;
href = 0.0;
sref = 0.0;
uref = 0.0;
rhoref = 0.0;
xbase = 0.0;
Tbase = 0.0;
};
double LiBrSolution::rho(double T, double p, double x){
checkTPX(T, p, x);
return 1./molarToSpecific(x, 1./rho_mix(T,massToMole(x)));
checkTPX(T, p, x);
return 1./molarToSpecific(x, 1./rho_mix(T,massToMole(x)));
}
double LiBrSolution::c(double T, double p, double x){
checkTPX(T, p, x);
return molarToSpecific(x, cp_mix(T,massToMole(x)));
checkTPX(T, p, x);
return molarToSpecific(x, cp_mix(T,massToMole(x)));
}
//double h(double T, double p, double x){
// return h_u(T,p,x);
// return h_u(T,p,x);
//}
double LiBrSolution::s(double T, double p, double x){
checkTPX(T, p, x);
return molarToSpecific(x, s_mix(T,massToMole(x)));
checkTPX(T, p, x);
return molarToSpecific(x, s_mix(T,massToMole(x)));
}
double LiBrSolution::visc(double T, double p, double x){
throw ValueError("Viscosity is not defined for LiBr-solutions.");
throw ValueError("Viscosity is not defined for LiBr-solutions.");
}
double LiBrSolution::cond(double T, double p, double x){
throw ValueError("Thermal conductivity is not defined for LiBr-solutions.");
throw ValueError("Thermal conductivity is not defined for LiBr-solutions.");
}
double LiBrSolution::u(double T, double p, double x){
checkTPX(T, p, x);
return molarToSpecific(x, h_mix(T,massToMole(x)));
checkTPX(T, p, x);
return molarToSpecific(x, h_mix(T,massToMole(x)));
}
double LiBrSolution::psat(double T, double x){
//checkT(T,p,x);
if (debug) throw ValueError(format("Your concentration is %f in kg/kg and %f in mol/mol.",x,massToMole(x)));
return ps_mix(T,massToMole(x));
//checkT(T,p,x);
if (debug) throw ValueError(format("Your concentration is %f in kg/kg and %f in mol/mol.",x,massToMole(x)));
return ps_mix(T,massToMole(x));
};
double LiBrSolution::Tfreeze(double p, double x){
if (debug) throw ValueError(format("No freezing point data available for Lithium-Bromide: p=%f, x=%f",p,x));
return Tmin;
if (debug) throw ValueError(format("No freezing point data available for Lithium-Bromide: p=%f, x=%f",p,x));
return Tmin;
}
/// Default constructor
JSONIncompressibleLibrary::JSONIncompressibleLibrary(){
_is_empty = true;
_is_empty = true;
// fluid_map.clear();
// name_vector.clear();
// string_to_index_map.clear();
@@ -344,168 +344,168 @@ JSONIncompressibleLibrary::JSONIncompressibleLibrary(){
/// Default destructor
JSONIncompressibleLibrary::~JSONIncompressibleLibrary(){
// freeClear(fluid_map);
// fluid_map.clear();
// freeClear(fluid_map);
// fluid_map.clear();
// name_vector.clear();
// string_to_index_map.clear();
};
/// A general function to parse the json files that hold the coefficient matrices
IncompressibleData JSONIncompressibleLibrary::parse_coefficients(rapidjson::Value &obj, std::string id, bool vital){
IncompressibleData fluidData;
if (obj.HasMember(id.c_str())) {
//rapidjson::Value value = obj[id.c_str()];
if (obj[id.c_str()].HasMember("type")){
if (obj[id.c_str()].HasMember("coeffs")){
std::string type = cpjson::get_string(obj[id.c_str()], "type");
if (!type.compare("polynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("logexponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exppolynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("polyoffset")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYOFFSET;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (vital){
throw ValueError(format("The type [%s] is not understood for [%s] of incompressible fluids. Please check your JSON file.", type.c_str(), id.c_str()));
}
else{
IncompressibleData fluidData;
if (obj.HasMember(id.c_str())) {
//rapidjson::Value value = obj[id.c_str()];
if (obj[id.c_str()].HasMember("type")){
if (obj[id.c_str()].HasMember("coeffs")){
std::string type = cpjson::get_string(obj[id.c_str()], "type");
if (!type.compare("polynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("logexponential")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_LOGEXPONENTIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("exppolynomial")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (!type.compare("polyoffset")){
fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYOFFSET;
fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
return fluidData;
}
else if (vital){
throw ValueError(format("The type [%s] is not understood for [%s] of incompressible fluids. Please check your JSON file.", type.c_str(), id.c_str()));
}
else{
//std::cout << format("The type [%s] is not understood for [%s] of incompressible fluids. Please check your JSON file.\n", type.c_str(), id.c_str());
}
}
else{
throw ValueError(format("Your file does not have an entry for \"coeffs\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
throw ValueError(format("Your file does not have an entry for \"type\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
}
return fluidData;
}
}
else{
throw ValueError(format("Your file does not have an entry for \"coeffs\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
throw ValueError(format("Your file does not have an entry for \"type\" in [%s], which is vital for this function.", id.c_str()));
}
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
}
return fluidData;
}
/// Get a double from the JSON storage if it is defined, otherwise return def
double JSONIncompressibleLibrary::parse_value(rapidjson::Value &obj, std::string id, bool vital, double def=0.0){
if (obj.HasMember(id.c_str())) {
if (obj.HasMember(id.c_str())) {
return cpjson::get_double(obj, id);
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
else{
return def;
}
}
else{
if (vital) {
throw ValueError(format("Your file does not have information for [%s], which is vital for an incompressible fluid.", id.c_str()));
}
else{
return def;
}
}
}
/// Get an integer from the JSON storage to identify the composition
composition_types JSONIncompressibleLibrary::parse_ifrac(rapidjson::Value &obj, std::string id){
std::string res = cpjson::get_string(obj, id);
if (!res.compare("mass")) return IFRAC_MASS;
if (!res.compare("mole")) return IFRAC_MOLE;
if (!res.compare("volume")) return IFRAC_VOLUME;
if (!res.compare("not defined")) return IFRAC_UNDEFINED;
if (!res.compare("pure")) return IFRAC_PURE;
std::string res = cpjson::get_string(obj, id);
if (!res.compare("mass")) return IFRAC_MASS;
if (!res.compare("mole")) return IFRAC_MOLE;
if (!res.compare("volume")) return IFRAC_VOLUME;
if (!res.compare("not defined")) return IFRAC_UNDEFINED;
if (!res.compare("pure")) return IFRAC_PURE;
throw ValueError(format("Cannot recognise the entry for [%s], [%s] is unknown for incompressible fluids.", id.c_str(), res.c_str()));
return IFRAC_UNDEFINED;
throw ValueError(format("Cannot recognise the entry for [%s], [%s] is unknown for incompressible fluids.", id.c_str(), res.c_str()));
return IFRAC_UNDEFINED;
}
/// Add all the fluid entries in the rapidjson::Value instance passed in
void JSONIncompressibleLibrary::add_many(rapidjson::Value &listing) {
for (rapidjson::Value::ValueIterator itr = listing.Begin();
itr != listing.End(); ++itr) {
add_one(*itr);
}
for (rapidjson::Value::ValueIterator itr = listing.Begin();
itr != listing.End(); ++itr) {
add_one(*itr);
}
};
void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
_is_empty = false;
_is_empty = false;
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
//fluid_map[index].reset(new IncompressibleFluid());
//fluid_map[index].reset(new IncompressibleFluid());
// Add index->fluid mapping
fluid_map[index] = IncompressibleFluid();
//fluid_map[index].reset(new IncompressibleFluid());
//fluid_map[index].reset(new IncompressibleFluid());
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.setName("unloaded");
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
fluid.setName("unloaded");
try
{
fluid.setName(cpjson::get_string(fluid_json, "name"));
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading base values for %s ",fluid.getName().c_str()) << std::endl;
fluid.setDescription(cpjson::get_string(fluid_json, "description"));
fluid.setReference(cpjson::get_string(fluid_json, "reference"));
fluid.setTmax( parse_value(fluid_json, "Tmax", true, 0.0));
fluid.setTmin( parse_value(fluid_json, "Tmin", true, 0.0));
fluid.setxmax( parse_value(fluid_json, "xmax", false, 1.0));
fluid.setxmin( parse_value(fluid_json, "xmin", false, 0.0));
fluid.setxid( parse_ifrac(fluid_json, "xid") );
fluid.setTminPsat(parse_value(fluid_json, "TminPsat", false, 0.0));
fluid.setName(cpjson::get_string(fluid_json, "name"));
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading base values for %s ",fluid.getName().c_str()) << std::endl;
fluid.setDescription(cpjson::get_string(fluid_json, "description"));
fluid.setReference(cpjson::get_string(fluid_json, "reference"));
fluid.setTmax( parse_value(fluid_json, "Tmax", true, 0.0));
fluid.setTmin( parse_value(fluid_json, "Tmin", true, 0.0));
fluid.setxmax( parse_value(fluid_json, "xmax", false, 1.0));
fluid.setxmin( parse_value(fluid_json, "xmin", false, 0.0));
fluid.setxid( parse_ifrac(fluid_json, "xid") );
fluid.setTminPsat(parse_value(fluid_json, "TminPsat", false, 0.0));
fluid.setTbase(parse_value(fluid_json, "Tbase", false, 0.0));
fluid.setxbase(parse_value(fluid_json, "xbase", false, 0.0));
fluid.setTbase(parse_value(fluid_json, "Tbase", false, 0.0));
fluid.setxbase(parse_value(fluid_json, "xbase", false, 0.0));
/// Setters for the coefficients
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading coefficients for %s ",fluid.getName().c_str()) << std::endl;
fluid.setDensity(parse_coefficients(fluid_json, "density", true));
fluid.setSpecificHeat(parse_coefficients(fluid_json, "specific_heat", true));
fluid.setViscosity(parse_coefficients(fluid_json, "viscosity", false));
fluid.setConductivity(parse_coefficients(fluid_json, "conductivity", false));
fluid.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
fluid.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
fluid.setMass2input(parse_coefficients(fluid_json, "mass2input", false));
fluid.setVolume2input(parse_coefficients(fluid_json, "volume2input", false));
fluid.setMole2input(parse_coefficients(fluid_json, "mole2input", false));
/// Setters for the coefficients
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading coefficients for %s ",fluid.getName().c_str()) << std::endl;
fluid.setDensity(parse_coefficients(fluid_json, "density", true));
fluid.setSpecificHeat(parse_coefficients(fluid_json, "specific_heat", true));
fluid.setViscosity(parse_coefficients(fluid_json, "viscosity", false));
fluid.setConductivity(parse_coefficients(fluid_json, "conductivity", false));
fluid.setPsat(parse_coefficients(fluid_json, "saturation_pressure", false));
fluid.setTfreeze(parse_coefficients(fluid_json, "T_freeze", false));
fluid.setMass2input(parse_coefficients(fluid_json, "mass2input", false));
fluid.setVolume2input(parse_coefficients(fluid_json, "volume2input", false));
fluid.setMole2input(parse_coefficients(fluid_json, "mole2input", false));
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading reference state for %s ",fluid.getName().c_str()) << std::endl;
fluid.set_reference_state(
parse_value(fluid_json, "Tref", false, 20+273.15) ,
parse_value(fluid_json, "pref", false, 1.01325e5) ,
parse_value(fluid_json, "xref", false, 0.0) ,
parse_value(fluid_json, "href", false, 0.0) ,
parse_value(fluid_json, "sref", false, 0.0)
);
if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading reference state for %s ",fluid.getName().c_str()) << std::endl;
fluid.set_reference_state(
parse_value(fluid_json, "Tref", false, 20+273.15) ,
parse_value(fluid_json, "pref", false, 1.01325e5) ,
parse_value(fluid_json, "xref", false, 0.0) ,
parse_value(fluid_json, "href", false, 0.0) ,
parse_value(fluid_json, "sref", false, 0.0)
);
/// A function to check coefficients and equation types.
fluid.validate();
/// A function to check coefficients and equation types.
fluid.validate();
// Add name->index mapping
string_to_index_map[fluid.getName()] = index;
// Add name->index mapping
string_to_index_map[fluid.getName()] = index;
// Add name to vector of names
if (fluid.is_pure()){
this->name_vector_pure.push_back(fluid.getName());
}
else{
this->name_vector_solution.push_back(fluid.getName());
}
if (fluid.is_pure()){
this->name_vector_pure.push_back(fluid.getName());
}
else{
this->name_vector_solution.push_back(fluid.getName());
}
}
catch(std::exception &e)
{
@@ -516,16 +516,16 @@ void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
};
void JSONIncompressibleLibrary::add_obj(IncompressibleFluid fluid_obj) {
_is_empty = false;
_is_empty = false;
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Get the next index for this fluid
std::size_t index = fluid_map.size();
// Add index->fluid mapping
fluid_map[index] = fluid_obj;
// Add index->fluid mapping
fluid_map[index] = fluid_obj;
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
// Create an instance of the fluid
IncompressibleFluid &fluid = fluid_map[index];
/// A function to check coefficients and equation types.
fluid.validate();
@@ -536,20 +536,20 @@ void JSONIncompressibleLibrary::add_obj(IncompressibleFluid fluid_obj) {
// Get an IncompressibleFluid instance stored in this library
IncompressibleFluid& JSONIncompressibleLibrary::get(std::string key) {
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()) {
return get(it->second);
} else {
throw ValueError(
format(
"key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",
key.c_str()
)
);
}
std::map<std::string, std::size_t>::iterator it;
// Try to find it
it = string_to_index_map.find(key);
// If it is found
if (it != string_to_index_map.end()) {
return get(it->second);
} else {
throw ValueError(
format(
"key [%s] was not found in string_to_index_map in JSONIncompressibleLibrary",
key.c_str()
)
);
}
};
/// Get a IncompressibleFluid instance stored in this library
@@ -557,16 +557,16 @@ IncompressibleFluid& JSONIncompressibleLibrary::get(std::string key) {
@param key The index of the fluid in the map
*/
IncompressibleFluid& JSONIncompressibleLibrary::get(std::size_t key) {
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()) {
return it->second;
} else {
throw ValueError(
format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
std::map<std::size_t, IncompressibleFluid>::iterator it;
// Try to find it
it = fluid_map.find(key);
// If it is found
if (it != fluid_map.end()) {
return it->second;
} else {
throw ValueError(
format("key [%d] was not found in JSONIncompressibleLibrary",key));
}
};
@@ -597,21 +597,21 @@ static JSONIncompressibleLibrary library;
void load_incompressible_library()
{
rapidjson::Document dd;
rapidjson::Document dd;
// This json formatted string comes from the all_incompressibles_JSON.h header which is a C++-escaped version of the JSON file
dd.Parse<0>(all_incompressibles_JSON.c_str());
if (dd.HasParseError()){
if (dd.HasParseError()){
throw ValueError("Unable to load all_incompressibles_JSON.json");
} else{
try{library.add_many(dd);}catch(std::exception &e){std::cout << e.what() << std::endl;}
}
// TODO: Implement LiBr in the source code!
//library.add_obj(LiBrSolution());
// TODO: Implement LiBr in the source code!
//library.add_obj(LiBrSolution());
}
JSONIncompressibleLibrary & get_incompressible_library(void){
if (library.is_empty()){ load_incompressible_library(); }
return library;
if (library.is_empty()){ load_incompressible_library(); }
return library;
}
IncompressibleFluid& get_incompressible_fluid(std::string fluid_string){

168
src/Backends/Incompressible/IncompressibleLibrary.h
View File

@@ -29,101 +29,101 @@ extern int get_debug_level();
class LiBrSolution : public IncompressibleFluid{
protected:
static double const M_H2O; /* kg/mol, molar mass of H2O */
static double const M_LiBr; /* kg/mol, molar mass of LiBr */
static double const T0; /* K, constant */
static double const M_H2O; /* kg/mol, molar mass of H2O */
static double const M_LiBr; /* kg/mol, molar mass of LiBr */
static double const T0; /* K, constant */
/* Critical point of H2O */
static double const Tc_H2O; /* K, temperature */
static double const pc_H2O; /* MPa, pressure */
static double const rhoc_H2O; /* mol/m^3 (322 kg/m^3), molar density */
static double const hc_H2O; /* J/mol, molar enthalpy */
static double const sc_H2O; /* J/(mol.K) molar entropy*/
/* Critical point of H2O */
static double const Tc_H2O; /* K, temperature */
static double const pc_H2O; /* MPa, pressure */
static double const rhoc_H2O; /* mol/m^3 (322 kg/m^3), molar density */
static double const hc_H2O; /* J/mol, molar enthalpy */
static double const sc_H2O; /* J/(mol.K) molar entropy*/
/*Triple point of H2O */
static double const Tt_H2O; /* K, temperature */
static double const cpt_H2O; /* J/(mol.K), molar isobaric heat capacity*/
/*Triple point of H2O */
static double const Tt_H2O; /* K, temperature */
static double const cpt_H2O; /* J/(mol.K), molar isobaric heat capacity*/
double ps_mix(double T, double x);
double rho_mix(double T, double x);
double cp_mix(double T, double x);
double h_mix(double T, double x);
double s_mix(double T, double x);
double ps_H2O(double T);
double rho_H2O(double T);
double cp_H2O(double T);
double h_H2O(double T);
double s_H2O(double T);
double ps_mix(double T, double x);
double rho_mix(double T, double x);
double cp_mix(double T, double x);
double h_mix(double T, double x);
double s_mix(double T, double x);
double ps_H2O(double T);
double rho_H2O(double T);
double cp_H2O(double T);
double h_H2O(double T);
double s_H2O(double T);
/** Finished with the code from the paper. Now we need to
* convert the molar values to mass-based units. */
double massToMole(double w);
double molarToSpecific(double w, double value);
/** Finished with the code from the paper. Now we need to
* convert the molar values to mass-based units. */
double massToMole(double w);
double molarToSpecific(double w, double value);
static const bool debug;
static const bool debug;
public:
LiBrSolution();
LiBrSolution();
double rho(double T, double p, double x);
double c(double T, double p, double x);
//double h(double T_K, double p, double x);
double s(double T, double p, double x);
double visc(double T, double p, double x);
double cond(double T, double p, double x);
double u(double T, double p, double x);
double psat(double T, double x);
double Tfreeze(double p, double x);
double rho(double T, double p, double x);
double c(double T, double p, double x);
//double h(double T_K, double p, double x);
double s(double T, double p, double x);
double visc(double T, double p, double x);
double cond(double T, double p, double x);
double u(double T, double p, double x);
double psat(double T, double x);
double Tfreeze(double p, double x);
/* Some functions can be inverted directly, those are listed
* here. It is also possible to solve for other quantities, but
* that involves some more sophisticated processing and is not
* done here, but in the backend, T(h,p) for example.
*/
/// Temperature as a function of density, pressure and composition.
double T_rho (double Dmass, double p, double x){throw NotImplementedError(format("%s (%d): T from density is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of heat capacities as a function of temperature, pressure and composition.
double T_c (double Cmass, double p, double x){throw NotImplementedError(format("%s (%d): T from heat capacity is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of entropy as a function of temperature, pressure and composition.
double T_s (double Smass, double p, double x){throw NotImplementedError(format("%s (%d): T from entropy is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of internal energy as a function of temperature, pressure and composition.
double T_u (double Umass, double p, double x){throw NotImplementedError(format("%s (%d): T from internal energy is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of enthalpy, pressure and composition.
//double T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));}
/// Viscosity as a function of temperature, pressure and composition.
double T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));}
/// Thermal conductivity as a function of temperature, pressure and composition.
double T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));}
/// Saturation pressure as a function of temperature and composition.
double T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));}
/// Composition as a function of freezing temperature and pressure.
double x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));}
/* Some functions can be inverted directly, those are listed
* here. It is also possible to solve for other quantities, but
* that involves some more sophisticated processing and is not
* done here, but in the backend, T(h,p) for example.
*/
/// Temperature as a function of density, pressure and composition.
double T_rho (double Dmass, double p, double x){throw NotImplementedError(format("%s (%d): T from density is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of heat capacities as a function of temperature, pressure and composition.
double T_c (double Cmass, double p, double x){throw NotImplementedError(format("%s (%d): T from heat capacity is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of entropy as a function of temperature, pressure and composition.
double T_s (double Smass, double p, double x){throw NotImplementedError(format("%s (%d): T from entropy is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of internal energy as a function of temperature, pressure and composition.
double T_u (double Umass, double p, double x){throw NotImplementedError(format("%s (%d): T from internal energy is not implemented for LiBr.",__FILE__,__LINE__));}
/// Temperature as a function of enthalpy, pressure and composition.
//double T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));}
/// Viscosity as a function of temperature, pressure and composition.
double T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));}
/// Thermal conductivity as a function of temperature, pressure and composition.
double T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));}
/// Saturation pressure as a function of temperature and composition.
double T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));}
/// Composition as a function of freezing temperature and pressure.
double x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));}
/// Overwrite some standard functions that cannot be used with LiBr
void setName(std::string name){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setDescription(std::string description){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setReference(std::string reference){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTmax(double Tmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTmin(double Tmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxmax(double xmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxmin(double xmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTminPsat(double TminPsat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
/// Overwrite some standard functions that cannot be used with LiBr
void setName(std::string name){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setDescription(std::string description){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setReference(std::string reference){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTmax(double Tmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTmin(double Tmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxmax(double xmax){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxmin(double xmin){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTminPsat(double TminPsat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTbase(double Tbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxbase(double xbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTbase(double Tbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setxbase(double xbase){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setDensity(IncompressibleData density){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setSpecificHeat(IncompressibleData specific_heat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setViscosity(IncompressibleData viscosity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setConductivity(IncompressibleData conductivity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setPsat(IncompressibleData p_sat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTfreeze(IncompressibleData T_freeze){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setVolToMass(IncompressibleData volToMass){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setMassToMole(IncompressibleData massToMole){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setDensity(IncompressibleData density){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setSpecificHeat(IncompressibleData specific_heat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setViscosity(IncompressibleData viscosity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setConductivity(IncompressibleData conductivity){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setPsat(IncompressibleData p_sat){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setTfreeze(IncompressibleData T_freeze){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setVolToMass(IncompressibleData volToMass){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
void setMassToMole(IncompressibleData massToMole){throw ValueError(format("%s (%d): Cannot change property of LiBr class",__FILE__,__LINE__));}
bool is_pure() {return false;};
bool is_pure() {return false;};
};
@@ -143,8 +143,8 @@ a rapidjson array of fluids to the add_many function.
class JSONIncompressibleLibrary
{
/// Map from CAS code to JSON instance.
/** This is not practical for the incompressibles, the CAS may not be
* defined for blends of heat transfer fluids and solutions.
/** This is not practical for the incompressibles, the CAS may not be
* defined for blends of heat transfer fluids and solutions.
*/
std::map<std::size_t, IncompressibleFluid> fluid_map;
std::vector<std::string> name_vector_pure, name_vector_solution;
@@ -183,8 +183,8 @@ public:
/// Return a comma-separated list of incompressible pure fluid names
std::string get_incompressible_list_pure(void){ return strjoin(name_vector_pure, ",");};
/// Return a comma-separated list of solution names
std::string get_incompressible_list_solution(void){ return strjoin(name_vector_solution, ",");};
/// Return a comma-separated list of solution names
std::string get_incompressible_list_solution(void){ return strjoin(name_vector_solution, ",");};
};
/// Get a reference to the library instance

20
src/Backends/REFPROP/REFPROPBackend.cpp
View File

@@ -23,21 +23,21 @@
namespace CoolProp {
REFPROPBackend::REFPROPBackend(const std::string & fluid_name) {
// Do the REFPROP instantiation for this fluid
// Do the REFPROP instantiation for this fluid
// Try to add this fluid to REFPROP - might want to think about making array of
// components and setting mole fractions if they change a lot.
std::vector<std::string> component_names(1,fluid_name);
set_REFPROP_fluids(component_names);
// Try to add this fluid to REFPROP - might want to think about making array of
// components and setting mole fractions if they change a lot.
std::vector<std::string> component_names(1,fluid_name);
set_REFPROP_fluids(component_names);
// Set the mole fraction to 1 in the base class (we can't set the mole fraction in this class,
// otherwise a NotImplementedError will be returned)
std::vector<long double> x(1, 1.0); // (one element with value of 1.0)
REFPROPMixtureBackend::set_mole_fractions(x);
// Set the mole fraction to 1 in the base class (we can't set the mole fraction in this class,
// otherwise a NotImplementedError will be returned)
std::vector<long double> x(1, 1.0); // (one element with value of 1.0)
REFPROPMixtureBackend::set_mole_fractions(x);
}
REFPROPBackend::~REFPROPBackend() {
// TODO Auto-generated destructor stub
// TODO Auto-generated destructor stub
}
} /* namespace CoolProp */

8
src/Backends/REFPROP/REFPROPBackend.h
View File

@@ -20,10 +20,10 @@ and exposes just the pure fluid interface.
class REFPROPBackend : public REFPROPMixtureBackend {
public:
REFPROPBackend();
REFPROPBackend(const std::string &fluid_name);
virtual ~REFPROPBackend();
REFPROPBackend();
REFPROPBackend(const std::string &fluid_name);
virtual ~REFPROPBackend();
};
} /* namespace CoolProp */

132
src/Backends/REFPROP/REFPROPMixtureBackend.cpp
View File

@@ -71,7 +71,7 @@ surface tension N/m
#define filepathlength 255
#define lengthofreference 3
#define errormessagelength 255
#define ncmax 20 // Note: ncmax is the max number of components
#define ncmax 20 // Note: ncmax is the max number of components
#define numparams 72
#define maxcoefs 50
@@ -82,7 +82,7 @@ std::string LoadedREFPROPRef;
#define filepathlength 255
#define lengthofreference 3
#define errormessagelength 255
#define ncmax 20 // Note: ncmax is the max number of components
#define ncmax 20 // Note: ncmax is the max number of components
#define numparams 72
#define maxcoefs 50
@@ -318,7 +318,7 @@ double setFunctionPointers()
SETMODdll = (SETMODdll_POINTER) getFunctionPointer((char *)SETMODdll_NAME);
SETREFdll = (SETREFdll_POINTER) getFunctionPointer((char *)SETREFdll_NAME);
SETUPdll = (SETUPdll_POINTER) getFunctionPointer((char *)SETUPdll_NAME);
// SPECGRdll = (SPECGRdll_POINTER) getFunctionPointer((char *)SPECGRdll_NAME); // not in library
// SPECGRdll = (SPECGRdll_POINTER) getFunctionPointer((char *)SPECGRdll_NAME); // not in library
SUBLPdll = (SUBLPdll_POINTER) getFunctionPointer((char *)SUBLPdll_NAME);
SUBLTdll = (SUBLTdll_POINTER) getFunctionPointer((char *)SUBLTdll_NAME);
SURFTdll = (SURFTdll_POINTER) getFunctionPointer((char *)SURFTdll_NAME);
@@ -595,81 +595,81 @@ void REFPROPMixtureBackend::check_status(void)
void REFPROPMixtureBackend::limits(double &Tmin, double &Tmax, double &rhomolarmax, double &pmax)
{
/*
*
subroutine LIMITS (htyp,x,tmin,tmax,Dmax,pmax)
c
c returns limits of a property model as a function of composition
c
c Pure fluid limits are read in from the .fld files; for mixtures, a
c simple mole fraction weighting in reduced variables is used.
c
c inputs:
c htyp--flag indicating which models are to be checked [character*3]
c 'EOS': equation of state for thermodynamic properties
c 'ETA': viscosity
c 'TCX': thermal conductivity
c 'STN': surface tension
c x--composition array [mol frac]
c outputs:
c tmin--minimum temperature for model specified by htyp [K]
c tmax--maximum temperature [K]
c Dmax--maximum density [mol/L]
c pmax--maximum pressure [kPa]
*
*/
double Dmax_mol_L,pmax_kPa;
/*
*
subroutine LIMITS (htyp,x,tmin,tmax,Dmax,pmax)
c
c returns limits of a property model as a function of composition
c
c Pure fluid limits are read in from the .fld files; for mixtures, a
c simple mole fraction weighting in reduced variables is used.
c
c inputs:
c htyp--flag indicating which models are to be checked [character*3]
c 'EOS': equation of state for thermodynamic properties
c 'ETA': viscosity
c 'TCX': thermal conductivity
c 'STN': surface tension
c x--composition array [mol frac]
c outputs:
c tmin--minimum temperature for model specified by htyp [K]
c tmax--maximum temperature [K]
c Dmax--maximum density [mol/L]
c pmax--maximum pressure [kPa]
*
*/
double Dmax_mol_L,pmax_kPa;
char htyp[] = "EOS";
LIMITSdll(htyp, &(mole_fractions[0]), &Tmin, &Tmax, &Dmax_mol_L, &pmax_kPa, 3);
pmax = pmax_kPa*1000;
rhomolarmax = Dmax_mol_L*1000;
LIMITSdll(htyp, &(mole_fractions[0]), &Tmin, &Tmax, &Dmax_mol_L, &pmax_kPa, 3);
pmax = pmax_kPa*1000;
rhomolarmax = Dmax_mol_L*1000;
}
long double REFPROPMixtureBackend::calc_pmax(void){
double Tmin, Tmax, rhomolarmax, pmax;
limits(Tmin, Tmax, rhomolarmax, pmax);
return static_cast<long double>(pmax);
double Tmin, Tmax, rhomolarmax, pmax;
limits(Tmin, Tmax, rhomolarmax, pmax);
return static_cast<long double>(pmax);
};
long double REFPROPMixtureBackend::calc_Tmax(void){
double Tmin, Tmax, rhomolarmax, pmax;
limits(Tmin, Tmax, rhomolarmax, pmax);
return static_cast<long double>(Tmax);
double Tmin, Tmax, rhomolarmax, pmax;
limits(Tmin, Tmax, rhomolarmax, pmax);
return static_cast<long double>(Tmax);
};
long double REFPROPMixtureBackend::calc_T_critical(){
long ierr;
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(Tcrit);
return static_cast<long double>(Tcrit);
};
long double REFPROPMixtureBackend::calc_p_critical(){
long ierr;
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(pcrit_kPa*1000);
return static_cast<long double>(pcrit_kPa*1000);
};
long double REFPROPMixtureBackend::calc_rhomolar_critical(){
long ierr;
char herr[255];
double Tcrit, pcrit_kPa, dcrit_mol_L;
CRITPdll(&(mole_fractions[0]),&Tcrit,&pcrit_kPa,&dcrit_mol_L,&ierr,herr,255); if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return static_cast<long double>(dcrit_mol_L*1000);
return static_cast<long double>(dcrit_mol_L*1000);
};
long double REFPROPMixtureBackend::calc_Ttriple(){
if (mole_fractions.size() != 1){throw ValueError("calc_Ttriple cannot be evaluated for mixtures");}
long icomp = 0;
double wmm, ttrp, tnbpt, tc, pc, Dc, Zc, acf, dip, Rgas;
INFOdll(&icomp, &wmm, &ttrp, &tnbpt, &tc, &pc, &Dc, &Zc, &acf, &dip, &Rgas);
return static_cast<long double>(ttrp);
return static_cast<long double>(ttrp);
};
long double REFPROPMixtureBackend::calc_molar_mass(void)
{
double wmm_kg_kmol;
WMOLdll(&(mole_fractions[0]), &wmm_kg_kmol); // returns mole mass in kg/kmol
_molar_mass = wmm_kg_kmol/1000; // kg/mol
return static_cast<long double>(_molar_mass.pt());
return static_cast<long double>(_molar_mass.pt());
};
double REFPROPMixtureBackend::calc_melt_Tmax()
{
long ierr;
@@ -687,33 +687,33 @@ double REFPROPMixtureBackend::calc_melt_Tmax()
}
long double REFPROPMixtureBackend::calc_melting_line(int param, int given, long double value)
{
long ierr;
long ierr;
char herr[255];
if (param == iP && given == iT){
double _T = static_cast<double>(value), p_kPa;
MELTTdll(&_T, &(mole_fractions[0]),
if (param == iP && given == iT){
double _T = static_cast<double>(value), p_kPa;
MELTTdll(&_T, &(mole_fractions[0]),
&p_kPa,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else if (param == iT && given == iP){
double p_kPa = static_cast<double>(value), _T;
MELTPdll(&p_kPa, &(mole_fractions[0]),
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else if (param == iT && given == iP){
double p_kPa = static_cast<double>(value), _T;
MELTPdll(&p_kPa, &(mole_fractions[0]),
&_T,
&ierr,herr,errormessagelength); // Error message
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else{
throw ValueError(format("calc_melting_line(%s,%s,%Lg) is an invalid set of inputs ",
get_parameter_information(param,"short").c_str(),
get_parameter_information(given,"short").c_str(),
value
)
);
}
if (ierr > 0) { throw ValueError(format("%s",herr).c_str()); } //else if (ierr < 0) {set_warning(format("%s",herr).c_str());}
return p_kPa*1000;
}
else{
throw ValueError(format("calc_melting_line(%s,%s,%Lg) is an invalid set of inputs ",
get_parameter_information(param,"short").c_str(),
get_parameter_information(given,"short").c_str(),
value
)
);
}
}
@@ -1294,13 +1294,13 @@ void REFPROPMixtureBackend::update(CoolProp::input_pairs input_pair, double valu
// Use flash routine to find properties
TQFLSHdll(&_T,&_Q,&(mole_fractions[0]),&kq,&p_kPa,&rho_mol_L,
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&rhoLmol_L,&rhoVmol_L,&(mole_fractions_liq[0]),&(mole_fractions_vap[0]), // Saturation terms
&emol,&hmol,&smol,&cvmol,&cpmol,&w, // Other thermodynamic terms
&ierr,herr,errormessagelength); // Error terms
if (ierr > 0) {
throw ValueError(format("TQ(%s): %s",LoadedREFPROPRef.c_str(), herr).c_str());
}// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());
throw ValueError(format("TQ(%s): %s",LoadedREFPROPRef.c_str(), herr).c_str());
}// TODO: else if (ierr < 0) {set_warning(format("%s",herr).c_str());
// Set all cache values that can be set with unit conversion to SI
_p = p_kPa*1000; // 1000 for conversion from kPa to Pa

114
src/Backends/REFPROP/REFPROPMixtureBackend.h
View File

@@ -18,96 +18,96 @@ namespace CoolProp {
class REFPROPMixtureBackend : public AbstractState {
protected:
std::size_t Ncomp;
bool _mole_fractions_set;
static bool _REFPROP_supported;
std::vector<double> mole_fractions, mass_fractions;
std::vector<double> mole_fractions_liq, mole_fractions_vap;
bool _mole_fractions_set;
static bool _REFPROP_supported;
std::vector<double> mole_fractions, mass_fractions;
std::vector<double> mole_fractions_liq, mole_fractions_vap;
public:
REFPROPMixtureBackend(){};
REFPROPMixtureBackend(){};
/// The instantiator
/// @param fluid_names The vector of strings of the fluid components, without file ending
REFPROPMixtureBackend(const std::vector<std::string>& fluid_names);
virtual ~REFPROPMixtureBackend();
/// The instantiator
/// @param fluid_names The vector of strings of the fluid components, without file ending
REFPROPMixtureBackend(const std::vector<std::string>& fluid_names);
virtual ~REFPROPMixtureBackend();
// REFPROP backend uses mole fractions
bool using_mole_fractions(){return true;}
bool using_mass_fractions(){return false;}
bool using_volu_fractions(){return false;}
/// Updating function for REFPROP
/**
In this function we take a pair of thermodynamic states, those defined in the input_pairs
enumeration and update all the internal variables that we can. REFPROP calculates
a lot of other state variables each time you use a flash routine so we cache all the
outputs that we can, which saves on computational time.
/// Updating function for REFPROP
/**
In this function we take a pair of thermodynamic states, those defined in the input_pairs
enumeration and update all the internal variables that we can. REFPROP calculates
a lot of other state variables each time you use a flash routine so we cache all the
outputs that we can, which saves on computational time.
@param input_pair Integer key from CoolProp::input_pairs to the two inputs that will be passed to the function
@param value1 First input value
@param value2 Second input value
*/
void update(CoolProp::input_pairs,
double value1,
double value2
);
@param input_pair Integer key from CoolProp::input_pairs to the two inputs that will be passed to the function
@param value1 First input value
@param value2 Second input value
*/
void update(CoolProp::input_pairs,
double value1,
double value2
);
long double calc_molar_mass(void);
/// Returns true if REFPROP is supported on this platform
bool REFPROP_supported(void);
/// Returns true if REFPROP is supported on this platform
bool REFPROP_supported(void);
long double calc_cpmolar_idealgas(void);
long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant);
/// Set the fluids in REFPROP DLL by calling the SETUPdll function
/**
@param fluid_names The vector of strings of the fluid components, without file ending
*/
void set_REFPROP_fluids(const std::vector<std::string> &fluid_names);
/// Set the fluids in REFPROP DLL by calling the SETUPdll function
/**
@param fluid_names The vector of strings of the fluid components, without file ending
*/
void set_REFPROP_fluids(const std::vector<std::string> &fluid_names);
/// Set the mole fractions
/**
@param mole_fractions The vector of mole fractions of the components
*/
void set_mole_fractions(const std::vector<long double> &mole_fractions);
/// Set the mole fractions
/**
@param mole_fractions The vector of mole fractions of the components
*/
void set_mole_fractions(const std::vector<long double> &mole_fractions);
/// Set the mass fractions
/**
@param mass_fractions The vector of mass fractions of the components
*/
void set_mass_fractions(const std::vector<long double> &mass_fractions);
/// Set the mass fractions
/**
@param mass_fractions The vector of mass fractions of the components
*/
void set_mass_fractions(const std::vector<long double> &mass_fractions);
void calc_phase_envelope(const std::string &type);
void calc_phase_envelope(const std::string &type);
std::vector<long double> calc_mole_fractions_liquid(void){return std::vector<long double>(mole_fractions_liq.begin(), mole_fractions_liq.end());}
std::vector<long double> calc_mole_fractions_vapor(void){return std::vector<long double>(mole_fractions_vap.begin(), mole_fractions_vap.end());}
/// Check if the mole fractions have been set, etc.
void check_status();
/// Check if the mole fractions have been set, etc.
void check_status();
/// Get the viscosity [Pa-s] (based on the temperature and density in the state class)
long double calc_viscosity(void);
/// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class)
long double calc_conductivity(void);
/// Get the surface tension [N/m] (based on the temperature in the state class). Invalid for temperatures above critical point or below triple point temperature
long double calc_surface_tension(void);
/// Get the viscosity [Pa-s] (based on the temperature and density in the state class)
long double calc_viscosity(void);
/// Get the thermal conductivity [W/m/K] (based on the temperature and density in the state class)
long double calc_conductivity(void);
/// Get the surface tension [N/m] (based on the temperature in the state class). Invalid for temperatures above critical point or below triple point temperature
long double calc_surface_tension(void);
long double calc_fugacity_coefficient(int i);
long double calc_melting_line(int param, int given, long double value);
long double calc_melting_line(int param, int given, long double value);
bool has_melting_line(){return true;};
double calc_melt_Tmax();
long double calc_T_critical(void);
long double calc_p_critical(void);
long double calc_rhomolar_critical(void);
long double calc_Ttriple(void);
/// A wrapper function to calculate the limits for the EOS
void limits(double &Tmin, double &Tmax, double &rhomolarmax, double &pmax);
/// Calculate the maximum pressure
long double calc_pmax(void);
/// Calculate the maximum temperature
long double calc_Tmax(void);
/// A wrapper function to calculate the limits for the EOS
void limits(double &Tmin, double &Tmax, double &rhomolarmax, double &pmax);
/// Calculate the maximum pressure
long double calc_pmax(void);
/// Calculate the maximum temperature
long double calc_Tmax(void);
};
} /* namespace CoolProp */