This paper presents two- and three-dimensional direct numerical simulations of the flow around a circular cylinder placed symmetrically in a plane channel. Results are presented in the Reynolds number range (based on the cylinder diameter and centerline velocity) of 10 to 390 for a blockage ratio (ratio of the cylinder diameter to the channel height) of 0.2. The aim of this work was to investigate in detail the confinement effect due to the channel’s stationary walls on the force coefficients and the associated Strouhal numbers, as well as on the generated flow regimes. Present results suggest a transition from a 2-D to a 3-D shedding flow regime between *Re* = 180 and *Re* = 210. This transition was found to be dominated by mode A and mode B three dimensional instabilities, similar to those observed in the case of an unconfined circular cylinder. This is the first time that the existence of the two modes, and of naturally occurring vortex dislocations, has been confirmed via full 3-D simulations for the case of a confined circular cylinder in a channel. A discontinuity in the variation of the Strouhal number *St*, and of the base pressure coefficient *C*_{pb}, with *Re* was also observed. This was found to be associated with the onset of mode A instability and the development of vortex dislocations, and parallels what occurs in the unconfined case, but previous studies could not confirm its existence in the confined case. Furthermore, by analyzing the mechanisms affecting the shape and evolution of these instabilities, it is demonstrated that they are significantly affected by the confinement only in the far wake.

## I. INTRODUCTION

While the flow over a circular cylinder represents one of the classical problems in fluid mechanics, the case of flow over a confined cylinder in a plane channel remains relatively unexplored. The extra confinement provided by the stationary no-slip walls of the channel affects the nature and stability of the flow. Understanding the wide variety of rich flow phenomena that ensue in this case is of intrinsic interest for the overall understanding of bluff body fluid dynamics. Even more importantly, such a flow configuration represents an idealization of several industrially important flows, where flow inserts can be used to enhance mixing and heat transfer; typical examples include flow past dividers in polymer processing, turbulence promoters in the liquid-metal blankets of fusion reactors, etc. Understanding the dynamics of a three-dimensional wake flow behind a cylinder can provide valuable knowledge with practical importance, with respect to its effect on heat and mass transfer.

Transition to three-dimensionality in the cylinder wake for the *unconfined* case is well understood and a lot of work has already been done.^{1–13} For this case, the flow is known to remain purely two-dimensional for values of the Reynolds number up to $Re\u2248189$.^{9} At even higher Re values, three-dimensional effects appear, and Roshko^{1} was the first to observe a transition regime in the wake of an unconfined circular cylinder. In this transition regime, three physically different instabilities are observed, referred to by Williamson^{7} as “mode A,” “mode B,” and “vortex dislocations,” Mode A is characterized by a discontinuous change in the wake formation, as the primary spanwise vortices become unstable and generate large-scale streamwise vortex loops, at a wavelength of around 3 to 4 cylinder diameters. Mode B corresponds to the appearance of small scale streamwise vortex structures, with a wavelength of approximately one cylinder diameter. Finally, vortex dislocations are generated between spanwise cells of different frequencies and evolve along the span. Williamson^{4} suggests that this phenomenon is mostly responsible for the low frequency fluctuations reported by Roshko^{1} to characterize the transition regime and the appearance of turbulent motion. Wake transition is defined by two discontinuities in the Strouhal frequency. At the first discontinuity, the Strouhal frequency drops from the laminar curve to one corresponding to a mode A three-dimensional shedding. As the Reynolds number is increased further, another discontinuity is observed related to mode B instabilities. This discontinuity is not hysteretic. Experiments have shown that the transition from mode A to mode B shedding takes place in the *Re* range 190 to 260. According to Williamson,^{7} wake transition is further characterized by velocity and pressure modifications. The appearance of mode-A instability and large scale dislocations are consistently accompanied by a reduction of base suction, a reduction of (two-dimensional) Reynolds stress level, and a growth in the size of the formation region.

In the case of a circular cylinder *confined* between parallel walls, studies investigating three dimensional effects at low Reynolds numbers, either experimentally or numerically, are relatively scarce. Almost all three-dimensional studies are at reasonably higher Reynolds numbers, in the range *O*(10^{4})–*O*(10^{5}).^{14–17} These studies investigate the flow characteristics in the wake of a circular cylinder placed near a plane wall at various gap heights, G (distance between the cylinder and the wall). There results show that for small gap ratios, *G*/*D* ≤ 0.3, vortex shedding is suppressed. For larger gap ratios, some flow quantities, such as the Strouhal number, remain remarkably constant and independent of the gap ratio. Moreover, Wang and Tan^{17} found that for the intermediate gap ratios 0.3 ≤ *G*/*D* ≤ 0.6, the influence of the wall is fairly strong, causing the flow to develop a distinct asymmetry about the cylinder centreline.

Only two studies, quite recently, focused on the development of three-dimensional effects for a confined circular cylinder. Rehimi *et al.*^{18} carried out experimental studies in the Reynolds number range 30–277 for a blockage ratio of 1/3. In their work, three-dimensional instabilities were observed having great similarities with mode A and mode B patterns of vortex shedding found in the unconfined case. However, in their experiments mode A appeared at *Re* = 159, which is considered premature relative to the critical value of $Re\u2248189$ for the onset of this instability in the unconfined case. They attributed this behavior to wall perturbations. The second study was that of Camarri and Giannetti^{19} where a three-dimensional Floquet stability analysis was used to investigate the three-dimensional stability of the wake behind a symmetrically confined circular cylinder for a different blockage ratio of 1/5 for *Re* up to 300. Their results showed that the transition to a three-dimensional state has the same space-time symmetries as the unconfined case. They also found that the critical Reynolds numbers for the onset of instabilities, when based on the centerline velocity, are similar to the ones found in Barkley and Henderson^{9} for the unconfined case. However, they observed that due to the inversion of the wake vortices, the linear unstable modes are significantly affected, leading to differences between the confined and the unconfined cases.

