# Writing Turbulence Solver, k-e model

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

 April 6, 2009, 12:11 Writing Turbulence Solver, k-e model #1 Member   Doug Hunsaker Join Date: Mar 2009 Location: Logan, UT Posts: 63 Rep Power: 8 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

 April 7, 2009, 02:58 #2 Senior Member   N/A Join Date: Mar 2009 Posts: 188 Rep Power: 8 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: 8 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.

 April 7, 2009, 17:00 #4 Senior Member   N/A Join Date: Mar 2009 Posts: 188 Rep Power: 8 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.

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

 April 7, 2009, 23:56 #6 Senior Member   N/A Join Date: Mar 2009 Posts: 188 Rep Power: 8 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: 8 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: 8 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

 October 21, 2009, 02:13 #9 Member   MAZI Join Date: Oct 2009 Posts: 80 Rep Power: 8 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..?

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

 October 21, 2009, 10:57 #11 Member   Doug Hunsaker Join Date: Mar 2009 Location: Logan, UT Posts: 63 Rep Power: 8 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 Member   MAZI Join Date: Oct 2009 Posts: 80 Rep Power: 8 Hi doug can you send me your code?

 October 25, 2009, 03:25 #13 New Member   Join Date: Mar 2009 Posts: 7 Rep Power: 8 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: 7 I am interested. Is it for un structure grid? I Will try to write one.

 Tags boundary conditions, k-e, turbulence model

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36 Travis FLUENT 4 January 7, 2009 02:18 kharati Main CFD Forum 1 July 29, 2008 02:12 FredPacheo FLUENT 0 July 24, 2008 11:06

All times are GMT -4. The time now is 09:56.