`ISRN Mathematical PhysicsVolume 2012 (2012), Article ID 630801, 16 pageshttp://dx.doi.org/10.5402/2012/630801`
Research Article

## Numerical Simulation of Viscous Flow over a Square Cylinder Using Lattice Boltzmann Method

1Department of Aeronautical Engineering, Noorul Islam Centre for Higher Education, Noorul Islam University, Kanyakumari 629180, India
2Department of Mechanical Engineering, Indian Institute of Technology Guwahati, Guwahati 781039, India

Received 6 April 2012; Accepted 9 July 2012

Academic Editors: A. Sanyal and F. Sugino

Copyright © 2012 D. Arumuga Perumal et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

This work is concerned with the lattice Boltzmann computation of two-dimensional incompressible viscous flow past a square cylinder confined in a channel. It is known that the nature of the flow past cylindrical obstacles is very complex. In the present work, computations are carried out both for steady and unsteady flows using lattice Boltzmann method. Effects of Reynolds number, blockage ratio, and channel length are studied in detail. As good care has been taken to include appropriate measures in the computational method, these results enjoy good credibility. To sum up, the present study reveals many interesting features of square cylinder problem and demonstrates the capability of the lattice Boltzmann method to capture these features.

#### 1. Introduction

Incompressible viscous flow around a cylinder with square cross-section confined in a channel is also one of the important problems in computational fluid dynamics (CFD) [14]. It is known that the presence of wall confinement is one of the common examples in many practical flow situations (e.g., road vehicles, heat exchangers, flow around tall buildings, suspension bridges, etc.), and the study of effect of wall confinement on flow field characteristics is of crucial importance [57]. The main advantage of this square cylinder problem is the relative simplicity of the geometry compared with circular cylinder. Another feature of the present square cylinder problem is that it has fixed separation points.

The present square cylinder confined in a channel is mainly characterized by flow separation, periodicity, and well-known vortex street formation. At very low Re, the flow is laminar and steady and remains attached to the square cylinder. As the Reynolds number increases, the flow separates from the trailing edge but remains laminar and steady up to Re of about 50. Beyond this Reynolds number (Re) unsteadiness develops, and the flow becomes periodic. At these Reynolds numbers localized regions of high vorticity are shed alternatively from either side of the cylinder and are convected downstream. At lower values of postcritical Reynolds numbers the flow is still laminar. Because of the popularity of the present problem in terms of applications and importance a number of experimental and numerical results are available in the literature for flow past a square cylinder mounted in a confined channel [813].

In most of the available literature, the square cylinder flow field calculation is based on solving the Navier-Stokes equation using a finite difference method (FDM), finite volume method (FVM), (or) finite element method (FEM). In 1984, Davis et al. [1] studied the confined flow past a square cylinders for a wide range of Re with different blockage ratios such as B = 4 and 6 experimentally and numerically. In 1990, Franke et al. [2] presented numerical calculations of laminar vortex-shedding flows for flow past circular and square cylinders. In 1992, Mukhopadhyay et al. [3] carried out two-dimensional numerical simulation of flow over a square cylinder with two different blockage ratios (B = 4.0 and 8.0). In 1993, Suzuki et al. [4] simulated flow past a square cylinder on a nonequidistant grid for different Reynolds number range and blockage ratios of 5.0 and 20.0. Sohankar et al. [5, 6] presented calculations of unsteady two-dimensional flow around square and rectangular cylinders at incidence for low Reynolds numbers. In 2000, Pavlov et al. [7] examined the results of numerical calculations of transient, 2-D incompressible flow past a square cylinder at different Reynolds numbers. In 2000, Breuer et al. [8] also investigated the confined flow past a square cylinder placed inside a plane channel (blockage ratio B = 8.0) by two different numerical techniques such as the lattice-Boltzmann automata (LBA) and the finite volume method (FVM). In 2002, Wan et al. [9] proposed a discrete singular convolution-finite subdomain method (DSC–FSM) for the analysis of incompressible viscous flows to multiply connected complex geometries, namely, the flow over a backward-facing step and the laminar flow past a square cylinder. In 2004, Roy and Bandyopadhyay [10] investigated flow past a square cylinder placed inside a channel with two different blockage ratios, namely, 4 and 8 at different Reynolds numbers. In 2005, Abide and Viazzo [11] presented the utility of the compact scheme projection decomposition method approach for two benchmark problems, namely, the flow over a backward-facing step and the flow past a square cylinder. In 2005, Hasan et al. [12] also presented an approach for handling outflow boundary condition based on the radial behaviour of the velocity field at large distances from the body in external incompressible viscous flows.

