We model a 3D turbulent fluid, evolving toward a statistical equilibrium, by adding to the equations for the mean field (v, p) a term like −α∇ • (ℓ(x)Dv t ). This is of the Kelvin-Voigt form, where the Prandtl mixing length ℓ is not constant and vanishes at the solid walls. We get estimates for velocity v inx , that allow us to prove the existence and uniqueness of a regular-weak solutions (v, p) to the resulting system, for a given fixed eddy viscosity. We then prove a structural compactness result that highlights the robustness of the model. This allows us to pass to the limit in the quadratic source term in the equation for the turbulent kinetic energy k, which yields the existence of a weak solution to the corresponding Reynolds Averaged Navier-Stokes system satisfied by (v, p, k).