How to implement a wall function?

AliE · May 20, 2018, 10:47

Hello guys,

This question is in my mind since a long time but now its time to solve it. I would like to add k-w models to my fvm colocated code. It is simple to me solve transport equation for k and w but at the wall I have not a clear idea on how to proceed from a practical point of view.
So here are my doubts:
-without wall function k shoud be zero at walls, but omega?
- what have I to do in order to implement a wall function? How boundary coeffs change? In my knowledge wall functions fix the value of k omega and u at the first cell, so I have not clear if I have to remove those unknown from the linear system or not.

Can anybody shed some light on this issue?

Thank,

AliE

sbaffini · May 27, 2018, 10:37

Dear Alie,

wall functions, as many other models, are not implemented in a unique way. Also, some specific models (mostly k-w family) have been developed within companies that, for obvious reasons, do not share any single detail of their implementation.

Now, going to the k-w family, that's probably one of the most obscure one, in the sense that most people tend to implement them in their way.

Let just focus on cell centered finite volume codes.

If no wall function is applied, k=0 is the most prominent correct bc for k. Yet, dk/dn = 0 is also valid as well. For w, which goes to infinity at walls, you can implement it by assigning w in the first few cells with the known analytical behavior near the walls (wilcox method) or, more naturally for a finite volume method, use the analytical w value in the first cell as bc on the cell face (menter method). For the wilcox method, those cells where w is fixed have to be excluded from the system of equations of course.

In case of wall functions, dk/dn=0 is the most common bc for k. Another possibility is to use k = u_tau^2/sqrt(Cmu), where u_tau comes from the velocity wall function (that is already solved at this stage). I have no idea how wall functions for w are implemented for w in case of wilcox method, but for the menter method the idea is to use the same approach as for the no wall function method, just that the analytical wall function is extended to include also the buffer and log layer.

On top of these differences, you also have some more options regarding the analytical wall behaviors for w, that depend from a velocity scale:

1) use u_tau as coming from the velocity wall function
2) use u_tau = u_star = sqrt(k sqrt(Cmu)), where k is in the first cell near the wall
3) a combination of the previous two

For what concerns the implementation in implicit codes, you typically use a deferred correction approach, where wall functions only enter explicitly.

Another topic of concern is if k-w models, like k-e models, require a special treatment of the k production term in the k equation when wall functions are in use. This is typically the case, because the velocity gradient on which the term is based is wrong if calculated naively from the velocity field. But such integration depends from the wall law used for the velocity wall function.

There would be much more to say, but this gives a first overview.

You can find some details in the Fluent manual https://www.sharcnet.ca/Software/Ans...near_wall.html

There are also a couple of reports from chalmers university:

http://www.tfd.chalmers.se/~lada/pos..._report_WF.pdf

http://www.tfd.chalmers.se/~lada/pos...ME_orlando.pdf

Popovac and Hanjalic also have a nice paper:
http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf

This paper also gives some insight on some details I cited above:
https://www.researchgate.net/publica...blicationTitle

Finally, there is a reference must that is from Kalitzin:
http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf

There really is a lot of material on the matter but, let me tell you, you won't find a definitive, satisfactory answer to all of your question. You will need to try several recipes until you find what satisfy the most both your code and your theoretical rigor.

Good Luck

AshwaniAssam · May 28, 2018, 02:02

You can look in this paper for more detail and also in the Menter Automatic wall treatment desciption in Menter's paper or CFX mannual.
http://fluidsengineering.asmedigital...icleid=2670682

AliE · May 28, 2018, 05:27

Dear Paolo,

Thank you for your very nice and correct answer. I am happy to know that I came up to the same conclusions a couple of days ago

Now, I will describe what I have learned, so that other people can take advantage of our discussion and also comment on

As Paolo said, this topic is not easy to cover and a lot of possibilities exists. Here some of my findings described in a "practical way".

Basically you have to solve a transport equation for turbulent kinetic enegy (k) and one for omega. You have a different form of the source term in these equtions as function of the model you choose (Wilcox, BSL, SST ecc...).

Firstly, you have to perform the following distinction:

- low Re (put at least 10 cells in the viscous sublayer)
- high Re (put the first cell in the log-law region (yplus> 30) )

In the first case, the boundary condition are:

- k=0
- omega=10*6*viscos/(Cbeta yplus^2) to be applied at cell center of the boundary cell.

Another possibilities which is very used for k is

- dk/dn=0 ( thus zero gradient) and fix the value of the production Pk (this term appear into the source term of k-equation) at the cell center according to known formulas (see Wilcox's book or Moukalled's book)

For the High Re:

- Solve iteratively, with Newthon-Raphson, the wall function for the velocity, i.e. (U/utau)=(1/kappa)*log(E*yplus), in order to extract utau and later yplus. You have to do this operation in each boundary cell adjacent to a no-slip wall. You can also find after some manipulation a nice formula to compute firstly yplus and then utau. In any case you will end up this step with a value of utau, one value for the wall shear stress (tw) coming from utau and one for yplus to be used later. Mind that U is the magnitude of the velocity parallel to the wall in the cell considered.

- Fix the values of k and omega at the cell center according to the wall function formulation. You can easily do this by setting the diagonal coeff =1 and the neighbour coeffs = 0 into the linear system matrix. The RHS of the corresponding cell is equal to the wanted value for k and omega.

- For the momentum equation: you can exitimate the diffusion flux through the boundary face either explicitly, either implicitly. In the first case just add tw*FaceArea to the RHS in the equation for the boundary cell. In the other case, you can conveniently modify the viscosity at wall (by defining a suitable eddy viscosity again see wilcox or Moukalled) and discretize implicitly as in the standard case.

- You can implement the so called "automatic wall treatment" to prevent the case your yplus falls below 11.46 (viscous sublayer limit) where the wall functions do not perform well. This method is slightly more complicated to implement but surerly more efficient. (I will go for this solution.)

If somebody has something else to say, to add or to correct, you are welcome!

Thak you all,

AliE

AliE · June 1, 2018, 06:08

Hello Guys,

I'm comingup with a simple question. In k-omega family models, which equation is to be solve first?

Is it:

- solve for velocity
- solve for k
- solve for omega

or:

- solve for velocity
- solve for omega
- solve for k

sbaffini · June 1, 2018, 08:23

I don't think there is a correct or wrong answer to this. Also, I think there is only a minor influence. k first seems legit to me.

Still, I solve them coupled, so I can't be of much help here.

I suggest you to take a look into the major open source codes (OpenFOAM, Code_Saturne, SU2) for such minor issues. While they are definitely not the bible (quite the contrary, actually), they can surely give you a hint.

AliE · June 7, 2018, 04:25

Hello Paolo,

My implementation is going well, thank you for your advices. I have another doubt to be solved reguarding the production Pk into the cells adjacent to a wall. I have two options that I find out looking at the literature.

- Peric in his caffa code fixes the production as Pk= abs(tw)*utau/(kappa*dist) where:
- dist is the cell center to wall distance
- utau is calculated as (Cmu^0.25)*sqrt(k)
- tw is the Laminar shear stress (this sounds me a bit strange)

In the another case the same formula for Pk is used but utau and tw come from the wall function updated at the current step (This is my guess).

What do you use in your code?

Thanks,

AliE

AshwaniAssam · June 8, 2018, 01:53

Both way works fine. The later one seemed simpler as, one needs to apply neumann condition for k at the wall with k at the wall calculated from the algebraic equation as seen here,
https://www.cfd-online.com/Wiki/Near...k-omega_models

"..- tw is the Laminar shear stress (this sounds me a bit strange)"
It might be updating the value of twusing wall function but storing it in laminar shear stress term so as keep consistency in the code. (Its just a guess however )

AliE · June 8, 2018, 02:16

Quote:

Originally Posted by AshwaniAssam

Both way works fine. The later one seemed simpler as, one needs to apply neumann condition for k at the wall with k at the wall calculated from the algebraic equation as seen here,
https://www.cfd-online.com/Wiki/Near...k-omega_models