The development of lattice Boltzmann method (LBM) is the promising method that uses different kind of nonconventional techniques for applications in CFD [1418]. It is now observed that an increasing number of researchers simply use the LBM as an alternative to conventional numerical methods for the Navier-Stokes equation. To demonstrate the ability of present lattice Boltzmann Method to handle flows over square geometries confined in a channel without much difficulty, in the present work, the influence of the locations of channel wall, that is, blockage effect, outlet boundary location, and different lengthwise cylinder locations are considered. Another objective of this present work is to capture the unsteady flow states characterized by the square cylinder at relatively higher Reynolds numbers.

#### 2. Lattice Boltzmann Method

The lattice Boltzmann method with single-relaxation-time (LBM-SRT) model, which is a commonly used LBM, is given by [14] the following: where and are the particle and equilibrium distribution functions at (), is the particle velocity along the th direction, and is the single-relaxation-time parameter that controls the rate of approach to equilibrium. For simulating two-dimensional flows, the two-dimensional nine-velocity square lattice model (D2Q9) with nine discrete velocities () is used (Figure 1).

Figure 1: Two-dimensional nine-velocity square lattice model.

The discrete velocity set is written as [14] In the D2Q9 square lattice, a suitable equilibrium distribution function that has been proposed is [15] where the lattice weights are given by and . The relaxation time that fixes the rate of approach to equilibrium is related to the viscosity by [15] where is the kinematic viscosity measured in lattice units. The macroscopic quantities such as density and momentum density are obtained as velocity moments of the distribution function as follows: where . The density is determined from the particle distribution function. The density and the velocities satisfy the Navier-Stokes equations in the low-Mach number limit. This can be demonstrated by using the Chapman-Enskog expansion [16].

##### 2.1. Code Validation

First, the developed LBM code is used to compute the flow past a single circular cylinder confined in a channel. For a blockage ratio B = H/D = 8 computations are carried out at Re = 60 with the help of lattice Boltzmann method. The lattice size of 500 × 80 is used for the present configuration. In Table 1 we present the coefficient of drag for different steady-flow Reynolds numbers. Expectedly as the Reynolds number increases, coefficient of drag () decreases. It is worth mentioning that the numerical simulations of our LBM are much closer to existing available results.

Table 1: Coefficient of drag () for different circular cylinder steady-flow Reynolds numbers by LBM.

##### 2.2. Problem and Associated Features of LBM

The geometry of the computational domain for flow past a square cylinder placed in a confined channel is shown in Figure 2. Here the cylinder is placed in the middle of the channel of height with blockage ratio B = H/D = 8, where is the side of the square cylinder. The channel length is fixed at L/D = 50 to reduce the influence of outflow boundary condition on accuracy. An upstream length of l = L/4 (or 12.5D) has been chosen. The Reynolds number is the main parameter which changes the flow behaviour; it can be defined by UD/, where is the kinematic viscosity. No-slip boundary conditions are used on the lower and upper boundaries and uniform flow velocity at the inlet.

Figure 2: Schematic diagram of the flow past a square cylinder confined in a channel.
##### 2.3. Boundary Conditions

