FLOW DYNAMICS OF LID-DRIVEN CAVITIES WITH OBSTACLES OF VARIOUS SHAPES AND CONFIGURATIONS USING THE LATTICE BOLTZMANN METHOD

This work implements the emerging computational technique namely the Lattice Boltzmann Method (LBM) to a fluid flow problem of single sided lid-driven cavities with various shapes of obstacles placed in it. The numerical methodology employs the Single-Relaxation-Time (SRT) model applicable to low Mach number hydrodynamic problem for incompressible flow regime. Three geometrical shapes of the obstacles considered are circular, square, and elliptic. Cavity with obstacles exhibited remarkable circulation zones and structures in contrast to the classical lid driven cavity. The flow mechanics and the vortex dynamics are studied for various values of Reynolds Number (Re = 100, 400, and 1000). Due to the introduction of the obstacles, a strong induced vortex forms close to the obstacles and its size changes interestingly with the variation of Reynolds number, which is captured by LBM. Further the study is extended to examine the vortex phenomena induced by changing the position of the obstacles within the cavity. It is observed that the flow structures change dramatically with little change in the position of obstacle inside the cavity which helps to identify position with enhanced mixing characteristics.


INTRODUCTION
Fluid dynamics simulation using numerical methods is a very hot topic in the past few decades [1]. Closed recirculating flow generated inside bounded wall actuated by a moving lid represents an interesting problem in the field of fluid mechanics. This type of cavity problem is an exemplar problem commonly studied and expounded while developing any numerical code [2]. With tremendous advancements in the domain of computational technologies, Computational Fluid Dynamics (CFD) utilizing the computational resources is emerging to solve the governing equation of fluid mechanics, due to which there is manifold development in numerical algorithms. The lid driven cavity problem underpins various problems of academic and industrial importance such as in roll-coating and drying techniques, melt spinning process, cut-outs, cavities in surfaces of aeroplane bodies, heat exchangers etc. [3]. Additionally, its research-oriented importance is because it illustrates the formation of corner vortices and eddies, longitudinal vortices, non-uniqueness, and transition to turbulence.
Initially, CFD simulations of two-dimensional lid driven cavities were performed by Burgraff [4], and Pan & Acrivos [5] on square and rectangular shaped cavities. Precise computations were performed by Ghia et al. [6] when he implemented a stream function -vorticity formulation along with a coupled implicit multigrid method to accurately compute flows in square cavities at various Reynolds Number. In recent past many articles have been dedicated to exploring various other computational techniques to compute incompressible fluid flow in lid driven cavities. Some of the well-known techniques are Gradient Smoothing Method (GSM) [7], Finite Element Method (FEM) [8], Finite Volume Method (FVM) [9] and Smooth Particle Hydrodynamics [10].
Apart from the above listed continuum-based methods, researchers have adopted the Lattice Boltzmann Method (LBM) which is a discrete computational method based on the Boltzmann equation. The LBM has proved its mettle in playing a crucial role as an essential simulation tool for understanding micro-and macro fluid flows [11]. The popularity that LBM has achieved is due to its simplicity, numerical stability, and high efficiency for 84 implementation in parallel computers. Thus, this has led to many articles on successful numerical simulation using LBM in square cavities [12], rectangular cavities [13]. By now it is clear that many of the researchers in the past carried out numerical investigation that are directed only on studying the lid-driven cavity with actuation by one or two lids.
In general, literature on lid driven cavity with obstacles is sparse. Oztop et al. [14] has carried out mixed convection heat transfer simulation of lid-driven air flow with different values of Reynolds Number such as 100, 500, and 1000 using SIMPLE algorithm. They observed that the block of obstacle placed inside the flow plays an important role in controlling the fluid flow and heat transfer characteristics. Few researchers have targeted heat transfer problems in cavities involving circular obstacles in general. Rahman et al. [15] has carried out the numerical analysis of conjugate joule heating with magnetic field in a lid driven cavity with circular obstruction using the Galerkin Finite Element formulation.
Gangawane [16] numerically studied the mixed convection and heat transfer analysis in square lid-driven cavity containing heated triangular block with constant heat flux condition. Then, Gangawane et al. [17] examined the influence of position of the heated triangular block along the vertical centreline of cavity on mixed convective characteristics. They assumed flow is steady, 2D, laminar and for Newtonian and incompressible. Numerical investigation of double-sided lid-driven cubical cavity induced by a cylindrical shape at centre by FVM using multigrid acceleration [18]. The pressure and velocity of LBM boundary conditions are discussed by Zou and He [19]. Abbassi et al. [20] numerically studied natural convection in an incinerator shaped domain with a localized heated source situated at the bottom. Hussein and Hussain [21] presented for mixed convection flow of air within a parallel motion two sided lid-driven parallelogrammic cavity in the presence of magnetic field. Taghikhani [22] investigated magnetic field effect on the heat transfer in a nanofluid filled lid driven cavity with Joule heating.
Based on the literature review, it can be concluded that most numerical studies involving obstacles in lid driven cavity flow have focused on solving problems only with circular obstacle using the continuum-based methods. Further, the inclusion of different shaped obstacles such as ellipse and square results in interesting feature of flow like induced vortices, bifurcating vortices, recirculation adjacent to the obstacles. These flow topologies depend largely on the shape and orientation of the obstacles placed inside the cavity. Variation in the position of the obstacles helps to identify the optimum region for enhanced fluid mixing and improved heat transfer coefficients.
The current work is organized into four sections. Section 2 talks about the numerical methodology with algorithm, its validation and the problem description. Section 3 discusses the results for each obstacle and the vortex positions. Streamline plots, vorticity contours, and centerline velocity profiles are presented. Section 4 concludes the work followed by acknowledgments and conclusions.