"..- tw is the Laminar shear stress (this sounds me a bit strange)"
It might be updating the value of twusing wall function but storing it in laminar shear stress term so as keep consistency in the code. (Its just a guess however )

Hello and thank you very much for your replay! Well, unfortunately it uses the value of laminar shear stress since tw is calculated as mu*(Up-Uw)/dist. That's why it sounds me strage

however I am convincing myself that this are implementation details that do not have a unique answer, thus thank you for sharing your experience!

sbaffini · June 8, 2018, 06:53

Dear Alie,

first of all, a disclaimer. While I'm currently spending part of my spare time working on my own wall function formulation (still unfinished work, see https://www.cfd-online.com/Forums/bl...on-part-1.html and following blog posts), my work is mostly theoretical, and I rely on previous works for the practical implementation (still, I'm working on that too). Basically, I'm a practitioner, just like you.

Now, the matter is that you can really find two main approaches in the field. The one where $u_{\tau}$ is computed from the

value in the adjacent fluid cell (typically denoted as

in this case) and the one where it comes from the already solved wall function for the velocity field.

Imagine, for example, LES or the Spalart-Allmaras turbulence models. They don't have a

equation to use. In this case, you are obliged to compute $u_{\tau}$ from the velocity wall function (in order to use it for $\tau_w$ at the wall).

Computing $u_{\tau}$ from

was basically introduced with the $k-\epsilon$ family of turbulence models.

The same is also true for the $\tau_w$ actual implementation. You can either store its value directly and later use it, or modify the turbulent viscosity at the walls and compute it as in the laminar case.

As a code developer and mantainer, as well as a theoretically rigorous guy, I consider unacceptable that two turbulence models require different wall function implementations (SA and LES vs. 2 eq. models), at least for what concerns the velocity wall function.

Also, directly storing $\tau_w$ , instead of a turbulent viscosity at walls (that has no sense for laminar cases), is much more simple and less error prone for me.

Thus, this is what I do concerning the momentum equations. This is all you need if using LES (with algebraic SGS models) or the SA model. I also choose to do this independently from the turbulence model, thus even for 2 eq. models.

However, 2 eq. models require the use of wall functions for the turbulence models themselves.

Here it comes another disclaimer. I work on an immersed boundary code. For most things you can imagine it just like an unstructured, cell centered, finite volume code. The only exception is that I need to handle both body fitted and immersed boundaries. This practically means that my

on the bodies can vary greatly, and ensuring it to be within a certain range is not a possibility. As a consequence, I only use 2 eq. models that can work on all

values. This practically means that I only use models of the $k-\omega$ family and the

model of Kalitzin.

Going back to the matter, we already discussed the

equation, which is quite simple to handle with wall functions. I use

without wall functions and $\frac{\partial k}{\partial n}=0$ with wall functions.

What remains is how to set the $\omega$ bc and the production term (if any at all), $P_{k}$ , in the

equation.

First of all, I use $u_{\tau}$ for both, but I am still experimenting on this, including mixing $u_{\tau}$ and

.

Second, for $P_{k}$ I should derive a formula from my specific wall function (see the blog post at the beginning), but I'm not there yet. So, for the moment, I am using the Fluent one for non equilibrium wall functions (eq. 4-321 here https://www.sharcnet.ca/Software/Ans...h_noneqwf.html). In any case, the $\tau_w$ you see in it is the global one coming from the wall function and not the laminar one, that only exists notionally, but not practically, in a turbulent computation (storing $\tau_w$ is of some help here in avoiding confusion).

As I said, my general approach on these things is: do it a la OpenFOAM until a better solution is found. Note that, in this particular case, the OF approach with the $k-\omega$ model is the same as for the $k-\epsilon$ model. Thus, in this very case, I'm not actually following my suggestion.

sbaffini · June 9, 2018, 01:42

Also note that both $\omega$ and

formulas might need to be rearranged properly in order to be used safely in the limit of $u_\tau$ going to 0.

