CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Wall function implementation in a FD scheme (https://www.cfd-online.com/Forums/main/157111-wall-function-implementation-fd-scheme.html)

salazardetroya July 20, 2015 10:53

Wall function implementation in a FD scheme
 
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 U_p and z_p 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 u_*. 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 u_* and subtract U_p. You will end up with something like f(u_*)=0, which needs to be solved iteratively for u_* (e.g., by Newton-Raphson).

salazardetroya July 20, 2015 11:33

Doing so won't I obtain U_p at the LHS? I forgot to mention that I'm implementing an implicit scheme, therefore U_p 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 U_p 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 (Post 556188)
Doing so won't I obtain U_p at the LHS? I forgot to mention that I'm implementing an implicit scheme, therefore U_p 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_*
U_p= \frac{u_*}{\kappa} log( \frac{z_p u_*}{\nu} )+ 5.5 u_*

Substract U_p

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

should I just take the previous value for U_p and calculate u_* ?

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 (Post 556373)
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.


All times are GMT -4. The time now is 19:52.