Lattice Boltzmann Formulation
In this present work, the Single-Relaxation-Time (SRT) model obtained by applying the Chapman -Enskog expansion is used. In kinetic theory treatment this discretized equation is formulated as [11], Here, fi is the particle distribution function, ci is the particle velocity in the ith direction. fi eq (x,t) is the equilibrium distribution function at x, t and τ is the relaxation time, For an incompressible flow, the relaxation time is calculated in relation with the viscosity of the fluid based on the continuum assumption. Since the problem involved is 2D, the lattice model chosen is D2Q9 square lattice. It is known that, compared with D2Q7 hexagonal lattice D2Q9 square lattice gives better results. It is as shown in Figure 1 having nine discrete velocities. In this lattice model, each node has eight neighbors connected by eight links. During the streaming process, particles residing on a node move to a neighboring lattice along these links in a time step.
A suitable equilibrium function has been proposed by Lim et al [11] as, After solving for discrete Boltzmann equation on lattices, the macroscopic properties are found out by taking the moments of the distribution functions. Density and momentum are defined and evaluated as: Here stands for the speed of sound and is taken as = √ . Moreover, the relation between the relaxation time and the kinematic viscosity defined as, Treatment and incorporation of boundary conditions play a very important role in any numerical solution. To elucidate the treatment of the boundary conditions using the interaction of particles at the bottom wall as shown in the Figure 2. Dirichlet boundary conditions are implemented by equating the lattice velocities at the moving wall to the fixed lid velocities. During initialization, at the stationary wall, both the velocity components are initialized to zero. Second order accurate bounce back boundary conditions are also used to model the no-slip boundary condition on the wall [11]. The boundary conditions given by Zou & He [19] are applied to the moving walls as the velocity is known already. In case of the top wall, the unknown particle distribution functions are computed as: Here, , , and stands for density, v -velocity, and u -velocity at the bottom wall respectively.

