Wall functions - questions about implementation
 

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.

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.
Q3)
What is db object? It appears from for example here:
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.

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 ;)

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

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?

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

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

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)

Thanks for the comments Anne.


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.

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+.


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

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?

Hey Robert,

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

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.
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.

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

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.

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

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

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 !

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! :-)

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 :D

function definitions
 

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

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

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

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

Hi all,
Could anyone explain where the "k" equation (let's say in the k-epsilon model) does the equivalent to
or (in k-omega)
, 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
The zeroGradient boundary condition only has a member function evaluate
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.

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
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

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?

Yes, most of the wall functions in OpenFOAM, with the exception of the kqRWallFunction, manipulate either the first cell value, or the wall value (in the case of the mut or nut wall functions) and this needs to be explicitly calculated and performed for every cell or cell face.

For more information on wall functions and how they are used you can also take a look at Pope, Turbulent Flows, 2000, pp 442-443 and Versteeg and Malalasekera, An introduction to computational fluid dynamics: The finite volume method, pp. 276-277

su and sp treatment when epsilon gets redined in the first cell
 

Hello

I have a question related to what is being discussed in this thread.
Say for example, I wish to set epsilon in the first cell as a function of k in the first cell. Is simply calling the appropriate wall function sufficient to implement that?

From what I understand,

^calls the relevant epsilon wall function.

Subsequent to this, the solving happens. But what I wonder is, shouldn't the sp and the su terms for the first cell also be modified so that the solver does not overwrite the newly calculated first cell values for epsilon?

For example, if we consider a 1D discretised equation:

ap*eps_p=ae*eps_e+aw*eps_w+su (1)
ap=ae+aw-sp

Now if I set sp and su such as:

su (@first_cell_next_to_wall)=1e10*epsilon_new
sp (@first_cell_next_to_wall)=-1e10

(where epsilon_new is the new value calculated as a function of k)

then,
sp>>aw, sp>>ae. su>>ae, su>>aw. Hence eq. 1 gives

-sp*eps_p=su => eps_p=su/-sp=1e10*epsilon_new/1e10=epsilon_new
Shouldn't such a treatment as above also be there?

Please let me know if I haven't made myself clear. Thanks.

Hi All,

I have a question about the wall functions implemented in Openfoam. Theoretically, when wall function is used, the velocity at the wall is not zero and instead the velocity there is the friction velocity which is solved from the log-law or linear law. (If the zero velocity is emposed there, then it is equivalent to the non-slip boundary condition.) Then the wall shear stress is obtained tau_w=rho*u_tau**2. Then this wall shear stress can be used for the discretization of the momentum equations in the finite volume appraoch. However, in openfoam's specifications, when wall function is applied, the velocity at the wall is still (0 0 0).

Why can we still use (0 0 0) for the wall velocity here?

Thank you in advance.

OFFO

You can use 0 for the velocity, because the shear velocity is not a velocity in this sense. It is calculated as

where is the the wall shear stress calculated from

As you can see from that it is not the velocity at the wall, but the velocity gradient at the wall that plays a role.

I think when you look into the code you can see that you don't need the newly calculated epsilon in the first cell anymore, after you have solved the system. http://foam.sourceforge.net/docs/cpp/a10359_source.html

The only part where you need to use epsilon again is in line 329 and following to calculate mut, however, in the case of wall functions, mut will be set in the first cell anyway and therefore it doesn't matter if your epsilon is wrong or not. So epsilon in the first cell only needs to be correct before we solve epsilon, not after.

Dear Roman,

Thank you for your reply. In the practical CFD implementations, if we use cell centered finite volume method as is done in the OF, the velocity gradient is needed to calculate the numerical flux for the diffusion term in the momentum equations. Is this understanding correct?

Is the following line of reasoning correct:
1, use the log law or linear linear law to get the friction velocity
2, use the equation you mentioned to predict the wall shear stress, which can be used for the momentum equation discretizations.

Is the above understanding correct? Thank you. I am working on the wall function implementations and thus would like to have these knowledge for my work. Thanks.

OFFO

I think your first assumption is correct.

For the other parts:
1. you use the log law or the linear law to find the velocity in the first point away from the wall away from the grid, not he friction velocity. This velocity can then be used to calculate the wall shear stress as you know the velocity at the wall (either no slip, partial slip or whatever you define).

I describe here http://www.cfd-online.com/Forums/ope...-function.html how the nut wall function works, which is the wall function for the velocity. Since setting the velocity in the first point is a bad choice (it will be overwritten), setting a source on the wall strong enough to make the velocity in the first grid point exactly what you need it much better.

Take a look at the references given in the mentioned post and also in this thread. These are:

• Pope, Turbulent Flows, 2000, pp 442-443
• Versteeg and Malalasekera, An introduction to computational fluid dynamics: The finite volume method, pp. 276-277

They are very helpful when it comes to understanding wall function implementations and how wall functions work.

Manan, I thought about the same thing. As far as I understand it, this is done with the "boundaryManipulate" function, called before solving.
I found it in "src/finiteVolume/lnInclude/fvMatrix.C". It calls a "manipulateMatrix"-function for all the boundary patches. This function can be found in the wall function file. Here I am lost, because I don't know which of the functions (from which .C file) "setValues(...)" and "manipulateMatrix" is called. But this pretty much looks like the right spot!
Cheers,
Philipp.

Dear Roman,

Thank you so much for your reply about the wall function implementations. Very helpful for me.

About the following:

My opinion is: the velocity at the first interior grid nodes (cell centroids if cell centered FVM is applied) is from the outer flow calculations like LES and RANS. The height of that nodes is always known. Having these two quantities and the law of wall (like log-law), we can get the friction velocity from iterative solutions of the law of wall.

Please directly point out if what I am saying is not correct.
Thank you.
OFFO

Of course the velocity in the first interior grid point is calculated through the solution of the velocity field. In the law of the wall / or any other wall function, you would like to set the velocity in the first grid point, because you know the velocity in that point based on theory. This is where the extra source term for the calculations comes in, in the case of openfoam, this source term is inserted into the momentum equation through nut/mut at the wall, as mut/nut at the wall are generally zero, due to the fact that k is zero.

I am not sure what you mean with calculating the friction velocity, because it results directly from the solution without having to know the law of the wall. As the friction velocity is

and

where the last part can just be calculated from


As you can see there is not need to calculate the shear velocity using the law of the wall and iterations

Ok, thanks. But in your line of reasoning, what did you mean by setting the velocity in the first grid point? Besides, how was the u_w obtained? I think if we calculate the tau_w directly from the third equations you provided, there will be some inaccuracies. Please comment if any. Thanks.

Wall functions predict the velocity in the first grid point based on theory. However, it is very hard to just prescribe the velocity in the first grid point and keep it fixed. In order to achieve the predicted velocity in the first grid point and not "fix" and thereby destabilize the differential equation (look at Pope, 2000 for this one), we generally introduce a source term, which in openfoam is introduced through the nut/mut wall function.

Thank you for your explanation, Roman.

In Openfoam, when we use wall function for both RANS and LES, the velocity at the walls are zero. So here the non-slip velocity is still applied. Can the slip velocity boundary conditions be used in Openfoam?

best regards,
OFFO

Hej,

when you take a look here http://www.openfoam.org/docs/cpp/ all the way at the bottom, there is a link to boundary conditions. It will take you to the boundary conditions that can be found in OpenFOAM togther (not the wall function, these are separate because they are part of the specific turbulence models). Take a look there. You can find wall boundary conditions and all the explanations how to use them, i.e. how to set them up in the variable files. Under wall boundaries there is something called a partial slip, and there should of course also be a complete slip (otherwise partial can be used as complete slip and no slip based on how large the partial slip is).

Wall Function Implementation in OpenFOAM
 

Hello OpenFOAM community,

I know I am reviving a very old thread, but I thought I'd try this before possibly starting a new thread.

Hey everyone,

Can someone please direct me to some resources which describe the theory behind implementation of this nutUSpaldingWallFunction?

Any assistance is greatly appreciated.

Thank you

Hi Ovi,
did you get any more information on the specifics of this wall function?

Regards,
Martin

Hey Martin,

I just used a combination of looking at the source code, reading the entry in the online OpenFOAM documentation and running appropriate 1D boundary layer profiles using incompressible/boundaryFoam/ to understand how the wall-function which results in a modified friction velocity can be incorporated into the near-wall cells to correct the shear stresses. However, I am not entirely clear about how this affects the eddy-viscosity terms.

Hi Foamers,

I have a question about the standard wall function.
I read that at high Reynolds number in the closest cell to the wall we should calculate the k by this equation:

k= sqr (utau)/Cmu25

but I have seen that everywhere it is recommended to use

kqRWallFunction

that in fact gives the zeroGradient boundary condition.

Maybe my question is completely wrong but could anybody help me to understand this?

Regards,
Amir

I think the above thinking is correct. But I have a question, how to apply the wall shear stress into the momentum equation discretizations as a boundary condition?