In LBM, implementing the boundary condition is a difficult task because of the necessity of imposing conditions in terms of particle distribution functions. Boundary conditions on different walls are as follows. (i) Top and bottom walls: on the top and bottom walls, we use the well-known bounce-back boundary condition [14] which indicates no-slip. At the inlet boundary we impose bounce-back condition with momentum proposed by Yu et al. [16] that can be written as (ii) Outlet boundary: at the outlet boundary, we employ a second-order accurate extrapolation boundary condition, which is written as [16] where is the number of lattices in the x-direction. (iii) On the solid curved boundary, that is, on the surface of the cylinder, a robust second-order accurate boundary treatment proposed by Bouzidi et al. [17] is used for the particle distribution function.

##### 2.4. Force Evaluation

The momentum exchange method [16] is used for the force evaluation in the present section. The total resultant fluid, force F, acting on a solid body by fluid can be written as [16] where is the number of nonzero lattice velocity vectors, and is an indicator, which is 0 at and 1 at . The two most important characteristic quantities of flow around a cylinder are the coefficient of drag and coefficient of lift. The coefficients are defined as [18]

where and are the and components of the total fluid force acting on the cylinder. In the present work, LBM computations for steady state are carried out till the following convergence criterion is satisfied: where is the fluid velocity, and is the iteration level.

##### 2.5. Themes of Studies Are as Follows:

(i)for a fixed blockage ratio to study the square cylinder flow field characteristics for various Reynolds numbers, (ii)for different blockage ratios to study the square cylinder flow field characteristics, (iii)for different lengthwise cylinder locations to study the square cylinder flow field characteristics, and(iv)for different outlet boundary locations to study the square cylinder flow field characteristics.

#### 3. Results and Discussion

##### 3.1. Fixed Bclokage Ratio

Computations are carried out for a blockage ratio B = H/D = 8 at different Reynolds numbers with the help of single-relaxation-time LBM method based on the D2Q9 model. The moderate lattice size of 800 × 128 is used for the present flow configuration.

Case i (steady flow). Figure 3 shows the streamline patterns around the square cylinder for different Reynolds numbers (Re = 4, 10, 20, and 30) with B = 8. At very low Reynolds number (Re = 4), the flow is laminar, steady, and slightly separated from the cylinder depicted in Figure 3(a). With the increase in Reynolds number (Re = 10) flow separates fully and two vortices symmetrically placed about the channel centreline as shown in Figure 3(b). It is seen (Figures 3(c) and 3(d)) that two vortices downstream of the cylinder, symmetrically placed about the channel centreline, develop and remain attached to the cylinder. From Figures 3(a)3(d), it is also seen that for steady flows the sizes of the vortices increase with Reynolds number.

Figure 3: Streamline patterns for steady flows past a square cylinder at different Reynolds numbers for a blockage ratio B = 8: (a) Re = 4 (b) Re=10 (c) Re = 20 and (d) Re = 30; lattice size: 800 × 128.

Case ii (unsteady flow). Beyond the critical Reynolds number (Re) unsteadiness develops, and the flow becomes periodic. In the present section, periodic flow is computed at Re = 60, 100, and 150. Figure 4 shows the streamline patterns at a certain instant for Re = 60, 100, and 150. From the streamline patterns, it is seen that the flow pattern changes at different Reynolds numbers.

Figure 4: Instantaneous streamline patterns for flows past a square cylinder at different Reynolds numbers for a blockage ratio B = 8: (a) (b) and (c) ; lattice size: 800 × 128.

Figure 5 depicts the vorticity contours for various Reynolds numbers with B = 8. In Figure 5(a), the solid lines represent the positive vorticity (anticlockwise sense), and the dashed lines represent the negative vorticity (clockwise sense). The positive vortex appears on the lower side of the square cylinder while the negative vortex is on the upper side of the square cylinder for all the Reynolds numbers. Steady flow vortices seem to be more and more elongated as the Reynolds number increases (Figures 5(a), 5(b), and 5(c)). The flow behind the square cylinder is symmetric for precritical Reynolds numbers. Figures 5(d) and 5(e) depicts the flow field at higher Reynolds number of 60 and 100. In these Reynolds numbers the flow becomes periodic, and the frequency of vortex shedding increases with Reynolds number.

