ENH: actuation disk modifications and rough Wall function for atm inlet

src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.C

@@ -66,6 +66,12 @@ void Foam::actuationDiskSource::checkData() const
            << "disk direction vector is approximately zero"
            << exit(FatalIOError);
     }
+    if (returnReduce(upstreamCellId_, maxOp<label>()) == -1)
+    {
+        FatalErrorIn("Foam::actuationDiskSource::checkData()")
+           << "upstream location " << upstreamPoint_  << " not found in mesh"
+           << exit(FatalIOError);
+    }
 }
 
 
@@ -83,7 +89,9 @@ Foam::actuationDiskSource::actuationDiskSource
     diskDir_(coeffs_.lookup("diskDir")),
     Cp_(readScalar(coeffs_.lookup("Cp"))),
     Ct_(readScalar(coeffs_.lookup("Ct"))),
-    diskArea_(readScalar(coeffs_.lookup("diskArea")))
+    diskArea_(readScalar(coeffs_.lookup("diskArea"))),
+    upstreamPoint_(coeffs_.lookup("upstreamPoint")),
+    upstreamCellId_(-1)
 {
     coeffs_.lookup("fieldNames") >> fieldNames_;
     applied_.setSize(fieldNames_.size(), false);
@@ -91,6 +99,8 @@ Foam::actuationDiskSource::actuationDiskSource
     Info<< "    - creating actuation disk zone: "
         << this->name() << endl;
 
+    upstreamCellId_ = mesh.findCell(upstreamPoint_);
+
     checkData();
 }
 

src/fieldSources/basicSource/actuationDiskSource/actuationDiskSource.H

@@ -48,6 +48,7 @@ Description
             Cp              0.1;        // power coefficient
             Ct              0.5;        // thrust coefficient
             diskArea        5.0;        // disk area
+            upstreamPoint   (0 0 0);    // upstream point
         }
 
 
@@ -92,6 +93,12 @@ protected:
         //- Disk area
         scalar diskArea_;
 
+        //- Upstream point sample
+        point upstreamPoint_;
+
+        //- Upstream cell ID
+        label upstreamCellId_;
+
 
 private:
 

src/fieldSources/basicSource/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -42,20 +42,27 @@ void Foam::actuationDiskSource::addActuationDiskAxialInertialResistance
 ) const
 {
     scalar a = 1.0 - Cp_/Ct_;
-    scalarField T(cells.size());
     vector uniDiskDir = diskDir_/mag(diskDir_);
     tensor E(tensor::zero);
     E.xx() = uniDiskDir.x();
     E.yy() = uniDiskDir.y();
     E.zz() = uniDiskDir.z();
 
-    forAll(cells, i)
+    vector upU = vector(VGREAT, VGREAT, VGREAT);
+    scalar upRho = VGREAT;
+    if (upstreamCellId_ != -1)
     {
-        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
+        upU =  U[upstreamCellId_];
+        upRho = rho[upstreamCellId_];
     }
+    reduce(upU, minOp<vector>());
+    reduce(upRho, minOp<scalar>());
+
+    scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1 - a);
+
     forAll(cells, i)
     {
-        Usource[cells[i]] -= ((Vcells[cells[i]]/V())*T[i]*E) & U[cells[i]];
+        Usource[cells[i]] += ((Vcells[cells[i]]/V())*T*E) & upU;
     }
 }
 

src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSource.H

@@ -53,6 +53,7 @@ Description
             Ct              0.5;        // thrust coefficient
             diskArea        5.0;        // disk area
             coeffs          (0.1 0.5 0.01); // radial distribution coefficients
+            upstreamPoint   (0 0 0);    // upstream point
         }
 
 

src/fieldSources/basicSource/radialActuationDiskSource/radialActuationDiskSourceTemplates.C

