From 2b1e0a1cc784507c9b2bfec50a76e786002e2b17 Mon Sep 17 00:00:00 2001 From: David Andrs Date: Thu, 8 Aug 2024 20:32:47 -0600 Subject: [PATCH] Adding methane properties --- docs/CMakeLists.txt | 1 + docs/classes/methane.rst | 5 + docs/fluids/methane.rst | 31 ++++ docs/pyplots/ch4.yml | 16 ++ docs/references.bib | 46 +++++ include/fprops/helmholtz.h | 1 - include/fprops/methane.h | 41 +++++ include/fprops/transport_models.h | 153 ++++++++++++++++ python/src/fprops.cpp | 8 + src/methane.cpp | 292 ++++++++++++++++++++++++++++++ test/src/methane_test.cpp | 105 +++++++++++ 11 files changed, 698 insertions(+), 1 deletion(-) create mode 100644 docs/classes/methane.rst create mode 100644 docs/fluids/methane.rst create mode 100644 docs/pyplots/ch4.yml create mode 100644 include/fprops/methane.h create mode 100644 src/methane.cpp create mode 100644 test/src/methane_test.cpp diff --git a/docs/CMakeLists.txt b/docs/CMakeLists.txt index c7c7329..c491fb6 100644 --- a/docs/CMakeLists.txt +++ b/docs/CMakeLists.txt @@ -49,6 +49,7 @@ if(DOXYGEN_FOUND) add_fluid(n2 n2.yml) add_fluid(co2 co2.yml) add_fluid(o2 o2.yml) + add_fluid(ch4 ch4.yml) add_custom_command( OUTPUT diff --git a/docs/classes/methane.rst b/docs/classes/methane.rst new file mode 100644 index 0000000..a3ce625 --- /dev/null +++ b/docs/classes/methane.rst @@ -0,0 +1,5 @@ +Methane +======== + +.. doxygenclass:: fprops::Methane + :members: diff --git a/docs/fluids/methane.rst b/docs/fluids/methane.rst new file mode 100644 index 0000000..930251f --- /dev/null +++ b/docs/fluids/methane.rst @@ -0,0 +1,31 @@ +Methane +======= + +Range of validity: + +- Temperature: `54 - 625 K` +- Pressure: to `1000 MPa` + +Computation of thermodynamic properties is based on :cite:t:`10.1063/1.555898`. +Viscosity and thermal conductivity is based on :cite:t:`10.1063/1.555828` and cite:t:`doi:10.1021/jp0618577`, respectively. + +Verification against CoolProp package +------------------------------------- + +.. image:: ../pyplots/ch4_err_rho.png + +.. image:: ../pyplots/ch4_err_u.png + +.. image:: ../pyplots/ch4_err_h.png + +.. image:: ../pyplots/ch4_err_s.png + +.. image:: ../pyplots/ch4_err_w.png + +.. image:: ../pyplots/ch4_err_c_p.png + +.. image:: ../pyplots/ch4_err_c_v.png + +.. image:: ../pyplots/ch4_err_mu.png + +.. image:: ../pyplots/ch4_err_k.png diff --git a/docs/pyplots/ch4.yml b/docs/pyplots/ch4.yml new file mode 100644 index 0000000..d4a2c44 --- /dev/null +++ b/docs/pyplots/ch4.yml @@ -0,0 +1,16 @@ +fprops: + name: ch4 + coolprop: methane + fprops: Methane + p: [ + 0.1, 0.2, + 0.5, 1, 2, 5, 10, 20, 50, 100, + 200, 500, 1000, 2000 + ] + T: [ + 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130, 135, 140, 145, 150, + 160, 170, 180, 190, + 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, + 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, + 400, 450, 500, 550, 600, 6250 + ] diff --git a/docs/references.bib b/docs/references.bib index e529bd9..d880e9c 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -127,3 +127,49 @@ @article{10.1063/1.555897 url = {https://doi.org/10.1063/1.555897}, eprint = {https://pubs.aip.org/aip/jpr/article-pdf/20/5/917/9766717/917\_1\_online.pdf}, } + +@article{10.1063/1.555898, + author = {Setzmann, U. and Wagner, W.}, + title = "{A New Equation of State and Tables of Thermodynamic Properties for Methane Covering the Range from the Melting Line to 625 K at Pressures up to 1000 MPa}", + journal = {Journal of Physical and Chemical Reference Data}, + volume = {20}, + number = {6}, + pages = {1061-1155}, + year = {1991}, + month = {11}, + abstract = "{This work reviews the data on thermodynamic properties of methane which were available up to the middle of 1991 and presents a new equation of state in the form of a fundamental equation explicit in the Helmholtz free energy. A new strategy for optimizing the structure of empirical thermodynamic correlation equations was used to determine the functional form of the equation. The Helmholtz function containing 40 fitted coefficients was fitted to selected experimental data of the following properties: (a) thermal properties of the single phase (pρT) and (b) of the liquid‐vapor saturation curve (psρ′ρ″) including the Maxwell criterion, (c) speed of sound w, (d) isochoric heat capacity cv, (e) isobaric heat capacity cp, (f) difference of enthalpy Δh, and (g) second virial coefficient B. Independent equations are also included for the vapor pressure, the saturated liquid and vapor densities, the isobaric ideal gas heat capacity and the melting pressure as functions of temperature. Tables for the thermodynamic properties of methane from 90 K to 620 K for pressures up to 1000 MPa are presented. For the density, uncertainties of ±0.03\\% for pressures below 12 MPa and temperatures below 350 K and ±0.03\\% to ±0.15\\% for higher pressures and temperatures are estimated. For the speed of sound, the uncertainty ranges from ±0.03\\% to ±0.3\\% depending on temperature and pressure. Heat capacities may be generally calculated within an uncertainty of ±1\\%. To verify the accuracy of the new formulation, the calculated property values are compared with selected experimental results and existing equations of state for methane. The new equation of state corresponds to the new International Temperature Scale of 1990 (ITS‐90) and is extrapolable up to pressures of 20000 MPa.}", + issn = {0047-2689}, + doi = {10.1063/1.555898}, + url = {https://doi.org/10.1063/1.555898}, + eprint = {https://pubs.aip.org/aip/jpr/article-pdf/20/6/1061/8183450/1061\_1\_online.pdf}, +} + +@article{10.1063/1.555828, + author = {Friend, Daniel G. and Ely, James F. and Ingham, Hepburn}, + title = "{Thermophysical Properties of Methane}", + journal = {Journal of Physical and Chemical Reference Data}, + volume = {18}, + number = {2}, + pages = {583-638}, + year = {1989}, + month = {04}, + abstract = "{New correlations for the thermophysical properties of fluid methane are presented. The correlations are based on a critical evaluation of the available experimental data and have been developed to represent these data over a broad range of the state variables. Estimates for the accuracy of the equations and comparisons with measured properties are given. The reasons for this new study of methane include significant new and more accurate data, and improvements in the correlation functions which allow increased accuracy of the correlations especially in the extended critical region. For the thermodynamic properties, a classical equation for the molar Helmholtz energy, which contains terms multiplied by the exponential of the quadratic and quartic powers of the system density, is used. The resulting equation of state is accurate from about 91 to 600 K for pressures \\<100 MPa and was developed by considering PVT, second virial coefficient, heat capacity, and sound speed data. Tables of coefficients and equations are presented to allow the calculation of these and other thermodynamic quantities. Ancillary equations for properties along the liquid–vapor phase boundary, which are consistent with the equation of state and lowest order scaling theory, are also given. For the viscosity of fluid methane, a low‐density contribution based on theory is combined with an empirical representation of the excess contribution. The approximate range of the resulting correlation is 91 to 400 K for pressures \\<55 MPa. The correlation for the thermal conductivity includes a theoretically based expression for the critical enhancement; the range for the resulting correlation is about 91 to 700 K for pressures below 100 MPa.}", + issn = {0047-2689}, + doi = {10.1063/1.555828}, + url = {https://doi.org/10.1063/1.555828}, + eprint = {https://pubs.aip.org/aip/jpr/article-pdf/18/2/583/11744744/583\_1\_online.pdf}, +} + +@article{doi:10.1021/jp0618577, + author = {Quiñones-Cisneros, Sergio E. and Deiters, Ulrich K.}, + title = {Generalization of the Friction Theory for Viscosity Modeling}, + journal = {The Journal of Physical Chemistry B}, + volume = {110}, + number = {25}, + pages = {12820-12834}, + year = {2006}, + doi = {10.1021/jp0618577}, + note = {PMID: 16800618}, + URL = {https://doi.org/10.1021/jp0618577}, + eprint = {https://doi.org/10.1021/jp0618577} +} diff --git a/include/fprops/helmholtz.h b/include/fprops/helmholtz.h index d16c210..e891147 100644 --- a/include/fprops/helmholtz.h +++ b/include/fprops/helmholtz.h @@ -1009,7 +1009,6 @@ class Helmholtz : public SinglePhaseFluidProperties { double delta(double rho) const; double tau(double T) const; -private: double temperature(double u, double tau, double da_dt) const; double pressure(double rho, double T, double delta, double da_dd) const; double internal_energy(double T, double tau, double da_dt) const; diff --git a/include/fprops/methane.h b/include/fprops/methane.h new file mode 100644 index 0000000..e016ab1 --- /dev/null +++ b/include/fprops/methane.h @@ -0,0 +1,41 @@ +// SPDX-FileCopyrightText: 2024 David Andrs +// SPDX-License-Identifier: MIT + +#pragma once + +#include "fprops/helmholtz.h" +#include "fprops/transport_models.h" + +namespace fprops { + +/// Methane fluid properties +/// +/// Reference: +/// +class Methane : public Helmholtz { +public: + Methane(); + +private: + [[nodiscard]] double alpha(double delta, double tau) const override; + [[nodiscard]] double dalpha_ddelta(double delta, double tau) const override; + [[nodiscard]] double dalpha_dtau(double delta, double tau) const override; + [[nodiscard]] double d2alpha_ddelta2(double delta, double tau) const override; + [[nodiscard]] double d2alpha_dtau2(double delta, double tau) const override; + [[nodiscard]] double d2alpha_ddeltatau(double delta, double tau) const override; + [[nodiscard]] double mu_from_rho_T(double rho, double T) const override; + [[nodiscard]] double k_from_rho_T(double rho, double T) const override; + + IdealGasLead lead; + IdealGasLogTau log_tau; + IdealGasPlanckEinsteinFunctionT pefnt; + IdealEnthalpyEntropyOffset offset; + ResidualPower power_r; + ResidualPowerExp power_exp_r; + ResidualGaussian gauss; + + PowersOfTemperature eta_0; + FrictionTheory eta_f; +}; + +} // namespace fprops diff --git a/include/fprops/transport_models.h b/include/fprops/transport_models.h index ac5afbe..2ef6889 100644 --- a/include/fprops/transport_models.h +++ b/include/fprops/transport_models.h @@ -518,6 +518,12 @@ class QuadraticGeneralFrictionTheory { { } + /// Evaluate the model + /// + /// @param tau + /// @param p_a Pressure [bar] + /// @param p_r Pressure [bar] + /// @param p_id Pressure [bar] TYPE value(double tau, double p_a, double p_r, double p_id) const { @@ -557,4 +563,151 @@ class QuadraticGeneralFrictionTheory { std::array C; }; +template +class PowersOfTr { +public: + /// + /// + /// @param T_reducing Reducing temperature [K] + /// @param a \f$a_i\f$ + /// @param t \f$t_i\f$ + PowersOfTr(double T_reducing, const std::vector & a, const std::vector & t) : + T_reducing(T_reducing), + a(a), + t(t) + { + } + + TYPE + value(TYPE T) const + { + TYPE Tr = T / this->T_reducing; + TYPE summer = 0; + for (std::size_t i = 0; i < a.size(); ++i) + summer += this->a[i] * pow(Tr, this->t[i]); + return summer; + } + +private: + double T_reducing; + std::vector a; + std::vector t; +}; + +/// +template +class FrictionTheory { +public: + FrictionTheory(double c1, + double c2, + const std::vector & Ai, + const std::vector & Aa, + double Na, + const std::vector & Ar, + double Nr, + const std::vector & Aaa, + double Naa, + const std::vector & Arr, + const std::vector & Adrdr, + double Nrr, + const std::vector & Aii, + double Nii, + const std::vector & Arrr, + double Nrrr, + const std::vector & Aaaa, + double Naaa) : + c1(c1), + c2(c2), + Ai(Ai), + Aa(Aa), + Na(Na), + Ar(Ar), + Nr(Nr), + Aaa(Aaa), + Naa(Naa), + Arr(Arr), + Adrdr(Adrdr), + Nrr(Nrr), + Aii(Aii), + Nii(Nii), + Arrr(Arrr), + Nrrr(Nrrr), + Aaaa(Aaaa), + Naaa(Naaa) + { + } + + TYPE + value(double tau, double p, double pr, double pid) const + { + TYPE kii = 0, krrr = 0, kaaa = 0, krr, kdrdr; + + TYPE psi1 = std::exp(tau) - this->c1; + TYPE psi2 = std::exp(math::pow(tau, 2)) - this->c2; + + TYPE ki = (this->Ai[0] + this->Ai[1] * psi1 + this->Ai[2] * psi2) * tau; + + TYPE ka = + (this->Aa[0] + this->Aa[1] * psi1 + this->Aa[2] * psi2) * math::pow(tau, this->Na); + TYPE kr = + (this->Ar[0] + this->Ar[1] * psi1 + this->Ar[2] * psi2) * math::pow(tau, this->Nr); + TYPE kaa = + (this->Aaa[0] + this->Aaa[1] * psi1 + this->Aaa[2] * psi2) * math::pow(tau, this->Naa); + if (this->Arr.empty()) { + krr = 0; + kdrdr = (this->Adrdr[0] + this->Adrdr[1] * psi1 + this->Adrdr[2] * psi2) * + math::pow(tau, this->Nrr); + } + else { + krr = (this->Arr[0] + this->Arr[1] * psi1 + this->Arr[2] * psi2) * + math::pow(tau, this->Nrr); + kdrdr = 0; + } + if (!this->Aii.empty()) { + kii = (this->Aii[0] + this->Aii[1] * psi1 + this->Aii[2] * psi2) * + math::pow(tau, this->Nii); + } + if (!this->Arrr.empty() && !this->Aaaa.empty()) { + krrr = (this->Arrr[0] + this->Arrr[1] * psi1 + this->Arrr[2] * psi2) * + math::pow(tau, this->Nrrr); + kaaa = (this->Aaaa[0] + this->Aaaa[1] * psi1 + this->Aaaa[2] * psi2) * + math::pow(tau, this->Naaa); + } + + TYPE pa = p - pr; + TYPE deltapr = pr - pid; + TYPE 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; + } + +private: + double c1; + double c2; + + std::vector Ai; + + std::vector Aa; + double Na; + + std::vector Ar; + double Nr; + + std::vector Aaa; + double Naa; + + std::vector Arr; + std::vector Adrdr; + double Nrr; + + std::vector Aii; + double Nii; + + std::vector Arrr; + double Nrrr; + + std::vector Aaaa; + double Naaa; +}; + } // namespace fprops diff --git a/python/src/fprops.cpp b/python/src/fprops.cpp index c439a95..67ffaec 100644 --- a/python/src/fprops.cpp +++ b/python/src/fprops.cpp @@ -10,6 +10,7 @@ #include "fprops/helium.h" #include "fprops/carbon_dioxide.h" #include "fprops/oxygen.h" +#include "fprops/methane.h" #include "version.h" using namespace fprops; @@ -86,4 +87,11 @@ PYBIND11_MODULE(fprops, m) .def("rho_p", &Oxygen::rho_p) .def("p_T", &Oxygen::p_T) .def("v_u", &Oxygen::v_u); + + py::class_(m, "Methane") + .def(py::init()) + .def("rho_T", &Methane::rho_T) + .def("rho_p", &Methane::rho_p) + .def("p_T", &Methane::p_T) + .def("v_u", &Methane::v_u); } diff --git a/src/methane.cpp b/src/methane.cpp new file mode 100644 index 0000000..498a2df --- /dev/null +++ b/src/methane.cpp @@ -0,0 +1,292 @@ +// SPDX-FileCopyrightText: 2024 David Andrs +// SPDX-License-Identifier: MIT + +#include "fprops/methane.h" + +namespace fprops { + +static const double GAS_CONSTANT = 8.314510; +static const double MOLAR_MASS = 16.0428e-3; +static const double T_CRIT = 190.564; +static const double RHO_MOLAR_CRIT = 10139.128; +static const double RHO_CRIT = RHO_MOLAR_CRIT * MOLAR_MASS; + +Methane::Methane() : + Helmholtz(GAS_CONSTANT, MOLAR_MASS, RHO_CRIT, T_CRIT), + lead(9.91243972, -6.33270087), + log_tau(3.0016), + pefnt(T_CRIT, { 0.008449, 4.6942, 3.4865, 1.6572, 1.4115 }, { 648, 1957, 3895, 5705, 15080 }), + offset(-12.8829893867948, 9.22344625310864), + power_r({ 0.04367901028, + 0.6709236199, + -1.765577859, + 0.8582330241, + -1.206513052, + 0.512046722, + -0.0004000010791, + -0.01247842423, + 0.03100269701, + 0.001754748522, + -3.171921605e-06, + -2.24034684e-06, + 2.947056156e-07 }, + { 1, 1, 1, 2, 2, 2, 2, 3, 4, 4, 8, 9, 10 }, + { -0.5, 0.5, 1, 0.5, 1, 1.5, 4.5, 0, 1, 3, 1, 3, 3 }), + power_exp_r({ 0.1830487909, 0.1511883679, -0.4289363877, 0.06894002446, -0.01408313996, + -0.0306305483, -0.02969906708, -0.01932040831, -0.1105739959, 0.09952548995, + 0.008548437825, -0.06150555662, -0.04291792423, -0.0181320729, 0.0344590476, + -0.00238591945, -0.01159094939, 0.06641693602, -0.0237154959, -0.03961624905, + -0.01387292044, 0.03389489599, -0.002927378753 }, + { 1, 1, 1, 2, 4, 5, 6, 1, 2, 3, 4, 4, 3, 5, 5, 8, 2, 3, 4, 4, 4, 5, 6 }, + { 0, 1, 2, 0, 0, 2, 2, 5, 5, 5, 2, 4, 12, 8, 10, 10, 10, 14, 12, 18, 22, 18, 14 }, + { 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4 }), + gauss({ 9.324799946e-05, -6.287171518, 12.71069467, -6.423953466 }, + { 2, 0, 0, 0 }, + { 2, 0, 1, 2 }, + { 20, 40, 40, 40 }, + { 1, 1, 1, 1 }, + { 200, 250, 250, 250 }, + { 1.07, 1.11, 1.11, 1.11 }), + // + eta_0({ 2.60536e-06, -1.85247e-05, 2.34216e-05, 0.0 }, { 0, 0.25, 0.5, 0.75 }), + eta_f(1, + 1, + { 3.49668e-08, -1.73176e-08, 0.0 }, + { -3.12118e-08, 1.99422e-10, 0 }, + 1, + { 5.98858e-08, -4.91143e-08, 0.0 }, + 1, + { -8.52992e-13, -3.58009e-13, 0.0 }, + 3, + {}, + { 1.60099e-11, 8.50221e-13, 0.0 }, + 3, + { -3.55631e-10, 2.80326e-10, 0.0 }, + 3, + {}, + 0, + {}, + 0) +{ +} + +double +Methane::alpha(double delta, double tau) const +{ + // clang-format off + return + this->lead.alpha(delta, tau) + + this->log_tau.alpha(delta, tau) + + this->pefnt.alpha(delta, tau) + + this->offset.alpha(delta, tau) + + this->power_r.alpha(delta, tau) + + this->power_exp_r.alpha(delta, tau) + + this->gauss.alpha(delta, tau); + // clang-format on +} + +double +Methane::dalpha_ddelta(double delta, double tau) const +{ + // clang-format off + return + this->lead.ddelta(delta, tau) + + this->power_r.ddelta(delta, tau) + + this->power_exp_r.ddelta(delta, tau) + + this->gauss.ddelta(delta, tau); + // clang-format on +} + +double +Methane::dalpha_dtau(double delta, double tau) const +{ + // clang-format off + return + this->lead.dtau(delta, tau) + + this->log_tau.dtau(delta, tau) + + this->pefnt.dtau(delta, tau) + + this->offset.dtau(delta, tau) + + this->power_r.dtau(delta, tau) + + this->power_exp_r.dtau(delta, tau) + + this->gauss.dtau(delta, tau); + // clang-format on +} + +double +Methane::d2alpha_ddelta2(double delta, double tau) const +{ + // clang-format off + return + this->lead.d2delta(delta, tau) + + this->power_r.d2delta(delta, tau) + + this->power_exp_r.d2delta(delta, tau) + + this->gauss.d2delta(delta, tau); + // clang-format on +} + +double +Methane::d2alpha_dtau2(double delta, double tau) const +{ + // clang-format off + return + this->log_tau.d2tau(delta, tau) + + this->pefnt.d2tau(delta, tau) + + this->power_r.d2tau(delta, tau) + + this->power_exp_r.d2tau(delta, tau) + + this->gauss.d2tau(delta, tau); + // clang-format on +} + +double +Methane::d2alpha_ddeltatau(double delta, double tau) const +{ + // clang-format off + return + this->power_r.d2deltatau(delta, tau) + + this->power_exp_r.d2deltatau(delta, tau) + + this->gauss.d2deltatau(delta, tau); + // clang-format on +} + +double +Methane::mu_from_rho_T(double rho, double T) const +{ + auto Tr = T / T_CRIT; + + auto d = delta(rho); + auto t = tau(T); + auto da_dd = dalpha_ddelta(d, t); + auto p = pressure(rho, T, d, da_dd); + // [bar]; 1e5 for conversion from Pa -> bar + auto p_bar = p / 1e5; + auto rho_molar = rho / MOLAR_MASS; + auto dp_dT = GAS_CONSTANT * rho_molar * d * da_dd; + auto p_r_bar = T * dp_dT / 1e5; + auto p_a_bar = p_bar - p_r_bar; + auto p_id_bar = rho_molar * GAS_CONSTANT * T / 1e5; + // clang-format off + return + this->eta_0.value(Tr) + + this->eta_f.value(t, p_a_bar, p_r_bar, p_id_bar); + // clang-format on +} + +double +Methane::k_from_rho_T(double rho, double T) const +{ + const double d = delta(rho); + const double t = tau(T); + + // NOTE: slightly different rho_crit and T_crit used in computation of + // `eta`s + double delta_eta = rho / 10139.0; + double tau_eta = 190.55 / T; + + // Viscosity formulation from Friend, JPCRD, 1989 + // Dilute + double C[] = { 0, + -3.0328138281, + 16.918880086, + -37.189364917, + 41.288861858, + -24.615921140, + 8.9488430959, + -1.8739245042, + 0.20966101390, + -9.6570437074e-3 }; + double OMEGA22_sum = 0; + double t_eta = T / 174.0; + for (int i = 1; i <= 9; ++i) + OMEGA22_sum += C[i] * math::pow(t_eta, (i - 1.0) / 3.0 - 1.0); + double eta_dilute = 10.50 * sqrt(t_eta) * OMEGA22_sum; + double re[] = { 0, 1, 1, 2, 2, 2, 3, 3, 4, 4, 1, 1 }; + double se[] = { 0, 0, 1, 0, 1, 1.5, 0, 2, 0, 1, 0, 1 }; + double ge[] = { 0, 0.41250137, -0.14390912, 0.10366993, 0.40287464, -0.24903524, + -0.12953131, 0.06575776, 0.02566628, -0.03716526, -0.38798341, 0.03533815 }; + double sum1 = 0; + for (int i = 1; i <= 9; ++i) + sum1 += ge[i] * math::pow(delta_eta, re[i]) * math::pow(tau_eta, se[i]); + double sum2 = 0; + for (int i = 10; i <= 11; ++i) + sum2 += ge[i] * math::pow(delta_eta, re[i]) * math::pow(tau_eta, se[i]); + + double eta_residual = 12.149 * sum1 / (1 + sum2); + double eta = eta_residual + eta_dilute; + + // Dilute + auto d2alpha0_dtau2 = this->log_tau.d2tau(d, t) + this->pefnt.d2tau(d, t); + + double f_int = 1.458850 - 0.4377162 / t_eta; + auto lambda_dilute = 0.51828 * eta_dilute * + (3.75 - f_int * (math::pow<2>(tau(T)) * d2alpha0_dtau2 + 1.5)); // [mW/m/K] + // Residual + double rl[] = { 0, 1, 3, 4, 4, 5, 5, 2 }; + double sl[] = { 0, 0, 0, 0, 1, 0, 1, 0 }; + double jl[] = { 0, 2.4149207, 0.55166331, -0.52837734, + 0.073809553, 0.24465507, -0.047613626, 1.5554612 }; + double sum = 0; + for (int i = 1; i <= 6; ++i) { + sum += jl[i] * math::pow(delta_eta, rl[i]) * math::pow(tau_eta, sl[i]); + } + double delta_sigma_star = 1.0; // Looks like a typo in Friend - should be 1 instead of 11 + // FIXME: uncomment + // if (HEOS.T() < HEOS.T_critical() && HEOS.rhomolar() < HEOS.rhomolar_critical()) { + // delta_sigma_star = HEOS.saturation_ancillary(iDmolar, 1, iT, HEOS.T()) / + // HEOS.keyed_output(CoolProp::irhomolar_critical); + // } + auto lambda_residual = + 6.29638 * (sum + jl[7] * math::pow<2>(delta_eta) / delta_sigma_star); // [mW/m/K] + + // FIXME: Ignoring critical region + // // Critical region + // double Tstar = 1 - 1 / tau; + // double rhostar = 1 - delta; + // double F_T = 2.646, F_rho = 2.678, F_A = -0.637; + // double F = std::exp(-F_T * std::sqrt(std::abs(Tstar)) - F_rho * math::pow<2>(rhostar) - F_A * + // rhostar); + // double CHI_T_star; + // if (std::abs(Tstar) < 0.03) { + // if (std::abs(rhostar) < 1e-16) { + // // Equation 26 + // const double LAMBDA = 0.0801, gamma = 1.190; + // CHI_T_star = LAMBDA * math::pow(std::abs(Tstar), -gamma); + // } + // else if (std::abs(rhostar) < 0.03) { + // // Equation 23 + // const double beta = 0.355, W = -1.401, S = -6.098, E = 0.287, a = 3.352, b = 0.732, + // R = 0.535, Q = 0.1133; + // double OMEGA = W * Tstar * math::pow(std::abs(rhostar), -1 / beta); + // double theta = 1; + // if (Tstar < -math::pow(std::abs(rhostar), -1 / beta) / S) { + // theta = 1 + E * math::pow(1 + S * Tstar * math::pow(std::abs(rhostar), -1 / + // beta), 2 * beta); + // } + // CHI_T_star = + // Q * math::pow(std::abs(rhostar), -a) * math::pow(theta, b) / (theta + OMEGA * + // (theta + R)); + // } + // else { + // // Equation 19a + // CHI_T_star = 0.28631 * delta * tau / + // (1 + 2 * delta * HEOS.dalphar_dDelta() + + // math::pow<2>(delta) * HEOS.d2alphar_dDelta2()); + // } + // } + // else { + // // Equation 19a + // CHI_T_star = + // 0.28631 * delta * tau / + // (1 + 2 * delta * HEOS.dalphar_dDelta() + math::pow<2>(delta) * + // HEOS.d2alphar_dDelta2()); + // } + + // auto lambda_critical = 91.855 / (eta * math::pow<2>(tau)) * + // math::pow<2>(1 + delta * HEOS.dalphar_dDelta() - + // delta * tau * HEOS.d2alphar_dDelta_dTau()) * + // math::pow(CHI_T_star, 0.4681) * F; //[mW/m/K] + // auto lambda = (lambda_dilute + lambda_residual + lambda_critical) * 0.001; + + return (lambda_dilute + lambda_residual) * 0.001; +} + +} // namespace fprops diff --git a/test/src/methane_test.cpp b/test/src/methane_test.cpp new file mode 100644 index 0000000..bb3b147 --- /dev/null +++ b/test/src/methane_test.cpp @@ -0,0 +1,105 @@ +#include "gtest/gtest.h" +#include +#include "exception_test_macros.h" +#include "fprops/methane.h" + +using namespace fprops; +using namespace testing; + +namespace { + +// T = 280 K, p = 1 MPa + +double w = 431.97873645972703; +double rho = 7.043156624753545; +double h = 859902.14463317941; +double s = 5323.7095834448064; +double T = 280; +double p = 1e6; +double cp = 2257.7903212900533; +double cv = 1680.8602830296586; +double u = 717920.35247972247; +double v = 1. / rho; +double mu = 1.0738075482060213e-05; +double k = 0.03232260677413022; + +State gold1 = { u, v, rho, p, T, mu, cp, cv, s, k, h, w }; + +} // namespace + +TEST(MethaneTest, rho_T) +{ + Methane fp; + auto state = fp.rho_T(rho, T); + + EXPECT_DOUBLE_EQ(state.rho, gold1.rho); + EXPECT_DOUBLE_EQ(state.T, gold1.T); + EXPECT_DOUBLE_EQ(state.p, gold1.p); + EXPECT_DOUBLE_EQ(state.u, gold1.u); + EXPECT_DOUBLE_EQ(state.cv, gold1.cv); + EXPECT_DOUBLE_EQ(state.cp, gold1.cp); + EXPECT_NEAR(state.mu, gold1.mu, 2e-7); + EXPECT_NEAR(state.k, gold1.k, 1e-3); + EXPECT_DOUBLE_EQ(state.v, gold1.v); + EXPECT_DOUBLE_EQ(state.s, gold1.s); + EXPECT_DOUBLE_EQ(state.h, gold1.h); + EXPECT_DOUBLE_EQ(state.w, gold1.w); +} + +TEST(MethaneTest, rho_p) +{ + Methane fp; + auto state = fp.rho_p(rho, p); + + EXPECT_DOUBLE_EQ(state.rho, gold1.rho); + EXPECT_DOUBLE_EQ(state.T, gold1.T); + EXPECT_DOUBLE_EQ(state.p, gold1.p); + EXPECT_DOUBLE_EQ(state.u, gold1.u); + EXPECT_DOUBLE_EQ(state.cv, gold1.cv); + EXPECT_DOUBLE_EQ(state.cp, gold1.cp); + EXPECT_NEAR(state.mu, gold1.mu, 2e-7); + EXPECT_NEAR(state.k, gold1.k, 1e-3); + EXPECT_DOUBLE_EQ(state.v, gold1.v); + EXPECT_DOUBLE_EQ(state.s, gold1.s); + EXPECT_DOUBLE_EQ(state.h, gold1.h); + EXPECT_DOUBLE_EQ(state.w, gold1.w); +} + +TEST(MethaneTest, p_T) +{ + Methane fp; + auto state = fp.p_T(p, T); + + EXPECT_DOUBLE_EQ(state.rho, gold1.rho); + EXPECT_DOUBLE_EQ(state.T, gold1.T); + EXPECT_DOUBLE_EQ(state.p, gold1.p); + EXPECT_DOUBLE_EQ(state.u, gold1.u); + EXPECT_DOUBLE_EQ(state.cv, gold1.cv); + EXPECT_DOUBLE_EQ(state.cp, gold1.cp); + EXPECT_NEAR(state.mu, gold1.mu, 2e-7); + EXPECT_NEAR(state.k, gold1.k, 1e-3); + EXPECT_DOUBLE_EQ(state.v, gold1.v); + EXPECT_DOUBLE_EQ(state.s, gold1.s); + EXPECT_DOUBLE_EQ(state.h, gold1.h); + EXPECT_DOUBLE_EQ(state.w, gold1.w); +} + +TEST(MethaneTest, v_u) +{ + Methane fp; + auto state0 = fp.p_T(p, T); + State state = fp.v_u(state0.v, state0.u); + + EXPECT_DOUBLE_EQ(state.rho, gold1.rho); + EXPECT_DOUBLE_EQ(state.T, gold1.T); + EXPECT_DOUBLE_EQ(state.p, gold1.p); + EXPECT_DOUBLE_EQ(state.u, gold1.u); + EXPECT_DOUBLE_EQ(state.cv, gold1.cv); + EXPECT_DOUBLE_EQ(state.cp, gold1.cp); + EXPECT_NEAR(state.mu, gold1.mu, 2e-7); + EXPECT_NEAR(state.k, gold1.k, 1e-3); + EXPECT_DOUBLE_EQ(state.v, gold1.v); + EXPECT_DOUBLE_EQ(state.s, gold1.s); + EXPECT_DOUBLE_EQ(state.h, gold1.h); + EXPECT_DOUBLE_EQ(state.w, gold1.w); +}