Adding methane properties

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

docs/classes/methane.rst

@@ -0,0 +1,5 @@
+Methane
+========
+
+.. doxygenclass:: fprops::Methane
+   :members:

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

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
+  ]

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}
+}

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;

include/fprops/methane.h

@@ -0,0 +1,41 @@
+// SPDX-FileCopyrightText: 2024 David Andrs <andrsd@gmail.com>
+// 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<double> lead;
+    IdealGasLogTau<double> log_tau;
+    IdealGasPlanckEinsteinFunctionT<double> pefnt;
+    IdealEnthalpyEntropyOffset<double> offset;
+    ResidualPower<double> power_r;
+    ResidualPowerExp<double, unsigned int> power_exp_r;
+    ResidualGaussian<double> gauss;
+
+    PowersOfTemperature<double> eta_0;
+    FrictionTheory<double> eta_f;
+};
+
+} // namespace fprops

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<double, 3> C;
 };
 
+template <typename TYPE>
+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<double> & a, const std::vector<double> & 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<double> a;
+    std::vector<double> t;
+};
+
+///
+template <typename TYPE>
+class FrictionTheory {
+public:
+    FrictionTheory(double c1,
+                   double c2,
+                   const std::vector<double> & Ai,
+                   const std::vector<double> & Aa,
+                   double Na,
+                   const std::vector<double> & Ar,
+                   double Nr,
+                   const std::vector<double> & Aaa,
+                   double Naa,
+                   const std::vector<double> & Arr,
+                   const std::vector<double> & Adrdr,
+                   double Nrr,
+                   const std::vector<double> & Aii,
+                   double Nii,
+                   const std::vector<double> & Arrr,
+                   double Nrrr,
+                   const std::vector<double> & 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<double> Ai;
+
+    std::vector<double> Aa;
+    double Na;
+
+    std::vector<double> Ar;
+    double Nr;
+
+    std::vector<double> Aaa;
+    double Naa;
+
+    std::vector<double> Arr;
+    std::vector<double> Adrdr;
+    double Nrr;
+
+    std::vector<double> Aii;
+    double Nii;
+
+    std::vector<double> Arrr;
+    double Nrrr;
+
+    std::vector<double> Aaaa;
+    double Naaa;
+};
+
 } // namespace fprops

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_<Methane>(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);
 }