CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

implementing a new wall function

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By romant
  • 1 Post By romant

Reply
 
LinkBack Thread Tools Display Modes
Old   August 7, 2013, 10:43
Default implementing a new wall function
  #1
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
I would like to implement a new wall function in OpenFOAM, however, I need to understand the existing wall function implementation first. I would like to know if I understand the current implementation.

at first we have the stresses near the wall

(nu + nut) dUdy_wall = tau_w / rho

where I can now rewrite to

nu_t = tau_w / rho * (dUdy)^-1 - nu= tau_w / rho (U_P / y_P)^-1 - nu

where U_P and y_P are the velocity in the first point from the wall from the field solution and the distance to the first point from the wall, respectively.

I think tau_w can then be replaced with (Versteeg and Malalasekera)

tau_w = rho (c_mu^1/4) (k_P^1/2) U_P kappa / (ln(E yPlus))

or in other version using yPlus=(c_mu^1/4) (k_P^1/2) y_P / nu this becomes

nu_t = nu (yPlus kappa / ln(E yPlus) - 1 )
(OF nutkWallFunction)

As far as I understand I only need to replace the tau_w in the above derivation and have implemented an alternative wall function. I know that I also need to apply different terms to the epsilonwall function and the k wall function if needed. However, I had a harder time understanding the implementation of the momentum term that in Versteeg and Malalasekera is just given as tau_w * A_cell = S and should be somehow introduced in the momentum equation.

The question is now, did I correctly understand the implementation and is it really just necessary to replace tau_w in the above expression in order to get a different wall function and is the replacement of dUdy = U_P/y_P correctly applied.
__________________
~roman
romant is offline   Reply With Quote

Old   April 8, 2014, 11:00
Default
  #2
Senior Member
 
Join Date: Jan 2013
Posts: 297
Rep Power: 6
openfoammaofnepo is on a distinguished road
Dear Roman,

Thank you for your help with my questions in another thread. Your above derivations are correct and helpful. At this stage, actually I have the same question with you:

Code:
How mutw (I am looking at the compressible RAS mutUWallFunction) affects the discretization of the momentum equations?
Do you have further understanding about this issue?

OFFO
openfoammaofnepo is offline   Reply With Quote

Old   April 8, 2014, 12:57
Default
  #3
Senior Member
 
Join Date: Jan 2013
Posts: 297
Rep Power: 6
openfoammaofnepo is on a distinguished road
Dear Roman,

About your last arguments, I think it is correct: the wall shear stress is calculated in different ways in different models. For example, in the following wall function:

Code:
https://github.com/OpenFOAM/OpenFOAM-2.1.x/blob/master/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutUSpaldingWallFunction/mutUSpaldingWallFunctionFvPatchScalarField.C
The wall shear stress is obtained from friction velocity, which is from the Newton iteration from Spalding's Law. About the wall function mutUSpaldingWallFunction in Openfoam, I have the following questions:

1, in the above source files I quoted, what does the following quantity magUp stand for?

Code:
scalarField magUp(mag(Uw.patchInternalField() - Uw));
2, About the velocity gradient normal to the wall, dU/dy, in OF it is equal to:

Code:
const scalarField magGradU(mag(Uw.snGrad()));
So in Openfoam, the dU/dy is set to be the normal gradient at the walls. If I am not making a mistake, this calculation very much depends on the near wall resolutions. If the first node is far from the wall (large yPlus, actually this always the case for most the wall functions in LES), there will be some inaccuracies if we still predict the dU/dy like that.

Any comments are welcome.

OFFO
openfoammaofnepo is offline   Reply With Quote

Old   April 9, 2014, 02:47
Default
  #4
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Quote:
Originally Posted by openfoammaofnepo View Post
Dear Roman,

Thank you for your help with my questions in another thread. Your above derivations are correct and helpful. At this stage, actually I have the same question with you:

Code:
How mutw (I am looking at the compressible RAS mutUWallFunction) affects the discretization of the momentum equations?
Do you have further understanding about this issue?

OFFO
I am not sure what you mean by this. However, when you look at the momentum equation in OpenFOAM, you can see that there is something like (mut+mu)=muEff, which is the effective viscosity. This is where mutw, will be inserted into the equation.
__________________
~roman
romant is offline   Reply With Quote

