
[Sponsors] 
No vortex shedding for flow over cylinder at Re = 59 using pimpleFoam 

LinkBack  Thread Tools  Display Modes 
May 11, 2016, 13:32 
No vortex shedding for flow over cylinder at Re = 59 using pimpleFoam

#1 
New Member
Katherine
Join Date: Feb 2016
Posts: 14
Rep Power: 3 
Hello,
I am using the solver pimpleFoam with a structured mesh generated in Cubit to solve for laminar flow over a cylinder, but I am having some problems at a Reynolds number of 59. The diameter of the cylinder is 1 m, and the velocity is 10 m/s in the xdirection. The value of kinematic viscosity I am using is 0.17 m^2/s; therefore, the Reynolds number is 59. I am using pimpleFoam with a time step of 0.1 seconds and I set my maximum Courant number to be 5 for this case. The mesh that I am using is structured and symmetric with 406,400 elements. The width of the domain is 100D and the length is 100D. Boundary conditions on the top and bottom of the domain are slip conditions for velocity. Here is the problem that I am having; after running my simulation for 160 seconds of flow time, the wake is still steady; there is no development of a vortex street behind the cylinder. I have included some pictures of the solution I have after 160 seconds of flow time. The vorticity magnitude is shown. I had obtained good results for the Strouhal number from simulations using the same mesh and setup for slightly higher Reynolds numbers (Re = 80 and Re = 100). I compared my values of the Strouhal number at these Reynolds numbers to the results from C. H. K. Williamson (1996) and my percent error is between 0.9  1.3%. These were obtained after 20 seconds of flow time. Upon suggestion from my advisor, I used the solution from the final timestep obtained for a Re of 100 and changed the Reynolds number to 59 to see what happens to the wake. After letting the simulation run for a flow time of 22 seconds, I obtained a Strouhal number of 0.133; C.H.K. Williamson found a Strouhal number of 0.135 at Re = 60, so I think this is fairly reasonable. However, I am very confused as to what went wrong with my original simulation for a Reynolds number of Re = 59. I have included my controlDict, fvSolution, and fvSchemes case files here to provide more information on my setup: Code:
FoamFile { version 2.0; format ascii; class dictionary; location "system"; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application pimpleFoam; startFrom startTime; startTime 121; // 0  121 s already obtained, start from last timestep stopAt endTime; endTime 160; deltaT 0.1; writeControl adjustableRunTime; writeInterval 0.05; purgeWrite 0; writeFormat ascii; writePrecision 12; writeCompression off; timeFormat general; timePrecision 12; runTimeModifiable true; adjustTimeStep yes; maxCo 5; //maxDeltaT 0.1; functions { forces { type forces; functionObjectLibs ("libforces.so"); outputControl timeStep; outputInterval 1; patches (fixedWalls); pName p; UName U; rhoName rhoInf; rhoInf 1.0; log true; CofR (0 0 0); writeFields yes; binData { nBin 90; //output data into 20 bins direction (1 0 0); // bin direction format gnuplot; cumulative yes; } } forceCoefficients { type forceCoeffs; functionObjectLibs ("libforces.so"); outputControl timeStep; outputInterval 1; patches (fixedWalls); writeFields yes; log true; patches (fixedWalls); // on cyl only pName p; UName U; dragDir (1 0 0); liftDir (0 1 0); pitchAxis (0 0 1); CofR (0 0 0); magUInf 10; lRef 1; Aref 0.79; // cross sectional area rhoName rhoInf; rhoInf 1; origin (0 0 0); coordinateRotation { type EulerRotation degrees true; rotation (0 0 0); } binData { nBin 90; //output data into 20 bins direction (1 0 0); // bin direction format gnuplot; cumulative yes; } } } Code:
FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p // linear equation system solver for p { solver GAMG; // very efficient multigrid solver tolerance 1e04; // solver finishes if either absolute relTol 0.0; // tolerance is reached or the relative // tolerance here minIter 3; // a minimum number of iterations maxIter 100; // limitation of iterions number smoother DIC; // setting for GAMG nPreSweeps 1; // 1 for p, set to 0 for all other! nPostSweeps 2; // 2 is fine nFinestSweeps 2; // 2 is fine scaleCorrection true; // true is fine directSolveCoarsestLevel false; // false is fine cacheAgglomeration true; // on is fine; set to off, if dynamic // mesh refinement is used! nCellsInCoarsestLevel 128; // 500 is fine, // otherwise sqrt(number of cells) agglomerator faceAreaPair; // faceAreaPair is fine mergeLevels 1; // 1 is fine } pFinal { solver GAMG; tolerance 1e04; relTol 0.0; smoother DICGaussSeidel; nPreSweeps 1; nPostSweeps 0; nFinestSweeps 0; cacheAgglomeration true; nCellsInCoarsestLevel 128; agglomerator faceAreaPair; mergeLevels 1; } U { solver smoothSolver; smoother symGaussSeidel; tolerance 1e4; relTol 0.0; } UFinal { solver smoothSolver; smoother symGaussSeidel; tolerance 1e04; relTol 0.0; } } PIMPLE // residual control is set to false by default { nNonOrthogonalCorrectors 0; // corrPISO nCorrectors 2; // nCorrPISO; solving pressure field twice nOuterCorrectors 50; // nCorrPIMPLE; default is 1, and if you use default, this is really operating in PISO mode..not PIMPLE. 10 is good after a while; start with 50 pRefCell 0; pRefValue 0; residualControl // add residual control after you reached a good solution during one time step; not before { U { tolerance 1e4; relTol 0; } p { tolerance 5e3; relTol 0; } } } relaxationFactors // 0 < relaxationFactors <= 1 { fields { p 0.7; pFinal 0.7; } equations { U 0.3; UFinal 0.3; } } Code:
FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default CrankNicolson 1; // psi=1 pure C.N, psi=0 pure Euler } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; } divSchemes // s/t more to add here? { default none; div(phi,U) Gauss limitedLinear 1; // changed from linear div(phi,k) Gauss limitedLinear 1; div(phi,B) Gauss limitedLinear 1; div(B) Gauss linear; div(phi,nuTilda) Gauss limitedLinear 1; div((nuEff*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default corrected; } fluxRequired { default no; p ; } wallDist { method meshWave; } Code:
FoamFile { version 2.0; format ascii; class volScalarField; location "0"; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 2 0 0 0 0]; internalField uniform 0; boundaryField { symmetry_4 { type zeroGradient; } velocityinlet_5 { type zeroGradient; } outflow_6 { type fixedValue; value uniform 0; } fixedWalls { type zeroGradient; } frontAndBackPlanes { type empty; } } Code:
FoamFile { version 2.0; format ascii; class volVectorField; location "0"; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { symmetry_4 { type slip; } velocityinlet_5 { type uniformFixedValue; uniformValue table 4 ( (0 (10 5 0)) (1 (10 0 0)) (15 (10 1 0)) (16 (10 0 0)) ) ; value uniform (10 5 0); } outflow_6 { type zeroGradient; } fixedWalls { type fixedValue; value uniform (0 0 0); } frontAndBackPlanes { type empty; } } I hope I have explained my problem; please let me know if there is anything I have been unclear on. I would really appreciate any and all advice as to what could be my error in my setup. Thank you very much in advance for your time! 

May 13, 2016, 00:25 

#2 
Senior Member
Mahdi Hosseinali
Join Date: Apr 2009
Location: NB, Canada
Posts: 271
Rep Power: 11 
I think what is happening is that the second order discretization error is diffusive for that Reynolds number and damps the shedding at the very starts, so it keeps steady.
Similar case happens if you want to solve a channel flow with unperturbed (or randomly perturbed) initial conditions. In this situations a workaround usually is to use a good initial condition (also the case for channel flow) or a forcing function. Using higher reynolds number is an initial condition that lets the shedding happen, and when it's strong enough not to be damped by the numerical diffusion then it takes over. However this applies to every transient case and you probably know it, but as a reminder, you should let the flow simulation runs for a few eddy turn over before you start averaging or sampling. 

May 17, 2016, 14:31 

#3 
New Member
Katherine
Join Date: Feb 2016
Posts: 14
Rep Power: 3 
Thank you anishtain4 for your reply! Your suggestion of starting with a good initial condition was very helpful; I was able to obtain a Strouhal number of 0.133 after letting the simulation run for 22 seconds of flow time and starting with an initial condition obtained from the last time step for a case where the Reynolds number was 100. Williamson provides a St # of about 0.135 for a Reynolds number of 60, so I am pretty satisfied with this result. However, I continued to let the simulation run after this result, and after another 20 seconds, the wake becomes steady and the vortex shedding has damped out. I posted a video of the results on dropbox if anyone is curious: https://www.dropbox.com/s/fi6dwhw5lh...re_60.avi?dl=0. The first snapshot of the video is the initial condition obtained from the Re=100 case, and the last snapshot is the steady wake after 42 seconds of flow time at Re=60.
At this point I have a acceptable Strouhal number, but I am just curious about why the second order discretization error is diffusive for this Reynolds number. Is this a problem inherent to using the pimpleFoam solver? I would be curious to know which solver you used for a channel flow simulation as you mentioned in your post. Thank you again for your help! 

May 18, 2016, 01:05 

#4 
New Member
Join Date: Jun 2014
Posts: 9
Rep Power: 5 
Hi kaszt
Try making these changes: 1) Code:
maxCo 2; // use 1 if 2 doesnt work 2) Code:
solvers { p { solver GAMG; tolerance 1e08; relTol 0.1; smoother GaussSeidel; nPreSweeps 1; nPostSweeps 2; nFinestSweeps 0; scaleCorrection true; cacheAgglomeration true; nCellsInCoarsestLevel 128; agglomerator faceAreaPair; mergeLevels 1; } pFinal { solver GAMG; tolerance 1e08; relTol 0.0; smoother DICGaussSeidel; nPreSweeps 1; nPostSweeps 2; nFinestSweeps 0; cacheAgglomeration true; nCellsInCoarsestLevel 128; agglomerator faceAreaPair; mergeLevels 1; } U { solver smoothSolver; smoother symGaussSeidel; tolerance 1e8; relTol 0.1; } UFinal { solver smoothSolver; smoother symGaussSeidel; tolerance 1e08; relTol 0.0; } } PIMPLE { nNonOrthogonalCorrectors 0; nCorrectors 2; nOuterCorrectors 50; pRefCell 0; pRefValue 0; } residualControl { U { tolerance 1e8; relTol 0; } p { tolerance 1e8; relTol 0; } } relaxationFactors { fields { p 0.7; pFinal 1; } equations { U 0.3; UFinal 1; } } 3) Code:
div(phi,U) Gauss limitedLinearV 1; 

