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

Merge branch 'master' of https://github.com/coolprop/coolprop

Browse Source
This commit is contained in:
Ian Bell
2015-02-10 18:04:15 -07:00
parent 8ae3205f68 b645f028b3
commit 102e2a8cb4
13 changed files with 291 additions and 177 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -90,11 +90,11 @@ void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double valu
_T = this->DmassP_flash(value1,value2);
break;
}
case PUmass_INPUTS: {
_p = value1;
_T = this->PUmass_flash(value1, value2);
break;
}
// case PUmass_INPUTS: {
// _p = value1;
// _T = this->PUmass_flash(value1, value2);
// break;
// }
case PSmass_INPUTS: {
_p = value1;
_T = this->PSmass_flash(value1, value2);
@@ -264,7 +264,35 @@ long double IncompressibleBackend::HmassP_flash(long double hmass, long double p
@returns T The temperature in K
*/
long double IncompressibleBackend::PSmass_flash(long double p, long double smass){
return fluid->T_s(smass, p, _fractions[0]);
class PSmass_residual : public FuncWrapper1D {
protected:
double p,x,s_in;
IncompressibleFluid* fluid;
protected:
PSmass_residual(){};
public:
PSmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &s_in){
this->p = p;
this->x = x;
this->s_in = s_in;
this->fluid = fluid;
}
virtual ~PSmass_residual(){};
double call(double target){
return fluid->s(target,p,x) - s_in;
}
};
PSmass_residual res = PSmass_residual(fluid, p, _fractions[0], smass);
std::string errstring;
double macheps = DBL_EPSILON;
double tol = DBL_EPSILON*1e3;
int maxiter = 10;
double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring);
//if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl;
return result;
}
/// Calculate T given pressure and internal energy
@@ -274,7 +302,35 @@ long double IncompressibleBackend::PSmass_flash(long double p, long double smass
@returns T The temperature in K
*/
long double IncompressibleBackend::PUmass_flash(long double p, long double umass){
return fluid->T_u(umass, p, _fractions[0]);
class PUmass_residual : public FuncWrapper1D {
protected:
double p,x,u_in;
IncompressibleFluid* fluid;
protected:
PUmass_residual(){};
public:
PUmass_residual(IncompressibleFluid* fluid, const double &p, const double &x, const double &u_in){
this->p = p;
this->x = x;
this->u_in = u_in;
this->fluid = fluid;
}
virtual ~PUmass_residual(){};
double call(double target){
return fluid->u(target,p,x) - u_in;
}
};
PUmass_residual res = PUmass_residual(fluid, p, _fractions[0], umass);
std::string errstring;
double macheps = DBL_EPSILON;
double tol = DBL_EPSILON*1e3;
int maxiter = 10;
double result = Brent(&res, fluid->getTmin(), fluid->getTmax(), macheps, tol, maxiter, errstring);
//if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl;
return result;
}
@@ -436,7 +492,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
CAPTURE(res);
CHECK( check_abs(val,res,acc) );
}
// Make sure the ruslt does not change -> reference state...
// Make sure the result does not change -> reference state...
val = backend.calc_smass();
backend.update( CoolProp::PT_INPUTS, p, T );
res = backend.calc_smass();

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

