Mesh Convergence for Turbulence Model
 

Greetings, everyone!

I'm testing the implementation of the k-ε-a-b turbulence model [For example decribed here]. For this purpose i've written a 1D CFD solver, which uses finite volume method with Godunov-type scheme for subset of hyperbolic equations and implicit linear scheme for subset of parabolic equations.

As a test case i had tried to run Poggi experiment on Richtmyer–Meshkov instability [Some examples of such calculus], but faced some troubles. When i'm trying to run such task with only artificial (but rather big, about 1000 times as much as normal for gases) viscosity, diffusion coefficient and heat conductivity, i.e. just a Navier-Stokes system of equations without tubulence model, solution coverges rather well with mesh refinement. Both in time dependence of mixing zone width and, for example, temperature at the end of the calculation (2.5 ms).

But when i'm trying to include turbulence, results begin to show great dependence on the mesh. Turbulent viscosity varies from case to case, and decrease significantly with mesh refinement. As a result mixing zone width after reshock also strongly depends on mesh and decrease when a try more detailed meshes.

How can i achieve convergence of the solution with turbulence model?

Thanks for any help!

Rough solution
 

The problem seems to be connected with deviatoric part of Reynolds stress tensor in source terms of k (and eps) equation. The only way to deal with it i've found so far - is simply to bound this source term (highlighted by red frame) by the maximum value of 2E8 J/m^3/s in k equation (and proportionally in eps equation). By this convergence of mixing zone width seems to be significantly improved.