From 792c466c3d2d6058f37248acab486b09e9eeda81 Mon Sep 17 00:00:00 2001 From: Jorrit Wronski Date: Tue, 27 Jan 2015 20:44:05 +0100 Subject: [PATCH] I hope this is the last commit to fix #429 and to fix #413 --- include/IncompressibleFluid.h | 8 +- include/PolyMath.h | 2 +- .../Incompressible/IncompressibleBackend.cpp | 48 +++++++++--- .../Incompressible/IncompressibleFluid.cpp | 78 ++++++++++++++----- 4 files changed, 105 insertions(+), 31 deletions(-) diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h index cbd3aa54..5445b91f 100644 --- a/include/IncompressibleFluid.h +++ b/include/IncompressibleFluid.h @@ -65,6 +65,10 @@ protected: double uref, rhoref; double xbase, Tbase; + double _xref, _Tref, _pref; + double _href, _sref; + double _uref, _rhoref; + /// These are the objects that hold the coefficients /** Note that all polynomials require a 2-dimensional array * of coefficients. This array may have only one row or @@ -272,7 +276,7 @@ protected: * pressure employing functions for internal energy and * density. Provides consistent formulations. */ double h_u(double T, double p, double x) { - return u(T,p,x)+(p-pref)/rhoref; + return u(T,p,x)+(p-_pref)/_rhoref; }; /// Internal energy from h, p and rho. @@ -280,7 +284,7 @@ protected: * and pressure employing functions for enthalpy and * density. Provides consistent formulations. */ double u_h(double T, double p, double x) { - return h(T,p,x)-(p-pref)/rhoref; + return h(T,p,x)-(p-_pref)/_rhoref; }; diff --git a/include/PolyMath.h b/include/PolyMath.h index 9ef18197..79fc500d 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -159,7 +159,7 @@ protected: double baseHorner(const std::vector &coefficients, double x); DEPRECATED(double baseHorner(const std::vector > &coefficients, double x, double y)); - bool do_debug(void){return get_debug_level()>=18;} + bool do_debug(void){return get_debug_level()>=500;} }; diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index 0b83405e..e353c1b9 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -27,7 +27,7 @@ IncompressibleBackend::IncompressibleBackend(IncompressibleFluid* fluid) { //this->_fractions_id = fluid->getxid(); this->fluid = fluid; if (this->fluid->is_pure()){ - this->set_fractions(std::vector(1,0)); + this->set_fractions(std::vector(1,1.0)); } } @@ -35,7 +35,7 @@ IncompressibleBackend::IncompressibleBackend(const std::string &fluid_name) { this->fluid = &get_incompressible_fluid(fluid_name); //this->_fractions_id = this->fluid->getxid(); if (this->fluid->is_pure()){ - this->set_fractions(std::vector(1,0)); + this->set_fractions(std::vector(1,1.0)); } } @@ -64,8 +64,8 @@ void IncompressibleBackend::update(CoolProp::input_pairs input_pair, double valu if (fluid->is_pure()){ this->_fluid_type = FLUID_TYPE_INCOMPRESSIBLE_LIQUID; if (get_debug_level()>=50) std::cout << format("Incompressible backend: Fluid type is %d ",this->_fluid_type) << std::endl; - if (_fractions[0]!=0.0){ - throw ValueError(format("%s is a pure fluid. The composition has to be set to a vector with one entry equal to 0.0 or nothing. %s is not valid.",this->name().c_str(),vec_to_string(_fractions).c_str())); + if (_fractions[0]!=1.0){ + throw ValueError(format("%s is a pure fluid. The composition has to be set to a vector with one entry equal to 1.0. %s is not valid.",this->name().c_str(),vec_to_string(_fractions).c_str())); } } else { this->_fluid_type = FLUID_TYPE_INCOMPRESSIBLE_SOLUTION; @@ -147,8 +147,8 @@ void IncompressibleBackend::set_fractions(const std::vector &fracti void IncompressibleBackend::set_mole_fractions(const std::vector &mole_fractions){ if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_mole_fractions with %s ",vec_to_string(mole_fractions).c_str()) << std::endl; if (mole_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mole fraction vector and not %d.",mole_fractions.size())); - if (fluid->getxid()==IFRAC_PURE) { - this->set_fractions(std::vector(1,0)); + if ((fluid->getxid()==IFRAC_PURE) && true ){//( this->_fractions[0]!=1.0 )){ + this->set_fractions(std::vector(1,1.0)); if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(mole_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl; } else if (fluid->getxid()==IFRAC_MOLE) { this->set_fractions(mole_fractions); @@ -168,8 +168,8 @@ void IncompressibleBackend::set_mole_fractions(const std::vector &m void IncompressibleBackend::set_mass_fractions(const std::vector &mass_fractions){ if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_mass_fractions with %s ",vec_to_string(mass_fractions).c_str()) << std::endl; if (mass_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mass fraction vector and not %d.",mass_fractions.size())); - if (fluid->getxid()==IFRAC_PURE) { - this->set_fractions(std::vector(1,0)); + if ((fluid->getxid()==IFRAC_PURE) && true ){// ( this->_fractions[0]!=1.0 )) { + this->set_fractions(std::vector(1,1.0)); if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(mass_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl; } else if (fluid->getxid()==IFRAC_MASS) { this->set_fractions(mass_fractions); @@ -189,8 +189,8 @@ void IncompressibleBackend::set_mass_fractions(const std::vector &m void IncompressibleBackend::set_volu_fractions(const std::vector &volu_fractions){ if (get_debug_level()>=10) std::cout << format("Incompressible backend: Called set_volu_fractions with %s ",vec_to_string(volu_fractions).c_str()) << std::endl; if (volu_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the volume fraction vector and not %d.",volu_fractions.size())); - if (fluid->getxid()==IFRAC_PURE) { - this->set_fractions(std::vector(1,0)); + if ((fluid->getxid()==IFRAC_PURE) && true ){// ( this->_fractions[0]!=1.0 )) { + this->set_fractions(std::vector(1,1.0)); if (get_debug_level()>=20) std::cout << format("Incompressible backend: Overwriting fractions for pure fluid with %s -> %s",vec_to_string(volu_fractions).c_str(),vec_to_string(this->_fractions).c_str()) << std::endl; } else if (fluid->getxid()==IFRAC_VOLUME) { this->set_fractions(volu_fractions); @@ -436,6 +436,18 @@ 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... + val = backend.calc_smass(); + backend.update( CoolProp::PT_INPUTS, p, T ); + res = backend.calc_smass(); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(val); + CAPTURE(res); + CHECK( check_abs(val,res,acc) ); + } val = fluid.u(T, p, x); res = backend.calc_umass(); @@ -535,6 +547,22 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(errmsg); CHECK( check_abs(expected,actual,acc) ); } + // entropy reference state problems + //CoolProp::set_debug_level(51); + expected = CoolProp::PropsSI("S","T",T,"P",p,"INCOMP::"+fluid+format("[%f]",x)); + actual = CoolProp::PropsSI("S","T",T,"P",p,"INCOMP::"+fluid+format("[%f]",x)); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(expected); + CAPTURE(actual); + std::string name = "INCOMP::"+fluid+format("[%f]",x); + CAPTURE(name); + std::string errmsg = CoolProp::get_global_param_string("errstring"); + CAPTURE(errmsg); + CHECK( check_abs(expected,actual,acc) ); + } } SECTION("SecCool example") { diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp index 346969ca..5e7a1ea8 100644 --- a/src/Backends/Incompressible/IncompressibleFluid.cpp +++ b/src/Backends/Incompressible/IncompressibleFluid.cpp @@ -19,15 +19,23 @@ and transport properties. //IncompressibleFluid::IncompressibleFluid(); void IncompressibleFluid::set_reference_state(double T0, double p0, double x0, double h0, double s0){ - this->Tref = T0; - this->rhoref = rho(T0,p0,x0); - this->pref = p0; - this->href = h0; - // Now we take care of the energy related values - this->uref = 0.0; - this->uref = u(T0,p0,x0) - h0; // (value without ref) - (desired ref) - this->sref = 0.0; - this->sref = s(T0,p0,x0); + 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; } void IncompressibleFluid::validate(){ @@ -130,7 +138,7 @@ double IncompressibleFluid::s (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, -1, 0, Tbase, xbase) - sref; + return poly.integral(specific_heat.coeffs, T, x, 0, -1, 0, Tbase, xbase) - _sref; 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)); @@ -147,7 +155,7 @@ 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; + 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)); @@ -422,7 +430,7 @@ double IncompressibleFluid::T_c (double Cmass, double p, double x){ } /// 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; + 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); @@ -438,7 +446,7 @@ double IncompressibleFluid::T_s (double Smass, double p, double x){ } /// 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; + 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); @@ -667,8 +675,8 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(Tref); CAPTURE(pref); CAPTURE(xref); - CAPTURE(href) - CAPTURE(res) + CAPTURE(href); + CAPTURE(res); CHECK( check_abs(href,res,acc) ); } { @@ -676,10 +684,26 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(Tref); CAPTURE(pref); CAPTURE(xref); - CAPTURE(sref) - CAPTURE(res) + CAPTURE(sref); + CAPTURE(res); CHECK( check_abs(sref,res,acc) ); - } + double res2 = XLT.s(Tref,pref,xref); + CAPTURE(Tref); + CAPTURE(pref); + CAPTURE(xref); + CAPTURE(res); + CAPTURE(res2); + CHECK( check_abs(res2,res,acc) ); + res = XLT.s(Tref,pref,xref); + CAPTURE(Tref); + CAPTURE(pref); + CAPTURE(xref); + CAPTURE(res); + CAPTURE(res2); + CHECK( check_abs(res2,res,acc) ); + } + + Tref = 25+273.15; pref = 0.0; @@ -849,6 +873,15 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(actual); CHECK( check_abs(expected,actual,acc) ); } + expected = CH3OH.s(T,p,x); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(expected); + CAPTURE(actual); + CHECK( check_abs(expected,actual,acc) ); + } expected = 0.0; actual = CH3OH.s(Tref,pref,xref); @@ -860,6 +893,15 @@ TEST_CASE("Internal consistency checks and example use cases for the incompressi CAPTURE(actual); CHECK( expected==actual ); } + expected = CH3OH.s(Tref,pref,xref); + { + CAPTURE(T); + CAPTURE(p); + CAPTURE(x); + CAPTURE(expected); + CAPTURE(actual); + CHECK( expected==actual ); + } // Compare u expected = -60043.78429641827;