Numerical Analysis of flow over a circular cylinder at different Reynolds number.

I. Aim: Simulate the vortex shedding around a circular cylinder

II. Introduction: Flow over a circular cylinder is an interesting research topic studied for several years. The dynamics of flow over circular cylinder has an application in understanding the heat transfer in cross-flow pipe heat exchangers, the vibration of industrial chimneys due to wind flow, drag forces on underwater pipelines, hissing noise of overhead power conductors (Aeolian vibration) etc. It still remains a challenging problem in fluid mechanics, where intensive investigations are continued even today to understand the complex unsteady dynamics of the cylinder wake flow. Cross-flow normal to the axis of a stationary circular cylinder and the associated problems of heat and mass transport are encountered in a wide variety of engineering applications.

The below figures show the characteristic flow regimes around a circular cylinder for different Reynolds numbers.

Figure 1: Influence of Reynolds number on the flow around a circular cylinder.

Figure 1(a) shows the flow past a cylinder for a very low Reynolds number(~1). The flow smoothly divides and reunites around the cylinder. At a Reynolds number of about 4, the flow (boundary layer) separates downstream of the cylinder and the wake is formed by two symmetric eddies. The eddies remain steady and symmetrical but grow in size up to a Reynolds number of about 40 as shown in Figure 1(b). At a Reynolds number above 40, oscillation in the wake induces asymmetry and finally, the wake starts shedding vortices into the stream. This situation is termed as the onset of periodicity as shown in Figure .1(c) and the wake keeps on undulating up to a Reynolds number of 90.

At a Reynolds number above 90, the eddies are shed alternately from the top and bottom of the cylinder and the regular pattern of alternately shed clockwise and counterclockwise vortices form Von Karman vortex street as in Figure 1(d).  Periodicity is eventually induced in the flow field with the vortex-shedding phenomenon.  Karman investigated the phenomenon and concluded that a non-staggered row of vortices is unstable, and a staggered row is stable only if the ratio of lateral distance between the vortices to their longitudinal distance is 0.28. Because of the similarity of the wake with footprints in a street, the staggered row of vortices behind a bluff body is called a Karman Vortex Street. The vortices move downstream at a speed smaller than the upstream velocity U.

In the range 40 < Re < 80, the vortex street does not interact with the pair of attached vortices. As Re is increased beyond 80 the vortex street forms closer to the cylinder, and the attached eddies themselves begin to oscillate. Finally, the attached eddies periodically break off alternately from the two sides of the cylinder.  While an eddy on one side is shed, that on the other side forms, resulting in an unsteady flow near the cylinder. As vortices of opposite circulations are shed off alternately from the two sides, the circulation around the cylinder changes sign, resulting in an oscillating "lift" or lateral force. If the frequency of vortex shedding is close to the natural frequency of some mode of vibration of the cylinder body, then an appreciable lateral vibration culminates.

The periodicity is characterized by the frequency of vortex shedding `f`.

In non-dimensional form, the vortex shedding frequency is expressed in terms of cylinder diameter 'D'  and upstream velocity 'U' as `S =(f*D)/U`known as the Strouhal number named after V. Strouhal, a German physicist who experimented with wires singing in the wind. The Strouhal number shows a slight but continuous variation with the Reynolds number around a value of 0.21. The boundary layer on the cylinder surface remains laminar and separation takes place at about 81⁰ from the forward stagnation point. At about Re = 500, multiple frequencies start showing up and the wake tends to become Chaotic. As the Reynolds number becomes higher, the boundary layer around the cylinder tends to become turbulent. The wake, of course, shows fully turbulent characters as shown in figure 1 (e).

For larger Reynolds numbers, the boundary layer becomes turbulent. A turbulent boundary layer offers greater resistance to separation than a laminar boundary layer. As a consequence, the separation point moves downstream and the separation angle is delayed to 110⁰ from the forward stagnation point figure 1(f).

Numerical analysis is a viable tool to study the complex flow phenomena around a circular cylinder. The governing equations of continuity, momentum equations can be solved numerically for different Reynolds numbers, to characterize the vortex shedding behind a cylinder.

III Objective:

1. Simulate the flow with the steady and unsteady case and calculate the Strouhal Number for Re= 100. 

2. Calculate the coefficient of drag and lift over a cylinder by setting the Reynolds number to 10,100,1000,10000 & 100000 using a steady-state solver

3. Discuss the effect of Reynolds number on the coefficient of drag. 

