# 2D cylinder simulation

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

April 5, 2024, 18:32
2D cylinder simulation
#1
New Member

aaron LLoyd
Join Date: Apr 2024
Posts: 3
Rep Power: 2
I am simulating an incompressible external flow 2D circular cylinder. I know that it is a “classical” example but I haven`t been able to find a similar case online to improve my results in order to obtain the most liable/precise simulation.
My Reynolds number is 2.4 x10^4 and I should get a Cd = 1.2. Instead I am obtaining a higher value (1,43 in average).
I used snappyHexMesh to generate the mesh and then extrudeMesh, I tried several trials changing layers (number, thickness, conformation) and mesh. In all cases the max value of y+ was <1.
I am attaching the results of a simulation in which there were 40 layers and average y+=0.17 and max y+ = 0.32.
The code below refers to the 0 folder boundary conditions k, nut, omega, p, U
k
Code:
```

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.0294;

boundaryField
{
inlet
{
type            fixedValue;
value           uniform 0.0294;
}
outlet
{
type            inletOutlet;
inletValue      uniform 0.0294;
value           uniform 0.0294;
}
left
{
type            symmetryPlane;
}
right
{
type            symmetryPlane;
}
top
{
type            empty;
}
bottom
{
type            empty;
}
cylinder
{
type            fixedValue;
value           uniform 1e-10;
}
}```
nut
Code:
```

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
inlet
{
type            calculated;
value           uniform 0;
}
outlet
{
type            calculated;
value           uniform 0;
}
left
{
type            symmetryPlane;
}
right
{
type            symmetryPlane;
}
top
{
type            empty;
}
bottom
{
type            empty;
}
cylinder
{
type            calculated;
value           uniform 0;
}
}```
omega
Code:
```

dimensions      [0 0 -1 0 0 0 0];

internalField   uniform 164;

boundaryField
{
inlet
{
type            fixedValue;
value           uniform 164;
}
outlet
{
type            inletOutlet;
inletValue      uniform 164;
value           uniform 164;
}
left
{
type            symmetryPlane;
}
right
{
type            symmetryPlane;
}
top
{
type            empty;
}
bottom
{
type            empty;
}
cylinder
{
type            omegaWallFunction;
value           uniform 164;
}
}```

p
Code:
```

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
inlet
{
}
outlet
{
type            fixedValue;
value           uniform 0;
}
left
{
type            symmetryPlane;
}
right
{
type            symmetryPlane;
}
top
{
type            empty;
}
bottom
{
type            empty;
}
cylinder
{
}
}```
and U

Code:
```

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (14 0 0);

boundaryField
{
inlet
{
type            fixedValue;
value           uniform (14 0 0);
}
outlet
{
type            inletOutlet;
inletValue      uniform (0 0 0);
value           uniform (14 0 0);
}
left
{
type            symmetryPlane;
}
right
{
type            symmetryPlane;
}
top
{
type            empty;
}
bottom
{
type            empty;
}
cylinder
{
type            fixedValue;
value           uniform (0 0 0);
}
}```

and controlDict
Code:
``` application     pimpleFoam;

startFrom       latestTime;

stopAt          endTime;

endTime         100;

deltaT          0.00001;

writeInterval   0.001;

purgeWrite      1;

writeFormat     binary;

writePrecision  7;

writeCompression no;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

maxCo           0.9;

functions
{
#include "forces"
#include "forceCoeffs"
#include "yPlus1"
#include "pressureCoefficient"
#include "wallShearStress"
#include "Q"
#include "wallBoundedStreamLine"
}```
and fvSchemes

Code:
``` ddtSchemes
{
default         CrankNicolson 0.9;
}

{
default         Gauss linear;

}

divSchemes
{
default         none;

div(phi,U)      Gauss linear;
div(phi,k)      Gauss upwind;
div(phi,omega)  Gauss upwind;

}

laplacianSchemes
{
default         Gauss linear limited 1.0;
}

interpolationSchemes
{
default         linear;
}

{
default         limited 1.0;
}

wallDist
{
method          meshWave;
}```
and fvSolution

Code:
```  solvers
{
p
{
solver          GAMG;
smoother        GaussSeidel;
tolerance       1e-6;
relTol          0.01;
}

pFinal
{
solver          GAMG;
smoother        GaussSeidel;
tolerance       1e-8;
relTol          0;
}

U
{
solver          smoothSolver;
smoother        symGaussSeidel;
tolerance       1e-8;
relTol          0.1;
}

UFinal
{
\$U;
relTol      0;
}

k
{
solver          smoothSolver;
smoother        symGaussSeidel;
tolerance       1e-8;
relTol          0.1;
}

kFinal
{
\$U;
relTol      0;
}

omega
{
solver          smoothSolver;
smoother        symGaussSeidel;
tolerance       1e-8;
relTol          0.1;
}

omegaFinal
{
\$U;
relTol      0;
}

}

PIMPLE
{

nCorrectors              2;

nNonOrthogonalCorrectors 1;

nOuterCorrectors    1;

}```
I am attaching a screenshot of the mesh and the images of the results: RESIDUAL(Ux, Uy, p) and Cd. The chart shows the final period as the simulation appears already for a while convergent.
Attached Images
 Cd.jpg (30.7 KB, 7 views) Residual.jpg (60.7 KB, 11 views) mesh.png (169.0 KB, 9 views)

 April 6, 2024, 06:21 #2 New Member   ROIN bin Join Date: Dec 2019 Posts: 11 Rep Power: 6 It looks very close, I think we can adjust the calculation area and boundary conditions (such as outlet), and I noticed that you used the built-in post-processing tool to output the coefficient, perhaps test it with force.

 April 6, 2024, 11:36 #3 New Member   aaron LLoyd Join Date: Apr 2024 Posts: 3 Rep Power: 2 Thanks for your reply and help: let me try to answer. In forceCoeffs lRef= 0.0272 which corresponds to the diameter of the cylinder, and Aref= 0.00272 since I set thikness= 0.1 in extrudeMesh. In the postProcessing folder, taking the Cd value at time 0.4, as an example, I have Cd= 1.5499957, and taking at the same time the value of the Total force = Pressure + Viscous in the x direction = 5.0612938e-01. F= 1/2*1.225*14^2*0.00272*1.5499957 which gives a slightly different value: 0.506116759, but this difference seems negligible to me . As for the outlet boundary conditions I just tried to replace in k and omega "type inletOutlet" with "type zeroGradient", but I haven`t seen any change.

 April 10, 2024, 09:21 #4 New Member   aaron LLoyd Join Date: Apr 2024 Posts: 3 Rep Power: 2 Hi, can anyone give me an advice? I would really appreciate it

 Tags cylinder, mesh, openfoam, simple