Coding the Launder-Sharma Model
 

We've coded the Launder-Sharma Lo-Reynolds number k-e model. However, the solution converges to negative values for k and epsilon. If we limit them so that k and epsilon are simply set to a very small positive number whenever they come out to be negative, the solution eventually converges to k=e=0 which is the laminar solution. However, if we don't limit k and epsilon, they eventually go negative and the solution diverges.

Has anyone experienced something like this? How did you get the model to converge to a viable turbulent solution?

Thanks.

Doug

There are a number of possibilties.
1. What convection scheme are you using? Unbounded schemes cause the result to diverge
2. How are you implementing the BC for the different boundaries? Can also cause divergence
3. Is it compressible or incompressible?
4. Staggered or collocated grid?

The case is for fully-developed channel flow for a channel with width, L.

1. We are coding it for 1-D fully-developed flow, so there are no convection terms. Only source terms and diffusion terms.
2. For boundary conditions we are using: U(0) = 0, k(0) = 0, e_tilda(0) = 0, U'(L/2) = 0, k'(L/2) = 0, e'(L/2) = 0. We are using second-order methods for estimating the second derivative of k on the wall.
3. The case is incompressible.
4. The grid is a collocated grid. Again, it is only 1-D, and we are storing the values at the cell faces. So, this can be seen as a finite-difference method rather than a finite volume method.

Any ideas what may be causing it to diverge?

From your BC it seems that you are using a low Re modification of the launder sharma model. Since it is incompressible, you may want to check for mass conservation. Another issue might be with the time marching if you are using explicit schemes/ diagonal dominance of implicit schmes.

1- You are not the first to report negative values for the dissipation, just google for (negative dissipation fluid mechanics), there are thousands of entries.
2- It depends on how and what equation you are using to calculate the dissipation
if you are using the basic definition
dissipation = nu * (rate of deformation tensor)^2
it is clear you can not get negative values even if the rate of deformation tensor is negative, but it is squared (s ij)^2
so I gather you are not using the basic definition.
What equation are you using, write it down and see how negative values could be produced
My guess is that there is coding problem. You could down load ISSAC (was written principally to test all turbulence models in Fortran) and check your code
Good luck

Thanks for your input. We've finally arrived at a converged solution. We found that the solution was very sensitive to the initial guess. If the initial guess doesn't have some of the basic properties of the solution, then it won't find the solution.

The results should not be so sensitive to the initial conditions. Try to look at some new implementations of the k epsilon model.