May 20, 2018, 10:47	How to implement a wall function?	#1
AliE Senior Member Join Date: Dec 2017 Posts: 153 Rep Power: 8	Hello guys, This question is in my mind since a long time but now its time to solve it. I would like to add k-w models to my fvm colocated code. It is simple to me solve transport equation for k and w but at the wall I have not a clear idea on how to proceed from a practical point of view. So here are my doubts: -without wall function k shoud be zero at walls, but omega? - what have I to do in order to implement a wall function? How boundary coeffs change? In my knowledge wall functions fix the value of k omega and u at the first cell, so I have not clear if I have to remove those unknown from the linear system or not. Can anybody shed some light on this issue? Thank, AliE kindstrike likes this.

May 27, 2018, 10:37		#2
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	Dear Alie, wall functions, as many other models, are not implemented in a unique way. Also, some specific models (mostly k-w family) have been developed within companies that, for obvious reasons, do not share any single detail of their implementation. Now, going to the k-w family, that's probably one of the most obscure one, in the sense that most people tend to implement them in their way. Let just focus on cell centered finite volume codes. If no wall function is applied, k=0 is the most prominent correct bc for k. Yet, dk/dn = 0 is also valid as well. For w, which goes to infinity at walls, you can implement it by assigning w in the first few cells with the known analytical behavior near the walls (wilcox method) or, more naturally for a finite volume method, use the analytical w value in the first cell as bc on the cell face (menter method). For the wilcox method, those cells where w is fixed have to be excluded from the system of equations of course. In case of wall functions, dk/dn=0 is the most common bc for k. Another possibility is to use k = u_tau^2/sqrt(Cmu), where u_tau comes from the velocity wall function (that is already solved at this stage). I have no idea how wall functions for w are implemented for w in case of wilcox method, but for the menter method the idea is to use the same approach as for the no wall function method, just that the analytical wall function is extended to include also the buffer and log layer. On top of these differences, you also have some more options regarding the analytical wall behaviors for w, that depend from a velocity scale: 1) use u_tau as coming from the velocity wall function 2) use u_tau = u_star = sqrt(k sqrt(Cmu)), where k is in the first cell near the wall 3) a combination of the previous two For what concerns the implementation in implicit codes, you typically use a deferred correction approach, where wall functions only enter explicitly. Another topic of concern is if k-w models, like k-e models, require a special treatment of the k production term in the k equation when wall functions are in use. This is typically the case, because the velocity gradient on which the term is based is wrong if calculated naively from the velocity field. But such integration depends from the wall law used for the velocity wall function. There would be much more to say, but this gives a first overview. You can find some details in the Fluent manual https://www.sharcnet.ca/Software/Ans...near_wall.html There are also a couple of reports from chalmers university: http://www.tfd.chalmers.se/~lada/pos..._report_WF.pdf http://www.tfd.chalmers.se/~lada/pos...ME_orlando.pdf Popovac and Hanjalic also have a nice paper: http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf This paper also gives some insight on some details I cited above: https://www.researchgate.net/publica...blicationTitle Finally, there is a reference must that is from Kalitzin: http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf There really is a lot of material on the matter but, let me tell you, you won't find a definitive, satisfactory answer to all of your question. You will need to try several recipes until you find what satisfy the most both your code and your theoretical rigor. Good Luck arjun and AliE like this.

May 28, 2018, 02:02		#3
AshwaniAssam Senior Member Ashwani Join Date: Sep 2013 Location: Hyderabad Posts: 154 Rep Power: 12	You can look in this paper for more detail and also in the Menter Automatic wall treatment desciption in Menter's paper or CFX mannual. http://fluidsengineering.asmedigital...icleid=2670682 phPatras likes this.

