Howto IncompressibleWallRoughness

From OpenFOAMWiki
Jump to: navigation, search

1 OBSOLETE OBSOLETE OBSOLETE

2 OBSOLOTE How to add wall roughness OBSOLETE

Im not sure if this is entirely correct, but the test I ran showed improved results on the pressure drop of a simple pipe flow.

The implementation follows the Fluent stuff, found at for instance here http://www.ipt.ntnu.no/manuals/Fluent.Inc/fluent6.3.26/help/html/ug//node253.htm



Go to the directory src/turbulenceModels/RAS/incompressible.

Open the file RASModel/RASModel.H and after the 'dimensionedScalar E_'-line add

       dimensionedScalar kappa_;
       dimensionedScalar E_;

       //- Roughness height
       dimensionedScalar Ks_;
       //- Roughness constant
       dimensionedScalar Cs_;

Next, open up the file RASModel/RASModel.C and after the E_ constructor add

   E_
   (
       dimensioned<scalar>::lookupOrAddToDict
       (
           "E",
           subDict("wallFunctionCoeffs"),
           9.0
       )
   ),

   Ks_
   (
       dimensioned<scalar>::lookupOrAddToDict
       (
           "Ks",
           subDict("wallFunctionCoeffs"),
           0.0
       )
   ),

   Cs_
   (
       dimensioned<scalar>::lookupOrAddToDict
       (
           "Cs",
           subDict("wallFunctionCoeffs"),
           0.0
       )
   ),

Next file to edit: wallFunctions/wallViscosityI.H Change the forAll-loop to what's below.

           forAll(curPatch, facei)
           {
               label faceCelli = curPatch.faceCells()[facei];

               scalar uStar = Cmu25*sqrt(k_[faceCelli]);
               scalar yPlus = y_[patchi][facei]*uStar/nuw[facei];
               scalar KsPlus = Ks_.value()*uStar/nuw[facei];

               // hydrodynamically smooth
               scalar logDB = 1.0;
               
               if (KsPlus > 2.25)
               {
                   //  fully rough
                   if (KsPlus > 90.0)
                   {
                       logDB = 1.0 + Cs_.value()*KsPlus;
                   }
                   else
                   // transitional
                   {
                       if (Ks_.value() > 0.0)
                       {
                           scalar n = ::sin(0.4258*(log(KsPlus) - 0.811));
                           scalar a = (KsPlus - 2.25)/87.75 + Cs_.value()*KsPlus;
                           logDB = ::pow(a, n);
                       }
                   }                    
               }

               if (yPlus > yPlusLam_ || Ks_.value() > 0)
               {
                   nutw[facei] =
                        nuw[facei]
                       *(yPlus*kappa_.value()/log(E_.value()*yPlus/logDB) - 1);
               }
               else
               {
                   nutw[facei] = 0.0;
               }
           }

Now you just need to run 'wmake libso' to recompile the library.

The law of the wall states that

\frac{u_p u^*}{\tau_w/\rho} = \frac{1}{\kappa}\log (E \frac{\rho u^* y_p}{\mu}) + \Delta B

where ΔB is the addition due to surface roughness.

Using the fact that log(a) − log(b) = log(a / b),log(a) * b = log(ab)

we simply just need to modify the constant E by the appropriate factor to incorporate the effect of ΔB

 \frac{1}{\kappa}\log(E Y) - \frac{1}{\kappa}\log(X) = \frac{1}{\kappa}\log(\frac{E}{X}Y)

The only doubt I have is the line

             if (yPlus > yPlusLam_ || Ks_.value() > 0)

which I modified so that if you have surface roughness, it always uses the log-law and it doesnt check if you're in the right region or not.

--Niklas 13:37, 16 December 2008 (CET)