Lattice Boltzmann Method Algorithm
The constituents necessary for a complete LBM algorithm are: an initialization process, the lattice Boltzmann equation (Equation 1), the velocity moment equations (Equation 5), and appropriate boundary conditions. After these constituents are specified the entire algorithm can be expounded. The flowchart for the LBM algorithm used to solve the problem is illustrated in Figure 3. A numerical code to implement the algorithm was developed using C++ programming language and thus the simulation results were obtained. A maximum lattice lid velocity was taken as U = 0.1. During initialization, the density, velocity and distribution function were initialized with the moving lid conditions and equilibrium distribution function (Equation 3) obtained from the initial conditions. The developed code consists of the steps illustrated in the flowchart. After the convergence criteria is satisfied, streamline patterns, vorticity contours, pressure contours and centerline velocity profiles were obtained.

Numerical Code Validation
For validation of the developed numerical code, a square lid-driven cavity is simulated. After solution, the streamline patterns and centerline velocity profiles were plotted for Re = 100. Figure 4 shows the schematic and the non-dimensionalized boundary condition imposed on the lid driven cavity. The streamline pattern is shown in Figure 5. At this Reynolds number, the formation of a strong primary vortex happens at the center of the cavity, along with two secondary vortices forming at the two bottom corners of the cavity. Figure 6 depicts the centerline velocity profile and its validation with the benchmark numerical results of Ghia et al. [6]. The number of lattices in both x and y direction are taken to be 129×129 as suggested by Ghia et al. [6].   The extension that has been done is the placement of obstacles of different shapes such as circular, elliptic obstacles at two orientations, square obstacles. All the obstacles are placed at the center of the cavity.
In the case of the elliptic cavity two orientations are considered. Firstly, the major axis of the ellipse is aligned with the direction of motion of the lid (α=0). Secondly, the major axis is aligned perpendicular to the motion of the lid (α=90). In the present study, the diameter of the circular is taken to be d = 0.4L. Here L is the length of the . Here U is the lid velocity, L is length of the cavity and the kinematic viscosity is taken as ν.

Lattice Convergence Test
To choose the optimum lattice (grid) resolution, lattice convergence test should be performed. So, the lattice convergence tests are performed for three lattice sizes ranging from 150 × 150, 300 × 300, and 450 × 450 with the circular obstacle at the center of the cavity with Re = 400. The location of the center of the vortices generated in the three grid resolution is enlisted in Table 1. From the table it can be observed that there is insignificant improvement in the vortex positions as the lattice size is increased from 300 × 300 to 450 × 450. Hence, an optimum lattice size of 300 × 300 is considered for further simulations, as lattice size of 450 × 450 would consume enormous computational time.

90
To specify the convergence criteria to stop the simulation after the convergence has been reached, it is a crucial step to find the errors after every finite number of iterations. To compute the residual error, a modified relative L 2 error is found which considers the effect of both the velocity components. The error norm is given by: Here, u n and u n-1 refer to the new and old horizontal velocity components respectively after every iteration to ascertain the residual error and check the convergence. Steady state simulations have been obtained with the criterion: Error ≤ 10 -6 .

RESULTS AND DISCUSSION
To study the fluid flow dynamics generated by actuation of a lid on a square cavity with obstacles, we use the present LBM model which also establishes its credibility to effectively model curved boundaries. It is thus applied to four pertinent problems: Cavity with circular, flat ellipse (α=0), vertical ellipse (α=90), and a square obstacle respectively. Later the code is extended to study the interesting flow features generated because of various position of obstacles in the cavity along the center-lines in both vertical and horizontal directions. To expound deeply on the flow dynamics generated the plots of streamline patterns and vorticity contours are depicted. Streamlines show the direction of the movement of the fluid particle. The streamlines are tangent to the resultant velocity vector at every point in the flow domain and thus it gives the sense of the motion of the fluid. Vorticity contours shows the tendency of the fluid to rotate. It is defined as the curl of the velocity vector. For a two-dimensional problem, it is given as, In the two-dimensional case, the curl is directed into or out of the plane. This phenomenon can be readily observed in Figure 8. This figure shows the streamline patterns and vorticity contours for four cases -Re = 10, 100, 400, and 1000 employing the SRT model of LBM. In case of Re = 10, an interesting phenomenon is observed. The induced vortex is divided into two sub vortices. Further it is also noticeable that both the sub vortices have a vertical centro-symmetry. From Figure 8(b), in case of Re = 100, a primary vortex can be seen encapsulating the obstacle and taking the form of the obstacle close to it while the two sub vortices coalesce and form a single induced vortex. An elongated vortex is induced right above the obstacle. Additionally, two secondary vortices are also observed opposite to each other on the left and right sides of the cavity in an asymmetric fashion.