To the best of our knowledge, no published direct numerical simulations of the three-dimensional flow in a channel exists in the Reynolds number range examined here (150 to 390 based on the centerline velocity). The main objective of this study is to investigate how the confinement provided by the channel walls affects the onset and development of three-dimensional instabilities in the wake of the flow at moderate Reynolds numbers. Comparisons are made with two-dimensional studies for the confined case, as well as with the extensive literature available for the unconfined case.

In the sections to follow, first the complete problem is formulated and presented together with a detailed grid-sensitivity analysis (Sec. II). The main results from the 3-D simulations together with a discussion of confinement effects are presented in Sec. III, followed by the conclusions of the present study.

## II. PROBLEM STATEMENT AND FORMULATION

### A. Flow configuration

The geometry considered in this study is shown in Figure 1. The geometry consists of a circular cylinder, of diameter *D*, symmetrically placed in a plane channel. The ratio of the cylinder diameter to the distance between the channel walls *H*, defines the blockage ratio, β = *D*/*H*. Results for β = 1/5 are presented in the present paper, allowing for a direct comparison with the linear stability analysis of Camarri and Gianneti.^{19} We decided not to use the higher blockage ratio (β = 1/3) that was used in the experiments of Rehimi *et al.*,^{18} even though that would have meant a smaller computational domain in favour of computational cost, “in order to avoid peculiar flow features related to a complex interaction between the wake and the confining walls,” as mentioned in Camarri and Giannetti.^{19} By choosing the lower blockage ratio, we can bring into focus the wall-blocking effects without these being obscured by the more complex interactions with the near-wall viscous regions.

The channel inlet is placed at a distance of *L*_{i} = 12.5*D* upstream of the circular cylinder, while the outlet is located at *L*_{o} = 35.5*D* behind the body. This choice of parameters ensures minimal distortion of the flow structure due to the boundary conditions,^{20} while maintaining a reasonable computational cost, and is in line with those used by others in the literature.^{19,21} Two spanwise lengths were considered, *W* = 8*D* and *W* = 12*D*. Based on observations from the transitional wake of open uniform flow past a circular cylinder, the value of 8*D* was considered sufficient for the development of both mode B and the larger mode A three-dimensional wake instabilities. The relatively large spanwise length of 12*D* has been used only in the analysis of natural vortex dislocations. In this case, the use of such a large span was motivated by the need to provide sufficient domain for the development of such irregularities in the wake, and it is in agreement with the practice followed by Braza *et al.*^{12}

### B. Mathematical formulation

The flow is completely described by the set of Navier-Stokes equations for an incompressible Newtonian fluid, of density ρ, dynamic viscosity μ, and kinematic viscosity ν = μ/ρ. Using the cylinder’s diameter, *D*, and the centerline inflow velocity, *U*_{c}, as the characteristic length and velocity scales, respectively, the non-dimensional continuity and momentum equations in a Cartesian coordinate system are given by

where *Re* = *U _{c}D*/ν is the Reynolds number.

The drag and lift coefficients were both determined by considering the viscous and the pressure forces on the cylinder surface,

where *F*_{D} and *F*_{L} are the drag and lift forces per unit length of the cylinder, defined as

Here, *n*_{j} denotes the unit normal vector pointing in the direction of *x*_{j}. The Strouhal number, characterizing the vortex shedding phenomenon, is based on the dominant frequency of the lift coefficient

The base pressure coefficient is defined as

where *p*_{b} is the spatiotemporal average, with respect to the time and the spanwise coordinate *z*, of pressure at the rear stagnation point of the cylinder (180° from the front), and *p*_{∞} is the free-stream pressure at the inlet boundary (*x* = −12.5*D*, *y* = 0).

### C. Numerical method and boundary conditions

The current computations have been performed using an unstructured collocated nodal-based finite-volume code. A second-order accurate centered-difference scheme, with skewness corrections, is applied to discretize the diffusive and nonlinear terms. The mass and momentum equations [Eqs. (1) and (2)] are coupled using a fractional-step method. The time integration of the flow equations is done using a Crank-Nicholson scheme for the diffusive terms. The non-linear terms are treated semi-implicitly. A very detailed description of the numerical techniques used by this code is reported in Ham *et al.*,^{22} and Moin *et al.*^{23}

Boundary conditions are given at the domain’s inlet by prescribing a Poiseuille, parabolic velocity profile,

while a no-slip boundary condition is imposed on the cylinder surface, the top, and the bottom walls. At the outlet, in order to minimize reflective effects and avoid the distortion of the flow structures leaving the domain, a convective boundary condition is applied,

where $Uconv\u2261\u222buds/\u222bds$. For the 3-D simulations, a Neumann boundary condition has been adopted for the velocity field in the spanwise direction,

which is considered appropriate and in accordance with experimental observations of the flow regime considered in the present study.^{11} A series of simulations performed with periodic spanwise boundary conditions suggested that, even though numerical results where similar, not all of the critical physical mechanisms could be captured as accurately as with Neumann conditions.

### D. Computational grids

The influence of the grid resolution on the computed physical characteristics of the flow was examined in detail in order to optimize the simulation in terms of accuracy and computational cost. For this reason, a series of two- and three-dimensional simulations were carried out at *Re* = 300. The parameters of the different grid configurations tested are summarized in Table I.

