An overview of the Lattice Boltzmann Method has been presented with an in house algorithm for the numerical simulation of fluid flow around a circular cylinder. The linearization of the collision operator has been discussed for distributions not close to t he local equilibrium state, and numerical simulation has been carried out for stable initial conditions up to a Reynolds number of 80. An overview of the lattice gas automata with regard to Boolean variables describing the particle occupation has also been defined. A comparison between the data obtained from the two dimensional fluid flow around the cylinder and previous experimentation has also been made. Keywords: Lattice Boltzmann Method, Gas Automata, Computational Fluid Dynamics, Mesh Formation 1.0 INTROD UCTION The lattice Boltzmann method (LBM) has been widely used as a numerical method to simulate fluid flows in the field of computational fluid dynamics. It strains from the lattice gas (LG) automata method, which makes use of a discrete lattice and time schedule wherein particle velocity in the distribution function, along with both space and time have been discretized. The lattice Boltzmann method reflects the employing of the Boltzmann equation without the consideration of a Eulerian led approach. This is why the lattice Boltzmann equation can be introduced as a reduction of the discrete kinetic equation for the particle distribution function (Chen et al 1998) which can also be regarded as a direct transcription from the LG automata. The fundamental idea behind the evolution of the LG automata into the lattice Boltzmann method is to simplify kinetic models by employing the consideration that the collocated macroscopic properties of a fluid system are the result of the evolution of particle distributions. This consideration takes into account that macroscopic dynamics are not influenced by microscopic interactions considerably, and hence avoids the over -complication of the kinetic equation by adopting a mesoscopic kinetic distribution instead of a microsco pic distribution. Unlike previous molecular dynamics simulations that took into account the trajectories and coherent properties of each particle, the averaged macroscopic distribution allows the implementation of boundary conditions to relatively complex boundaries without the accumulation of stochastic noise. Before understanding the numerical flow around an arbitrary body using boundary conditions imposed by the LBM, a brief overview of the LG automata is imperative to be able to shift from discrete tim e behaviour to an algorithm involving continuous space and time. Two further models exist wherein particle velocity is discretized and space and time are continuous. These include the single speed particle velocity and the multi speed discrete velocity mod el (Inamuro and Sturtevant 1990). One of the first models that considered both space and time to be discretized was given by Hardy et al (1976). Fig 1: Computational grid for the numerical simulation of flow around the cylinder 2 1.1 Lattice Gas Automata The LG Automata is often affected by fluctuations, making the presence of stochastic noise a troubling factor whilst simulating nodal fu nctions and their inherent values. Since the LG automata involves discrete particle velocity, space as well as time, it is relatively simple to derive further, a model that rests on statistical classical mechanics, by deriving procedures and values that ar e averaged over a particular subspace. As explained by Chen & Doolen (1998), the lattice gas automata consists of a set of Boolean variables ni ( x, t )( i = 1, … , M) describing the particle occupation, where M is the number of directions of particle velocities at every node. As indicated in [1], the evolution equation of the lattice gas automata can be represented as follows: ni ( x + ei , t + 1 ) = ni ( x; t ) + Wi ( n ( x, t )) where ( i = 0, 1, … , M) The lattice Boltzmann method is mesoscopic in nature since it replaces the occupation variable with single distribution functions that address groups of particles acting in coherence in place of single particles, ther eby reducing computational complexity and corresponding statistical noise. From a theoretical perspective, the lattice Boltzmann equation is not dependent on the LG automata and can be derived directly from the continuous Boltzmann equation as shown by He & Luo (1997). 1.2 Use of Interpolation to Increase the Numerical Accuracy of the Lattice Boltzmann Method Setting of a hydrodynamic limit for the Boltzmann equation was necessary to derive the detailed function spaces of the macroscopic quantities. This can be idealistically carried out by applying the conservation principle to the given system that is being a nalysed. The solution to the Boltzmann equation is always continuous since the equation does not depend on any spatial variables. This was found by introducing a set of Lagrangian multipliers and then applying the conservation principles. As observed by Bhatnagar, Gross and Krook (BGK), the main effect of the collision term is to bring the velocity distribution function closer to the equilibrium distribution. Another adjustment was made by Higuera & Jimenez (1989) to linearize the collision operator by usin g this assumption set forth in the BGK model. The BGK model with regard to the lattice Boltzmann method involves two steps: relaxation and streaming. In the former, the particle distributions at nodal points settle to their equilibrium state, and in the s treaming process, particle distributions exhibit their inbound velocities and migrate to another lattice site. Since the nodal addresses of migration points cannot be determined beforehand, the new algorithm made use of an interpolation based strategy wher ein the probability of migration of particles at both nodal points, as well as intermediate points is determined. 2.0 METHODOLOGY 2.1 Setting up of the Computational Grid The setting of a hydrodynamic limit made the computational analysis through lattice Boltzmann method analytically efficient. To test this rudimentary set up of nodal grids, we consider a two dimensional cylinder bound in a rectangular boundary set up withi n the limits of incompressible flow. Since we wanted to establish a basis of simulation at the local equilibrium state, the two dimensional cylinder was simulated for a range of Reynolds numbers up to 80. The simulation starts with an initial state of unif orm flow, and simulation begins with the introduction of the cylinder with diameter D being suddenly introduced in the uniform flow. The development of initial acoustic waves was followed by the recirculation of the flow. The recirculation was seen to be i nitially very scattered however settled to a steady degree after about a thousand time steps. Fig 2: Schematic representation of the ( x, y) plane of the cylinder in a two dimensional Cartesian coordinate system The flow consideration outside the circular cylinder can be transformed into a rectangular strip as s hown by Fig 2 and Fig 5. Since the boundaries should all have a finite value for the numerical simulation to take place, the outer boundary in this case is assumed to be located at 𝝃 = 𝝃 ∞ . y a x θ x 3 Fig 3: The introduction of the body in a uniform stream wherein acoustic waves develop and begin to propagate throughout the fluid medium at a Reynolds number of 53.4 Fig 4: A periodic shedding flow begins to form at a Reynolds number of 77.8. The shedding intensifies and becomes more uniform as the number of time steps increases The computational mesh can be seen as a set of concentric circles spaced uniformly in the ( x, y) plane. The mesh considered has NX= 520 and NY = 180. There were two advantages in si mulating such a grid structure. First, the no slip boundary condition could be directly implemented on the cylinder surface and not on arbitrary lattice intersection points. Second, the density packing of the grid was considerable high, reducing the redund ancy of the algorithm. 2.2 Boundary and Initial Conditions Since the initial condition was symmetric, the lattice Boltzmann method application to the flow can be assumed to be much faster in comparison to pseudo -spectral methods. It was assumed that no symmetry breaking would arise at any given point during the nu merical simulation since there was no symmetry requirement with regard to the structure of the computational grid. The lattice Boltzmann method is heavily dependent on the initial definition of density distribution at every boundary point on the lattice. Three boundary conditions were defined, namely: the upstream axis, cylinder wall and the infinite flow filed vector considerations. 3.0 RESULTS AND DISCUSSION The simulation results were compared to previous studies such as Higuera & Succi (1989) and Noble et al (1996). The agreement of the simulation results was seen up to a R eynolds number of 60 as seen in Fig 3. A probable reason of variations at higher Reynolds numbers was found to be due to the difference in area covered by the nodal functions or the computational domain, versus that of the fluid flow region. This issue tha t could not be addressed despite successive iterations due to 4 Table 1 : Comparison of flow parameters with previous experimental and numerical simulation Re Authors L/a Α, 10 Dennis and Chang (1970) 0.53 29.6 Nieuwstadt and Keller (1973) 0.434 27.96 Countanceau and Bouard (1977a) 0.68 32.5 Fornberg (1980) - - He and Doolen (1997) 0.474 26.89 Present 0.498 27.05 20 Dennis and Chang (1970) 1.88 43.7 Nieuwstadt and Keller (1973) 1.786 43.37 Countanceau and Bouard (1977a) 1.86 44.8 Fornberg (1980) 1.82 - He and Doolen (1997) 1.842 42.96 Present 1.851 43.21 40 Dennis and Chang (1970) 4.69 53.8 Nieuwstadt and Keller (1973) 4.357 53.34 Countanceau and Bouard (1977a) 4.26 53.5 Fornberg (1980) 4.48 - He and Doolen (1997) 4.490 52.84 Present 4.487 51.29 Note : L: wake length. Α,: separation angle limitations in the computing power available. It was observed that the flow around the cylinder repeated after a certain number of time steps successively at a Reynolds number greater than 53.4 after a long transient phase. The incumbent Reynolds number of the flow has an effect on the pressure distribution over the cylinder. This was a phenomenon that was previously studied by Wagner (1994). The in house algorithm to implement the simulation was different in terms of Strouhal number from other simulations by 2.7%, and hence all results obtained were considered to be extremely accurate. To ensure high experimental accuracy, the no slip boundary layer condition was imposed on the cylinder boundary directly and not the jagged mesh points which could lead to er roneous values. One of the major criteria in evaluating any computational algorithm employing the LBM for numerical simulation is the definition of the boundary and the application of the no slip condition as explained above. The mesh is considered to be u niform across the computational domain, and hence the easiest method of implementation is to draw the boundary and mark all the links that the boundary line intersects and then defining the velocity index midway between two corresponding intersection point s. A periodic shedding flow was observed to commence at a Reynolds number of 77.8 as shown in Fig 4, in coherence with experimental results obtained by Tritton (1959) and Gerrard (1978). Fig 5: Schematic representation of the (𝝃, n) plane of the cylinder in a two dimensional plane 4.0 CONCLUSIONS Averaging properties on the microscopic scale for macroscopic visualization has been numerically presented using a bottom down approach that is often characterized by distributions close to the local equilibrium state. This generalization was then applied simulate the dynamics of low Reynolds number fluid flow around a two dimensional cylinder, the only added requirement being the increase in ή 1 0 𝝃 ∞ 𝝃 -1 5 Fig 6 : The building up of the transient which lasts for an order of around six cycles at a Reynolds number of 77.8 Fig 7 : The flow visualization one cycle before the formation of the initial Karman vortex sheet at a Reynolds number of 77.8 size of the computational domain using advances in parallel computing. The computational efficiency was greatly enhanced using the in house algorithm while leaving the basic grid structure same as that of previous simulations in literature. It was noted that the vortex shedding Reynolds number, although computationally obtained, differed as compared to current values. This could be due to difference in Strouhal number, as well as implementation of a no slip condition. The method of interpolation may also be applied to the immersed boundary method to analyse flow over dynamic flow systems.