Figure 5: Vorticity contours for the flow past a square cylinder for different Reynolds numbers with B = 8: (a) , (b) , (c) , (d) , and (e) Re = 100. lattice size: 800 × 128.

Based on the numerical investigation a lower value of the Reynolds number for periodic vortex motion () was determined by Kelkar and Patankar [13]. In the present work, Figure 6 depicts the streamline pattern and vorticity contour for Re = 52. In our present LBM computations, periodicity develops for Reynolds number 52 at a particular blockage ratio B = 8, and periodic shedding starts without any triggering. Figure 7 shows the instantaneous coefficient of lift values with respect to time demonstrating the periodicity of flow at Re = 52.

Figure 6: Instantaneous streamline pattern and vorticity contours for flows past a square cylinder at Re = 52.
Figure 7: Instantaneous coefficient of lift for flow past a square cylinder at Re = 52.

Figure 8 shows the periodic time variation of the lift coefficient () and drag coefficient () for the flow past a square cylinder at Re = 60 and 100. It is seen that, for each time period, has two crests and two troughs of unequal amplitudes, which are a consequence of the periodic vortex shedding from the top and bottom surfaces. The same figure shows that for each time period has just one trough and one crest; this is because is not influenced by the pressure distribution on the rear side of the cylinder, which is strongly influenced by the vortex shedding at the top and bottom sides of the cylinder. The value of the lift force fluctuation is directly connected to the formation and shedding of the eddy and: therefore, its value varies between a positive maximum and a negative maximum of equal magnitude. Table 2 gives the mean drag coefficient values for different Reynolds numbers. The presented drag coeffcients are close to the results of [8] which were obtained by using the lattice Boltzmann and finite volumes discretizations. One more thing we observed here is as Reynolds number increases, coefficient of drag decreases even for the vortex-shedding flows.

Table 2: Mean drag coefficients for the flow past a square cylinder at different Reynolds numbers.
Figure 8: Flow over a square cylinder. Time-dependent lift and drag coefficients ( and ) at (a) Re = 60 and (b) Re = 100 on a lattice size of 800 × 128.

In order to study the flow field characteristics, velocity variation along the channel centerline is presented for Re = 60. Figure 9 shows the instantaneous results of streamwise velocity (u) and cross-stream velocity (v) distribution along the channel centreline for a Reynolds number 60. It is observed that fluctuations of velocity at the exit boundary are small. At the surface of the cylinder velocity is zero, and at upstream locations flow velocity varies smoothly without much oscillation. Figure 10 shows velocity profiles of and at three different axial positions, x = 200, 300, and 500. From the figures we can clearly observe that velocity profiles, typical of the wake region, do not exhibit any orderly trend.

Figure 9: Instantaneous results at a certain moment for the flow past a square cylinder: (a) streamwise (u) and (b) cross-stream (v) velocities along the centerline (y = 64); Re = 60, lattice size: 800 × 128.
Figure 10: Instantaneous results at a certain moment for the flow past a square cylinder: (a) streamwise (u) and (b) cross-stream (v) velocities at three different positions in the flow field, center of cylinder (x = 200), near-wake (x = 300), and far-wake (x = 500) at Re = 60 and lattice size of 800 × 128.
##### 3.2. Different Blockage Ratios

To study the effect of different blockage ratios Re = 51, cylinder location l = L/4, and a channel length of L/D = 50 are considered. Figure 11 shows the vorticity contours for confined flow over a square cylinder at different blockage ratios. Five different blockage ratios B = 8, 10, 15, 20, and 25 are studied in the present work. From our LBM simulations we observe that the moderate Blockage ratio B = 10 onwards unsteadiness develops in the flow field for Re = 51. We can safely say that this will happen also for Re > 51. The flow appears to be steady in the low blockage ratios. It can be clearly observed that at low blockage ratios wall boundary layer has strong influence over around the cylinder so that it inhibits vortex shedding. As the blockage ratio increases computation, time required for developing unsteadiness increases.

