[KYPCK444]

Hi all,

I am trying to simulate a water flow against elliptical body, adopting kEpsilon model and simpleFoam as solver. The body is inside a pipe, and the whole domain is "split" in a half (I am considering a half of the pipe+body). I am adopting "freestream approach" for che cylindrical surface of the pipe (but treating it as wall does not lead to different results).

My problem is with p residuals, which show very high oscillating behaviour (see attached image), while the other quantities show good convergence. Run stops either after few timesteps, but I think I cannot accept the results.

I think the mesh (3,4 billions cells) should be ok.

I searched the forum to find out what I could change but I cannot get to a solution.

These are the case files:

PRESSURE------------------------

internalField uniform 0;
boundaryField
{
inlet
{ type zeroGradient;}
outlet
{ type fixedValue;
value $internalField;}
periphery
{ type freestreamPressure;}
symmetryPlane
{ type symmetry; }
body
{ type zeroGradient; }
}


VELOCITY--------------------------------------

internalField uniform (0 0 1.5);
boundaryField
{
inlet
{ type fixedValue;
value uniform (0 0 1.5);}
outlet
{ type zeroGradient;}
periphery
{ type freestream;
freestreamValue $internalField;}
symmetryPlane
{ type symmetry;}
body
{type noSlip;}
}


K---------------------------------------------

internalField uniform 0.0030;

boundaryField
{ inlet
{ type fixedValue;
value $internalField; }
outlet
{ type zeroGradient;}
periphery
{ type freestream;
freestreamValue $internalField;}
symmetryPlane
{ type symmetry;}
body
{ type kqRWallFunction;
value uniform 0;}
}


EPSILON-------------------------------------

internalField uniform 0.0003;
boundaryField
{ inlet
{ type zeroGradient;}
outlet
{ type zeroGradient;}
periphery
{ type freestream;
freestreamValue $internalField;}
symmetryPlane
{ type symmetry;}
body
{ type epsilonWallFunction;
value $internalField;}
}


NUt---------------------------------------------

internalField uniform 0.0026;

boundaryField
{ inlet
{ type calculated;
value $internalField;}
outlet
{ type calculated;
value $internalField;}
periphery
{ type freestream;
freestreamValue $internalField;}
symmetryPlane
{ type symmetry;}
body
{ type nutkWallFunction;
value uniform 0; }
}


fvSchemes--------------------------------------------

ddtSchemes
{
default steadyState;
}

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

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

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}


fvSolution----------------------------------------------

solvers
{
p
{
solver GAMG;
maxIter 500;
smoother GaussSeidel;

tolerance 1e-10;
relTol 0.01;
}

pFinal
{
$p;
tolerance 1e-10;
relTol 0;
}

"(U|k|epsilon)"
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0.1;
}

"(U|k|epsilon)Final"
{
$U;
tolerance 1e-05;
relTol 0;
}
}

PIMPLE
{
nOuterCorrectors 4;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}

SIMPLE
{

pRefCell 0;
pRefValue 0;

nNonOrthogonalCorrectors 3;
nCorrector 1;
residualControl
{
p 1e-2;
U 1e-3;
"(k|epsilon)" 1e-3;
}
}

relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.7;
k 0.7;
"epsilon.*" 0.7;
}
}

cache
{
grad(U);
}

-----------------------------

Any hint would be really appreciated, thanks.

[piu58]

What is your time step. How large is the max Courant number? If you have problems firts try the get maxCo < 0.1.

[KYPCK444]

Thanks for your reply.

Time step is set to 1, max Courant number (calculated as dt*U/d, where minimum cell dimension is calculated as the square root of minimum cell area, as stated in checkMesh) should be 885,
but I was not taking into account this, because I knew Co is not involved in steady simulations (Smallest element, simpleFoam and CFL condition).

Have I misunderstood something?

Thanks again

[piu58]

Even in steady simulations the Co number says something about your simulation. The larges change in the values calculated spans about 900 cells! That is much to high.

[akidess]

Originally Posted by piu58

Even in steady simulations the Co number says something about your simulation. The larges change in the values calculated spans about 900 cells! That is much to high.

Courant number is 0, since dt is 0. Your reasoning is odd.

[KYPCK444]

Well, I obtain same risults with delta t = 0,001 (large Co) and = 1.

Anyway, I have got better residuals by changing fvSchemes as below:

fvSchemes ---------------------------------------------------------------

ddtSchemes
{
default steadyState;
}

gradSchemes
{
default Gauss linear;
}

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

laplacianSchemes
{
default Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}
--------------------------------------------------------------------------

Now, the very problem is that in spite of better residuals behaviour, my run does not converge, I see it from drag coefficient plot over time not converging to a stable value.

I already adopt initialization with potentialFoam (adding potentialFLow line fvSolution), and I'm monitoring yPlus (in my case yPlus is between 2,5 and 65; I know RAS yPlus should not be less than 30, I don't know if this could be the problem).
I have also tried to increase domain lenght past the submersed pipe, but things do not change.

The "strange" feature is that the simulation stops after few timesteps (ca 100), saying it has converged ..

Attached my case files (as you can see in bc, I am now treating the peripherical surface of the domain as a wall).

Any hint? Thanks