Grids . | . | N_{cyl}
. | Δn_{walls}
. | Δz
. | N_{z}
. | W
. | Total nodes . |
---|---|---|---|---|---|---|---|

G1 | 2-D | 120 | 0.020 | — | — | — | 25 531 |

G2 | 2-D | 240 | 0.005 | — | — | — | 66 989 |

G3 | 2-D | 320 | 0.002 | — | — | — | 112 657 |

G2C | 3-D | 240 | 0.005 | 0.4 | 21 | 8 | 1 406 769 |

G2M | 3-D | 240 | 0.005 | 0.2 | 41 | 8 | 2 726 109 |

G2F | 3-D | 240 | 0.005 | 0.1 | 81 | 8 | 5 426 109 |

G2W | 3-D | 240 | 0.005 | 0.12 | 101 | 12 | 6 765 889 |

Grids . | . | N_{cyl}
. | Δn_{walls}
. | Δz
. | N_{z}
. | W
. | Total nodes . |
---|---|---|---|---|---|---|---|

G1 | 2-D | 120 | 0.020 | — | — | — | 25 531 |

G2 | 2-D | 240 | 0.005 | — | — | — | 66 989 |

G3 | 2-D | 320 | 0.002 | — | — | — | 112 657 |

G2C | 3-D | 240 | 0.005 | 0.4 | 21 | 8 | 1 406 769 |

G2M | 3-D | 240 | 0.005 | 0.2 | 41 | 8 | 2 726 109 |

G2F | 3-D | 240 | 0.005 | 0.1 | 81 | 8 | 5 426 109 |

G2W | 3-D | 240 | 0.005 | 0.12 | 101 | 12 | 6 765 889 |

First, a series of two-dimensional simulations were performed using three non-uniform grids (G1, G2, G3), mainly differing in the spatial resolution in the vicinity of the cylinder and the channel walls. The hyperbolic tangent function was used for stretching the cell sizes in a clustered region close to the cylinder, and linear grid stretching was applied in the direction normal to the channel walls, as shown in Figure 2. Further upstream and downstream of the cylinder, a uniform grid was used. The discrepancy between the values of mean drag coefficient $C\u2212$ _{D}, rms value of the lift coefficient $C'L$, Strouhal number *St*, and the base pressure coefficient *C*_{pb} resulting from the use of different grids are shown in Table II. Results show that a grid independent solution can be achieved with the grids considered. For example, the percentage difference between the values predicted on the coarsest grid with respect to the ones obtained on the finest grid is below 2.1%. This discrepancy is further reduced to less than 0.5% when the two finest grids are compared. Hence, one can conclude that the intermediate grid G2 is sufficiently fine to resolve the flow.

Grids . | . | Re
. | $C\u2212D$ . | $C'L$ . | St
. | C_{pb}
. |
---|---|---|---|---|---|---|

G1 | 2-D | 300 | 1.210 | 0.551 | 0.2021 | 1.0440 |

G2 | 2-D | 300 | 1.232 | 0.559 | 0.2023 | 1.0472 |

G3 | 2-D | 300 | 1.236 | 0.562 | 0.2024 | 1.0488 |

G2C | 3-D | 300 | 1.157 | 0.405 | 0.1915 | 0.929 |

G2M | 3-D | 300 | 1.167 | 0.388 | 0.1967 | 0.937 |

G2F | 3-D | 300 | 1.172 | 0.376 | 0.1989 | 0.942 |

Grids . | . | Re
. | $C\u2212D$ . | $C'L$ . | St
. | C_{pb}
. |
---|---|---|---|---|---|---|

G1 | 2-D | 300 | 1.210 | 0.551 | 0.2021 | 1.0440 |

G2 | 2-D | 300 | 1.232 | 0.559 | 0.2023 | 1.0472 |

G3 | 2-D | 300 | 1.236 | 0.562 | 0.2024 | 1.0488 |

G2C | 3-D | 300 | 1.157 | 0.405 | 0.1915 | 0.929 |

G2M | 3-D | 300 | 1.167 | 0.388 | 0.1967 | 0.937 |

G2F | 3-D | 300 | 1.172 | 0.376 | 0.1989 | 0.942 |

To study the effect of the spanwise grid spacing in the numerical solution for the 3-D cases, several three-dimensional grids were generated by repeating the grid G2 along the spanwise direction. Three grids differing in the spanwise resolution were tested, namely G2C, G2M, and G2F with Δ*z* grid spacings of 0.4*D*, 0.2*D*, and 0.1*D*, respectively. The spanwise dimension of 8*D* was adopted for these meshes. Simulations were performed over periods of at least 900 dimensionless time units, corresponding to about 180 vortex shedding cycles, with time averaging of results performed over the last 100 shedding cycles, when the flow had reached a “fully developed” state.

Results obtained with the three grids are listed in Table II and overall show good grid convergence. For example, replacing the coarsest grid G2C with the finest grid G2F, resulted in only 3.7% change in *St* and 1.4% in *C*_{pb}, whereas the corresponding changes were only 1.1% and 0.6% when replacing the intermediate grid G2M by the finest grid. Based on these results, one could have opted for the intermediate grid G2M. However, because we were particularly interested in capturing the details of the vortical structures downstream of the cylinder, we decided to use the finest grid G2F, despite the increased computational cost. For example, a typical simulation with grid G2F required a total of 27 days of computation on a 32 node (64 processor) Linux cluster to complete 180 shedding cycles. Each node has dual Opteron 244 (1.8 GHz) processors with 4 GB of RAM.

