nutURoughWallFunction: Updated for input consistency with nutkRoughWallFunction

Now both the nutkRoughWallFunction and nutURoughWallFunction us the same field
based Ks and Cs input to support non-uniform surface roughness.

Note that nutURoughWallFunction is not exactly consistent with the
nutUWallFunction in the limit of roughness height = 0 and also not consistent
with nutkRoughWallFunction, particularly if the near-wall cell is in the laminar
sub-layer for which nutURoughWallFunction is currently incorrect.

Significant further work on the nutURoughWallFunction is needed to make it
consistent with both the nutUWallFunction and nutkRoughWallFunction BCs.

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.C

@@ -97,84 +97,78 @@ tmp<scalarField> nutURoughWallFunctionFvPatchScalarField::calcYPlus
     tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
     scalarField& yPlus = tyPlus.ref();
 
-    if (roughnessHeight_ > 0.0)
+    static const scalar c_2 = 2.25/(90 - 2.25);
-    {
+    static const scalar c_3 = 2.0*atan(1.0)/log(90/2.25);
-        // Rough Walls
+    static const scalar c_4 = c_3*log(2.25);
-        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)
+    {
+        if (Ks_[facei] > 0.0)
         {
-            // If KsPlus is based on YPlus the extra term added to the law
+            // Rough Walls
-            // of the wall will depend on yPlus
+
-            forAll(yPlus, facei)
+            const scalar c_1 = 1/(90 - 2.25) + Cs_[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;
+
+            const scalar dKsPlusdYPlus = Ks_[facei]/y[facei];
+
+            do
             {
-                const scalar magUpara = magUp[facei];
+                yPlusLast = yp;
-                const scalar Re = magUpara*y[facei]/nuw[facei];
-                const scalar kappaRe = kappa_*Re;
 
-                scalar yp = yPlusLam_;
+                // The non-dimensional roughness height
-                const scalar ryPlusLam = 1.0/yp;
+                const scalar KsPlus = yp*dKsPlusdYPlus;
 
-                int iter = 0;
+                // The extra term in the law-of-the-wall
-                scalar yPlusLast = 0.0;
+                scalar G = 0.0;
-                scalar dKsPlusdYPlus = roughnessHeight_/y[facei];
 
-                // Additional tuning parameter - nominally = 1
+                scalar yPlusGPrime = 0.0;
-                dKsPlusdYPlus *= roughnessFactor_;
 
-                do
+                if (KsPlus >= 90)
                 {
-                    yPlusLast = yp;
+                    const scalar t_1 = 1 + Cs_[facei]*KsPlus;
+                    G = log(t_1);
+                    yPlusGPrime = Cs_[facei]*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));
+                }
 
-                    // The non-dimensional roughness height
+                scalar denom = 1.0 + log(E_*yp) - G - yPlusGPrime;
-                    scalar KsPlus = yp*dKsPlusdYPlus;
+                if (mag(denom) > vSmall)
+                {
+                    yp = (kappaRe + yp*(1 - yPlusGPrime))/denom;
+                }
+            } while
+            (
+                mag(ryPlusLam*(yp - yPlusLast)) > 0.0001
+             && ++iter < 10
+             && yp > vSmall
+            );
 
-                    // The extra term in the law-of-the-wall
+            yPlus[facei] = max(0.0, yp);
-                    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] = max(0.0, yp);
-            }
         }