@@ -44,11 +44,9 @@ addRadialActuationDiskAxialInertialResistance
 ) const
 {
     scalar a = 1.0 - Cp_/Ct_;
-    scalarField T(cells.size());
     scalarField Tr(cells.size());
     const vector uniDiskDir = diskDir_/mag(diskDir_);
 
-
     tensor E(tensor::zero);
     E.xx() = uniDiskDir.x();
     E.yy() = uniDiskDir.y();
@@ -65,21 +63,27 @@ addRadialActuationDiskAxialInertialResistance
       + radialCoeffs_[1]*sqr(maxR)/2.0
       + radialCoeffs_[2]*pow4(maxR)/3.0;
 
+    vector upU = vector(VGREAT, VGREAT, VGREAT);
+    scalar upRho = VGREAT;
+    if (upstreamCellId_ != -1)
+    {
+        upU =  U[upstreamCellId_];
+        upRho = rho[upstreamCellId_];
+    }
+    reduce(upU, minOp<vector>());
+    reduce(upRho, minOp<scalar>());
+
+    scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1.0 - a);
     forAll(cells, i)
     {
-        T[i] = 2.0*rho[cells[i]]*diskArea_*mag(U[cells[i]])*a/(1.0 - a);
-
         scalar r2 = magSqr(mesh().cellCentres()[cells[i]] - avgCentre);
 
         Tr[i] =
-            T[i]
+            T
            *(radialCoeffs_[0] + radialCoeffs_[1]*r2 + radialCoeffs_[2]*sqr(r2))
            /intCoeffs;
-    }
 
-    forAll(cells, i)
+        Usource[cells[i]] += ((Vcells[cells[i]]/V_)*Tr[i]*E) & upU;
-    {
-        Usource[cells[i]] -= ((Vcells[cells[i]]/V_)*Tr[i]*E) & U[cells[i]];
     }
 
     if (debug)

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.C

@@ -51,6 +51,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     z_(pTraits<vector>::zero),
     z0_(0),
     kappa_(0.41),
+    Uref_(0),
+    Href_(0),
     zGround_(0)
 {}
 
@@ -69,6 +71,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     z_(ptf.z_),
     z0_(ptf.z0_),
     kappa_(ptf.kappa_),
+    Uref_(ptf.Uref_),
+    Href_(ptf.Href_),
     zGround_(ptf.zGround_)
 {}
 
@@ -82,10 +86,12 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
-    Ustar_(readScalar(dict.lookup("Ustar"))),
+    Ustar_(p.size()),
     z_(dict.lookup("z")),
     z0_(readScalar(dict.lookup("z0"))),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
+    Uref_(readScalar(dict.lookup("Uref"))),
+    Href_(readScalar(dict.lookup("Href"))),
     zGround_("zGround", dict, p.size())
 {
     if (mag(z_) < SMALL)
@@ -103,6 +109,11 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
             << abort(FatalError);
     }
 
+    forAll (Ustar_, i)
+    {
+        Ustar_[i] = kappa_*Uref_/(log((Href_  + z0_[i])/max(z0_[i] , 0.001)));
+    }
+
     z_ /= mag(z_);
 
     evaluate();
@@ -121,6 +132,8 @@ atmBoundaryLayerInletEpsilonFvPatchScalarField
     z_(blpsf.z_),
     z0_(blpsf.z0_),
     kappa_(blpsf.kappa_),
+    Uref_(blpsf.Uref_),
+    Href_(blpsf.Href_),
     zGround_(blpsf.zGround_)
 {}
 
@@ -138,14 +151,16 @@ void atmBoundaryLayerInletEpsilonFvPatchScalarField::updateCoeffs()
 void atmBoundaryLayerInletEpsilonFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchScalarField::write(os);
-    os.writeKeyword("Ustar")
-        << Ustar_ << token::END_STATEMENT << nl;
     os.writeKeyword("z")
         << z_ << token::END_STATEMENT << nl;
     os.writeKeyword("z0")
         << z0_ << token::END_STATEMENT << nl;
     os.writeKeyword("kappa")
         << kappa_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Uref")
