Merge branch 'issue-1838-yPlus' into 'develop'

BUG: fix yPlus predictions of nut wall functions

See merge request !484

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 #include "ShihQuadraticKE.H"
 #include "bound.H"
 #include "wallFvPatch.H"
-#include "nutkWallFunctionFvPatchScalarField.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,11 +116,16 @@ yPlus() const
         )
     );
     const scalarField& y = turbModel.y()[patchi];
+
     const tmp<scalarField> tnuw = turbModel.nu(patchi);
     const scalarField& nuw = tnuw();
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
 
-    return y*sqrt(nuw*mag(Uw.snGrad()))/nuw;
+    tmp<scalarField> tnuEff = turbModel.nuEff(patchi);
+    const scalarField& nuEff = tnuEff();
+
+    return y*sqrt(nuEff*mag(Uw.snGrad()))/nuw;
 }
 
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -183,7 +183,8 @@ Foam::tmp<Foam::scalarField>
 Foam::nutUWallFunctionFvPatchScalarField::yPlus() const
 {
     const label patchi = patch().index();
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -191,10 +192,33 @@ Foam::nutUWallFunctionFvPatchScalarField::yPlus() const
             internalField().group()
         )
     );
+
+    const scalarField& y = turbModel.y()[patchi];
+
+    tmp<scalarField> tnuw = turbModel.nu(patchi);
+    const scalarField& nuw = tnuw();
+
+    tmp<scalarField> tnuEff = turbModel.nuEff(patchi);
+    const scalarField& nuEff = tnuEff();
+
     const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
     const scalarField magUp(mag(Uw.patchInternalField() - Uw));
+    const scalarField magGradUw(mag(Uw.snGrad()));
 
-    return calcYPlus(magUp);
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus.ref();
+
+    forAll(yPlus, facei)
+    {
+        if (yPlusLam_ > yPlus[facei])
+        {
+            // viscous sublayer
+            yPlus[facei] =
+                y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nuw[facei];
+        }
+    }
+
+    return tyPlus;
 }
 
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016, 2019 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,7 +33,6 @@ License
 #include "wallFvPatch.H"
 #include "addToRunTimeSelectionTable.H"
 
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 Foam::tmp<Foam::scalarField> Foam::nutkWallFunctionFvPatchScalarField::
@@ -142,7 +141,7 @@ yPlus() const
 {
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
+    const auto& turbModel = db().lookupObject<turbulenceModel>
     (
         IOobject::groupName
         (
@@ -153,12 +152,39 @@ yPlus() const
 
     const scalarField& y = turbModel.y()[patchi];
 
-    const tmp<volScalarField> tk = turbModel.k();
+    tmp<volScalarField> tk = turbModel.k();
     const volScalarField& k = tk();
-    tmp<scalarField> kwc = k.boundaryField()[patchi].patchInternalField();
-    tmp<scalarField> nuw = turbModel.nu(patchi);
+    tmp<scalarField> tkwc = k.boundaryField()[patchi].patchInternalField();
+    const scalarField& kwc = tkwc();
+
+    tmp<scalarField> tnuw = turbModel.nu(patchi);
+    const scalarField& nuw = tnuw();
+
+    tmp<scalarField> tnuEff = turbModel.nuEff(patchi);
+    const scalarField& nuEff = tnuEff();
+
+    const fvPatchVectorField& Uw = U(turbModel).boundaryField()[patchi];
+    const scalarField magGradUw(mag(Uw.snGrad()));
+
+    const scalar Cmu25 = pow025(Cmu_);
+
+    auto tyPlus = tmp<scalarField>::New(patch().size(), Zero);
+    auto& yPlus = tyPlus.ref();
+
+    forAll(yPlus, facei)
+    {
+        // inertial sublayer
+        yPlus[facei] = Cmu25*y[facei]*sqrt(kwc[facei])/nuw[facei];
+
+        if (yPlusLam_ > yPlus[facei])
+        {
+            // viscous sublayer
+            yPlus[facei] =
+                y[facei]*sqrt(nuEff[facei]*magGradUw[facei])/nuw[facei];
+        }
+    }
 
-    return pow025(Cmu_)*y*sqrt(kwc)/nuw;
+    return tyPlus;
 }
 
 

src/functionObjects/field/yPlus/yPlus.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -70,7 +70,8 @@ Foam::functionObjects::yPlus::yPlus
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
-    writeFile(obr_, name, typeName, dict)
+    writeFile(obr_, name, typeName, dict),
+    useWallFunction_(true)
 {
     read(dict);
 
@@ -101,10 +102,14 @@ Foam::functionObjects::yPlus::yPlus
 
 bool Foam::functionObjects::yPlus::read(const dictionary& dict)
 {
-    fvMeshFunctionObject::read(dict);
-    writeFile::read(dict);
+    if (fvMeshFunctionObject::read(dict) && writeFile::read(dict))
+    {
+        useWallFunction_ = dict.getOrDefault("useWallFunction", true);
 
-    return true;
+        return true;
+    }
+
+    return false;
 }
 
 
@@ -139,7 +144,11 @@ bool Foam::functionObjects::yPlus::execute()
         {
             const fvPatch& patch = patches[patchi];
 
-            if (isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi]))
+            if
+            (
+                isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi])
+             && useWallFunction_
+            )
             {
                 const nutWallFunctionFvPatchScalarField& nutPf =
                     dynamic_cast<const nutWallFunctionFvPatchScalarField&>
@@ -166,6 +175,14 @@ bool Foam::functionObjects::yPlus::execute()
         WarningInFunction
             << "Unable to find turbulence model in the "
             << "database: yPlus will not be calculated" << endl;
+
+        if (postProcess)
+        {
+            WarningInFunction
+                << "Please try to use the solver option -postProcess, e.g.:"
+                << " <solver> -postProcess -func yPlus" << endl;
+        }
+
         return false;
     }
 

src/functionObjects/field/yPlus/yPlus.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,9 @@ Usage
         type        yPlus;
         libs        (fieldFunctionObjects);
 
+        // Optional entries
+        useWallFunction    true;
+
         // Optional (inherited) entries
         ...
     }
@@ -60,6 +63,10 @@ Usage
       Property   | Description                        | Type | Req'd | Dflt
       type       | Type name: yPlus                   | word |  yes  | -
       libs       | Library name: fieldFunctionObjects | word |  yes  | -
+      useWallFunction | Flag to use the local expressions of <!--
+                      --> the selected nut wall function to compute yPlus, <!--
+                      --> otherwise directly compute yPlus from flow field <!--
+                      -->                             | bool | no    | true
     \endtable
 
     The inherited entries are elaborated in:
@@ -106,6 +113,12 @@ class yPlus
     public fvMeshFunctionObject,
     public writeFile
 {
+    // Private Data
+
+        //- Flag to use the wall function expressions to compute yPlus
+        bool useWallFunction_;
+
+
     // Private Member Functions
 
         //- File header information