CFD Online Discussion Forums - SimpleFoam bounded k and epsilon always diverging

Hello!

For some time now I'm having trouble with my simulations. I was hoping someone would be able to advise me on something I might be doing wrong. Full problem statement at the end.

OF version:
My version of OpenFOAM is the OpenFOAM for Windows from CFD-support. The version is 3.0.x.

My simulation in short:
I'm looking into pressure drop in a certain pipe configuration (single pipe, splits to double pipe and then single pipe again) using simpleFoam. So basically the geometry consists of an inlet patch, outlet patch and the walls. Since I'm looking for a pressure drop from constant inlet conditions, I'm looking for the steady state result.

I've attached my fvSchemes, fvSolver files, below I stated some I.C. (0). I've only put k and omega here, but pretty much all turbulence models give divergent results.

epsilon:

Code:

dimensions      [0 2 -3 0 0 0 0];



internalField   uniform 0.322;



boundaryField

{

    inlet

    {

        type            fixedValue;

        value           uniform 0.322;

    }

    outlet

    {

        type            zeroGradient;

    }

    walls

    {

        type            epsilonWallFunction;

        value           uniform 0.322;

    }

}

Code:

dimensions      [0 2 -2 0 0 0 0];



internalField   uniform 0.0139;



boundaryField

{

    inlet

    {

        type            fixedValue;

        value           uniform 0.0139;

    }

    outlet

    {

        type            zeroGradient;

    }

    walls

    {

        type            kqRWallFunction;

        value           uniform 0.0139;

    }

}

Code:

dimensions      [0 2 -2 0 0 0 0];



internalField   uniform 0;



boundaryField

{

    inlet

    {

        type            zeroGradient;

    }



    outlet

    {

        type            fixedValue;

        value           uniform 0;

    }



    walls

    {

        type            zeroGradient;

    }



}

Code:

dimensions      [0 1 -1 0 0 0 0];



internalField   uniform (0 0 0);



boundaryField

{

    inlet

    {

        type            fixedValue;

        value           uniform (1 0 0);

    }



    outlet

    {

        type            zeroGradient;

    }



    walls

    {

        type            fixedValue;

        value           uniform (0 0 0);

    }



}

The problem
The problem is that both k and epsilon (or omega for example, when using kOmega turbulence scheme) diverge after less than 10 iterations. No matter the timestep. The strange thing is that some time ago I did more or less the same simulations and these got results (but were quite off from actual measurements). I didn't change anything in any source code, so this seems really strange to me. The error message, as expected with divergence, is:

Code:

Time = 0.0007



smoothSolver:  Solving for Ux, Initial residual = 0.0092005991, Final residual = 0.00036478284, No Iterations 1

smoothSolver:  Solving for Uy, Initial residual = 0.32864692, Final residual = 0.0107253, No Iterations 1

smoothSolver:  Solving for Uz, Initial residual = 0.013299208, Final residual = 0.00017554984, No Iterations 1

GAMG:  Solving for p, Initial residual = 1, Final residual = 0.083184461, No Iterations 3

time step continuity errors : sum local = 2.5183777e+063, global = -3.1488013e+062, cumulative = -3.1488013e+062

smoothSolver:  Solving for epsilon, Initial residual = 4.5122245e-007, Final residual = 4.5122245e-007, No Iterations 0

bounding epsilon, min: 2.2586618e-018 max: 4.3673913e+082 average: 3.8104317e+077

smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.098475985, No Iterations 11

bounding k, min: -1.8065296e+098 max: 1.5304793e+103 average: 3.7462759e+097

ExecutionTime = 18.737 s  ClockTime = 19 s



Time = 0.0008



smoothSolver:  Solving for Ux, Initial residual = 0.064686933, Final residual = 0.0037937872, No Iterations 1

smoothSolver:  Solving for Uy, Initial residual = 0.10249877, Final residual = 0.0076579062, No Iterations 1

smoothSolver:  Solving for Uz, Initial residual = 0.62092711, Final residual = 0.027387706, No Iterations 1

GAMG:  Solving for p, Initial residual = 4.3415891e-024, Final residual = 4.3415891e-024, No Iterations 0

time step continuity errors : sum local = 2.7715949e+065, global = 1.163728e+064, cumulative = 1.13224e+064

smoothSolver:  Solving for epsilon, Initial residual = nan, Final residual = nan, No Iterations 1000

bounding epsilon, min: 8.4530623e-020 max: 1.0590742e+157 average: nan

smoothSolver:  Solving for k, Initial residual = nan, Final residual = nan, No Iterations 1000

ExecutionTime = 66.668 s  ClockTime = 67 s







--> FOAM FATAL IO ERROR:

wrong token type - expected Scalar, found on line 0 the word 'nan'



file: C:/OpenFOAM/z003nhwx-3.0.x/run/DN25simple/system/data.solverPerformance.epsilon at line 0.



    From function operator>>(Istream&, Scalar&)

    in file lnInclude/Scalar.C at line 93.



FOAM exiting

Steps taken so far
Since I have this problem for some time now, I have tried several thing:

1. Changing initial and boundary conditions
I've changed more or less most initial conditions to previous values, small values, large values and even zero. Only the case for no flow (u inlet = 0) and no stated pressure, only pressureGradient = 0) gave a result, the expected result of nothing happening. Then changing the inlet velocity to 1 m/s gave again an exploding behavior for the turbulence parameters.

2. Laminar flow
To look if the problem truly resides with the turbulence model I performed the case above for laminar flow with icoFoam. This gave a result which was more or less in line of the results I got when all was working.

3. Changing fvSchemes
I've changed the divergent schemes, especially for epsilon and k, but all resulted in exploding behavior. This also means having bounded or not bounded schemes.

4. Changing solver and relaxation factors
Changing the solver (used both GAMG and PCG) gave the same result. Also changing the relaxation factor (0.5, 0.7, 0.9, 0.98, 1) didn't led to any converging behavior.

5. Adapting mesh
Regarding mesh: when running checkMesh gives an OK mesh. Despite this, I tried running the simulation with a mesh with and without BL refinement. No difference here.

6. Changing timestep
As mentioned, changing the timestep didn't do anything to the divergent behavior. Ranging from 1s to 0.0000001s gave massive values within 10 iterations.

7. Running both in parallel as series
To rule out an error in parallel running, I ran the simulations in parallel and without MPI. No change whatsoever.

Any help would be very very welcome, since I'm getting a little frustrated with this and my time on this project is limited. Thank you very much in advance!

Dennis