CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   No vortex shedding for flow over cylinder at Re = 59 using pimpleFoam (https://www.cfd-online.com/Forums/openfoam-solving/171586-no-vortex-shedding-flow-over-cylinder-re-59-using-pimplefoam.html)

 kaszt May 11, 2016 14:32

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

3 Attachment(s)
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 x-direction. 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      1e-04;          // 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      1e-04;         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      1e-4;         relTol          0.0;     }         UFinal     {         solver          smoothSolver;         smoother        symGaussSeidel;         tolerance      1e-04;         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  1e-4;             relTol      0;         }         p         {             tolerance  5e-3;             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; }```
My boundary condition case files for pressure and velocity are as follows:

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;     }     velocity-inlet_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;     }     velocity-inlet_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;     } }```
For the inlet velocity boundary condition, I used a ramp inlet condition to provide a perturbation in the y-direction after receiving some advice that this perturbation might induce a vortex street if you had trouble with the simulation using a symmetric and structured mesh.

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!

 anishtain4 May 13, 2016 01:25

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.

 kaszt May 17, 2016 15:31

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!

 Kire May 18, 2016 02:05

Hi kaszt

Try making these changes:

1)
Code:

`maxCo  2; // use 1 if 2 doesnt work`

2)
Code:

```solvers {     p        {         solver            GAMG;          tolerance        1e-08;            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        1e-08;         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    1e-8;         relTol          0.1;     }         UFinal     {         solver          smoothSolver;         smoother    symGaussSeidel;         tolerance    1e-08;         relTol          0.0;     } } PIMPLE {     nNonOrthogonalCorrectors 0;     nCorrectors    2;     nOuterCorrectors 50;     pRefCell        0;     pRefValue      0; }     residualControl     {         U         {             tolerance  1e-8;             relTol      0;         }         p         {             tolerance  1e-8;             relTol      0;         }     } relaxationFactors {     fields     {         p      0.7;         pFinal  1;     }     equations     {         U      0.3;         UFinal  1;        } }```

3)
Code:

`    div(phi,U)      Gauss limitedLinearV 1;`

 RodriguezFatz May 20, 2016 04:17

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

 kaszt May 31, 2016 13:57

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 half-way 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 re-meshing 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).

 Kire June 1, 2016 01:40

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 steady-state result by running simpleFoam with/without turbulence off.

I can't remember how long it took me to get the results.

 RodriguezFatz June 1, 2016 03:19

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?

 navidmt November 20, 2018 08:08

Hi all,

I am trying to simulate flow passing over an airfoil with a Reynolds number 60000 , I am modeling it using RANS with SA turbulence model. I am running it using pisoFoam solver and by checking the results, so far it's given good results but the problem is that it is very slow. I wanted to increase the speed of simulations and I read in one of the topics that using pimpleFoam since we can go for higher Courant number is an option. It gave the results very fast but when I check it , there was no vortex shedding and it looked like an average results.

Do you have any idea how I can set the pimpleFoam in order to capture vortex shedding?

 All times are GMT -4. The time now is 15:04.