[parthiv1991]

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:

    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:

    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:

/*--------------------------------*- 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:

/*--------------------------------*- 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:

  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