May 20, 2016, 03:17 

#5 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 20 
Katherine, I think what Mahdi ment is "second order discretization error is too diffusive". So every discretization scheme (except linear on an orthogonal mesh) introduces an error that behaves like additional diffusion (=higher viscosity). This error can damp the oscillations in a way that they don't occour. If this is actually the problem, you will have the same problem with any other solver.
In my experience, sometimes the vortex street takes a very long time to develop, but best way is to force the start as Mahdi wrote. Also, in my simulations it looked like the oscillations will start to occur at the outlet (on the right hand side in your pictures) and propagate upstream until they will occur at the bluff body. Have a look at this video here: https://www.youtube.com/watch?v=pqBATR2_kfg I can see your mesh is very fine at the cylinder but not that fine at the outlet, so maybe it is too coarse for this (at the outlet). Philipp.
__________________
The skeleton ran out of shampoo in the shower. 

May 31, 2016, 12:57 

#6 
New Member
Katherine
Join Date: Feb 2016
Posts: 14
Rep Power: 3 
Hi Philipp,
Thank you for sharing a video of your simulation; I can see that the formation of the vortex street starts in the wake about halfway near the outflow boundary, and propagates upstream. Since I am in the moderate Reynolds number range, there is an upstream influence that affects what happens at the cylinder. I did not make my mesh very fine at the outlet in order to conserve computational time, but I think that remeshing would be a good test to see what the difference is. Thank you for your help! Also thank you, Kire, for your suggestions; I tried using them exactly as you said, but I did not make much progress because advancing in time by 0.05 seconds took my computer about 23 hours...could I ask if you had happened to try this, and how long it took you to advance by one time step? (Whatever the size of your time step may be; I was using convective Courant number control to calculate my time step, and writing every 0.05 seconds of flow time). 

