epsilonWallFunction: Updated to work with both low- and high-Reynolds number turbulence models

This boundary condition provides a turbulence dissipation wall constraint
for low- and high-Reynolds number turbulence models.

The condition can be applied to wall boundaries for which it
- calculates \c epsilon and \c G
- specifies the near-wall epsilon value

where

\vartable
    epsilon | turblence dissipation field
    G       | turblence generation field
\endvartable

The model switches between laminar and turbulent functions based on the
laminar-to-turbulent y+ value derived from kappa and E.

Recent tests have shown that this formulation is more accurate than
the standard high-Reynolds number form for 10 < y+ < 30 with both
standard and continuous wall-functions.

Replaces epsilonLowReWallFunction and should be used for all
low-Reynolds number models for which the epsilonLowReWallFunction BC was
recommended.

src/TurbulenceModels/turbulenceModels/Make/files

@@ -43,7 +43,6 @@ $(nutWallFunctions)/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarFiel
 
 epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
 $(epsilonWallFunctions)/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C
-$(epsilonWallFunctions)/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C
 
 omegaWallFunctions = $(wallFunctions)/omegaWallFunctions
 $(omegaWallFunctions)/omegaWallFunction/omegaWallFunctionFvPatchScalarField.C

src/TurbulenceModels/turbulenceModels/RAS/v2f/v2f.H

@@ -41,7 +41,7 @@ Description
     Wall boundary conditions are:
 
         k       = kLowReWallFunction
-        epsilon = epsilonLowReWallFunction
+        epsilon = epsilonWallFunction
         v2      = v2WallFunction
         f       = fWallFunction
 
@@ -90,7 +90,7 @@ See also
     Foam::RASModels::v2fBase
     Foam::RASModels::kEpsilon
     Foam::kLowReWallFunctionFvPatchScalarField
-    Foam::epsilonLowReWallFunctionFvPatchScalarField
+    Foam::epsilonWallFunctionFvPatchScalarField
     Foam::v2WallFunctionFvPatchScalarField
     Foam::fWallFunctionFvPatchScalarField
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.C