Circular Obstacle in Lid Driven Cavity
As the Reynolds number is increased to 400 (Fig 8(c)), the induced vortex becomes larger in size and move to the immediate vicinity of the obstacle. Again, the size of the induced vortex increases leading to a shift in the vortex center towards the direction of the lid movement and this also cause it to move downwards due to the smooth slope of the circular obstacle. It is also to be noticed that the size of the secondary vortices increases, and the right secondary vortex becomes larger that the left one. Additionally, they also move upward and toward the center of the cavity. Further increasing the Reynolds number to 1000 results in the induced vortex becoming smaller due to the high velocity of the fluid actuated by the fast-moving lid. Remarkably, the secondary vortices become much larger in size as compared to the previous cases.
The vorticity contours are also illustrated on the right side of Figure 8. The results show that the moving wall generates vorticity that cause the fluid to rotate inside the cavity which eventually diffuses inside, and this diffusion mechanism acts as the driving mechanism of the flow. It is observed that strong vorticity gradients are formed close to the top wall and this increases with the increase in the Reynolds number. When the Reynolds number increases the induced vortex core has no viscous motion which is deduced by the absence of vorticity gradients. This phenomenon is typically observed in such close bounded lid driven cavity flows [1]. But close to the walls, at high Re, these gradient slides the fluid past the wall which is as opposed in the case of low Reynolds number. The vorticity contour at Re = 1000 matches well with the results given by Oztop [14]. Oztop has also observed the formation of such huge secondary vortices close to the bottom corner of the cavity. Table 2 presents the locations of vortices obtained using LBM for the circular obstacle for Re = 100, 400, and 1000. It is observed that, as the Reynolds number increases the primary and secondary vortex move towards the geometric center of the cavity.   Table 3 for both the elliptic configurations. It is observed that the position of vortices changes with position of obstacles. Furthermore, the growth of the two secondary vortices is more predominant and faster than the case of the circular obstacle. On the other hand, vertical ellipse is on the other side of the coin. The elongated portion plays an important role in controlling the size of the vortex. The vortex is constricted and thus as the Reynolds number increases there is a larger tendency for the induced vortex to fall due to the high slope of the ellipse. Moreover, it is observed that the size and position of the vortex can be controlled by the shape of the obstacle. Thus, in the cases of mixing problem and in heat transfer enhancements in electronic cooling, such geometries entail in enhancing and suppressing the induced vortex which controls the process parameters of the system.  Figure 11. Again, the phenomenon of vortex formation is quite like that of the circular obstacle. An interesting feature of the square geometry is that the sharp edge causes the induced vortex to divide. As the Reynolds number increases from 100 to 1000, the induced vortex, due to presence of the sharp corner, splits into two induced vortices as observable in Figure 11(c). It is clearly deductible that geometries of sharp and smooth features can be used effectively to control the flow. Such features when used in cooling systems play a crucial role in heat transfer. Actuating the opposite side in case of Re = 1000 will lead to two more induced vortex on the size exactly opposite about the center of the cavity thus resulting in effective flow control. The position of the vortex cores inside the cavity is given in Table 4.

