[nikosb]

I am trying to run DNS of a sphere at Re=3700. I will be comparing the results of OpenFOAM against other codes and results in the literature.

I created an O grid for the mesh around a sphere in a cylindrical domain. The inlet is 5D from the sphere, the outlet is 20D and the radius of the cylindrical domain is 10D. Here are some photos of the mesh near the sphere. Around the sphere the grid is slightly stretched in the radial direction to cluster points near the surface of the sphere. The resolution is kept fine in a region behind the sphere up to 5D to resolve the near wake.

This photo showing a cross section cut showing the grid around the sphere and the wake



Another photo showing a cross section cut in the wake with crinkle clip




The simulation blows up after 20 iterations. I ran the simulation with both pisoFOAM and simpleFOAM without any success. Simulationtype is laminar. Here is the output of the pisoFOAM run.

Time = 0.036

Courant Number mean: 0.0422403 max: 1.36537
smoothSolver:  Solving for Ux, Initial residual = 0.0587041, Final residual = 8.42802e-06, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 0.0587602, Final residual = 7.79945e-06, No Iterations 5
smoothSolver:  Solving for Uz, Initial residual = 0.00212533, Final residual = 2.55993e-06, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.99465, Final residual = 0.0430159, No Iterations 2
time step continuity errors : sum local = 1.02157e-06, global = -6.00863e-12, cumulative = -3.02265e-08
GAMG:  Solving for p, Initial residual = 0.979142, Final residual = 8.77358e-07, No Iterations 21
time step continuity errors : sum local = 2.42567e-11, global = -3.09213e-12, cumulative = -3.02296e-08
ExecutionTime = 401.43 s  ClockTime = 403 s

forceCoeffs forces write:
    Cm    = -0.0331513
    Cd    = 0.697653
    Cl    = 0.0638549
    Cl(f) = -0.00122385
    Cl(r) = 0.0650787

Time = 0.038

Courant Number mean: 0.0422442 max: 2.85164
smoothSolver:  Solving for Ux, Initial residual = 0.0998439, Final residual = 8.97566e+216, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.100005, Final residual = 2.91442e+216, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.00432983, Final residual = 2.49464e+215, No Iterations 1000

The simulation blows up suddenly from one time step to the next. The time step is 2e-3. I ran with smaller time step of 2e-4 and I still get the same problem.

Here is the fvScheme file

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         Gauss linear;
}

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

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

and the fvSolutions

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

    pFinal
    {
        $p;
        smoother        DICGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|k|B|nuTilda|s)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
}

and the controDict

application     pisoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         300;

deltaT          2e-03;

adjustTimeStep  yes;

maxCo           1.0;

writeControl    adjustableRunTime;

writeInterval   10;

purgeWrite      0;

writeFormat     binary;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

Any help with why it is blowing up and what to try next?

[joshwilliams]

blows up after 20 iterations

This would suggest to me that the flow begins to develop (like flow separation behind the sphere). And then it is not being dissipated. I guess maybe you need a longer domain or finer mesh to promote dissipation of the turbulence by molecular viscosity. To see if the lack of viscous dissipation is an issue, you may try use an LES model (just as a test to check this is the issue).


Alternatively, it may be a numerical issue in your solution or boundary conditions. It seems your tolerances are quite low. We typically go for like 1e-10 for pressure and 1e-8 for velocity. Also why using two different solvers for pressure (GAMG and DICGaussSiedel)? I am unsure what pressure condition you use at the outlet but if there is still recirculation at the outlet you may try a stabilising BC (totalPressure for p, inletOutlet for U).