# pimpleFoam blows:calculating laminar flow field

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

 April 25, 2014, 16:09 pimpleFoam blows:calculating laminar flow field #1 Member   India Join Date: Oct 2012 Posts: 84 Rep Power: 13 Hi Foamers, I am trying to calculate transient flow field in a single bifurcating tube(using OF-2.1.1 and swak4Foam) with inlet boundary condition varying with space and time: I will start illustrating my case by giving 'checkMesh' report: Create time Create polyMesh for time = 0 Time = 0 Mesh stats points: 140311 faces: 1127971 internal faces: 1083988 cells: 521531 boundary patches: 4 point zones: 0 face zones: 1 cell zones: 1 Overall number of cells of each type: hexahedra: 0 prisms: 125835 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 395696 polyhedra: 0 Checking topology... Boundary definition OK. Cell to face addressing OK. Point usage OK. Upper triangular ordering OK. Face vertices OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ... Patch Faces Points Surface topology INLET 910 612 ok (non-closed singly connected) OUTLET1 562 390 ok (non-closed singly connected) OUTLET2 566 394 ok (non-closed singly connected) WALL 41945 21065 ok (non-closed singly connected) Checking geometry... Overall domain bounding box (-0.108847 -0.161672 -0.0190498) (0.108851 0.127 0.0190498) Mesh (non-empty, non-wedge) directions (1 1 1) Mesh (non-empty) directions (1 1 1) Boundary openness (2.61084e-17 -1.48172e-16 3.15641e-17) OK. Max cell openness = 2.91142e-16 OK. Max aspect ratio = 7.00156 OK. Minumum face area = 1.85341e-07. Maximum face area = 3.64501e-06. Face area magnitudes OK. Min volume = 5.00262e-11. Max volume = 2.11586e-09. Total volume = 0.000347759. Cell volumes OK. Mesh non-orthogonality Max: 63.6578 average: 14.781 Non-orthogonality check OK. Face pyramids OK. Max skewness = 2.94579 OK. Coupled point location match (average 0) OK. Mesh OK. End Boundary Conditions: 0/U Code: ```dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { INLET { type groovyBC; variables "vavg=(-0.000163-0.00336*cos(time()*0.3254)-0.3416*sin(time()*0.3254));r2=(pow(pos().x,2)+pow(pos().z,2));R2=sum(area())/pi;para=-((R2-r2)/R2)*normal();"; valueExpression "2*vavg*para"; value uniform (0 0 0); } OUTLET1 { type zeroGradient; } OUTLET2 { type zeroGradient; } WALL { type fixedValue; value uniform (0 0 0); } }``` 0/p Code: ```dimensions [0 2 -2 0 0 0 0]; internalField uniform 0; boundaryField { INLET { type zeroGradient; } OUTLET1 { type fixedValue; value uniform 0; } OUTLET2 { type fixedValue; value uniform 0; } WALL { type zeroGradient; } }``` controlDict Code: ```libs ( "libOpenFOAM.so" "libsimpleSwakFunctionObjects.so" "libswakFunctionObjects.so" "libgroovyBC.so" ) ; application pimpleFoam; startFrom startTime; startTime 0; stopAt endTime; endTime 20; deltaT 0.02; writeControl adjustableRunTime; writeInterval 0.4; purgeWrite 0; writeFormat ascii; writePrecision 6; writeCompression uncompressed; timeFormat general; timePrecision 6; runTimeModifiable true; adjustTimeStep no; maxCo 100; functions { initSwak { // needed to allow DataEntry to work type initSwakFunctionObject; region region0; } }; functions { volFlowInlet { type swakExpression; valueType patch; patchName INLET; verbose true; expression "phi"; accumulations ( sum ); } volFlowOutlet1 { type swakExpression; valueType patch; patchName OUTLET1; verbose true; expression "phi"; accumulations ( sum ); } volFlowOutlet2 { type swakExpression; valueType patch; patchName OUTLET2; verbose true; expression "phi"; accumulations ( sum ); } }``` fvSchemes Code: ```ddtSchemes { default CrankNicholson 1; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(p_rgh) Gauss linear; grad(U) Gauss linear; } divSchemes { default none; div(phi,U) Gauss linearUpwind grad(U); div(phi,k) Gauss limitedLinear 1; div(phi,epsilon) Gauss limitedLinear 1; div(phi,R) Gauss limitedLinear 1; div(R) Gauss linear; div(phi,nuTilda) Gauss limitedLinear 1; div((nuEff*dev(T(grad(U))))) Gauss linear; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear corrected; laplacian((1|A(U)),p) Gauss linear corrected; laplacian((1|A(U)),p_rgh) Gauss linear corrected; laplacian(rAU,p_rgh) Gauss linear corrected; laplacian(DkEff,k) Gauss linear corrected; laplacian(DepsilonEff,epsilon) Gauss linear corrected; laplacian(DREff,R) Gauss linear corrected; laplacian(DnuTildaEff,nuTilda) Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; p_rgh ; }``` fvSolution Code: ```solvers { p { solver GAMG; tolerance 1e-06; relTol 0.01; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } pFinal { solver GAMG; tolerance 1e-06; relTol 0.01; smoother GaussSeidel; cacheAgglomeration true; nCellsInCoarsestLevel 10; agglomerator faceAreaPair; mergeLevels 1; } "(U|kl|kt|omega)" { solver PBiCG; preconditioner DILU; tolerance 1e-06; relTol 0; } "(U|kl|kt|omega)Final" { \$U; tolerance 1e-06; relTol 0; } } PIMPLE { nOuterCorrectors 50; nCorrectors 1; nNonOrthogonalCorrectors 2; pRefCell 0; pRefValue 0; residualControl { "(p|U)" { tolerance 1e-4; relTol 0; } } } relaxationFactors { fields { } equations { "U.*" 0.5; "kl.*" 1; "kt.*" 1; "omega.*" 1; } }``` The simulation blows blows at time=0.9s or little higher. The pulse lasts 20s and the reynolds number peaks(Re=2200) at 5s and 15s. If I choose adjustable time step then the it blows by showing time step in range of 10^-13 or so. If I choose a fixed time step of 0.01 or 0.02 or 0.005 it blows showing high courant number or unusually high flow rates. I have employed meshes of various kinds.....even pure tetrahedral mesh doesn't seem to work. Here is the link to download the case directory along with the mesh: https://www.dropbox.com/sh/u7sqlfvqtc6uynq/1MMKNryk7Z Please suggest me how to make it work .

 May 1, 2014, 15:38 #2 Member   India Join Date: Oct 2012 Posts: 84 Rep Power: 13 Somebody please help. Output at the time of explosion looks like this: Expression volFlowInlet : sum=0.000155198 Expression volFlowOutlet1 : sum=-1.01936e+12 Expression volFlowOutlet2 : sum=1.01938e+12 Courant Number mean: 4.53531e+15 max: 9.33679e+17 Time = 1.23 Also the time step continuity errors goes exponentially high: PIMPLE: iteration 50 DILUPBiCG: Solving for Ux, Initial residual = 0.0124649, Final residual = 3.77817e-07, No Iterations 5 DILUPBiCG: Solving for Uy, Initial residual = 0.00973289, Final residual = 3.76167e-07, No Iterations 5 DILUPBiCG: Solving for Uz, Initial residual = 0.0198487, Final residual = 4.315e-08, No Iterations 6 GAMG: Solving for p, Initial residual = 0.0537055, Final residual = 0.00027507, No Iterations 3 GAMG: Solving for p, Initial residual = 0.00744695, Final residual = 6.84879e-05, No Iterations 3 GAMG: Solving for p, Initial residual = 0.00143221, Final residual = 8.22524e-06, No Iterations 3 time step continuity errors : sum local = 3.47333e+09, global = -7.64618e+08, cumulative = -2.18106e+10 PIMPLE: not converged within 50 iterations ExecutionTime = 55306.7 s ClockTime = 88055 s Velocity fields have exploded and there is a flow reversal too. I tried avoiding flow reversal using outletinlet boudary condition since my outlets should behave as inlets in the beginning: Code: ```type outletInlet; outletValue uniform 0; value uniform 0;``` but it didn't help either.