In order to study the natural occurrence of vortex dislocations, a case at *Re* = 240 was computed in a wider domain with a spanwise dimension of 12*D*. To accommodate this case, an additional grid (G2W) was generated, which in view of the increased computational requirements had a spanwise grid spacing of 0.12*D*, which is slightly coarser than G2F, but still finer than G2M.

The dimensionless time step was kept constant during each simulation but was determined independently for each case in order to satisfy the Courand-Friedrichs-Lewy stability criterion, *CFL* ≤ 1. This yielded values in the range 7.5 × 10^{−3}*D*/*U*_{c} to 10 × 10^{−3}*D*/*U*_{c}.

## III. RESULTS

### A. Validation

To the best of our knowledge, the only 3-D results reported in the literature for the case of flow over a confined circular cylinder in a channel are from Rehimi *et al.*,^{18} who, however, considered a different blockage ratio of β = 1/3. For the blockage ratio examined in the present paper (β = 1/5), the only previously available results had been obtained from 2-D simulations. Therefore, 2-D numerical simulations were performed first, allowing the validation of the numerical code and the chosen parameters through comparisons with existing literature.

The critical Reynolds number *Re*_{cr} for the transition from steady to unsteady flow, the mean drag coefficient $C\u2212$ _{D}, and the Strouhal number *St* were computed and compared against previous studies. Our results for *Re*_{cr}, and the corresponding Strouhal number *St*_{cr}, are listed in Table III. When compared with the values reported in Refs. 21, 24, and 25, these show an excellent agreement. Figure 3 displays a comparison of both the drag coefficient and the Strouhal number versus *Re*. As shown, the agreement is again very satisfactory, confirming the accuracy of the numerical code and the choice of appropriate numerical parameters.

. | Present . | Sahin and Owens (Ref. 24) . | Zovatto and Pedrizzetti (Ref. 21) . | Chen et al. (Ref. 25)
. |
---|---|---|---|---|

Re_{cr} | 69.5 | 69.9 | 68.9 | 69.3 |

St_{cr} | 0.1567 | 0.1567 | — | 0.1559 |

### B. Transition to three-dimensionality

In order to carefully investigate transitional effects, a series of nine three-dimensional simulations have been carried out for the range 150 ≤ *Re* ≤ 390, in steps of 30. A summary of the results obtained is listed in Table IV.

Cases . | Re
. | W
. | $C\u2212D$ . | $C'L$ . | St
. | C_{pb}
. |
---|---|---|---|---|---|---|

1 | 150 | 8 | 1.2389 | 0.271 | 0.1850 | 0.8487 |

2 | 180 | 8 | 1.2253 | 0.342 | 0.1893 | 0.8845 |

3 | 210 | 8 | 1.166 | 0.25 | 0.1852 | 0.832 |

4 | 240 | 12 | 1.169 | 0.30 | 0.1895 | 0.880 |

5 | 270 | 8 | 1.170 | 0.37 | 0.1949 | 0.913 |

6 | 300 | 8 | 1.172 | 0.37 | 0.1989 | 0.942 |

7 | 330 | 8 | 1.139 | 0.33 | 0.2015 | 0.918 |

8 | 360 | 8 | 1.120 | 0.30 | 0.2035 | 0.895 |

9 | 390 | 8 | 1.099 | 0.26 | 0.2053 | 0.883 |

Cases . | Re
. | W
. | $C\u2212D$ . | $C'L$ . | St
. | C_{pb}
. |
---|---|---|---|---|---|---|

1 | 150 | 8 | 1.2389 | 0.271 | 0.1850 | 0.8487 |

2 | 180 | 8 | 1.2253 | 0.342 | 0.1893 | 0.8845 |

3 | 210 | 8 | 1.166 | 0.25 | 0.1852 | 0.832 |

4 | 240 | 12 | 1.169 | 0.30 | 0.1895 | 0.880 |

5 | 270 | 8 | 1.170 | 0.37 | 0.1949 | 0.913 |

6 | 300 | 8 | 1.172 | 0.37 | 0.1989 | 0.942 |

7 | 330 | 8 | 1.139 | 0.33 | 0.2015 | 0.918 |

8 | 360 | 8 | 1.120 | 0.30 | 0.2035 | 0.895 |

9 | 390 | 8 | 1.099 | 0.26 | 0.2053 | 0.883 |

It is important to note that the vortex shedding in the present 3-D simulations is initiated without imposing or forcing any artificial or external flow-disturbances. Instead, the round-off and truncation errors, which are uniformly distributed over the whole computational domain, were allowed to generate the self-excitation needed for the flow to naturally develop three-dimensionalities. This choice was associated with long computational times that became even longer for the cases that were closer to the transitional regime, as we discuss further below. However, it was motivated from the lack of previously reported data on this flow configuration.

The first steps of transition to three-dimensionality are reflected in the amplification of the spanwise velocity-component, *u*_{z}, in the near wake. Figure 4 shows the time evolution of *u*_{z} in the near-wake, at Reynolds numbers ranging from 180 to 300. For the sake of clarity, each signal is truncated at a different evolution time to avoid overlaps. For the lower Reynolds numbers considered here, namely *Re* = 180 (and *Re* = 150, which is not shown in Figure 4), the flow did not exhibit any sign of spanwise fluctuations, an indication that it remained completely two-dimensional. In contrast, for higher values of the Reynolds number, *Re* ≥ 210, the flow showed a non-linear growth of *u*_{z}, indicating the inception of three-dimensionality. As a general trend, with increasing Reynolds number, the *u*_{z} amplification was found to initiate earlier. Thus, the present simulations indicate that the transition in the cylinder’s wake occurs within the interval between *Re* = 180 and *Re* = 210. This is consistent with the results of Camarri and Giannetti,^{19} who carried out a Floquet stability analysis for the same configuration and found a critical value of *Re*_{crA} ≈ 201 for the transition to a three-dimensional state.

