diff --git a/dumux/material/components/h2o.hh b/dumux/material/components/h2o.hh index 0d2ccdc20a13c8d656fc3f0b9875548d044954c3..0fd5d831b3b5900a90f1a518202b0760d7bc568b 100644 --- a/dumux/material/components/h2o.hh +++ b/dumux/material/components/h2o.hh @@ -771,7 +771,7 @@ public: Scalar rho = gasDensity(temperature, pressure); return Common::viscosity(temperature, rho); - }; + } /*! * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of pure water. @@ -795,7 +795,7 @@ public: Scalar rho = liquidDensity(temperature, pressure); return Common::viscosity(temperature, rho); - }; + } /*! * \brief Thermal conductivity \f$\mathrm{[[W/(m K)]}\f$ of water (IAPWS) . @@ -814,10 +814,13 @@ public: { // Thermal conductivity of water is empirically fit. // Evaluating that fitting-function outside the area of validity does not make sense. - assert( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) ) + if( not ( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) ) or (pressure <= 200e6 and ((398.15<temperature) and (temperature<=523.15)) ) or (pressure <= 150e6 and ((523.15<temperature) and (temperature<=673.15)) ) - or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) ); + or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) ) ){ + DUNE_THROW(NumericalProblem, + "Evaluating the IAPWS fit function for thermal conductivity outside range of applicability: p= " << pressure << "T= " << temperature); + } Scalar rho = liquidDensity(temperature, pressure); return Common::thermalConductivityIAPWS(temperature, rho); @@ -840,10 +843,13 @@ public: { // Thermal conductivity of water is empirically fit. // Evaluating that fitting-function outside the area of validity does not make sense. - assert( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) ) + if( not ( (pressure <= 400e6 and ((273.15<=temperature) and (temperature<=398.15)) ) or (pressure <= 200e6 and ((398.15<temperature) and (temperature<=523.15)) ) or (pressure <= 150e6 and ((523.15<temperature) and (temperature<=673.15)) ) - or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) ); + or (pressure <= 100e6 and ((673.15<temperature) and (temperature<=1073.15)) ) ) ){ + DUNE_THROW(NumericalProblem, + "Evaluating the IAPWS fit function for thermal conductivity outside range of applicability: p= " << pressure << " T= " << temperature); + } Scalar rho = gasDensity(temperature, pressure); return Common::thermalConductivityIAPWS(temperature, rho); @@ -857,7 +863,7 @@ private: Region1::tau(temperature) * Region1::dgamma_dtau(temperature, pressure) * Rs*temperature; - }; + } // the unregularized specific isobaric heat capacity static Scalar heatCap_p_Region1_(Scalar temperature, Scalar pressure) @@ -866,7 +872,7 @@ private: - pow(Region1::tau(temperature), 2 ) * Region1::ddgamma_ddtau(temperature, pressure) * Rs; - }; + } // the unregularized specific isochoric heat capacity static Scalar heatCap_v_Region1_(Scalar temperature, Scalar pressure) @@ -879,7 +885,7 @@ private: - pow(tau, 2 ) * Region1::ddgamma_ddtau(temperature, pressure) * Rs + diff; - }; + } // the unregularized specific internal energy for liquid water static Scalar internalEnergyRegion1_(Scalar temperature, Scalar pressure) @@ -888,7 +894,7 @@ private: Rs * temperature * ( Region1::tau(temperature)*Region1::dgamma_dtau(temperature, pressure) - Region1::pi(pressure)*Region1::dgamma_dpi(temperature, pressure)); - }; + } // the unregularized specific volume for liquid water static Scalar volumeRegion1_(Scalar temperature, Scalar pressure) @@ -897,7 +903,7 @@ private: Region1::pi(pressure)* Region1::dgamma_dpi(temperature, pressure) * Rs * temperature / pressure; - }; + } // the unregularized specific enthalpy for steam static Scalar enthalpyRegion2_(Scalar temperature, Scalar pressure) @@ -906,7 +912,7 @@ private: Region2::tau(temperature) * Region2::dgamma_dtau(temperature, pressure) * Rs*temperature; - }; + } // the unregularized specific internal energy for steam static Scalar internalEnergyRegion2_(Scalar temperature, Scalar pressure) @@ -915,7 +921,7 @@ private: Rs * temperature * ( Region2::tau(temperature)*Region2::dgamma_dtau(temperature, pressure) - Region2::pi(pressure)*Region2::dgamma_dpi(temperature, pressure)); - }; + } // the unregularized specific isobaric heat capacity static Scalar heatCap_p_Region2_(Scalar temperature, Scalar pressure) @@ -924,7 +930,7 @@ private: - pow(Region2::tau(temperature), 2 ) * Region2::ddgamma_ddtau(temperature, pressure) * Rs; - }; + } // the unregularized specific isochoric heat capacity static Scalar heatCap_v_Region2_(Scalar temperature, Scalar pressure) @@ -937,7 +943,7 @@ private: - pow(tau, 2 ) * Region2::ddgamma_ddtau(temperature, pressure) * Rs - diff; - }; + } // the unregularized specific volume for steam static Scalar volumeRegion2_(Scalar temperature, Scalar pressure) @@ -946,7 +952,7 @@ private: Region2::pi(pressure)* Region2::dgamma_dpi(temperature, pressure) * Rs * temperature / pressure; - }; + } }; // end class template <class Scalar> diff --git a/dumux/material/components/mesitylene.hh b/dumux/material/components/mesitylene.hh index 35afc8678dfdc7483e76730c263f411d28989922..a0aa8ffd1b22d18a55d9df75f623999618df2560 100644 --- a/dumux/material/components/mesitylene.hh +++ b/dumux/material/components/mesitylene.hh @@ -108,20 +108,29 @@ public: * \param temperature temperature of component in \f$\mathrm{[K]}\f$ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ - static Scalar liquidEnthalpy(Scalar temperature, Scalar pressure) + static Scalar liquidEnthalpy(const Scalar temperature, + const Scalar pressure) { - // Gauss quadrature rule: - // Interval: [0K; temperature (K)] - // Gauss-Legendre-Integration with variable transformation: - // \int_a^b f(T) dT \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 ) - // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1 - // here: a=0, b=actual temperature in Kelvin - // \leadsto h(T) = \int_0^T c_p(T) dT - // \approx 0.5 T * (cp( (0.5-0.5*\sqrt(1/3)) T) + cp((0.5+0.5*\sqrt(1/3)) T)) - // = 0.5 T * (cp(0.2113 T) + cp(0.7887 T) ) - - // enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input - return 0.5*temperature*(spHeatCapLiquidPhase_(0.2113*temperature) + spHeatCapLiquidPhase_(0.7887*temperature)); + // Gauss quadrature rule: + // Interval: [0K; temperature (K)] + // Gauss-Legendre-Integration with variable transformation: + // \int_a^b f(T) dT \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 ) + // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1 + // here: a=273.15K, b=actual temperature in Kelvin + // \leadsto h(T) = \int_273.15^T c_p(T) dT + // \approx 0.5 (T-273.15) * (cp( 0.5(temperature-273.15)sqrt(1/3) ) + cp(0.5(temperature-273.15)(-1)sqrt(1/3)) + + // Enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input and is + // fit over a certain temperature range. This suggests choosing an interval of integration being in the actual fit range. + // I.e. choosing T=273.15K as reference point for liquid enthalpy. + + constexpr Scalar sqrt1over3 = std::sqrt(1./3.); + const Scalar TEval1 = 0.5*(temperature-273.15)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration + const Scalar TEval2 = 0.5*(temperature-273.15)* (-1)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration + + const Scalar h_n = 0.5 * (temperature-273.15) * ( liquidHeatCapacity(TEval1, pressure) + liquidHeatCapacity(TEval2, pressure) ) ; + + return h_n; } /*! diff --git a/dumux/material/components/xylene.hh b/dumux/material/components/xylene.hh index cd4e7faced4511bf0a205d1508908d82f0bce272..eb97e4e18e9b03d5770cae0b70619acabc450f44 100644 --- a/dumux/material/components/xylene.hh +++ b/dumux/material/components/xylene.hh @@ -151,22 +151,29 @@ public: * \param temp temperature of component in \f$\mathrm{[K]}\f$ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$ */ - static Scalar liquidEnthalpy(Scalar temp, - Scalar pressure) + static Scalar liquidEnthalpy(const Scalar temperature, + const Scalar pressure) { - // Gauss quadrature rule: - // Interval: [0K; temperature (K)] - // Gauss-Legendre-Integration with variable transformation: - // \int_a^b f(T) dT \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 ) - // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1 - // here: a=0, b=actual temperature in Kelvin - // \leadsto h(T) = \int_0^T c_p(T) dT - // \approx 0.5 T * (cp( (0.5-0.5*\sqrt(1/3)) T) + cp((0.5+0.5*\sqrt(1/3)) T)) - // = 0.5 T * (cp(0.2113 T) + cp(0.7887 T) ) - - // enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input - return 0.5*temp*(spHeatCapLiquidPhase(0.2113*temp,pressure) - + spHeatCapLiquidPhase(0.7887*temp,pressure)); + // Gauss quadrature rule: + // Interval: [0K; temperature (K)] + // Gauss-Legendre-Integration with variable transformation: + // \int_a^b f(T) dT \approx (b-a)/2 \sum_i=1^n \alpha_i f( (b-a)/2 x_i + (a+b)/2 ) + // with: n=2, legendre -> x_i = +/- \sqrt(1/3), \apha_i=1 + // here: a=273.15K, b=actual temperature in Kelvin + // \leadsto h(T) = \int_273.15^T c_p(T) dT + // \approx 0.5 (T-273.15) * (cp( 0.5(temperature-273.15)sqrt(1/3) ) + cp(0.5(temperature-273.15)(-1)sqrt(1/3)) + + // Enthalpy may have arbitrary reference state, but the empirical/fitted heatCapacity function needs Kelvin as input and is + // fit over a certain temperature range. This suggests choosing an interval of integration being in the actual fit range. + // I.e. choosing T=273.15K as reference point for liquid enthalpy. + + constexpr Scalar sqrt1over3 = std::sqrt(1./3.); + const Scalar TEval1 = 0.5*(temperature-273.15)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration + const Scalar TEval2 = 0.5*(temperature-273.15)* (-1)* sqrt1over3 + 0.5*(273.15+temperature) ; // evaluation points according to Gauss-Legendre integration + + const Scalar h_n = 0.5 * (temperature-273.15) * ( liquidHeatCapacity(TEval1, pressure) + liquidHeatCapacity(TEval2, pressure) ) ; + + return h_n; } /*! diff --git a/dumux/material/fluidsystems/h2oairfluidsystem.hh b/dumux/material/fluidsystems/h2oairfluidsystem.hh index 3c16c7481e1567b513ed0aae13448a1c57bc3913..9f890c7c853d8d1e6be33f0064b57f9e2c6ff4cd 100644 --- a/dumux/material/fluidsystems/h2oairfluidsystem.hh +++ b/dumux/material/fluidsystems/h2oairfluidsystem.hh @@ -305,7 +305,7 @@ public: init(/*tempMin=*/273.15, /*tempMax=*/623.15, /*numTemp=*/100, - /*pMin=*/-10, + /*pMin=*/-10., /*pMax=*/20e6, /*numP=*/200); } @@ -683,7 +683,7 @@ public: } }; -} // end namepace FluidSystems +} // end namespace FluidSystems #ifdef DUMUX_PROPERTIES_HH // forward defintions of the property tags