+        << Uref_ << token::END_STATEMENT << nl;
+    os.writeKeyword("Href")
+        << Href_ << token::END_STATEMENT << nl;
     os.writeKeyword("zGround")
         << zGround_ << token::END_STATEMENT << nl;
     writeEntry("value", os);

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletEpsilon/atmBoundaryLayerInletEpsilonFvPatchScalarField.H

@@ -41,6 +41,16 @@ Description
         z0 is the surface roughness lenght
         zGround minium vlaue in z direction
 
+
+    and:
+
+        Ustar = K Uref/ln((Zref + z0)/z0)
+
+    where:
+
+        Uref is the reference velocity at Zref
+        Zref is the reference height.
+
     \endverbatim
 
     Reference:
@@ -78,17 +88,23 @@ class atmBoundaryLayerInletEpsilonFvPatchScalarField
     // Private data
 
         //- Frictional velocity
-        const scalar Ustar_;
+        scalarField Ustar_;
 
         //- Direction of the z-coordinate
         vector z_;
 
         //- Surface roughness lenght
-        const scalar z0_;
+        scalarField z0_;
 
         //- Von Karman constant
         const scalar kappa_;
 
+        //- Reference velocity
+        const scalar Uref_;
+
+        //- Reference hight
+        const scalar Href_;
+
         //- Minimum corrdinate value in z direction
         const scalarField zGround_;
 
@@ -158,7 +174,7 @@ public:
     // Member functions
 
         //- Return max value
-        scalar Ustar() const
+        scalarField Ustar() const
         {
             return Ustar_;
         }

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.C

@@ -88,16 +88,16 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
 )
 :
     fixedValueFvPatchVectorField(p, iF),
-    Ustar_(0),
+    Ustar_(p.size()),
     n_(dict.lookup("n")),
     z_(dict.lookup("z")),
-    z0_(readScalar(dict.lookup("z0"))),
+    z0_("z0", dict, p.size()),
     kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
     Uref_(readScalar(dict.lookup("Uref"))),
     Href_(readScalar(dict.lookup("Href"))),
     zGround_("zGround", dict, p.size())
 {
-    if (mag(n_) < SMALL || mag(z_) < SMALL || mag(z0_) < SMALL)
+    if (mag(n_) < SMALL || mag(z_) < SMALL)
     {
         FatalErrorIn
         (
@@ -108,14 +108,17 @@ atmBoundaryLayerInletVelocityFvPatchVectorField
                 "onst dictionary&"
             ")"
         )
-            << "magnitude of n, z and z0 vectors must be greater than zero"
+            << "magnitude of n or z must be greater than zero"
             << abort(FatalError);
     }
 
     n_ /= mag(n_);
     z_ /= mag(z_);
 
-    Ustar_ = kappa_*Uref_/(log((Href_  + z0_)/max(z0_ , 0.001)));
+    forAll (Ustar_, i)
+    {
+        Ustar_[i] = kappa_*Uref_/(log((Href_  + z0_[i])/max(z0_[i] , 0.001)));
+    }
 
     evaluate();
 }
@@ -153,12 +156,12 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
         if ((coord[i] - zGround_[i]) < Href_)
         {
             Un[i] =
-                (Ustar_/kappa_)
+                (Ustar_[i]/kappa_)
-              * log((coord[i] - zGround_[i] + z0_)/max(z0_, 0.001));
+              * log((coord[i] - zGround_[i] + z0_[i])/max(z0_[i], 0.001));
         }
         else
         {
-            Un[i] = (Uref_);
+            Un[i] = Uref_;
         }
     }
 