Transition to three-dimensionality is also depicted on the time history of the lift, *C*_{L}, and drag, *C*_{D}, coefficients as shown in Figure 5. After a transient period, during which the flow remains in a two-dimensional state and periodic vortex shedding is observed, the signal eventually losses its coherence and a drastic reduction in the lift and drag forces occurs, corresponding to a three-dimensional state of the flow. The flow develops these three-dimensional effects soon after the amplification of the spanwise velocity *u*_{z}. The length of the transient period was particularly long for *Re* = 210, approaching 500 time units, but it got shorter as the Reynolds number was further increased. For all other cases up to *Re* = 390, the transient period was ranging between 50 and 200 time units. This behavior bears strong resemblance to the simulations of Mittal and Balachandar^{28} for the case of an unconfined circular cylinder.

### C. Effect of Re on Strouhal number St and base pressure coefficient *C*_{pb}

As the *Re* number is increased, modifications in the dynamics of the wake structure produce distinct changes in the shedding frequency and the pressure field. The transient nature of the flow necessitates long simulation times before reliable statistics can be collected, especially for *Re* values close to the wake transition critical points. In our case, simulation times were made even longer because of relying on self excitation to trigger transition. Once a statistically stationary state was reached in the computation, the calculation was continued for another 100 to 150 vortex shedding cycles to obtain the time-averaged values. Strouhal numbers have been obtained from the lift coefficient signals using Welch’s averaged periodogram method.^{29} A Hamming window was applied to each overlapping segment of data.

The variation of the Strouhal number, *St*, as a function of the Reynolds number is shown in Figure 6, where results from the present 2-D and 3-D simulations are compared to the experimental results of Williamson^{2} for the *unconfined* case. Dashed vertical lines indicate the critical Reynolds numbers, which, according to Camarri and Giannetti,^{19} mark the onset of different wake instabilities for the case of a confined circular cylinder with the same blockage ratio as used in the present study. As shown in Sec. III B, the flow remains completely two-dimensional up to *Re* = 180. Not surprisingly then, results from the 2-D flow simulations are identical to those from the 3-D simulations in this range of *Re* values. However, at *Re* = 210, where three-dimensional effects start to develop, a significant difference between the 2-D and 3-D cases is observed. In the 3-D case, *St* undergoes a sudden drop of approximately 4% that persists up to *Re* = 240. As *Re* is further increased, differences are observed to be less significant, and another discontinuity in the *St–Re* relationship shows up, which this time is not hysteretic. At *Re* = 270, the difference between 2-D and 3D cases is around 2.5%, while for *Re* = 300–390, the difference decreases to 1.5%. Interestingly, results from the three-dimensional simulations compare well to the experimental results of Williamson for the case of an unconfined cylinder, but with a small delay in the critical *Re* values. This offset, towards higher Re in our case, seems to be attributed to the additional confinement from the channel walls, which presumably stabilizes the flow and produces higher transitional Reynolds numbers.^{30}

Williamson^{7} associated these two discontinuous changes in the *St–Re* curve, with the development of different instabilities in the wake; mode A instability combined with intermittent vortex dislocations and mode B instability. In accordance with these observations, current results indicate the existence of two instability regions: the first occurs around *Re* = 210 and corresponds to the onset of three-dimensionalities in the flow, while the second occurs around *Re* = 270. We will discuss about these instabilities in greater detail in Sec. III D. These findings are consistent with the results of Camarri and Giannetti,^{19} who carried out a linear stability analysis for the case of a confined circular cylinder having the same blockage ratio as in the present study and also found the existence of two instability regions. They found the critical value for the onset of the second instability to be *Re*_{crB} ≈ 256. According to their work, these instabilities have the same symmetries as the mode A and mode B instabilities found in the unconfined case. However, no drop in the Strouhal number is reported in their study. This is not surprising taking into consideration the limitations of the linear stability approach, and the fact that this transition feature, observed at *Re* ≈ *Re*_{crA} and above, is not a linear effect of the instability but is a result of strongly nonlinear phenomena.^{9,10}

Despite its value in assessing flow instabilities, the evolution of the base pressure coefficient, *C*_{pb}, as a function of the Reynolds number has not been previously discussed in the literature for the confined case. Figure 7 shows this evolution for the present 2-D and 3-D simulations together with the experimental results of Williamson and Roshko^{3} and the data of Dennis and Chang^{31} for comparison with the unconfined case. Again, as expected, results from the 2-D and 3-D simulated flows are identical up to *Re* = 180, where the flow remains two-dimensional. For *Re* > 210, where three-dimensional effects start to take place, the results from the 3-D simulation show a marked drop in the level of *C*_{pb}, which remains undetected in the 2-D simulations. Once again, the *C*_{pb} predictions from the current 3-D simulations exhibit a very similar evolution with *Re* as found experimentally by Williamson and Roshko^{3} for the unconfined case.

