From a1091a254aadfe8d971cf21c5e339e9232ac3b51 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Fri, 7 Jun 2019 23:24:55 +0100
Subject: [PATCH] nutkRoughWallFunction: Re-derived and implemented to improve
 stability and range of applicability

This wall function is now stable even if the near-wall cell centre distance is
less than the roughness height although it is unlikely to be very accurate if
used in this way.

The sudden change limit stabilisation has been re-implemented to avoid the
near-wall viscosity being lower limited to half the laminar viscosity which has
no physical meaning.
---
 .../nutkRoughWallFunctionFvPatchScalarField.C | 37 +++++++++----------
 1 file changed, 18 insertions(+), 19 deletions(-)

diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
index 21b1666150..7c7aab9675 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkRoughWallFunction/nutkRoughWallFunctionFvPatchScalarField.C
@@ -43,8 +43,11 @@ scalar nutkRoughWallFunctionFvPatchScalarField::fnRough
 ) const
 {
     // Return fn based on non-dimensional roughness height
-
-    if (KsPlus < 90.0)
+    if (KsPlus < 2.25)
+    {
+        return 1;
+    }
+    else if (KsPlus < 90)
     {
         return pow
         (
@@ -54,7 +57,7 @@ scalar nutkRoughWallFunctionFvPatchScalarField::fnRough
     }
     else
     {
-        return (1.0 + Cs*KsPlus);
+        return (1 + Cs*KsPlus);
     }
 }
 
@@ -84,31 +87,25 @@ tmp<scalarField> nutkRoughWallFunctionFvPatchScalarField::calcNut() const
 
     forAll(nutw, facei)
     {
-        label celli = patch().faceCells()[facei];
+        const label celli = patch().faceCells()[facei];
 
-        scalar uStar = Cmu25*sqrt(k[celli]);
-        scalar yPlus = uStar*y[facei]/nuw[facei];
-        scalar KsPlus = uStar*Ks_[facei]/nuw[facei];
-
-        scalar Edash = E_;
-        if (KsPlus > 2.25)
-        {
-            Edash /= fnRough(KsPlus, Cs_[facei]);
-        }
-
-        scalar limitingNutw = max(nutw[facei], nuw[facei]);
+        const scalar uStar = Cmu25*sqrt(k[celli]);
+        const scalar KsPlus = uStar*Ks_[facei]/nuw[facei];
+        const scalar Edash = E_/fnRough(KsPlus, Cs_[facei]);
+        const scalar yPlusMin = constant::mathematical::e/Edash;
+        const scalar yPlus = max(uStar*y[facei]/nuw[facei], yPlusMin);
 
         // To avoid oscillations limit the change in the wall viscosity
-        // which is particularly important if it temporarily becomes zero
         nutw[facei] =
             max
             (
                 min
                 (
                     nuw[facei]
-                   *(yPlus*kappa_/log(max(Edash*yPlus, 1+1e-4)) - 1),
-                    2*limitingNutw
-                ), 0.5*limitingNutw
+                   *max(yPlus*kappa_/log(max(Edash*yPlus, 1)) - 1, 0),
+                    max(2*nutw[facei], nuw[facei])
+                ),
+                0.5*nutw[facei]
             );
 
         if (debug)
@@ -116,6 +113,8 @@ tmp<scalarField> nutkRoughWallFunctionFvPatchScalarField::calcNut() const
             Info<< "yPlus = " << yPlus
                 << ", KsPlus = " << KsPlus
                 << ", Edash = " << Edash
+                << ", yPlusMin " << yPlusMin
+                << ", yPlusLam " << yPlusLam(kappa_, Edash)
                 << ", nutw = " << nutw[facei]
                 << endl;
         }