[AlmostSurelyRob]

Dear Foam Users and Developers,

I am really interested in understanding wall functions well. The main reason for this interest is that I would like to know precisely which RAS model is compatible with each of those functions.

Therefore, I embarked on the quest to understand the implantation in OF. To paraphrase one sketch, I would like to emerge one day from this journey, probably a little bruised, but hopefully a little wiser. So here I am after a few days of literature review and perusing the code with a set of compiled questions - not yet wiser, but not too bruised as well. I am hoping to get both with this post.

I am looking at OF 2.0

Q1)
nutUWallFunction stands for Standard Wall function. yplus is estimated with at most 10 Newton iterations under given U and y and then turbulent viscosity is modified.

According to the literature e.g. Wilcox having u_\tau allows to adjust k and epsilon/omega. But in OF I don't see any equivalent wall boundary conditions for k and omega. Instead there is an implementation which changes omega/epsilon and adjusts the production term of k. In fact k and omega/epsilon have only one wall function each. Is nutUWallFunction compatible with epsilon and kqR wall functions?


Q2)
G is k production i.e. \bar{uv} dU/dy. Please look at the comments, the code is taken from omega wall function.

G[faceCellI] =
//This is an estimate of shear stress \bar uv
            (nutw[faceI] + nuw[faceI])
           *magGradUw[faceI]
//This is an estimate of dU/dy from the logarithmic velocity profile
           *Cmu25*sqrt(k[faceCellI])
           /(kappa_*y[faceI]);

Are these comments right and if they are what is the point of having laminar viscosity in shear stress calculation and then calculating the gradient from log-law anyway?

Q2') When is this actually invoked? k equation is always solved second in RAS. omega/epsilon is always solved before and it has this suspiciously looking method boundaryManipulate.

    omegaEqn().boundaryManipulate(omega_.boundaryField());

Q3)
What is db object? It appears from for example here:

    volScalarField& G =
        const_cast<volScalarField&>(db().lookupObject<volScalarField>(GName_));

Q4) Am I right that k wall function in OF is just zeroGradient?

Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re models? Should I put calculated or should I put nutLowReWallFunction or... is there no difference? The latter seems to have only constructors and cloning functions.


I will appreciate any comments or hints on the above questions. I find OF implementation very neat but as I outlined above I have still some doubts and gaps in my knowledge.

[Anne Lincke]

Hey Robert,

I am also investigating in wall functions and how they are implemented in OpemFOAM. Maybe, if we merge our knowledge, we'll become wiser

Originally Posted by AlmostSurelyRob

Q1)
nutUWallFunction stands for Standard Wall function. yplus is estimated with at most 10 Newton iterations under given U and y and then turbulent viscosity is modified.

According to the literature e.g. Wilcox having u_\tau allows to adjust k and epsilon/omega. But in OF I don't see any equivalent wall boundary conditions for k and omega. Instead there is an implementation which changes omega/epsilon and adjusts the production term of k. In fact k and omega/epsilon have only one wall function each. Is nutUWallFunction compatible with epsilon and kqR wall functions?

I do not know nutUWallFunction but I use nutWallFunction.
This wall function is compatible with both, kqRWallFunction and omegaWallFunction/epsilonWallFunction.
It uses the log law to derive the boundary condition for nut, like can be read in this thread

http://www.cfd-online.com/Forums/ope...lfunction.html

Originally Posted by AlmostSurelyRob

Q2)
G is k production i.e. \bar{uv} dU/dy. Please look at the comments, the code is taken from omega wall function.

G[faceCellI] =
//This is an estimate of shear stress \bar uv
            (nutw[faceI] + nuw[faceI])
           *magGradUw[faceI]
//This is an estimate of dU/dy from the logarithmic velocity profile
           *Cmu25*sqrt(k[faceCellI])
           /(kappa_*y[faceI]);

Are these comments right and if they are what is the point of having laminar viscosity in shear stress calculation and then calculating the gradient from log-law anyway?

The comments are right. The production term is \bar uv dU/dy.
\bar uv is equal to (u_tau)^2 = tau_w= (nu + nut) dU/dy = (nutw[faceI] + nuw[faceI])) *magGradUw[faceI]
and dU/dy = u_tau/ (kappa*y) = c_mu^(0.25) * sqrt(k)
But I ask myself, what happens when using a low Reynolds turbulence model?
u_tau is computed according to the log law. This is only valid for a high Reynolds turbulence model using wall functions. Is this ansatz right for a computation without wall functions and y+ ~1?

Originally Posted by AlmostSurelyRob

