[sinatahmooresi]

Hello dear Foammers,



I have difficulties finding out about the yPlus definition in OF. Does openFoam calculate the y in y+ formula from the first node after above the wall boundary or the cell centre (half of the distance from the wall to the first node)? I could see from the source code (yPlus.C):

yPlusBf[patchi] =
d[patchi]

*sqrt

(

model.nuEff(patchi)

*mag(UBf[patchi].snGrad())

)/model.nu(patchi);


where, d in d[patchi]:
const nearWallDist nwd(mesh_);
const volScalarField::Boundary& d = nwd.y();


where nearWallDist returns: Distance calculation for cells with face on a wall. Searches pointNeighbours to find closest (nearWallDist.H).





Based on the above mentioned, my take is on is that OF accounts for the first node distance (not the cell centre), which is more reasonable to me given that the cell centres (volume cells) are for scalars while velocity cell centre is on scalar face.

Thank you,

Sina

[Tobermory]

Your terminology is confusing and doesn't match OpenFOAM (are you used to using CFX, by any chance?). The only "nodes" in OpenFOAM are at the cell centre ... this is where all of the primitive variables, including U, are stored/defined. This is a necessary condition for a finite volume, non-staggered grid solver like OpenFOAM. Fluxes, such as the mass flux phi, are calculated at the centres of the cell surfaces. No data is calculated at the cell vertices.

Coming back to your question, the wall distance calculation is the distance from the wall to the cell centre, which is the only metric that makes sense in OF. Note that for highly skewed cells, this distance varies depending on the wallDist method that you have specified (i.e. whether you choose a method that calculates the wall normal component or the absolute distance). Hence, the yPlus value is based on the "half-height" of the cell.

[sinatahmooresi]

Originally Posted by Tobermory

Your terminology is confusing and doesn't match OpenFOAM (are you used to using CFX, by any chance?). The only "nodes" in OpenFOAM are at the cell centre ... this is where all of the primitive variables, including U, are stored/defined. This is a necessary condition for a finite volume, non-staggered grid solver like OpenFOAM. Fluxes, such as the mass flux phi, are calculated at the centres of the cell surfaces. No data is calculated at the cell vertices.

Coming back to your question, the wall distance calculation is the distance from the wall to the cell centre, which is the only metric that makes sense in OF. Note that for highly skewed cells, this distance varies depending on the wallDist method that you have specified (i.e. whether you choose a method that calculates the wall normal component or the absolute distance). Hence, the yPlus value is based on the "half-height" of the cell.

Thank you for your reply. I was confused with the staggered and non-staggered approaches and the one hired in OF. Thank you for the clarification. I was wrong using the terminology "nodes" when trying to convey my understanding. I know that finite volume, in a conventional staggered approach, has the pressure and scalars on the cell centres while velocity (which is expressed in terms of u*A; flux) is defined on centres of the cell surface [please correct if I am wrong]

My confusion with yplus, apart from the wrong way of expression in which I used "nodes" in a wrong way, comes from my calculation to see whether I can get to the "u_tau" from two different way.

OF returns both wallShearStress and yPlus. Regrading the definitions of these two parameters, u_tau can be derived either from the wallShearStress (as sqrt(wallShearStress/rho), where rho would be 1 in an incompressible solver) or from the yPlus as yplus*nu/y, where y would be the distance of the wall to the first cell centre. However, using these two methods, I can't verify one from another, and it puzzles me with the calculation procedure.


Back to the "node" term I used, can you please clarify on the "nearWallDist" in OF, as it "Searches pointNeighbours to find closest" as stated by the OF documentation. where the half-height you mentioned is indicated ?
Thank you for time and valuable reply

[Tobermory]

No problem - just remember that OpenFOAM is based on a non-staggered grid approach.