Figure 11: Instantaneous vorticity contours for flows past a square cylinder at different blockage ratios: (a) B = 8, (b) B = 10, (c) B = 15, (d) B = 20, and (e) B = 25.
##### 3.3. Different Cylinder Locations

To study the effect of different cylinder locations Re = 51 is chosen in the present work. As usual a blockage ratio of B = 8 and L/D = 50 are considered. All the cases studied so far take l = 12.5D. In the present case, l = 25.0D, 16.0D, 12.5D, 10.0D, and 8.5D are considered. Figure 12 shows the instantaneous vorticity contours for confined flow over a square cylinder for various cylinder locations. It is seen that, for a Reynolds number (here 51), if the cylinder is closer to the inflow boundary, it is more likely to develop unsteadiness. In Figure 12, except for case (a), all the other cases show unsteadiness. For the unsteady cases it is observed that, nearer the cylinder to the entrance, the lower is the time for unsteadiness to develop.

Figure 12: Instantaneous vorticity contours for flows past a square cylinder at different cylinder locations for B = 8, Re = 51: (a) l = 25.0D, (b) l = 16.0D, (c) l = 12.5D, (d) l = 10.0D, and (e) l = 8.5D; lattice size of 500 × 80.
##### 3.4. Different Outlet Boundary Locations

To study the effect of different outlet boundary locations Re = 51 is taken. In this case, blockage ratio B = 8 and the upstream location l = L/4 are chosen. Outflow boundary locations given by L/D = 20, 30, 40, and 50 are considered. Figure 13 shows the instantaneous vorticity contours for confined flow past a square cylinder at different outlet boundary locations. It is clearly observed from Figure 13 that the outlet boundary location is also playing important role in the vortex shedding process for a particular Reynolds number. For the extrapolation boundary condition for the particle distribution function at the outflow boundary to hold accurately, the boundary should be far enough downstream of the cylinder. In our present study, results reveal the fact that L/D = 30 onwards periodic solutions are obtained for Re = 51.

Figure 13: Instantaneous vorticity contours for flows past a square cylinder at different outlet boundary locations for Re = 51, l = 12.5D: (a) L/D = 20, (b) L/D = 30, (c) L/D = 40, and (d) L/D = 50.

#### 4. Conclusion

A lack of accurate LBM computations was found in the literature for the confined flow past square cylinder. From the available literature of numerical methods, it is observed that when confined walls are present in the flow, the behavior of vortex shedding behind the body is affected appreciably. Therefore, we presented accurate and reliable LBM computations in the present work. In this work to understand the flow behaviour of the problem for the square cylinder, we investigate four different cases. First, we examine the flow field characteristics at different Reynolds numbers. Our LBM computations compare well with existing results. For a blockage ratio (B = 8) we found that periodicity starts at Re = 52. Next, we study different blockage ratios at Re = 52. We conclude that, at low blockage ratios, wall boundary layer has strong influence over around the cylinder so that it inhibits vortex shedding. Then, we study different upstream locations. We observe that as cylinder is located nearer to the entrance, periodicity develops faster. Next, we examined different outlet boundary locations. It is seen that as the length of the downstream portion of the channel increases, the flow gets periodic behaviour faster. The present computations demonstrate the capability of lattice Boltzmann method to capture vortex shedding phenomena and efficient handling of instantaneous flows.

#### References