-    }
+        else
-    else
-    {
-        // Smooth Walls
-        forAll(yPlus, facei)
         {
+            // Smooth Walls
             const scalar magUpara = magUp[facei];
             const scalar Re = magUpara*y[facei]/nuw[facei];
             const scalar kappaRe = kappa_*Re;
@@ -208,10 +202,22 @@ nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    nutWallFunctionFvPatchScalarField(p, iF),
+    nutUWallFunctionFvPatchScalarField(p, iF),
-    roughnessHeight_(Zero),
+    Ks_(p.size(), 0.0),
-    roughnessConstant_(Zero),
+    Cs_(p.size(), 0.0)
-    roughnessFactor_(Zero)
+{}
+
+
+nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    nutUWallFunctionFvPatchScalarField(p, iF, dict),
+    Ks_("Ks", dict, p.size()),
+    Cs_("Cs", dict, p.size())
 {}
 
 
@@ -223,24 +229,9 @@ nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    nutWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
+    nutUWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    roughnessHeight_(ptf.roughnessHeight_),
+    Ks_(mapper(ptf.Ks_)),
-    roughnessConstant_(ptf.roughnessConstant_),
+    Cs_(mapper(ptf.Cs_))
-    roughnessFactor_(ptf.roughnessFactor_)
-{}
-
-
-nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary& dict
-)
-:
-    nutWallFunctionFvPatchScalarField(p, iF, dict),
-    roughnessHeight_(readScalar(dict.lookup("roughnessHeight"))),
-    roughnessConstant_(readScalar(dict.lookup("roughnessConstant"))),
-    roughnessFactor_(readScalar(dict.lookup("roughnessFactor")))
 {}
 
 
@@ -249,10 +240,9 @@ nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
     const nutURoughWallFunctionFvPatchScalarField& rwfpsf
 )
 :
-    nutWallFunctionFvPatchScalarField(rwfpsf),
+    nutUWallFunctionFvPatchScalarField(rwfpsf),
-    roughnessHeight_(rwfpsf.roughnessHeight_),
+    Ks_(rwfpsf.Ks_),
-    roughnessConstant_(rwfpsf.roughnessConstant_),
+    Cs_(rwfpsf.Cs_)
-    roughnessFactor_(rwfpsf.roughnessFactor_)
 {}
 
 
@@ -262,31 +252,38 @@ nutURoughWallFunctionFvPatchScalarField::nutURoughWallFunctionFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    nutWallFunctionFvPatchScalarField(rwfpsf, iF),
+    nutUWallFunctionFvPatchScalarField(rwfpsf, iF),
-    roughnessHeight_(rwfpsf.roughnessHeight_),
+    Ks_(rwfpsf.Ks_),
-    roughnessConstant_(rwfpsf.roughnessConstant_),
+    Cs_(rwfpsf.Cs_)
-    roughnessFactor_(rwfpsf.roughnessFactor_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-tmp<scalarField> nutURoughWallFunctionFvPatchScalarField::yPlus() const
+void nutURoughWallFunctionFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
 {
-    const label patchi = patch().index();
+    nutUWallFunctionFvPatchScalarField::autoMap(m);
+    m(Ks_, Ks_);
+    m(Cs_, Cs_);
+}
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
-        (
-            turbulenceModel::propertiesName,
-            internalField().group()
-        )
-    );
-    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
-    tmp<scalarField> magUp = mag(Uw.patchInternalField() - Uw);
 
-    return calcYPlus(magUp());
+void nutURoughWallFunctionFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    nutUWallFunctionFvPatchScalarField::rmap(ptf, addr);
+
+    const nutURoughWallFunctionFvPatchScalarField& nrwfpsf =
+        refCast<const nutURoughWallFunctionFvPatchScalarField>(ptf);
+
+    Ks_.rmap(nrwfpsf.Ks_, addr);
+    Cs_.rmap(nrwfpsf.Cs_, addr);
 }
 
 
@@ -294,9 +291,8 @@ void nutURoughWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
     writeLocalEntries(os);
-    writeEntry(os, "roughnessHeight", roughnessHeight_);
+    writeEntry(os, "Cs", Cs_);
-    writeEntry(os, "roughnessConstant", roughnessConstant_);
+    writeEntry(os, "Ks", Ks_);
-    writeEntry(os, "roughnessFactor", roughnessFactor_);
     writeEntry(os, "value", *this);
 }
 

