diff --git a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C b/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
index cd1c39f1fe1d436953d670eff73bbd382455a655..53cebf634b42e540671c85f5b0811d4ce915cf5e 100644
--- a/src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
+++ b/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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C
index cbbed4fdc3accd3afe76b8ba0ee940fd70169686..9f4cf26bfa69849ffbab763c7223cd1a26369b40 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutLowReWallFunction/nutLowReWallFunctionFvPatchScalarField.C
+++ b/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;
 }
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C
index 003bbc48764469f76b79d6ffb3129b31c825e68b..b513f54afa615a14f6b52e80ba7ddc690baa183b 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUWallFunction/nutUWallFunctionFvPatchScalarField.C
+++ b/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;
 }
 
 
diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
index 826a8f354fc53e2c4b5cb63e510039cced908433..cf8159e4282305b1278052cf9eb4d0c5ec6ec617 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkWallFunction/nutkWallFunctionFvPatchScalarField.C
+++ b/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;
 }
 
 
diff --git a/src/functionObjects/field/yPlus/yPlus.C b/src/functionObjects/field/yPlus/yPlus.C
index 8fc394738e729afe08a3e00706121e00f3a12b40..99895924eb003ebd2182265d0201c05e9273cdac 100644
--- a/src/functionObjects/field/yPlus/yPlus.C
+++ b/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;
     }
 
diff --git a/src/functionObjects/field/yPlus/yPlus.H b/src/functionObjects/field/yPlus/yPlus.H
index 6614ef33ff20097366a8b2e7bc5d8b95218de2e1..1e70952e535be2ad837e03d84b328c36dd204dc2 100644
--- a/src/functionObjects/field/yPlus/yPlus.H
+++ b/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