ThermophysicalTransportModels: Added temperature gradient based heat flux models

Fourier, eddyDiffusivity and nonUnityLewisEddyDiffusivity thermophysical
transport models now apply an implicit energy correction to a temperature
gradient based heat-flux to provide computational stability and efficiency while
converging to temperature gradient based solution.  This ensures consistent heat
exchange between fluid and solid regions in CHT cases and with heat-flux
boundaries.

The Fourier and eddyDiffusivity models support single specie systems only
whereas nonUnityLewisEddyDiffusivity supports specie diffusion with independent
specification of turbulent Prandtl and Schmidt numbers, i.e. non-unity Lewis
number.

The unityLewisFourier and unityLewisEddyDiffusivity thermophysical transport
models use an implicit energy gradient based heat-flux which is optimal for
numerical stability and convergence but does not guarantee consistent heat
exchange between fluid and solid regions and heat-flux boundaries in the
presence of gradients of heat capacity.  Both of these models support specie
diffusion with the restriction that the laminar and turbulent Prandtl and
Schmidt numbers are equal, i.e. unity Lewis number.

The thermophysical transport model is specified in the optional
thermophysicalTransport dictionary; if this file is not present the
unityLewisFourier model is selected for laminar and unityLewisEddyDiffusivity
for turbulent cases for backward compatibility.

The chtMultiRegionFoam tutorial cases have been updated to use the most
appropriate of the new thermophysical transport models.

src/ThermophysicalTransportModels/fluidReactionThermo/fluidReactionThermophysicalTransportModels.C

@@ -42,6 +42,9 @@ makeThermophysicalTransportModels
 #include "Fourier.H"
 makeLaminarThermophysicalTransportModel(Fourier);
 
+#include "unityLewisFourier.H"
+makeLaminarThermophysicalTransportModel(unityLewisFourier);
+
 
 // -------------------------------------------------------------------------- //
 // RAS models
@@ -50,6 +53,9 @@ makeLaminarThermophysicalTransportModel(Fourier);
 #include "eddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(RAS, eddyDiffusivity);
 
+#include "unityLewisEddyDiffusivity.H"
+makeRASLESThermophysicalTransportModel(RAS, unityLewisEddyDiffusivity);
+
 #include "nonUnityLewisEddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(RAS, nonUnityLewisEddyDiffusivity);
 
@@ -61,6 +67,9 @@ makeRASLESThermophysicalTransportModel(RAS, nonUnityLewisEddyDiffusivity);
 #include "eddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(LES, eddyDiffusivity);
 
+#include "unityLewisEddyDiffusivity.H"
+makeRASLESThermophysicalTransportModel(LES, unityLewisEddyDiffusivity);
+
 #include "nonUnityLewisEddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(LES, nonUnityLewisEddyDiffusivity);
 

src/ThermophysicalTransportModels/fluidThermo/fluidThermophysicalTransportModels.C

@@ -42,6 +42,9 @@ makeThermophysicalTransportModels
 #include "Fourier.H"
 makeLaminarThermophysicalTransportModel(Fourier);
 
+#include "unityLewisFourier.H"
+makeLaminarThermophysicalTransportModel(unityLewisFourier);
+
 
 // -------------------------------------------------------------------------- //
 // RAS models
@@ -50,6 +53,9 @@ makeLaminarThermophysicalTransportModel(Fourier);
 #include "eddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(RAS, eddyDiffusivity);
 
+#include "unityLewisEddyDiffusivity.H"
+makeRASLESThermophysicalTransportModel(RAS, unityLewisEddyDiffusivity);
+
 
 // -------------------------------------------------------------------------- //
 // LES models
@@ -58,5 +64,8 @@ makeRASLESThermophysicalTransportModel(RAS, eddyDiffusivity);
 #include "eddyDiffusivity.H"
 makeRASLESThermophysicalTransportModel(LES, eddyDiffusivity);
 
+#include "unityLewisEddyDiffusivity.H"
+makeRASLESThermophysicalTransportModel(LES, unityLewisEddyDiffusivity);
+
 
 // ************************************************************************* //

src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C

@@ -25,6 +25,7 @@ License
 
 #include "Fourier.H"
 #include "fvmLaplacian.H"
+#include "fvcLaplacian.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -68,9 +69,46 @@ bool Fourier<BasicThermophysicalTransportModel>::read()
 }
 
 
