From 81f6cf91d019b6c112ba013165b2f7530d3d185b Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Fri, 29 Oct 2021 12:49:43 +0100
Subject: [PATCH] liquidProperties: Added properties for NH3 (ammonia)

---
 applications/test/liquid/Test-liquid.C        |  64 +++++-
 .../NSRDS/NSRDS0/NSRDS0.C                     |  20 +-
 .../NSRDS/NSRDS0/NSRDS0.H                     |   6 +-
 .../thermophysicalProperties/Make/files       |   1 +
 .../liquidProperties/NH3/NH3.C                | 138 +++++++++++++
 .../liquidProperties/NH3/NH3.H                | 194 ++++++++++++++++++
 .../liquidProperties/NH3/NH3I.H               | 122 +++++++++++
 7 files changed, 535 insertions(+), 10 deletions(-)
 create mode 100644 src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.C
 create mode 100644 src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.H
 create mode 100644 src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3I.H

diff --git a/applications/test/liquid/Test-liquid.C b/applications/test/liquid/Test-liquid.C
index 43f82fb2de..005dd335e1 100644
--- a/applications/test/liquid/Test-liquid.C
+++ b/applications/test/liquid/Test-liquid.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,20 +25,68 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "C7H8.H"
+#include "argList.H"
+#include "IFstream.H"
+#include "OFstream.H"
+#include "liquidProperties.H"
+
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Main program:
 
-int main()
+int main(int argc, char *argv[])
 {
-    C7H8 fuel;
+    argList::validArgs.append("liquidName");
+    argList::validArgs.append("pMin");
+    argList::validArgs.append("pMax");
+    argList::validArgs.append("nP");
+    argList::validArgs.append("Tmin");
+    argList::validArgs.append("Tmax");
+    argList::validArgs.append("nT");
+    argList args(argc, argv);
+    const word liquidName(args[1]);
 
-    Info<< fuel.rho(1e5, 300) << endl;
-    Info<< fuel << endl;
+    const scalar pMin = args.argRead<scalar>(2);
+    const scalar pMax = args.argRead<scalar>(3);
+    const scalar nP = args.argRead<label>(4);
+    const scalar Tmin = args.argRead<scalar>(5);
+    const scalar Tmax = args.argRead<scalar>(6);
+    const scalar nT = args.argRead<label>(7);
 
-    Info<< "End\n" << endl;
+    autoPtr<liquidProperties> liquidPtr = liquidProperties::New(liquidName);
+
+    OFstream plotFile(liquidName + ".dat");
+
+    plotFile << "# p T rho Cp Hs Ha pv hl Cpg mu mug kappa kappag sigma" << nl;
+
+    for (label pi = 0; pi < nP; ++ pi)
+    {
+        const scalar p = pMin + (pMax - pMin)*pi/(nP - 1);
+
+        for (label Ti = 0; Ti < nT; ++ Ti)
+        {
+            const scalar T = Tmin + (Tmax - Tmin)*Ti/(nT - 1);
+
+            plotFile
+                << p << ' '
+                << T << ' '
+                << liquidPtr->rho(p, T) << ' '
+                << liquidPtr->Cp(p, T) << ' '
+                << liquidPtr->Hs(p, T) << ' '
+                << liquidPtr->Ha(p, T) << ' '
+                << liquidPtr->pv(p, T) << ' '
+                << liquidPtr->hl(p, T) << ' '
+                << liquidPtr->Cpg(p, T) << ' '
+                << liquidPtr->mu(p, T) << ' '
+                << liquidPtr->mug(p, T) << ' '
+                << liquidPtr->kappa(p, T) << ' '
+                << liquidPtr->kappag(p, T) << ' '
+                << liquidPtr->sigma(p, T)
+                << nl;
+        }
+    }
+
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.C b/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.C
index 11bc9c16c2..0827736b5e 100644
--- a/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.C
+++ b/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,24 @@ Foam::scalar Foam::Function1s::NSRDS0::integral
 }
 
 