Comparison of Velocity Profile for Various Obstacles
To compare the various obstacles velocity profiles of lid-driven cavity is discussed in this section. Figure 12  (a), (b), (c), and (d) shows the variation of u -velocity and v -velocity along the vertical and horizontal center-lines inside the cavity. As the Reynolds number increase the magnitude of velocity above and below and on both sides of the cavity increases. This shows a measure of the strength of the vortex formed. It can also be noticed that as the dimension of the obstacle becomes larger in a specific direction, the fluid flow region gets constricted and hence the velocity increases to satisfy continuity. From this aspect, we can observe that when the Reynolds number increases the induced vortex becomes smaller and thus the mixing phenomena happening close to the obstacle is reduced despite the higher velocities attained. Hence, lower Reynolds number entails in better and thorough mixing due to large induced vortex. In this further study, the analysis is extended to study the nature of flow dynamics generated when the position of the obstacle is varied. For this study the square obstacle is chosen because of the sharp corner present which enhances control over the flow characteristics. The dimension of the square obstacle taken is l = 0.2L. Here nine position of the obstacle is considered. The obstacle is moved along the horizontal and vertical center-lines of the cavity. The Reynolds number considered is 400. The streamline plots showing the flow behavior is shown in Figure  13. From the figure, it can be observed that the flow generated is very interesting at different positions. When the obstacle is placed at the center in the bottom of the cavity (Fig 13(a)), a large primary vortex is generated which occupies a significant portion of the cavity. Due to presence of the obstacle, two secondary vortices are generated. This results in enhanced mixing at the vicinity of the cavity. As the obstacle is moved up one step further (Fig 13(b)), the formation of large secondary vortices vanishes and only one vortex forms at the bottom right corner of the cavity. Also, the size of the induced vortex has reduced. Further, when the obstacle is moved to the center of the cavity (Fig  13(c)), the flow structure is same as that discussed in previous section. Since the size of the obstacle has decreased, the size of induced vortex is larger. As the obstacle is moved further up (Fig 13(d)), the induced vortex attains the lateral shape of the obstacle as it is pushed downwards due to strong momentum transfer from the lid. It can also be observed that the size of the secondary vortices has increases as compared to the previous case. As the obstacle is moved to the topmost center of the cavity (Fig 13(e)), the flow structure is quite like inverted image of Fig 13(a). The secondary vortex created is enormous in this configuration, while the primary vortex is much smaller. The mixing is this case is much that in first case because of the strong vortex regions at the vicinity of the obstacle.
In the second category, the obstacles are moved along the horizontal centerline. In the configuration shown in Fig 13(f), the presence of the obstacle results in a large primary vortex accompanied with three secondary vortices. As the obstacle is moved rightward (Fig 13(g)), an interesting phenomenon happens -two tertiary vortices are formed. Further when the obstacle is moved (Fig 13(h)), a tertiary vortex is formed right below the obstacle. As the obstacle touches the right wall (Fig 13(i)), two tertiary vortices are formed. And the secondary vortex is enormous in size. The interesting features of all this configuration is that the flow features change quite drastically when the obstacle is moved by a very little distance. This study should help choosing the optimum position of the tube or pipes in electronic cooling or nuclear reactor for maximum heat transfer.

CONCLUSIONS
The fluid flow dynamics in lid driven cavity with various shapes of obstacles -circular, flat ellipse, vertical ellipse, and square -in terms of Reynolds number have been studied in detail using the Single Relaxation Time (SRT) model of Lattice Boltzmann Method (LBM). The streamline patterns, vorticity contours and the velocity profiles in horizontal and vertical centerline show that obstructed cavity exhibits remarkable circulation zones and structures in contrast to the conventional lid driven cavity. The following conclusions can be made from the present work: • The size of the induced vortex reduces when the Reynolds number of the flow increases, and the secondary vortices increase in size. • Sharp edges present in the cavity helps in bifurcating the induced vortex at higher Reynolds number. This to effectively suppress and enhance the vortex structures formed at the obstacle's vicinity. • As the position of the obstacle is moved inside the cavity a tremendous change in the flow structures take place with small changes in position of the obstacle. • Our future work is to extend the present study in deep and shallow lid driven cavities.