Now, as to why you cannot reproduce yplus using the wall shear stress - this is interesting and a little complex! Essentially, it boils down to the wall function. In the earliest CFD codes, the wall function was simply based on the friction velocity, Utau. This is fine for flows with no pressure gradient, but consider what happens to U+ at a separation point where the wall shear stress (and therefore Utau) goes to zero. Hence, modern practice is to use a more complex wall function approach, and to define a wall scaling velocity based on the TKE level in the wall adjacent cell, U* = Cmu^0.25 k^0.5. With this definition, we define the normalised wall distance as y*, and the wall function is more readily applicable to a wider range of flows.

So what does this mean for the yPlus postprocessing function? Well, the value returned depends on the boundary patch type. Let's look at the coding. In yPlus.C we have (I have stripped out most of the lines, for the sake of brevity):

 Foam::tmp<Foam::volScalarField> Foam::functionObjects::yPlus::calcYPlus
 (
     const momentumTransportModel& turbModel
 )
 {
    ...
 
     forAll(patches, patchi)
     {
         ...
 
         if (isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi]))
         {
             const nutWallFunctionFvPatchScalarField& nutPf =
                 dynamic_cast<const nutWallFunctionFvPatchScalarField&>
                 (
                     nutBf[patchi]
                 );
 
             yPlusBf[patchi] = nutPf.yPlus();
         }
         else if (isA<wallFvPatch>(patch))
         {
             yPlusBf[patchi] =
                 d[patchi]
                *sqrt
                 (
                     nuEffBf[patchi]
                    *mag(turbModel.U().boundaryField()[patchi].snGrad())
                 )/nuBf[patchi];
         }
     }
 
     return tyPlus;
 }

So, if the patch is a wall function patch (first part of the if statement), then it sets yPlus equal to the return from the wall function yPlus() function. e.g. if we assume that that's nutkWallFunctionFvPatchScalarField.C, then we have:

 tmp<scalarField> nutkWallFunctionFvPatchScalarField::yPlus() const
 {
     ...
     return pow025(Cmu_)*y*sqrt(kwc)/nuw;
 }

which is the y* calculation that I listed above.

If there's no wall function applied on the wall boundary (second part of the if statement), then it will calculate the wall shear stress from the velocity gradient and effective viscosity, and will then use that to calculate y+, and returns this for yPlus.

So the yPlus utility returns either y* or y+, depending on the patch type, and if you want to check that it is returning the right value then for your application you probably need to be comparing against a calculated y* and not a calculated y+. As a demo, here are some numbers for two wall cells from the simpleFoam pitzDaily tutorial:

variable face 1 face 3
y 0.00015942 0.00015942
nu 1.00E-05 1.00E-05

tau/rho 0.613171 0.594907
uTau 0.783052361 0.771302146

kCell 0.362828 0.327753
u* 0.329921809 0.313569609

y+ 12.48342 12.29610
y* 5.25961 4.99893
yPlus 5.25962 4.99894

Back to the "node" term I used, can you please clarify on the "nearWallDist" in OF, as it "Searches pointNeighbours to find closest" as stated by the OF documentation. where the half-height you mentioned is indicated ?
Thank you for time and valuable reply

Agreed that it is perhaps misleading, but it is inferred that the distance is to the cell centre, since that is where all of the primary variable data is stored. You can dive into the coding if you wish.

[sinatahmooresi]

Originally Posted by Tobermory

No problem - just remember that OpenFOAM is based on a non-staggered grid approach.

Now, as to why you cannot reproduce yplus using the wall shear stress - this is interesting and a little complex! Essentially, it boils down to the wall function. In the earliest CFD codes, the wall function was simply based on the friction velocity, Utau. This is fine for flows with no pressure gradient, but consider what happens to U+ at a separation point where the wall shear stress (and therefore Utau) goes to zero. Hence, modern practice is to use a more complex wall function approach, and to define a wall scaling velocity based on the TKE level in the wall adjacent cell, U* = Cmu^0.25 k^0.5. With this definition, we define the normalised wall distance as y*, and the wall function is more readily applicable to a wider range of flows.

