[kaszt]

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:

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;
        }
    }
}

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;    
    }
}

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:

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;
    }
}

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]

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]

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]

Hi kaszt


Try making these changes:

1)

maxCo   2; // use 1 if 2 doesnt work

2)

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)

    div(phi,U)      Gauss limitedLinearV 1;

[RodriguezFatz]

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.

[kaszt]

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]

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]

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]

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?

[sandrobrad]

Same here,
Did you eventually solve that problem?. I've spent a lot of time only getting 0.001 of my simulation using pisoFoam.
Thanks in advance.