@@ -19,23 +19,11 @@ and transport properties.
//IncompressibleFluid::IncompressibleFluid();
void IncompressibleFluid::set_reference_state(double T0, double p0, double x0, double h0, double s0){
this->_href = h0;
this->_pref = p0;
this->_Tref = T0;
this->_xref = x0;
this->_rhoref = rho(T0,p0,x0);
this->_uref = h0; // Reset the old reference value
this->_uref = u(T0,p0,x0); // uses _uref
this->_sref = s0; // Reset the old reference value
this->_sref = s(T0,p0,x0); // uses _sref
// Save the values for later calls at composition changes
this->href = h0;
this->pref = p0;
this->Tref = T0;
this->xref = x0;
this->rhoref = this->_rhoref;
this->uref = this->_uref;
this->sref = s0;
this->Tref = T0;
this->pref = p0;
this->xref = x0;
this->href = h0;
this->sref = s0;
}
void IncompressibleFluid::validate(){
@@ -88,6 +76,7 @@ double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, do
return _HUGE;
}
/// Density as a function of temperature, pressure and composition.
double IncompressibleFluid::rho (double T, double p, double x){
switch (density.type) {
@@ -135,10 +124,15 @@ double IncompressibleFluid::c (double T, double p, double x){
/// Entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::s (double T, double p, double x){
double dsdTatp = _HUGE;
double dsdpatT = _HUGE;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - _sref;
//TODO: Optimise this, add a cached value or direct calculation of definite integral
dsdTatp = poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, -1, 0, Tbase, xbase);
dsdpatT = this->rho(T,p,x);
dsdpatT = this->drhodTatPx(T,p,x)/dsdpatT/dsdpatT * (p-pref);
return sref + dsdTatp + dsdpatT;
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
@@ -151,24 +145,28 @@ double IncompressibleFluid::s (double T, double p, double x){
}
/// Internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::u (double T, double p, double x){
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//throw NotImplementedError("Here you should implement the polynomial.");
return poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - _uref;
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for internal energy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
double IncompressibleFluid::u (double T, double p, double x){return u_h(T,p,x);};
/// Enthalpy as a function of temperature, pressure and composition.
double IncompressibleFluid::h (double T, double p, double x){return h_u(T,p,x);};
double IncompressibleFluid::h (double T, double p, double x){
double dhdTatP = _HUGE;
double dhdpatT = _HUGE;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
//TODO: Optimise this, add a cached value or direct calculation of definite integral
dhdTatP = poly.integral(specific_heat.coeffs, T, x, 0, 0, 0, Tbase, xbase) - poly.integral(specific_heat.coeffs, Tref, xref, 0, 0, 0, Tbase, xbase);
dhdpatT = this->dhdpatTx(T, p, x) * (p-pref);
return href + dhdTatP + dhdpatT;
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for internal energy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
/// Viscosity as a function of temperature, pressure and composition.
double IncompressibleFluid::visc(double T, double p, double x){
@@ -280,6 +278,43 @@ double IncompressibleFluid::Tfreeze( double p, double x){
return _HUGE;
}
/* Below are direct calculations of the derivatives. Nothing
* special is going on, we simply use the polynomial class to
* derive the different functions with respect to temperature.
*/
/// Partial derivative of density with respect to temperature at constant pressure and composition
double IncompressibleFluid::drhodTatPx (double T, double p, double x){
switch (density.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.derivative(density.coeffs, T, x, 0, 0, 0, Tbase, xbase);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,density.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for density.",__FILE__,__LINE__,density.type));
break;
}
return _HUGE;
}
/* Other useful derivatives
*/
/// Partial derivative of enthalpy with respect to temperature at constant pressure and composition
double IncompressibleFluid::dhdpatTx (double T, double p, double x){
double rho = 0;
rho = this->rho(T,p,x);
return 1/rho * ( 1 + T/rho * this->drhodTatPx(T,p,x));
}
/// Partial derivative of entropy with respect to temperature at constant pressure and composition
double IncompressibleFluid::dsdpatTx (double T, double p, double x){
double rho = 0;
rho = this->rho(T,p,x);
return 1/rho/rho * this->drhodTatPx(T,p,x);
}
/// Mass fraction conversion function
/** If the fluid type is mass-based, it does not do anything. Otherwise,
* it converts the mass fraction to the required input. */
@@ -292,7 +327,7 @@ double IncompressibleFluid::inputFromMass (double T, double x){
throw NotImplementedError("Mass composition conversion has not been implemented.");
switch (mass2input.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
return poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
return baseExponential(mass2input, x, 0.0);
@@ -301,7 +336,7 @@ double IncompressibleFluid::inputFromMass (double T, double x){
return baseLogexponential(mass2input, x, 0.0);
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
return exp(poly.evaluate(mass2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
return basePolyOffset(mass2input, T, x);
@@ -329,7 +364,7 @@ double IncompressibleFluid::inputFromVolume (double T, double x){
throw NotImplementedError("Volume composition conversion has not been implemented.");
switch (volume2input.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
return poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
return baseExponential(volume2input, x, 0.0);
@@ -338,7 +373,7 @@ double IncompressibleFluid::inputFromVolume (double T, double x){
return baseLogexponential(volume2input, x, 0.0);
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
return exp(poly.evaluate(volume2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
return basePolyOffset(volume2input, T, x);
@@ -366,7 +401,7 @@ double IncompressibleFluid::inputFromMole (double T, double x){
throw NotImplementedError("Mole composition conversion has not been implemented.");
switch (mole2input.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase is defined in the correct way
return poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_EXPONENTIAL:
return baseExponential(mole2input, x, 0.0);
@@ -375,7 +410,7 @@ double IncompressibleFluid::inputFromMole (double T, double x){
return baseLogexponential(mole2input, x, 0.0);
break;
case IncompressibleData::INCOMPRESSIBLE_EXPPOLYNOMIAL:
return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase is defined in the correct way
return exp(poly.evaluate(mole2input.coeffs, T, x, 0, 0, 0.0, 0.0)); // TODO: make sure Tbase and xbase are defined in the correct way
break;
case IncompressibleData::INCOMPRESSIBLE_POLYOFFSET:
return basePolyOffset(mole2input, T, x);
@@ -428,49 +463,6 @@ double IncompressibleFluid::T_c (double Cmass, double p, double x){
}
return _HUGE;
}
/// Temperature as a function of entropy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_s (double Smass, double p, double x){
double s_raw = Smass + _sref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, s_raw, Tmin, Tmax, 0, -1, 0, Tbase, xbase, 0);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
/// Temperature as a function of internal energy as a function of temperature, pressure and composition.
double IncompressibleFluid::T_u (double Umass, double p, double x){
double u_raw = Umass + _uref;
switch (specific_heat.type) {
case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL:
return poly.solve_limitsInt(specific_heat.coeffs, x, u_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase, 0);
break;
case IncompressibleData::INCOMPRESSIBLE_NOT_SET:
throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type));
break;
default:
throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type));
break;
}
return _HUGE;
}
///// Temperature as a function of enthalpy, pressure and composition.
//double IncompressibleFluid::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 IncompressibleFluid::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 IncompressibleFluid::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 IncompressibleFluid::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 IncompressibleFluid::x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));}
/*
* Some more functions to provide a single implementation
@@ -737,7 +729,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare s
val = 145.59157247249246;
val = 144.08;
res = XLT.s(T,p,x);
{
CAPTURE(T);
@@ -756,7 +748,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare u
val = 45212.407309106304;
val = 44724.1;
res = XLT.u(T,p,x);
{
CAPTURE(T);
@@ -775,7 +767,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare h
val = 46388.7;
val = 45937;
res = XLT.h(T,p,x);
{
CAPTURE(T);
@@ -863,7 +855,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare s
expected = -206.62646783739274;
expected = -207.027;
actual = CH3OH.s(T,p,x);
{
CAPTURE(T);
@@ -904,7 +896,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare u
expected = -60043.78429641827;
expected = -60157.1;
actual = CH3OH.u(T,p,x);
{
CAPTURE(T);
@@ -927,7 +919,7 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi
}
// Compare h
expected = -58999.1;
expected = -59119;
actual = CH3OH.h(T,p,x);
{
CAPTURE(T);

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

@@ -288,8 +288,6 @@ LiBrSolution::LiBrSolution():IncompressibleFluid(){
pref = 0.0;
href = 0.0;
sref = 0.0;
uref = 0.0;
rhoref = 0.0;
xbase = 0.0;
Tbase = 0.0;

5
src/PolyMath.cpp
View File

@@ -596,7 +596,8 @@ double Polynomial2DFrac::derivative(const Eigen::MatrixXd &coefficients, const d
/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension
/// @param x_base double value that represents the base value for a centred fit in the 1st dimension
/// @param y_base double value that represents the base value for a centred fit in the 2nd dimension
double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){
/// @param ax_val double value that represents the base value for the definite integral on the chosen axis.
double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const double &ax_val){
Eigen::MatrixXd newCoefficients;
int int_exp,other_exp;
@@ -628,6 +629,8 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou
}
if (int_exp<-1) throw NotImplementedError(format("%s (%d): This function is only implemented for lowest exponents >= -1, int_exp=%d ",__FILE__,__LINE__,int_exp));
// TODO: Fix this and allow the direct calculation of integrals
if (std::abs(ax_val)>DBL_EPSILON) throw NotImplementedError(format("%s (%d): This function is only implemented for indefinite integrals, ax_val=%d ",__FILE__,__LINE__,ax_val));
size_t r = newCoefficients.rows();
size_t c = newCoefficients.cols();