[jrr75166]

Hi OpenFOAM users!

I'm very very very new to CFD, and have been trying my hardest to learn OpenFOAM. Thus far, the '3 weeks' series of tutorials on the wiki has been a godsend, but I've encountered a problem with one of them.

Here's the pdf with all the information about the case setup: https://www.fluiddynamics.at/downloa...xample_Ten.pdf

I've followed all the steps up to running simpleFoam, but when I do it I get a floating point exception:

[stack trace]
=============
#1  Foam::sigFpe::sigHandler(int) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#2  ? in /lib/x86_64-linux-gnu/libc.so.6
#3  Foam::multiply(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so
#4  ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libincompressibleTurbulenceModels.so
#5  ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/lib/libincompressibleTurbulenceModels.so
#6  ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/bin/simpleFoam
#7  ? in /lib/x86_64-linux-gnu/libc.so.6
#8  __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#9  ? in /usr/lib/openfoam/openfoam2312/platforms/linux64GccDPInt32Opt/bin/simpleFoam
=============
Floating point exception

I know the error means I'm likely dividing by 0, and if you look at the last 2 iterations before it crashes:

Time = 16

smoothSolver:  Solving for Ux, Initial residual = 0.00672444, Final residual = 0.000187917, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.00620943, Final residual = 9.26625e-05, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 0.0153594, Final residual = 0.000618596, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.0183642, Final residual = 0.00045666, No Iterations 2
time step continuity errors : sum local = 6.28939e+78, global = -1.49885e+65, cumulative = -1.49885e+65
smoothSolver:  Solving for epsilon, Initial residual = 0.0864303, Final residual = 3.64626e-20, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 0.920353, Final residual = 0.0777223, No Iterations 1
bounding k, min: -2.2312e+102 max: 3.2152e+132 average: 1.51071e+129
ExecutionTime = 0.24 s  ClockTime = 1 s

Time = 17

smoothSolver:  Solving for Ux, Initial residual = 2.68236e-13, Final residual = 4.81414e-14, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 5.77589e-13, Final residual = 8.05644e-14, No Iterations 1
smoothSolver:  Solving for Uz, Initial residual = 2.57082e-14, Final residual = 2.85844e-15, No Iterations 1
GAMG:  Solving for p, Initial residual = 0.0256865, Final residual = 5.91851e-05, No Iterations 1
time step continuity errors : sum local = 3.03522e+100, global = -6.86672e+84, cumulative = -6.86672e+84

k and epsilon are clearly blowing up and the time step continuity errors are huge.

I've done some reaserch and tried as many fixes that have been suggested to others as I could, but nothing much has changed. Changing my fvSchemes managed to get it from 9 to 17 iterations, but that's still very little. At this point my best guess is that, since the pdf doesn't really specify the boundary condition values of k and epsilon, I must have chosen incorrectly, but I'm really not sure how to even begin working out what they should be. I've tried a few different ones, but none have led to anything.

Here are what I think are the relevant files, let me know if any others would be helpful!:

0/U:

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

internalField   uniform (0 0 0);

boundaryField
{
    inlet_one
    {
        type            fixedValue;
        value           uniform (0.1 0 0);
    }

    inlet_two
    {
        type            fixedValue;
        value           uniform (-0.2 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    walls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

}

0/p:

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

internalField   uniform 1e5;

boundaryField
{
    inlet_one
    {
        type            zeroGradient;
    }

    inlet_two
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    walls
    {
        type            zeroGradient;
    }
}

0/epsilon:

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

internalField   uniform 14.855;

boundaryField
{
    inlet_one
    {
        type            fixedValue;
        value           uniform 14.855;
    }

    inlet_two
    {
        type            fixedValue;
        value           uniform 14.855;
    }

    outlet
    {
        type            zeroGradient;
    }

    walls
    {
        type            epsilonWallFunction;
        value           uniform 14.855;
    }

}

0/k:

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

internalField   uniform 0.375;

boundaryField
{
    inlet_one
    {
        type            fixedValue;
        value           uniform 0.375;
    }

    inlet_two
    {
        type            fixedValue;
        value           uniform 0.375;
    }

    outlet
    {
        type            zeroGradient;
    }

    walls
    {
        type            kqRWallFunction;
        value           uniform 0.375;
    }
}

0/nut:

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

internalField   uniform 0;

boundaryField
{
    inlet_one
    {
        type            calculated;
        value           uniform 0;
    }

    inlet_two
    {
        type            calculated;
        value           uniform 0;
    }

    outlet
    {
        type            calculated;
        value           uniform 0;
    }

    walls
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
}

fvSchemes:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwind grad(U);

    turbulence      bounded Gauss limitedLinear 1;
    div(phi,k)      bounded Gauss linearUpwind grad(k);
    div(phi,epsilon) bounded Gauss linearUpwind grad(epsilon);
    div(phi,omega)  $turbulence;

    div(nonlinearStress) Gauss linear;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method          meshWave;
}

fvSolution:

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

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

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      yes;

    residualControl
    {
        p               1e-2;
        U               1e-3;
        "(k|epsilon|omega|f|v2)" 1e-3;
    }
}

relaxationFactors
{
    equations
    {
        U               0.95; // 0.9 is more stable but 0.95 more convergent
        ".*"            0.95; // 0.9 is more stable but 0.95 more convergent
    }
}

I can't post my blockMeshDict, but I've run checkMesh whenever I've changed anything and it always comes out saying there are no issues.

Thank you in advance for your help!