June 8, 2018, 06:53		#10
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	Dear Alie, first of all, a disclaimer. While I'm currently spending part of my spare time working on my own wall function formulation (still unfinished work, see https://www.cfd-online.com/Forums/bl...on-part-1.html and following blog posts), my work is mostly theoretical, and I rely on previous works for the practical implementation (still, I'm working on that too). Basically, I'm a practitioner, just like you. Now, the matter is that you can really find two main approaches in the field. The one where $u_{\tau}$ is computed from the value in the adjacent fluid cell (typically denoted as in this case) and the one where it comes from the already solved wall function for the velocity field. Imagine, for example, LES or the Spalart-Allmaras turbulence models. They don't have a equation to use. In this case, you are obliged to compute $u_{\tau}$ from the velocity wall function (in order to use it for $\tau_w$ at the wall). Computing $u_{\tau}$ from was basically introduced with the $k-\epsilon$ family of turbulence models. The same is also true for the $\tau_w$ actual implementation. You can either store its value directly and later use it, or modify the turbulent viscosity at the walls and compute it as in the laminar case. As a code developer and mantainer, as well as a theoretically rigorous guy, I consider unacceptable that two turbulence models require different wall function implementations (SA and LES vs. 2 eq. models), at least for what concerns the velocity wall function. Also, directly storing $\tau_w$ , instead of a turbulent viscosity at walls (that has no sense for laminar cases), is much more simple and less error prone for me. Thus, this is what I do concerning the momentum equations. This is all you need if using LES (with algebraic SGS models) or the SA model. I also choose to do this independently from the turbulence model, thus even for 2 eq. models. However, 2 eq. models require the use of wall functions for the turbulence models themselves. Here it comes another disclaimer. I work on an immersed boundary code. For most things you can imagine it just like an unstructured, cell centered, finite volume code. The only exception is that I need to handle both body fitted and immersed boundaries. This practically means that my on the bodies can vary greatly, and ensuring it to be within a certain range is not a possibility. As a consequence, I only use 2 eq. models that can work on all values. This practically means that I only use models of the $k-\omega$ family and the model of Kalitzin. Going back to the matter, we already discussed the equation, which is quite simple to handle with wall functions. I use without wall functions and $\frac{\partial k}{\partial n}=0$ with wall functions. What remains is how to set the $\omega$ bc and the production term (if any at all), $P_{k}$ , in the equation. First of all, I use $u_{\tau}$ for both, but I am still experimenting on this, including mixing $u_{\tau}$ and . Second, for $P_{k}$ I should derive a formula from my specific wall function (see the blog post at the beginning), but I'm not there yet. So, for the moment, I am using the Fluent one for non equilibrium wall functions (eq. 4-321 here https://www.sharcnet.ca/Software/Ans...h_noneqwf.html). In any case, the $\tau_w$ you see in it is the global one coming from the wall function and not the laminar one, that only exists notionally, but not practically, in a turbulent computation (storing $\tau_w$ is of some help here in avoiding confusion). As I said, my general approach on these things is: do it a la OpenFOAM until a better solution is found. Note that, in this particular case, the OF approach with the $k-\omega$ model is the same as for the $k-\epsilon$ model. Thus, in this very case, I'm not actually following my suggestion. AliE likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
decomposePar problem: Cell 0contains face labels out of range	vaina74	OpenFOAM Pre-Processing	37	July 20, 2020 05:38
How to implement wall function in my code?	dli	Main CFD Forum	5	May 21, 2015 03:57
[blockMesh] non-orthogonal faces and incorrect orientation?	nennbs	OpenFOAM Meshing & Mesh Conversion	7	April 17, 2013 05:42
[blockMesh] Axisymmetrical mesh	Rasmus Gjesing (Gjesing)	OpenFOAM Meshing & Mesh Conversion	10	April 2, 2007 14:00
Droplet Evaporation	Christian	Main CFD Forum	2	February 27, 2007 06:27

