[cnsidero]

I'm having some difficulties getting a pimpleDyMFoam calculation to run. It will run for a long simulation time (as in thousands of time steps) but the calculation always seems to suddenly blow up. However, it needs to run longer and with second order accurate schemes. A quick summary of what I'm trying to do (I can't show a picture of the actual geometry):

- running 1.6-ext, pimpleDyMFoam
- 2D incompressible flow (V_inf = 10 m/s, L_ref ~= 0.5 m, air)
- uniform, horizontal flow with a cylindrical rotating section (omega = 3000 RPM, R_ref ~= 0.053 m => V_tip ~= 2.65 m/s)
- turbulent flow, using standard k-e
- using a GGI between the non-rotating/rotating sections.
- initial condition is an MRFSimpleFoam solution.
- running adaptive time-stepping with CFL = 2.0

The mesh is a mix of quad and tri cells (quad in BL regions, tri's elsewhere). Here's the last bit of the checkMesh output:

Checking geometry...
    This is a 2-D mesh
    Overall domain bounding box (-3.175 -3.175 0) (3.175 3.175 0.00127)
    Mesh (non-empty, non-wedge) directions (1 1 0)
    Mesh (non-empty) directions (1 1 0)
    Mesh (non-empty, non-wedge) dimensions 2
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (1.547876833e-20 -1.278680862e-20 -2.545349965e-23) Threshold = 1e-06 OK.
    Max cell openness = 1.030359725e-15 OK.
    Max aspect ratio = 222.5044955 OK.
    Minumum face area = 2.787829408e-10. Maximum face area = 0.06796272984.  Face area magnitudes OK.
    Min volume = 3.540543349e-13. Max volume = 8.63126669e-05.  Total volume = 0.05112695654.  Cell volumes OK.
    Mesh non-orthogonality Max: 58.16212739 average: 8.861513366 Threshold = 70
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.91991223 OK.

Mesh OK.

so the mesh appears to be good. (the aspect ratio is due to the extrusion of the 2D mesh in the 3rd dim)

I've attached my controlDict, fvSchemes, fvSolution and dynamicMeshDict for completeness but here's a short summary for each:

fvSchemes:

ddt: Euler
grad: cellMDLimited Gauss linear 1.0
div: upwind

As you can see, pretty much first order for everything (I'm just trying to get it to run!)

fvSolution:

p: GAMG (tolerance = 1e-8, relTol = 0)
everything else: PBiCG (tolerance = 1e-12, relTol = 0)
nOuterCorrectors: 2
nCorrectors: 2
nNonOrthogonalCorrectors: 0

Here's the output from the second last time before it blows up:

Initializing the GGI interpolator between master/shadow patches: outerGGI/innerGGI
DILUPBiCG:  Solving for Ux, Initial residual = 2.019988529e-07, Final residual = 8.180846808e-15, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 1.146793419e-06, Final residual = 3.853642016e-13, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.1681149053, Final residual = 0.001629595763, No Iterations 6
time step continuity errors : sum local = 5.338602432e-11, global = 3.312247682e-13, cumulative = 3.86768347e-10
GAMG:  Solving for p, Initial residual = 0.1122173235, Final residual = 0.001055831713, No Iterations 7
time step continuity errors : sum local = 3.632168291e-11, global = 2.825140759e-13, cumulative = 3.870508611e-10
DILUPBiCG:  Solving for epsilon, Initial residual = 5.015908983e-08, Final residual = 3.506292395e-13, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 1.795206307e-07, Final residual = 1.091405791e-13, No Iterations 3
DILUPBiCG:  Solving for Ux, Initial residual = 9.373482943e-08, Final residual = 1.385869522e-13, No Iterations 3
DILUPBiCG:  Solving for Uy, Initial residual = 5.330058884e-07, Final residual = 6.507134342e-14, No Iterations 8
GAMG:  Solving for p, Initial residual = 0.05525130139, Final residual = 0.0004889471911, No Iterations 5
time step continuity errors : sum local = 2.027277872e-11, global = -1.753070309e-15, cumulative = 3.87049108e-10
GAMG:  Solving for p, Initial residual = 0.003209023427, Final residual = 9.919113054e-09, No Iterations 300
time step continuity errors : sum local = 3.890095611e-16, global = -1.955326102e-17, cumulative = 3.870490885e-10
DILUPBiCG:  Solving for epsilon, Initial residual = 1.597098642e-08, Final residual = 2.082955361e-13, No Iterations 2
DILUPBiCG:  Solving for k, Initial residual = 5.561266417e-08, Final residual = 1.619021747e-14, No Iterations 3
ExecutionTime = 2036.12 s  ClockTime = 2040 s

Courant Number mean: 0.001213163235 max: 1.990485424 velocity magnitude: 21.07635725
deltaT = 1.553405431e-06
GGI pair (outerGGI, innerGGI) : 0.001002355386 0.001002444145 Diff = -3.192354117e-08 or 0.00318485256 %

To me every looks pretty good except maybe the pressure. The last "Initial residual" only gets to 0.0032. The other thing I noticed is the number of iterations in the last pressure update goes up dramatically - throughout most of the simulation it only needs 30-40 but the second last time step needs 300, the last one goes to 1000. Another indicator pressure might be the guilty culprit.

I've checked the mechanics of everything and all seems to be working (e.g. the ggi mass flux error is small). I looked at several instances of the solution prior to it blowing up with Paraview: the rotating zone is moving correctly and there are no spikes in velocity, pressure or the turbulence quantities.

This all leads me to believe it's has to do with the pressure equation. I've tried different # of nCorrectors, outerCorrectorsn NonOrthogonalCorrectors and numerical schemes (I'm using the most stable ones already!) but it doesn't seem help. As it appears to be the pressure equation blowing up, what should I tweak to get better control/convergence of the pressure?

Sorry for the long post but I'm running out of things to try so I was hoping some fresh eyes and someone with more knowledgeable about the inner workings of pimple* type solvers might be able to send me down a better rabbit hole. Let me know if you need more information.

[flowris]

Hello,

I have a similar issue. I run turbDyMFoam in 1.5-dev, and I also use GAMG for my p calculation. This also blows up my number of iterations, and my speed increases very suddenly to unphysically high values.

I am running the same case with PBiCG for p as well as for the other variables, and until now it seems stable (but is much slower). Maybe there is a difficulty in combining GGI with grid coarsening?

[cnsidero]

Thanks for the tip flowris. I tried that but unfortunately that didn't make a difference for me.

I took a look at how the pimple*Foam works. It appears the outer corrector updates all the equations whereas the inner corrector updates only the pressure equation. Here's the code from pimpleFoam:

for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
        {
            if (nOuterCorr != 1)
            {
                p.storePrevIter();
            }

#           include "UEqn.H"

            // --- PISO loop
            for (int corr=0; corr<nCorr; corr++)
            {
#               include "pEqn.H"
            }

            turbulence->correct();
        }

Since it seemed the pressure equation was giving me difficulties, I just turned up the number of nCorrectors (inner corrections). It seems 8 nCorrectors is working. My relTol for p is only 2 orders of magnitude so the extra overhead of the additional corrector steps (I was only doing 2 before) isn't too bad. Actually it might even being doing fewer total iterations because the final pressure correction is not needing as many to converge (I have relTol for pFinal = 0).

Futhermore, I've finally been able to switch over to second-order accurate div and ddt schemes. X'ing my fingers ...

[cnsidero]

Another issue I've discovered with pimple*Foam is that you still can not use large CFL numbers. At least for my applications I've found anything over a CFL of 2 and I start having stability issues. I suppose it's better than a PISO algorithm but you still have keep the CFL within reason.