What about solving with a different variant of the standard kEpsilon? eg. realizable kEpsilon, takes additonal information from local flow - in case that stabilizes the solution process. I can only speak from my own experience here: when I monitored k epsilon with 1st and 2nd order accuracy, the Upwind showed a stable, mean flow solution, while linearUpwind solution showed a wavering field around this mean. This was obviously more realistic, but an important takeaway was that 2nd order introduced oscillations. It could be that the standard kEps model is missing some key aspect needed to prevent blow up of these.