From 864df6d590eaa43cddeaedda6e2c54c9f6bae706 Mon Sep 17 00:00:00 2001 From: jowr Date: Fri, 20 Jun 2014 17:02:44 +0200 Subject: [PATCH] Double-checked entropy equations --- include/PolyMath.h | 3 - .../Helmholtz/Fluids/Incompressible.cpp | 250 +++++++++--------- .../Helmholtz/Fluids/Incompressible.h | 80 +++--- 3 files changed, 156 insertions(+), 177 deletions(-) diff --git a/include/PolyMath.h b/include/PolyMath.h index 697b7a55..f407c691 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -8,9 +8,6 @@ #include #include #include "Solvers.h" -//#include // inner_product -//#include -//#include "float.h" #include "MatrixMath.h" #include diff --git a/src/Backends/Helmholtz/Fluids/Incompressible.cpp b/src/Backends/Helmholtz/Fluids/Incompressible.cpp index 5dc2cb56..d176a7ba 100644 --- a/src/Backends/Helmholtz/Fluids/Incompressible.cpp +++ b/src/Backends/Helmholtz/Fluids/Incompressible.cpp @@ -1,130 +1,120 @@ -/* - * Incompressible.cpp - * - * Created on: 20 Dec 2013 - * Author: jowr - */ - -#include "Incompressible.h" - - -namespace CoolProp { - -std::vector > Incompressible::changeAxis(const std::vector &input){ - std::vector > output; - std::vector tmp; - size_t sizeX = input.size(); - for (size_t i = 0; i < sizeX; ++i){ - tmp.clear(); - tmp.push_back(input[i]); - output.push_back(tmp); - } - return output; -} - -// /* All functions need T and p as input. Might not -// * be necessary, but gives a clearer structure. -// */ -// /// Density as a function of temperature, pressure and composition. -// double Incompressible::rho(double T_K, double p) { -// return poly.polyval(cRho, getxInput(x), getTInput(T_K)); -// } -// /// Heat capacities as a function of temperature, pressure and composition. -// double Incompressible::c(double T_K, double p) { -// return poly.polyval(cHeat, getxInput(x), getTInput(T_K)); -// } -// /// Enthalpy as a function of temperature, pressure and composition. -// double Incompressible::h(double T_K, double p) { -// return h_u(T_K, p); -// } -// /// Entropy as a function of temperature, pressure and composition. -// double Incompressible::s(double T_K, double p) { -// return poly.polyfracintcentral(cHeat, getxInput(x), T_K, Tbase) -// - poly.polyfracintcentral(cHeat, getxInput(x), Tref, Tbase); -// } -// /// Viscosity as a function of temperature, pressure and composition. -// double Incompressible::visc(double T_K, double p) { -// return expo.expval(cVisc, getxInput(x), getTInput(T_K), 2) / 1e3; -// } -// /// Thermal conductivity as a function of temperature, pressure and composition. -// double Incompressible::cond(double T_K, double p) { -// return poly.polyval(cCond, getxInput(x), getTInput(T_K)); -// } -// /// Internal energy as a function of temperature, pressure and composition. -// double Incompressible::u(double T_K, double p) { -// return poly.polyint(cHeat, getxInput(x), getTInput(T_K)) -// - poly.polyint(cHeat, getxInput(x), getTInput(Tref)); -// } - -/// Saturation pressure as a function of temperature and composition. -double Incompressible::psat(double T_K ){throw NotImplementedError("Psat is not available");}; -/// Freezing temperature as a function of pressure and composition. -double Incompressible::Tfreeze( double p){throw NotImplementedError("Tfreeze is not available");}; - - -/* - * Some more functions to provide a single implementation - * of important routines. - * We start with the check functions that can validate input - * in terms of pressure p, temperature T and composition x. - */ - -/// Check validity of temperature input. -/** Compares the given temperature T to the result of a - * freezing point calculation. This is not necessarily - * defined for all fluids, default values do not - * cause errors. */ -bool Incompressible::checkT(double T_K, double p){ - if( Tmin < 0. ) { - throw ValueError("Please specify the minimum temperature."); - } else if( Tmax < 0.) { - throw ValueError("Please specify the maximum temperature."); - } else if ( (Tmin>T_K) || (T_K>Tmax) ) { - throw ValueError(format("Your temperature %f is not between %f and %f.",T_K,Tmin,Tmax)); - } else if (T_K < Tfreeze(p)) { - throw ValueError(format("Your temperature %f is below the freezing point of %f.",T_K,Tfreeze(p))); - } else { - return true; - } - return false; -} - -/// Check validity of pressure input. -/** Compares the given pressure p to the saturation pressure at - * temperature T and throws and exception if p is lower than - * the saturation conditions. - * The default value for psat is -1 yielding true if psat - * is not redefined in the subclass. - * */ -bool Incompressible::checkP(double T_K, double p) { - double ps = psat(T_K); - if (px) || (x>xmax) ) { - throw ValueError(format("Your composition %f is not between %f and %f.",x,xmin,xmax)); - } else { - return true; - } - return false; -} - -/// Check validity of temperature, pressure and composition input. -bool Incompressible::checkTPX(double T, double p, double x) { - return (checkT(T,p) && checkP(T,p) && checkX(x)); -} - -} /* namespace CoolProp */ +///* +// * Incompressible.cpp +// * +// * Created on: 20 Dec 2013 +// * Author: jowr +// */ +// +//#include "Incompressible.h" +// +// +//namespace CoolProp { +// +// +// +//// /* All functions need T and p as input. Might not +//// * be necessary, but gives a clearer structure. +//// */ +//// /// Density as a function of temperature, pressure and composition. +//// double Incompressible::rho(double T_K, double p) { +//// return poly.polyval(cRho, getxInput(x), getTInput(T_K)); +//// } +//// /// Heat capacities as a function of temperature, pressure and composition. +//// double Incompressible::c(double T_K, double p) { +//// return poly.polyval(cHeat, getxInput(x), getTInput(T_K)); +//// } +//// /// Enthalpy as a function of temperature, pressure and composition. +//// double Incompressible::h(double T_K, double p) { +//// return h_u(T_K, p); +//// } +//// /// Entropy as a function of temperature, pressure and composition. +//// double Incompressible::s(double T_K, double p) { +//// return poly.polyfracintcentral(cHeat, getxInput(x), T_K, Tbase) +//// - poly.polyfracintcentral(cHeat, getxInput(x), Tref, Tbase); +//// } +//// /// Viscosity as a function of temperature, pressure and composition. +//// double Incompressible::visc(double T_K, double p) { +//// return expo.expval(cVisc, getxInput(x), getTInput(T_K), 2) / 1e3; +//// } +//// /// Thermal conductivity as a function of temperature, pressure and composition. +//// double Incompressible::cond(double T_K, double p) { +//// return poly.polyval(cCond, getxInput(x), getTInput(T_K)); +//// } +//// /// Internal energy as a function of temperature, pressure and composition. +//// double Incompressible::u(double T_K, double p) { +//// return poly.polyint(cHeat, getxInput(x), getTInput(T_K)) +//// - poly.polyint(cHeat, getxInput(x), getTInput(Tref)); +//// } +// +///// Saturation pressure as a function of temperature and composition. +//double Incompressible::psat(double T_K ){throw NotImplementedError("Psat is not available");}; +///// Freezing temperature as a function of pressure and composition. +//double Incompressible::Tfreeze( double p){throw NotImplementedError("Tfreeze is not available");}; +// +// +///* +// * Some more functions to provide a single implementation +// * of important routines. +// * We start with the check functions that can validate input +// * in terms of pressure p, temperature T and composition x. +// */ +// +///// Check validity of temperature input. +///** Compares the given temperature T to the result of a +// * freezing point calculation. This is not necessarily +// * defined for all fluids, default values do not +// * cause errors. */ +//bool Incompressible::checkT(double T_K, double p){ +// if( Tmin < 0. ) { +// throw ValueError("Please specify the minimum temperature."); +// } else if( Tmax < 0.) { +// throw ValueError("Please specify the maximum temperature."); +// } else if ( (Tmin>T_K) || (T_K>Tmax) ) { +// throw ValueError(format("Your temperature %f is not between %f and %f.",T_K,Tmin,Tmax)); +// } else if (T_K < Tfreeze(p)) { +// throw ValueError(format("Your temperature %f is below the freezing point of %f.",T_K,Tfreeze(p))); +// } else { +// return true; +// } +// return false; +//} +// +///// Check validity of pressure input. +///** Compares the given pressure p to the saturation pressure at +// * temperature T and throws and exception if p is lower than +// * the saturation conditions. +// * The default value for psat is -1 yielding true if psat +// * is not redefined in the subclass. +// * */ +//bool Incompressible::checkP(double T_K, double p) { +// double ps = psat(T_K); +// if (px) || (x>xmax) ) { +// throw ValueError(format("Your composition %f is not between %f and %f.",x,xmin,xmax)); +// } else { +// return true; +// } +// return false; +//} +// +///// Check validity of temperature, pressure and composition input. +//bool Incompressible::checkTPX(double T, double p, double x) { +// return (checkT(T,p) && checkP(T,p) && checkX(x)); +//} +// +//} /* namespace CoolProp */ diff --git a/src/Backends/Helmholtz/Fluids/Incompressible.h b/src/Backends/Helmholtz/Fluids/Incompressible.h index 5b4a55ca..675228ef 100644 --- a/src/Backends/Helmholtz/Fluids/Incompressible.h +++ b/src/Backends/Helmholtz/Fluids/Incompressible.h @@ -8,7 +8,9 @@ #ifndef INCOMPRESSIBLE_H_ #define INCOMPRESSIBLE_H_ +#include #include "PolyMath.h" +#include "MatrixMath.h" namespace CoolProp { @@ -25,20 +27,16 @@ protected: double TminPsat; double xref, Tref; double xbase, Tbase; - double x; - std::vector< std::vector > cRho; - std::vector< std::vector > cHeat; - std::vector< std::vector > cVisc; - std::vector< std::vector > cCond; - std::vector< std::vector > cPsat; - std::vector cTfreeze; + Eigen::MatrixXd cRho; + Eigen::MatrixXd cHeat; + Eigen::MatrixXd cVisc; + Eigen::MatrixXd cCond; + Eigen::MatrixXd cPsat; + Eigen::MatrixXd cTfreeze; + Eigen::MatrixXd cV2M; - std::vector > changeAxis(const std::vector &input); - - //BasePolynomial poly; - //PolynomialSolver solver; - //BaseExponential expo; + Polynomial2DFrac poly; public: Incompressible(); @@ -58,12 +56,10 @@ public: double getxref() const {return xref;} double getTbase() const {return Tbase;} double getxbase() const {return xbase;} - double getx() const {return x;} void setName(std::string name) {this->name = name;} void setDescription(std::string description) {this->description = description;} void setReference(std::string reference) {this->reference = reference;} - void setTmax(double Tmax) {this->Tmax = Tmax;} void setTmin(double Tmin) {this->Tmin = Tmin;} void setxmax(double xmax) {this->xmax = xmax;} @@ -73,21 +69,15 @@ public: void setxref(double xref) {this->xref = xref;} void setTbase(double Tbase) {this->Tbase = Tbase;} void setxbase(double xbase) {this->xbase = xbase;} - void setx(double x) {this->x = x;} // Setters for the coefficients - void setcRho(std::vector > cRho){this->cRho = cRho;} - void setcHeat(std::vector > cHeat){this->cHeat = cHeat;} - void setcVisc(std::vector > cVisc){this->cVisc = cVisc;} - void setcCond(std::vector > cCond){this->cCond = cCond;} - void setcPsat(std::vector > cPsat){this->cPsat = cPsat;} - void setcTfreeze(std::vector cTfreeze){this->cTfreeze = cTfreeze;} - - void setcRho(std::vector cRho){this->cRho = changeAxis(cRho);} - void setcHeat(std::vector cHeat){this->cHeat = changeAxis(cHeat);} - void setcVisc(std::vector cVisc){this->cVisc = changeAxis(cVisc);} - void setcCond(std::vector cCond){this->cCond = changeAxis(cCond);} - void setcPsat(std::vector cPsat){this->cPsat = changeAxis(cPsat);} + void setcRho(Eigen::MatrixXd cRho){this->cRho = cRho;} + void setcHeat(Eigen::MatrixXd cHeat){this->cHeat = cHeat;} + void setcVisc(Eigen::MatrixXd cVisc){this->cVisc = cVisc;} + void setcCond(Eigen::MatrixXd cCond){this->cCond = cCond;} + void setcPsat(Eigen::MatrixXd cPsat){this->cPsat = cPsat;} + void setcTfreeze(Eigen::MatrixXd cTfreeze){this->cTfreeze = cTfreeze;} + void setcV2M(Eigen::MatrixXd cV2M){this->cV2M = cV2M;} double getTInput(double curTValue){return curTValue-Tbase;} double getxInput(double curxValue){return (curxValue-xbase)*100.0;} @@ -96,25 +86,27 @@ public: * be necessary, but gives a clearer structure. */ /// Density as a function of temperature, pressure and composition. - virtual double rho (double T_K, double p); + virtual double rho (double T_K, double p, double x); /// Heat capacities as a function of temperature, pressure and composition. - virtual double c (double T_K, double p); - virtual double cp (double T_K, double p){return c(T_K,p);}; - virtual double cv (double T_K, double p){return c(T_K,p);}; + virtual double c (double T_K, double p, double x); + virtual double cp (double T_K, double p, double x){return c(T_K,p,x);}; + virtual double cv (double T_K, double p, double x){return c(T_K,p,x);}; /// Entropy as a function of temperature, pressure and composition. - virtual double s (double T_K, double p); + virtual double s (double T_K, double p, double x); /// Internal energy as a function of temperature, pressure and composition. - virtual double u (double T_K, double p); + virtual double u (double T_K, double p, double x); /// Enthalpy as a function of temperature, pressure and composition. - virtual double h (double T_K, double p); + virtual double h (double T_K, double p, double x); /// Viscosity as a function of temperature, pressure and composition. - virtual double visc(double T_K, double p); + virtual double visc(double T_K, double p, double x); /// Thermal conductivity as a function of temperature, pressure and composition. - virtual double cond(double T_K, double p); + virtual double cond(double T_K, double p, double x); /// Saturation pressure as a function of temperature and composition. - virtual double psat(double T_K ); + virtual double psat(double T_K, double x ); /// Freezing temperature as a function of pressure and composition. - virtual double Tfreeze( double p); + virtual double Tfreeze( double p, double x); + /// Conversion from volume-based to mass-based composition. + virtual double V2M( double x); protected: @@ -127,16 +119,16 @@ protected: /** Calculate enthalpy as a function of temperature and * pressure employing functions for internal energy and * density. Provides consistent formulations. */ - double h_u(double T_K, double p) { - return u(T_K,p)+p/rho(T_K,p); + double h_u(double T_K, double p, double x) { + return u(T_K,p,x)+p/rho(T_K,p,x); }; /// Internal energy from h, p and rho. /** Calculate internal energy as a function of temperature * and pressure employing functions for enthalpy and * density. Provides consistent formulations. */ - double u_h(double T_K, double p) { - return h(T_K,p)-p/rho(T_K,p); + double u_h(double T_K, double p, double x) { + return h(T_K,p,x)-p/rho(T_K,p,x); }; @@ -150,7 +142,7 @@ protected: /** Compares the given temperature T to the result of a * freezing point calculation. This is not necessarily * defined for all fluids, default values do not cause errors. */ - bool checkT(double T_K, double p); + bool checkT(double T_K, double p, double x); /// Check validity of pressure input. /** Compares the given pressure p to the saturation pressure at @@ -159,7 +151,7 @@ protected: * The default value for psat is -1 yielding true if psat * is not redefined in the subclass. * */ - bool checkP(double T_K, double p); + bool checkP(double T_K, double p, double x); /// Check validity of composition input. /** Compares the given composition x to a stored minimum and