Q4) Am I right that k wall function in OF is just zeroGradient?

Yes you are right, the k wall function implements a zero Gradient boundary condition.

Originally Posted by AlmostSurelyRob

Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re models? Should I put calculated or should I put nutLowReWallFunction or... is there no difference? The latter seems to have only constructors and cloning functions.

I would use "calculated" for nut when using a low-Re model.


I am concerning about the production term, especially what happens without wall functions for the boundary condition. According to the implementation, G is computed using the log law, no matter how big y+ is.
Isn't this wrong?

Best,
Anne

[romant]

Q5) Finally, what is the correct boundary condition for nut in k-epsilon low-Re models? Should I put calculated or should I put nutLowReWallFunction or... is there no difference? The latter seems to have only constructors and cloning functions.

The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum)

[AlmostSurelyRob]

Thanks for the comments Anne.

Originally Posted by Anne Lincke

I do not know nutUWallFunction but I use nutWallFunction.
This wall function is compatible with both, kqRWallFunction and omegaWallFunction/epsilonWallFunction.

Which version of OF are you using? I am on OF2.1.x and I do not have nutWallFunction at all. Only:

nutUWallFunction
nutkWallFunction

and their rough-wall equivalents.

Have you seen this:
http://http://www.tfd.chalmers.se/~lada/postscript_files/jonas_report_WF.pdf
Please look at Standard wall function (p 11) and Launder-Spalding methodology (incidentally, also p 11. ;-))

This is I believe the difference between nutU* and nutk* wall function. Interestingly, the latter were implemented first and this is why I think omegaWallFunction is and kqRWallFunction were also tailed to this methodology.

Originally Posted by Anne Lincke

But I ask myself, what happens when using a low Reynolds turbulence model?
u_tau is computed according to the log law. This is only valid for a high Reynolds turbulence model using wall functions. Is this ansatz right for a computation without wall functions and y+ ~1?

No, I don't think it's valid. In fact when you look at some tutorials that use low-Re models you will see that that omega/epsilon and k have different BC. BC for k is fixedValue with value=0 and epsilon depending on the model is zero gradient or a value. Wilcox book has a good chapter about this.

Wilcox, Turbulence Modeling for CFD, Chapter 4 section 9 pp 138 Low Reynolds-Number Effects

Have a look if you haven't done so already.



An interesting aspect is that some models employ certain blending procedures. I don't think this is at present implemented in OF.

Two comments here:

1) I once expanded the article in CFD-wiki about this (feel free to correct it you find anything!)
http://www.cfd-online.com/Wiki/Near-...k-omega_models

2) Read the description to kOmegaSST in OF. For some reason OF developers do not believe this sort of tricks.

Nevertheless, I think we could probably implement these low to high switching techniques in OF. Are you interested in doing that? It's just that they might not be valid, but they should at least provide correct limiting behaviours for very low wall y+ and very high y+.

Originally Posted by Anne Lincke

I would use "calculated" for nut when using a low-Re model.
Anne

I was doing the same but romant comment... made me see the light!

Originally Posted by romant

The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum)
Anne

Many thanks romant.

Originally Posted by Anne Lincke

I am concerning about the production term, especially what happens without wall functions for the boundary condition. According to the implementation, G is computed using the log law, no matter how big y+ is.
Isn't this wrong?

I believe you're right here. Shall we run some simple case with y+~1 to see what is the difference between omegaWallFunction with and without this altered production term? That should be fairly straightforward and should settle the matter, shouldn't it?

[Anne Lincke]

Hey Robert,

Originally Posted by AlmostSurelyRob

Which version of OF are you using? I am on OF2.1.x and I do not have nutWallFunction at all. Only:

nutUWallFunction
nutkWallFunction

and their rough-wall equivalents.

I am using OpenFOAM-1.7.0. The nutWallFunction implementd there equals nutkWallFunction.

Originally Posted by AlmostSurelyRob

Have you seen this:
http://http://www.tfd.chalmers.se/~lada/postscript_files/jonas_report_WF.pdf
Please look at Standard wall function (p 11) and Launder-Spalding methodology (incidentally, also p 11. ;-))

This is I believe the difference between nutU* and nutk* wall function. Interestingly, the latter were implemented first and this is why I think omegaWallFunction is and kqRWallFunction were also tailed to this methodology.

I think you are right. The basic difference between nutUWallFunction and nutkWallfunction is the computation of y+. In nutkWallFunction it is computed from the turbulent kinetic energy and the constant c_mu, whereas in nutUWallFunction y+ is simply read in.

