adding yPlus calcs for incompressible S-A models

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.C

@@ -39,6 +39,134 @@ namespace incompressible
 namespace RASModels
 {
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+tmp<scalarField>
+nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::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 scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
+    scalarField& yPlus = tyPlus();
+
+    if (roughnessHeight_ > 0.0)
+    {
+        // 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);
+        static const scalar c_4 = c_3*log(2.25);
+
+        //if (KsPlusBasedOnYPlus_)
+        {
+            // If KsPlus is based on YPlus the extra term added to the law
+            // of the wall will depend on yPlus
+            forAll(yPlus, facei)
+            {
+                const scalar magUpara = magUp[facei];
+                const scalar Re = magUpara*y[facei]/nuw[facei];
+                const scalar kappaRe = kappa_*Re;
+
+                scalar yp = yPlusLam;
+                const scalar ryPlusLam = 1.0/yp;
+
+                int iter = 0;
+                scalar yPlusLast = 0.0;
+                scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
+
+                // Enforce the roughnessHeight to be less than the distance to
+                // the first cell centre
+                if (dKsPlusdYPlus > 1)
+                {
+                    dKsPlusdYPlus = 1;
+                }
+
+                // Additional tuning parameter (fudge factor) - nominally = 1
+                dKsPlusdYPlus *= roughnessFudgeFactor_;
+
+                do
+                {
+                    yPlusLast = yp;
+
+                    // The non-dimensional roughness height
+                    scalar KsPlus = yp*dKsPlusdYPlus;
+
+                    // The extra term in the law-of-the-wall
+                    scalar G = 0.0;
+
+                    scalar yPlusGPrime = 0.0;
+
+                    if (KsPlus >= 90)
+                    {
+                        const scalar t_1 = 1 + roughnessConstant_*KsPlus;
+                        G = log(t_1);
+                        yPlusGPrime = roughnessConstant_*KsPlus/t_1;
+                    }
+                    else if (KsPlus > 2.25)
+                    {
+                        const scalar t_1 = c_1*KsPlus - c_2;
+                        const scalar t_2 = c_3*log(KsPlus) - c_4;
+                        const scalar sint_2 = sin(t_2);
+                        const scalar logt_1 = log(t_1);
+                        G = logt_1*sint_2;
+                        yPlusGPrime =
+                            (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
+                    }
+
+                    scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
+                    if (mag(denom) > VSMALL)
+                    {
+                        yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
+                    }
+                } while
+                (
+                    mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
+                 && ++iter < 10
+                 && yp > VSMALL
+                );
+
+                yPlus[facei] = yp;
+            }
+        }
+    }
+    else
+    {
+        // Smooth Walls
+        forAll(yPlus, facei)
+        {
+            const scalar magUpara = magUp[facei];
+            const scalar Re = magUpara*y[facei]/nuw[facei];
+            const scalar kappaRe = kappa_*Re;
+
+            scalar yp = yPlusLam;
+            const scalar ryPlusLam = 1.0/yp;
+
+            int iter = 0;
+            scalar yPlusLast = 0.0;
+
+            do
+            {
+                yPlusLast = yp;
+                yPlus = (kappaRe + yp)/(1.0 + log(E_*yp));
+
+            } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.0001 && ++iter < 10);
+
+            yPlus[facei] = yp;
+        }
+    }
+
+    return tyPlus;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::
@@ -123,131 +251,24 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() 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& nuw = rasModel.nu().boundaryField()[patchI];
 
     // The flow velocity at the adjacent cell centre
-    scalarField magUp = mag(Uw.patchInternalField() - Uw);
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
+
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
 
     tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
     scalarField& nutw = tnutw();
 
-    if (roughnessHeight_ > 0.0)
+    forAll(yPlus, facei)
     {
-        // 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);
-        static const scalar c_4 = c_3*log(2.25);
-
-        //if (KsPlusBasedOnYPlus_)
+        if (yPlus[facei] > yPlusLam)
         {
-            // If KsPlus is based on YPlus the extra term added to the law
-            // of the wall will depend on yPlus.
-            forAll(nutw, facei)
-            {
-                const scalar magUpara = magUp[facei];
-                const scalar Re = magUpara*y[facei]/nuw[facei];
-                const scalar kappaRe = kappa_*Re;
-
-                scalar yPlus = yPlusLam;
-                const scalar ryPlusLam = 1.0/yPlus;
-
-                int iter = 0;
-                scalar yPlusLast = 0.0;
-                scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
-
-                // Enforce the roughnessHeight to be less than the distance to
-                // the first cell centre
-                if (dKsPlusdYPlus > 1)
-                {
-                    dKsPlusdYPlus = 1;
-                }
-
-                // Additional tuning parameter (fudge factor) - nominally = 1
-                dKsPlusdYPlus *= roughnessFudgeFactor_;
-
-                do
-                {
-                    yPlusLast = yPlus;
-
-                    // The non-dimensional roughness height.
-                    scalar KsPlus = yPlus*dKsPlusdYPlus;
-
-                    // The extra term in the law-of-the-wall.
-                    scalar G = 0.0;
-
-                    scalar yPlusGPrime = 0.0;
-
-                    if (KsPlus >= 90)
-                    {
-                        const scalar t_1 = 1 + roughnessConstant_*KsPlus;
-                        G = log(t_1);
-                        yPlusGPrime = roughnessConstant_*KsPlus/t_1;
-                    }
-                    else if (KsPlus > 2.25)
-                    {
-                        const scalar t_1 = c_1*KsPlus - c_2;
-                        const scalar t_2 = c_3*log(KsPlus) - c_4;
-                        const scalar sint_2 = sin(t_2);
-                        const scalar logt_1 = log(t_1);
-                        G = logt_1*sint_2;
-                        yPlusGPrime =
-                            (c_1*sint_2*KsPlus/t_1) + (c_3*logt_1*cos(t_2));
-                    }
-
-                    scalar denom = 1.0 + log(E_*yPlus) - G - yPlusGPrime;
-                    if (mag(denom) > VSMALL)
-                    {
-                        yPlus = (kappaRe + yPlus*(1 - yPlusGPrime))/denom;
-                    }
-                    else
-                    {
-                        // Ensure immediate end and nutw = 0.
-                        yPlus = 0;
-                    }
-
-                } while
-                (
-                    mag(ryPlusLam*(yPlus - yPlusLast)) > 0.0001
-                 && ++iter < 10
-                 && yPlus > VSMALL
-                );
-
-                if (yPlus > yPlusLam)
-                {
-                    nutw[facei] = nuw[facei]*(yPlus*yPlus/Re - 1);
-                }
-            }
-        }
-    }
-    else
-    {
-        // Smooth Walls.
-        forAll(nutw, facei)
-        {
-            const scalar magUpara = magUp[facei];
-            const scalar Re = magUpara*y[facei]/nuw[facei];
-            const scalar kappaRe = kappa_*Re;
-
-            scalar yPlus = yPlusLam;
-            const 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.0001 && ++iter < 10);
-
-            if (yPlus > yPlusLam)
-            {
-                nutw[facei] = nuw[facei]*(yPlus*yPlus/Re - 1);
-            }
+            const scalar Re = magUp[facei]*y[facei]/nuw[facei];
+            nutw[facei] = nuw[facei]*(sqr(yPlus[facei])/Re - 1);
         }
     }
 