@@ -1,191 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "epsilonLowReWallFunctionFvPatchScalarField.H"
-#include "turbulenceModel.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::epsilonLowReWallFunctionFvPatchScalarField::yPlusLam
-(
-    const scalar kappa,
-    const scalar E
-)
-{
-    scalar ypl = 11.0;
-
-    for (int i=0; i<10; i++)
-    {
-        ypl = log(max(E*ypl, 1))/kappa;
-    }
-
-    return ypl;
-}
-
-
-void Foam::epsilonLowReWallFunctionFvPatchScalarField::calculate
-(
-    const turbulenceModel& turbModel,
-    const List<scalar>& cornerWeights,
-    const fvPatch& patch,
-    scalarField& G0,
-    scalarField& epsilon0
-)
-{
-    const label patchi = patch.index();
-
-    const scalarField& y = turbModel.y()[patchi];
-
-    const scalar Cmu25 = pow025(Cmu_);
-    const scalar Cmu75 = pow(Cmu_, 0.75);
-
-    const tmp<volScalarField> tk = turbModel.k();
-    const volScalarField& k = tk();
-
-    const tmp<scalarField> tnuw = turbModel.nu(patchi);
-    const scalarField& nuw = tnuw();
-
-    const tmp<scalarField> tnutw = turbModel.nut(patchi);
-    const scalarField& nutw = tnutw();
-
-    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
-
-    const scalarField magGradUw(mag(Uw.snGrad()));
-
-    const DimensionedField<scalar, volMesh>& G =
-        db().lookupObject<DimensionedField<scalar, volMesh>>
-        (
-            turbModel.GName()
-        );
-
-    // Set epsilon and G
-    forAll(nutw, facei)
-    {
-        label celli = patch.faceCells()[facei];
-
-        scalar yPlus = Cmu25*sqrt(k[celli])*y[facei]/nuw[facei];
-
-        scalar w = cornerWeights[facei];
-
-        if (yPlus > yPlusLam_)
-        {
-            epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
-
-            G0[celli] +=
-                w
-               *(nutw[facei] + nuw[facei])
-               *magGradUw[facei]
-               *Cmu25*sqrt(k[celli])
-               /(kappa_*y[facei]);
-        }
-        else
-        {
-            epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
-            G0[celli] += G[celli];
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::epsilonLowReWallFunctionFvPatchScalarField::
-epsilonLowReWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    epsilonWallFunctionFvPatchScalarField(p, iF),
-    yPlusLam_(yPlusLam(kappa_, E_))
-{}
-
-
-Foam::epsilonLowReWallFunctionFvPatchScalarField::
-epsilonLowReWallFunctionFvPatchScalarField
-(
-    const epsilonLowReWallFunctionFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    epsilonWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    yPlusLam_(ptf.yPlusLam_)
-{}
-
-
-Foam::epsilonLowReWallFunctionFvPatchScalarField::
-epsilonLowReWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    epsilonWallFunctionFvPatchScalarField(p, iF, dict),
-    yPlusLam_(yPlusLam(kappa_, E_))
-{}
-
-
-Foam::epsilonLowReWallFunctionFvPatchScalarField::
-epsilonLowReWallFunctionFvPatchScalarField
-(
-    const epsilonLowReWallFunctionFvPatchScalarField& ewfpsf
-)
-:
-    epsilonWallFunctionFvPatchScalarField(ewfpsf),
-    yPlusLam_(ewfpsf.yPlusLam_)
-{}
-
-
-Foam::epsilonLowReWallFunctionFvPatchScalarField::
-epsilonLowReWallFunctionFvPatchScalarField
-(
-    const epsilonLowReWallFunctionFvPatchScalarField& ewfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    epsilonWallFunctionFvPatchScalarField(ewfpsf, iF),
-    yPlusLam_(ewfpsf.yPlusLam_)
-{}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    makePatchTypeField
-    (
-        fvPatchScalarField,
-        epsilonLowReWallFunctionFvPatchScalarField
-    );
-}
-
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonLowReWallFunction/epsilonLowReWallFunctionFvPatchScalarField.H

@@ -1,190 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::epsilonLowReWallFunctionFvPatchScalarField
-
-Group
-    grpWallFunctions
-
-Description
-    This boundary condition provides a turbulence dissipation wall function
-    condition for low- and high-Reynolds number turbulent flow cases.
-
-    The condition can be applied to wall boundaries, whereby it inserts near
-    wall epsilon values directly into the epsilon equation to act as a
-    constraint.
-
-    The model operates in two modes, based on the computed laminar-to-turbulent
-    switch-over y+ value derived from kappa and E.
-
-Usage
-    \table
-        Property     | Description             | Required    | Default value
-        Cmu          | model coefficient       | no          | 0.09
-        kappa        | Von Karman constant     | no          | 0.41
-        E            | model coefficient       | no          | 9.8
-    \endtable
-
-    Example of the boundary condition specification:
-    \verbatim
-    <patchName>
-    {
-        type            epsilonLowReWallFunction;
-    }
-    \endverbatim
-
-See also
-    Foam::epsilonWallFunctionFvPatchScalarField
-
-SourceFiles
-    epsilonLowReWallFunctionFvPatchScalarField.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef epsilonLowReWallFunctionFvPatchScalarField_H
-#define epsilonLowReWallFunctionFvPatchScalarField_H
-
-#include "epsilonWallFunctionFvPatchScalarField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-           Class epsilonLowReWallFunctionFvPatchScalarField Declaration
-\*---------------------------------------------------------------------------*/
-
-class epsilonLowReWallFunctionFvPatchScalarField
-:
-    public epsilonWallFunctionFvPatchScalarField
-{
-
-protected:
-
-    // Protected data
-
-        //- Y+ at the edge of the laminar sublayer
-        scalar yPlusLam_;
-
-
-    // Protected Member Functions
-
-        //- Calculate the Y+ at the edge of the laminar sublayer
-        scalar yPlusLam(const scalar kappa, const scalar E);
-
-        //- Calculate the epsilon and G
-        virtual void calculate
-        (
-            const turbulenceModel& turbulence,
-            const List<scalar>& cornerWeights,
-            const fvPatch& patch,
-            scalarField& G,
-            scalarField& epsilon
-        );
-
-
-public:
-
-    //- Runtime type information
-    TypeName("epsilonLowReWallFunction");
-
-
-    // Constructors
-
-        //- Construct from patch and internal field
-        epsilonLowReWallFunctionFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&
-        );
-
-        //- Construct from patch, internal field and dictionary
-        epsilonLowReWallFunctionFvPatchScalarField
-        (
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const dictionary&
-        );
-
-        //- Construct by mapping given
-        //  epsilonLowReWallFunctionFvPatchScalarField
-        //  onto a new patch
-        epsilonLowReWallFunctionFvPatchScalarField
-        (
-            const epsilonLowReWallFunctionFvPatchScalarField&,
-            const fvPatch&,
-            const DimensionedField<scalar, volMesh>&,
-            const fvPatchFieldMapper&
-        );
-
-        //- Construct as copy
-        epsilonLowReWallFunctionFvPatchScalarField
-        (
-            const epsilonLowReWallFunctionFvPatchScalarField&
-        );
-
-        //- Construct and return a clone
-        virtual tmp<fvPatchScalarField> clone() const
-        {
-            return tmp<fvPatchScalarField>
-            (
-                new epsilonLowReWallFunctionFvPatchScalarField(*this)
-            );
-        }
-
-        //- Construct as copy setting internal field reference
-        epsilonLowReWallFunctionFvPatchScalarField
-        (
-            const epsilonLowReWallFunctionFvPatchScalarField&,
-            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 epsilonLowReWallFunctionFvPatchScalarField(*this, iF)
-            );
-        }
-
-    //- Destructor
-    virtual ~epsilonLowReWallFunctionFvPatchScalarField()
-    {}
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "epsilonWallFunctionFvPatchScalarField.H"
+#include "nutWallFunctionFvPatchScalarField.H"
 #include "turbulenceModel.H"
 #include "fvPatchFieldMapper.H"
 #include "fvMatrix.H"
@@ -173,7 +174,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
     scalarField& epsilon0
 )
 {
-    // accumulate all of the G and epsilon contributions
+    // Accumulate all of the G and epsilon contributions
     forAll(cornerWeights_, patchi)
     {
         if (!cornerWeights_[patchi].empty())
@@ -186,7 +187,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
         }
     }
 
-    // apply zero-gradient condition for epsilon
+    // Apply zero-gradient condition for epsilon
     forAll(cornerWeights_, patchi)
     {
         if (!cornerWeights_[patchi].empty())
@@ -201,49 +202,63 @@ void Foam::epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
 
 void Foam::epsilonWallFunctionFvPatchScalarField::calculate
 (
-    const turbulenceModel& turbulence,
+    const turbulenceModel& turbModel,
     const List<scalar>& cornerWeights,
     const fvPatch& patch,
-    scalarField& G,
+    scalarField& G0,
-    scalarField& epsilon
+    scalarField& epsilon0
 )
 {
     const label patchi = patch.index();
 
-    const scalarField& y = turbulence.y()[patchi];
+    const scalarField& y = turbModel.y()[patchi];
 
     const scalar Cmu25 = pow025(Cmu_);
     const scalar Cmu75 = pow(Cmu_, 0.75);
 
-    const tmp<volScalarField> tk = turbulence.k();
+    const tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
 
-    const tmp<scalarField> tnuw = turbulence.nu(patchi);
+    const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
 
-    const tmp<scalarField> tnutw = turbulence.nut(patchi);
+    const tmp<scalarField> tnutw = turbModel.nut(patchi);
     const scalarField& nutw = tnutw();
 
-    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchi];
+    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
 
     const scalarField magGradUw(mag(Uw.snGrad()));
 
+    typedef DimensionedField<scalar, volMesh> FieldType;
+    const FieldType& G = db().lookupObject<FieldType>(turbModel.GName());
+
     // Set epsilon and G
     forAll(nutw, facei)
     {
-        label celli = patch.faceCells()[facei];
+        const label celli = patch.faceCells()[facei];
 
-        scalar w = cornerWeights[facei];
+        const scalar yPlus = Cmu25*y[facei]*sqrt(k[celli])/nuw[facei];
 
-        epsilon[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
+        const scalar w = cornerWeights[facei];
 
-        G[celli] +=
+        if (yPlus > yPlusLam_)
+        {
+            epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);
+
+            G0[celli] +=
                 w
                *(nutw[facei] + nuw[facei])
                *magGradUw[facei]
                *Cmu25*sqrt(k[celli])
                /(kappa_*y[facei]);
         }
+        else
+        {
+            epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
+
+            G0[celli] += w*G[celli];
+        }
+    }
 }
 
 
@@ -260,6 +275,7 @@ epsilonWallFunctionFvPatchScalarField
     Cmu_(0.09),
     kappa_(0.41),
     E_(9.8),
+    yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
     G_(),
     epsilon_(),
     initialised_(false),
@@ -283,6 +299,7 @@ epsilonWallFunctionFvPatchScalarField
     Cmu_(ptf.Cmu_),
     kappa_(ptf.kappa_),
     E_(ptf.E_),
+    yPlusLam_(ptf.yPlusLam_),
     G_(),
     epsilon_(),
     initialised_(false),
@@ -305,6 +322,7 @@ epsilonWallFunctionFvPatchScalarField
     Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     E_(dict.lookupOrDefault<scalar>("E", 9.8)),
+    yPlusLam_(nutWallFunctionFvPatchScalarField::yPlusLam(kappa_, E_)),
     G_(),
     epsilon_(),
     initialised_(false),
@@ -313,7 +331,7 @@ epsilonWallFunctionFvPatchScalarField
 {
     checkType();
 
-    // apply zero-gradient condition on start-up
+    // Apply zero-gradient condition on start-up
     this->operator==(patchInternalField());
 }
 
@@ -328,6 +346,7 @@ epsilonWallFunctionFvPatchScalarField
     Cmu_(ewfpsf.Cmu_),
     kappa_(ewfpsf.kappa_),
     E_(ewfpsf.E_),
+    yPlusLam_(ewfpsf.yPlusLam_),
     G_(),
     epsilon_(),
     initialised_(false),
@@ -349,6 +368,7 @@ epsilonWallFunctionFvPatchScalarField
     Cmu_(ewfpsf.Cmu_),
     kappa_(ewfpsf.kappa_),
     E_(ewfpsf.E_),
+    yPlusLam_(ewfpsf.yPlusLam_),
     G_(),
     epsilon_(),
     initialised_(false),
@@ -487,7 +507,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::updateWeightedCoeffs
 
     scalarField& epsilonf = *this;
 
-    // only set the values if the weights are > tolerance
+    // Only set the values if the weights are > tolerance
     forAll(weights, facei)
     {
         scalar w = weights[facei];
@@ -545,7 +565,7 @@ void Foam::epsilonWallFunctionFvPatchScalarField::manipulateMatrix
 
     forAll(weights, facei)
     {
-        // only set the values if the weights are > tolerance
+        // Anly set the values if the weights are > tolerance
         if (weights[facei] > tolerance_)
         {
             nConstrainedCells++;

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.H

@@ -28,13 +28,12 @@ Group
     grpWallFunctions
 
 Description
-    This boundary condition provides a turbulence dissipation wall function
+    This boundary condition provides a turbulence dissipation wall constraint
-    condition for high Reynolds number, turbulent flow cases.
+    for low- and high-Reynolds number turbulence models.
 
-    The condition can be applied to wall boundaries, whereby it
+    The condition can be applied to wall boundaries for which it
     - calculates \c epsilon and \c G
-    - inserts near wall epsilon values directly into the epsilon equation
+    - specifies the near-wall epsilon value
-        to act as a constraint
 
     where
 
@@ -43,6 +42,9 @@ Description
         G       | turblence generation field
     \endvartable
 
+    The model switches between laminar and turbulent functions based on the
+    laminar-to-turbulent y+ value derived from kappa and E.
+
 Usage
     \table
         Property     | Description             | Required    | Default value
@@ -61,6 +63,7 @@ Usage
 
 See also
     Foam::fixedInternalValueFvPatchField
+    Foam::omegaWallFunctionFvPatchScalarField
 
 SourceFiles
     epsilonWallFunctionFvPatchScalarField.C
@@ -103,6 +106,9 @@ protected:
         //- E coefficient
         scalar E_;
 
+        //- y+ at the edge of the laminar sublayer
+        scalar yPlusLam_;
+
         //- Local copy of turbulence G field
         scalarField G_;