[AmiN.D]

Dear Cfd-Online!


I want to simulate the vortices around [back ] of a 2D cylinder in openfoam! But i haven't been successful yet!
I should be able to see vortices when 60<Re<200
I put the cylinder in a channel with a size of 30cm*10cm and the diameter of cylinder is 2cm. The inlet velocity is 0.075.


I use picoFoam [I have to use it, because i've been asked to] but after a few steps [say around time=22s] either the Courant number rockets high OR i got Time Step continuity errors.


here is in case of Time Step continuity errors:

Courant Number mean: 0.234093 max: 5.92467
smoothSolver:  Solving for Ux, Initial residual = 0.804164, Final residual = 9.99732e-06, No Iterations 888
smoothSolver:  Solving for Uy, Initial residual = 0.816256, Final residual = 9.97816e-06, No Iterations 533
GAMG:  Solving for p, Initial residual = 0.00591756, Final residual = 0.00050205, No Iterations 5
time step continuity errors : sum local = 0.0182318, global = 0.00146394, cumulative = 0.00126626
GAMG:  Solving for p, Initial residual = 9.51768e-05, Final residual = 8.79025e-07, No Iterations 5
time step continuity errors : sum local = 0.0170765, global = 0.000947447, cumulative = 0.0022137
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib64/libc.so.6
#3  Foam::GaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
#4  Foam::GaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
#5  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
#8  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#9  Foam::SolverPerformance<double> Foam::solve<double>(Foam::tmp<Foam::fvMatrix<double> > const&) at ??:?
#10  Foam::RASModels::kEpsilon<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::correct() at ??:?
#11  ? at ??:?
#12  __libc_start_main in /lib64/libc.so.6
#13  ? at ??:?

if i set the Delta T= 0.000125 the courant number rockets up to 12e+9!




Here is the case summery:

• Mesh via Gmsh , fine around cylinder and a little bit bigger in other zone.
• Not a structured mesh - prisms: 9682
• The checkMesh also shows no error
• Max aspect ratio = 2.32685 OK
  Mesh non-orthogonality Max: 30.1949 average: 4.04066
  Non-orthogonality check OK.
  Max skewness = 0.436476 OK.
• Delta T= 0.00125 and endTime=30s
• Solver for P is GAMG

• I tried to change the Guass linear to Upwind and it didn't help]
• I tried to address high Courant number with decreasing delta T [around 0.000125] but it happens with a delay, but happens again and again.
• I'm using k-epsilon model and the turbulence is on!
• Please be aware that setting maxCo in Controldict is not working with pisoFoam

U boundary :

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.075 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    walls
    {
         type            fixedValue;
        value           uniform (0 0 0);
    }

    

    frontAndback
    {
        type            empty;
    }
    

    cylinder
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

p:

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    frontAndback
    {
        type            empty;
    }

    walls
    {
        type            zeroGradient;
    }


    cylinder
    {
         type            zeroGradient;
    }
}

k

internalField   uniform 7.65e-06;


boundaryField
{
  inlet
      {
               type         zeroGradient;
      }

      outlet
        {
          type            zeroGradient;
        }

   frontAndback
    {
        type            empty;
    }

        walls
        {
          type      zeroGradient;
        }

    
        cylinder
        {
    type   zeroGradient;
      }



}

epsilon:

internalField   uniform 1.72e-07;


boundaryField
{
  inlet
    {
            type         zeroGradient;
    }

       outlet
      {
        type         zeroGradient;
      }


    frontAndback
    {
        type            empty;
    }

      walls
      {
        type    zeroGradient;
      }



      cylinder
      {
         type   zeroGradient;

}
}

nut

internalField   uniform 0;


boundaryField
{
    inlet
    {
       type            zeroGradient;
    }
    outlet
    {       
    type            zeroGradient;
    }

    walls
   
 {
           type            zeroGradient;
  }


    cylinder
    {
             type            zeroGradient;     
    }

 frontAndback
    {
        type            empty;
    }
}

solvers

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.1;
        smoother        GaussSeidel;
    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }

    "(U|k|epsilon|omega|R|nuTilda)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

fvSchemes [openfoam default]

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss limitedLinearV 1;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div(phi,omega)  Gauss limitedLinear 1;
    div(phi,R)      Gauss limitedLinear 1;
    div(R)          Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

So any help would be appreciated !

[ano]

A good initial solution could help you:

1. Use potentialFoam to get the laminar inviscid flow around the cylinder.
2. Then use pisoFoam without a turbulence model.
3. Then switch on the turbulence model.



Increase nNonOrthogonalCorrectors (your mesh around the cylinder is probably unorthogonal).