+template<class TurbulenceThermophysicalTransportModel>
+tmp<volScalarField>
+Fourier<TurbulenceThermophysicalTransportModel>::DEff
+(
+    const volScalarField& Yi
+) const
+{
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " unityLewisFourier"
+        << exit(FatalError);
+
+    return tmp<volScalarField>(nullptr);
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<scalarField>
+Fourier<TurbulenceThermophysicalTransportModel>::DEff
+(
+    const volScalarField& Yi,
+    const label patchi
+) const
+{
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " unityLewisFourier"
+        << exit(FatalError);
+
+    return tmp<scalarField>(nullptr);
+}
+
+
 template<class BasicThermophysicalTransportModel>
 tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::q() const
 {
+    const thermoModel& thermo = this->thermo();
+
     return volVectorField::New
     (
         IOobject::groupName
@@ -78,7 +116,7 @@ tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::q() const
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->thermo().alpha()*this->alpha()*fvc::grad(this->thermo().he())
+       -(this->alpha()*thermo.kappa())*fvc::grad(thermo.T())
     );
 }
 
@@ -87,7 +125,13 @@ template<class BasicThermophysicalTransportModel>
 tmp<fvScalarMatrix>
 Fourier<BasicThermophysicalTransportModel>::divq(volScalarField& he) const
 {
-    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), he);
+    const thermoModel& thermo = this->thermo();
+
+    // Return heat flux source as an implicit energy correction
+    // to the temperature gradient flux
+    return
+        -correction(fvm::laplacian(this->alpha()*thermo.alpha(), he))
+        -fvc::laplacian(this->alpha()*thermo.kappa(), thermo.T());
 }
 
 
@@ -97,15 +141,13 @@ tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::j
     const volScalarField& Yi
 ) const
 {
-    return volVectorField::New
+    FatalErrorInFunction
-    (
+        << type() << " supports single component systems only, " << nl
-        IOobject::groupName
+        << "    for multi-component transport select"
-        (
+           " unityLewisFourier"
-            "j(" + Yi.name() + ')',
+        << exit(FatalError);
-            this->momentumTransport().alphaRhoPhi().group()
+
-        ),
+    return tmp<volVectorField>(nullptr);
-       -this->thermo().alpha()*this->alpha()*fvc::grad(Yi)
-    );
 }
 
 
@@ -113,7 +155,13 @@ template<class BasicThermophysicalTransportModel>
 tmp<fvScalarMatrix>
 Fourier<BasicThermophysicalTransportModel>::divj(volScalarField& Yi) const
 {
-    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), Yi);
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " unityLewisFourier"
+        << exit(FatalError);
+
+    return tmp<fvScalarMatrix>(nullptr);
 }
 
 

src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H

@@ -25,7 +25,12 @@ Class
     Foam::laminarThermophysicalTransportModels::Fourier
 
 Description
-    Fourier's gradient heat flux model for laminar flow.
+    Fourier's temperature gradient heat flux model
+    for single specie laminar flow.
+
+    The heat flux source is implemented as an implicit energy correction to the
+    temperature gradient based flux source.  At convergence the energy
+    correction is 0.
 
 SourceFiles
     Fourier.C
@@ -96,6 +101,17 @@ public:
         //- Read thermophysicalTransport dictionary
         virtual bool read();
 
+        //- Effective mass diffusivity for a given specie mass-fraction [kg/m/s]
+        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const;
+
+        //- Effective mass diffusivity for a given specie mass-fraction
+        //  for patch [kg/m/s]
+        virtual tmp<scalarField> DEff
+        (
+            const volScalarField& Yi,
+            const label patchi
+        ) const;
+
         //- Return the heat flux
         virtual tmp<volVectorField> q() const;
 

src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "laminarThermophysicalTransportModel.H"
-#include "Fourier.H"
+#include "unityLewisFourier.H"
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -126,12 +126,12 @@ Foam::laminarThermophysicalTransportModel
     else
     {
         Info<< "Selecting default laminar thermophysical transport model "
-            << laminarThermophysicalTransportModels::Fourier<
+            << laminarThermophysicalTransportModels::unityLewisFourier<
                BasicThermophysicalTransportModel>::typeName << endl;
 
         return autoPtr<laminarThermophysicalTransportModel>
         (
-            new laminarThermophysicalTransportModels::Fourier
+            new laminarThermophysicalTransportModels::unityLewisFourier
             <
                 BasicThermophysicalTransportModel
             >(momentumTransport, thermo)

src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.H

@@ -182,14 +182,7 @@ public:
         }
 
         //- Effective mass diffusivity for a given specie mass-fraction [kg/m/s]
-        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const
+        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const = 0;
-        {
-            return volScalarField::New
-            (
-                "DEff",
-                alphaEff()
-            );
-        }
 
         //- Effective mass diffusivity for a given specie mass-fraction
         //  for patch [kg/m/s]
@@ -197,10 +190,7 @@ public:
         (
             const volScalarField& Yi,
             const label patchi
-        ) const
+        ) const = 0;
-        {
-            return alphaEff(patchi);
-        }
 
         //- Correct the laminar transport
         virtual void correct();

