chtMultiRegionFoam: Changed solid energy from enthalpy to internal energy

The solid is currently assumed incompressible (the solid pressure is not
updated) and in general would be near incompressible so internal energy is a
more appropriate energy choice than enthalpy which would require a pressure work
term currently not implemented.  Additionally due to the way in which the
conduction is handled in terms of the gradient of energy the accuracy of the
current enthalpy implementation is sensitive to the pressure distribution as
this introduces an enthalpy gradient from the p/rho term which would need to be
corrected; this issue is avoided by solving for internal energy instead.

This improvement requires the scheme and solver settings for the solids in
chtMultiRegionFoam cases to be changed from "h" to "e" and the thermo-physical
properties in <solid>/thermophysicalProperties to be set to the corresponding
internal energy forms, e.g.:

    thermo          eConst;
    .
    .
    .
    energy          sensibleInternalEnergy;

All tutorials have be updated to reflect this and provide guidance when updating
cases.

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H

@@ -34,8 +34,8 @@ forAll(solidRegions, i)
             coordinateSystem::New(solidRegions[i], thermos[i])
         );
 
-        tmp<volVectorField> tkappaByCp =
+        tmp<volVectorField> tkappaByCv =
-            thermos[i].Kappa()/thermos[i].Cp();
+            thermos[i].Kappa()/thermos[i].Cv();
 
         aniAlphas.set
         (
@@ -54,7 +54,7 @@ forAll(solidRegions, i)
                 dimensionedSymmTensor
                 (
                     "zero",
-                    tkappaByCp().dimensions(),
+                    tkappaByCv().dimensions(),
                     Zero
                 ),
                 zeroGradientFvPatchSymmTensorField::typeName
@@ -62,7 +62,7 @@ forAll(solidRegions, i)
         );
 
         aniAlphas[i].primitiveFieldRef() =
-            coordinates[i].R().transformVector(tkappaByCp());
+            coordinates[i].R().transformVector(tkappaByCv());
         aniAlphas[i].correctBoundaryConditions();
 
     }

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/setRegionSolidFields.H

@@ -3,18 +3,18 @@ solidThermo& thermo = thermos[i];
 tmp<volScalarField> trho = thermo.rho();
 const volScalarField& rho = trho();
 
-tmp<volScalarField> tcp = thermo.Cp();
+tmp<volScalarField> tCv = thermo.Cv();
-const volScalarField& cp = tcp();
+const volScalarField& Cv = tCv();
 
 tmp<volSymmTensorField> taniAlpha;
 if (!thermo.isotropic())
 {
     volSymmTensorField& aniAlpha = aniAlphas[i];
-    tmp<volVectorField> tkappaByCp = thermo.Kappa()/cp;
+    tmp<volVectorField> tkappaByCv = thermo.Kappa()/Cv;
     const coordinateSystem& coodSys = coordinates[i];
 
     aniAlpha.primitiveFieldRef() =
-        coodSys.R().transformVector(tkappaByCp());
+        coodSys.R().transformVector(tkappaByCv());
     aniAlpha.correctBoundaryConditions();
 
     taniAlpha = tmp<volSymmTensorField>
@@ -24,7 +24,7 @@ if (!thermo.isotropic())
 }
 
 
-volScalarField& h = thermo.he();
+volScalarField& e = thermo.he();
 
 const volScalarField& betav = betavSolid[i];
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H

@@ -1,25 +1,25 @@
 {
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix hEqn
+        fvScalarMatrix eEqn
         (
-            fvm::ddt(betav*rho, h)
+            fvm::ddt(betav*rho, e)
           - (
                 thermo.isotropic()
-              ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+              ? fvm::laplacian(betav*thermo.alpha(), e, "laplacian(alpha,e)")
-              : fvm::laplacian(betav*taniAlpha(), h, "laplacian(alpha,h)")
+              : fvm::laplacian(betav*taniAlpha(), e, "laplacian(alpha,e)")
             )
           ==
-            fvOptions(rho, h)
+            fvOptions(rho, e)
         );
 
-        hEqn.relax();
+        eEqn.relax();
 
-        fvOptions.constrain(hEqn);
+        fvOptions.constrain(eEqn);
 
-        hEqn.solve();
+        eEqn.solve();
 
-        fvOptions.correct(h);
+        fvOptions.correct(e);
     }
 }
 

etc/templates/singleFluidCHT/templates/materials/aluminium/thermophysicalProperties

@@ -19,10 +19,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -41,7 +41,7 @@ mixture
     thermodynamics
     {
         Hf      0;
-        Cp      896;    // [J/kg/K]
+        Cv      896;    // [J/kg/K]
     }
 
     equationOfState

etc/templates/singleFluidCHT/templates/materials/copper/thermophysicalProperties

@@ -19,10 +19,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -41,7 +41,7 @@ mixture
     thermodynamics
     {
         Hf      0;
-        Cp      385;        // [J/kg/K]
+        Cv      385;        // [J/kg/K]
     }
 
     equationOfState

