Wall function implementation in a FD scheme

salazardetroya · July 20, 2015, 10:53

Hello everybody

I'm trying to implement a k- $\varepsilon$ model using Finite Difference in a simple open channel vertical profile. I'm having issues with the wall function. This is what I'm doing. The equation corresponding to the node near the wall is replaced by

$\frac{U_p}{u_*} = \frac{1}{\kappa} log( \frac{z_p u_*}{\nu} )+ 5.5$

Where

and

are the velocity and distance to the wall of the near wall node. Now as you can see, I have a new unknown, the shear velocity

. What I do is to add another equation to have as many unknowns as equations. I impose the following:

$\kappa z_p \frac{\partial u}{\partial z} = u_*^2$

at the first node near the wall. Because I'm using a FD scheme. to evaluate the gradient at the first node, I use

$\frac{\partial u}{\partial z} = \frac{U_{p+1} - U_p}{\Delta z}$

where $U_{p+1}$ is the node above the first node. This implementation is giving me excessive gradients in the first node. I was wondering if there were a differente implementation for FD schemes.

Thanks
Miguel

sbaffini · July 20, 2015, 11:30

Multiply your first equation by

and subtract

. You will end up with something like

, which needs to be solved iteratively for

(e.g., by Newton-Raphson).

salazardetroya · July 20, 2015, 11:33

Doing so won't I obtain

at the LHS? I forgot to mention that I'm implementing an implicit scheme, therefore

is also unknown.

Hamidzoka · July 20, 2015, 15:25

hi;
First of all, since the application of wall function is intended you should make sure that width of near wall cells are selected such that y+ >11.6 (the exact value can be found in reference books). unless a linear correlation must be used instead of the logarithmic one as you wrote. remember that standard wall function has two parts: first linear and then logarithmic.
Second, even if you have an implicit approach you have always access to a "Up" value from an initial guess or the previous step of iteration. use these values to make your system of equations linear.
Third, wall function approach by nature predicts more or less steep changes in profile compared with low Re models.So, how you judge the calculated data? "excessive gradients in the first node" is just a qualitative phrase and needs to be assessed by some experimental data. Do you have such experimental results?

regards

salazardetroya · July 20, 2015, 15:35

They are excessive in the sense that are much larger than gradients above the first node. The velocity profile is clearly steep in that node, which is what I thought the wall function would avoid. Hence why I thought my error was in trying to impose a gradient in that node.

Given that I have

from the previous iteration, I can use it to calculate the shear velocity with the log law, but what do I use it for? Imposing the gradient as I'm doing?

EDIT: the distance from the wall satisfies the criteria y+ > 11.6

sbaffini · July 21, 2015, 03:52

Quote:

Originally Posted by salazardetroya

Doing so won't I obtain

at the LHS? I forgot to mention that I'm implementing an implicit scheme, therefore

is also unknown.

1) No

2) The typical approach for implicit cases is to discretize implicitly as if no wall function is used and explicitly for the difference between the wall function and the implicit discretization. This is also known as deferred correction.

3) A wall function always return a a wall shear stress (rho*u_tau^2), how to use it is up to you. You should use it as diffusive part on the given wall.

4) As mentioned by others, your function will not work properly if somewhere you get a better resolution (y+ < 11)

salazardetroya · July 21, 2015, 09:44

Thanks for your response, although it's still not clear to me. This is what I obtain doing the algebra operations you mention:

Multiply by

$U_p= \frac{u_*}{\kappa} log( \frac{z_p u_*}{\nu} )+ 5.5 u_*$

Substract

$0 = \frac{u_*}{\kappa} log( \frac{z_p u_*}{\nu} )+ 5.5 u_* - U_p$

should I just take the previous value for

and calculate

?

With respect to including the term in the diffusive part, how could I do so in my equation? I have a vertical profile, with only the streamwise components, my equation is simply something like this:

$\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] + \frac{\rho_s - \rho_w}{\rho_w} c g_x$

In order to evalue the gradient of the diffusive flux, which is:
$\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}$

I should evaluate it using a simple first order FD scheme, first using the value at the node, and at the other side, the shear stress? Something like this?