Old   April 9, 2014, 02:56
Default
  #5
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Quote:
Originally Posted by openfoammaofnepo View Post
Dear Roman,

About your last arguments, I think it is correct: the wall shear stress is calculated in different ways in different models. For example, in the following wall function:

Code:
https://github.com/OpenFOAM/OpenFOAM-2.1.x/blob/master/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/wallFunctions/mutWallFunctions/mutUSpaldingWallFunction/mutUSpaldingWallFunctionFvPatchScalarField.C
The wall shear stress is obtained from friction velocity, which is from the Newton iteration from Spalding's Law. About the wall function mutUSpaldingWallFunction in Openfoam, I have the following questions:

1, in the above source files I quoted, what does the following quantity magUp stand for?

Code:
scalarField magUp(mag(Uw.patchInternalField() - Uw));
2, About the velocity gradient normal to the wall, dU/dy, in OF it is equal to:

Code:
const scalarField magGradU(mag(Uw.snGrad()));
So in Openfoam, the dU/dy is set to be the normal gradient at the walls. If I am not making a mistake, this calculation very much depends on the near wall resolutions. If the first node is far from the wall (large yPlus, actually this always the case for most the wall functions in LES), there will be some inaccuracies if we still predict the dU/dy like that.

Any comments are welcome.

