CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Simulating Flow past Circular Cylinder

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

Like Tree1Likes
  • 1 Post By piu58

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 19, 2018, 12:25
Default Simulating Flow past Circular Cylinder
  #1
New Member
 
PLD
Join Date: Jun 2017
Location: Braunschweig, Germany
Posts: 13
Rep Power: 5
parthiv1991 is on a distinguished road
Deal All,

I have been trying to simulate the Schäfer and Turek circular cylinder case in OpenFOAM using two custom solvers using higher order Runge-Kutta time integration schemes: 1. DIRK 2. ROSW methods.

So, I have simulated the Steady-state 2D cylinder case (stated as test case 2D-1 in Schäfer and Turek's paper). Inlet U = 0.3 m/s, and Re = 20 (Laminar flow)

I have used two different meshes (as seen below).
The problem is that while computing the force Coefficients, the drag value obtained is within the bounds mentioned in the paper, but the lift coeff value is below the bounds. Essentially, I am getting correct drag coeff. values but incorrect lift coeff. values.
I cannot determine from where the error might generate (due to the mesh, or Aref & Lref values, or the various schemes).

The coarse mesh: No. of cells = 4800, No. of points = 9980
anf the checkMesh for this is as below:
Code:
    Overall domain bounding box (0 0 -0.0005) (2.2 0.41 0.0005)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-5.28852e-20 1.11437e-19 -3.76119e-16) OK.
    Max cell openness = 2.11208e-16 OK.
    Max aspect ratio = 4.5 OK.
    Minimum face area = 3.86507e-06. Maximum face area = 0.000495.  Face area magnitudes OK.
    Min volume = 2.50828e-08. Max volume = 4.95e-07.  Total volume = 0.000894154.  Cell volumes OK.
    Mesh non-orthogonality Max: 32.295 average: 6.50871
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.576444 OK.
    Coupled point location match (average 0) OK.
Fine mesh: No. of cells = 22620, No. of points = 46056
and the checkMesh for this shows:
Code:
    Overall domain bounding box (0 0 -0.0005) (2.2 0.41 0.0005)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-6.66972e-20 5.66631e-20 5.41013e-15) OK.
    Max cell openness = 2.33164e-16 OK.
    Max aspect ratio = 4.98629 OK.
    Minimum face area = 1.3752e-06. Maximum face area = 9.13442e-05.  Face area magnitudes OK.
    Min volume = 2.16552e-09. Max volume = 9.13442e-08.  Total volume = 0.000894147.  Cell volumes OK.
    Mesh non-orthogonality Max: 33.1296 average: 5.87691
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.853518 OK.
    Coupled point location match (average 0) OK.

Mesh OK.
The fvSolution file looks like this:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }
    kihelp
     {
      solver          GAMG;
      smoother        DILU;
      tolerance       1e-06;
     }  


     slope
     {
      solver          PBiCGStab;
      smoother        DILU;
      tolerance       1e-06;
     }


    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }
}

PISO
{
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}
And the fvSchemes file:

Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         leastSquares;
    grad(p)         leastSquares;
}

divSchemes
{
    default         Gauss linear;
    div(phi,U)      Gauss linearUpwind cellLimited Gauss linear 1;
    div(U)          Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}
My forceCoeffs code snippet in controlDict look like this:

Code:
  forces
    {
       type    forces;
       functionObjectLibs ("libforces.so");

       patches ( cylinder );
       pName p;
       UName U;

       log   true;
       CofR (0.2 0.2 0);
       rho   rhoInf;
       rhoInf 1;

       outputControl timeStep;
       outputInterval 10;

    } 

    forceCoefficients
    {
       type    forceCoeffs;
       functionObjectLibs ("libforces.so");
       log yes;

       outputControl timeStep;
       outputInterval 10;

       patches ( cylinder );
       pName p;
       UName U;
       rho   rhoInf;
       rhoInf 1;
    //   origin (0 0 0);

     //  porosity   no;

       liftDir (0 -1 0);
       dragDir (1 0 0);
       CofR    (0.2 0.2 0);  // Centre of cylinder
       pitchAxis (0 0 1);

       magUInf 0.2;  // 2*U/3 , where U = 0.3 m/s (inflow velocity)
       rhoInf 1;
       lRef 0.1;  // Dia of cylinder 0.1m
       Aref 0.0001; // (Dia = 0.1) * (z-axis cylinder length = 0.001m)
And I have attached the Lift and Drag Coeff plots below for the two meshes.
Acc. to Schäfer and Turek, the bound for Drag coeff : [5.57, 5.59] and Lift: [0.0104, 0.0110]

Could anyone suggest where the error might come from??
Also I have notices that a finer mesh gives me a lower value of Lift Coeff (which is far lower than the bound, than for a coarse mesh which gives a Lift coeff value a bit closer to the bound).

Any help would be appreciated.

Thank You
Attached Images
File Type: jpg mesh.jpg (114.1 KB, 37 views)
File Type: jpg mesh_fine.jpg (110.1 KB, 32 views)
File Type: jpg LIFT.jpg (34.6 KB, 21 views)
File Type: jpg DRAG.jpg (26.5 KB, 20 views)
parthiv1991 is offline   Reply With Quote

Old   February 19, 2018, 12:54
Default
  #2
Senior Member
 
piu58's Avatar
 
Uwe Pilz
Join Date: Feb 2017
Location: Leipzig, Germany
Posts: 588
Rep Power: 9
piu58 is on a distinguished road
I would avoid the sudden changes in mesh density. It is far better to change the mesh gradually. blockMesh gives a good support for that.
parthiv1991 likes this.
__________________
Uwe Pilz
--
Die der Hauptbewegung überlagerte Schwankungsbewegung ist in ihren Einzelheiten so hoffnungslos kompliziert, daß ihre theoretische Berechnung aussichtslos erscheint. (Hermann Schlichting, 1950)
piu58 is offline   Reply With Quote

Old   February 19, 2018, 13:41
Default
  #3
Senior Member
 
Santiago Lopez Castano
Join Date: Nov 2012
Posts: 312
Rep Power: 12
Santiago is on a distinguished road
Quote:
Originally Posted by parthiv1991 View Post
Deal All,

I have been trying to simulate the Schäfer and Turek circular cylinder case in OpenFOAM using two custom solvers using higher order Runge-Kutta time integration schemes: 1. DIRK 2. ROSW methods.

So, I have simulated the Steady-state 2D cylinder case (stated as test case 2D-1 in Schäfer and Turek's paper). Inlet U = 0.3 m/s, and Re = 20 (Laminar flow)

I have used two different meshes (as seen below).
The problem is that while computing the force Coefficients, the drag value obtained is within the bounds mentioned in the paper, but the lift coeff value is below the bounds. Essentially, I am getting correct drag coeff. values but incorrect lift coeff. values.
I cannot determine from where the error might generate (due to the mesh, or Aref & Lref values, or the various schemes).

The coarse mesh: No. of cells = 4800, No. of points = 9980
anf the checkMesh for this is as below:
Code:
    Overall domain bounding box (0 0 -0.0005) (2.2 0.41 0.0005)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-5.28852e-20 1.11437e-19 -3.76119e-16) OK.
    Max cell openness = 2.11208e-16 OK.
    Max aspect ratio = 4.5 OK.
    Minimum face area = 3.86507e-06. Maximum face area = 0.000495.  Face area magnitudes OK.
    Min volume = 2.50828e-08. Max volume = 4.95e-07.  Total volume = 0.000894154.  Cell volumes OK.
    Mesh non-orthogonality Max: 32.295 average: 6.50871
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.576444 OK.
    Coupled point location match (average 0) OK.
Fine mesh: No. of cells = 22620, No. of points = 46056
and the checkMesh for this shows:
Code:
    Overall domain bounding box (0 0 -0.0005) (2.2 0.41 0.0005)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (-6.66972e-20 5.66631e-20 5.41013e-15) OK.
    Max cell openness = 2.33164e-16 OK.
    Max aspect ratio = 4.98629 OK.
    Minimum face area = 1.3752e-06. Maximum face area = 9.13442e-05.  Face area magnitudes OK.
    Min volume = 2.16552e-09. Max volume = 9.13442e-08.  Total volume = 0.000894147.  Cell volumes OK.
    Mesh non-orthogonality Max: 33.1296 average: 5.87691
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.853518 OK.
    Coupled point location match (average 0) OK.

Mesh OK.
The fvSolution file looks like this:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }
    kihelp
     {
      solver          GAMG;
      smoother        DILU;
      tolerance       1e-06;
     }  


     slope
     {
      solver          PBiCGStab;
      smoother        DILU;
      tolerance       1e-06;
     }


    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }
}

PISO
{
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}
And the fvSchemes file:

Code:
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         leastSquares;
    grad(p)         leastSquares;
}

divSchemes
{
    default         Gauss linear;
    div(phi,U)      Gauss linearUpwind cellLimited Gauss linear 1;
    div(U)          Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}
My forceCoeffs code snippet in controlDict look like this:

Code:
  forces
    {
       type    forces;
       functionObjectLibs ("libforces.so");

       patches ( cylinder );
       pName p;
       UName U;

       log   true;
       CofR (0.2 0.2 0);
       rho   rhoInf;
       rhoInf 1;

       outputControl timeStep;
       outputInterval 10;

    } 

    forceCoefficients
    {
       type    forceCoeffs;
       functionObjectLibs ("libforces.so");
       log yes;

       outputControl timeStep;
       outputInterval 10;

       patches ( cylinder );
       pName p;
       UName U;
       rho   rhoInf;
       rhoInf 1;
    //   origin (0 0 0);

     //  porosity   no;

       liftDir (0 -1 0);
       dragDir (1 0 0);
       CofR    (0.2 0.2 0);  // Centre of cylinder
       pitchAxis (0 0 1);

       magUInf 0.2;  // 2*U/3 , where U = 0.3 m/s (inflow velocity)
       rhoInf 1;
       lRef 0.1;  // Dia of cylinder 0.1m
       Aref 0.0001; // (Dia = 0.1) * (z-axis cylinder length = 0.001m)
And I have attached the Lift and Drag Coeff plots below for the two meshes.
Acc. to Schäfer and Turek, the bound for Drag coeff : [5.57, 5.59] and Lift: [0.0104, 0.0110]

Could anyone suggest where the error might come from??
Also I have notices that a finer mesh gives me a lower value of Lift Coeff (which is far lower than the bound, than for a coarse mesh which gives a Lift coeff value a bit closer to the bound).

Any help would be appreciated.

Thank You
I have tested a similar case, that of Muzaferija et al. (1995) with Re=100, using a variant of PISO. I got decent results, within bounds, and a relative error on Cd & Cf of 3% on the Richardsons' extrapolated values when confronted with reference ones. One question: you said you use RK ad time scheme, but I see Euler as your ddt, could it be the problem?

Why you use such a scheme for div(phi,U)? Gauss linear is enough.

What is your maxCo?

Finally, is there any reference on your particular code? Is difficult to know if the algorithm is not clear.
Santiago is offline   Reply With Quote

Old   February 20, 2018, 08:14
Default
  #4
New Member
 
PLD
Join Date: Jun 2017
Location: Braunschweig, Germany
Posts: 13
Rep Power: 5
parthiv1991 is on a distinguished road
Hello Santiago,

Thank you for your reply.

My max Courant number is 0.181081

Yes I am using a custom RK scheme for time integration which I have put in the code, which is a modified version of icoFoam solver. So my entire time integration is done inside the solver itself. Then should I set the ddt scheme as "none" and run? I will try this once and let you know.

But I think that it would give me an error because for example, the initial guess for the velocity slope "k" I have chosen as ddt(U). Then I guess I should not use this ddt(U) if I have to turn off the ddt schemes in fvSchemes.

But I don't think that the ddtScheme is causing the problem. But I will try your suggestion.

Thank you
parthiv1991 is offline   Reply With Quote

Old   June 3, 2018, 06:19
Default
  #5
Member
 
Johan Lorentzon
Join Date: Mar 2009
Location: Lunds University, Sweden
Posts: 76
Rep Power: 19
pi06jl6 will become famous soon enough
Greetings,

I tried Turek's case on OpenFOAM a few years ago. It works fine for FSI1 but not FSI2 and FSI3. I would appreciate for a solution, it doesn't work on for at least OF2.3.X

The problem is the following: Mesh motion solver and the boundary condition, you need to create your own point field distribution and move the whole field of points on the walls that define your 3D case, since you run the case in 2D frozen mode.

This is solved however for extended Foam, using their tetFem mesh motion solver. But I have not adapted my FSI solver to that distribution.

If you need to validate with using 3D mesh, I do recommend to use the 3D pipe problem, used by many, for example, Degroote/IQN-ILS.

I wish you best luck with the problem

/J

PS. I am still curious if there is a solution to this problem. However, it is of minor importance for my work since I have validated my solver on other cases.

Last edited by pi06jl6; June 5, 2018 at 15:21.
pi06jl6 is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Discrepancy in the Strouhal number of a flow past circular cylinder HectorRedal Main CFD Forum 13 April 6, 2017 19:20
flow around a circular cylinder with velocity inlet and outflow outlet shuoxue OpenFOAM 1 March 3, 2014 11:42
flow around a circular cylinder with velocity inlet and outflow outlet shuoxue OpenFOAM Running, Solving & CFD 0 November 2, 2013 05:32
3D FLOW OVER A CIRCULAR CYLINDER Srinivas Mettu FLUENT 2 April 4, 2010 23:11
Simulation of the Flow past a circular cylinder using STAR-CD M. S. GUEROUACHE Main CFD Forum 0 October 1, 1998 11:51


All times are GMT -4. The time now is 14:59.