Additional wall functions: k/epsilonNonEqWallFunctions

Non-equilibrium wall functions for k and epsilon.

Author: Filip Sutalo, Merge: Vuko Vukcevic

src/turbulenceModels/incompressible/RAS/Make/files

@@ -37,6 +37,7 @@ $(nutWallFunctions)/nutMEWTWallFunction/nutMEWTWallFunctionFvPatchScalarField.C
 
 epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
 $(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
+$(epsilonWallFunctions)/epsilonNonEqWallFunction/epsilonNonEqWallFunctionFvPatchScalarField.C
 
 omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
 $(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
@@ -45,6 +46,7 @@ $(omegaWallFunctions)/omegaMEWTWallFunction/omegaMEWTWallFunctionFvPatchScalarFi
 
 kqRWallFunctions = $(wallFunctions)/kqRWallFunctions
 $(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C
+$(kqRWallFunctions)/kNonEqWallFunction/kNonEqWallFunctionFvPatchScalarField.C
 
 RWallFunctions = $(wallFunctions)/RWallFunctions
 $(RWallFunctions)/RWallFunction/RWallFunctionFvPatchSymmTensorField.C

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/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
+
+// ************************************************************************* //