$\frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] \approx \frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{i} - \rho u_*^2}{\Delta z}$

If this were the case, which $\Delta z$ should I use? The distance to the wall?

sbaffini · July 21, 2015, 11:17

This is a non linearity, just like for nu_t. You either perform outer iterations or simply use previous values (for Up). As i said, you can use an implicit approach with no wall function and add the difference between the wall function and the no wall function case as an explicit source term.

For what concern the implementation:

1) rho is not needed in your last equation

2) I always use central schemes for diffusion (to follow the actual physics)

3) At point 2 (1 being on the wall) i would write the full diffusion term as difference between 2+1/2 and 1+1/2. At 2+1/2 you can evaluate du/dz between 2 and 3 while at 1+1/2 you can use u_tau^2 (being careful to the signs) as the wall function model implicitly assumes that the full diffusion tern is not changing with z.

HOWEVER, i think that using a wall function model in an overall 1D problem might simply be unacceptable. In practice the wall function can be seen as the solution of your 1D problem without the unsteady and gravity terms (and nu_t given by a mixing length model). Consider that there are wall function approaches using your full 1D problem as wall function itself.

Why are you considering wall functions for a 1D problem?

salazardetroya · July 21, 2015, 11:38

I'm implementing a $k-\varepsilon$ turbulence model in a vertical profile. I want to use the wall function to avoid having to refine the mesh near the wall. Near the wall $\nu_t$ is given a mixing length model, which is different than for the rest of the profile ( $\nu_t$ given by the turbulence model) Why could this be wrong?

Edit: I'm writing down the discretization you mentioned so it's clearer to other people and do not confuse it with the one I suggested first (which is wrong). Also I follow this document, page 8 http://www.staff.city.ac.uk/~ra600/M...CFDlecture.pdf

$\left.\frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] \right|_{2}\approx$
$\frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{2+1/2} - \left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{1+1/2}}{\frac{1}{2} (z_{3} - z_{1})} = \frac{\left.\left[\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\right]\right|_{2+1/2} - u^2_*}{\frac{1}{2} (z_{3} - z_{1})} =$
$\frac{\left( \nu + \nu_t \right) \rvert_{2+1/2} \frac{U_{3} - U_{2}}{z_{3} - z_2} - u^2_*}{\frac{1}{2} (z_{3} - z_{1})}$

sbaffini · July 21, 2015, 11:59

This is not wrong per se (also, the k-eps model actually requires the wall function as it cannot be integrated down to the wall).

What is problematic, in my opinion, is that your whole model is a 1D profile. That is, your whole physics is the 1D profile. If you cut your 1D profile to use a wall function (not including the additional terms), i think you are strongly affecting the results.

However, for the k-eps model, the interpretation is clearly different (because you cannot simply increase the resolution) and the wall function can be tought as part of the k-eps model.

Nonetheless, i see some sense in doing this only if this is a sort of exercise. Otherwise i would just use a different turbulent model (k-w) and integrate down to the wall. I don't see the problem in using up to 1000 points in 1D. Your integration should be nearly immediate on any decent processor.

salazardetroya · July 21, 2015, 12:13

I'm assuming my flow is fully developed, therefore I neglect the convective terms in the longitudinal direction. From the continuity equation, I obtain that the vertical velocity component is zero, being left out with just the streamwise component.

What did you mean with "(not including the additional terms),"? You meant the gravity terms?

sbaffini · July 21, 2015, 15:06

Yes, i meant the gravity and unsteady terms.

Consider a full 3D computation; wall functions are used on a minimum part of the domain. Still, they are gratly influencing your results. When i look at your case, i can see that the fraction of the domain which is approximated by a wall function is much higher and i do not expect any better result than in 3D.

Also, i repeat, i do not see any point in using a wall function in a 1D problem, besides pure testing or as an exercise.

However, coming back to your implementation, your last formula is exactly what i meant. Nonetheless, i suggest the following approach to implement it implicitly:

1) Discretize implicitly the velocity equation without concern for the wall function (even for point 2). Which means using:

$\frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{2+1/2} - \left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{1+1/2}}{\frac{1}{2} (z_{3} - z_{1})}$

