Since, you are using Gauss-Seidel point to point iteration, then it does matter which one of the parameters calculated first. You may need to reduce the time step or you may use the Gauss-Seidel line relaxation method with an under-relaxation factor less than one.

Farouk is right on track... also don't start with zeroes. While you're solving the governing equations with G-S, you're also solving nonlinear contributions (the turbulence model) with feedback. This often takes finesse (you may have to try several things, so don't become discouraged). What works for one problem may not work for another. Keep at it and don't give up!

Dear All,

Thank you for your answers. I encountered a paper titled "Implementing a k-epsilon Turbulence Model in the FEniCS Finite Element Programming Environment". Although it uses finite elements, I saw that they are calculating the coefficients and variables in an order (at least this was what I deduced). Therefore I implement following procedure and had a solution:

1. Treat source terms half implicit half explicit: This means generations are left explicit and dissipation moved into main coefficient (implicit or source term linearization as Patankar named).

2. First calculate coefficients of k, then calculate k, then calculate coefficients of epsilon, then calculate epsilon, then calculate turbulent viscosity and repeat.

3. When calculating turbulent viscosity, limit its maximum value to for instance 1000 and if there is a limited turbulent viscosity point, calculate epsilon at that point by using limited turbulent viscosity and calculated k.

This procedure ensures that epsilon always have a proportionality with k.

I should also mention that in a recent book (Moukalled et al. - The Finite Volume Method in Computational Fluid Dynamics), i saw that they start solving epsilon first.

I should also mention that, the dissipation of k is linearized by writing e/k* where k* is the old k value of previous iteration.

I hope this information will be useful to someone.

Best Regards,

Eyub