
[Sponsors] 
January 26, 2011, 22:22 
pimpleDyMFoam stability problems

#1 
Senior Member
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 406
Rep Power: 15 
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.6ext, 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 ke  using a GGI between the nonrotating/rotating sections.  initial condition is an MRFSimpleFoam solution.  running adaptive timestepping 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 2D mesh Overall domain bounding box (3.175 3.175 0) (3.175 3.175 0.00127) Mesh (nonempty, nonwedge) directions (1 1 0) Mesh (nonempty) directions (1 1 0) Mesh (nonempty, nonwedge) dimensions 2 All edges aligned with or perpendicular to nonempty directions. Boundary openness (1.547876833e20 1.278680862e20 2.545349965e23) Threshold = 1e06 OK. Max cell openness = 1.030359725e15 OK. Max aspect ratio = 222.5044955 OK. Minumum face area = 2.787829408e10. Maximum face area = 0.06796272984. Face area magnitudes OK. Min volume = 3.540543349e13. Max volume = 8.63126669e05. Total volume = 0.05112695654. Cell volumes OK. Mesh nonorthogonality Max: 58.16212739 average: 8.861513366 Threshold = 70 Nonorthogonality check OK. Face pyramids OK. Max skewness = 1.91991223 OK. Mesh OK. 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 = 1e8, relTol = 0) everything else: PBiCG (tolerance = 1e12, 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.019988529e07, Final residual = 8.180846808e15, No Iterations 3 DILUPBiCG: Solving for Uy, Initial residual = 1.146793419e06, Final residual = 3.853642016e13, 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.338602432e11, global = 3.312247682e13, cumulative = 3.86768347e10 GAMG: Solving for p, Initial residual = 0.1122173235, Final residual = 0.001055831713, No Iterations 7 time step continuity errors : sum local = 3.632168291e11, global = 2.825140759e13, cumulative = 3.870508611e10 DILUPBiCG: Solving for epsilon, Initial residual = 5.015908983e08, Final residual = 3.506292395e13, No Iterations 2 DILUPBiCG: Solving for k, Initial residual = 1.795206307e07, Final residual = 1.091405791e13, No Iterations 3 DILUPBiCG: Solving for Ux, Initial residual = 9.373482943e08, Final residual = 1.385869522e13, No Iterations 3 DILUPBiCG: Solving for Uy, Initial residual = 5.330058884e07, Final residual = 6.507134342e14, 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.027277872e11, global = 1.753070309e15, cumulative = 3.87049108e10 GAMG: Solving for p, Initial residual = 0.003209023427, Final residual = 9.919113054e09, No Iterations 300 time step continuity errors : sum local = 3.890095611e16, global = 1.955326102e17, cumulative = 3.870490885e10 DILUPBiCG: Solving for epsilon, Initial residual = 1.597098642e08, Final residual = 2.082955361e13, No Iterations 2 DILUPBiCG: Solving for k, Initial residual = 5.561266417e08, Final residual = 1.619021747e14, No Iterations 3 ExecutionTime = 2036.12 s ClockTime = 2040 s Courant Number mean: 0.001213163235 max: 1.990485424 velocity magnitude: 21.07635725 deltaT = 1.553405431e06 GGI pair (outerGGI, innerGGI) : 0.001002355386 0.001002444145 Diff = 3.192354117e08 or 0.00318485256 % 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. 

January 28, 2011, 06:36 

#2 
Senior Member
Join Date: Apr 2010
Posts: 151
Rep Power: 9 
Hello,
I have a similar issue. I run turbDyMFoam in 1.5dev, 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? 

January 29, 2011, 13:27 

#3 
Senior Member
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 406
Rep Power: 15 
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(); } Futhermore, I've finally been able to switch over to secondorder accurate div and ddt schemes. X'ing my fingers ... 

January 29, 2011, 13:36 

#4 
Senior Member
Chris Sideroff
Join Date: Mar 2009
Location: Ottawa, ON, CAN
Posts: 406
Rep Power: 15 
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.


Thread Tools  
Display Modes  


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 
Semiimplicit RungeKutta  Bobby  Main CFD Forum  10  August 5, 2005 15:43 