[doug]

All,

I'm writing my own turbulence solver. I've decided to start with the k-e model using the Lam and Bremhorst damping functions as I think this is one of the simplest models.

I'm using deps/dn = 0 as the bc for epsilon, and k = 0 for the bc for k. However, my code is blowing up. I think the problem stems from the fact that f_mu is supposed to be small near the wall, but when you get a small k value near the wall, Re_t gets very small and f_mu gets large. Do I need to constrain f_mu somehow?

Thanks for your help.

Doug

[harishg]

What numerical scheme do you employ for the convection term ? Use of central differences can sometimes cause oscillation and you might need to add a bit of dissipation with upwinding or artificial dissipation. What is your viscosity ? The damping term has been calibrated based on asymptotic expansion near the wall and might not be the real cause of problem. Also deps/dn =0 is not the exact BC. Replace it with e=nu d^2k/dy^2. Check the book of Wilcox for more details.

[doug]

I'm using a first-order upwinding scheme for the convection terms. I'm actually using deferred correction so I can blend between first-order upwinding and central differencing, but right now I'm setting it to use solely first-order upwinding. I'm running at Re=10000 with V = 100 and mu = .01 and rho = 1.0. I've implemented the correct BC for eps like you mentioned and the problem still persists.

Is there some certain grid refinement that I have to have before this model will work? More specifically, I believe the problem originates with f_mu. By definition,

f_mu = (1-exp(-.0165*Ry))^2 * (1 + 20.5/Rt)

As you approach the wall, Rt becomes small because k is small and Rt = k^2/(nu*eps). If you look at the asymptotic behavior of f_mu, you see that it should approach f_mu = (.0165^2)*20.5*2.0 = .01116225 at the wall. This can be derived by looking at the taylor series expansion for k at the wall and realizing that at the wall:
k = 0
dk/dy = 0
eps = nu*d^2k/dy^2

However, numerically, I get huge values of f_mu. Do I need to fix my closest f_mu to the wall to some value? That would seem odd because I thought the purpose of the damping functions was to allow you to integrate straight to the wall.

Any advice would be helpful. Thanks in advance.

[harishg]

You can set the values of k and epsilon at the first grid cell above the wall to the asymptotic value. This has been commonly employed for k-omega models even when they are integrated straight to the wall. Wilcox has a very good discussion on this issue in chapter 7.

[doug]

After looking through Wilcox, it looks like if my cell is within a certain y+ distance from the wall, I'll have to use asymptotic values for a few of the equations. Is this right? I guess I'm surprised because I thought wall functions were supposed to do something like this, and the beauty of damping functions is that you don't have to modify the equations anywhere in the domain... you can just set the bcs and let the thing iterate.

Wilcox gives modifications near a wall for the k-w model. Does anyone have suggestions for what to use for the k-e model? I've searched the orig Lam Bremhorst (1981) paper and the Patel, Rodi, Scheuerer (1985) paper and don't see any suggestions for what to use in this "sublayer" region.

Tks.

[harishg]

This is a numerical problem and it has nothing to do with the with damping function. You can use the modifications of Wilcox. Use the expression for omega to calculate the dissipation. This is obviously not exactly suited for your model but it would work since in the region yplus < 1, the viscous stresses are more important and they determine your dissipation.

[doug]

Thanks for your help! I'm looking at this more closely and will post my results when I get a little further.

In the meantime, does anyone know what Fluent uses very near the wall in this "numerically troublesome" region?

Thanks.

[Sammy]

Dear sir
Dose anyone have "Lam and Bremhorst, 1981, A modified form of the k-epsilon model for predicting wall turbulence? I have a question about epsilon wall boundary condition. They mention with epsilon=nu*d^2k/dy^2 . Some numerical model apply this equation with e(wall) = 4*nu*k1/y1^2 - e1 by obtain subscript "1" denotes the first grid point from the wall. What should be initial value and source term for epsilon boundary condition?

Sam

[mazdak]

Hi Dear Doug.

do you have any benchmarks for comparing results. I write a code in FORTRAN 90 for solving turbulent flow in a 2D duct using Standar K-E model. when the first nodes near the wall (that wall function imposed on it) is in a place that y+ is about 10 my code converges but when is greater than 10 it oscilates.

do you know how much is the value of yplus for first nodes?

greater than 13

about 30 and..?

[doug]

Mazdak,

Are you using wall functions or wall-damping functions? I've only ever written code that will solve the k-e model with wall-damping functions. To get a converged solution with wall-damping functions you have to have y+<1.0 generally. I believe if you're using wall functions your y+ should be between 10 and 60 or so.

Doug

[doug]

Sammy,

The "boundary condition" for epsilon at a wall is simply the second derivative of k at the wall (y+=0). So, depending on how your grid is set up and how you have discretized the domain, you must approximate the second derivative of k at the wall and set the wall value of epsilon to that value. Make sense? The exact implementation of this is dependent on if you're using finite difference, finite volume, or finite element methods.

Doug

[mazdak]

Hi doug

can you send me your code?

[christina]

can you email me the code as well please?

[shirazbj]

I am interested. Is it for un structure grid? I Will try to write one.