So what does this mean for the yPlus postprocessing function? Well, the value returned depends on the boundary patch type. Let's look at the coding. In yPlus.C we have (I have stripped out most of the lines, for the sake of brevity):

 Foam::tmp<Foam::volScalarField> Foam::functionObjects::yPlus::calcYPlus
 (
     const momentumTransportModel& turbModel
 )
 {
    ...
 
     forAll(patches, patchi)
     {
         ...
 
         if (isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi]))
         {
             const nutWallFunctionFvPatchScalarField& nutPf =
                 dynamic_cast<const nutWallFunctionFvPatchScalarField&>
                 (
                     nutBf[patchi]
                 );
 
             yPlusBf[patchi] = nutPf.yPlus();
         }
         else if (isA<wallFvPatch>(patch))
         {
             yPlusBf[patchi] =
                 d[patchi]
                *sqrt
                 (
                     nuEffBf[patchi]
                    *mag(turbModel.U().boundaryField()[patchi].snGrad())
                 )/nuBf[patchi];
         }
     }
 
     return tyPlus;
 }

So, if the patch is a wall function patch (first part of the if statement), then it sets yPlus equal to the return from the wall function yPlus() function. e.g. if we assume that that's nutkWallFunctionFvPatchScalarField.C, then we have:

 tmp<scalarField> nutkWallFunctionFvPatchScalarField::yPlus() const
 {
     ...
     return pow025(Cmu_)*y*sqrt(kwc)/nuw;
 }

which is the y* calculation that I listed above.

If there's no wall function applied on the wall boundary (second part of the if statement), then it will calculate the wall shear stress from the velocity gradient and effective viscosity, and will then use that to calculate y+, and returns this for yPlus.

So the yPlus utility returns either y* or y+, depending on the patch type, and if you want to check that it is returning the right value then for your application you probably need to be comparing against a calculated y* and not a calculated y+. As a demo, here are some numbers for two wall cells from the simpleFoam pitzDaily tutorial:

variable face 1 face 3
y 0.00015942 0.00015942
nu 1.00E-05 1.00E-05

tau/rho 0.613171 0.594907
uTau 0.783052361 0.771302146

kCell 0.362828 0.327753
u* 0.329921809 0.313569609

y+ 12.48342 12.29610
y* 5.25961 4.99893
yPlus 5.25962 4.99894




Agreed that it is perhaps misleading, but it is inferred that the distance is to the cell centre, since that is where all of the primary variable data is stored. You can dive into the coding if you wish.

Hello Tobermory,

Thank you for this clarifying explanation.
Speaking of wall function , I am using WALE model in OF for sgs model, and as WALE, it accounts for an adaptive nu_sgs at the vicinity of the wall boundary. Based on what I can see in OF documentation, it looks like that "nutLowReWallFunction" would serve fine when having a wall resolved approach by LES in OF. It might seem controversial, but using this wall function, OF puts the nut = 0 at the wall and calculate the y+. It is another source of confusion in cfd codes for me. I know that OF calculates first "yplus_laminar", which is a decisive parameter to see whether it switch on which formula for the wall cell treatment. As far as I know, the production of TKE is zero at the wall not the TKE itself. Therefore, putting the nu_sgs to zero directly at the wall boundary in a wall resolved LES, does not sounds correct to me. The latter, is a practice I've seen it in some threads of conversation on the wall treatment in OF. At the end of the day, what is the fate of the nu_sgs on the wall?



Could you share thoughts with me on this matter as well.



Regards,

Sina

[Tobermory]

Within the viscous sublayer, i.e. y+ < 5, viscosity dominates and there is negligible turbulence (i.e. TKE ~ 0). So putting nut or nu_sgs equal to zero at the wall seems like a good idea to me, in a boundary layer resolved simulation.