
[Sponsors] 
April 6, 2009, 12:11 
Writing Turbulence Solver, ke model

#1 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
All,
I'm writing my own turbulence solver. I've decided to start with the ke 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 

April 7, 2009, 02:58 

#2 
Senior Member
N/A
Join Date: Mar 2009
Posts: 189
Rep Power: 10 
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.


April 7, 2009, 16:03 

#3 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
I'm using a firstorder upwinding scheme for the convection terms. I'm actually using deferred correction so I can blend between firstorder upwinding and central differencing, but right now I'm setting it to use solely firstorder 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 = (1exp(.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. 

April 7, 2009, 17:00 

#4 
Senior Member
N/A
Join Date: Mar 2009
Posts: 189
Rep Power: 10 
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 komega models even when they are integrated straight to the wall. Wilcox has a very good discussion on this issue in chapter 7.


April 7, 2009, 20:20 

#5 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
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 kw model. Does anyone have suggestions for what to use for the ke 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. 

April 7, 2009, 23:56 

#6 
Senior Member
N/A
Join Date: Mar 2009
Posts: 189
Rep Power: 10 
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.


April 13, 2009, 11:14 

#7 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
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. 

October 6, 2009, 05:26 
epsilon boundary condition

#8 
New Member
Rungroj W.
Join Date: Sep 2009
Posts: 4
Rep Power: 10 
Dear sir
Dose anyone have "Lam and Bremhorst, 1981, A modified form of the kepsilon 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 

October 21, 2009, 02:13 

#9 
Senior Member
MAZI
Join Date: Oct 2009
Posts: 100
Rep Power: 10 
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 KE 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..? 

October 21, 2009, 10:53 

#10 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
Mazdak,
Are you using wall functions or walldamping functions? I've only ever written code that will solve the ke model with walldamping functions. To get a converged solution with walldamping 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 

October 21, 2009, 10:57 

#11 
Member
Doug Hunsaker
Join Date: Mar 2009
Location: Logan, UT
Posts: 63
Rep Power: 10 
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 

October 24, 2009, 01:14 

#12 
Senior Member
MAZI
Join Date: Oct 2009
Posts: 100
Rep Power: 10 
Hi doug
can you send me your code? 

October 25, 2009, 03:25 

#13 
New Member
Join Date: Mar 2009
Posts: 7
Rep Power: 10 
can you email me the code as well please?


July 31, 2010, 22:56 

#14 
Senior Member
Cean
Join Date: Feb 2010
Posts: 128
Rep Power: 9 
I am interested. Is it for un structure grid? I Will try to write one.


Tags 
boundary conditions, ke, turbulence model 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Superlinear speedup in OpenFOAM 13  msrinath80  OpenFOAM Running, Solving & CFD  18  March 3, 2015 06:36 
v2f Turbulence Model (and others?)  Travis  FLUENT  4  January 7, 2009 02:18 
baldwin lomax turbulence model  kharati  Main CFD Forum  1  July 29, 2008 02:12 
how could i define a custom turbulence model  FredPacheo  FLUENT  0  July 24, 2008 11:06 