src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutURoughWallFunction/nutURoughWallFunctionFvPatchScalarField.H

@@ -26,30 +26,33 @@ Class
 
 Description
     This boundary condition provides a turbulent kinematic viscosity condition
-    when using wall functions for rough walls, based on velocity.
+    when using wall functions for rough walls, based on velocity.  The
+    condition manipulates the E parameter to account for roughness effects.
 
 Usage
     \table
         Property     | Description             | Required    | Default value
-        roughnessHeight | roughness height     | yes         |
+        Ks           | sand-grain roughness height | yes     |
-        roughnessConstant | roughness constanr | yes         |
+        Cs           | roughness constant      | yes         |
-        roughnessFactor | scaling factor       | yes         |
     \endtable
 
+    Parameter ranges
+    - roughness height (Ks) = sand-grain roughness (0 for smooth walls)
+    - roughness constant (Cs) = 0.5 - 1.0
+
     Example of the boundary condition specification:
     \verbatim
     <patchName>
     {
         type            nutURoughWallFunction;
-        roughnessHeight 1e-5;
+        Ks              uniform 1e-5;
-        roughnessConstant 0.5;
+        Cs              uniform 0.5;
-        roughnessFactor 1;
         value           uniform 0;
     }
     \endverbatim
 
 See also
-    Foam::nutWallFunctionFvPatchScalarField
+    Foam::nutUWallFunctionFvPatchScalarField
 
 SourceFiles
     nutURoughWallFunctionFvPatchScalarField.C
@@ -59,7 +62,7 @@ SourceFiles
 #ifndef nutURoughWallFunctionFvPatchScalarField_H
 #define nutURoughWallFunctionFvPatchScalarField_H
 
-#include "nutWallFunctionFvPatchScalarField.H"
+#include "nutUWallFunctionFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,20 +75,15 @@ namespace Foam
 
 class nutURoughWallFunctionFvPatchScalarField
 :
-    public nutWallFunctionFvPatchScalarField
+    public nutUWallFunctionFvPatchScalarField
 {
     // Private data
 
-        // Roughness model parameters
+        //- Roughness height
+        scalarField Ks_;
 
-            //- Height
+        //- Roughness constant
-            scalar roughnessHeight_;
+        scalarField Cs_;
-
-            //- Constant
-            scalar roughnessConstant_;
-
-            //- Scale factor
-            scalar roughnessFactor_;
 
 
     // Protected Member Functions
@@ -168,53 +166,45 @@ public:
 
     // Member functions
 
-        // Access
+        // Access functions
 
             //- Return the roughness height
-            scalar roughnessHeight() const
+            const scalarField& Ks() const
             {
-                return roughnessHeight_;
+                return Ks_;
             }
 
             //- Return reference to the roughness height to allow adjustment
-            scalar& roughnessHeight()
+            scalarField& Ks()
             {
-                return roughnessHeight_;
+                return Ks_;
             }
 
 
             //- Return the roughness constant scale
-            scalar roughnessConstant() const
+            const scalarField& Cs() const
             {
-                return roughnessConstant_;
+                return Cs_;
             }
 
             //- Return reference to the roughness constant to allow adjustment
-            scalar& roughnessConstant()
+            scalarField& Cs()
             {
-                return roughnessConstant_;
+                return Cs_;
-            }
-
-            //- Return the roughness scale factor
-            scalar roughnessFactor() const
-            {
-                return roughnessFactor_;
-            }
-
-            //- Return reference to the roughness scale factor to allow
-            //  adjustment
-            scalar& roughnessFactor()
-            {
-                return roughnessFactor_;
             }
 
 
-    // I-O
+        // Mapping functions
 
-        // Evaluation functions
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap(const fvPatchFieldMapper&);
 
-            //- Calculate and return the yPlus at the boundary
+            //- Reverse map the given fvPatchField onto this fvPatchField
-            virtual tmp<scalarField> yPlus() const;
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
 
 
         // I-O