2) On the right hand side, at point 2, add as source term (evaluated explicitly):

$-\frac{\left.\left[\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\right]\right|_{1+1/2} - u^2_*}{\frac{1}{2} (z_{3} - z_{1})}$

sbaffini · July 21, 2015, 15:13

If you want to add additional terms to the wall function formulation, you can consider this simple model:

http://www.cfd-online.com/Forums/blo...t-effects.html

The Reichardt law is a well known profile covering the whole y+ range starting from 0.
It has been extended (approximately) to include pressure gradient effects at no additional costs. You can then add both the unsteady and the gravity term, assuming them to be constant with z.

salazardetroya · July 21, 2015, 15:32

Quote:

Originally Posted by sbaffini

Also, i repeat, i do not see any point in using a wall function in a 1D problem, besides pure testing or as an exercise.

What do you suggest as an alternative? Implementing a $k-\omega$ model and integrate the entire domain?

With regards to neglecting the unsteady term, the standard wall function assumes this in order to derive the log law, but the flow is usually unsteady, why doesn't the standard wall function include it then? I thought we could assume an steady solution near the wall.

sbaffini · July 22, 2015, 04:09

At this point, i suggest using both approaches (k-eps+WF and k-w) to evaluate the differences, unless you have a very specific reason to use k-eps. I insist on this because,
nowadays, there are codes which use your full 1D equation as wall function (not solving for turbulence, however). Computational resources should not be an issue here.

Standard wall functions are not derived, usually, from an equation (so that you can include terms at your wish). They are just simple (approximate) representations of the near wall velocity profile in very specific conditions (mostly those in steady, fully developed, turbulent channel flows or equilibrium boundary layers). That's just the best you can usually do.

In the link above you can find a trick to include additional terms, as long as they are treated as constant along z. But that's just a trick. And the original Reichardt law is also a trick. I mean, it's just a model, the fact that it is your best option doesn't mean that its absolute merit is high. It depends from the context.

July 20, 2015, 10:53	Wall function implementation in a FD scheme	#1
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	Hello everybody I'm trying to implement a k- $\varepsilon$ model using Finite Difference in a simple open channel vertical profile. I'm having issues with the wall function. This is what I'm doing. The equation corresponding to the node near the wall is replaced by $\frac{U_p}{u_} = \frac{1}{\kappa} log( \frac{z_p u_}{\nu} )+ 5.5$ Where and are the velocity and distance to the wall of the near wall node. Now as you can see, I have a new unknown, the shear velocity . What I do is to add another equation to have as many unknowns as equations. I impose the following: $\kappa z_p \frac{\partial u}{\partial z} = u_*^2$ at the first node near the wall. Because I'm using a FD scheme. to evaluate the gradient at the first node, I use $\frac{\partial u}{\partial z} = \frac{U_{p+1} - U_p}{\Delta z}$ where $U_{p+1}$ is the node above the first node. This implementation is giving me excessive gradients in the first node. I was wondering if there were a differente implementation for FD schemes. Thanks Miguel

Thread Tools	Search this Thread
Show Printable Version Email this Page	Search this Thread: Advanced Search
Display Modes
Linear Mode Switch to Hybrid Mode Switch to Threaded Mode

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[Commercial meshers] Fluent3DMeshToFoam	simvun	OpenFOAM Meshing & Mesh Conversion	50	January 19, 2020 15:33
Natural convection in a closed domain STILL NEEDING help!	Yr0gErG	FLUENT	4	December 2, 2019 00:04
[snappyHexMesh] How to define to right point for locationInMesh	Mirage12	OpenFOAM Meshing & Mesh Conversion	7	March 13, 2016 14:07
Wall Function	aja1345	FLUENT	0	May 27, 2015 11:01
Wall function in adverse pressure gradients	stephane baralon	Main CFD Forum	11	September 2, 1999 04:05

July 20, 2015, 11:30		#2
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	Multiply your first equation by and subtract . You will end up with something like , which needs to be solved iteratively for (e.g., by Newton-Raphson).