IV: Solving & Modelling approach:

2D CFD analysis is carried out in Ansys workbench using Fluent solver to simulate the flow across a circular cylinder for different Reynolds numbers. The cases considered are as below:

Case1: 2D flow over a cylinder with steady and unsteady solver in Ansys Fluent for a Reynolds number, Re =100.

Case 2: 2D flow over a cylinder with steady solver in Ansys Fluent for Re = 10,100,1000,10000,100000.

V: Pre-processing and solver setting:

The first step in a CFD analysis is preprocessing. This includes preparing the geometry/ computational domain, discretization of the domain (meshing), assigning material properties, boundary conditions, solver setup.

1. Geometry/Domain:

The computational domain is prepared using Ansys space-claim design modeler. Since the cylinder has a symmetric cross-section along its length, the flow across it can be approximated as a 2D flow. Hence, the domain is a 2D surface. The domain is generated by sketching a plane surface with a circle in it. A cylinder of unit length with 2m diameter (circle) is considered for the analysis. The details about geometry and fluid domain are given in the below figures:

Geometry	Fluid domain

 

 2. Domain discretization /Meshing:

The 2d fluid domain is discretized into a finite number of control volumes using the Ansys meshing component. A 2D all triangular mesh, with a global element size of 0.25 m is used for the analysis, with the cells near to the cylinder wall are inflated using the inflation option to create boundary layer cells. A total of 6 boundary layer cells are created with a first-layer thickness of 5 mm and a growth rate of 1.2. The cylinder wall is divided into 36 divisions. The generated 2D mesh with a closeup of boundary layer cells is shown below are shown below:

Mesh	Boundary cells near cylinder wall
Global element size = 0.25 m, Cells = 40248 Average Element quality = 0.95

 The boundaries of the domain are also defined in the Ansys mesher using the named selection option. The boundaries are defined as follows:

Boundary/Zone	Boundary/Zone name
Inlet	inlet
Outlet	outlet
Cylinder wall	cylinder_wall
Symmetry	symmetry
Fluid domain	fluid_domain

3. Solver setup:

The 2D analysis is performed on the discretized domain(meshed model) in Ansys Fluent CFD software. A double-precision parallel solver with 4 cores is used. The double-precision solver uses 64-bit precision, which reduces the round-off error in the results.

The above settings are defined at the fluent launcher setup as shown below:

 

 

Once the mesh is opened in fluent, mesh quality is checked using the check mesh, quality -report quality options in the domain tab. 

The properties of the fluid considered for the analysis are given in the below table:

Fluid properties
Fluid	Density - `rho` (kg/m^3)	Dynamic viscosity - `mu`(Pa.s)
User-defined	1	0.02

 

The inlet velocities are varied to get the desired Reynolds number (`Re = (rhoUD)/mu`) as given below :

Fluid inlet velocity (m/s)	Reynolds number `Re`
0.1	10
1	100
10	1000
100	10000
1000	100000

From the above table, it can be inferred that the `Re` is in the laminar region (<10^5 over surfaces ). Hence the laminar viscous model is considered. The pressure-based solver is considered for the analysis. A steady and unsteady state analysis is carried for a Re = 100 as an initial study.

The type of solver (steady/unsteady), operating conditions, viscous model, material properties, zones, and boundary conditions are defined in the respective section of the physics tab as in the below figures.  The operating conditions define the atmosphere in which the analysis is carried out. This sets the reference pressure. In this analysis, a reference pressure of 1 Bar(101325 Pa) is considered to simulate atmospheric conditions. The laminar model is selected in the viscous option of the models section. The properties of the fluid material are defined in the materials section by creating a new user material. Cell zones are defined to apply the materials to the domains. In the current problem, we have only one zone ie. fluid zone with material as 'user-material'.

 

 

The inlet velocity is defined as values corresponding to Reynold's number. The outlet BC is a pressure type with gauze pressure equal to zero, as it is open to the operating condition(atmosphere). These along with wall boundary conditions of no-slip for cylinder wall and symmetry BC are assigned in the boundaries option of the zones section. The below images represent the BC assignment in Fluent.

 

 

Reference Values:

In order to calculate the drag forces on the cylinder, the reference frontal area of the cylinder facing the fluid flow direction along with fluid flow properties are defined in the reference values section as shown below:

 

Solution setup:

After setting up the operating conditions, material properties, zones, and boundary conditions, the solution method, controls, type of reports required, initialization type, and a number of iterations to solve are defined in the solution tab.

The default simple algorithm is used for both steady and unsteady cases. The simple algorithm solves the momentum and continuity equation sequentially to compute velocity and pressure in the flow field. The first-order transient implicit scheme is used for the unsteady case. The under-relaxation factors for better convergence are also defined in control methods.  The discretization schemes chosen for flow variables and gradients are as shown below. Second-order upwind schemes are used for flow variables to improve the solution accuracy.

 

 

The solution is run iteratively in the steady-state case and with respect to time in the unsteady case with parameters as given below:

Steady-state	Unsteady case

 

Reports:

The convergence criteria for the governing equations can be set up using residual monitors in the report section. The reports of drag coefficient, lift coefficient, and velocity at a point in the wake of the cylinder are defined in the reports tab and plotted to monitor them during the solution.

 

 

Initialization:

After setting all the above setup conditions, initialization of the domain has to be done necessary, because the solution is an iterative process that can only be executed if there is an initial guess of the solution variables. standard initialization is employed in the current project, computing the initial values in the domain from the inlet as shown below:

 

VI. Solution Run/ Solving:

After initialization, the iterative solution process is executed using the run calculation section of the solution tab. We have to provide the number of iterations (steady case), timesteps and iterations for time step (unsteady case) to be performed and click calculate to start the iterative solution process.

VII. Postprocessing/ Results extraction:

After the solution process is completed, post-processing has to be done to extract and analyze the results. Post-processing can be done in Fluent as well as Ansys dedicated postprocessor CFD-Post.

VIII. Results and Discussion:

1. Flow over a cylinder for Re = 100

As an initial case, the flow over the cylinder for Re = 100 is analyzed using steady-state and unsteady-state solvers. The steady-state solution is sufficient to compute the lift and drag coefficient on the cylinder and to visualize the vortex shedding phenomenon. But in order to compute the Strouhal frequency which is a time-dependent characteristic of the vortex shedding phenomenon, an unsteady solution is necessary. Hence a comparison is made between steady and unsteady results and also the Strouhal number is computed from the unsteady solution.

Steady-state solution:

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

It can be seen from the above plots, that the residuals, drag coefficient (Cd), become constant and velocity in the wake, lift coefficient (Cl) show a repetitive pattern from around 500 iterations. Hence the solution is said to be converged after 500 iterations. The cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, indicates the vortex shedding phenomenon with the shedding of periodic vortices from the top and bottom side of the cylinder alternatively. The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

It can be seen from the above velocity contour plot, that on the backside of the cylinder a wake is formed with periodic vortices shedding from the top and bottom of the cylinder. The path formed by the periodic vortices is nothing but the Von-Karman vortex street. 

The pressure contour shows the pressure variation on the front side and backside of the cylinder. It can be seen that the pressure increases at the cylinder front end due to stagnation of the fluid. Also, the flow streams separate at the top and bottom of the cylinder creating awake behind the cylinder. This leads to a negative pressure gradient behind the cylinder(blue region) which leads to flow recirculation initiating the vortex development.

The Cd, Cl can be computed from the steady solution, but to compute the Strouhal number a transient/unsteady solution is required.

Unsteady state solution:

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

In a transient solution, the iterations are carried for each timestep in the time marching process. In this analysis, a timestep of 0.1 sec is considered and the solution is run until the steady-state is reached. The number of iterations for each time step is set to 20.  It can be seen from the above plots, that the drag coefficient (Cd), becomes constant and velocity in the wake, lift coefficient (Cl) show a repetitive pattern after around 140 seconds ( 28000 iterations). Hence the solution is said to be converged after 140 seconds. The cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, indicates the vortex shedding phenomenon with the shedding of periodic vortices from the top and bottom side of the cylinder alternatively. The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

The contours are for the time = 160 seconds. The contours are similar to that of a steady-state solution. 

The lift coefficient (Cl ) report file in the time domain is converted into the frequency domain using the (Fast Fourier transform)FFT option in Fluent software to compute the Strouhal number. The Strouhal number value for the flow is the value of non-dimensional frequency at which the Cl is maximum. The FFT plot of Cl is shown below:

 

 

The comparison of the results along with the Strouhal number for the steady and unsteady solution for the flow with Re = 100 is given in the below table.

Re = 100	Cd	Cl	Strouhal Number
Steady State solution	1.338	0.153	-
Unsteady state solution	1.365	0.212	0.162
Rajani et.al	1.335	0.179	0.156