May 28, 2018, 05:27		#4
AliE Senior Member Join Date: Dec 2017 Posts: 153 Rep Power: 8	Dear Paolo, Thank you for your very nice and correct answer. I am happy to know that I came up to the same conclusions a couple of days ago Now, I will describe what I have learned, so that other people can take advantage of our discussion and also comment on As Paolo said, this topic is not easy to cover and a lot of possibilities exists. Here some of my findings described in a "practical way". Basically you have to solve a transport equation for turbulent kinetic enegy (k) and one for omega. You have a different form of the source term in these equtions as function of the model you choose (Wilcox, BSL, SST ecc...). Firstly, you have to perform the following distinction: - low Re (put at least 10 cells in the viscous sublayer) - high Re (put the first cell in the log-law region (yplus> 30) ) In the first case, the boundary condition are: - k=0 - omega=106viscos/(Cbeta yplus^2) to be applied at cell center of the boundary cell. Another possibilities which is very used for k is - dk/dn=0 ( thus zero gradient) and fix the value of the production Pk (this term appear into the source term of k-equation) at the cell center according to known formulas (see Wilcox's book or Moukalled's book) For the High Re: - Solve iteratively, with Newthon-Raphson, the wall function for the velocity, i.e. (U/utau)=(1/kappa)log(Eyplus), in order to extract utau and later yplus. You have to do this operation in each boundary cell adjacent to a no-slip wall. You can also find after some manipulation a nice formula to compute firstly yplus and then utau. In any case you will end up this step with a value of utau, one value for the wall shear stress (tw) coming from utau and one for yplus to be used later. Mind that U is the magnitude of the velocity parallel to the wall in the cell considered. - Fix the values of k and omega at the cell center according to the wall function formulation. You can easily do this by setting the diagonal coeff =1 and the neighbour coeffs = 0 into the linear system matrix. The RHS of the corresponding cell is equal to the wanted value for k and omega. - For the momentum equation: you can exitimate the diffusion flux through the boundary face either explicitly, either implicitly. In the first case just add tw*FaceArea to the RHS in the equation for the boundary cell. In the other case, you can conveniently modify the viscosity at wall (by defining a suitable eddy viscosity again see wilcox or Moukalled) and discretize implicitly as in the standard case. - You can implement the so called "automatic wall treatment" to prevent the case your yplus falls below 11.46 (viscous sublayer limit) where the wall functions do not perform well. This method is slightly more complicated to implement but surerly more efficient. (I will go for this solution.) If somebody has something else to say, to add or to correct, you are welcome! Thak you all, AliE

June 1, 2018, 06:08		#5
AliE Senior Member Join Date: Dec 2017 Posts: 153 Rep Power: 8	Hello Guys, I'm comingup with a simple question. In k-omega family models, which equation is to be solve first? Is it: - solve for velocity - solve for k - solve for omega or: - solve for velocity - solve for omega - solve for k

June 1, 2018, 08:23		#6
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	I don't think there is a correct or wrong answer to this. Also, I think there is only a minor influence. k first seems legit to me. Still, I solve them coupled, so I can't be of much help here. I suggest you to take a look into the major open source codes (OpenFOAM, Code_Saturne, SU2) for such minor issues. While they are definitely not the bible (quite the contrary, actually), they can surely give you a hint.

June 7, 2018, 04:25		#7
AliE Senior Member Join Date: Dec 2017 Posts: 153 Rep Power: 8	Hello Paolo, My implementation is going well, thank you for your advices. I have another doubt to be solved reguarding the production Pk into the cells adjacent to a wall. I have two options that I find out looking at the literature. - Peric in his caffa code fixes the production as Pk= abs(tw)utau/(kappadist) where: - dist is the cell center to wall distance - utau is calculated as (Cmu^0.25)*sqrt(k) - tw is the Laminar shear stress (this sounds me a bit strange) In the another case the same formula for Pk is used but utau and tw come from the wall function updated at the current step (This is my guess). What do you use in your code? Thanks, AliE

June 8, 2018, 01:53		#8
AshwaniAssam Senior Member Ashwani Join Date: Sep 2013 Location: Hyderabad Posts: 154 Rep Power: 12	Both way works fine. The later one seemed simpler as, one needs to apply neumann condition for k at the wall with k at the wall calculated from the algebraic equation as seen here, https://www.cfd-online.com/Wiki/Near...k-omega_models "..- tw is the Laminar shear stress (this sounds me a bit strange)" It might be updating the value of twusing wall function but storing it in laminar shear stress term so as keep consistency in the code. (Its just a guess however )

June 9, 2018, 01:42		#11
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,151 Blog Entries: 29 Rep Power: 39	Also note that both $\omega$ and formulas might need to be rearranged properly in order to be used safely in the limit of $u_\tau$ going to 0.