@@ -258,13 +279,13 @@ nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::calcNut() const
 tmp<scalarField>
 nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus() const
 {
-    notImplemented
-    (
-        "nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField::yPlus()"
-        "const"
-    );
+    const label patchI = patch().index();
 
-    return tmp<scalarField>(NULL);
+    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
+    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
+    const scalarField magUp = mag(Uw.patchInternalField() - Uw);
+
+    return calcYPlus(magUp);
 }
 
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardRoughWallFunction/nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField.H

@@ -72,6 +72,9 @@ class nutSpalartAllmarasStandardRoughWallFunctionFvPatchScalarField
 
     // Protected member functions
 
+        //- Calculate yPLus
+        virtual tmp<scalarField> calcYPlus(const scalarField& magUp) const;
+
         //- Calculate the turbulence viscosity
         virtual tmp<scalarField> calcNut() const;
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.C

@@ -39,6 +39,48 @@ namespace incompressible
 namespace RASModels
 {
 
+// * * * * * * * * * * * * * Private 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 scalarField& nuw = rasModel.nu().boundaryField()[patchI];
+
+    tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
+    scalarField& yPlus = tyPlus();
+
+    forAll(yPlus, facei)
+    {
+        scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei];
+
+        scalar yp = yPlusLam;
+        scalar ryPlusLam = 1.0/yp;
+
+        int iter = 0;
+        scalar yPlusLast = 0.0;
+
+        do
+        {
+            yPlusLast = yp;
+            yPlus = (kappaRe + yp)/(1.0 + log(E_*yp));
+
+        } while(mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 );
+
+        yPlus[facei] = yp;
+    }
+
+    return tyPlus;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::
@@ -107,37 +149,22 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() 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 scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
+    tmp<scalarField> tyPlus = calcYPlus(magUp);
+    scalarField& yPlus = tyPlus();
+
     tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
     scalarField& nutw = tnutw();
 
-    forAll(nutw, facei)
+    forAll(yPlus, facei)
     {
-        scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei];
-
-        scalar yPlus = yPlusLam;
-        scalar ryPlusLam = 1.0/yPlus;
-
-        int iter = 0;
-        scalar yPlusLast = 0.0;
-
-        do
+        if (yPlus[facei] > yPlusLam)
         {
-            yPlusLast = yPlus;
-            yPlus = (kappaRe + yPlus)/(1.0 + log(E_*yPlus));
-
-        } while(mag(ryPlusLam*(yPlus - yPlusLast)) > 0.01 && ++iter < 10 );
-
-        if (yPlus > yPlusLam)
-        {
-            nutw[facei] = nuw[facei]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
+            nutw[facei] =
+                nuw[facei]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0);
         }
     }
 
@@ -148,13 +175,12 @@ nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::calcNut() const
 tmp<scalarField>
 nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::yPlus() const
 {
-    notImplemented
-    (
-        "nutSpalartAllmarasStandardWallFunctionFvPatchScalarField::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);
 }
 
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasStandardWallFunction/nutSpalartAllmarasStandardWallFunctionFvPatchScalarField.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> calcNut() const;
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutSpalartAllmarasWallFunction/nutSpalartAllmarasWallFunctionFvPatchScalarField.C

@@ -167,9 +167,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::calcNut() const
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-
     const scalarField magGradU = mag(Uw.snGrad());
-
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
     return max(0.0, sqr(calcUTau(magGradU))/magGradU - nuw);
@@ -183,9 +181,7 @@ nutSpalartAllmarasWallFunctionFvPatchScalarField::yPlus() const
 
     const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
     const scalarField& y = rasModel.y()[patchI];
-
     const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];
-
     const scalarField& nuw = rasModel.nu().boundaryField()[patchI];
 
     return y*calcUTau(mag(Uw.snGrad()))/nuw;