Improved turbulence modelling and wall functions. Vuko Vukcevic

src/finiteVolume/fields/fvPatchFields/derived/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C

@@ -48,6 +48,7 @@ turbulentIntensityKineticEnergyInletFvPatchScalarField
     this->valueFraction() = 0.0;
 }
 
+
 Foam::turbulentIntensityKineticEnergyInletFvPatchScalarField::
 turbulentIntensityKineticEnergyInletFvPatchScalarField
 (
@@ -101,6 +102,7 @@ turbulentIntensityKineticEnergyInletFvPatchScalarField
     this->valueFraction() = 0.0;
 }
 
+
 Foam::turbulentIntensityKineticEnergyInletFvPatchScalarField::
 turbulentIntensityKineticEnergyInletFvPatchScalarField
 (

src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C

@@ -192,10 +192,12 @@ tmp<volSymmTensorField> LamBremhorstKE::devReff() const
 
 tmp<fvVectorMatrix> LamBremhorstKE::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C

@@ -198,10 +198,12 @@ tmp<volSymmTensorField> LaunderSharmaKE::devReff() const
 
 tmp<fvVectorMatrix> LaunderSharmaKE::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C

@@ -276,11 +276,13 @@ tmp<volSymmTensorField> LienCubicKE::devReff() const
 
 tmp<fvVectorMatrix> LienCubicKE::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
         fvc::div(nonlinearStress_)
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/LienCubicKELowRe/LienCubicKELowRe.C

@@ -346,11 +346,13 @@ tmp<volSymmTensorField> LienCubicKELowRe::devReff() const
 
 tmp<fvVectorMatrix> LienCubicKELowRe::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
         fvc::div(nonlinearStress_)
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C

@@ -234,11 +234,12 @@ tmp<volSymmTensorField> LienLeschzinerLowRe::devReff() const
 
 tmp<fvVectorMatrix> LienLeschzinerLowRe::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-    //- (fvc::grad(U_) & fvc::grad(nuEff()))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
-      - fvc::div(nuEff()*T(fvc::grad(U_)))
     );
 }
 

src/turbulenceModels/incompressible/RAS/Make/files

@@ -32,24 +32,35 @@ $(nutWallFunctions)/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctio
 $(nutWallFunctions)/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C
