[Arya_the_Rabbit]

Hi guys!
I have some problems with the stability and convergence of a case with 70 millions cells and complicated automotive geometry.

-Mesh quality (tetra with prism layer):
skew <3.5
nonOrtho <60

I'm using a k-omega SST turbulence model with an high reynolds number approach.

During the first 600 iterations I use the first order scheme below:
ddtSchemes
{
default steadyState;
}

gradSchemes
{
default cellLimited leastSquares 1;
grad(U) cellLimited leastSquares 1;
}

divSchemes
{
default none;
/*--------------Schemi numerici modello STD KOMEGASST-----------------------*/
div(phi,U) bounded Gauss upwind grad(U);
div(phi,k) bounded Gauss upwind;
div(phi,omega) bounded Gauss upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
div(div(phi,U)) Gauss linear;
}

laplacianSchemes
{
default Gauss linear limited 0.5;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default limited 0.5;
}

fluxRequired
{
default no;
p;
}

wallDist
{
method meshWave;

and coupled solver:
solvers
{
p
{
solver GAMG;
smoother GaussSeidel;
tolerance 1e-7;
relTol 0.01;
nPreSweep 0;
nPostSweep 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
}

Phi
{
$p;
}



U
{
type coupled;
solver PBiCICG;
preconditioner DILU;
tolerance (1e-8 1e-8 1e-8);
relTol (0 0 0);
}

"(k|omega)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-8;
relTol 0.1;
}
/*---------------Aggiunti per il modello termico-------------*/
"(h|e)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.01;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 0;
/*--------------Aggiunti per il modello termico---------------*/
rhoMin rhoMin [1 -3 0 0 0] 0.5;
rhoMax rhoMax [1 -3 0 0 0] 2.0;
}

potentialFlow
{
nNonOrthogonalCorrectors 10;
}

relaxationFactors
{
fields
{
p 0.2;
rho 0.3;
}
equations
{
U 0.7;
"(k|omega)" 0.5;
"(h|e)" 0.7;
}
}

cache
{
grad(U);
}

Then a second order scheme:
ddtSchemes
{
default steadyState;
}

gradSchemes
{
default cellMDLimited Gauss linear 1;
grad(U) cellMDLimited Gauss linear 1;
}

divSchemes
{
default none;
div(phi,U) bounded Gauss SuperBeeV grad(U);
div(phi,k) bounded Gauss upwind;
div(phi,omega) bounded Gauss upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
div(div(phi,U)) Gauss linear;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}

wallDist
{
method meshWave;
}

and solver:
solvers
{
p
{
solver GAMG;
smoother GaussSeidel;
tolerance 1e-7;
relTol 0.01;
nPreSweep 0;
nPostSweep 2;
cacheAgglomeration on;
agglomerator faceAreaPair;
nCellsInCoarsestLevel 10;
mergeLevels 1;
}

Phi
{
$p;
}



U
{
type coupled;
solver PBiCICG;
preconditioner DILU;
tolerance (1e-8 1e-8 1e-8);
relTol (0 0 0);
minIter 3;
}

"(k|omega)"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-8;
relTol 0.1;
}
/*---------------Aggiunti per il modello termico-------------*/
"(h|e)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-8;
relTol 0.01;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 1;
/*--------------Aggiunti per il modello termico---------------*/
rhoMin rhoMin [1 -3 0 0 0] 0.5;
rhoMax rhoMax [1 -3 0 0 0] 2.0;
}

potentialFlow
{
nNonOrthogonalCorrectors 10;
}

relaxationFactors
{
fields
{
p 0.3;
rho 0.3;
}
equations
{
U 0.7;
"(k|omega)" 0.7;
"(h|e)" 0.7;
}
}

cache
{
grad(U);
}




I tried to use a coupled solver for a steady simulation for the better performance which the coupled solver should have. I've already tried with a segregated approach but the results are the same.

My problem is that the residual of the pressure never go under 10e-2 and the aerodynamic coefficients are unstable without a good convergence.
When I switch to second order approach the residual are stable but always high, especially the pressure, whereas the coefficients always unstable.

Is there some problems with the numerical schemes?
How can I improve my solver?

Thank you all for your advices.
Regards