[andrewburns]

I have a fully converged solution (residuals between 1e-6 and 1e-8) using first order upwind discretisation. I want to switch to second order to get more accurate drag results however every time I try switching to a second order scheme my k and epsilon values go unbounded (negative) within a few iterations, stay that way and my forces go dodgy.

My discretization schemes are currently:

gradSchemes
{
default Gauss linear;
grad(p) Gauss midPoint;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss linear;
div(R) Gauss linear;
div(phi,nuTilda) Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian (nuEff,U) Gauss linear corrected;
laplacian ((1|A(U)),p) Gauss linear limited 1;
laplacian (DkEff,k) Gauss linear limited 1;
laplacian (DepsilonEff,epsilon) Gauss linear limited 1;
laplacian (DREff,R) Gauss linear corrected;
laplacian (DnuTildaEff,nuTilda) Gauss linear corrected;
}

Running a simple steady-state solution in simpleFoam. My current relaxation factors are

p 0.3
U 0.5
k 0.3
epsilon 0.2
R 0.7
nuTilda 0.7

I've also tried:

QUICK
linearUpwind (which threw me an error for some reason)
limitedVanLeer (strictly bound to be positive for k and epsilon, which it appeared to ignore)
SFCD

Of all of those the most stable one appeared to be SFCD but even so k and epsilon unbounded fairly rapidly.

Can someone help me out? I know other people must be able to perform successful second order solutions in openfoam but I just can't make it work.

[andrewburns]

Can anyone help me here? It seems that if I take any fully converged first order solution and then switch to second order my k and epsilon (and omega when using komegaSST) go unbounded and stay there.