The steady and unsteady results of Cd and Cl are almost identical. The results are in good agreement with the unsteady 2d analysis done by Rajani et.al .  

Further, steady-state analysis is conducted  for the Reynolds numbers 10, 1000, 10000 and 100000 and the results are discussed below:

2. Flow over a cylinder for Re = 10

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

 It can be seen from the above plots, that the drag coefficient (Cd),  velocity in the wake, and lift coefficient (Cl) become constant from around 800 iterations. Hence the solution is said to be converged after 800 iterations. There is no cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, which indicates the vortex shedding phenomenon is not seen in this case (Re = 10), since the boundary layer although separates and forms a wake behind the cylinder, eddies are not formed.  The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

 

3. Flow over a cylinder for Re = 1000

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

 It can be seen from the above plots, that the drag coefficient (Cd),  velocity in the wake, and lift coefficient (Cl) become constant/repetitive from around 400 iterations. Hence the solution is said to be converged after 400 iterations. The cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, indicates the vortex shedding phenomenon with the shedding of periodic vortices from the top and bottom side of the cylinder alternatively. The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

The vortex shedding is more diffused in this case due to higher flow velocity (Re). To have a clear visualization of flow the velocity streamline plotted in CFD post is shown below:

 

 

4. Flow over a cylinder for Re = 10000

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

 It can be seen from the above plots, that the drag coefficient (Cd),  velocity in the wake, and lift coefficient (Cl) become constant/repetitive from around 300 iterations. Hence the solution is said to be converged after 300 iterations. The cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, indicates the vortex shedding phenomenon with the shedding of periodic vortices from the top and bottom side of the cylinder alternatively. The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

The vortex shedding is more diffused in this case due to higher flow velocity (Re). To have a clear visualization of flow the velocity streamline plotted in CFD post is shown below:

 

 

 

5. Flow over a cylinder for Re = 100000

Convergence plots:

Residuals	Velocity at a point in the wake of the cylinder (x = 10m from the cylinder)
Drag Coefficient	Lift Coefficient

 It can be seen from the above plots, that the drag coefficient (Cd),  velocity in the wake, and lift coefficient (Cl) become constant/repetitive from around 500 iterations. Hence the solution is said to be converged after 500 iterations. The cyclic repetitive pattern of velocity in the wake, as well as lift coefficient, indicates the vortex shedding phenomenon with the shedding of periodic vortices from the top and bottom side of the cylinder alternatively. The same can be inferred from the velocity contour plot as shown below. 

Velocity contour :	Pressure contour

The vortex shedding is more diffused in this case due to higher flow velocity (Re). To have a clear visualization of flow the velocity streamlines plotted in the CFD post is shown below:

 

 

 

Validation:

The drag coefficient (Cd) results obtained from the steady-state numerical analysis are compared and validated with the work done by Rajani et.al and Wieselsberger’s (1921) work obtained from work published by Anatol Roshko. The details are given below:

Reynolds Number	Cd, Numerical analysis	Cd,  Rajani et.al	Cd, Wieselsberger’s (1921)	Erorr % Rajani et.al	Erorr %, Wieselsberger’s (1921)
10	2.92	-	2.720	-	7.35
100	1.338	1.335	1.425	0.22	6.10
1000	0.725	-	0.956	-	24.16
10000	0.850	-	1.100	-	22.72
100000	0.952	-	1.2	-	20.66

 

Conclusion:

2D numercial analysis of flow across a circular cylinder is investigated for different Reynolds numbers. There was no vortex shedding observed although wake is formed for Re = 10, for Re = 100, a sustained periodic vortex shedding is observed with forming the Von Karman vortex street. The vortex shedding is more diffused for Re =1000 and there is no sustained vortex street observed. For Re = 10000 and 100000, the flow in the wake becomes unstable and the wake becomes narrower and disorganized. The drag coefficient is higher at Re = 10, thereafter it shows a decreasing trend up to Re = 1000 and increases at a higher Reynolds number.

References:

1. https://nptel.ac.in/content/storage2/courses/112104118/lecture-31/31-3_mechanics.htm

2. https://www.sciencedirect.com/science/article/pii/S0307904X08000243

3. https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/experiments-on-the-flow-past-a-circular-cylinder-at-very-high-reynolds-number/7859A6C46BF4B0F43F11F52AE1C60150