Merge pull request #86 from andrsd/ammonia

Adding ammonia

docs/CMakeLists.txt

@@ -50,6 +50,7 @@ if(DOXYGEN_FOUND)
     add_fluid(co2 co2.yml)
     add_fluid(o2 o2.yml)
     add_fluid(ch4 ch4.yml)
+    add_fluid(nh3 nh3.yml)
 
     add_custom_command(
         OUTPUT

docs/classes/ammonia.rst

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

docs/fluids/ammonia.rst

@@ -0,0 +1,31 @@
+Ammonia
+=======
+
+Range of validity:
+
+- Temperature: `melting line - 725 K`
+- Pressure: to `1000 MPa`
+
+Computation of thermodynamic properties is based on :cite:t:`gao2020thermodynamic`.
+Viscosity and thermal conductivity is based on :cite:t:`fenghour1995viscosity` and :cite:t:`tufeu1984thermal`.
+
+Verification against CoolProp package
+-------------------------------------
+
+.. image:: ../pyplots/nh3_err_rho.png
+
+.. image:: ../pyplots/nh3_err_u.png
+
+.. image:: ../pyplots/nh3_err_h.png
+
+.. image:: ../pyplots/nh3_err_s.png
+
+.. image:: ../pyplots/nh3_err_w.png
+
+.. image:: ../pyplots/nh3_err_c_p.png
+
+.. image:: ../pyplots/nh3_err_c_v.png
+
+.. image:: ../pyplots/nh3_err_mu.png
+
+.. image:: ../pyplots/nh3_err_k.png

docs/pyplots/nh3.yml

@@ -0,0 +1,16 @@
+fprops:
+  name: nh3
+  coolprop: ammonia
+  fprops: Ammonia
+  p: [
+    0.1, 0.2,
+    0.5, 1, 2, 5, 10, 20, 50, 100,
+    200, 500, 1000
+  ]
+  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, 650, 700
+  ]

docs/references.bib