According to Williamson^{7} and Roshko,^{32} the variations of the base pressure coefficient correspond to the presence of instabilities in the flow. In the unconfined case, the drop in the base pressure coefficient at *Re* = 180 is thought to be to related to the presence of mode-A instability and vortex dislocations.^{7} At *Re* = 260, on the other hand, there is a local maximum that corresponds to a saturation of the primary instability growth.^{32} Around that point, the secondary spanwise structure changes to one with smaller scale, mode B instability.^{7} As Re is further increased, three-dimensional structures become more disordered and the base suction begins to decrease. Our computations show the same characteristics in the variation of *C*_{pb} with *Re*, the only difference being again a small delay in the initiation of 3-D effects in the confined case relative to the unconfined results (*Re* = 210 instead of *Re* = 180) and similar delay in the saturation effects (*Re* = 270–300 instead of *Re* = 260), which are attributed to the additional confinement of the channel walls. Based on these results, one would expect the overall structure of the confined wake to be similar to that of the unconfined case. This expectation is confirmed by the findings presented in the following sections.

### D. Instabilities in the wake

In order to identify three-dimensional vortex structures in free shear flows, the vorticity magnitude is usually used. However, in our case due to the existence of vorticity at the channel walls, using the vorticity magnitude to identify vortices would result in deformed vorticity structures and difficulties in visualizing them. The λ_{2} criterion is by far more appropriate for boundary layer type of flows and it is defined as the second eigenvalue of *S*^{2} + *X*^{2}, where *S* and *X* denote the symmetric and antisymmetric parts of the velocity gradient tensor, respectively.^{33} Thus, iso-surfaces of λ_{2} were used in order to exclude the wall shear region and focus on the swirling motion of the primary and induced vortices of the cylinder wake. Indeed, through the application of the λ_{2} criterion, we were able to clearly detect the vortical motions of interest for such a highly three-dimensional wall-bounded shear flow.

A visual impression of the three-dimensional structures related to the different instabilities found in the flow can be obtained from Figure 8. There, we show snapshots of iso-surfaces of λ_{2} normalized by its absolute minimum, λ_{2,min}, for different Reynolds numbers in increasing order. Iso-surfaces are colored by the streamwise vorticity component, ω_{x}, to reveal the streamwise rotation direction of each vortical structure. The spanwise rollers essentially identify the primary vortex cores.

At *Re* = 240 [Fig. 8(a)], a Reynolds number corresponding to a regime well after the onset of the first instability, a spanwise waviness of the primary vortex cores is observed in the cylinder’s wake, along with the formation of counter-rotating streamwise vortex pairs. Over successive half cycles of vortex shedding, vortex pairs of opposite-sign vorticity are formed. The flow displays a dominant spanwise wavelength of around 4*D*, in agreement with the value of 4.65*D* reported by Camarri and Gianneti^{19} for a similar configuration. This periodic, out-of-phase, three-dimensional flow pattern is topologically similar to that of mode A,^{7} at least in the near wake region (as it will be discussed in the Sec. III E). Mode A instability appears in the wake of an unconfined circular cylinder with a spanwise wavelength of approximately 4*D*.^{6,7,9}

For *Re* = 300 and above, finer scale streamwise vortex pairs are observed [Figs. 8(d) and 8(e)], with a distinctly smaller spanwise wavelength of approximately 1*D*. This time, successive vortex pairs from one braid layer to the next, have the same orientation, forming an in-phase streamwise vortex pattern. This instability, in contrast to mode A, is restricted to the near wake only and it is not found in the far downstream locations. With increasing Re, the number of streamwise vortices is increased, and the flow becomes more distorted. This flow pattern is analogous to that of mode B vortex shedding described by Williamson^{7} for the case of an unconfined circular cylinder. It also compares well with the predicted critical wavelength of 0.86*D* found by Camarri and Gianneti^{19} for the confined case.

In the case of *Re* = 270, which is close to the instability threshold, $ReB\u2248256$, found by Camarri and Gianneti,^{19} modes of vortex shedding similar to both mode A [Fig. 8(b)] and mode B [Fig. 8(c)] were observed at different instances of the flow. This intermittent nature of the flow is in line with the experimental observations of Williamson,^{7} the direct numerical simulations of Henderson,^{10} and the stability analysis of Barkley^{34} for the unconfined case. They have shown that the transition from mode A to mode B is associated with a gradual transfer of energy from one mode to the other, resulting in a mixed-mode state approximately in the Re range 230–265. Williamson attributes this to the intermittent swapping between the two modes, rather than the coexistence of both modes. Also, Behara and Mittal^{13} in their numerical investigations demonstrated the swapping between modes up to *Re* = 275 for the unconfined case.

The present three-dimensional simulations are the first to reveal the presence of natural vortex dislocations in the confined wake of a circular cylinder. Williamson^{4,7} has shown that large-scale spot-like vortex dislocations are an intrinsic phenomenon of wake transition. These irregularities occur spontaneously along the span as a natural feature of the wake flow and are associated with the presence of mode A instability. However, few three-dimensional simulations have captured them. Zhang *et al.*^{6} reproduced numerically Williamson’s vortex-dislocations, after applying strong localized spanwise inhomogeneity in the initial conditions. Braza *et al.*^{12} were the first to obtain the vortex dislocations naturally by means of a complete Navier-Stokes simulation. In the present work, clear observations of the existence of natural vortex dislocations are presented for the case of a confined circular cylinder.

The occurrence of such dislocations can clearly be seen in Figure 9, where snapshots of iso-surfaces of λ_{2} = −0.4 and pressure *p* = −0.4 are presented for *Re* = 240 at *t* = 870*D*/*U*_{c}. Iso-surfaces of λ_{2} are colored by streamwise vorticity component, ω_{x}. At the occurrence of dislocation, in the vicinity of *z*/*D* = −2, the span-wise coherence of the primary vortex core is lost and a break in the continuity of the vortex tube can be seen. Moreover, during this phase of the flow, dislocations affect the shedding phenomenon, which as a result becomes irregular.