src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.C

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "unityLewisFourier.H"
+#include "fvmLaplacian.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarThermophysicalTransportModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicThermophysicalTransportModel>
+unityLewisFourier<BasicThermophysicalTransportModel>::unityLewisFourier
+(
+    const momentumTransportModel& momentumTransport,
+    const thermoModel& thermo
+)
+:
+    laminarThermophysicalTransportModel<BasicThermophysicalTransportModel>
+    (
+        typeName,
+        momentumTransport,
+        thermo
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicThermophysicalTransportModel>
+const dictionary&
+unityLewisFourier<BasicThermophysicalTransportModel>::coeffDict() const
+{
+    return dictionary::null;
+}
+
+
+template<class BasicThermophysicalTransportModel>
+bool unityLewisFourier<BasicThermophysicalTransportModel>::read()
+{
+    return true;
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<volVectorField>
+unityLewisFourier<BasicThermophysicalTransportModel>::q() const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "q",
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+       -this->thermo().alpha()*this->alpha()*fvc::grad(this->thermo().he())
+    );
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+unityLewisFourier<BasicThermophysicalTransportModel>::
+divq(volScalarField& he) const
+{
+    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), he);
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<volVectorField>unityLewisFourier<BasicThermophysicalTransportModel>::j
+(
+    const volScalarField& Yi
+) const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "j(" + Yi.name() + ')',
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+       -this->thermo().alpha()*this->alpha()*fvc::grad(Yi)
+    );
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+unityLewisFourier<BasicThermophysicalTransportModel>::
+divj(volScalarField& Yi) const
+{
+    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), Yi);
+}
+
+
+template<class BasicThermophysicalTransportModel>
+void unityLewisFourier<BasicThermophysicalTransportModel>::correct()
+{
+    laminarThermophysicalTransportModel
+    <
+        BasicThermophysicalTransportModel
+    >::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarThermophysicalTransportModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::laminarThermophysicalTransportModels::unityLewisFourier
+
+Description
+    unityLewisFourier's energy gradient heat flux model for laminar flow.
+    Specie fluxes are computed assuming a unity turbulent Lewis number.
+
+SourceFiles
+    unityLewisFourier.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef unityLewisFourier_H
+#define unityLewisFourier_H
+
+#include "laminarThermophysicalTransportModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarThermophysicalTransportModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class unityLewisFourier Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicThermophysicalTransportModel>
+class unityLewisFourier
+:
+    public laminarThermophysicalTransportModel
+    <
+        BasicThermophysicalTransportModel
+    >
+{
+
+public:
+
+    typedef typename BasicThermophysicalTransportModel::alphaField
+        alphaField;
+
+    typedef typename BasicThermophysicalTransportModel::momentumTransportModel
+        momentumTransportModel;
+
+    typedef typename BasicThermophysicalTransportModel::thermoModel
+        thermoModel;
+
+
+    //- Runtime type information
+    TypeName("unityLewisFourier");
+
+
+    // Constructors
+
+        //- Construct from components
+        unityLewisFourier
+        (
+            const momentumTransportModel& momentumTransport,
+            const thermoModel& thermo
+        );
+
+
+    //- Destructor
+    virtual ~unityLewisFourier()
+    {}
+
+
+    // Member Functions
+
+        //- Const access to the coefficients dictionary
+        virtual const dictionary& coeffDict() const;
+
+        //- Read thermophysicalTransport dictionary
+        virtual bool read();
+
+        //- Effective mass diffusivity for a given specie mass-fraction [kg/m/s]
+        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const
+        {
+            return volScalarField::New
+            (
+                "DEff",
+                this->alphaEff()
+            );
+        }
+
+        //- Effective mass diffusivity for a given specie mass-fraction
+        //  for patch [kg/m/s]
+        virtual tmp<scalarField> DEff
+        (
+            const volScalarField& Yi,
+            const label patchi
+        ) const
+        {
+            return this->alphaEff(patchi);
+        }
+
+        //- Return the heat flux
+        virtual tmp<volVectorField> q() const;
+
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
+
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
+        //- Correct the unityLewisFourier viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarThermophysicalTransportModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "unityLewisFourier.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/turbulence/LES/LESThermophysicalTransportModel/LESThermophysicalTransportModel.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LESThermophysicalTransportModel.H"
-#include "eddyDiffusivity.H"
+#include "unityLewisEddyDiffusivity.H"
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -126,22 +126,22 @@ Foam::LESThermophysicalTransportModel
     else
     {
         typedef
-            turbulenceThermophysicalTransportModels::eddyDiffusivity
+            turbulenceThermophysicalTransportModels::unityLewisEddyDiffusivity
             <
                 LESThermophysicalTransportModel
                 <
                     BasicThermophysicalTransportModel
                 >
-            > LESeddyDiffusivity;
+            > LESunityLewisEddyDiffusivity;
 
         Info<< "Selecting default LES thermophysical transport model "
-            <<  LESeddyDiffusivity::typeName << endl;
+            <<  LESunityLewisEddyDiffusivity::typeName << endl;
 
         return autoPtr<LESThermophysicalTransportModel>
         (
-            new LESeddyDiffusivity
+            new LESunityLewisEddyDiffusivity
             (
-                LESeddyDiffusivity::typeName,
+                LESunityLewisEddyDiffusivity::typeName,
                 momentumTransport,
                 thermo,
                 true

src/ThermophysicalTransportModels/turbulence/RAS/RASThermophysicalTransportModel/RASThermophysicalTransportModel.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "RASThermophysicalTransportModel.H"
-#include "eddyDiffusivity.H"
+#include "unityLewisEddyDiffusivity.H"
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -126,22 +126,22 @@ Foam::RASThermophysicalTransportModel
     else
     {
         typedef
-            turbulenceThermophysicalTransportModels::eddyDiffusivity
+            turbulenceThermophysicalTransportModels::unityLewisEddyDiffusivity
             <
                 RASThermophysicalTransportModel
                 <
                     BasicThermophysicalTransportModel
                 >
-            > RASeddyDiffusivity;
+            > RASunityLewisEddyDiffusivity;
 
         Info<< "Selecting default RAS thermophysical transport model "
-            <<  RASeddyDiffusivity::typeName << endl;
+            <<  RASunityLewisEddyDiffusivity::typeName << endl;
 
         return autoPtr<RASThermophysicalTransportModel>
         (
-            new RASeddyDiffusivity
+            new RASunityLewisEddyDiffusivity
             (
-                RASeddyDiffusivity::typeName,
+                RASunityLewisEddyDiffusivity::typeName,
                 momentumTransport,
                 thermo,
                 true

src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.C

@@ -25,6 +25,7 @@ License
 
 #include "eddyDiffusivity.H"
 #include "fvmLaplacian.H"
+#include "fvcLaplacian.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,25 +54,6 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::eddyDiffusivity
     const momentumTransportModel& momentumTransport,
     const thermoModel& thermo
 )
-:
-    eddyDiffusivity
-    (
-        typeName,
-        momentumTransport,
-        thermo,
-        false
-    )
-{}
-
-
-template<class TurbulenceThermophysicalTransportModel>
-eddyDiffusivity<TurbulenceThermophysicalTransportModel>::eddyDiffusivity
-(
-    const word& type,
-    const momentumTransportModel& momentumTransport,
-    const thermoModel& thermo,
-    const bool allowDefaultPrt
-)
 :
     TurbulenceThermophysicalTransportModel
     (
@@ -80,22 +62,7 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::eddyDiffusivity
         thermo
     ),
 
-    Prt_
+    Prt_(dimensioned<scalar>("Prt", dimless, this->coeffDict_)),
-    (
-        allowDefaultPrt
-      ? dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Prt",
-            this->coeffDict_,
-            1
-        )
-      : dimensioned<scalar>
-        (
-            "Prt",
-            dimless,
-            this->coeffDict_
-        )
-    ),
 
     alphat_
     (
@@ -123,7 +90,7 @@ bool eddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
 {
     if (TurbulenceThermophysicalTransportModel::read())
     {
-        Prt_.readIfPresent(this->coeffDict());
+        Prt_.read(this->coeffDict());
 
         return true;
     }
@@ -134,6 +101,41 @@ bool eddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
 }
 
 
+template<class TurbulenceThermophysicalTransportModel>
+tmp<volScalarField>
+eddyDiffusivity<TurbulenceThermophysicalTransportModel>::DEff
+(
+    const volScalarField& Yi
+) const
+{
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " nonUnityLewisEddyDiffusivity or unityLewisEddyDiffusivity"
+        << exit(FatalError);
+
+    return tmp<volScalarField>(nullptr);
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<scalarField>
+eddyDiffusivity<TurbulenceThermophysicalTransportModel>::DEff
+(
+    const volScalarField& Yi,
+    const label patchi
+) const
+{
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " nonUnityLewisEddyDiffusivity or unityLewisEddyDiffusivity"
+        << exit(FatalError);
+
+    return tmp<scalarField>(nullptr);
+}
+
+
 template<class TurbulenceThermophysicalTransportModel>
 tmp<volVectorField>
 eddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
@@ -145,7 +147,7 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->alphaEff()*this->alpha()*fvc::grad(this->thermo().he())
+       -(this->alpha()*this->kappaEff()*fvc::grad(this->thermo().T()))
     );
 }
 
@@ -157,7 +159,11 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     volScalarField& he
 ) const
 {
-    return -fvm::laplacian(this->alpha()*this->alphaEff(), he);
+    // Return heat flux source as an implicit energy correction
+    // to the temperature gradient flux
+    return
+        -correction(fvm::laplacian(this->alpha()*this->alphaEff(), he))
+        -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T());
 }
 
 
@@ -168,15 +174,13 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::j
     const volScalarField& Yi
 ) const
 {
-    return volVectorField::New
+    FatalErrorInFunction
-    (
+        << type() << " supports single component systems only, " << nl
-        IOobject::groupName
+        << "    for multi-component transport select"
-        (
+           " nonUnityLewisEddyDiffusivity or unityLewisEddyDiffusivity"
-            "j(" + Yi.name() + ')',
+        << exit(FatalError);
-            this->momentumTransport().alphaRhoPhi().group()
+
-        ),
+    return tmp<volVectorField>(nullptr);
-       -this->DEff(Yi)*this->alpha()*fvc::grad(Yi)
-    );
 }
 
 
@@ -187,7 +191,13 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divj
     volScalarField& Yi
 ) const
 {
-    return -fvm::laplacian(this->alpha()*this->DEff(Yi), Yi);
+    FatalErrorInFunction
+        << type() << " supports single component systems only, " << nl
+        << "    for multi-component transport select"
+           " nonUnityLewisEddyDiffusivity or unityLewisEddyDiffusivity"
+        << exit(FatalError);
+
+    return tmp<fvScalarMatrix>(nullptr);
 }
 
 

src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.H

@@ -25,9 +25,12 @@ Class
     Foam::turbulenceThermophysicalTransportModels::eddyDiffusivity
 
 Description
-    Eddy-diffusivity based gradient heat flux model for RAS or LES
+    Eddy-diffusivity based temperature gradient heat flux model
-    of turbulent flow. Specie fluxes are computed assuming a unity turbulent
+    for single specie RAS or LES of turbulent flow.
-    Lewis number.
+
+    The heat flux source is implemented as an implicit energy correction to the
+    temperature gradient based flux source.  At convergence the energy
+    correction is 0.
 
 Usage
     \verbatim
@@ -109,17 +112,6 @@ public:
             const thermoModel& thermo
         );
 
-        //- Construct from a type name, a momentum transport model and a thermo
-        //  model, and whether to default the turbulent Prandtl number to one
-        //  if it is not specified
-        eddyDiffusivity
-        (
-            const word& type,
-            const momentumTransportModel& momentumTransport,
-            const thermoModel& thermo,
-            const bool allowDefaultPrt
-        );
-
 
     //- Destructor
     virtual ~eddyDiffusivity()
@@ -179,14 +171,7 @@ public:
         }
 
         //- Effective mass diffusivity for a given specie mass-fraction [kg/m/s]