June 1, 2016, 00:40 

#7 
New Member
Join Date: Jun 2014
Posts: 9
Rep Power: 5 
Hi kaszt
I would use PISO and limit the time step so that Courant number is less than 1. You can also consider initialising the case with a steadystate result by running simpleFoam with/without turbulence off. I can't remember how long it took me to get the results. 

June 1, 2016, 02:19 

#8 
Senior Member
Philipp
Join Date: Jun 2011
Location: Germany
Posts: 1,297
Rep Power: 20 
Katherine, the video isn't mine, I just took some random hit from the search
Before trying the PISO stuff, please post some log output of your simulation. I don't think that you actually need to resolve Co<1 but you must resolve the shedding frequency properly, such as 20 time steps per cycle. So normally, this works with PIMPLE and pretty large time steps. Also, when I see your grid, this looks insanely fine. You will get a vortex street on a much coarser grid. Sometimes these instabilities grow much faster on a coarser grid. So if I were you, I would coarsen the grid (this is 2d, right) to something like less than 100000 cells. How many do you have right now?
__________________
The skeleton ran out of shampoo in the shower. 

Tags 
cylinder flow, pimplefoam, structured mesh 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Vortex Shedding in a 2D cylinder Using buoyantBoussinesqPimpleFoam  veerag96  OpenFOAM  0  April 16, 2016 22:22 
Cylinder vortex shedding (3D)  Apocolapse  STARCCM+  2  April 3, 2014 20:33 
Need advices for simulating vortex shedding of moving cylinder  quyetthang  CFX  0  October 1, 2011 05:39 
3D cylinder and Vortex Street shedding issues  josik_1982  FLUENT  2  July 17, 2010 10:19 
Motion of cylinder due to vortex shedding  ACHILLEAS TSOMPANOS  FLUENT  3  November 24, 2006 04:06 