One can identify the occurrence of dislocations through their signature in instantaneous velocity signals and the corresponding spectra at different spanwise positions. Figure 10 displays the crossflow velocity component signal from a series of probes spaced in the spanwise direction and placed at the intersection of the horizontal plane y/D = 0.5 and the vertical plane x/D = 1.5. This is a location slightly downstream from the cylinder. At certain spanwise positions, the normal Karman vortex shedding gives way to a pronounced modulation in the time history of fluctuating velocities. These can be seen as the marked areas in Figure 10(a), or more clearly in Figure 10(b), where two samples of velocity signals are displayed in an interval corresponding to the time instant of the snapshot shown in Figure 9. These irregularities or glitches in the velocity signal correspond directly to the passage of a vortex dislocation structure past the measuring probe.^{4} The spontaneous occurrence of such dislocations at different spanwise locations was found to be regular in time.

Another way to visualize dislocations is by looking at the corresponding power spectra density of the velocity signal (Fig. 11). Braza *et al.*^{12} found that in the regions where a vortex dislocation occurs, the spectral energy of the fundamental frequency is reduced considerably. This behavior is clearly seen in Figure 11. In the case of *Re* = 240, where the flow exhibits mode A vortex shedding, there is a significant decrease in the spectral energy in the vicinity of *z*/*D* = 4 and *z*/*D* = −3. The drop, at certain spanwise positions, reaches up to 50% from the maximum value. On the other hand, for the cases of *Re* = 270 and *Re* = 300, which correspond to a mixed or pure mode B vortex shedding, changes in the amplitude of the spectral density are much smaller, approximately 15% and 20%, respectively. In these cases, the flow is devoid of vortex dislocations, in agreement with the experiments of Williamson^{4} and the numerical simulations of Behara and Mittal.^{13}

### E. Effects of confinement

In order to identify differences between the wake characteristics of confined and unconfined cylinders, two more simulations were performed, at *Re* = 240 and *Re* = 300, this time for an unconfined cylinder. To accommodate the unconfined simulations, grids G2W and G2F were laterally extended from 5*D* to 50*D*, in order to minimize blockage effects,^{9,11} while providing a comparable resolution to the corresponding confined cases.

The effects of confinement on mode A flow structures are examined first. A comparison between the confined and the unconfined case is shown in Figure 12, where snapshots of iso-surfaces of the λ_{2} criterion are plotted at *Re* = 240. Although mode A is identified in the confined case, the downstream evolution of the cylinder wake is strongly affected by confinement. In the near wake region, i.e., for *x*/*D* < 5, mode A is clearly identifiable in the confined case, and one cannot discern clear differences between the vortical patterns generated by the confined and unconfined cases. Moving further downstream though, one finds that the structure of the confined wake departs from that of a standard unconfined cylinder wake. For example, at 5 < *x*/*D* < 17, the hairpin structures of the braid shear layer that are a standard feature of the unconfined wake, can still be identified in the confined case, but their motion and shape is significantly modified. These differences may be attributed to the inversion of the Von Karman vortices caused by wall interactions in the confined case, as described by Camarri and Gianneti.^{19} For *x*/*D* > 17, the confined wake is much more fragmented, without a clear presence of primary vortex cores, which in the unconfined case are seen to persist this far downstream. In the same region, one finds a stronger streamwise alignment of vortices in the confined case compared to the unconfined case. In the confined case, hairpin vortices sustain their coherency and persist for larger distances, up to *x*/*D* ≈ 30.

The breakdown of the primary vortex cores and the fragmented nature of the flow depicted by iso-surfaces of λ_{2} = −0.3 (0.5% of its absolute minimum) is represented over successive time instants at *Re* = 240 in Figure 13. In order to have a more clear view of the spanwise vortex cores, rather than the streamwise vortices, iso-surfaces are colored by the spanwise vorticity component ω_{z}. At *t* = 1111*D*/*U*_{c} one can observe that from approximately *x*/*D* ≈ 17 and above, the primary vortex cores lose there coherence and break into smaller structures that become aligned in the streamwise direction as they are advected further downstream. If we focus our attention on one such a pair of counter-rotating spanwise vortex cores [Fig. 13(a), box A], and follow them in time as they move further downstream, the initial small waviness grows, and they seem to be gradually pulled backwards and towards the channel walls. When they reach the vicinity of *x*/*D* = 17 [Fig. 13(b), box B], the front counterclockwise vortex roller (top roller, red online) starts losing its continuity and breaks. At *x*/*D* ≈ 22, eventually both vortices break down and they are transformed to U-shaped vortices [Fig. 13(c), box C]. The legs of these structures are later separated and transformed to streamwise rollers as they are advected downstream [Fig. 13(d), box D]. Interestingly, this pattern of behavior is systematic and all primary vortex cores are stretched and eventually break down in a similar manner and around the same spanwise positions.