+$(nutWallFunctions)/nutCWTWallFunction/nutCWTWallFunctionFvPatchScalarField.C
+$(nutWallFunctions)/nutMEWTWallFunction/nutMEWTWallFunctionFvPatchScalarField.C
 
 epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
 $(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+$(epsilonWallFunctions)/epsilonNonEqWallFunction/epsilonNonEqWallFunctionFvPatchScalarField.C
 
 omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
 $(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
+$(omegaWallFunctions)/omegaCWTWallFunction/omegaCWTWallFunctionFvPatchScalarField.C
+$(omegaWallFunctions)/omegaMEWTWallFunction/omegaMEWTWallFunctionFvPatchScalarField.C
 
 kqRWallFunctions = $(wallFunctions)/kqRWallFunctions
 $(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C
+$(kqRWallFunctions)/kNonEqWallFunction/kNonEqWallFunctionFvPatchScalarField.C
 
 RWallFunctions = $(wallFunctions)/RWallFunctions
 $(RWallFunctions)/RWallFunction/RWallFunctionFvPatchSymmTensorField.C
 
+wallFunctionFunctionObjects = $(wallFunctions)/functionObjects
+$(wallFunctionFunctionObjects)/pressureGradient/pressureGradientFunctionObject.C
+$(wallFunctionFunctionObjects)/velocityConvection/velocityConvectionFunctionObject.C
 
 /* Patch fields */
 derivedFvPatchFields/turbulentHeatFluxTemperature/turbulentHeatFluxTemperatureFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthDissipationRateInlet/turbulentMixingLengthDissipationRateInletFvPatchScalarField.C
 derivedFvPatchFields/turbulentMixingLengthFrequencyInlet/turbulentMixingLengthFrequencyInletFvPatchScalarField.C
+derivedFvPatchFields/nuRatioEpsilon/nuRatioEpsilonFvPatchScalarField.C
+derivedFvPatchFields/nuRatioOmega/nuRatioOmegaFvPatchScalarField.C
 
 backwardsCompatibility/wallFunctions/backwardsCompatibilityWallFunctions.C
 

src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C

@@ -264,11 +264,13 @@ tmp<volSymmTensorField> NonlinearKEShih::devReff() const
 
 tmp<fvVectorMatrix> NonlinearKEShih::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
         fvc::div(nonlinearStress_)
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -208,10 +208,12 @@ tmp<volSymmTensorField> RNGkEpsilon::devReff() const
 
 tmp<fvVectorMatrix> RNGkEpsilon::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -312,12 +312,12 @@ tmp<volSymmTensorField> SpalartAllmaras::devReff() const
 
 tmp<fvVectorMatrix> SpalartAllmaras::divDevReff() const
 {
-    volScalarField nuEff_ = nuEff();
+    const volScalarField nuEffective = nuEff();
 
     return
     (
-      - fvm::laplacian(nuEff_, U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff_*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/blockKOmegaSST/blockKOmegaSST.C

@@ -355,10 +355,12 @@ tmp<volSymmTensorField> blockKOmegaSST::devReff() const
 
 tmp<fvVectorMatrix> blockKOmegaSST::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C

@@ -194,10 +194,12 @@ tmp<volSymmTensorField> coupledKEpsilon::devReff() const
 
 tmp<fvVectorMatrix> coupledKEpsilon::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.C

@@ -355,10 +355,12 @@ tmp<volSymmTensorField> coupledKOmegaSST::devReff() const
 
 tmp<fvVectorMatrix> coupledKOmegaSST::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/nuRatioEpsilon/nuRatioEpsilonFvPatchScalarField.C

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nuRatioEpsilonFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    nuRatio_(1),
+    Cmu_(0.09),
+    kName_("undefined-k"),
+    phiName_("undefined-phi")
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField
+(
+    const nuRatioEpsilonFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    inletOutletFvPatchScalarField(ptf, p, iF, mapper),
+    nuRatio_(ptf.nuRatio_),
+    Cmu_(ptf.Cmu_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+Foam::nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    nuRatio_(readScalar(dict.lookup("nuRatio"))),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kName_(dict.lookupOrDefault<word>("kName", "k")),
+    phiName_(dict.lookupOrDefault<word>("phiName", "phi"))
+{
+    if (nuRatio_< SMALL)
+    {
+        FatalErrorIn
+        (
+            "nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField"
+            "(const fvPatch& p, const DimensionedField<scalar, volMesh>& iF, "
+            "const dictionary& dict)"
+        )   << "Invalid eddy viscosity ratio (nuRatio) specified: " << nuRatio_
+            << "\n    on patch " << this->patch().name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    if (Cmu_ < SMALL)
+    {
+        FatalErrorIn
+        (
+            "nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField"
+            "(const fvPatch& p, const DimensionedField<scalar, volMesh>& iF, "
+            "const dictionary& dict)"
+        )   << "Invalid Cmu_ constant specified: " << Cmu_
+            << "\n    on patch " << this->patch().name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField
+(
+    const nuRatioEpsilonFvPatchScalarField& ptf
+)
+:
+    inletOutletFvPatchScalarField(ptf),
+    nuRatio_(ptf.nuRatio_),
+    Cmu_(ptf.Cmu_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+Foam::nuRatioEpsilonFvPatchScalarField::nuRatioEpsilonFvPatchScalarField
+(
+    const nuRatioEpsilonFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(ptf, iF),
+    nuRatio_(ptf.nuRatio_),
+    Cmu_(ptf.Cmu_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::nuRatioEpsilonFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Get turbulent kinetic energy field for this patch
+    const fvPatchScalarField& kp =
+        lookupPatchField<volScalarField, vector>(kName_);
+
+    // Get flux field for this patch
+    const fvsPatchScalarField& phip =
+        lookupPatchField<surfaceScalarField, scalar>(phiName_);
+
+    // Get RASModel
+    const incompressible::RASModel& rasModel =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            incompressible::RASModel
+        >("RASProperties");
+
+    // Get laminar viscosity for this patch
+    const fvPatchScalarField& nup =
+        rasModel.nu().boundaryField()[this->patch().index()];
+
+    this->refValue() = Cmu_*sqrt(kp)/nuRatio_/nup;
+    this->valueFraction() = neg(phip);
+
+    inletOutletFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::nuRatioEpsilonFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("nuRatio") << nuRatio_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kName") << kName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("phiName") << phiName_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        nuRatioEpsilonFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/nuRatioEpsilon/nuRatioEpsilonFvPatchScalarField.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::nuRatioEpsilonFvPatchScalarField
+
+Description
+    Calculate turbulent dissipation rate (epsilon) from the provided eddy
+    viscosity ratio. Note: derived from inletOutlet so it can readily be used
+    for both inlets and outlets.
+
+    Example of the boundary condition specification:
+    \verbatim
+        inlet
+        {
+            type        nuRatioEpsilon;
+            nuRatio     10;             // Eddy viscosity ratio
+            Cmu         0.09;           // Default value 0.09
+            kName       k;              // Name of the turbulent energy field
+                                        // (k by default)
+            phiName     phi;            // Name of the flux field
+                                        // (phi by default)
+            value       uniform 1;      // placeholder
+        }
+    \endverbatim
+
+SourceFiles
+    nuRatioEpsilonFvPatchScalarField.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nuRatioEpsilonFvPatchScalarField_H
+#define nuRatioEpsilonFvPatchScalarField_H
+
+#include "inletOutletFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+     Class eddyViscosityRatioTurbulenceDissipationRateFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class nuRatioEpsilonFvPatchScalarField
+:
+    public inletOutletFvPatchScalarField
+{
+    // Private data
+
+        //- Target viscosity ratio
+        scalar nuRatio_;
+
+        //- Cmu constant
+        scalar Cmu_;
+
+        //- Name of the turbulent energy field
+        word kName_;
+
+        //- Name of the flux field
+        word phiName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nuRatioEpsilon");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        nuRatioEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        nuRatioEpsilonFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  nuRatioEpsilonFvPatchScalarField
+        //  onto a new patch
+        nuRatioEpsilonFvPatchScalarField
+        (
+            const nuRatioEpsilonFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        nuRatioEpsilonFvPatchScalarField
+        (
+            const nuRatioEpsilonFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nuRatioEpsilonFvPatchScalarField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        nuRatioEpsilonFvPatchScalarField
+        (
+            const nuRatioEpsilonFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nuRatioEpsilonFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/nuRatioOmega/nuRatioOmegaFvPatchScalarField.C

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nuRatioOmegaFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    nuRatio_(1),
+    kName_("undefined-k"),
+    phiName_("undefined-phi")
+{
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField
+(
+    const nuRatioOmegaFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    inletOutletFvPatchScalarField(ptf, p, iF, mapper),
+    nuRatio_(ptf.nuRatio_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+Foam::nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    inletOutletFvPatchScalarField(p, iF),
+    nuRatio_(readScalar(dict.lookup("nuRatio"))),
+    kName_(dict.lookupOrDefault<word>("kName", "k")),
+    phiName_(dict.lookupOrDefault<word>("phiName", "phi"))
+{
+    if (nuRatio_< SMALL)
+    {
+        FatalErrorIn
+        (
+            "nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField"
+            "(const fvPatch& p, const DimensionedField<scalar, volMesh>& iF, "
+            "const dictionary& dict)"
+        )   << "Invalid eddy viscosity ratio (nuRatio) specified: " << nuRatio_
+            << "\n    on patch " << this->patch().name()
+            << " of field " << this->dimensionedInternalField().name()
+            << " in file " << this->dimensionedInternalField().objectPath()
+            << exit(FatalError);
+    }
+
+    fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
+
+    this->refValue() = 0.0;
+    this->refGrad() = 0.0;
+    this->valueFraction() = 0.0;
+}
+
+
+Foam::nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField
+(
+    const nuRatioOmegaFvPatchScalarField& ptf
+)
+:
+    inletOutletFvPatchScalarField(ptf),
+    nuRatio_(ptf.nuRatio_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+Foam::nuRatioOmegaFvPatchScalarField::nuRatioOmegaFvPatchScalarField
+(
+    const nuRatioOmegaFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    inletOutletFvPatchScalarField(ptf, iF),
+    nuRatio_(ptf.nuRatio_),
+    kName_(ptf.kName_),
+    phiName_(ptf.phiName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::nuRatioOmegaFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Get turbulent kinetic energy field for this patch
+    const fvPatchScalarField& kp =
+        lookupPatchField<volScalarField, vector>(kName_);
+
+    // Get flux field for this patch
+    const fvsPatchScalarField& phip =
+        lookupPatchField<surfaceScalarField, scalar>(phiName_);
+
+    // Get RASModel
+    const incompressible::RASModel& rasModel =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            incompressible::RASModel
+        >("RASProperties");
+
+    // Get laminar viscosity for this patch
+    const fvPatchScalarField& nup =
+        rasModel.nu().boundaryField()[this->patch().index()];
+
+    this->refValue() = kp/(nuRatio_*nup);
+    this->valueFraction() = neg(phip);
+
+    inletOutletFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::nuRatioOmegaFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("nuRatio") << nuRatio_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kName") << kName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("phiName") << phiName_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        nuRatioOmegaFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/nuRatioOmega/nuRatioOmegaFvPatchScalarField.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::nuRatioOmegaFvPatchScalarField
+
+Description
+    Calculate turbulent dissipation rate (epsilon) from the provided eddy
+    viscosity ratio. Note: derived from inletOutlet so it can readily be used
+    for both inlets and outlets.
+
+    Example of the boundary condition specification:
+    \verbatim
+        inlet
+        {
+            type        nuRatioOmega;
+            nuRatio     10;             // Eddy viscosity ratio
+            kName       k;              // Name of the turbulent energy field
+                                        // (k by default)
+            phiName     phi;            // Name of the flux field
+                                        // (phi by default)
+            value       uniform 1;      // placeholder
+        }
+    \endverbatim
+
+SourceFiles
+    nuRatioOmegaFvPatchScalarField.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nuRatioOmegaFvPatchScalarField_H
+#define nuRatioOmegaFvPatchScalarField_H
+
+#include "inletOutletFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+     Class eddyViscosityRatioTurbulenceDissipationRateFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class nuRatioOmegaFvPatchScalarField
+:
+    public inletOutletFvPatchScalarField
+{
+    // Private data
+
+        //- Target viscosity ratio
+        scalar nuRatio_;
+
+        //- Name of the turbulent energy field
+        word kName_;
+
+        //- Name of the flux field
+        word phiName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nuRatioOmega");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        nuRatioOmegaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        nuRatioOmegaFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  nuRatioOmegaFvPatchScalarField
+        //  onto a new patch
+        nuRatioOmegaFvPatchScalarField
+        (
+            const nuRatioOmegaFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        nuRatioOmegaFvPatchScalarField
+        (
+            const nuRatioOmegaFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nuRatioOmegaFvPatchScalarField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        nuRatioOmegaFvPatchScalarField
+        (
+            const nuRatioOmegaFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nuRatioOmegaFvPatchScalarField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+
+    // Member functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonNonEqWallFunction/epsilonNonEqWallFunctionFvPatchScalarField.C

@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "epsilonNonEqWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void epsilonNonEqWallFunctionFvPatchScalarField::checkType()
+{
+    if (!this->patch().isWall())
+    {
+        FatalErrorIn("epsilonNonEqWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+epsilonNonEqWallFunctionFvPatchScalarField::
+epsilonNonEqWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchScalarField(p, iF),
+    UName_("U"),
+    kName_("k"),
+    nuName_("nu"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41)
+{
+    checkType();
+}
+
+
+epsilonNonEqWallFunctionFvPatchScalarField::
+epsilonNonEqWallFunctionFvPatchScalarField
+(
+    const epsilonNonEqWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedInternalValueFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    kName_(ptf.kName_),
+    nuName_(ptf.nuName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_)
+{
+    checkType();
+}
+
+
+epsilonNonEqWallFunctionFvPatchScalarField::
+epsilonNonEqWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedInternalValueFvPatchScalarField(p, iF, dict),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    kName_(dict.lookupOrDefault<word>("k", "k")),
+    nuName_(dict.lookupOrDefault<word>("nu", "nu")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+{
+    checkType();
+}
+
+
+epsilonNonEqWallFunctionFvPatchScalarField::
+epsilonNonEqWallFunctionFvPatchScalarField
+(
+    const epsilonNonEqWallFunctionFvPatchScalarField& ewfpsf
+)
+:
+    fixedInternalValueFvPatchScalarField(ewfpsf),
+    UName_(ewfpsf.UName_),
+    kName_(ewfpsf.kName_),
+    nuName_(ewfpsf.nuName_),
+    nutName_(ewfpsf.nutName_),
+    Cmu_(ewfpsf.Cmu_),
+    kappa_(ewfpsf.kappa_)
+{
+    checkType();
+}
+
+
+epsilonNonEqWallFunctionFvPatchScalarField::
+epsilonNonEqWallFunctionFvPatchScalarField
+(
+    const epsilonNonEqWallFunctionFvPatchScalarField& ewfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchScalarField(ewfpsf, iF),
+    UName_(ewfpsf.UName_),
+    kName_(ewfpsf.kName_),
+    nuName_(ewfpsf.nuName_),
+    nutName_(ewfpsf.nutName_),
+    Cmu_(ewfpsf.Cmu_),
+    kappa_(ewfpsf.kappa_)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void epsilonNonEqWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patch().index()];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu75 = pow(Cmu_, 0.75);
+
+    // Epsilon is now a refValue (derived from fixedInternalValueFvPatchField)
+    scalarField& epsilon = refValue();
+
+    const volScalarField& k = db().lookupObject<volScalarField>(kName_);
+
+    const scalarField& nuw =
+        lookupPatchField<volScalarField, scalar>(nuName_);
+
+     const scalarField& nutw =
+        lookupPatchField<volScalarField, scalar>(nutName_);
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    // Set epsilon
+    forAll(nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+
+        // Viscous sublayer thickness
+        const scalar yVis = nuw[faceI]*11.225/(Cmu25*sqrt(k[faceCellI]));
+
+        if (y[faceI] > yVis)
+        {
+           epsilon[faceI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
+        }
+        else
+        {
+           epsilon[faceI] = 2.0*nuw[faceI]*k[faceCellI]/sqr(y[faceI]);
+        }
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+
+    fixedInternalValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void epsilonNonEqWallFunctionFvPatchScalarField::evaluate
+(
+    const Pstream::commsTypes commsType
+)
+{
+    fixedInternalValueFvPatchScalarField::evaluate(commsType);
+}
+
+
+void epsilonNonEqWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fixedInternalValueFvPatchScalarField::write(os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    epsilonNonEqWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonNonEqWallFunction/epsilonNonEqWallFunctionFvPatchScalarField.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::epsilonNonEqWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for epsilon when using wall functions
+    - calculates epsilon and G
+    - epsilon values added directly into the matrix to act as a constraint
+    Takes into acount non-equilibrium effects
+
+    Reference (bibtex entry):
+
+    @article{sutaloThesis2017,
+        author = {\v{S}utalo, F.}
+        title = {{Evaluation of Variants of Enhanced Wall Treatment Wall
+                  Function in Turbulent Flow Simulations}},
+        school = {Faculty of Mechanical Engineering and Naval Architecture,
+                  University of Zagreb},
+        year = {2017},
+    }
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+SourceFiles
+    epsilonNonEqWallFunctionFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef epsilonNonEqWallFunctionFvPatchScalarField_H
+#define epsilonNonEqWallFunctionFvPatchScalarField_H
+
+#include "fixedInternalValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class epsilonNonEqWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class epsilonNonEqWallFunctionFvPatchScalarField
+:
+    public fixedInternalValueFvPatchScalarField
+{
+    // Private data
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of turbulence kinetic energy field
+        word kName_;
+
+        //- Name of turbulence generation field
+        word GName_;
+
+        //- Name of laminar viscosity field
+        word nuName_;
+
+        //- Name of turbulent viscosity field
+        word nutName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("epsilonNonEqWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        epsilonNonEqWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        epsilonNonEqWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  epsilonNonEqWallFunctionFvPatchScalarField
+        //  onto a new patch
+        epsilonNonEqWallFunctionFvPatchScalarField
+        (
+            const epsilonNonEqWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        epsilonNonEqWallFunctionFvPatchScalarField
+        (
+            const epsilonNonEqWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new epsilonNonEqWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        epsilonNonEqWallFunctionFvPatchScalarField
+        (
+            const epsilonNonEqWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new epsilonNonEqWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+            //- Evaluate the patchField
+            virtual void evaluate
+            (
+                const Pstream::commsTypes  = Pstream::blocking
+            );
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/functionObjects/pressureGradient/pressureGradientFunctionObject.C

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "pressureGradientFunctionObject.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(pressureGradientFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        pressureGradientFunctionObject,
+        dictionary
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pressureGradientFunctionObject::pressureGradientFunctionObject
+(
+    const word& name,
+    const Time& t,
+    const dictionary& dict
+)
+:
+    functionObject(name),
+    name_(name),
+    time_(t),
+    regionName_(dict.lookupOrDefault<word>("region", polyMesh::defaultRegion)),
+    pName_(dict.lookup("pName")),
+    gradp_
+    (
+        IOobject
+        (
+            "pressureGradient",
+            time_.timeName(),
+            time_.lookupObject<fvMesh>(regionName_),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        time_.lookupObject<fvMesh>(regionName_),
+        // Note: dimensions reset in the constructor body
+        dimensionedVector("zero", dimless, vector::zero)
+    )
+{
+    Info<< "Creating pressureGradientFunctionObject for "
+        << pName_ << " field." << endl;
+
+    // Reset dimensions
+    const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName_);
+
+    if (mesh.foundObject<volScalarField>(pName_))
+    {
+        // Found pressure, reset dimensions for pressure gradient
+        const volScalarField& p =
+            mesh.lookupObject<volScalarField>(pName_);
+
+        gradp_.dimensions().reset(p.dimensions()/dimLength);
+    }
+    else
+    {
+        WarningIn
+        (
+            "pressureGradientFunctionObject::"
+            "pressureGradientFunctionObject"
+            "\n("
+            "\n    const word& name,"
+            "\n    const Time& t,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Could not find " << pName_ << " field."
+            << " Possible dimensions mismatch."
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::pressureGradientFunctionObject::start()
+{
+    return true;
+}
+
+
+bool Foam::pressureGradientFunctionObject::execute()
+{
+    const fvMesh& mesh =
+        time_.lookupObject<fvMesh>(regionName_);
+
+    if (mesh.foundObject<volScalarField>(pName_))
+    {
+        const volScalarField& p =
+            mesh.lookupObject<volScalarField>(pName_);
+
+        gradp_ = fvc::grad(p);
+
+        return true;
+    }
+    else
+    {
+        InfoIn("bool pressureGradientFunctionObject::execute()")
+            << "Field " << pName_ << " not found.  Returning."
+            << endl;
+
+        return false;
+    }
+}
+
+
+bool Foam::pressureGradientFunctionObject::read(const dictionary& dict)
+{
+    pName_ = word(dict.lookup("pName"));
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/functionObjects/pressureGradient/pressureGradientFunctionObject.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    pressureGradientFunctionObject
+
+Description
+    FunctionObject that calculates the pressure gradient given a name for the
+    pressure field and registers it into the database. Useful for pressure
+    sensitised wall functions (see omega/nutCWT/MEWT wall functions).
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+SourceFiles
+    pressureGradientFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pressureGradientFunctionObject_H
+#define pressureGradientFunctionObject_H
+
+#include "functionObject.H"
+#include "dictionary.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class pressureGradientFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class pressureGradientFunctionObject
+:
+    public functionObject
+{
+    // Private data
+
+        //- Name
+        const word name_;
+
+        //- Reference to main object registry
+        const Time& time_;
+
+        //- Region name
+        word regionName_;
+
+        //- Name of the pressure field
+        word pName_;
+
+        //- Pressure gradient field
+        volVectorField gradp_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        pressureGradientFunctionObject
+        (
+            const pressureGradientFunctionObject&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const pressureGradientFunctionObject&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pressureGradient");
+
+
+    // Constructors
+
+        //- Construct from components
+        pressureGradientFunctionObject
+        (
+            const word& name,
+            const Time&,
+            const dictionary&
+        );
+
+
+    // Member Functions
+
+        //- start is called at the start of the time-loop
+        virtual bool start();
+
+        //- execute is called at each ++ or += of the time-loop
+        virtual bool execute();
+
+        //- Read and set the function object if its data has changed
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/functionObjects/velocityConvection/velocityConvectionFunctionObject.C

@@ -0,0 +1,170 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+\*---------------------------------------------------------------------------*/
+
+#include "velocityConvectionFunctionObject.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvcDiv.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(velocityConvectionFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        velocityConvectionFunctionObject,
+        dictionary
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::velocityConvectionFunctionObject::velocityConvectionFunctionObject
+(
+    const word& name,
+    const Time& t,
+    const dictionary& dict
+)
+:
+    functionObject(name),
+    name_(name),
+    time_(t),
+    regionName_(dict.lookupOrDefault<word>("region", polyMesh::defaultRegion)),
+    UName_(dict.lookup("UName")),
+    phiName_(dict.lookup("phiName")),
+    convection_
+    (
+        IOobject
+        (
+            "velocityConvection",
+            time_.timeName(),
+            time_.lookupObject<fvMesh>(regionName_),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        time_.lookupObject<fvMesh>(regionName_),
+        // Note: dimensions will be overwritten on start
+        dimensionedVector("zero", dimless, vector::zero)
+    )
+{
+    Info<< "Creating velocityConvectionFunctionObject for "
+        << UName_ << " and " << phiName_ << " field." << endl;
+
+    // Reset dimensions
+    const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName_);
+
+    if
+    (
+        mesh.foundObject<volVectorField>(UName_)
+     && mesh.foundObject<surfaceScalarField>(phiName_)
+    )
+    {
+        // Found velocity and flux, reset dimensions for convection
+        const volVectorField& U =
+            mesh.lookupObject<volVectorField>(UName_);
+
+        const surfaceScalarField& phi =
+            mesh.lookupObject<surfaceScalarField>(phiName_);
+
+        convection_.dimensions().reset
+        (
+            phi.dimensions()*U.dimensions()/dimArea/dimLength
+        );
+    }
+    else
+    {
+        WarningIn
+        (
+            "velocityConvectionFunctionObject::"
+            "velocityConvectionFunctionObject"
+            "\n("
+            "\n    const word& name,"
+            "\n    const Time& t,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Could not find " << UName_ << " and " << phiName_
+            << " fields. Possible dimensions mismatch."
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::velocityConvectionFunctionObject::start()
+{
+    return true;
+}
+
+
+bool Foam::velocityConvectionFunctionObject::execute()
+{
+    const fvMesh& mesh = time_.lookupObject<fvMesh>(regionName_);
+
+    if
+    (
+        mesh.foundObject<volVectorField>(UName_)
+     && mesh.foundObject<surfaceScalarField>(phiName_)
+    )
+    {
+        // Found velocity and flux, reset dimensions for convection
+        const volVectorField& U =
+            mesh.lookupObject<volVectorField>(UName_);
+
+        const surfaceScalarField& phi =
+            mesh.lookupObject<surfaceScalarField>(phiName_);
+
+        convection_ = fvc::div(phi, U);
+
+        return true;
+    }
+    else
+    {
+        InfoIn("bool velocityConvectionFunctionObject::execute()")
+            << "Field " << UName_ << " or " << phiName_ << " not found."
+            << " Returning."
+            << endl;
+
+        return false;
+    }
+}
+
+
+bool Foam::velocityConvectionFunctionObject::read(const dictionary& dict)
+{
+    UName_ = word(dict.lookup("UName"));
+    phiName_ = word(dict.lookup("phiName"));
+
+    return false;
+}
+
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/functionObjects/velocityConvection/velocityConvectionFunctionObject.H

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    velocityConvectionFunctionObject
+
+Description
+    FunctionObject that calculates the velocity convection given names for the
+    velocity field and the flux field and registers it into the database. Useful
+    for convection sensitised wall functions (see omega/nutCWT/MEWT wall
+    functions).
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+SourceFiles
+    velocityConvectionFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef velocityConvectionFunctionObject_H
+#define velocityConvectionFunctionObject_H
+
+#include "functionObject.H"
+#include "dictionary.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class velocityConvectionFunctionObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class velocityConvectionFunctionObject
+:
+    public functionObject
+{
+    // Private data
+
+        //- Name
+        const word name_;
+
+        //- Reference to main object registry
+        const Time& time_;
+
+        //- Region name
+        word regionName_;
+
+        //- Name of the velocity field
+        word UName_;
+
+        //- Name of the flux field
+        word phiName_;
+
+        //- Convection field
+        volVectorField convection_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        velocityConvectionFunctionObject
+        (
+            const velocityConvectionFunctionObject&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const velocityConvectionFunctionObject&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("velocityConvection");
+
+
+    // Constructors
+
+        //- Construct from components
+        velocityConvectionFunctionObject
+        (
+            const word& name,
+            const Time&,
+            const dictionary&
+        );
+
+
+    // Member Functions
+
+        //- start is called at the start of the time-loop
+        virtual bool start();
+
+        //- execute is called at each ++ or += of the time-loop
+        virtual bool execute();
+
+        //- Read and set the function object if its data has changed
+        virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kNonEqWallFunction/kNonEqWallFunctionFvPatchScalarField.C

@@ -0,0 +1,292 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kNonEqWallFunctionFvPatchScalarField.H"
+#include "fvPatchFieldMapper.H"
+#include "addToRunTimeSelectionTable.H"
+#include "wallFvPatch.H"
+#include "RASModel.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+void kNonEqWallFunctionFvPatchScalarField::checkType()
+{
+    if (!this->patch().isWall())
+    {
+        FatalErrorIn("kNonEqWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << this->patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << this->patch().type()
+            << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+kNonEqWallFunctionFvPatchScalarField::kNonEqWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+
+)
+:
+    zeroGradientFvPatchScalarField(p, iF),
+    UName_("U"),
+    kName_("k"),
+    epsilonName_("epsilon"),
+    GName_("RASModel::G"),
+    nuName_("nu"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41)
+{
+    checkType();
+}
+
+
+kNonEqWallFunctionFvPatchScalarField::kNonEqWallFunctionFvPatchScalarField
+(
+    const kNonEqWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+
+)
+:
+    zeroGradientFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    kName_(ptf.kName_),
+    epsilonName_(ptf.epsilonName_),
+    GName_(ptf.GName_),
+    nuName_(ptf.nuName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_)
+
+{
+    checkType();
+}
+
+
+kNonEqWallFunctionFvPatchScalarField::kNonEqWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    zeroGradientFvPatchScalarField(p, iF, dict),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    kName_(dict.lookupOrDefault<word>("k", "k")),
+    epsilonName_(dict.lookupOrDefault<word>("epsilon", "epsilon")),
+    GName_(dict.lookupOrDefault<word>("G", "RASModel::G")),
+    nuName_(dict.lookupOrDefault<word>("nu", "nu")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41))
+
+{
+    checkType();
+}
+
+
+kNonEqWallFunctionFvPatchScalarField::kNonEqWallFunctionFvPatchScalarField
+(
+    const kNonEqWallFunctionFvPatchScalarField& tkqrwfpf
+)
+:
+    zeroGradientFvPatchScalarField(tkqrwfpf),
+    UName_(tkqrwfpf.UName_),
+    kName_(tkqrwfpf.kName_),
+    epsilonName_(tkqrwfpf.epsilonName_),
+    GName_(tkqrwfpf.GName_),
+    nuName_(tkqrwfpf.nuName_),
+    nutName_(tkqrwfpf.nutName_),
+    Cmu_(tkqrwfpf.Cmu_),
+    kappa_(tkqrwfpf.kappa_)
+{
+    checkType();
+}
+
+
+kNonEqWallFunctionFvPatchScalarField::kNonEqWallFunctionFvPatchScalarField
+(
+    const kNonEqWallFunctionFvPatchScalarField& tkqrwfpf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    zeroGradientFvPatchScalarField(tkqrwfpf, iF),
+    UName_(tkqrwfpf.UName_),
+    kName_(tkqrwfpf.kName_),
+    epsilonName_(tkqrwfpf.epsilonName_),
+    GName_(tkqrwfpf.GName_),
+    nuName_(tkqrwfpf.nuName_),
+    nutName_(tkqrwfpf.nutName_),
+    Cmu_(tkqrwfpf.Cmu_),
+    kappa_(tkqrwfpf.kappa_)
+{
+    checkType();
+}
+
+
+void kNonEqWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    // If G field is not present, execute zero gradient evaluation
+    // HJ, 20/Mar/2011
+    if ( !db().foundObject<volScalarField>(GName_) )
+
+    {
+        InfoIn("void kNonEqWallFunctionFvPatchScalarField::updateCoeffs()")
+            << "Cannot access " << GName_ << " field for patch "
+            << patch().name() << ".  Evaluating as zeroGradient"
+            << endl;
+
+        fvPatchScalarField::updateCoeffs();
+
+        zeroGradientFvPatchScalarField::evaluate();
+
+
+        return;
+    }
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patch().index()];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+    const scalar Cmu75 = pow(Cmu_, 0.75);
+
+    volScalarField& G = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>(GName_));
+    scalarField& GIn = G.internalField();
+
+    volScalarField& epsilon = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>(epsilonName_));
+    scalarField& epsilonIn = epsilon.internalField();
+
+    const volScalarField& k = db().lookupObject<volScalarField>(kName_);
+    const scalarField& kIn = k.internalField();
+
+    const scalarField& nuw =
+        lookupPatchField<volScalarField, scalar>(nuName_);
+
+    const scalarField& nutw =
+        lookupPatchField<volScalarField, scalar>(nutName_);
+
+    const fvPatchVectorField& Uw =
+        lookupPatchField<volVectorField, vector>(UName_);
+
+    const scalarField magGradUw = mag(Uw.snGrad());
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    // Averaged G and epsilon
+    forAll (nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+
+        // Viscous sublayer thickness
+        const scalar yVis = nuw[faceI]*11.225/(Cmu25*sqrt(k[faceCellI]));
+
+        // Height of the wall adjacent volume
+        const scalar yn = 2.0*y[faceI];  
+
+        if (y[faceI] > yVis)
+        {
+            GIn[faceCellI] =
+                sqr
+                (
+                    (nutw[faceI] + nuw[faceI])*magGradUw[faceI]
+                )
+               *log(yn/yVis)
+                /(Cmu25*sqrt(kIn[faceCellI])*kappa_*yn);
+
+            epsilonIn[faceCellI] =
+                (
+                    2.0*nuw[faceI]*kIn[faceCellI]/yVis+pow(kIn[faceCellI], 1.5)
+                   *Cmu75*log(yn/yVis)/kappa_
+                )/yn;
+        }
+        else
+        {
+            GIn[faceCellI] = 0.0;
+            epsilonIn[faceCellI] = 2.0*nuw[faceI]*kIn[faceCellI]/sqr(yVis);
+        }
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+
+    zeroGradientFvPatchScalarField::updateCoeffs();
+}
+
+void kNonEqWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    zeroGradientFvPatchScalarField::write(os);
+    this->writeEntry("value", os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "epsilon", "epsilon", epsilonName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    kNonEqWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kqRWallFunctions/kNonEqWallFunction/kNonEqWallFunctionFvPatchScalarField.H

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::kNonEqWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulence k, Q, and R when using wall functions.
+    Acts as a zero gradient condition and takes into account non equilibrium
+    effects.
+
+    Reference (bibtex entry):
+
+    @article{sutaloThesis2017,
+        author = {\v{S}utalo, F.}
+        title = {{Evaluation of Variants of Enhanced Wall Treatment Wall
+                  Function in Turbulent Flow Simulations}},
+        school = {Faculty of Mechanical Engineering and Naval Architecture,
+                  University of Zagreb},
+        year = {2017},
+    }
+
+SourceFiles
+    kNonEqWallFunctionFvPatchScalarField.C
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kNonEqWallFunctionFvPatchScalarField_H
+#define kNonEqWallFunctionFvPatchScalarField_H
+
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+          Class kNonEqWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class kNonEqWallFunctionFvPatchScalarField
+:
+    public zeroGradientFvPatchScalarField
+{
+    // Private data
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of turbulence kinetic energy field
+        word kName_;
+
+        //- Name of dissipation rate field
+        word epsilonName_;
+
+        //- Name of turbulence generation field
+        word GName_;
+
+        //- Name of laminar viscosity field
+        word nuName_;
+
+        //- Name of turbulent viscosity field
+        word nutName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kNonEqWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        kNonEqWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        kNonEqWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  kNonEqWallFunctionFvPatchScalarField
+        //  onto a new patch
+        kNonEqWallFunctionFvPatchScalarField
+        (
+            const kNonEqWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        kNonEqWallFunctionFvPatchScalarField
+        (
+            const kNonEqWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new kNonEqWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        kNonEqWallFunctionFvPatchScalarField
+        (
+            const kNonEqWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new kNonEqWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutCWTWallFunction/nutCWTWallFunctionFvPatchScalarField.C

@@ -0,0 +1,361 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nutCWTWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void nutCWTWallFunctionFvPatchScalarField::checkType()
+{
+    if (!patch().isWall())
+    {
+        FatalErrorIn("nutCWTWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+scalar nutCWTWallFunctionFvPatchScalarField::calcYPlusLam
+(
+    const scalar kappa,
+    const scalar E
+) const
+{
+    scalar ypl = 11.0;
+
+    for (int i = 0; i < 10; i++)
+    {
+        ypl = log(E*ypl)/kappa;
+    }
+
+    return ypl;
+}
+
+
+tmp<scalarField> nutCWTWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    // Get normals
+    const vectorField n = patch().nf();
+
+    // Patch velocity field at this wall
+    const fvPatchVectorField& Uw =
+          lookupPatchField<volVectorField, vector>(UName_);
+
+    const scalarField magGradUw = mag(Uw.snGrad());
+    const vectorField UwIn = Uw.patchInternalField();
+
+    // Patch internal velocity field tangential to the wall
+    const vectorField UwInTang = UwIn - (UwIn & n)*n;
+    const scalarField magUwInTang = mag(UwInTang);
+
+    // Calculate tangential direction for patch cells
+    const vectorField tDir = UwInTang/(magUwInTang + SMALL);
+
+    // Wall-velocity vector field tangential to the wall
+    const vectorField UwTang = Uw - (Uw & n)*n;
+    const scalarField magUwTang = mag(UwTang);
+
+
+    // Pressure effects. Lookup the pressure gradient field from the registry
+    // (created with pressureGradient function object)
+
+    // Pressure gradient projected on the wall parallel velocity direction
+    scalarField gradpTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("pressureGradient")
+    )
+    {
+        const volVectorField& gradP =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("pressureGradient");
+
+        // Update pressure gradient projected on the wall parallel velocity
+        gradpTang =
+            gradP.boundaryField()[this->patch().index()].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field pressureGradient not found. Neglecting pressure gradient "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Convective terms. Lookup the convection field from the registry (created
+    // with velocityConvection function object)
+
+    // Velocity convection projected on the wall parallel velocity direction
+    scalarField convectionTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("velocityConvection")
+    )
+    {
+        const volVectorField& convection =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("velocityConvection");
+
+        // Upate convection term projected on the wall parallel velocity
+        convectionTang =
+            convection.boundaryField()
+            [
+                this->patch().index()
+            ].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field velocityConvection not found. Neglecting convection "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), SMALL));
+    scalarField& nutw = tnutw();
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    forAll(nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+        const scalar uStar = Cmu25*sqrt(k[faceCellI]);
+
+	    // Note: here yPlus is actually yStar
+        const scalar yPlus = uStar*y[faceI]/nuw[faceI];
+
+        // Relative tangential velocity
+        const scalar magUrel = magUwInTang[faceI] - magUwTang[faceI];
+
+        const scalar Cu =  convectionTang[faceI] + gradpTang[faceI];
+        const scalar Psi = 1.0 - Cu/(kappa_*uStar*magGradUw[faceI]);
+
+        const scalar tauwVis = nuw[faceI]*magGradUw[faceI];
+        const scalar tauwLog = kappa_*uStar*magUrel*Psi/log(E_*yPlus);
+
+        // Kader blending
+        const scalar gamma = -0.01*pow(yPlus, 4)/(1.0 + 5.0*yPlus);
+        const scalar tauw = tauwVis*exp(gamma) + tauwLog*exp(1.0/gamma);
+
+        nutw[faceI] = tauw/(magGradUw[faceI] + SMALL) - nuw[faceI];
+    }
+
+    return tnutw;
+}
+
+
+void nutCWTWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "p", "p", pName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+nutCWTWallFunctionFvPatchScalarField::nutCWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    UName_("U"),
+    pName_("p"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
+{
+    checkType();
+}
+
+
+nutCWTWallFunctionFvPatchScalarField::nutCWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
+{
+    checkType();
+}
+
+
+nutCWTWallFunctionFvPatchScalarField::nutCWTWallFunctionFvPatchScalarField
+(
+    const nutCWTWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    pName_(ptf.pName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    E_(ptf.E_),
+    yPlusLam_(ptf.yPlusLam_)
+{
+    checkType();
+}
+
+
+nutCWTWallFunctionFvPatchScalarField::nutCWTWallFunctionFvPatchScalarField
+(
+    const nutCWTWallFunctionFvPatchScalarField& wfpsf
+)
+:
+    fixedValueFvPatchScalarField(wfpsf),
+    UName_(wfpsf.UName_),
+    pName_(wfpsf.pName_),
+    nutName_(wfpsf.nutName_),
+    Cmu_(wfpsf.Cmu_),
+    kappa_(wfpsf.kappa_),
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
+{
+    checkType();
+}
+
+
+nutCWTWallFunctionFvPatchScalarField::nutCWTWallFunctionFvPatchScalarField
+(
+    const nutCWTWallFunctionFvPatchScalarField& wfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(wfpsf, iF),
+    UName_(wfpsf.UName_),
+    pName_(wfpsf.pName_),
+    nutName_(wfpsf.nutName_),
+    Cmu_(wfpsf.Cmu_),
+    kappa_(wfpsf.kappa_),
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void nutCWTWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    operator==(calcNut());
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+tmp<scalarField> nutCWTWallFunctionFvPatchScalarField::yPlus() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField kwc = k.boundaryField()[patchI].patchInternalField();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    return pow(Cmu_, 0.25)*y*sqrt(kwc)/nuw;
+}
+
+
+void nutCWTWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    writeLocalEntries(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, nutCWTWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutCWTWallFunction/nutCWTWallFunctionFvPatchScalarField.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::nutCWTWallFunctionFvPatchScalarField
+
+Description
+    Improved boundary condition for turbulent (kinematic) viscosity, taking into
+    account non-equilibrium effects and are sensitive to flow unsteadiness and
+    presence of the pressure gradient.
+
+    Reference (bibtex entry):
+
+    @article{popovacHanjalic2007,
+        author = {Popovac, M. and Hanjali\'{c}, K.,}
+        title = {{Compound Wall Treatment for RANS Computation of Complex
+                  Turbulent Flows and Heat Transfer}},
+        journal = {Flow Turbulence Combust},
+        year = {2007},
+        pages = {78--177}
+        DOI = {10.1007/s10494-006-9067-x}
+    }
+
+SourceFiles
+    nutCWTWallFunctionFvPatchScalarField.C
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nutCWTWallFunctionFvPatchScalarField_H
+#define nutCWTWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+             Class nutCWTWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class nutCWTWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+protected:
+
+    // Protected data
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of pressure field
+        word pName_;
+
+        //- Name of eddy viscosity field
+        word nutName_;
+
+        //- Name of omega field
+        word omegaName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+        //- Y+ at the edge of the laminar sublayer
+        scalar yPlusLam_;
+
+
+    // Protected member functions
+
+        //- Check the type of the patch
+        virtual void checkType();
+
+        //- Calculate the Y+ at the edge of the laminar sublayer
+        virtual scalar calcYPlusLam(const scalar kappa, const scalar E) const;
+
+        //- Calculate the turbulence viscosity
+        virtual tmp<scalarField> calcNut() const;
+
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nutCWTWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        nutCWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        nutCWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  nutCWTWallFunctionFvPatchScalarField
+        //  onto a new patch
+        nutCWTWallFunctionFvPatchScalarField
+        (
+            const nutCWTWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        nutCWTWallFunctionFvPatchScalarField
+        (
+            const nutCWTWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nutCWTWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        nutCWTWallFunctionFvPatchScalarField
+        (
+            const nutCWTWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nutCWTWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Calculate and return the yPlus at the boundary
+            virtual tmp<scalarField> yPlus() const;
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutMEWTWallFunction/nutMEWTWallFunctionFvPatchScalarField.C

@@ -0,0 +1,409 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "nutMEWTWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::debug::tolerancesSwitch
+Foam::incompressible::RASModels::
+nutMEWTWallFunctionFvPatchScalarField::dimlessAFactorTol_
+(
+    "dimlessAFactorMEWTTolerance",
+    1e-6
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void nutMEWTWallFunctionFvPatchScalarField::checkType()
+{
+    if (!patch().isWall())
+    {
+        FatalErrorIn("nutMEWTWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+scalar nutMEWTWallFunctionFvPatchScalarField::calcYPlusLam
+(
+    const scalar kappa,
+    const scalar E
+) const
+{
+    scalar ypl = 11.0;
+
+    for (int i = 0; i < 10; i++)
+    {
+        ypl = log(E*ypl)/kappa;
+    }
+
+    return ypl;
+}
+
+
+tmp<scalarField> nutMEWTWallFunctionFvPatchScalarField::calcNut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    // Get normals
+    const vectorField n = patch().nf();
+
+    // Patch velocity field at this wall
+    const fvPatchVectorField& Uw =
+          lookupPatchField<volVectorField, vector>(UName_);
+
+    const scalarField magGradUw = mag(Uw.snGrad());
+    const vectorField UwIn = Uw.patchInternalField();
+
+    // Patch internal velocity field tangential to the wall
+    const vectorField UwInTang = UwIn - (UwIn & n)*n;
+    const scalarField magUwInTang = mag(UwInTang);
+
+    // Calculate tangential direction for patch cells
+    const vectorField tDir = UwInTang/(magUwInTang + SMALL);
+
+    // Wall-velocity vector field tangential to the wall
+    const vectorField UwTang = Uw - (Uw & n)*n;
+    const scalarField magUwTang = mag(UwTang);
+
+
+    // Pressure effects. Lookup the pressure gradient field from the registry
+    // (created with pressureGradient function object)
+
+    // Pressure gradient projected on the wall parallel velocity direction
+    scalarField gradpTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("pressureGradient")
+    )
+    {
+        const volVectorField& gradP =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("pressureGradient");
+
+        // Update pressure gradient projected on the wall parallel velocity
+        gradpTang =
+            gradP.boundaryField()[this->patch().index()].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field pressureGradient not found. Neglecting pressure gradient "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Convective terms. Lookup the convection field from the registry (created
+    // with velocityConvection function object)
+
+    // Velocity convection projected on the wall parallel velocity direction
+    scalarField convectionTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("velocityConvection")
+    )
+    {
+        const volVectorField& convection =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("velocityConvection");
+
+        // Upate convection term projected on the wall parallel velocity
+        convectionTang =
+            convection.boundaryField()
+            [
+                this->patch().index()
+            ].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field velocityConvection not found. Neglecting convection "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+    // Needed to calculate yPlus
+    const scalarField& eddyVis =
+        lookupPatchField<volScalarField, scalar>(nutName_);
+
+    tmp<scalarField> tnutw(new scalarField(patch().size(), SMALL));
+    scalarField& nutw = tnutw();
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    forAll(nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+        const scalar uStar = Cmu25*sqrt(k[faceCellI]);
+
+        // Calculate yPlus
+        const scalar yPlus =
+            sqrt
+            (
+                (eddyVis[faceI]+nuw[faceI])*magGradUw[faceI]
+            )*
+            y[faceI]/
+           (nuw[faceI] + SMALL);
+
+        // Relative tangential velocity
+        const scalar magUrel = magUwInTang[faceI] - magUwTang[faceI];
+
+        // Dimless A factor
+        scalar A = nuw[faceI]*
+            (gradpTang[faceI] + convectionTang[faceI])/(pow(uStar, 3) + SMALL);
+
+        // Numerical stabilisation of the A factor
+        if (A < SMALL)
+        {
+            A += dimlessAFactorTol_();
+        }
+
+        // Helper variables
+        const scalar S1 = sqrt(max(SMALL, 1.0 + A*yPlus));
+        const scalar p1 = sqrt(max(SMALL, 1.0 +6.0*A));
+        const scalar uPlusT =
+            (1.0/kappa_)*log(6.0*E_)
+           -(1.0/kappa_)*(2.0*p1 + log(mag(p1 - 1.0)) + log(p1 + 1.0));
+
+        // Dimless velocity in log layer
+        const scalar uLogPlus =
+            (1.0/kappa_)*(2.0*S1 + log(mag(S1 - 1.0)) + log(S1 + 1.0)) + uPlusT;
+
+        // Friction velocity in viscous sublayer
+        const scalar uTauVis = sqrt(nuw[faceI]*magUrel/(y[faceI] + SMALL));
+
+        // Friction velocity in log layer
+        const scalar uTauLog = magUrel/(uLogPlus + SMALL);
+
+        // Kader blending for friction velocity
+        const scalar gamma = -0.01*pow(yPlus, 4)/(1.0 + 5.0*yPlus);
+        const scalar uTau = uTauVis*exp(gamma) + uTauLog*exp(1.0/gamma);
+
+        // Need to limit nutw for stability reasons since at some point, yPlus
+        // is 0 and Kader blending is not defined
+        nutw[faceI] =
+            max
+            (
+                SMALL,
+                sqr(uTau)/(magGradUw[faceI] + SMALL) - nuw[faceI]
+            );
+    }
+
+    return tnutw;
+}
+
+
+void nutMEWTWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
+{
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "p", "p", pName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+nutMEWTWallFunctionFvPatchScalarField::nutMEWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    UName_("U"),
+    pName_("p"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
+{
+    checkType();
+}
+
+
+nutMEWTWallFunctionFvPatchScalarField::nutMEWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    yPlusLam_(calcYPlusLam(kappa_, E_))
+{
+    checkType();
+}
+
+
+nutMEWTWallFunctionFvPatchScalarField::nutMEWTWallFunctionFvPatchScalarField
+(
+    const nutMEWTWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    UName_(ptf.UName_),
+    pName_(ptf.pName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    E_(ptf.E_),
+    yPlusLam_(ptf.yPlusLam_)
+{
+    checkType();
+}
+
+
+nutMEWTWallFunctionFvPatchScalarField::nutMEWTWallFunctionFvPatchScalarField
+(
+    const nutMEWTWallFunctionFvPatchScalarField& wfpsf
+)
+:
+    fixedValueFvPatchScalarField(wfpsf),
+    UName_(wfpsf.UName_),
+    pName_(wfpsf.pName_),
+    nutName_(wfpsf.nutName_),
+    Cmu_(wfpsf.Cmu_),
+    kappa_(wfpsf.kappa_),
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
+{
+    checkType();
+}
+
+
+nutMEWTWallFunctionFvPatchScalarField::nutMEWTWallFunctionFvPatchScalarField
+(
+    const nutMEWTWallFunctionFvPatchScalarField& wfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(wfpsf, iF),
+    UName_(wfpsf.UName_),
+    pName_(wfpsf.pName_),
+    nutName_(wfpsf.nutName_),
+    Cmu_(wfpsf.Cmu_),
+    kappa_(wfpsf.kappa_),
+    E_(wfpsf.E_),
+    yPlusLam_(wfpsf.yPlusLam_)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void nutMEWTWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    operator==(calcNut());
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+tmp<scalarField> nutMEWTWallFunctionFvPatchScalarField::yPlus() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patchI];
+
+    const tmp<volScalarField> tk = rasModel.k();
+    const volScalarField& k = tk();
+    const scalarField kwc = k.boundaryField()[patchI].patchInternalField();
+    const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    return pow(Cmu_, 0.25)*y*sqrt(kwc)/nuw;
+}
+
+
+void nutMEWTWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    writeLocalEntries(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, nutMEWTWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutMEWTWallFunction/nutMEWTWallFunctionFvPatchScalarField.H

@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::nutMEWTWallFunctionFvPatchScalarField
+
+Description
+    Improved boundary condition for turbulent (kinematic) viscosity, taking into
+    account convection and pressure gradient terms.
+
+    Reference (bibtex entry):
+
+    @article{sutaloThesis2017,
+        author = {\v{S}utalo, F.}
+        title = {{Evaluation of Variants of Enhanced Wall Treatment Wall
+                  Function in Turbulent Flow Simulations}},
+        school = {Faculty of Mechanical Engineering and Naval Architecture,
+                  University of Zagreb},
+        year = {2017},
+    }
+
+SourceFiles
+    nutMEWTWallFunctionFvPatchScalarField.C
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef nutMEWTWallFunctionFvPatchScalarField_H
+#define nutMEWTWallFunctionFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+#include "tolerancesSwitch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+             Class nutMEWTWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class nutMEWTWallFunctionFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+
+    // Private static data member
+
+        //- Tolerance for the A dimensionless parameter
+        static const debug::tolerancesSwitch dimlessAFactorTol_;
+
+
+protected:
+
+    // Protected data
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of pressure field
+        word pName_;
+
+        //- Name of eddy viscosity field
+        word nutName_;
+
+        //- Name of omega field
+        word omegaName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+        //- Y+ at the edge of the laminar sublayer
+        scalar yPlusLam_;
+
+
+    // Protected member functions
+
+        //- Check the type of the patch
+        virtual void checkType();
+
+        //- Calculate the Y+ at the edge of the laminar sublayer
+        virtual scalar calcYPlusLam(const scalar kappa, const scalar E) const;
+
+        //- Calculate the turbulence viscosity
+        virtual tmp<scalarField> calcNut() const;
+
+        //- Write local wall function variables
+        virtual void writeLocalEntries(Ostream&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("nutMEWTWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        nutMEWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        nutMEWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  nutMEWTWallFunctionFvPatchScalarField
+        //  onto a new patch
+        nutMEWTWallFunctionFvPatchScalarField
+        (
+            const nutMEWTWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        nutMEWTWallFunctionFvPatchScalarField
+        (
+            const nutMEWTWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nutMEWTWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        nutMEWTWallFunctionFvPatchScalarField
+        (
+            const nutMEWTWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new nutMEWTWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Calculate and return the yPlus at the boundary
+            virtual tmp<scalarField> yPlus() const;
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaCWTWallFunction/omegaCWTWallFunctionFvPatchScalarField.C

@@ -0,0 +1,367 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "omegaCWTWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void omegaCWTWallFunctionFvPatchScalarField::checkType()
+{
+    if (!patch().isWall())
+    {
+        FatalErrorIn("omegaCWTWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF),
+    pName_("p"),
+    UName_("U"),
+    kName_("k"),
+    GName_("RASModel::G"),
+    nuName_("nu"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    beta1_(0.075)
+{
+    checkType();
+}
+
+
+omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
+(
+    const omegaCWTWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    pName_(ptf.pName_),
+    UName_(ptf.UName_),
+    kName_(ptf.kName_),
+    GName_(ptf.GName_),
+    nuName_(ptf.nuName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    E_(ptf.E_),
+    beta1_(ptf.beta1_)
+{
+    checkType();
+}
+
+
+omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF, dict),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    kName_(dict.lookupOrDefault<word>("k", "k")),
+    GName_(dict.lookupOrDefault<word>("G", "RASModel::G")),
+    nuName_(dict.lookupOrDefault<word>("nu", "nu")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075))
+{
+    checkType();
+}
+
+
+omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
+(
+    const omegaCWTWallFunctionFvPatchScalarField& owfpsf
+)
+:
+    fixedInternalValueFvPatchField<scalar>(owfpsf),
+    pName_(owfpsf.pName_),
+    UName_(owfpsf.UName_),
+    kName_(owfpsf.kName_),
+    GName_(owfpsf.GName_),
+    nuName_(owfpsf.nuName_),
+    nutName_(owfpsf.nutName_),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_)
+{
+    checkType();
+}
+
+
+omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
+(
+    const omegaCWTWallFunctionFvPatchScalarField& owfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
+    pName_(owfpsf.pName_),
+    UName_(owfpsf.UName_),
+    kName_(owfpsf.kName_),
+    GName_(owfpsf.GName_),
+    nuName_(owfpsf.nuName_),
+    nutName_(owfpsf.nutName_),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void omegaCWTWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // If G field is not present, execute zero gradient evaluation
+    // HJ, 20/Mar/2011
+    if (!db().foundObject<volScalarField>(GName_))
+    {
+        InfoIn("void omegaCWTWallFunctionFvPatchScalarField::updateCoeffs()")
+            << "Cannot access " << GName_ << " field for patch "
+            << patch().name() << ".  Evaluating as zeroGradient"
+            << endl;
+
+        fvPatchScalarField::updateCoeffs();
+        zeroGradientFvPatchScalarField::evaluate();
+
+        return;
+    }
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patch().index()];
+
+    const scalar Cmu25 = pow(Cmu_, 0.25);
+
+    volScalarField& G = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>(GName_));
+
+    // Note: omega is now a refValue and set in fixedInternalValueFvPatchField
+    // HJ, 3/Aug/2011
+    scalarField& omega = refValue();
+
+    const scalarField& k = db().lookupObject<volScalarField>(kName_);
+
+    const scalarField& nuw =
+        lookupPatchField<volScalarField, scalar>(nuName_);
+
+    const scalarField& nutw =
+        lookupPatchField<volScalarField, scalar>(nutName_);
+
+    // Velocity field at the patch
+    const fvPatchVectorField& Uw =
+        lookupPatchField<volVectorField, vector>(UName_);
+
+    // Patch normals
+    const vectorField n = patch().nf();
+
+    // Velocity patch internal field
+    const vectorField UwIn = Uw.patchInternalField();
+
+    // Velocity vector tangential to the wall
+    const vectorField UwInTang = UwIn - (UwIn & n)*n;
+    const scalarField magUwInTang = mag(UwInTang);
+
+    // Calculate tangential direction for patch cells
+    const vectorField tDir = UwInTang/(magUwInTang + SMALL);
+
+    // Magnitude of the surface normal gradient velocity
+    const scalarField magGradUw = mag(Uw.snGrad());
+
+
+    // Pressure effects. Lookup the pressure gradient field from the registry
+    // (created with pressureGradient function object)
+
+    // Pressure gradient projected on the wall parallel velocity direction
+    scalarField gradpTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("pressureGradient")
+    )
+    {
+        const volVectorField& gradP =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("pressureGradient");
+
+        // Update pressure gradient projected on the wall parallel velocity
+        gradpTang =
+            gradP.boundaryField()[this->patch().index()].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field pressureGradient not found. Neglecting pressure gradient "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Convective terms. Lookup the convection field from the registry (created
+    // with velocityConvection function object)
+
+    // Velocity convection projected on the wall parallel velocity direction
+    scalarField convectionTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("velocityConvection")
+    )
+    {
+        const volVectorField& convection =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("velocityConvection");
+
+        // Upate convection term projected on the wall parallel velocity
+        convectionTang =
+            convection.boundaryField()
+            [
+                this->patch().index()
+            ].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field velocityConvection not found. Neglecting convection "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    // Set omega and G
+    forAll(nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+        const scalar uStar= Cmu25*sqrt(k[faceCellI]);
+
+        // Note: here yPlus is actually yStar
+        const scalar yPlus = uStar*y[faceI]/nuw[faceI];
+
+        const scalar Cu =  convectionTang[faceI] + gradpTang[faceI];
+        const scalar Psi = 1.0 - Cu/(kappa_*uStar*magGradUw[faceI]);
+
+        // Kader blend for velocity gradient
+        const scalar gamma = -0.01*pow(yPlus, 4)/(1.0 + 5.0*yPlus);
+        const scalar gammaEps = - 0.001*pow(yPlus, 4)/(1.0 + yPlus);
+
+        const scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+        const scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
+
+        // Menter blend for omega
+        omega[faceI] = omegaVis*exp(gammaEps) + omegaLog*exp(1.0/gammaEps);
+
+        const scalar Gvis = k[faceCellI]*sqr(magGradUw[faceI])/omega[faceI];
+        const scalar Glog = pow(uStar, 3)/(Psi*kappa_*y[faceI]);
+
+        G[faceCellI]= Gvis*exp(gamma) + Glog*exp(1.0/gamma);
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+
+    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+}
+
+void omegaCWTWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fixedInternalValueFvPatchField<scalar>::write(os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "p", "p", pName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    omegaCWTWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaCWTWallFunction/omegaCWTWallFunctionFvPatchScalarField.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     4.0
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::omegaCWTWallFunctionFvPatchScalarField
+
+Description
+    Improved boundary condition for specific dissipation, taking into account
+    non-equilibrium effects and are sensitive to flow unsteadiness and presence
+    of the pressure gradient.
+
+    Reference (bibtex entry):
+
+    @article{popovacHanjalic2007,
+        author = {Popovac, M. and Hanjali\'{c}, K.,}
+        title = {{Compound Wall Treatment for RANS Computation of Complex
+                  Turbulent Flows and Heat Transfer}},
+        journal = {Flow Turbulence Combust},
+        year = {2007},
+        pages = {78--177}
+        DOI = {10.1007/s10494-006-9067-x}
+    }
+
+SourceFiles
+    omegaCWTWallFunctionFvPatchScalarField.C
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef omegaCWTWallFunctionFvPatchScalarField_H
+#define omegaCWTWallFunctionFvPatchScalarField_H
+
+#include "fixedInternalValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class omegaCWTWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class omegaCWTWallFunctionFvPatchScalarField
+:
+    public fixedInternalValueFvPatchScalarField
+{
+    // Private data
+
+         //- Name of pressure field
+        word pName_;
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of turbulence kinetic energy field
+        word kName_;
+
+        //- Name of turbulence generation field
+        word GName_;
+
+        //- Name of laminar viscosity field
+        word nuName_;
+
+        //- Name of turbulent viscosity field
+        word nutName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+        //- beta1 coefficient
+        scalar beta1_;
+
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("omegaCWTWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        omegaCWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        omegaCWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  omegaCWTWallFunctionFvPatchScalarField
+        //  onto a new patch
+        omegaCWTWallFunctionFvPatchScalarField
+        (
+            const omegaCWTWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        omegaCWTWallFunctionFvPatchScalarField
+        (
+            const omegaCWTWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaCWTWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        omegaCWTWallFunctionFvPatchScalarField
+        (
+            const omegaCWTWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaCWTWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaMEWTWallFunction/omegaMEWTWallFunctionFvPatchScalarField.C

@@ -0,0 +1,373 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "omegaMEWTWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void omegaMEWTWallFunctionFvPatchScalarField::checkType()
+{
+    if (!patch().isWall())
+    {
+        FatalErrorIn("omegaMEWTWallFunctionFvPatchScalarField::checkType()")
+            << "Invalid wall function specification" << nl
+            << "    Patch type for patch " << patch().name()
+            << " must be wall" << nl
+            << "    Current patch type is " << patch().type() << nl << endl
+            << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+omegaMEWTWallFunctionFvPatchScalarField::omegaMEWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF),
+    pName_("p"),
+    UName_("U"),
+    kName_("k"),
+    GName_("RASModel::G"),
+    nuName_("nu"),
+    nutName_("nut"),
+    Cmu_(0.09),
+    kappa_(0.41),
+    E_(9.8),
+    beta1_(0.075)
+{
+    checkType();
+}
+
+
+omegaMEWTWallFunctionFvPatchScalarField::omegaMEWTWallFunctionFvPatchScalarField
+(
+    const omegaMEWTWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedInternalValueFvPatchField<scalar>(ptf, p, iF, mapper),
+    pName_(ptf.pName_),
+    UName_(ptf.UName_),
+    kName_(ptf.kName_),
+    GName_(ptf.GName_),
+    nuName_(ptf.nuName_),
+    nutName_(ptf.nutName_),
+    Cmu_(ptf.Cmu_),
+    kappa_(ptf.kappa_),
+    E_(ptf.E_),
+    beta1_(ptf.beta1_)
+{
+    checkType();
+}
+
+
+omegaMEWTWallFunctionFvPatchScalarField::omegaMEWTWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedInternalValueFvPatchField<scalar>(p, iF, dict),
+    pName_(dict.lookupOrDefault<word>("p", "p")),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    kName_(dict.lookupOrDefault<word>("k", "k")),
+    GName_(dict.lookupOrDefault<word>("G", "RASModel::G")),
+    nuName_(dict.lookupOrDefault<word>("nu", "nu")),
+    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
+    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
+    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    beta1_(dict.lookupOrDefault<scalar>("beta1", 0.075))
+{
+    checkType();
+}
+
+
+omegaMEWTWallFunctionFvPatchScalarField::omegaMEWTWallFunctionFvPatchScalarField
+(
+    const omegaMEWTWallFunctionFvPatchScalarField& owfpsf
+)
+:
+    fixedInternalValueFvPatchField<scalar>(owfpsf),
+    pName_(owfpsf.pName_),
+    UName_(owfpsf.UName_),
+    kName_(owfpsf.kName_),
+    GName_(owfpsf.GName_),
+    nuName_(owfpsf.nuName_),
+    nutName_(owfpsf.nutName_),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_)
+{
+    checkType();
+}
+
+
+omegaMEWTWallFunctionFvPatchScalarField::omegaMEWTWallFunctionFvPatchScalarField
+(
+    const omegaMEWTWallFunctionFvPatchScalarField& owfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedInternalValueFvPatchField<scalar>(owfpsf, iF),
+    pName_(owfpsf.pName_),
+    UName_(owfpsf.UName_),
+    kName_(owfpsf.kName_),
+    GName_(owfpsf.GName_),
+    nuName_(owfpsf.nuName_),
+    nutName_(owfpsf.nutName_),
+    Cmu_(owfpsf.Cmu_),
+    kappa_(owfpsf.kappa_),
+    E_(owfpsf.E_),
+    beta1_(owfpsf.beta1_)
+{
+    checkType();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void omegaMEWTWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // If G field is not present, execute zero gradient evaluation
+    // HJ, 20/Mar/2011
+    if (!db().foundObject<volScalarField>(GName_))
+    {
+        InfoIn("void omegaMEWTWallFunctionFvPatchScalarField::updateCoeffs()")
+            << "Cannot access " << GName_ << " field for patch "
+            << patch().name() << ".  Evaluating as zeroGradient"
+            << endl;
+
+        fvPatchScalarField::updateCoeffs();
+        zeroGradientFvPatchScalarField::evaluate();
+
+        return;
+    }
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalarField& y = rasModel.y()[patch().index()];
+
+    volScalarField& G = const_cast<volScalarField&>
+        (db().lookupObject<volScalarField>(GName_));
+
+    // Note: omega is now a refValue and set in fixedInternalValueFvPatchField
+    // HJ, 3/Aug/2011
+    scalarField& omega = refValue();
+
+    const scalarField& nuw =
+        lookupPatchField<volScalarField, scalar>(nuName_);
+
+    const scalarField& nutw =
+        lookupPatchField<volScalarField, scalar>(nutName_);
+
+    // Velocity field at the patch
+    const fvPatchVectorField& Uw =
+        lookupPatchField<volVectorField, vector>(UName_);
+
+    // Patch normals
+    const vectorField n = patch().nf();
+
+    // Velocity patch internal field
+    const vectorField UwIn = Uw.patchInternalField();
+
+    // Velocity vector tangential to the wall
+    const vectorField UwInTang = UwIn - (UwIn & n)*n;
+    const scalarField magUwInTang = mag(UwInTang);
+
+    // Calculate tangential direction for patch cells
+    const vectorField tDir = UwInTang/magUwInTang;
+
+    // Magnitude of the surface normal gradient velocity
+    const scalarField magGradUw = mag(Uw.snGrad());
+
+
+    // Pressure effects. Lookup the pressure gradient field from the registry
+    // (created with pressureGradient function object)
+
+    // Pressure gradient projected on the wall parallel velocity direction
+    scalarField gradpTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("pressureGradient")
+    )
+    {
+        const volVectorField& gradP =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("pressureGradient");
+
+        // Update pressure gradient projected on the wall parallel velocity
+        gradpTang =
+            gradP.boundaryField()[this->patch().index()].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field pressureGradient not found. Neglecting pressure gradient "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Convective terms. Lookup the convection field from the registry (created
+    // with velocityConvection function object)
+
+    // Velocity convection projected on the wall parallel velocity direction
+    scalarField convectionTang(magGradUw.size(), 0.0);
+
+    if
+    (
+        this->dimensionedInternalField().mesh().foundObject
+        <
+            volVectorField
+        >("velocityConvection")
+    )
+    {
+        const volVectorField& convection =
+            this->dimensionedInternalField().mesh().lookupObject
+            <volVectorField>("velocityConvection");
+
+        // Upate convection term projected on the wall parallel velocity
+        convectionTang =
+            convection.boundaryField()
+            [
+                this->patch().index()
+            ].patchInternalField()
+          & tDir;
+    }
+    else
+    {
+        Info<< "Field velocityConvection not found. Neglecting convection "
+            << "effects for wall functions at patch: " << patch().name()
+            << endl;
+    }
+
+
+    // Get face cells
+    const unallocLabelList& fc = patch().faceCells();
+
+    // Set omega and G
+    forAll(nutw, faceI)
+    {
+        const label faceCellI = fc[faceI];
+
+        // Tangential stress at the wall
+        const scalar tauw = (nuw[faceI] + nutw[faceI])*magGradUw[faceI];
+
+        const scalar yPlus = sqrt(tauw)*y[faceI]/nuw[faceI];
+
+        // Velocity gradient for viscous sublayer
+		const scalar dudyVis= magGradUw[faceI];
+
+        // Velocity gradient for log layer
+        const scalar dudyLog =
+            sqrt
+            (
+                max
+                (
+                    SMALL,
+                    (gradpTang[faceI] + convectionTang[faceI])*y[faceI] +tauw
+                )
+            )/(kappa_*y[faceI]);
+
+        // Kader blending for velocity gradient
+        const scalar gamma = -0.01*pow(yPlus, 4)/(1.0 + 5.0*yPlus);
+        const scalar dudy = dudyVis*exp(gamma) + dudyLog*exp(1.0/gamma);
+
+        const scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
+        const scalar omegaLog = dudyLog/(sqrt(Cmu_));
+
+        // Menter blend for omega
+        omega[faceI] = sqrt(sqr(omegaVis) + sqr(omegaLog));
+
+        // Production term
+        G[faceCellI]= tauw*dudy;
+    }
+
+    // TODO: perform averaging for cells sharing more than one boundary face
+
+    fixedInternalValueFvPatchField<scalar>::updateCoeffs();
+}
+
+void omegaMEWTWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fixedInternalValueFvPatchField<scalar>::write(os);
+    writeEntryIfDifferent<word>(os, "U", "U", UName_);
+    writeEntryIfDifferent<word>(os, "p", "p", pName_);
+    writeEntryIfDifferent<word>(os, "k", "k", kName_);
+    writeEntryIfDifferent<word>(os, "G", "RASModel::G", GName_);
+    writeEntryIfDifferent<word>(os, "nu", "nu", nuName_);
+    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
+    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
+    os.writeKeyword("beta1") << beta1_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchScalarField,
+    omegaMEWTWallFunctionFvPatchScalarField
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaMEWTWallFunction/omegaMEWTWallFunctionFvPatchScalarField.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | foam-extend: Open Source CFD
+   \\    /   O peration     | Version:     3.2
+    \\  /    A nd           | Web:         http://www.foam-extend.org
+     \\/     M anipulation  | For copyright notice see file Copyright
+-------------------------------------------------------------------------------
+License
+    This file is part of foam-extend.
+
+    foam-extend 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.
+
+    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::incompressible::RASModels::omegaMEWTWallFunctionFvPatchScalarField
+
+Description
+    Improved boundary condition for specific dissipation, taking into account
+    convection and pressure gradient terms.
+
+    Reference (bibtex entry):
+
+    @article{sutaloThesis2017,
+        author = {\v{S}utalo, F.}
+        title = {{Evaluation of Variants of Enhanced Wall Treatment Wall
+                  Function in Turbulent Flow Simulations}},
+        school = {Faculty of Mechanical Engineering and Naval Architecture,
+                  University of Zagreb},
+        year = {2017},
+    }
+
+SourceFiles
+    omegaMEWTWallFunctionFvPatchScalarField.C
+
+Author
+    Filip Sutalo, FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef omegaMEWTWallFunctionFvPatchScalarField_H
+#define omegaMEWTWallFunctionFvPatchScalarField_H
+
+#include "fixedInternalValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace incompressible
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+           Class omegaMEWTWallFunctionFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class omegaMEWTWallFunctionFvPatchScalarField
+:
+    public fixedInternalValueFvPatchScalarField
+{
+    // Private data
+
+        //- Name of pressure field
+        word pName_;
+
+        //- Name of velocity field
+        word UName_;
+
+        //- Name of turbulence kinetic energy field
+        word kName_;
+
+        //- Name of turbulence generation field
+        word GName_;
+
+        //- Name of laminar viscosity field
+        word nuName_;
+
+        //- Name of turbulent viscosity field
+        word nutName_;
+
+        //- Cmu coefficient
+        scalar Cmu_;
+
+        //- Von Karman constant
+        scalar kappa_;
+
+        //- E coefficient
+        scalar E_;
+
+        //- beta1 coefficient
+        scalar beta1_;
+
+
+    // Private member functions
+
+        //- Check the type of the patch
+        void checkType();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("omegaMEWTWallFunction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        omegaMEWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        omegaMEWTWallFunctionFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        //  omegaMEWTWallFunctionFvPatchScalarField
+        //  onto a new patch
+        omegaMEWTWallFunctionFvPatchScalarField
+        (
+            const omegaMEWTWallFunctionFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        omegaMEWTWallFunctionFvPatchScalarField
+        (
+            const omegaMEWTWallFunctionFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaMEWTWallFunctionFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        omegaMEWTWallFunctionFvPatchScalarField
+        (
+            const omegaMEWTWallFunctionFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new omegaMEWTWallFunctionFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace incompressible
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C

@@ -180,10 +180,12 @@ tmp<volSymmTensorField> kEpsilon::devReff() const
 
 tmp<fvVectorMatrix> kEpsilon::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/kOmega/kOmega.C

@@ -189,10 +189,12 @@ tmp<volSymmTensorField> kOmega::devReff() const
 
 tmp<fvVectorMatrix> kOmega::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C

@@ -342,10 +342,12 @@ tmp<volSymmTensorField> kOmegaSST::devReff() const
 
 tmp<fvVectorMatrix> kOmegaSST::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/laminar/laminar.C

@@ -167,10 +167,12 @@ tmp<volSymmTensorField> laminar::devReff() const
 
 tmp<fvVectorMatrix> laminar::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C

@@ -245,10 +245,12 @@ tmp<volSymmTensorField> qZeta::devReff() const
 
 tmp<fvVectorMatrix> qZeta::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }
 

src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C

@@ -233,10 +233,12 @@ tmp<volSymmTensorField> realizableKE::devReff() const
 
 tmp<fvVectorMatrix> realizableKE::divDevReff() const
 {
+    const volScalarField nuEffective = nuEff();
+
     return
     (
-      - fvm::laplacian(nuEff(), U_)
+      - fvm::laplacian(nuEffective, U_)
-      - fvc::div(nuEff()*dev(T(fvc::grad(U_))))
+      - (fvc::grad(U_) & fvc::grad(nuEffective))
     );
 }