OFFO
  1. magUp stands for the magnitude of the velocity in point P, the first point away from the wall. It uses the relative velocity difference between the velocity in P ( mag(Uw.patchInternalField() ) and the velocity at the wall ( Uw ). This is done in case there is a moving wall, or partial slip at the wall.
  2. This is correct, a larger y+ leads to larger errors in the predictions (but it is not really bad if your y+ is within limits of applicability). I am not sure what the applicability for the Spalding wall function is, as I am usually calculating things in RANS or U-RANS, where the applicability for the law of the wall is about y+=30 to 150, something even a little bit higher. The gradient is actually only used to set the "correction" source to get the right velocity in point P, as you can see in post #1 . It shouldn't have such a large influence. However, as previously said, you always also have to look at the applicability range of the wall function.
__________________
~roman
romant is offline   Reply With Quote

Old   April 9, 2014, 16:46
Question
  #6
Senior Member
 
Join Date: Jan 2013
Posts: 297
Rep Power: 6
openfoammaofnepo is on a distinguished road
Dear Romant,

Thank you for your continuous help. In my understanding, the reason why we need to model the wall stress lies in the following implementations:

if we use cell-centered finite volume discretization just used in OF, we will integrate the NS equation over the all the control volumes. For the diffusion terms (molecular diffusion + turbulent diffusion) in momentum equations, the contribution of each face from diffusive fluxes to each cell (the quantities are stored at the cell cnetroids) is linked to the gradient at the inter-cell face (and of course also the mu and mut at the inter-cell faces). If this face just coincides with the wall, so this gradient is the wall shear stress (or drag).

So the following problem is how to model the wall shear stress. From the above line of reasoning, the two terms (molecular diffusion + turbulent diffusion) will have the numerical flux contributions to the near-wall cell. The following issues are: how to predict the gradient corresponding to the molecular diffusion + turbulent diffusion.

My ultimate question is: why in openfoam, the wall function only applies to the mut, which corresponding to the turbulent diffusion?

Thank you so much. Please correct if what I am saying has some problems.
OFFO


Quote:
Originally Posted by romant View Post
I am not sure what you mean by this. However, when you look at the momentum equation in OpenFOAM, you can see that there is something like (mut+mu)=muEff, which is the effective viscosity. This is where mutw, will be inserted into the equation.
openfoammaofnepo is offline   Reply With Quote

Old   April 9, 2014, 16:49
Default
  #7
Senior Member
 
Join Date: Jan 2013
Posts: 297
Rep Power: 6
openfoammaofnepo is on a distinguished road
For Spalding wall function, they claimed that that equation can be used for the whole boundary layers: viscous sublayer, buffer layer and turbulence regions.

Quote:
Originally Posted by romant View Post
  1. magUp stands for the magnitude of the velocity in point P, the first point away from the wall. It uses the relative velocity difference between the velocity in P ( mag(Uw.patchInternalField() ) and the velocity at the wall ( Uw ). This is done in case there is a moving wall, or partial slip at the wall.
  2. This is correct, a larger y+ leads to larger errors in the predictions (but it is not really bad if your y+ is within limits of applicability). I am not sure what the applicability for the Spalding wall function is, as I am usually calculating things in RANS or U-RANS, where the applicability for the law of the wall is about y+=30 to 150, something even a little bit higher. The gradient is actually only used to set the "correction" source to get the right velocity in point P, as you can see in post #1 . It shouldn't have such a large influence. However, as previously said, you always also have to look at the applicability range of the wall function.
openfoammaofnepo is offline   Reply With Quote

Old   April 10, 2014, 02:54
Default
  #8
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Quote:
Originally Posted by openfoammaofnepo View Post
Dear Romant,

My ultimate question is: why in openfoam, the wall function only applies to the mut, which corresponding to the turbulent diffusion?
OFFO
It is only used for mut, because mut is 0 at the wall, which is why it can be used for something else entirely, like inserting a source term. As far as I can see this is the whole reason behind this. Of course, one could have probably solved in a different way, but in my opinion, this is an elegant solution. It re-purposes something that already exists and does not introduce another term that might not be used in all equations.
For RANS, there are also wall functions for alphat (temperature wall function) and also to set the turbulent production and dissipation (epsilon wall function).

Quote:
Originally Posted by openfoammaofnepo View Post
For Spalding wall function, they claimed that that equation can be used for the whole boundary layers: viscous sublayer, buffer layer and turbulence regions.
Then I think it shouldn't be a problem to use it for the whole range.
JasonWang3 likes this.
__________________
~roman
romant is offline   Reply With Quote

Old   April 10, 2014, 12:11
Question
  #9
Senior Member
 
Join Date: Jan 2013
Posts: 297
Rep Power: 6
openfoammaofnepo is on a distinguished road
Hi Romant,

Thank you for your very helpful suggestions and comments. Did you have any experience in using/implementing LES wall function in OF? How are their performances?

best regards,
OFFO

Quote:
Originally Posted by romant View Post
It is only used for mut, because mut is 0 at the wall, which is why it can be used for something else entirely, like inserting a source term. As far as I can see this is the whole reason behind this. Of course, one could have probably solved in a different way, but in my opinion, this is an elegant solution. It re-purposes something that already exists and does not introduce another term that might not be used in all equations.
For RANS, there are also wall functions for alphat (temperature wall function) and also to set the turbulent production and dissipation (epsilon wall function).



Then I think it shouldn't be a problem to use it for the whole range.
openfoammaofnepo is offline   Reply With Quote

Old   April 11, 2014, 04:54
Default
  #10
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Unfortunately, I don't have any experience with wall functions for LES. I have only worked with wall functions for RANS models and also implemented some for RANS.
__________________
~roman
romant is offline   Reply With Quote

Old   August 18, 2015, 11:51
Default
  #11
New Member
 
Nihar
Join Date: Oct 2014
Posts: 8
Rep Power: 3
Nero_CMU is on a distinguished road
Dear Roman

I too am trying to implement a new wall function in OpenFOAM and have similar doubt as you. Were you able to find if it's necessary to introduce the new tau_w term in the momentum equation for the new wall function to take effect?

Regards
Nihar
Nero_CMU is offline   Reply With Quote

Old   August 19, 2015, 08:49
Default
  #12
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Quote:
Originally Posted by Nero_CMU View Post
Dear Roman

I too am trying to implement a new wall function in OpenFOAM and have similar doubt as you. Were you able to find if it's necessary to introduce the new tau_w term in the momentum equation for the new wall function to take effect?

Regards
Nihar
There is no need to implement the new tau_w into the momentum equation by any means. Boundary conditions mut...WallFunction and nut...WallFunction take care of this. At the wall nut and mut are 0, therefore, you can reuse these values at the wall for other purposes, for example for introducing the tau_w wall shear stress into the momentum equation. So, if you now look at post #1, you see that it is fairly easy to implement a new wall function. Just take any of the simpler wall function models and replace tau_w in that model with the one that you calculate and you obtain a new wall function.
JasonWang3 likes this.
__________________
~roman
romant is offline   Reply With Quote

Old   September 7, 2015, 11:53
Default
  #13
Member
 
Xinguang Wang
Join Date: Feb 2015
Posts: 35
Rep Power: 3
JasonWang3 is on a distinguished road
Quote:
Originally Posted by romant View Post
It is only used for mut, because mut is 0 at the wall, which is why it can be used for something else entirely, like inserting a source term. As far as I can see this is the whole reason behind this. Of course, one could have probably solved in a different way, but in my opinion, this is an elegant solution. It re-purposes something that already exists and does not introduce another term that might not be used in all equations.
For RANS, there are also wall functions for alphat (temperature wall function) and also to set the turbulent production and dissipation (epsilon wall function).



Then I think it shouldn't be a problem to use it for the whole range.
Hi Romant

When I look into the code of kEpsilon.C, I find that only alphat, mut, epsilon have the code of correct the value using wall function, like this:mut_.correctBoundaries(). However, it is hard to find any clue about k, which has different types in the kWallFunction folder. Does it mean that k wall function is not used in solving kEpsilon turbulence modeling?

I know that the epsilon and production term in the epsilon wall function have the influence of solving k equations, but the code in kWallFunction do not include in the program, and in the 0 folder, we have to choose a wall function for k. It looks so strange.

Thanks
JasonWang3 is offline   Reply With Quote

Old   September 8, 2015, 03:49
Default
  #14
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Hi Jason,

the k wall function sets a value at the boundary, which means that the value in the first cell is untouched. It also does not use any values from the first cell to calculate the value k, which by the way is 0. However, nut and epsilon use values from the first cell, and or set values in the first cells, which is why the correction functions from the wallFunction class need to be called. It has to do with the order in which we calculate different values. For example,

1. Update values of epsilon and G in the first cell
2. Solve the epsilon balance equation
3. Solve the k balance equation
4. calculate nut everywhere
5. Correct nut at the wall based on velocity, k and epsilon in the first cell.


Quote:
Originally Posted by JasonWang3 View Post
Hi Romant

When I look into the code of kEpsilon.C, I find that only alphat, mut, epsilon have the code of correct the value using wall function, like this:mut_.correctBoundaries(). However, it is hard to find any clue about k, which has different types in the kWallFunction folder. Does it mean that k wall function is not used in solving kEpsilon turbulence modeling?

I know that the epsilon and production term in the epsilon wall function have the influence of solving k equations, but the code in kWallFunction do not include in the program, and in the 0 folder, we have to choose a wall function for k. It looks so strange.

Thanks
__________________
~roman
romant is offline   Reply With Quote

Old   September 8, 2015, 06:01
Default
  #15
Member
 
Xinguang Wang
Join Date: Feb 2015
Posts: 35
Rep Power: 3
JasonWang3 is on a distinguished road
Quote:
Originally Posted by romant View Post
Hi Jason,

the k wall function sets a value at the boundary, which means that the value in the first cell is untouched. It also does not use any values from the first cell to calculate the value k, which by the way is 0. However, nut and epsilon use values from the first cell, and or set values in the first cells, which is why the correction functions from the wallFunction class need to be called. It has to do with the order in which we calculate different values. For example,

1. Update values of epsilon and G in the first cell
2. Solve the epsilon balance equation
3. Solve the k balance equation
4. calculate nut everywhere
5. Correct nut at the wall based on velocity, k and epsilon in the first cell.
Hi Romant

Thanks for your reply. In the standard wall function, k is assumed to be constant across the near-wall fully turbulence region, and k_v = k_P( the first node, this formula is used in the epsilon wall function). It is sure that k=0 at the wall, however in some more advanced wall functions, k has a distribution, e.g. two-layer model, k=k_v(y/y_v)^2 when y<y_v, and k=by+a when y>y_v.
If I want to add this type of wall function in the code, do I need to add k_correctBoundary in the code, or I just need to change the epsilon wall function, where epsilon and production are modified near the wall, which then have influence to the k-equation?
I think the second way is better, because, as you mentioned above, k=0 at the wall and the distribution of k in the first cell result in the change of epsilon and production term.
Hopefully, you can understand what I said.
Regards,
Jason
JasonWang3 is offline   Reply With Quote

Old   September 8, 2015, 06:42
Default
  #16
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Hi Jason,

I understand what you mean. I haven't looked at it in a while, so I am not sure if you have to run this. I suggest you look into the velocity equations, maybe k boundaries are updated there? Otherwise, you might have to make a small diagram and see where which boundary is updated. It shouldn't take too long.

If this doesn't work out, there is always the option to program it and you add an echo or similar to the location where the update should occur and see if the part of the code is reached at any point.
__________________
~roman

Last edited by romant; September 8, 2015 at 06:43. Reason: spelling
romant is offline   Reply With Quote

Old   September 8, 2015, 10:47
Default
  #17
Member
 
Xinguang Wang
Join Date: Feb 2015
Posts: 35
Rep Power: 3
JasonWang3 is on a distinguished road
Quote:
Originally Posted by romant View Post
Hi Jason,

I understand what you mean. I haven't looked at it in a while, so I am not sure if you have to run this. I suggest you look into the velocity equations, maybe k boundaries are updated there? Otherwise, you might have to make a small diagram and see where which boundary is updated. It shouldn't take too long.

If this doesn't work out, there is always the option to program it and you add an echo or similar to the location where the update should occur and see if the part of the code is reached at any point.
Hi Romant

what do you mean velocity equations? Do you mean momentum equations?

When I look into the code, after solving the momentum equations, there is such code: U.correctBoundaryConditions, then this refer to fixedValueFvPatchField at the wall. However, in this code, there is no clue of k.

For my understanding, if it refer to a wall boundary condition, and wall function is used, the first cell of velocity should be modified. But I can not find anything in the code about this.

Thanks a lot for your suggestion.

Jason
JasonWang3 is offline   Reply With Quote

Old   September 8, 2015, 15:21
Default
  #18
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 366
Rep Power: 13
romant is on a distinguished road
Hej Jason,

the velocity is not modified in the first cell, because this might introduce instabilities in the equations. instead we use mut or nut to introduce a shear force at the wall which is large enough to created the value of U which we want to introduce in the first cell. please take a look at http://www.cfd-online.com/Forums/ope...tml#post452744 for some information on this. The whole thread could be of interest to you.
__________________
~roman
romant is offline   Reply With Quote

Old   September 9, 2015, 07:22
Default
  #19
Member
 
Xinguang Wang
Join Date: Feb 2015
Posts: 35
Rep Power: 3
JasonWang3 is on a distinguished road
Quote:
Originally Posted by romant View Post
Hej Jason,

the velocity is not modified in the first cell, because this might introduce instabilities in the equations. instead we use mut or nut to introduce a shear force at the wall which is large enough to created the value of U which we want to introduce in the first cell. please take a look at http://www.cfd-online.com/Forums/ope...tml#post452744 for some information on this. The whole thread could be of interest to you.
Hi Romant

Thank you for your reply. I did some numerical experiment, in which only the k BC is changed, the numerical results shows that the choice of the k wall function do influence the results.

Jason
JasonWang3 is offline   Reply With Quote

Old   September 9, 2015, 22:34
Default
  #20
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 277
Rep Power: 5
huangxianbei is on a distinguished road
Hi,Roman:
Could you please take a look at my wall function? I implemented a wall function to modify the velocity near-wall. However, no change is seen during the calculation, I don't know why. Could you give me some advice?

Xianbei
huangxianbei is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Radiation interface hinca CFX 15 January 26, 2014 18:11
Compile problem ivanyao OpenFOAM Running, Solving & CFD 1 October 12, 2012 09:31
channelFoam for a 3D pipe AlmostSurelyRob OpenFOAM 3 June 24, 2011 13:06
Axisymmetrical mesh Rasmus Gjesing (Gjesing) OpenFOAM Native Meshers: blockMesh 10 April 2, 2007 14:00
Quick Question - Wall Function D.Tandra Main CFD Forum 2 March 16, 2004 05:29


All times are GMT -4. The time now is 08:32.