@@ -173,3 +173,32 @@ @article{10.1021/jp0618577
     URL = {https://doi.org/10.1021/jp0618577},
     eprint = {https://doi.org/10.1021/jp0618577}
 }
+
+@incollection{gao2020thermodynamic,
+    title={Thermodynamic properties of ammonia for temperatures from the melting line to 725 K and pressures to 1000 MPa},
+    author={Gao, K and Wu, J and Bell, IH and Lemmon, EW},
+    booktitle={J. Phys. Chem. Ref. Data},
+    year={2020}
+}
+
+@article{tufeu1984thermal,
+    title={Thermal conductivity of ammonia in a large temperature and pressure range including the critical region},
+    author={Tufeu, R and Ivanov, DY and Garrabos, Y and Le Neindre, B},
+    journal={Berichte der Bunsengesellschaft f{\"u}r physikalische Chemie},
+    volume={88},
+    number={4},
+    pages={422--427},
+    year={1984},
+    publisher={Wiley Online Library}
+}
+
+@article{fenghour1995viscosity,
+    title={The viscosity of ammonia},
+    author={Fenghour, A and Wakeham, William A and Vesovic, V and Watson, JTR and Millat, J and Vogel, E},
+    journal={Journal of Physical and Chemical Reference Data},
+    volume={24},
+    number={5},
+    pages={1649--1667},
+    year={1995},
+    publisher={American Institute of Physics for the National Institute of Standards and~…}
+}

include/fprops/ammonia.h

@@ -0,0 +1,44 @@
+// 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 {
+
+/// Ammonia fluid properties
+///
+class Ammonia : public Helmholtz {
+public:
+    Ammonia();
+
+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;
+
+    double lambda_crit(double t, double d) const;
+
+    IdealGasLead<double> lead;
+    IdealGasLogTau<double> log_tau;
+    IdealGasPlanckEinstein<double> pe;
+    ResidualPower<double> power_r;
+    ResidualPowerExp<double, unsigned int> power_exp_r;
+    ResidualGaussian<double> gauss;
+    ResidualGaoB<double> gaob;
+
+    CollisionIntegral<double> eta_0;
+    RainwaterFriend<double> rf;
+    ModifiedBatshinskiHildebrand<double> eta_ho;
+    PolynomialRatio<double> lambda_0;
+    Polynomial<double> lambda_r;
+};
+
+} // namespace fprops

include/fprops/helmholtz.h

@@ -1005,6 +1005,169 @@ class Helmholtz : public SinglePhaseFluidProperties {
         std::vector<T> D;
     };
 
+    template <typename TYPE>
+    class ResidualGaoB {
+    public:
+        ResidualGaoB(const std::vector<double> & n,
+                     const std::vector<double> & d,
+                     const std::vector<double> & t,
+                     const std::vector<double> & eta,
+                     const std::vector<double> & beta,
+                     const std::vector<double> & gamma,
+                     const std::vector<double> & epsilon,
+                     const std::vector<double> & b) :
+            n(n),
+            d(d),
+            t(t),
+            eta(eta),
+            beta(beta),
+            gamma(gamma),
+            epsilon(epsilon),
+            b(b)
+        {
+        }
+
+        TYPE
+        alpha(TYPE delta, TYPE tau) const
+        {
+            TYPE a = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++)
+                a += this->n[i] * f_tau(i, tau) * f_delta(i, delta);
+            return a;
+        }
+
+        TYPE
+        ddelta(TYPE delta, TYPE tau) const
+        {
+            TYPE da_ddelta = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++) {
+                da_ddelta += this->n[i] * f_tau(i, tau) * delta_df_delta_ddelta(i, delta) / delta;
+            }
+            return da_ddelta;
+        }
+
+        TYPE
+        dtau(TYPE delta, TYPE tau) const
+        {
+            TYPE da_dtau = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++) {
+                da_dtau += this->n[i] * f_delta(i, delta) * tau_df_tau_dtau(i, tau) / tau;
+            }
+            return da_dtau;
+        }
+
+        TYPE
+        d2delta(TYPE delta, TYPE tau) const
+        {
+            TYPE d2a_ddelta2 = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++)
+                d2a_ddelta2 += this->n[i] * f_tau(i, tau) * delta2_d2f_delta_ddelta2(i, delta) /
+                               math::pow<2>(delta);
+            return d2a_ddelta2;
+        }
+
+        TYPE
+        d2tau(TYPE delta, TYPE tau) const
+        {
+            TYPE d2a_dtau2 = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++)
+                d2a_dtau2 +=
+                    this->n[i] * f_delta(i, delta) * tau2_d2f_tau_dtau2(i, tau) / math::pow<2>(tau);
+            return d2a_dtau2;
+        }
+
+        TYPE
+        d2deltatau(TYPE delta, TYPE tau) const
+        {
+            TYPE d2a_ddeltatau = 0.;
+            for (std::size_t i = 0; i < this->n.size(); i++)
+                d2a_ddeltatau += this->n[i] * tau_df_tau_dtau(i, tau) *
+                                 delta_df_delta_ddelta(i, delta) / tau / delta;
+            return d2a_ddeltatau;
+        }
+
+    private:
+        TYPE
+        f_tau(std::size_t i, double tau) const
+        {
+            return math::pow(tau, this->t[i]) *
+                   std::exp(1.0 /
+                            (this->b[i] + this->beta[i] * math::pow<2>(-this->gamma[i] + tau)));
+        }
+
+        TYPE
+        f_delta(std::size_t i, double delta) const
+        {
+            return math::pow(delta, this->d[i]) *
+                   std::exp(this->eta[i] * math::pow<2>(delta - this->epsilon[i]));
+        }
+
+        TYPE
+        delta_df_delta_ddelta(std::size_t i, double delta) const
+        {
+            return (this->d[i] * math::pow(delta, this->d[i]) +
+                    2. * math::pow(delta, this->d[i] + 1) * this->eta[i] *
+                        (delta - this->epsilon[i])) *
+                   std::exp(this->eta[i] * math::pow<2>(delta - this->epsilon[i]));
+        }
+
+        TYPE
+        tau_df_tau_dtau(std::size_t i, double tau) const
+        {
+            return (2. * this->beta[i] * math::pow(tau, this->t[i] + 1) * (this->gamma[i] - tau) +
+                    this->t[i] * math::pow(tau, this->t[i]) *
+                        math::pow<2>(this->b[i] +
+                                     this->beta[i] * math::pow<2>(this->gamma[i] - tau))) *
+                   std::exp(1.0 /
+                            (this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau))) /
+                   math::pow<2>(this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau));
+        }
+
+        TYPE
+        delta2_d2f_delta_ddelta2(std::size_t i, double delta) const
+        {
+            return math::pow(delta, this->d[i]) *
+                   (4. * this->d[i] * delta * this->eta[i] * (delta - this->epsilon[i]) +
+                    this->d[i] * (this->d[i] - 1) +
+                    2. * math::pow<2>(delta) * this->eta[i] *
+                        (2. * this->eta[i] * math::pow<2>(delta - this->epsilon[i]) + 1)) *
+                   std::exp(this->eta[i] * math::pow<2>(delta - this->epsilon[i]));
+        }
+
+        TYPE
+        tau2_d2f_tau_dtau2(std::size_t i, double tau) const
+        {
+            return math::pow(tau, this->t[i]) *
+                   (4 * this->beta[i] * this->t[i] * tau *
+                        math::pow<2>(this->b[i] +
+                                     this->beta[i] * math::pow<2>(this->gamma[i] - tau)) *
+                        (this->gamma[i] - tau) +
+                    2 * this->beta[i] * math::pow<2>(tau) *
+                        (4 * this->beta[i] *
+                             (this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau)) *
+                             math::pow<2>(this->gamma[i] - tau) +
+                         2 * this->beta[i] * math::pow<2>(this->gamma[i] - tau) -
+                         math::pow<2>(this->b[i] +
+                                      this->beta[i] * math::pow<2>(this->gamma[i] - tau))) +
+                    this->t[i] *
+                        math::pow(this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau),
+                                  4) *
+                        (this->t[i] - 1)) *
+                   std::exp(1.0 /
+                            (this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau))) /
+                   math::pow(this->b[i] + this->beta[i] * math::pow<2>(this->gamma[i] - tau), 4);
+        }
+
+        std::vector<double> n;
+        std::vector<double> d;
+        std::vector<double> t;
+        std::vector<double> eta;
+        std::vector<double> beta;
+        std::vector<double> gamma;
+        std::vector<double> epsilon;
+        std::vector<double> b;
+    };
+
 protected:
     double delta(double rho) const;
     double tau(double T) const;

python/src/fprops.cpp

@@ -5,6 +5,7 @@
 
 #include "fprops/single_phase_fluid_properties.h"
 #include "fprops/air.h"
+#include "fprops/ammonia.h"
 #include "fprops/ideal_gas.h"
 #include "fprops/nitrogen.h"
 #include "fprops/helium.h"
@@ -48,6 +49,13 @@ PYBIND11_MODULE(fprops, m)
         .def("p_T", &Air::p_T)
         .def("v_u", &Air::v_u);
 
+    py::class_<Ammonia>(m, "Ammonia")
+        .def(py::init())
+        .def("rho_T", &Ammonia::rho_T)
+        .def("rho_p", &Ammonia::rho_p)
+        .def("p_T", &Ammonia::p_T)
+        .def("v_u", &Ammonia::v_u);
+
     py::class_<IdealGas>(m, "IdealGas")
         .def(py::init<double, double>())
         .def("gamma", &IdealGas::gamma)