From 01903e5e8f5e9397467b4f6693fefb772e730a61 Mon Sep 17 00:00:00 2001
From: andy <a.heather@opencfd.co.uk>
Date: Thu, 23 Jul 2009 15:18:00 +0100
Subject: [PATCH] further updates to wall function calcs

---
 ...ndardRoughWallFunctionFvPatchScalarField.C | 267 ++++++++++--------
 ...ndardRoughWallFunctionFvPatchScalarField.H |   3 +
 ...asStandardWallFunctionFvPatchScalarField.C |  86 ++++--
 ...asStandardWallFunctionFvPatchScalarField.H |   3 +
 ...ndardRoughWallFunctionFvPatchScalarField.C |   2 +-
 ...asStandardWallFunctionFvPatchScalarField.C |   6 +-
 6 files changed, 208 insertions(+), 159 deletions(-)

diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
index 172cfb0db56..351f4ca80d2 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
@@ -39,104 +39,25 @@ namespace compressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mutWallFunctionFvPatchScalarField(p, iF),
-    roughnessHeight_(pTraits<scalar>::zero),
-    roughnessConstant_(pTraits<scalar>::zero),
-    roughnessFudgeFactor_(pTraits<scalar>::zero)
-{}
-
-
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
-(
-    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    mutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    roughnessHeight_(ptf.roughnessHeight_),
-    roughnessConstant_(ptf.roughnessConstant_),
-    roughnessFudgeFactor_(ptf.roughnessFudgeFactor_)
-{}
-
-
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    mutWallFunctionFvPatchScalarField(p, iF, dict),
-    roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
-    roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
-    roughnessFudgeFactor_(readScalar(dict.lookup("roughnessFudgeFactor")))
-{}
-
-
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
-(
-    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf
-)
-:
-    mutWallFunctionFvPatchScalarField(rwfpsf),
-    roughnessHeight_(rwfpsf.roughnessHeight_),
-    roughnessConstant_(rwfpsf.roughnessConstant_),
-    roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
-{}
-
-
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
-(
-    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    mutWallFunctionFvPatchScalarField(rwfpsf, iF),
-    roughnessHeight_(rwfpsf.roughnessHeight_),
-    roughnessConstant_(rwfpsf.roughnessConstant_),
-    roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 tmp<scalarField>
-mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
+(
+    const scalarField& magUp
+) const
 {
     const label patchI = patch().index();
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
-    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
-
-    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
     const scalarField& muw = rasModel.mu().boundaryField()[patchI];
-
-    scalarField magUp = mag(Uw.patchInternalField() - Uw);
-
     const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
 
-    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
-    scalarField& mutw = tmutw();
+    tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
+    scalarField& yPlus = tyPlus();
 
     if (roughnessHeight_ > 0.0)
     {
-        // Rough Walls.
+        // Rough Walls
         const scalar c_1 = 1/(90 - 2.25) + roughnessConstant_;
         static const scalar c_2 = 2.25/(90 - 2.25);
         static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
@@ -145,15 +66,15 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
         //if (KsPlusBasedOnYPlus_)
         {
             // If KsPlus is based on YPlus the extra term added to the law
-            // of the wall will depend on yPlus.
-            forAll(mutw, facei)
+            // of the wall will depend on yPlus
+            forAll(yPlus, facei)
             {
                 const scalar magUpara = magUp[facei];
                 const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
                 const scalar kappaRe = kappa_*Re;
 
-                scalar yPlus = yPlusLam;
-                const scalar ryPlusLam = 1.0/yPlus;
+                scalar yp = yPlusLam;
+                const scalar ryPlusLam = 1.0/yp;
 
                 int iter = 0;
                 scalar yPlusLast = 0.0;
@@ -171,10 +92,10 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
 
                 do
                 {
-                    yPlusLast = yPlus;
+                    yPlusLast = yp;
 
                     // The non-dimensional roughness height
-                    scalar KsPlus = yPlus*dKsPlusdYPlus;
+                    scalar KsPlus = yp*dKsPlusdYPlus;
 
                     // The extra term in the law-of-the-wall
                     scalar G = 0.0;
@@ -198,65 +119,159 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
                             (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
                     }
 
-                    scalar denom = 1.0 + log(E_*yPlus) - G - yPlusGPrime;
+                    scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
                     if (mag(denom) > VSMALL)
                     {
-                        yPlus = (kappaRe + yPlus*(1 - yPlusGPrime))/denom;
-                        if( yPlus < 0 )
-                        {
-                            yPlus = 0;
-                        }
-                    }
-                    else
-                    {
-                        // Ensure immediate end and mutw = 0
-                        yPlus = 0;
+                        yp = kappaRe + yp*(1 - yPlusGPrime))/denom;
                     }
 
                 } while
                 (
-                    mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001
+                    mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
                  && ++iter < 10
-                 && yPlus > VSMALL
+                 && yp > VSMALL
                 );
 
-                if (yPlus > yPlusLam)
-                {
-                    mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1);
-                }
+                yPlus[facei] = max(0.0, yp);
             }
         }
     }
     else
     {
         // Smooth Walls
-        forAll(mutw, facei)
+        forAll(yPlus, facei)
         {
             const scalar magUpara = magUp[facei];
             const scalar Re = rho[facei]*magUpara*y[facei]/muw[facei];
             const scalar kappaRe = kappa_*Re;
 
-            scalar yPlus = yPlusLam;
-            const scalar ryPlusLam = 1.0/yPlus;
+            scalar yp = yPlusLam;
+            const scalar ryPlusLam = 1.0/yp;
 
             int iter = 0;
             scalar yPlusLast = 0.0;
 
             do
             {
-                yPlusLast = yPlus;
-                yPlus = (kappaRe + yPlus)/(1.0 + log(E_*yPlus));
+                yPlusLast = yp;
+                yp = (kappaRe + yp)/(1.0 + log(E_*yp));
 
             } while
             (
-                mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001
+                mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
              && ++iter < 10
             );
 
-            if (yPlus > yPlusLam)
-            {
-                mutw[facei] = muw[facei]*(yPlus*yPlus/Re - 1);
-            }
+            yPlus[facei] = max(0.0, yp);
+        }
+    }
+
+    return tyPlus;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutWallFunctionFvPatchScalarField(p, iF),
+    roughnessHeight_(pTraits<scalar>::zero),
+    roughnessConstant_(pTraits<scalar>::zero),
+    roughnessFudgeFactor_(pTraits<scalar>::zero)
+{}
+
+
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    mutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    roughnessHeight_(ptf.roughnessHeight_),
+    roughnessConstant_(ptf.roughnessConstant_),
+    roughnessFudgeFactor_(ptf.roughnessFudgeFactor_)
+{}
+
+
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    mutWallFunctionFvPatchScalarField(p, iF, dict),
+    roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
+    roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
+    roughnessFudgeFactor_(readScalar(dict.lookup("roughnessFudgeFactor")))
+{}
+
+
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf
+)
+:
+    mutWallFunctionFvPatchScalarField(rwfpsf),
+    roughnessHeight_(rwfpsf.roughnessHeight_),
+    roughnessConstant_(rwfpsf.roughnessConstant_),
+    roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
+{}
+
+
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
+(
+    const mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField& rwfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    mutWallFunctionFvPatchScalarField(rwfpsf, iF),
+    roughnessHeight_(rwfpsf.roughnessHeight_),
+    roughnessConstant_(rwfpsf.roughnessConstant_),
+    roughnessFudgeFactor_(rwfpsf.roughnessFudgeFactor_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<scalarField>
+mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField& muw = rasModel.mu().boundaryField()[patchI];
+    const fvPatchScalarField& rho = rasModel.rho().boundaryField()[patchI];
+
+    scalarField magUp = mag(Uw.patchInternalField() - Uw);
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
+    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
+    scalarField& mutw = tmutw();
+
+    forAll(yPlus, facei)
+    {
+        if (yPlus[facei] > yPlusLam)
+        {
+            const scalar Re = rho[facei]*magUp[facei]*y[facei]/muw[facei];
+            mutw[facei] = muw[facei]*(sqr(yPlus[facei])/Re - 1);
         }
     }
 
@@ -267,13 +282,13 @@ mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcMut() const
 tmp<scalarField>
 mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
 {
-    notImplemented
-    (
-        "mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus()"
-        "const"
-    );
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
 
-    return tmp<scalarField>(NULL);
+    return calcYPlus(magUp);
 }
 
 
@@ -283,6 +298,8 @@ void mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::write
 ) const
 {
     fixedValueFvPatchScalarField::write(os);
+    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
     os.writeKeyword("roughnessHeight")
         << roughnessHeight_ << token::END_STATEMENT << nl;
     os.writeKeyword("roughnessConstant")
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H
index e2434c8183e..2414825d213 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardRoughWallFunction/mutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H
@@ -74,6 +74,9 @@ protected:
 
     // Protected member functions
 
+        //- Calculate yPLus
+        virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
+
         //- Calculate the turbulence viscosity
         virtual tmp<scalarField> calcMut() const;
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
index f2337dab471..493cf1598d4 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
@@ -39,6 +39,50 @@ namespace compressible
 namespace RASModels
 {
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField>
+nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
+(
+    const scalarField& magUp
+) const
+{
+    const label patchI = patch().index();
+
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
+    const scalarField& y = rasModel.y()[patchI];
+    const fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI];
+    const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
+
+    tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
+    scalarField& yPlus = tyPlus();
+
+    forAll(yPlus, faceI)
+    {
+        scalar kappaRe = kappa_*magUp[facei]*y[faceI]/(muw[faceI]/rhow[faceI]);
+
+        scalar yp = yPlusLam;
+        scalar ryPlusLam = 1.0/yp;
+
+        int iter = 0;
+        scalar yPlusLast = 0.0;
+
+        do
+        {
+            yPlusLast = yp;
+            yp = (kappaRe + yp)/(1.0 + log(E_*yp));
+
+        } while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10);
+
+        yPlus[facei] = max(0.0, yp);
+    }
+
+    return tyPlus;
+
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
@@ -108,40 +152,23 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const scalar yPlusLam = rasModel.yPlusLam(kappa_, E_);
     const scalarField& y = rasModel.y()[patchI];
-
     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-
     const scalarField magUp = mag(Uw.patchInternalField() - Uw);
-
     const fvPatchScalarField& rhow = rasModel.rho().boundaryField()[patchI];
-
     const fvPatchScalarField& muw = rasModel.mu().boundaryField()[patchI];
 
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
     tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
     scalarField& mutw = tmutw();
 
-    forAll(mutw, faceI)
+    forAll(yPlus, faceI)
     {
-        scalar magUpara = magUp[faceI];
-
-        scalar kappaRe = kappa_*magUpara*y[faceI]/(muw[faceI]/rhow[faceI]);
-
-        scalar yPlus = yPlusLam;
-        scalar ryPlusLam = 1.0/yPlus;
-
-        int iter = 0;
-        scalar yPlusLast = 0.0;
-
-        do
-        {
-            yPlusLast = yPlus;
-            yPlus = (kappaRe + yPlus)/(1.0 + log(E_*yPlus));
-
-        } while (mag(ryPlusLam*(yPlus - yPlusLast)) > 0.01 && ++iter < 10);
-
-        if (yPlus > yPlusLam)
+        if (yPlus[facei] > yPlusLam)
         {
-            mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
+            mutw[faceI] =
+                muw[faceI]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0);
         }
     }
 
@@ -152,13 +179,12 @@ mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcMut() const
 tmp<scalarField>
 mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
 {
-    notImplemented
-    (
-        "mutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() "
-        "const"
-    );
+    const label patchI = patch().index();
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
 
-    return tmp<scalarField>(NULL);
+    return calcYPlus(magUp);
 }
 
 
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H
index 6813e59667c..09ba112ce4f 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutSpalartAllmarasStandardWallFunction/mutSpalartAllmarasStandardWallFunctionFvPatchScalarField.H
@@ -60,6 +60,9 @@ protected:
 
     // Protected member functions
 
+        //- Calculate yPLus
+        virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
+
         //- Calculate the turbulence viscosity
         virtual tmp<scalarField> calcMut() const;
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
index d6b96ed5b45..f44da462b9c 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C
@@ -155,7 +155,7 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcYPlus
             do
             {
                 yPlusLast = yp;
-                yPlus = (kappaRe + yp)/(1.0 + log(E_*yp));
+                yp = (kappaRe + yp)/(1.0 + log(E_*yp));
 
             } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10);
 
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
index 25c9b87d98e..450c857a756 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C
@@ -39,7 +39,7 @@ namespace incompressible
 namespace RASModels
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 tmp<scalarField>
 nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
@@ -70,11 +70,11 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcYPlus
         do
         {
             yPlusLast = yp;
-            yPlus = (kappaRe + yp)/(1.0 + log(E_*yp));
+            yp = (kappaRe + yp)/(1.0 + log(E_*yp));
 
         } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 );
 
-        yPlus[facei] = yp;
+        yPlus[facei] = max(0.0, yp);
     }
 
     return tyPlus;
-- 
GitLab