+Foam::Function1s::NSRDS0 Foam::Function1s::NSRDS0::integral
+(
+    const word& name,
+    const scalar a
+) const
+{
+    if (f_ != 0)
+    {
+        FatalErrorInFunction
+            << "Integral function of " << typeName << " function requested"
+            << " but the \"f\" coefficient is not zero"
+            << exit(FatalError);
+    }
+
+    return NSRDS0(name, a, a_, b_/2, c_/3, d_/4, e_/5);
+}
+
+
 void Foam::Function1s::NSRDS0::write(Ostream& os) const
 {
     writeEntry(os, "a", a_);
diff --git a/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.H b/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.H
index 95ac99702f..249efb2dcb 100644
--- a/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.H
+++ b/src/thermophysicalModels/specie/thermophysicalFunctions/NSRDS/NSRDS0/NSRDS0.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -122,6 +122,10 @@ public:
         //- Integrate between two scalar values
         virtual scalar integral(const scalar x1, const scalar x2) const;
 
+        //- Return the integral as a function. Note that this will error if the
+        //  "f" coefficient is not zero.
+        NSRDS0 integral(const word& name, const scalar a) const;
+
         //- Write the function coefficients
         virtual void write(Ostream& os) const;
 };
diff --git a/src/thermophysicalModels/thermophysicalProperties/Make/files b/src/thermophysicalModels/thermophysicalProperties/Make/files
index 58ea96bef0..0be048910a 100644
--- a/src/thermophysicalModels/thermophysicalProperties/Make/files
+++ b/src/thermophysicalModels/thermophysicalProperties/Make/files
@@ -35,6 +35,7 @@ liquidProperties/MB/MB.C
 liquidProperties/CH4N2O/CH4N2O.C
 liquidProperties/nC3H8O/nC3H8O.C
 liquidProperties/iC3H8O/iC3H8O.C
+liquidProperties/NH3/NH3.C
 
 
 solidProperties/solidProperties/solidProperties.C
diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.C b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.C
new file mode 100644
index 0000000000..1b28f9ef74
--- /dev/null
+++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.C
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "NH3.H"
+#include "thermodynamicConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant::thermodynamic;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(NH3, 0);
+    addToRunTimeSelectionTable(liquidProperties, NH3,);
+    addToRunTimeSelectionTable(liquidProperties, NH3, dictionary);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::NH3::NH3()
+:
+    liquidProperties
+    (
+        typeName,
+        17.030,
+        405.65,
+        1.1278e+07,
+        0.07247,
+        0.242,
+        195.41,
+        6.1177e+03,
+        239.72,
+        4.9034e-30,
+        0.2520,
+        2.9217e+04
+    ),
+    rho_("rho", 3.5430e+00*W(), 2.5471e-01, Tc(), 2.8870e-01),
+    pv_("pv", 9.0451e+01, -4.6690e+03, -1.1601e+01, 1.7183e-02, 1),
+    hl_("hl", Tc(), 3.1523e+07/W(), 3.9140e-01, -2.2890e-01, 2.3090e-01, 0),
+    Cp_
+    (
+        "Cp",
+        1.0827e+06/W(),
+        -1.5541e+04/W(),
+        8.9011e+01/W(),
+        -2.2513e-01/W(),
+        2.1336e-04/W(),
+        0
+    ),
+    h_(Cp_.integral("h", - Cp_.integral("", 0).value(Tstd) - 46190000/W())),
+    Cpg_("Cpg", 3.3190e+04/W(), 7.4230e+04/W(), 5.4040e+02, 8.8950e-01),
+    B_("B", 1.5600e-02, -1.9900e+01, -5.0500e+06, -2.5330e+18, 3.8700e+20),
+    mu_("mu", -1.6430e+00, 4.5560e+02, -1.5637e+00, 0, 0),
+    mug_("mug", 4.1855e-08, 9.8060e-01, 3.0800e+01, 0),
+    kappa_("kappa", 1.1606e+00, -2.2840e-03, 0, 0, 0, 0),
+    kappag_("kappag", -4.5900e-02, 1.6520e-01, -1.7078e+03, 0),
+    sigma_("sigma", 9.1200e-02, 1.1028e+00, 0, 0, 0, 0),
+    D_("D", 14.9, 20.1, W(), 28),
+    Hf_(h_.value(Tstd))
+{}
+
+
+Foam::NH3::NH3
+(
+    const liquidProperties& l,
+    const Function1s::NSRDS5& density,
+    const Function1s::NSRDS1& vapourPressure,
+    const Function1s::NSRDS6& heatOfVapourisation,
+    const Function1s::NSRDS0& heatCapacity,
+    const Function1s::NSRDS0& enthalpy,
+    const Function1s::NSRDS3& idealGasHeatCapacity,
+    const Function1s::NSRDS4& secondVirialCoeff,
+    const Function1s::NSRDS1& dynamicViscosity,
+    const Function1s::NSRDS2& vapourDynamicViscosity,
+    const Function1s::NSRDS0& thermalConductivity,
+    const Function1s::NSRDS2& vapourThermalConductivity,
+    const Function1s::NSRDS6& surfaceTension,
+    const Function2s::APIdiffCoef& vapourDiffusivity
+)
+:
+    liquidProperties(l),
+    rho_(density),
+    pv_(vapourPressure),
+    hl_(heatOfVapourisation),
+    Cp_(heatCapacity),
+    h_(enthalpy),
+    Cpg_(idealGasHeatCapacity),
+    B_(secondVirialCoeff),
+    mu_(dynamicViscosity),
+    mug_(vapourDynamicViscosity),
+    kappa_(thermalConductivity),
+    kappag_(vapourThermalConductivity),
+    sigma_(surfaceTension),
+    D_(vapourDiffusivity),
+    Hf_(h_.value(Tstd))
+{}
+
+
+Foam::NH3::NH3(const dictionary& dict)
+:
+    NH3()
+{
+    readIfPresent(*this, dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::NH3::write(Ostream& os) const
+{
+    liquidProperties::write(*this, os);
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.H b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.H
new file mode 100644
index 0000000000..704f08a7be
--- /dev/null
+++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3.H
@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::NH3
+
+Description
+    Liquid ammonia
+
+SourceFiles
+    NH3.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NH3_H
+#define NH3_H
+
+#include "liquidProperties.H"
+#include "NSRDS0.H"
+#include "NSRDS1.H"
+#include "NSRDS2.H"
+#include "NSRDS3.H"
+#include "NSRDS4.H"
+#include "NSRDS5.H"
+#include "NSRDS6.H"
+#include "NSRDS7.H"
+#include "APIdiffCoef.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                             Class NH3 Declaration
+\*---------------------------------------------------------------------------*/
+
+class NH3
+:
+    public liquidProperties
+{
+    // Private Data
+
+        Function1s::NSRDS5 rho_;
+        Function1s::NSRDS1 pv_;
+        Function1s::NSRDS6 hl_;
+        Function1s::NSRDS0 Cp_;
+        Function1s::NSRDS0 h_;
+        Function1s::NSRDS3 Cpg_;
+        Function1s::NSRDS4 B_;
+        Function1s::NSRDS1 mu_;
+        Function1s::NSRDS2 mug_;
+        Function1s::NSRDS0 kappa_;
+        Function1s::NSRDS2 kappag_;
+        Function1s::NSRDS6 sigma_;
+        Function2s::APIdiffCoef D_;
+
+        //- Liquid heat of formation [J/kg]
+        scalar Hf_;
+
+
+public:
+
+    friend class liquidProperties;
+
+    //- Runtime type information
+    TypeName("NH3");
+
+
+    // Constructors
+
+        //- Construct null
+        NH3();
+
+        //- Construct from components
+        NH3
+        (
+            const liquidProperties& l,
+            const Function1s::NSRDS5& density,
+            const Function1s::NSRDS1& vapourPressure,
+            const Function1s::NSRDS6& heatOfVapourisation,
+            const Function1s::NSRDS0& heatCapacity,
+            const Function1s::NSRDS0& enthalpy,
+            const Function1s::NSRDS3& idealGasHeatCapacity,
+            const Function1s::NSRDS4& secondVirialCoeff,
+            const Function1s::NSRDS1& dynamicViscosity,
+            const Function1s::NSRDS2& vapourDynamicViscosity,
+            const Function1s::NSRDS0& thermalConductivity,
+            const Function1s::NSRDS2& vapourThermalConductivity,
+            const Function1s::NSRDS6& surfaceTension,
+            const Function2s::APIdiffCoef& vapourDiffusivity
+        );
+
+        //- Construct from dictionary
+        NH3(const dictionary& dict);
+
+        //- Construct and return clone
+        virtual autoPtr<liquidProperties> clone() const
+        {
+            return autoPtr<liquidProperties>(new NH3(*this));
+        }
+
+
+    // Member Functions
+
+        //- Liquid density [kg/m^3]
+        inline scalar rho(scalar p, scalar T) const;
+
+        //- Vapour pressure [Pa]
+        inline scalar pv(scalar p, scalar T) const;
+
+        //- Heat of vapourisation [J/kg]
+        inline scalar hl(scalar p, scalar T) const;
+
+        //- Liquid heat capacity [J/kg/K]
+        inline scalar Cp(scalar p, scalar T) const;
+
+        //- Liquid sensible enthalpy [J/kg]
+        inline scalar Hs(scalar p, scalar T) const;
+
+        //- Liquid heat of formation [J/kg]
+        inline scalar Hf() const;
+
+        //- Liquid absolute enthalpy [J/kg]
+        inline scalar Ha(scalar p, scalar T) const;
+
+        //- Ideal gas heat capacity [J/kg/K]
+        inline scalar Cpg(scalar p, scalar T) const;
+
+        //- Second Virial Coefficient [m^3/kg]
+        inline scalar B(scalar p, scalar T) const;
+
+        //- Liquid viscosity [Pa s]
+        inline scalar mu(scalar p, scalar T) const;
+
+        //- Vapour viscosity [Pa s]
+        inline scalar mug(scalar p, scalar T) const;
+
+        //- Liquid thermal conductivity  [W/m/K]
+        inline scalar kappa(scalar p, scalar T) const;
+
+        //- Vapour thermal conductivity  [W/m/K]
+        inline scalar kappag(scalar p, scalar T) const;
+
+        //- Surface tension [N/m]
+        inline scalar sigma(scalar p, scalar T) const;
+
+        //- Vapour diffusivity [m^2/s]
+        inline scalar D(scalar p, scalar T) const;
+
+        //- Vapour diffusivity [m^2/s] with specified binary pair
+        inline scalar D(scalar p, scalar T, scalar Wb) const;
+
+
+    // I-O
+
+        //- Write the function coefficients
+        void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NH3I.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3I.H b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3I.H
new file mode 100644
index 0000000000..d0bd0f1acc
--- /dev/null
+++ b/src/thermophysicalModels/thermophysicalProperties/liquidProperties/NH3/NH3I.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+inline Foam::scalar Foam::NH3::rho(scalar p, scalar T) const
+{
+    return rho_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::pv(scalar p, scalar T) const
+{
+    return pv_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::hl(scalar p, scalar T) const
+{
+    return hl_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::Cp(scalar p, scalar T) const
+{
+    return Cp_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::Hs(scalar p, scalar T) const
+{
+    return Ha(p, T) - Hf();
+}
+
+
+inline Foam::scalar Foam::NH3::Hf() const
+{
+    return Hf_;
+}
+
+
+inline Foam::scalar Foam::NH3::Ha(scalar p, scalar T) const
+{
+    return h_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::Cpg(scalar p, scalar T) const
+{
+    return Cpg_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::B(scalar p, scalar T) const
+{
+    return B_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::mu(scalar p, scalar T) const
+{
+    return mu_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::mug(scalar p, scalar T) const
+{
+    return mug_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::kappa(scalar p, scalar T) const
+{
+    return kappa_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::kappag(scalar p, scalar T) const
+{
+    return kappag_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::sigma(scalar p, scalar T) const
+{
+    return sigma_.value(T);
+}
+
+
+inline Foam::scalar Foam::NH3::D(scalar p, scalar T) const
+{
+    return D_.value(p, T);
+}
+
+
+inline Foam::scalar Foam::NH3::D(scalar p, scalar T, scalar Wb) const
+{
+    return D_.value(p, T, Wb);
+}
+
+
+// ************************************************************************* //