src/ThermophysicalTransportModels/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarFields.C

@@ -24,9 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "thermalBaffle1DFvPatchScalarField.H"
-
 #include "forSolids.H"
-
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -39,22 +37,22 @@ typedef
     <
         species::thermo
         <
-            hConstThermo
+            eConstThermo
             <
                 rhoConst<specie>
             >,
-            sensibleEnthalpy
+            sensibleInternalEnergy
         >
-    > hConstSolidThermoPhysics;
+    > eConstSolidThermoPhysics;
 
 typedef
-    compressible::thermalBaffle1DFvPatchScalarField<hConstSolidThermoPhysics>
+    compressible::thermalBaffle1DFvPatchScalarField<eConstSolidThermoPhysics>
     thermalBaffle1DHConstSolidThermoPhysicsFvPatchScalarField;
 
 defineTemplateTypeNameAndDebugWithName
 (
     thermalBaffle1DHConstSolidThermoPhysicsFvPatchScalarField,
-    "compressible::thermalBaffle1D<hConstSolidThermoPhysics>",
+    "compressible::thermalBaffle1D<eConstSolidThermoPhysics>",
     0
 );
 
@@ -69,22 +67,22 @@ typedef
     <
         species::thermo
         <
-            hPowerThermo
+            ePowerThermo
             <
                 rhoConst<specie>
             >,
-            sensibleEnthalpy
+            sensibleInternalEnergy
         >
-    > hPowerSolidThermoPhysics;
+    > ePowerSolidThermoPhysics;
 
 typedef
-    compressible::thermalBaffle1DFvPatchScalarField<hPowerSolidThermoPhysics>
+    compressible::thermalBaffle1DFvPatchScalarField<ePowerSolidThermoPhysics>
     thermalBaffle1DHPowerSolidThermoPhysicsFvPatchScalarField;
 
 defineTemplateTypeNameAndDebugWithName
 (
     thermalBaffle1DHPowerSolidThermoPhysicsFvPatchScalarField,
-    "compressible::thermalBaffle1D<hPowerSolidThermoPhysics>",
+    "compressible::thermalBaffle1D<ePowerSolidThermoPhysics>",
     0
 );
 

src/thermophysicalModels/solidSpecie/include/forSolids.H

@@ -30,11 +30,11 @@ License
 
 #include "rhoConst.H"
 
-#include "hConstThermo.H"
+#include "eConstThermo.H"
-#include "hPolynomialThermo.H"
+#include "ePolynomialThermo.H"
-#include "hPowerThermo.H"
+#include "ePowerThermo.H"
 
-#include "sensibleEnthalpy.H"
+#include "sensibleInternalEnergy.H"
 
 #include "constIsoSolidTransport.H"
 #include "constAnIsoSolidTransport.H"
@@ -51,9 +51,10 @@ License
     forThermo(Mu, He, Cp, rhoConst, specie, Macro, Args)
 
 #define forSolidEnergiesAndThermos(Mu, Macro, Args...)                         \
-    forSolidEquations(Mu, sensibleEnthalpy, hConstThermo, Macro, Args);        \
+    forSolidEquations(Mu, sensibleInternalEnergy, eConstThermo, Macro, Args);  \
-    forSolidEquations(Mu, sensibleEnthalpy, hPolynomialThermo, Macro, Args);   \
+    forSolidEquations                                                          \
-    forSolidEquations(Mu, sensibleEnthalpy, hPowerThermo, Macro, Args)
+        (Mu, sensibleInternalEnergy, ePolynomialThermo, Macro, Args);          \
+    forSolidEquations(Mu, sensibleInternalEnergy, ePowerThermo, Macro, Args)
 
 #define forSolidTransports(Macro, Args...)                                     \
     forSolidEnergiesAndThermos(constIsoSolidTransport, Macro, Args);           \

src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H

@@ -141,7 +141,8 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::S
     const scalar T
 ) const
 {
-    return Cp(p, T)*log(T/Tstd) + EquationOfState::S(p, T);
+    return Cp(p, T)*log(T/Tstd)
+      + EquationOfState::S(p, T); // Requires EquationOfState::S for Cv
 }
 
 
@@ -153,7 +154,7 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::Gstd
 {
     return
         Cv_*(T - Tref_) + Esref_ + Hf() + Pstd/EquationOfState::rho(Pstd, T)
-      - Cp(Pstd, T)*T*log(T/Tstd);
+      - S(Pstd, T)*T;
 }
 
 
@@ -164,7 +165,7 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::dCpdT
     const scalar T
 ) const
 {
-    return 0;
+    return 0; // EquationOfState::dCpdT
 }
 
 