For the implementation of the production term have a look at this thread

http://www.cfd-online.com/Forums/ope...tml#post361420

There is an interesting remark, which states that the computation of the production term is correct for low Reynolds number, as it is related to turbulent kinetic energy, which is very small at the viscous sublayer. So the production term also goes to zero in this region.

Originally Posted by AlmostSurelyRob

An interesting aspect is that some models employ certain blending procedures. I don't think this is at present implemented in OF.

Two comments here:

1) I once expanded the article in CFD-wiki about this (feel free to correct it you find anything!)
http://www.cfd-online.com/Wiki/Near-...k-omega_models

2) Read the description to kOmegaSST in OF. For some reason OF developers do not believe this sort of tricks.

There is a blending procedure in OpenFOAM, which is equal to the technique in the thread you posted. It switches between lowRe and highRe, so the omegaWallFunction can be used for both, lowRe and highRe.
Look at the source code of OF-2.0.0.

 scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));

        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);

        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));

This is exactly the same as in the thread which was mentioned before.

[malaboss]

Originally Posted by romant

The nutLowReWallFunction is only their to provide an access to the yPlus calculation tool (which in the original OF implementation gives weird values when using lowRe numbers, I suggest using plusPostRans, which is also available in the forum)

Hi,
I am also wondering what is hidden behind all thos wall conditions, as I want to know which range of y+ is allowed for a given wall condition.

based on what you said romant, the nut file in a case does not define the wall function that will be applied to the case.
But where is it defined if you have for example a k-epsilon turbulence model ?
Indeed kWall and epsilonWall functions are only zeroGradient BC.

Thx for your reply.

[Anne Lincke]

Hey,

the nutWallFunction implements the wall function via the boundary condition for nut.

Look at this thread

http://www.cfd-online.com/Forums/ope...lfunction.html

Kind Regards
Anne

[malaboss]

Thanks a lot, this thread answered clearly my questions and even more !

[malaboss]

Originally Posted by Anne Lincke

Yes you are right, the k wall function implements a zero Gradient boundary condition.

Hi,
I tried to understand why this wall-boundary-condition was used in the k epsilon model, and actually it seems weird to me that the gradient of k near the wall could be zero.

Tell me if I am wrong but on the wall :

k = utau²/sqrt(Cmu)
dk/dn = d[utau²/sqrt(Cmu)]/dn

with utau = sqrt(nu*du/dn) we have

dk/dn = d[ (nu*du/dn)/sqrt(Cmu) ]/dn
dk/dn = [nu/sqrt(Cmu)] * d²u/dn²

if dk/dn =0 then this means that on the wall there is no viscosity constraint which is false.

I must be wrong somewhere but I can't find where...
If you know, let me know !
Thx for all ... to all !

[AlmostSurelyRob]

This near-wall modelling business is so exciting. :-)

My understanding is the following. From the point of view of physics k=0 at the wall. Why? Because by definition

k=1/2 u_i'u_i'

and U=0 and therefore u'=0 at the wall. Meaning no fluctuations at the wall.

But when you do wall functions we approximate the behaviour of the near wall region by using various experimental and theoretical results accounting for turbulence boundary layer structure. You must look up viscous, buffer and logarithmic regions in the literature.

The expression you give
k = utau²/sqrt(Cmu)
has been established for logarithmic subregion only (300>y+>30). This is also called constant shear layer (correct me if I am wrong!). Experiments show that it is an okay approximation for log region.

Note that wall functions impose the value of k in the centre of the cell (rather than the face). dk/dn is just a way of saying - the wall will have no further influence on the value in the cell and will approximate any gradients as if shear in the cell was constant. This is a good for log-region, but if you have a much finer mesh (y+~1) it is simply not true.

Hope this description is a good approximation of what is really happening and hope it helps you on your quest! :-)

[malaboss]

Hi Robert, thanks for the answer.

I agree with you k = 0 at the wall.

For the relation I gave at first,

k = utau²/sqrt(Cmu)

it seems to be only applicable near the wall because in this region the generation of turbulence is equal to its dissipation see (4.2.18 in the link below):

http://staffweb.cms.gre.ac.uk/~ct02/...is/node54.html

If I got it, when we say zeroGradient for k on the wall, it is a way to avoid any more calculation as we already use the wall functions, which is enough ?
But then we should do it for every variable...
Actually, there must be something I still don't understand

[asoltoon]

Hi everybody,
Is there any way to see the functions implementation in Openfoam? I want to see the definition of some turbulence-related functions but don't know how?!