Figures 14(a) and 14(b) show iso-surfaces of the streamwise component of the velocity field for *u*_{x} = 0.5 and 1.5, respectively, for *Re* = 240 at *t* = 1130*D*/*U*_{c}. Looking at Figure 14(a), one can observe close to the channel walls “streaks” of low velocity fluid that protrude towards the center of the channel. These “streaks” are aligned in the streamwise direction and occur on both sides of the channel walls at similar spanwise positions. In between these low velocity fields, regions of high velocity fluid are observed that are also aligned in the streamwise direction, as seen in Figures 14(b) and 14(d).

This organized, “symmetric” mixing of low and high velocity fields results from the interaction of the streamwise vortex pairs from mode A instability and the channel walls. Figure 14(c) shows a contour plot of *u*_{x} in the plane *x*/*D* = 7.2 (plane A). Contour lines of streamwise vorticity, ω_{x} = ±0.6, ±0.8, are superimposed to give information on the location and rotation direction of the streamwise vortex pairs. The induced velocity due to the streamwise vortex pairs drives low velocity fluid towards the centreplane of the wake and high velocity fluid towards the channel walls. Because mode A instability gives an out-of-phase streamwise vortex pattern at particular spanwise locations, streamwise vortex pairs of opposite sign reside on opposing sides of channel walls. As a result, the deformation and stretching of low and high velocity fluid towards and away from the centreplane of the wake occur at the same spanwise positions on both sides. The self-sustaining nature of mode A, that gives an array of streamwise vortex pairs that travel further downstream at the same spanwise positions, explains why these low and high speed regions are aligned in the streamwise direction.

Figures 14(e) and 14(f) show contour plots of *u*_{x} in the planes *z*/*D* = 2 (plane C) and *z*/*D* = 4 (plane D) respectively, which correspond to areas in the wake where high and low velocity fields are more prominent accordingly. In order to investigate the effect of the velocity field on the primary vortex cores, lines representing values of spanwise vorticity, ω_{z} = ±0.8, are displayed on top. As seen in Figure 14(e), the primary vortex cores strongly interact with the high speed velocity field and are forced to follow a trajectory close to the centerline as they shed downstream. On the other hand, in the plane *z*/*D* = 4 [Fig. 14(f)], as segments of the vortex cores are caught in the low speed region, the primary vortex cores follow an oblique trajectory closer to the channel walls and progressively slow down (see vortex cores *P*1 and *P*2). This accelerating and slowing down of segments of the primary vortex cores at different spanwise locations causes the enhanced waviness observed. Around *x*/*D* = 17, the stretching of the vortex core becomes so strong that forces it to break down.

So far, the discussion of confinement effects has been limited to mode A instability. In the case of mode B, when compared with the unconfined case, the downstream evolution of the cylinder wake also seems to be strongly affected by confinement in a similar manner with mode A as shown in Fig. 15 for *Re* = 300. For higher values of *Re*, the break down of the primary vortex cores ensues further upstream [up to *x*/*D* = 12 for *Re* = 390, as shown in Fig. 8(e)] and the flow becomes more fragmented. However, in the presence of mode B instability, it was not possible to observe a regular pattern of vortex breakdown.

There are mainly two features of the flow responsible for this behavior. Unlike mode A, mode B instability leads to an in-phase streamwise vortex array with a smaller spanwise wavelength [Figs. 8 and 16(c)]. Also, a spanwise wandering of the streamwise vortices from one half-cycle to the next is observed, in line with the experimental visualizations of Williamson^{7} for the unconfined case. As a result, induced velocity from the streamwise vortices was found to generate irregular streams of low velocity, as shown in Figure 16(a). In contrast with mode A, these slow moving regions, were neither aligned along the streamwise direction nor observed at symmetric locations on the channel walls. In addition, no clustered regions of high velocity extending in the streamwise direction were identified [Fig. 16(b)]. These features promote the amplified distortion observed for mode B instability and result in the irregular break down of the primary vortex cores. With increasing *Re*, these effects become even more pronounced leading to the breakdown of the primary vortex cores earlier downstream.

## IV. CONCLUSIONS

Direct numerical simulations of the two- and three-dimensional flow around a circular cylinder placed in a plane channel have been performed. The blockage ratio was kept constant at 1/5, and the Reynolds number was varied between 10 and 390.

Present results indicate that up to *Re* = 180, the flow remains two-dimensional. For higher values of the Reynolds number, *Re* ≥ 210, the flow develops three-dimensional effects which are depicted on the time history of the lift, *C*_{L}, and drag, *C*_{D}, coefficients. Two discontinuous changes were detected in the *St*–*Re* curve corresponding to different spanwise instabilities in the wake, an effect that is linked to inherently nonlinear mechanisms and which previous studies failed to capture. We have also confirmed for the first time, that the critical points are also reflected in the *C*_{pb}–*Re* relationship, as was also found by Williamson for the unconfined case. Similar to the case of an unconfined circular cylinder, mode A 3-D shedding was observed for *Re* = 210 and 240. For *Re* ≥ 300, mode B vortex structures were detected. At *Re* = 270, the flow exhibited an intermittent swapping between the two modes. The intermittent presence of naturally occurring vortex dislocations, as a fundamental feature of wake transition, was also demonstrated. This is the first time that the existence of these instabilities has been confirmed via full 3-D simulations for the confined circular cylinder in a channel.

The present work leads to a clarification of how the shape and evolution of mode A and mode B instabilities are affected downstream by the confinement of the channel walls. In case of mode A, organized, “symmetric” mixing of low and high velocity fields is observed downstream, which eventually forces the primary vortex cores to break in a consistent and systematic way. In the case of mode B, irregular streams of low and high velocity are observed, which result to the irregular breakdown of the vortex cores. Understanding the hydrodynamic field and knowledge of the velocity profiles is important to properly predict local phenomena such as heat transfer and corrosion effects.

## ACKNOWLEDGMENTS

This work has been performed under the UCY-CompSci project, a Marie Curie Transfer of Knowledge (TOK-DEV) grant (Contract No. MTKD-CT-2004-014199), funded by the CEC under the 6th Framework Program and under the contract of association ERB 5005 CT 99 0100 between the European Atomic Energy Community and the Hellenic Republic. Partial support has also been received through a Center of Excellence grant from the Norwegian Research Council to the Center for Biomedical Computing.