src/thermophysicalModels/specie/thermo/ePolynomial/ePolynomialThermo.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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 "ePolynomialThermo.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+Foam::ePolynomialThermo<EquationOfState, PolySize>::ePolynomialThermo
+(
+    const dictionary& dict
+)
+:
+    EquationOfState(dict),
+    Hf_(dict.subDict("thermodynamics").lookup<scalar>("Hf")),
+    Sf_(dict.subDict("thermodynamics").lookup<scalar>("Sf")),
+    CvCoeffs_
+    (
+        dict.subDict("thermodynamics").lookup
+        (
+            "CvCoeffs<" + Foam::name(PolySize) + '>'
+        )
+    ),
+    eCoeffs_(),
+    sCoeffs_()
+{
+    eCoeffs_ = CvCoeffs_.integral();
+    sCoeffs_ = CvCoeffs_.integralMinus1();
+
+    // Offset e poly so that it is relative to the enthalpy at Tstd
+    eCoeffs_[0] -= eCoeffs_.value(Tstd);
+
+    // Offset s poly so that it is relative to the entropy at Tstd
+    sCoeffs_[0] -= sCoeffs_.value(Tstd);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+void Foam::ePolynomialThermo<EquationOfState, PolySize>::write
+(
+    Ostream& os
+) const
+{
+    EquationOfState::write(os);
+
+    dictionary dict("thermodynamics");
+    dict.add("Hf", Hf_);
+    dict.add("Sf", Sf_);
+    dict.add
+    (
+        word("CvCoeffs<" + Foam::name(PolySize) + '>'),
+        CvCoeffs_
+    );
+    os  << indent << dict.dictName() << dict;
+}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const ePolynomialThermo<EquationOfState, PolySize>& pt
+)
+{
+    pt.write(os);
+    return os;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/thermo/ePolynomial/ePolynomialThermo.H

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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::ePolynomialThermo
+
+Description
+    Thermodynamics package templated on the equation of state, using polynomial
+    functions for \c cv, \c e and \c s.
+
+    Polynomials for \c e and \c s derived from \c cv.
+
+Usage
+
+    \table
+        Property     | Description
+        Hf           | Heat of formation
+        Sf           | Standard entropy
+        CvCoeffs<8>  | Specific heat at constant volume polynomial coeffs
+    \endtable
+
+    Example of the specification of the thermodynamic properties:
+    \verbatim
+    thermodynamics
+    {
+        Hf              0;
+        Sf              0;
+        CvCoeffs<8>     ( 1000 -0.05 0.003 0 0 0 0 0 );
+    }
+    \endverbatim
+
+    The polynomial expression is evaluated as so:
+
+        \f[
+            Cv = 1000 - 0.05 T + 0.003 T^2
+        \f]
+
+Note
+    - Specific heat at constant volume polynomial coefficients evaluate to an
+      expression in [J/kg/K].
+
+SourceFiles
+    ePolynomialThermoI.H
+    ePolynomialThermo.C
+
+See also
+    Foam::Polynomial
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ePolynomialThermo_H
+#define ePolynomialThermo_H
+
+#include "scalar.H"
+#include "Polynomial.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class EquationOfState, int PolySize>
+class ePolynomialThermo;
+
+template<class EquationOfState, int PolySize>
+inline ePolynomialThermo<EquationOfState, PolySize> operator+
+(
+    const ePolynomialThermo<EquationOfState, PolySize>&,
+    const ePolynomialThermo<EquationOfState, PolySize>&
+);
+
+template<class EquationOfState, int PolySize>
+inline ePolynomialThermo<EquationOfState, PolySize> operator*
+(
+    const scalar,
+    const ePolynomialThermo<EquationOfState, PolySize>&
+);
+
+template<class EquationOfState, int PolySize>
+inline ePolynomialThermo<EquationOfState, PolySize> operator==
+(
+    const ePolynomialThermo<EquationOfState, PolySize>&,
+    const ePolynomialThermo<EquationOfState, PolySize>&
+);
+
+template<class EquationOfState, int PolySize>
+Ostream& operator<<
+(
+    Ostream&,
+    const ePolynomialThermo<EquationOfState, PolySize>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                      Class ePolynomialThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class EquationOfState, int PolySize=8>
+class ePolynomialThermo
+:
+    public EquationOfState
+{
+    // Private Data
+
+        //- Heat of formation
+        scalar Hf_;
+
+        //- Standard entropy
+        scalar Sf_;
+
+        //- Specific heat at constant volume polynomial coeffs [J/kg/K/K^i]
+        Polynomial<PolySize> CvCoeffs_;
+
+        //- Internal energy polynomial coeffs [J/kg/K^i]
+        //  Derived from Cv coeffs. Relative to Tstd.
+        typename Polynomial<PolySize>::intPolyType eCoeffs_;
+
+        //- Entropy polynomial coeffs [J/kg/K/K^i]
+        //  Derived from Cv coeffs. Relative to Tstd.
+        Polynomial<PolySize> sCoeffs_;
+
+
+    // Private Member Functions
+
+        //- Construct from components
+        inline ePolynomialThermo
+        (
+            const EquationOfState& pt,
+            const scalar Hf,
+            const scalar Sf,
+            const Polynomial<PolySize>& CvCoeffs,
+            const typename Polynomial<PolySize>::intPolyType& eCoeffs,
+            const Polynomial<PolySize>& sCoeffs
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        ePolynomialThermo(const dictionary& dict);
+
+        //- Construct as a named copy
+        inline ePolynomialThermo(const word&, const ePolynomialThermo&);
+
+
+    // Member Functions
+
+        //- Return the instantiated type name
+        static word typeName()
+        {
+            return "ePolynomial<" + EquationOfState::typeName() + '>';
+        }
+
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar) const;
+
+
+        // Fundamental properties
+
+            //- Heat capacity at constant volume [J/kg/K]
+            inline scalar Cv(const scalar p, const scalar T) const;
+
+            //- Sensible internal energy [J/kg]
+            inline scalar Es(const scalar p, const scalar T) const;
+
+            //- Absolute internal energy [J/kg]
+            inline scalar Ea(const scalar p, const scalar T) const;
+
+            //- Enthalpy of formation [J/kg]
+            inline scalar Hf() const;
+
+            //- Entropy [J/kg/K]
+            inline scalar S(const scalar p, const scalar T) const;
+
+            //- Gibbs free energy of the mixture in the standard state [J/kg]
+            inline scalar Gstd(const scalar T) const;
+
+            #include "EtoHthermo.H"
+
+
+        // Derivative term used for Jacobian
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
+        // I-O
+
+            //- Write to Ostream
+            void write(Ostream& os) const;
+
+
+    // Member Operators
+
+        inline void operator+=(const ePolynomialThermo&);
+        inline void operator*=(const scalar);
+
+
+    // Friend operators
+
+        friend ePolynomialThermo operator+ <EquationOfState, PolySize>
+        (
+            const ePolynomialThermo&,
+            const ePolynomialThermo&
+        );
+
+        friend ePolynomialThermo operator* <EquationOfState, PolySize>
+        (
+            const scalar,
+            const ePolynomialThermo&
+        );
+
+        friend ePolynomialThermo operator== <EquationOfState, PolySize>
+        (
+            const ePolynomialThermo&,
+            const ePolynomialThermo&
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <EquationOfState, PolySize>
+        (
+            Ostream&,
+            const ePolynomialThermo&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ePolynomialThermoI.H"
+
+#ifdef NoRepository
+    #include "ePolynomialThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/thermo/ePolynomial/ePolynomialThermoI.H

@@ -0,0 +1,280 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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 "ePolynomialThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+inline Foam::ePolynomialThermo<EquationOfState, PolySize>::ePolynomialThermo
+(
+    const EquationOfState& pt,
+    const scalar Hf,
+    const scalar Sf,
+    const Polynomial<PolySize>& CvCoeffs,
+    const typename Polynomial<PolySize>::intPolyType& eCoeffs,
+    const Polynomial<PolySize>& sCoeffs
+)
+:
+    EquationOfState(pt),
+    Hf_(Hf),
+    Sf_(Sf),
+    CvCoeffs_(CvCoeffs),
+    eCoeffs_(eCoeffs),
+    sCoeffs_(sCoeffs)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+inline Foam::ePolynomialThermo<EquationOfState, PolySize>::ePolynomialThermo
+(
+    const word& name,
+    const ePolynomialThermo& pt
+)
+:
+    EquationOfState(name, pt),
+    Hf_(pt.Hf_),
+    Sf_(pt.Sf_),
+    CvCoeffs_(pt.CvCoeffs_),
+    eCoeffs_(pt.eCoeffs_),
+    sCoeffs_(pt.sCoeffs_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::Cv
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return CvCoeffs_.value(T) + EquationOfState::Cv(p, T);
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::Es
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return eCoeffs_.value(T) + EquationOfState::E(p, T);
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::Ea
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return Es(p, T) + Hf();
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::Hf()
+const
+{
+    return Hf_;
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::S
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return sCoeffs_.value(T)
+      + EquationOfState::S(p, T); // Requires EquationOfState::S for Cv
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::Gstd
+(
+    const scalar T
+) const
+{
+    return
+        eCoeffs_.value(T) + Hf() + Pstd/EquationOfState::rho(Pstd, T)
+      - S(Pstd, T)*T;
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::ePolynomialThermo<EquationOfState, PolySize>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return CvCoeffs_.derivative(T); // + EquationOfState::dCpdT
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+inline void Foam::ePolynomialThermo<EquationOfState, PolySize>::operator+=
+(
+    const ePolynomialThermo<EquationOfState, PolySize>& pt
+)
+{
+    scalar Y1 = this->Y();
+
+    EquationOfState::operator+=(pt);
+
+    if (mag(this->Y()) > small)
+    {
+        Y1 /= this->Y();
+        const scalar Y2 = pt.Y()/this->Y();
+
+        Hf_ = Y1*Hf_ + Y2*pt.Hf_;
+        Sf_ = Y1*Sf_ + Y2*pt.Sf_;
+        CvCoeffs_ = Y1*CvCoeffs_ + Y2*pt.CvCoeffs_;
+        eCoeffs_ = Y1*eCoeffs_ + Y2*pt.eCoeffs_;
+        sCoeffs_ = Y1*sCoeffs_ + Y2*pt.sCoeffs_;
+    }
+}
+
+
+template<class EquationOfState, int PolySize>
+inline void Foam::ePolynomialThermo<EquationOfState, PolySize>::operator*=
+(
+    const scalar s
+)
+{
+    EquationOfState::operator*=(s);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class EquationOfState, int PolySize>
+inline Foam::ePolynomialThermo<EquationOfState, PolySize> Foam::operator+
+(
+    const ePolynomialThermo<EquationOfState, PolySize>& pt1,
+    const ePolynomialThermo<EquationOfState, PolySize>& pt2
+)
+{
+    EquationOfState eofs = pt1;
+    eofs += pt2;
+
+    if (mag(eofs.Y()) < small)
+    {
+        return ePolynomialThermo<EquationOfState, PolySize>
+        (
+            eofs,
+            pt1.Hf_,
+            pt1.Sf_,
+            pt1.CvCoeffs_,
+            pt1.eCoeffs_,
+            pt1.sCoeffs_
+        );
+    }
+    {
+        const scalar Y1 = pt1.Y()/eofs.Y();
+        const scalar Y2 = pt2.Y()/eofs.Y();
+
+        return ePolynomialThermo<EquationOfState, PolySize>
+        (
+            eofs,
+            Y1*pt1.Hf_ + Y2*pt2.Hf_,
+            Y1*pt1.Sf_ + Y2*pt2.Sf_,
+            Y1*pt1.CvCoeffs_ + Y2*pt2.CvCoeffs_,
+            Y1*pt1.eCoeffs_ + Y2*pt2.eCoeffs_,
+            Y1*pt1.sCoeffs_ + Y2*pt2.sCoeffs_
+        );
+    }
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::ePolynomialThermo<EquationOfState, PolySize> Foam::operator*
+(
+    const scalar s,
+    const ePolynomialThermo<EquationOfState, PolySize>& pt
+)
+{
+    return ePolynomialThermo<EquationOfState, PolySize>
+    (
+        s*static_cast<const EquationOfState&>(pt),
+        pt.Hf_,
+        pt.Sf_,
+        pt.CvCoeffs_,
+        pt.eCoeffs_,
+        pt.sCoeffs_
+    );
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::ePolynomialThermo<EquationOfState, PolySize> Foam::operator==
+(
+    const ePolynomialThermo<EquationOfState, PolySize>& pt1,
+    const ePolynomialThermo<EquationOfState, PolySize>& pt2
+)
+{
+    EquationOfState eofs
+    (
+        static_cast<const EquationOfState&>(pt1)
+     == static_cast<const EquationOfState&>(pt2)
+    );
+
+    const scalar Y1 = pt1.Y()/eofs.Y();
+    const scalar Y2 = pt2.Y()/eofs.Y();
+
+    return ePolynomialThermo<EquationOfState, PolySize>
+    (
+        eofs,
+        Y2*pt2.Hf_       - Y1*pt1.Hf_,
+        Y2*pt2.Sf_       - Y1*pt1.Sf_,
+        Y2*pt2.CvCoeffs_ - Y1*pt1.CvCoeffs_,
+        Y2*pt2.eCoeffs_  - Y1*pt1.eCoeffs_,
+        Y2*pt2.sCoeffs_  - Y1*pt1.sCoeffs_
+    );
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/thermo/ePower/ePowerThermo.C

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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 "ePowerThermo.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+Foam::ePowerThermo<EquationOfState>::ePowerThermo
+(
+    const dictionary& dict
+)
+:
+    EquationOfState(dict),
+    c0_(dict.subDict("thermodynamics").lookup<scalar>("C0")),
+    n0_(dict.subDict("thermodynamics").lookup<scalar>("n0")),
+    Tref_(dict.subDict("thermodynamics").lookup<scalar>("Tref")),
+    Hf_(dict.subDict("thermodynamics").lookup<scalar>("Hf"))
+{}
+
+
+// * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const ePowerThermo<EquationOfState>& et
+)
+{
+    et.write(os);
+    return os;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/thermo/ePower/ePowerThermo.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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::ePowerThermo
+
+Description
+    Power-function based thermodynamics package templated on EquationOfState.
+
+    In this thermodynamics package the heat capacity is a simple power of
+    temperature:
+
+        Cv(T) = c0*(T/Tref)^n0;
+
+    which is particularly suitable for solids.
+
+SourceFiles
+    ePowerThermoI.H
+    ePowerThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ePowerThermo_H
+#define ePowerThermo_H
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class EquationOfState> class ePowerThermo;
+
+template<class EquationOfState>
+inline ePowerThermo<EquationOfState> operator+
+(
+    const ePowerThermo<EquationOfState>&,
+    const ePowerThermo<EquationOfState>&
+);
+
+template<class EquationOfState>
+inline ePowerThermo<EquationOfState> operator*
+(
+    const scalar,
+    const ePowerThermo<EquationOfState>&
+);
+
+
+template<class EquationOfState>
+inline ePowerThermo<EquationOfState> operator==
+(
+    const ePowerThermo<EquationOfState>&,
+    const ePowerThermo<EquationOfState>&
+);
+
+
+template<class EquationOfState>
+Ostream& operator<<
+(
+    Ostream&,
+    const ePowerThermo<EquationOfState>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class ePowerThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class EquationOfState>
+class ePowerThermo
+:
+    public EquationOfState
+{
+    // Private Data
+
+        scalar c0_;
+        scalar n0_;
+        scalar Tref_;
+        scalar Hf_;
+
+
+    // Private Member Functions
+
+        //- Check given temperature is within the range of the fitted coeffs
+        inline void checkT(const scalar T) const;
+
+        //- Construct from components
+        inline ePowerThermo
+        (
+            const EquationOfState& st,
+            const scalar c0,
+            const scalar n0,
+            const scalar Tref,
+            const scalar Hf
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        ePowerThermo(const dictionary&);
+
+        //- Construct as a named copy
+        inline ePowerThermo
+        (
+            const word&,
+            const ePowerThermo&
+        );
+
+         //- Construct and return a clone
+        inline autoPtr<ePowerThermo> clone() const;
+
+        //- Selector from dictionary
+        inline static autoPtr<ePowerThermo> New(const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Return the instantiated type name
+        static word typeName()
+        {
+            return "ePower<" + EquationOfState::typeName() + '>';
+        }
+
+        //- Limit the temperature to be in the range Tlow_ to Thigh_
+        inline scalar limit(const scalar T) const;
+
+
+        // Fundamental properties
+
+            //- Heat capacity at constant volume [J/kg/K]
+            inline scalar Cv(const scalar p, const scalar T) const;
+
+            //- Sensible internal energy [J/kg]
+            inline scalar Es(const scalar p, const scalar T) const;
+
+            //- Absolute internal energy [J/kg]
+            inline scalar Ea(const scalar p, const scalar T) const;
+
+            //- Enthalpy of formation [J/kg]
+            inline scalar Hf() const;
+
+            //- Entropy [J/kg/K]
+            inline scalar S(const scalar p, const scalar T) const;
+
+            //- Gibbs free energy of the mixture in the standard state [J/kg]
+            inline scalar Gstd(const scalar T) const;
+
+            #include "EtoHthermo.H"
+
+
+        // Derivative term used for Jacobian
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
+    // Member Operators
+
+        inline void operator+=(const ePowerThermo&);
+
+
+    // Friend operators
+
+        friend ePowerThermo operator+ <EquationOfState>
+        (
+            const ePowerThermo&,
+            const ePowerThermo&
+        );
+
+        friend ePowerThermo operator* <EquationOfState>
+        (
+            const scalar,
+            const ePowerThermo&
+        );
+
+
+        friend ePowerThermo operator== <EquationOfState>
+        (
+            const ePowerThermo&,
+            const ePowerThermo&
+        );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <EquationOfState>
+        (
+            Ostream&,
+            const ePowerThermo&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ePowerThermoI.H"
+    #include "ePowerThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/thermo/ePower/ePowerThermoI.H

@@ -0,0 +1,312 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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 "ePowerThermo.H"
+#include "specie.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class EquationOfState>
+inline void Foam::ePowerThermo<EquationOfState>::checkT
+(
+    const scalar T
+) const
+{
+    if (T < 0)
+    {
+        FatalErrorInFunction
+            << "attempt to evaluate ePowerThermo<EquationOfState>"
+               " for negative temperature " << T
+            << abort(FatalError);
+    }
+}
+
+
+template<class EquationOfState>
+inline Foam::ePowerThermo<EquationOfState>::ePowerThermo
+(
+    const word& name,
+    const ePowerThermo& jt
+)
+:
+    EquationOfState(name, jt),
+    c0_(jt.c0_),
+    n0_(jt.n0_),
+    Tref_(jt.Tref_),
+    Hf_(jt.Hf_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+inline Foam::ePowerThermo<EquationOfState>::ePowerThermo
+(
+    const EquationOfState& st,
+    const scalar c0,
+    const scalar n0,
+    const scalar Tref,
+    const scalar Hf
+)
+:
+    EquationOfState(st),
+    c0_(c0),
+    n0_(n0),
+    Tref_(Tref),
+    Hf_(Hf)
+{}
+
+
+template<class EquationOfState>
+inline Foam::autoPtr<Foam::ePowerThermo<EquationOfState>>
+Foam::ePowerThermo<EquationOfState>::clone() const
+{
+    return autoPtr<ePowerThermo<EquationOfState>>
+    (
+        new ePowerThermo<EquationOfState>(*this)
+    );
+}
+
+
+template<class EquationOfState>
+inline Foam::autoPtr<Foam::ePowerThermo<EquationOfState>>
+Foam::ePowerThermo<EquationOfState>::New(const dictionary& dict)
+{
+    return autoPtr<ePowerThermo<EquationOfState>>
+    (
+        new ePowerThermo<EquationOfState>(dict)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::limit
+(
+    const scalar T
+) const
+{
+    return T;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::Cv
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return c0_*pow(T/Tref_, n0_) + EquationOfState::Cv(p, T);
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::Es
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+        c0_*(pow(T, n0_ + 1) - pow(Tstd, n0_ + 1))/(pow(Tref_, n0_)*(n0_ + 1))
+      + EquationOfState::E(p, T);
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::Ea
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return Es(p, T) + Hf();
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::Hf() const
+{
+    return Hf_;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::S
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+        c0_*(pow(T, n0_) - pow(Tstd, n0_))/(pow(Tref_, n0_)*n0_)
+      + EquationOfState::S(p, T); // Requires EquationOfState::S for Cv
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::Gstd
+(
+    const scalar T
+) const
+{
+    return
+        c0_*(pow(T, n0_ + 1) - pow(Tstd, n0_ + 1))/(pow(Tref_, n0_)*(n0_ + 1))
+      + Hf() + Pstd/EquationOfState::rho(Pstd, T)
+      - S(Pstd, T)*T;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::ePowerThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0; // + EquationOfState::dCpdT
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+inline void Foam::ePowerThermo<EquationOfState>::operator+=
+(
+    const ePowerThermo<EquationOfState>& ct
+)
+{
+    scalar Y1 = this->Y();
+
+    EquationOfState::operator+=(ct);
+
+    if (mag(this->Y()) > small)
+    {
+        Y1 /= this->Y();
+        const scalar Y2 = ct.Y()/this->Y();
+
+        Hf_ = Y1*Hf_ + Y2*ct.Hf_;
+        c0_ = Y1*c0_ + Y2*ct.c0_;
+        n0_ = Y1*n0_ + Y2*ct.n0_;
+        Tref_ = Y1*Tref_ + Y2*ct.Tref_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+template<class EquationOfState>
+inline Foam::ePowerThermo<EquationOfState> Foam::operator+
+(
+    const ePowerThermo<EquationOfState>& ct1,
+    const ePowerThermo<EquationOfState>& ct2
+)
+{
+    EquationOfState eofs
+    (
+        static_cast<const EquationOfState&>(ct1)
+      + static_cast<const EquationOfState&>(ct2)
+    );
+
+    if (mag(eofs.Y()) < small)
+    {
+        return ePowerThermo<EquationOfState>
+        (
+            eofs,
+            ct1.c0_,
+            ct1.n0_,
+            ct1.Tref_,
+            ct1.Hf_
+        );
+    }
+    else
+    {
+        return ePowerThermo<EquationOfState>
+        (
+            eofs,
+            ct1.Y()/eofs.Y()*ct1.c0_
+          + ct2.Y()/eofs.Y()*ct2.c0_,
+            ct1.Y()/eofs.Y()*ct1.n0_
+          + ct2.Y()/eofs.Y()*ct2.n0_,
+            ct1.Y()/eofs.Y()*ct1.Tref_
+          + ct2.Y()/eofs.Y()*ct2.Tref_,
+            ct1.Y()/eofs.Y()*ct1.Hf_
+          + ct2.Y()/eofs.Y()*ct2.Hf_
+        );
+    }
+}
+
+
+template<class EquationOfState>
+inline Foam::ePowerThermo<EquationOfState> Foam::operator*
+(
+    const scalar s,
+    const ePowerThermo<EquationOfState>& ct
+)
+{
+    return ePowerThermo<EquationOfState>
+    (
+        s*static_cast<const EquationOfState&>(ct),
+        ct.c0_,
+        ct.n0_,
+        ct.Tref_,
+        ct.Hf_
+    );
+}
+
+
+template<class EquationOfState>
+inline Foam::ePowerThermo<EquationOfState> Foam::operator==
+(
+    const ePowerThermo<EquationOfState>& ct1,
+    const ePowerThermo<EquationOfState>& ct2
+)
+{
+    EquationOfState eofs
+    (
+        static_cast<const EquationOfState&>(ct1)
+     == static_cast<const EquationOfState&>(ct2)
+    );
+
+    return ePowerThermo<EquationOfState>
+    (
+        eofs,
+        ct2.Y()/eofs.Y()*ct2.c0_
+      - ct1.Y()/eofs.Y()*ct1.c0_,
+        ct2.Y()/eofs.Y()*ct2.n0_
+      - ct1.Y()/eofs.Y()*ct1.n0_,
+        ct2.Y()/eofs.Y()*ct2.Tref_
+      - ct1.Y()/eofs.Y()*ct1.Tref_,
+        ct2.Y()/eofs.Y()*ct2.Hf_
+      - ct1.Y()/eofs.Y()*ct1.Hf_
+    );
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/templates/materials/aluminium/thermophysicalProperties

@@ -19,10 +19,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -41,7 +41,7 @@ mixture
     thermodynamics
     {
         Hf      0;
-        Cp      896;    // [J/kg/K]
+        Cv      896;    // [J/kg/K]
     }
 
     equationOfState

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/templates/materials/copper/thermophysicalProperties

@@ -19,10 +19,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -41,7 +41,7 @@ mixture
     thermodynamics
     {
         Hf      0;
-        Cp      385;        // [J/kg/K]
+        Cv      385;        // [J/kg/K]
     }
 
     equationOfState

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/templates/system/solid/fvOptions

@@ -37,7 +37,7 @@ fixedPower
 
     sources
     {
-        h
+        e
         {
             explicit    $power;
             implicit    0;

tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/constant/porous/fvOptions

@@ -24,7 +24,7 @@ porousToair
     master          true;
 
     nbrModel        airToporous;
-    fields          (h);
+    fields          (e);
     semiImplicit    no;
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/constant/heater/fvOptions

@@ -28,7 +28,7 @@ options
 
         sources
         {
-            h
+            e
             {
                 explicit    1e7; // W/m^3 == kg/m/s^3
                 implicit    0;

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/constant/heater/thermophysicalProperties

@@ -20,10 +20,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -45,7 +45,7 @@ mixture
     thermodynamics
     {
         Hf              0;
-        Cp              450;
+        Cv              450;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/constant/metal/thermophysicalProperties

@@ -20,10 +20,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -45,7 +45,7 @@ mixture
     thermodynamics
     {
         Hf              0;
-        Cp              900;
+        Cv              900;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/system/heater/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(alpha,h) Gauss linear corrected;
+    laplacian(alpha,e) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/system/heater/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    h
+    e
     {
         solver           GAMG;
         smoother         symGaussSeidel;
@@ -25,9 +25,9 @@ solvers
         relTol           0.1;
     }
 
-    hFinal
+    eFinal
     {
-        $h;
+        $e;
         relTol           0;
     }
 }

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/system/metal/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(alpha,h) Gauss linear corrected;
+    laplacian(alpha,e) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/heatTransfer/chtMultiRegionFoam/heatedDuct/system/metal/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    h
+    e
     {
         solver           GAMG;
         smoother         symGaussSeidel;
@@ -25,9 +25,9 @@ solvers
         relTol           0.1;
     }
 
-    hFinal
+    eFinal
     {
-        $h;
+        $e;
         relTol           0;
     }
 }

tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/solid/thermophysicalProperties

@@ -19,10 +19,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -38,7 +38,7 @@ mixture
     thermodynamics
     {
         Hf              0;
-        Cp              450;
+        Cv              450;
     }
     transport
     {

tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/system/solid/fvSchemes

@@ -32,7 +32,7 @@ divSchemes
 laplacianSchemes
 {
     default             none;
-    laplacian(alpha,h)  Gauss linear corrected;
+    laplacian(alpha,e)  Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/system/solid/fvSolution

@@ -16,7 +16,7 @@ FoamFile
 
 solvers
 {
-    h
+    e
     {
         solver           PCG;
         preconditioner   DIC;
@@ -24,9 +24,9 @@ solvers
         relTol           0.1;
     }
 
-    hFinal
+    eFinal
     {
-        $h;
+        $e;
         tolerance        1e-06;
         relTol           0;
     }
@@ -41,7 +41,7 @@ relaxationFactors
 {
     equations
     {
-        h               0.7;
+        e               0.7;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/shellAndTubeHeatExchanger/constant/solid/thermophysicalProperties

@@ -20,10 +20,10 @@ thermoType
     type            heSolidThermo;
     mixture         pureMixture;
     transport       constIso;
-    thermo          hConst;
+    thermo          eConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
@@ -45,7 +45,7 @@ mixture
     thermodynamics
     {
         Hf              0;
-        Cp              900;
+        Cv              900;
     }
 }
 

tutorials/heatTransfer/chtMultiRegionFoam/shellAndTubeHeatExchanger/system/solid/fvSchemes

@@ -33,7 +33,7 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(alpha,h) Gauss linear corrected;
+    laplacian(alpha,e) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/heatTransfer/chtMultiRegionFoam/shellAndTubeHeatExchanger/system/solid/fvSolution

@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    "h.*"
+    "e.*"
     {
         solver           GAMG;
         smoother         symGaussSeidel;