Thanks in advance

[David*]

Hello Ali,
they are hidden at
/opt/OpenFOAM-2.1.0/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions
(your path may differ.)
Regards,
David

[asoltoon]

Vielen Dank David,

But is there any easier way to find the implementation of some functions that are declared in *.C files in the Openfoam website?
For instance if I want to see the implementation of the function DepsilonEff() at line 233 of kEpsilon.C file at link below what should I do?

http://foam.sourceforge.net/docs/cpp/a08551_source.html

[David*]

Just follow the links, I guess:
"Return the effective diffusivity for epsilon.
Definition at line 129 of file kEpsilon.H."
(Note: better to have this kind of conversations in a separate thread.)
Regards, David

[RodriguezFatz]

Hi all,
Could anyone explain where the "k" equation (let's say in the k-epsilon model) does the equivalent to

epsilon_.boundaryField().updateCoeffs();
...
 epsEqn().boundaryManipulate(epsilon_.boundaryField());

or (in k-omega)

omega_.boundaryField().updateCoeffs();
...
 omegaEqn().boundaryManipulate(omega_.boundaryField());

, i.e. where in the code is the selected k-wall function actually called?

[romant]

Originally Posted by RodriguezFatz

Hi all,
Could anyone explain where the "k" equation (let's say in the k-epsilon model) does the equivalent to

epsilon_.boundaryField().updateCoeffs();
...
 epsEqn().boundaryManipulate(epsilon_.boundaryField());

or (in k-omega)

omega_.boundaryField().updateCoeffs();
...
 omegaEqn().boundaryManipulate(omega_.boundaryField());

, i.e. where in the code is the selected k-wall function actually called?

If you take a look at the implementation of the kqRWallFunction, which is the corresponding wall function for k, you will notice, that k does not have a function updateCoeffs. The wall function also describes what it does, which is to set a zeroGradient boundary condition.

I am not completely sure, but I think the zeroBoundary condition is evaluated while solving kEqn in

    kEqn().relax();
    solve(kEqn);
    bound(k_, kMin_);

The zeroGradient boundary condition only has a member function evaluate

template<class Type>
void kqRWallFunctionFvPatchField<Type>::evaluate
(
    const Pstream::commsTypes commsType
)
{
    zeroGradientFvPatchField<Type>::evaluate(commsType);
}

The zeroGradient condition can be applied to many different solution variables, like U, p and T, which all do not use an updateCoeff() function before or after there solution.

[RodriguezFatz]

Ok, for the "k=0" boundary condition, we don't change the equations at all, since only a zero is added.
Now for dk/dn=0, why is the "evaluate" function for k called. I am no professional programmer, but can read the code. Why is the "omega" b.c. explicitly called by the turbulence model, where "k" b.c (you were showing in your post) isn't?

[romant]

Originally Posted by RodriguezFatz

Ok, for the "k=0" boundary condition, we don't change the equations at all, since only a zero is added.
Now for dk/dn=0, why is the "evaluate" function for k called. I am no professional programmer, but can read the code. Why is the "omega" b.c. explicitly called by the turbulence model, where "k" b.c (you were showing in your post) isn't?

The omega or epsilon wall function do not just set a boundary condition for the differential equations in the common sense, like dk/dy=0 or k=const. But they manipulate the value at the first point or on the wall depending on the properties or solution variables close to the wall, which is different at every point of a wall. For example, take a look at omegaWallFunction
The code reads

forAll(mutw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]));
        scalar omegaLog = sqrt(k[faceCellI])/(Cmu25*kappa_*y[faceI]);
        omega[faceCellI] = sqrt(sqr(omegaVis) + sqr(omegaLog));

        G[faceCellI] =
            (mutw[faceI] + muw[faceI])
           *magGradUw[faceI]
           *Cmu25*sqrt(k[faceCellI])
           /(kappa_*y[faceI]);
    }

where the indices faceI and faceCellI refer to the location of the center of the face on the wall and the locate of the center of the cell next to the wall, respectively.

As you can see, not only is the value for omega in the first cell manipulated, but also for G in the first cell, where G is the average turbulence production term sometimes also called

[RodriguezFatz]

So what you are saying is that omega and epsilon wall functions have to change the cell value itself, so thay need to be explitcitly called and "overwrite" the calculated cell value. Whereas the k boundary condition defines the value at the cell's boundary and thus just changes theboundary condition of the wall cell but not the value in the first cell. Thus, k boundary condition changes something that is more hidden in the solving routine.
Is that correct?