[andreachr]

Hello!

I'm simulating a geometry that resembles a nozzle, using the k-epsilon model for turbulence. I have experimental data with which I can validate the results of the simulations. I define a mass flow rate at inlet, static pressure at outlet, fixed values of k and epsilon at inlet, temperature at inlet.

With rhoSimpleFoam, all my simulations converge (residuals less than 10^-5) and the values are very similar to the experimental values. But with sonicFoam, even though the residuals reach values lower than 10^-4, the results obtained are very different from the experimental values.

As an example, the calculated pressure at inlet varies around 25% between rhoSimpleFoam and sonicFoam, even though ALL boundary conditions are exactly the same. Does anyone know how the numerical schemes or solvers used might affect so much the results?

Here are my selected settings:

• For rhoSimpleFoam

fvSchemes

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;

div(phi,U) Gauss upwind;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,e) Gauss limitedLinear 1;
div(phi,epsilon) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,K) Gauss upwind;
div(phid,p) Gauss Linear;
div(phi,Ekp) Gauss upwind;
div((phi|interpolate(rho)),p) Gauss upwind;


}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}

fvSolution

solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.05;
smoother GaussSeidel;
cacheAgglomeration off;
nCellsInCoarsestLevel 20;
agglomerator faceAreaPair;
mergeLevels 1;
}

U
{
solver smoothSolver;
smoother GaussSeidel;
nSweeps 2;
tolerance 1e-06;
relTol 0.1;
}

e
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-06;
relTol 0.1;
}

"(k|epsilon)"
{
$U;
tolerance 1e-07;
relTol 0.1;
}
}

SIMPLE
{
nNonOrthogonalCorrectors 8;
nCorrectors 1;
pMin pMin [1 -1 -2 0 0 ] 15000;
rhoMin rhoMin [ 1 -3 0 0 0 ] 0.5;
rhoMax rhoMax [ 1 -3 0 0 0 ] 2;

residualControl
{
p 1e-5;
U 1e-5;
T 1e-5;
e 1e-4;

// possibly check turbulence fields
"(k|epsilon|omega)" 1e-5;
}
}

relaxationFactors
{
fields
{
p 0.1;
rho 0.001;
}
equations
{
U 0.3;
k 0.3;
epsilon 0.3;
e 0.3;
}
}

• For sonicFoam

fvSchemes


ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phid,p) Gauss limitedLinear 1;
div(phi,K) Gauss limitedLinear 1;
div(phi,e) Gauss limitedLinear 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss limitedLinear 1;
div(phiv,p) Gauss limitedLinear 1;
div(phi,omega) Gauss limitedLinear 1;
div((rho*R)) Gauss linear;
div(R) Gauss linear;
div(U) Gauss linear;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}

fvSolution

solvers
{
"rho.*"
{
solver diagonal;
}

"p.*"
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-12;
relTol 0;
}

"(U|e).*"
{
$p;
tolerance 1e-9;
}

"(k|epsilon).*"
{
$p;
tolerance 1e-10;
}
}

PIMPLE
{
nOuterCorrectors 200;
nCorrectors 2;
nNonOrthogonalCorrectors 4;
rhoMin rhoMin [ 1 -3 0 0 0 ] 0.2;
rhoMax rhoMax [ 1 -3 0 0 0 ] 20.0;

residualControl
{
"(U)"
{
relTol 0;
tolerance 1e-5;
}
"(k|epsilon|p)"
{
relTol 0;
tolerance 1e-5;
}

}
}

relaxationFactors
{
fields
{
"p.*" 0.3;
"rho.*" 0.01;
}
equations
{
"(U|e|k|epsilon|omega).*" 0.5;
}
}

The boundary conditions are exactly the same in both cases. The mesh is also ok (Mesh non-orthogonality Max: 60.7909 average: 7.41921).

I would really appreciate if someone can help with this problem!! Thanks!!!