July 20, 2015, 11:33		#3
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	Doing so won't I obtain at the LHS? I forgot to mention that I'm implementing an implicit scheme, therefore is also unknown.

July 20, 2015, 15:25		#4
Hamidzoka Senior Member Hamid Zoka Join Date: Nov 2009 Posts: 282 Rep Power: 18	hi; First of all, since the application of wall function is intended you should make sure that width of near wall cells are selected such that y+ >11.6 (the exact value can be found in reference books). unless a linear correlation must be used instead of the logarithmic one as you wrote. remember that standard wall function has two parts: first linear and then logarithmic. Second, even if you have an implicit approach you have always access to a "Up" value from an initial guess or the previous step of iteration. use these values to make your system of equations linear. Third, wall function approach by nature predicts more or less steep changes in profile compared with low Re models.So, how you judge the calculated data? "excessive gradients in the first node" is just a qualitative phrase and needs to be assessed by some experimental data. Do you have such experimental results? regards

July 20, 2015, 15:35		#5
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	They are excessive in the sense that are much larger than gradients above the first node. The velocity profile is clearly steep in that node, which is what I thought the wall function would avoid. Hence why I thought my error was in trying to impose a gradient in that node. Given that I have from the previous iteration, I can use it to calculate the shear velocity with the log law, but what do I use it for? Imposing the gradient as I'm doing? EDIT: the distance from the wall satisfies the criteria y+ > 11.6

July 21, 2015, 09:44		#7
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	Thanks for your response, although it's still not clear to me. This is what I obtain doing the algebra operations you mention: Multiply by $U_p= \frac{u_}{\kappa} log( \frac{z_p u_}{\nu} )+ 5.5 u_$ Substract $0 = \frac{u_}{\kappa} log( \frac{z_p u_}{\nu} )+ 5.5 u_ - U_p$ should I just take the previous value for and calculate ? With respect to including the term in the diffusive part, how could I do so in my equation? I have a vertical profile, with only the streamwise components, my equation is simply something like this: $\frac{\partial u }{\partial t} = \frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] + \frac{\rho_s - \rho_w}{\rho_w} c g_x$ In order to evalue the gradient of the diffusive flux, which is: $\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}$ I should evaluate it using a simple first order FD scheme, first using the value at the node, and at the other side, the shear stress? Something like this? $\frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] \approx \frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{i} - \rho u_*^2}{\Delta z}$ If this were the case, which $\Delta z$ should I use? The distance to the wall?

July 21, 2015, 11:17		#8
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	This is a non linearity, just like for nu_t. You either perform outer iterations or simply use previous values (for Up). As i said, you can use an implicit approach with no wall function and add the difference between the wall function and the no wall function case as an explicit source term. For what concern the implementation: 1) rho is not needed in your last equation 2) I always use central schemes for diffusion (to follow the actual physics) 3) At point 2 (1 being on the wall) i would write the full diffusion term as difference between 2+1/2 and 1+1/2. At 2+1/2 you can evaluate du/dz between 2 and 3 while at 1+1/2 you can use u_tau^2 (being careful to the signs) as the wall function model implicitly assumes that the full diffusion tern is not changing with z. HOWEVER, i think that using a wall function model in an overall 1D problem might simply be unacceptable. In practice the wall function can be seen as the solution of your 1D problem without the unsteady and gravity terms (and nu_t given by a mixing length model). Consider that there are wall function approaches using your full 1D problem as wall function itself. Why are you considering wall functions for a 1D problem?

July 21, 2015, 11:38		#9
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	I'm implementing a $k-\varepsilon$ turbulence model in a vertical profile. I want to use the wall function to avoid having to refine the mesh near the wall. Near the wall $\nu_t$ is given a mixing length model, which is different than for the rest of the profile ( $\nu_t$ given by the turbulence model) Why could this be wrong? Edit: I'm writing down the discretization you mentioned so it's clearer to other people and do not confuse it with the one I suggested first (which is wrong). Also I follow this document, page 8 http://www.staff.city.ac.uk/~ra600/M...CFDlecture.pdf $\left.\frac{\partial}{\partial z} \left[ \left( \nu + \nu_t \right) \frac{\partial u }{\partial z} \right] \right\|_{2}\approx$ $\frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{2+1/2} - \left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{1+1/2}}{\frac{1}{2} (z_{3} - z_{1})} = \frac{\left.\left[\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\right]\right\|_{2+1/2} - u^2_}{\frac{1}{2} (z_{3} - z_{1})} =$ $\frac{\left( \nu + \nu_t \right) \rvert_{2+1/2} \frac{U_{3} - U_{2}}{z_{3} - z_2} - u^2_}{\frac{1}{2} (z_{3} - z_{1})}$