1. R. W. Davis, E. F. Moore, and L. P. Purtell, “A numerical-experimental study of confined flow around rectangular cylinders,” Physics of Fluids, vol. 27, no. 1, pp. 46–59, 1984.
2. R. Franke, W. Rodi, and B. Schönung, “Numerical calculation of laminar vortex shedding flow past cylinders,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 35, pp. 237–257, 1990.
3. A. Mukhopadhyay, G. Biswas, and T. Sundararajan, “Numerical investigation of confined wakes behind a square cylinder in a channel,” International Journal for Numerical Methods in Fluids, vol. 14, no. 12, pp. 1473–1484, 1992.
4. H. Suzuki, K. Fukutani, T. Takishita, and K. Suzuki, “Unsteady flow in a channel obstructed by a square rod (crisscross motion of vortex),” International Journal of Heat Fluid Flow, vol. 14, no. 1, pp. 2–9, 1993.
5. A. Sohankar, C. Norberg, and L. Davidson, “Numerical simulation of unsteady low-reynolds number flow around rectangular cylinders at incidence,” Journal of Wind Engineering and Industrial Aerodynamics, vol. 69, pp. 189–201, 1997.
6. A. Sohankar, C. Norberg, and L. Davidson, “Low-Reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet 21 boundary condition,” International Journal for Numerical Methods in Fuids, vol. 26, pp. 39–56, 1998.
7. A. N. Pavlov, S. S. Sazhin, R. P. Fedorenko, and M. R. Heikal, “A conservative finite difference method and its application for the analysis of a transient flow around a square prism,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 10, no. 1, pp. 6–46, 2000.
8. M. Breuer, J. Bernsdorf, T. Zeiser, and F. Durst, “Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 21, no. 2, pp. 186–196, 2000.
9. D. C. Wan, B. S. V. Patnaik, and G. W. Wei, “Discrete singular convolution—finite subdomain method for the solution of incompressible viscous flows,” Journal of Computational Physics, vol. 180, no. 1, pp. 229–255, 2002.
10. A. Roy and G. Bandyopadhyay, “Numerical investigation of confined flow past a square cylinder placed in a channel,” in Proceedings of the All-India Seminar on Aircrafts and Trans-atmospheric Vehicles: Missions, Challenges and Perspectives, pp. 28–29, Kolkata, May 2004.
11. S. Abide and S. Viazzo, “A 2D compact fourth-order projection decomposition method,” Journal of Computational Physics, vol. 206, no. 1, pp. 252–276, 2005.
12. N. Hasan, S. F. Anwer, and S. Sanghi, “On the outflow boundary condition for external incompressible flows: a new approach,” Journal of Computational Physics, vol. 206, no. 2, pp. 661–683, 2005.
13. K. M. Kelkar and S. V. Patankar, “Numerical prediction of vortex shedding behind a square cylinder,” International Journal for Numerical Methods in Fluids, vol. 14, no. 3, pp. 327–341, 1992.
14. S. Hou, Q. Zou, S. Chen, G. Doolen, and A. C. Cogley, “Simulation of cavity flow by the Lattice Boltzmann method,” Journal of Computational Physics, vol. 118, pp. 329–347, 1995.
15. D. A. Perumal and A. K. Dass, “Multiplicity of steady solutions in two-dimensions lid-driven cavity flows by Lattice Boltzmann method,” Computers & Mathematics with Applications, vol. 61, no. 12, pp. 3711–3721, 2011.
16. D. Yu, R. Mei, L. S. Luo, and W. Shyy, “Viscous flow computations with the method of Lattice Boltzmann equation,” Progress in Aerospace Sciences, vol. 39, no. 5, pp. 329–367, 2003.
17. M. Bouzidi, M. Firdaouss, and P. Lallamand, “Momentum transfer of a Lattice Boltzmann fluid with boundaries,” Physics of Fluids, vol. 13, no. 11, pp. 3452–3459, 2001.
18. G. V. S. Kumar, D. A. Perumal, and A. K. Dass, “Numerical simulation of viscous flow over a circular cylinder using Lattice Boltzmann method,” in Proceedings of the International Conference on Fluid Mechanics and Fluid Power, IIT Madras, India, December 2010.
19. D. J. Tritton, “Experiments on the flow past a circular cylinder at low reynolds numbers,” Journal of Fluid Mechanics, vol. 6, no. 4, pp. 547–567, 1959.
20. B. A. Fornberg, “A numerical study of steady viscous flow past a circular cylinder,” Journal of Fluid Mechanics, vol. 98, no. 4, pp. 819–855, 1980.
21. D. Calhoun, “A Cartesian grid method for solving the two-dimensional streamfunction-vorticity equations in irregular regions,” Journal of Computational Physics, vol. 176, no. 2, pp. 231–275, 2002.