@@ -171,8 +174,7 @@ void atmBoundaryLayerInletVelocityFvPatchVectorField::updateCoeffs()
 void atmBoundaryLayerInletVelocityFvPatchVectorField::write(Ostream& os) const
 {
     fvPatchVectorField::write(os);
-    os.writeKeyword("z0")
+    zGround_.writeEntry("z0", os) ;
-        << z0_ << token::END_STATEMENT << nl;
     os.writeKeyword("n")
         << n_ << token::END_STATEMENT << nl;
     os.writeKeyword("z")

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/atmBoundaryLayerInletVelocity/atmBoundaryLayerInletVelocityFvPatchVectorField.H

@@ -92,7 +92,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
     // Private data
 
         //- Frictional velocity
-        scalar Ustar_;
+        scalarField Ustar_;
 
         //- Flow direction
         vector n_;
@@ -101,7 +101,7 @@ class atmBoundaryLayerInletVelocityFvPatchVectorField
         vector z_;
 
         //- Surface roughness lenght
-        const scalar z0_;
+        scalarField z0_;
 
         //- Von Karman constant
         const scalar kappa_;

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.C

@@ -65,22 +65,10 @@ tmp<scalarField> nutkAtmRoughWallFunctionFvPatchScalarField::calcNut() const
         scalar uStar = Cmu25*sqrt(k[faceCellI]);
         scalar yPlus = uStar*y[faceI]/nuw[faceI];
 
-        scalar Edash = (y[faceI] + z0_[faceI] - zGround_[faceI])/z0_[faceI];
+        scalar Edash = (y[faceI] + z0_[faceI])/z0_[faceI];
 
-        scalar limitingNutw = max(nutw[faceI], nuw[faceI]);
-
-        // To avoid oscillations limit the change in the wall viscosity
-        // which is particularly important if it temporarily becomes zero
         nutw[faceI] =
-            max
+            nuw[faceI]*(yPlus*kappa_/log(max(Edash, 1+1e-4)) - 1);
-            (
-                min
-                (
-                    nuw[faceI]
-                   *(yPlus*kappa_/log(max(Edash, 1+1e-4)) - 1),
-                    2*limitingNutw
-                ), 0.5*limitingNutw
-            );
 
         if (debug)
         {
@@ -105,8 +93,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(p, iF),
-    z0_(p.size(), 0.0),
+    z0_(p.size(), 0.0)
-    zGround_(p.size(), 0.0)
 {}
 
 
@@ -120,8 +107,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(ptf, p, iF, mapper),
-    z0_(ptf.z0_, mapper),
+    z0_(ptf.z0_, mapper)
-    zGround_(ptf.zGround_, mapper)
 {}
 
 
@@ -134,8 +120,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(p, iF, dict),
-    z0_("z0", dict, p.size()),
+    z0_("z0", dict, p.size())
-    zGround_("zGround", dict, p.size())
 {}
 
 
@@ -146,8 +131,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(rwfpsf),
-    z0_(rwfpsf.z0_),
+    z0_(rwfpsf.z0_)
-    zGround_(rwfpsf.zGround_)
 {}
 
 
@@ -159,8 +143,7 @@ nutkAtmRoughWallFunctionFvPatchScalarField
 )
 :
     nutkWallFunctionFvPatchScalarField(rwfpsf, iF),
-    z0_(rwfpsf.z0_),
+    z0_(rwfpsf.z0_)
-    zGround_(rwfpsf.zGround_)
 {}
 
 
@@ -173,7 +156,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::autoMap
 {
     nutkWallFunctionFvPatchScalarField::autoMap(m);
     z0_.autoMap(m);
-    zGround_.autoMap(m);
 }
 
 
@@ -189,7 +171,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::rmap
         refCast<const nutkAtmRoughWallFunctionFvPatchScalarField>(ptf);
 
     z0_.rmap(nrwfpsf.z0_, addr);
-    zGround_.rmap(nrwfpsf.zGround_, addr);
 }
 
 
@@ -198,7 +179,6 @@ void nutkAtmRoughWallFunctionFvPatchScalarField::write(Ostream& os) const
     fvPatchField<scalar>::write(os);
     writeLocalEntries(os);
     z0_.writeEntry("z0", os);
-    zGround_.writeEntry("zGround", os);
     writeEntry("value", os);
 }
 

src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutkAtmRoughWallFunction/nutkAtmRoughWallFunctionFvPatchScalarField.H

@@ -74,9 +74,6 @@ protected:
         //- Surface roughness lenght
         scalarField z0_;
 
-        //- Minimum corrdinate value in z direction
-        scalarField zGround_;
-
 
     // Protected Member Functions