July 21, 2015, 11:59		#10
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	This is not wrong per se (also, the k-eps model actually requires the wall function as it cannot be integrated down to the wall). What is problematic, in my opinion, is that your whole model is a 1D profile. That is, your whole physics is the 1D profile. If you cut your 1D profile to use a wall function (not including the additional terms), i think you are strongly affecting the results. However, for the k-eps model, the interpretation is clearly different (because you cannot simply increase the resolution) and the wall function can be tought as part of the k-eps model. Nonetheless, i see some sense in doing this only if this is a sort of exercise. Otherwise i would just use a different turbulent model (k-w) and integrate down to the wall. I don't see the problem in using up to 1000 points in 1D. Your integration should be nearly immediate on any decent processor.

July 21, 2015, 12:13		#11
salazardetroya New Member Miguel Salazar Join Date: Jul 2015 Posts: 9 Rep Power: 10	I'm assuming my flow is fully developed, therefore I neglect the convective terms in the longitudinal direction. From the continuity equation, I obtain that the vertical velocity component is zero, being left out with just the streamwise component. What did you mean with "(not including the additional terms),"? You meant the gravity terms?

July 21, 2015, 15:06		#12
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	Yes, i meant the gravity and unsteady terms. Consider a full 3D computation; wall functions are used on a minimum part of the domain. Still, they are gratly influencing your results. When i look at your case, i can see that the fraction of the domain which is approximated by a wall function is much higher and i do not expect any better result than in 3D. Also, i repeat, i do not see any point in using a wall function in a 1D problem, besides pure testing or as an exercise. However, coming back to your implementation, your last formula is exactly what i meant. Nonetheless, i suggest the following approach to implement it implicitly: 1) Discretize implicitly the velocity equation without concern for the wall function (even for point 2). Which means using: $\frac{\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{2+1/2} - \left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\rvert_{1+1/2}}{\frac{1}{2} (z_{3} - z_{1})}$ 2) On the right hand side, at point 2, add as source term (evaluated explicitly): $-\frac{\left.\left[\left( \nu + \nu_t \right) \frac{\partial u }{\partial z}\right]\right\|_{1+1/2} - u^2_*}{\frac{1}{2} (z_{3} - z_{1})}$

July 21, 2015, 15:13		#13
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	If you want to add additional terms to the wall function formulation, you can consider this simple model: http://www.cfd-online.com/Forums/blo...t-effects.html The Reichardt law is a well known profile covering the whole y+ range starting from 0. It has been extended (approximately) to include pressure gradient effects at no additional costs. You can then add both the unsteady and the gravity term, assuming them to be constant with z.

July 22, 2015, 04:09		#15
sbaffini Senior Member Paolo Lampitella Join Date: Mar 2009 Location: Italy Posts: 2,152 Blog Entries: 29 Rep Power: 39	At this point, i suggest using both approaches (k-eps+WF and k-w) to evaluate the differences, unless you have a very specific reason to use k-eps. I insist on this because, nowadays, there are codes which use your full 1D equation as wall function (not solving for turbulence, however). Computational resources should not be an issue here. Standard wall functions are not derived, usually, from an equation (so that you can include terms at your wish). They are just simple (approximate) representations of the near wall velocity profile in very specific conditions (mostly those in steady, fully developed, turbulent channel flows or equilibrium boundary layers). That's just the best you can usually do. In the link above you can find a trick to include additional terms, as long as they are treated as constant along z. But that's just a trick. And the original Reichardt law is also a trick. I mean, it's just a model, the fact that it is your best option doesn't mean that its absolute merit is high. It depends from the context.