CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

pimpleDyMFoam stability problems

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   January 26, 2011, 22:22
Default pimpleDyMFoam stability problems
  #1
Senior Member
 
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 369
Rep Power: 12
cnsidero is on a distinguished road
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:

Code:
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:

Code:
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.
Attached Files
File Type: txt dynamicMeshDict.txt (1.2 KB, 58 views)
File Type: txt controlDict.txt (2.2 KB, 58 views)
File Type: txt fvSchemes.txt (1.5 KB, 89 views)
File Type: txt fvSolution.txt (1.9 KB, 103 views)
cnsidero is offline   Reply With Quote

Old   January 28, 2011, 06:36
Default
  #2
Senior Member
 
Join Date: Apr 2010
Posts: 151
Rep Power: 7
flowris is on a distinguished road
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?
flowris is offline   Reply With Quote

Old   January 29, 2011, 13:27
Default
  #3
Senior Member
 
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 369
Rep Power: 12
cnsidero is on a distinguished road
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:

Code:
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 is offline   Reply With Quote

Old   January 29, 2011, 13:36
Default
  #4
Senior Member
 
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 369
Rep Power: 12
cnsidero is on a distinguished road
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.
cnsidero is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
SLTS / CoEuler stability problems deepblue17 OpenFOAM Running, Solving & CFD 1 February 24, 2013 16:23
rhoPisoFoam stability problems at low Mach ivan_cozza OpenFOAM 1 October 1, 2010 11:03
Stability problems with kepsilon in external aero edreed OpenFOAM Running, Solving & CFD 21 July 16, 2008 15:00
Stability startup problems with rhoSimpleFoam olesen OpenFOAM Running, Solving & CFD 1 July 18, 2006 08:09
Semi-implicit Runge-Kutta Bobby Main CFD Forum 10 August 5, 2005 15:43


All times are GMT -4. The time now is 18:57.