Additional wall functions: omega/nutMEWTWallFunctions

Modified enhanced wall treatment by Sutalo, wall functions are
sensitive to convection and pressure gradient effects.

Author: Filip Sutalo, Merge: Vuko Vukcevic

src/turbulenceModels/incompressible/RAS/Make/files

@@ -33,6 +33,7 @@ $(nutWallFunctions)/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasSta
 $(nutWallFunctions)/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C
 $(nutWallFunctions)/nutCWTWallFunction/nutCWTWallFunctionFvPatchScalarField.C
+$(nutWallFunctions)/nutMEWTWallFunction/nutMEWTWallFunctionFvPatchScalarField.C
 
 epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
 $(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
@@ -40,6 +41,7 @@ $(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarFiel
 omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
 $(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C
 $(omegaWallFunctions)/omegaCWTWallFunction/omegaCWTWallFunctionFvPatchScalarField.C
+$(omegaWallFunctions)/omegaMEWTWallFunction/omegaMEWTWallFunctionFvPatchScalarField.C
 
 kqRWallFunctions = $(wallFunctions)/kqRWallFunctions
 $(kqRWallFunctions)/kqRWallFunction/kqRWallFunctionFvPatchFields.C

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

@@ -118,7 +118,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("nutCWTWallTreatment");
+    TypeName("nutCWTWallFunction");
 
 
     // Constructors

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

@@ -0,0 +1,384 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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;
+
+    // Wall-velocity vector field tangential to the wall
+    const vectorField UwTang = Uw - (Uw & n)*n;
+    const scalarField magUwTang = mag(UwTang);
+
+
+    // Pressure terms
+    const volScalarField& p =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            volScalarField
+        >(pName_);
+
+    // Pressure gradient
+    const volVectorField gradp = fvc::grad(p);
+
+    // Pressure gradient in wall adjacent cell
+    const vectorField gradPIn =
+        gradp.boundaryField()[this->patch().index()].patchInternalField();
+
+    // Pressure gradient projected on the wall parallel velocity
+    const scalarField gradpTang= gradPIn & tDir;
+
+
+    // Convective terms
+    const volVectorField& U =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            volVectorField
+        >(UName_);
+
+    const surfaceScalarField& phi =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            surfaceScalarField
+        >("phi");
+
+    const volVectorField convection = fvc::div(phi, U);
+
+    const vectorField convectionIn =
+        convection.boundaryField()[this->patch().index()].patchInternalField();
+
+    // Convection term projected on the wall parallel velocity
+    const scalarField convectionTang = convectionIn & tDir;
+
+    // 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

@@ -131,7 +131,7 @@ omegaCWTWallFunctionFvPatchScalarField::omegaCWTWallFunctionFvPatchScalarField
 )
 :
     fixedInternalValueFvPatchField<scalar>(owfpsf),
-     pName_(owfpsf.pName_),
+    pName_(owfpsf.pName_),
     UName_(owfpsf.UName_),
     kName_(owfpsf.kName_),
     GName_(owfpsf.GName_),
@@ -232,6 +232,7 @@ void omegaCWTWallFunctionFvPatchScalarField::updateCoeffs()
     // Magnitude of the surface normal gradient velocity
     const scalarField magGradUw = mag(Uw.snGrad());
 
+
     // Pressure effects
     const volScalarField& p =
         this->dimensionedInternalField().mesh().lookupObject
@@ -244,9 +245,10 @@ void omegaCWTWallFunctionFvPatchScalarField::updateCoeffs()
     const vectorField gradPIn =
         gradp.boundaryField()[this->patch().index()].patchInternalField();
 
-    // Pressure gradient projected on the wall parallel velocity
+    // Pressure gradient projected on the wall parallel velocity direction
     const scalarField gradpTang= gradPIn & tDir;
 
+
     // Convective terms
     const volVectorField& U =
         this->dimensionedInternalField().mesh().lookupObject

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

@@ -113,7 +113,7 @@ class omegaCWTWallFunctionFvPatchScalarField
 public:
 
     //- Runtime type information
-    TypeName("omegaCWTWallTreatment");
+    TypeName("omegaCWTWallFunction");
 
 
     // Constructors

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

@@ -0,0 +1,345 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    const volScalarField& p =
+        this->dimensionedInternalField().mesh().lookupObject
+        <volScalarField>(pName_);
+
+    // Pressure gradient
+    const volVectorField gradp = fvc::grad(p);
+
+    // Pressure gradient in wall adjacent cell
+    const vectorField gradPIn =
+        gradp.boundaryField()[this->patch().index()].patchInternalField();
+
+    // Pressure gradient projected on the wall parallel velocity direction
+    const scalarField gradpTang= gradPIn & tDir;
+
+
+    // Convective terms
+    const volVectorField& U =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            volVectorField
+        >(UName_);
+
+    const surfaceScalarField& phi =
+        this->dimensionedInternalField().mesh().lookupObject
+        <
+            surfaceScalarField
+        >("phi");
+
+    const volVectorField convection = fvc::div(phi, U);
+
+    const vectorField convectionIn =
+        convection.boundaryField()[this->patch().index()].patchInternalField();
+
+    // Convection term projected on the wall parallel velocity
+    const scalarField convectionTang = convectionIn & tDir;
+
+    // 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
+
+// ************************************************************************* //