-        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const
+        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const;
-        {
-            return volScalarField::New
-            (
-                "DEff",
-                alphaEff()
-            );
-        }
 
         //- Effective mass diffusivity for a given specie mass-fraction
         //  for patch [kg/m/s]
@@ -194,10 +179,7 @@ public:
         (
             const volScalarField& Yi,
             const label patchi
-        ) const
+        ) const;
-        {
-            return alphaEff(patchi);
-        }
 
         //- Return the heat flux
         virtual tmp<volVectorField> q() const;

src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C

@@ -24,7 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "nonUnityLewisEddyDiffusivity.H"
+#include "fvcDiv.H"
 #include "fvcLaplacian.H"
+#include "fvcSnGrad.H"
+#include "fvmSup.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,7 +47,7 @@ nonUnityLewisEddyDiffusivity
     const thermoModel& thermo
 )
 :
-    eddyDiffusivity<TurbulenceThermophysicalTransportModel>
+    unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>
     (
         typeName,
         momentumTransport,
@@ -69,7 +73,11 @@ template<class TurbulenceThermophysicalTransportModel>
 bool
 nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
 {
-    if (eddyDiffusivity<TurbulenceThermophysicalTransportModel>::read())
+    if
+    (
+        unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>
+        ::read()
+    )
     {
         Sct_.readIfPresent(this->coeffDict());
 
@@ -86,26 +94,28 @@ template<class TurbulenceThermophysicalTransportModel>
 tmp<volVectorField>
 nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
 {
-    tmp<volVectorField> tmpq =
+    tmp<volVectorField> tmpq
-        eddyDiffusivity<TurbulenceThermophysicalTransportModel>::q();
+    (
-
+        volVectorField::New
-    if (mag(this->Prt_ - Sct_).value() > small)
-    {
-        const basicSpecieMixture& composition = this->thermo().composition();
-
-        const PtrList<volScalarField>& Y = composition.Y();
-
-        volScalarField alphatMinusDt
         (
-            "alphatMinusDt",
+            IOobject::groupName
-            this->alphat()*(1 - this->Prt_/Sct_)
+            (
-        );
+                "q",
+                this->momentumTransport().alphaRhoPhi().group()
+            ),
+           -(this->alpha()*this->kappaEff()*fvc::grad(this->thermo().T()))
+        )
+    );
 
+    const basicSpecieMixture& composition = this->thermo().composition();
+    const PtrList<volScalarField>& Y = composition.Y();
+
+    if (Y.size())
+    {
         forAll(Y, i)
         {
-            tmpq.ref() +=
+            tmpq.ref() -=
-                this->alpha()
+                this->alpha()*DEff(Y[i])
-               *alphatMinusDt
                *composition.HE(i, this->thermo().p(), this->thermo().T())
                *fvc::grad(Y[i]);
         }
@@ -122,32 +132,70 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     volScalarField& he
 ) const
 {
-    tmp<fvScalarMatrix> tmpDivq =
+    tmp<fvScalarMatrix> tmpDivq
-        eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq(he);
+    (
-
+        fvm::Su
-    if (mag(this->Prt_ - Sct_).value() > small)
-    {
-        const basicSpecieMixture& composition = this->thermo().composition();
-
-        const PtrList<volScalarField>& Y = composition.Y();
-
-        volScalarField alphatMinusDt
         (
-            "alphatMinusDt",
+            -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T()),
-            this->alphat()*(1 - this->Prt_/Sct_)
+            he
+        )
+    );
+
+    const basicSpecieMixture& composition = this->thermo().composition();
+    const PtrList<volScalarField>& Y = composition.Y();
+
+    if (!Y.size())
+    {
+        tmpDivq.ref() -=
+            correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+    }
+    else
+    {
+        tmpDivq.ref() -= fvm::laplacian(this->alpha()*this->alphaEff(), he);
+
+        volScalarField heNew
+        (
+            volScalarField::New
+            (
+                "he",
+                he.mesh(),
+                dimensionedScalar(he.dimensions(), 0)
+            )
+        );
+
+        surfaceScalarField hGradY
+        (
+            surfaceScalarField::New
+            (
+                "hGradY",
+                he.mesh(),
+                dimensionedScalar(he.dimensions()/dimLength, 0)
+            )
         );
 
         forAll(Y, i)
         {
-            tmpDivq.ref() +=
+            const volScalarField hi
-                fvc::laplacian
+            (
+                composition.HE(i, this->thermo().p(), this->thermo().T())
+            );
+
+            heNew += Y[i]*hi;
+            hGradY += fvc::interpolate(hi)*fvc::snGrad(Y[i]);
+        }
+
+        tmpDivq.ref() +=
+            fvc::laplacian(this->alpha()*this->alphaEff(), heNew);
+
+        tmpDivq.ref() -=
+            fvc::div
+            (
+                fvc::interpolate
                 (
                     this->alpha()
-                   *alphatMinusDt
+                   *this->thermo().alphaEff((this->Prt_/Sct_)*this->alphat())
-                   *composition.HE(i, this->thermo().p(), this->thermo().T()),
+                )*hGradY*he.mesh().magSf()
-                    Y[i]
+            );
-                );
-        }
     }
 
     return tmpDivq;

src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.H

@@ -25,10 +25,14 @@ Class
     Foam::turbulenceThermophysicalTransportModels::nonUnityLewisEddyDiffusivity
 
 Description
-    Non-unity-Lewis-Eddy-diffusivity based gradient heat flux model for RAS or
+    Non-unity-Lewis-Eddy-diffusivity based temperature gradient heat flux model
-    LES of turbulent flow. Allows independent specification of turbulent
+    for RAS or LES of turbulent flow. Allows independent specification of
-    Prandtl and Schmidt numbers. The laminar Lewis number is still assumed to
+    turbulent Prandtl and Schmidt numbers.  Unity laminar Lewis number is
-    equal one.
+    assumed.
+
+    The heat flux source is implemented as an implicit energy correction to the
+    temperature gradient based flux source.  At convergence the energy
+    correction is 0.
 
 Usage
     \verbatim
@@ -45,7 +49,7 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#include "eddyDiffusivity.H"
+#include "unityLewisEddyDiffusivity.H"
 
 #ifndef nonUnityLewisEddyDiffusivity_H
 #define nonUnityLewisEddyDiffusivity_H
@@ -64,7 +68,7 @@ namespace turbulenceThermophysicalTransportModels
 template<class TurbulenceThermophysicalTransportModel>
 class nonUnityLewisEddyDiffusivity
 :
-    public eddyDiffusivity<TurbulenceThermophysicalTransportModel>
+    public unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>
 {
 
 protected:
@@ -120,7 +124,7 @@ public:
             return volScalarField::New
             (
                 "DEff",
-                this->thermo().alphaEff(this->Prt_/Sct_*this->alphat())
+                this->thermo().alphaEff((this->Prt_/Sct_)*this->alphat())
             );
         }
 

src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.C

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "unityLewisEddyDiffusivity.H"
+#include "fvmLaplacian.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace turbulenceThermophysicalTransportModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class TurbulenceThermophysicalTransportModel>
+void unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::
+correctAlphat()
+{
+    alphat_ =
+        this->momentumTransport().rho()
+       *this->momentumTransport().nut()/Prt_;
+    alphat_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class TurbulenceThermophysicalTransportModel>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::
+unityLewisEddyDiffusivity
+(
+    const momentumTransportModel& momentumTransport,
+    const thermoModel& thermo
+)
+:
+    unityLewisEddyDiffusivity
+    (
+        typeName,
+        momentumTransport,
+        thermo,
+        false
+    )
+{}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::
+unityLewisEddyDiffusivity
+(
+    const word& type,
+    const momentumTransportModel& momentumTransport,
+    const thermoModel& thermo,
+    const bool allowDefaultPrt
+)
+:
+    TurbulenceThermophysicalTransportModel
+    (
+        type,
+        momentumTransport,
+        thermo
+    ),
+
+    Prt_
+    (
+        allowDefaultPrt
+      ? dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Prt",
+            this->coeffDict_,
+            1
+        )
+      : dimensioned<scalar>
+        (
+            "Prt",
+            dimless,
+            this->coeffDict_
+        )
+    ),
+
+    alphat_
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                "alphat",
+                this->momentumTransport().alphaRhoPhi().group()
+            ),
+            momentumTransport.time().timeName(),
+            momentumTransport.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        momentumTransport.mesh()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class TurbulenceThermophysicalTransportModel>
+bool unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
+{
+    if (TurbulenceThermophysicalTransportModel::read())
+    {
+        Prt_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<volVectorField>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "q",
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+       -this->alphaEff()*this->alpha()*fvc::grad(this->thermo().he())
+    );
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
+(
+    volScalarField& he
+) const
+{
+    return -fvm::laplacian(this->alpha()*this->alphaEff(), he);
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<volVectorField>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::j
+(
+    const volScalarField& Yi
+) const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "j(" + Yi.name() + ')',
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+       -this->DEff(Yi)*this->alpha()*fvc::grad(Yi)
+    );
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divj
+(
+    volScalarField& Yi
+) const
+{
+    return -fvm::laplacian(this->alpha()*this->DEff(Yi), Yi);
+}
+
+
+template<class TurbulenceThermophysicalTransportModel>
+void unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::
+correct()
+{
+    TurbulenceThermophysicalTransportModel::correct();
+    correctAlphat();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace turbulenceThermophysicalTransportModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.H

@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::turbulenceThermophysicalTransportModels::unityLewisEddyDiffusivity
+
+Description
+    Eddy-diffusivity based energy gradient heat flux model for RAS or LES
+    of turbulent flow.  Specie fluxes are computed assuming a unity turbulent
+    Lewis number.
+
+Usage
+    \verbatim
+    LES
+    {
+        model           unityLewisEddyDiffusivity;
+        Prt             0.85;
+    }
+    \endverbatim
+
+SourceFiles
+    unityLewisEddyDiffusivity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef unityLewisEddyDiffusivity_H
+#define unityLewisEddyDiffusivity_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace turbulenceThermophysicalTransportModels
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class unityLewisEddyDiffusivity Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class TurbulenceThermophysicalTransportModel>
+class unityLewisEddyDiffusivity
+:
+    public TurbulenceThermophysicalTransportModel
+{
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            //- Turbulent Prandtl number []
+            dimensionedScalar Prt_;
+
+        // Fields
+
+            //- Turbulent thermal diffusivity of enthalpy [kg/m/s]
+            volScalarField alphat_;
+
+
+    // Protected Member Functions
+
+        virtual void correctAlphat();
+
+
+public:
+
+    typedef typename TurbulenceThermophysicalTransportModel::alphaField
+        alphaField;
+
+    typedef typename
+        TurbulenceThermophysicalTransportModel::momentumTransportModel
+        momentumTransportModel;
+
+    typedef typename TurbulenceThermophysicalTransportModel::thermoModel
+        thermoModel;
+
+
+    //- Runtime type information
+    TypeName("unityLewisEddyDiffusivity");
+
+
+    // Constructors
+
+        //- Construct from a momentum transport model and a thermo model
+        unityLewisEddyDiffusivity
+        (
+            const momentumTransportModel& momentumTransport,
+            const thermoModel& thermo
+        );
+
+        //- Construct from a type name, a momentum transport model and a thermo
+        //  model, and whether to default the turbulent Prandtl number to one
+        //  if it is not specified
+        unityLewisEddyDiffusivity
+        (
+            const word& type,
+            const momentumTransportModel& momentumTransport,
+            const thermoModel& thermo,
+            const bool allowDefaultPrt
+        );
+
+
+    //- Destructor
+    virtual ~unityLewisEddyDiffusivity()
+    {}
+
+
+    // Member Functions
+
+        //- Read thermophysicalTransport dictionary
+        virtual bool read();
+
+        //- Turbulent thermal diffusivity for enthalpy [kg/m/s]
+        virtual tmp<volScalarField> alphat() const
+        {
+            return alphat_;
+        }
+
+        //- Turbulent thermal diffusivity for enthalpy for a patch [kg/m/s]
+        virtual tmp<scalarField> alphat(const label patchi) const
+        {
+            return alphat()().boundaryField()[patchi];
+        }
+
+        //- Effective thermal turbulent diffusivity for temperature
+        //  of mixture [W/m/K]
+        virtual tmp<volScalarField> kappaEff() const
+        {
+            return this->thermo().kappaEff(alphat());
+        }
+
+        //- Effective thermal turbulent diffusivity for temperature
+        //  of mixture for patch [W/m/K]
+        virtual tmp<scalarField> kappaEff(const label patchi) const
+        {
+            return this->thermo().kappaEff
+            (
+                alphat(patchi),
+                patchi
+            );
+        }
+
+        //- Effective thermal turbulent diffusivity of mixture [kg/m/s]
+        virtual tmp<volScalarField> alphaEff() const
+        {
+            return this->thermo().alphaEff(alphat());
+        }
+
+        //- Effective thermal turbulent diffusivity of mixture
+        //  for patch [kg/m/s]
+        virtual tmp<scalarField> alphaEff(const label patchi) const
+        {
+            return this->thermo().alphaEff
+            (
+                alphat(patchi),
+                patchi
+            );
+        }
+
+        //- Effective mass diffusivity for a given specie mass-fraction [kg/m/s]
+        virtual tmp<volScalarField> DEff(const volScalarField& Yi) const
+        {
+            return volScalarField::New
+            (
+                "DEff",
+                alphaEff()
+            );
+        }
+
+        //- Effective mass diffusivity for a given specie mass-fraction
+        //  for patch [kg/m/s]
+        virtual tmp<scalarField> DEff
+        (
+            const volScalarField& Yi,
+            const label patchi
+        ) const
+        {
+            return alphaEff(patchi);
+        }
+
+        //- Return the heat flux
+        virtual tmp<volVectorField> q() const;
+
+        //- Return the source term for the energy equation
+        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
+
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
+        //- Correct the unityLewisEddyDiffusivity viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace turbulenceThermophysicalTransportModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "unityLewisEddyDiffusivity.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/constant/thermophysicalTransport

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RAS
+{
+    model           eddyDiffusivity;
+
+    Prt             0.85;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/templates/materials/air/thermophysicalTransport

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RAS
+{
+    model           eddyDiffusivity;
+
+    Prt             0.85;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/coolingSphere/templates/materials/water/thermophysicalTransport

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RAS
+{
+    model           eddyDiffusivity;
+
+    Prt             0.85;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/heatExchanger/constant/porous/thermophysicalTransport

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+laminar
+{
+    model           Fourier;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/reverseBurner/constant/gas/thermophysicalTransport

@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RAS
+{
+    model           nonUnityLewisEddyDiffusivity;
+
+    Prt             0.85;
+    Sct             0.7;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/shellAndTubeHeatExchanger/constant/shell/thermophysicalTransport

@@ -0,0 +1,26 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalTransport;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RAS
+{
+    model           eddyDiffusivity;
+
+    Prt             0.85;
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/chtMultiRegionFoam/shellAndTubeHeatExchanger/constant/tube/thermophysicalTransport

@@ -0,0 +1 @@
+../shell/thermophysicalTransport