A lattice Boltzmann study of 2D steady and unsteady flows around a confined cylinder

In this work, the lattice Boltzmann (LB) method was applied to simulate incompressible steady and unsteady low Reynolds number (Re) flows around a confined cylinder. In the LB method, different collision models (Bhatnagar–Gross–Krook model, two-relaxation-time model, multi-relaxation-time model, and entropic lattice Boltzmann model) and a regularization model were used, and the results were compared. Numerical results pertaining to a two-dimensional flow around a cylinder are reported and compared with numerical and experimental data available in the literature. The results agree with the predictions made from the literature. A correlation for Strouhal number (St) for 55 ≤ Re ≤ 300 is suggested.


Introduction
The lattice Boltzmann (LB) method is a numerical tool based around a discrete form of statistical-mechanical Boltzmann equation [1]. One of LB's primary uses is hydrodynamics computations. Several collision models have been proposed in the past, that include the Bhatnagar-Gross-Krook model-BGK [2], the two-relaxation-time model-TRT [3], the multi-relaxation-time model-MRT [4], the entropic lattice Boltzmann model-ELB [5], and, although technically not a collision model, the regularized lattice Boltzmann model-RLB [6]. By applying the right parameters, all of these models can reproduce the Navier-Stokes equations. However, each collision model has a specific window of application, which for some is wider, for others narrower. The "better" collision models usually carry a greater computational cost, and different problems require different approaches in terms of collision model choice. Different collision models have been put head-to-head before, e.g., their performance in finding a steady-state solution of a 2D lid-driven cavity flow benchmark was compared [7]; however, the ELB implementation therein received a follow-up comment [8].
Another popular computational fluid dynamics (CFD) benchmark is fluid flow around a cylinder [9]. The flow around a cylinder is often considered for testing different solid boundary conditions. Curved walls of the cylinder represent a challenge to be modeled properly, especially in LB computations, as these are usually performed on square lattices [10][11][12][13][14][15][16]. This benchmark serves well in comparing different CFD approaches. However, flow around a cylinder is also a good starting point for testing the solver's performance in solving flow through complex geometries [17,18]. With the onset of an unsteady solution at Reynolds numbers (Re) above its critical value (Re crit ), where a vortex street (von Kármán vortex street) appears, flow around a cylinder is also popular for testing the performance in solving unsteady flows [19][20][21] and turbulent flows [22].
The aim of this work is to: first verify and validate an in-house-developed LB solver, utilizing the BGK, TRT, MRT, ELB, and RLB models, via the CFD benchmark for flow around a cylinder in both states-steady and unsteady [9]; and second to study the flow around a confined cylinder using the verified model. The first part is also meant to Technical Editor: Edson José Soares, Ph.D.
compare the different models' performance in solving the benchmark problem. In the second part Re crit , as well as the correlation between the Strouhal number (St) and Re in the laminar and transition vortex shedding modes [23], will be investigated.

The lattice Boltzmann method
The LB method is performed on a square lattice where at each node there is a discrete number of directions in which the fluid particles can move. In this study, a 2D model with nine lattice velocities-the D2Q9, will be discussed and applied.
A general form of the LB equation can be written as follows [1]: where f i is a discrete distribution function at position x , with average particle lattice velocity e i , pointing in direction i (see Fig. 1), at time t. coll is the collision operator, and will be discussed later. e i = (e x,i , e y,i ) is the directional vector and it takes the following values: The local fluid density is computed as the sum of all f i 's on site. The expression for D2Q9 lattice model is the following: The local fluid velocity u = u x , u y can be computed from the expression for local momentum:

Bhatnagar-Gross-Krook model
The Bhatnagar-Gross-Krook (BGK) model has the simplest form of the presented collision operators, which reads as follows [1,2]: f eq i is the local equilibrium distribution function, and is computed from the local macroscopic velocity u , and local density : where w i is the equilibrium weight factor in the ith direction. Its values add up to 1, and in D2Q9 they equal: is the relaxation rate and is the inverse of the relaxation time = 1 . From the relaxation rate, the fluid's kinematic viscosity can be calculated: + and − are relaxation rates of the model. They are linked via the "magic parameter" : For almost all of the benchmarking, = 1 4 was used as this was proven to be its most stable value [24]. Indeed, this gave the best results, but it was numerically unstable in some cases. There = 10 −5 was used. The positive relaxation rate + plays an identical role in TRT as does in BGK, and therefore directly correlates with the kinematic viscosity of the fluid :

Multi-relaxation-time model
The multi-relaxation-time (MRT) model was developed almost simultaneously with the LB method itself [4]. The form of MRT presented and used in this study is of the Gram-Schmidt type [25,26]. The algorithm of MRT collision is as follows: the f i 's at a node in the population space are transformed into the moment space to moments m i 's via the transformation matrix M , where the collision is carried out through the relaxation matrix S , and finally m i 's are transformed back to f i 's via the inverse of the transformation matrix M −1 . In general, the MRT collision operator can be written as: where f = f 0 , f 1 , … , f 8 , and f eq = f eq 0 , f eq 1 , … , f eq 8 . The streaming step is carried out with f i 's in the population space. Gram-Schmidt M for the D2Q9 model takes the following form: The inverse matrix M −1 can be looked up in the literature [25,26] or it can be computed with the help of a computer solver. After the transformation to moment space, the coll i moments relax toward the equilibrium moments m eq i which are: The relaxation is done through S , which in Gram-Schmidt procedure takes the form of a diagonal matrix, and can thus be written as: where is of particular interest, as it correlates with the fluid's kinematic viscosity in the same way as in BGK: is correlated with the fluid's bulk viscosity. However, it, along with the other relaxation rates (except ) is usually set to a fixed value, which requires some tuning on the user's part to find the optimal values. For the purpose of this study, the following values were used S = diag(0, 1.95, 1.95, 0, 1.4, 0, 1.4, , ).

Entropic lattice Boltzmann model
The entropic lattice Boltzmann (ELB) model was derived by finding the entropy functions through which the Navier-Stokes equations could be recovered [5,27]. For ELB, the collision operation can be written as: At the first glance, Eq. 19 seems very similar to BGK equation (Eq. 5), and indeed, if = 2 it becomes the BGK collision operator, as = 1 2 . is therefore directly correlated with the fluid's kinematic viscosity : However, f eq in ELB take a slightly different form: .
To understand the parameter , the entropy function H should be considered first. In D2Q9, H is defined as: should then be numerically determined as the non-trivial solution of the following equation, using the f eq calculated with Eq. 21:

Regularized lattice Boltzmann model
The regularized lattice Boltzmann (RLB) model actually incorporates the BGK collision, but it includes a regularization of the f i 's before the collision step itself [6]. There are several related regularization procedures in the literature [6,28,29]. In this study, the high-order regularization with filtered central moments by Mattila et al. was used [29]. First, let there be a non-equilibrium distribution function f neq i introduced: Then by combining Eqs. 1, 5, and 24, the LB equation can be rearranged as follows: The high-order f eq i is computed as: The actual regularization step occurs by computing f neq i as: where:

Boundary conditions
As a cylinder (or its 2D projection-a circle) requires the use of curved boundaries, the multi-reflection (MR) boundary condition was used for solid walls [30]. Unlike other approaches to curved boundaries, MR formulation is independent of the distance of the wall from the wet node q [31]. This property makes the MR fairly easy to implement and does not greatly increase the overall computational costs, especially in its reduced form [31,32]: q takes a value between 0 and 1, and represents the relative distance between the actual wall position and the boundary node. The positional parameters B and N represent the boundary, and neighboring nodes, respectively, N = B + e i t . Straight walls had the value of q set to 0.5, which reduces Eq. 29 to the standard no-slip bounce-back boundary condition. The inlet boundary condition can be a simple Dirichlet boundary condition, where the multi-reflection boundary is slightly modified: where the u xin is the inlet velocity, a function of y-coordinate, as in computations in this study a parabolic inlet velocity profile was applied. in is the density at inlet, and it can be calculated as [33]: . For the outlet boundary conditions, the non-equilibrium boundary was utilized [34]. This one is formulated in the following manner: where the parameters for f eq i indicate the values that are used in its computation. 0 takes a user-defined value and represents the constant density (pressure) at the outlet.

Benchmarking the code
The code was first tested in a simple setup-steady and unsteady 2D flow around a cylinder. The benchmark is such as described by Schäfer et al. [9]. A 2D domain, which is bounded on the north and south sides by static no-slip walls, was designed. A cylinder was positioned slightly off-midstream (Fig. 2); the inlet boundary condition is a fully developed parabolic profile, and the outlet boundary condition can be chosen freely. In this study, the boundary conditions were chosen as described above. For MR, at straight walls q equaled 0.5.
The benchmarking is done at two different benchmark Reynolds numbers, Re bm : where ū x is the mean flow speed in the x-direction at the inlet, D is the cylinder diameter, and is the fluid's kinematic viscosity. The benchmark for the steady flow is set at Re bm = 20 , while the unsteady case is studied at Re bm = 100.
The tested values were: the drag coefficient c D , the lift coefficient c L , the recirculation length L a at Re bm = 20 , the Strouhal number St bm at Re bm = 100 , and the downstream pressure drop at the cylinder P: .
where F D and F L are the drag and lift components of the force acting on the cylinder ( F = F D , F L ), respectively, and is the fluid density. x r is the x-coordinate downstream from the cylinder, where u x equals 0, x e is the x-coordinate on the edge of the cylinder that is the furthest downstream, equaling 0.25, and x a lies on the opposing side of the cylinder upstream. f sh is the vortex shedding frequency, and it was determined as the frequency at which there was a maximum of the Fourier transformed data of c L .
In the steady case, the resulting c D and c L were the steady solutions, whereas in the unsteady case both parameters in question were their maximum values, c max D and c max L , when a periodic steady state was reached (the von Kármán vortex street). F D and F L were calculated as the change of momentum in a time period t in a given direction, also accounting for the curved wall position q [30]: B represents all solid cylinder boundary nodes, N represents B + e i t if the said point is fluid, and N + 1 represents B + 2e i t.

Simulation setup
Again a cylinder is put between two walls in a 2D system. Unlike in the benchmark, here a channel with a cylinder placed midstream and a fifth of the channel length downstream was set up. The channel width was 8 times the cylinder diameter D, and channel length was 4 times the channel width. The channel-width-to-cylinder-diameter of 8 (blockage ratio of 0.125) was found to reproduce the results of an unbounded cylinder [35]. These dimensions also made sure that the inlet and outlet boundary conditions did not interfere with the computations. The setup is represented in Fig. 3. The boundary conditions used were the same as above.
The computations logged the density difference through the middle of the cylinder in the y-direction. These measurements are proportional to the pressure drop in the same direction. By processing these data, the frequency of vortex shedding f sh could be obtained. The simulations were performed in a range of Reynolds numbers Re between 50 and 300. Re is defined as follows: The flow speed ū ⋆ is calculated as the average x-component of the velocity, u x , at inlet over the obstacle's projection onto the inlet [36], in this case: According to Chen et al., this formulation of average flow speed allows for comparison of bounded and unbounded cylinders [36]. Additionally, the critical Reynolds number Re crit was determined by increasing Re between 50 and 60 in increments of 1. Re crit is the Re at which an instability in flow occurs.
To compare the present results to the ones found in other published work, f sh was converted to the Strouhal number St. This was then plotted against Re. St is defined as: The simulations ran for a total of 200 000 time steps per value of Re. The output data were then transformed with the Fourier transformation to obtain the vortex shedding frequency spectra. The simulations for Re crit ran for 2 000 000 time steps, as it took longer for instabilities to develop at lower Re.
The cylinder diameter chosen for the simulations was D = 32 . Computations were performed with the MRT collision model.

Hardware and software
The computations were carried out with an in-housedeveloped LB solver, written in CUDA C++. The code ran on two separate machines, namely the benchmarking was performed on an Intel ® Core™ i7-6700HQ, 16 GB RAM, NVIDIA ® GeForce ® GTX 960M 4 GB laptop, and the simulations ran on an Intel ® Core™ i7-8700K, 8 GB RAM, NVIDIA ® GeForce ® GTX 1060 6 GB desktop computer. Both machines were using Ubuntu 18.04.1 LTS operating system, NVIDIA ® CUDA ® Compiler Driver (nvcc) from NVIDIA ® Toolkit version 9.1, and GNU C++ (g++) compiler version 7.3.0. Data analysis was performed with ParaView 5.4.1, and Python 3.6.5, using the numpy, scipy, and matplotlib libraries. The Fourier transform utilized was the scipy.fftpack.fft. To smooth out the datasets computed in the benchmark (especially for the c max L and c max D determination), the Savitzky-Golay filter (scipy.signal.savgol_filter) was applied. This reduced the errors significantly.

Benchmarks
The benchmarking results are given in Table 1 From these results, it is evident that the emerging flow pattern remains pretty much unaffected by the size of D (results for St bm ) and the choice of the collision operator. However, MRT and RLB were the most stable, but MRT required shorter computation times (data not shown), so further computations were performed using it. D = 32 appears to be a good enough resolution for further computations, as D = 64 does not greatly improve the results, and smaller    D would result in greater errors. All the presented models performed better in the unsteady (Re bm = 100 ) case, which is good for this study, as the flow in the transition range of Re that was investigated later does not give a steady-state solution. The overall benchmarking results suggest that the presented approach is suitable for simulating 2D fluid flow around a cylinder. Overall the computation time of the BGK, TRT, and MRT models was similar ( ±6% ), whereas ELB and RLB computation time was about twice as long (data not shown). It should be noted, however, that RLB implementation here was not as far optimized as the other models'.

Re crit determination
Blockage ratio of the channel has an effect on Re crit , which is the lowest Re at which a disturbance in the wake of the cylinder can be noticed. In general, Re crit increases with greater blockage ratio. Sahin and Owens investigated this phenomenon [37]; however, in that study Re crit was not determined for blockage ratio of 0.125, which was used here. By interpolating their data, it was determined that in this case Re crit should be about Re = 54.6. The results of the numerical experiments, which were used to determine Re crit of the investigated system, are displayed in Fig. 8. In Fig. 8a, the maximum peak size of the Fourier transform of the data log was plotted against Re. The peak size is affected by the frequency and the amplitude of the signal. A more intense (higher amplitude) wave produces a higher peak size. Because the disturbance appears sooner in flows with higher Re, only the last 150,000 logged points were used in the Fourier transform. At this point, the disturbances were already fully developed, and the comparison between the signals could be made. Figure 8b displays the u y at the channel mid-point along the x-axis after 2,000,000 time steps. The plots shown are for Re = 53, 54, 55. This shows that practically no change in flow pattern occurs between Re 53 (a) Maximum peak size plotted vs Re.
(b) u y vs x. Fig. 8 Re crit determination. a Maximum peak size plotted against Re. From this graph, it is evident that a disturbance in the flow becomes significant at Re = 55. Plotted is also a best-fit line across 55 ≤ Re ≤ 57. b u y at the channel mid-point plotted along the x-axis for Re = 53, 54, 55, all after 2 000 000 time steps. The lines for Re= 53, 54 appear identical, whereas a significant disturbance is seen at Re = 55  [38], and Norberg [39] are based on experimental results, whereas the results of Zhang [40], Behara [41], Mittal [42], and Singha [35] are from computer simulations and 54, but a significant disturbance appears at Re = 55. Both plots mentioned above suggest that Re crit = 55.0 ± 0.5 , which favors the prediction made above (Figs. 6, 7).

Correlation between St and Re
The correlation between St and Re was investigated and compared to some results found in the literature, which were both experimental and computational [23,[38][39][40][41][42].
The results of these simulations are presented in Fig. 9, and an example of a velocity profile of fully developed flow at Re = 250 is displayed in Fig. 10. The present results agree well with the others. A general trend of St increasing with higher Re can be observed. At lower Re, St increases faster, but the increase slows down at higher Re. In particular, the agreement is good between the present results and [38,39] in the range 55 = Re crit ≤ Re ≤ 100 . A function to fit the data was also calculated:

Conclusions
A brief overview of the LB collision models was presented, along with the benchmark model for 2D flow around a cylinder [9]. The presented LB collision models were applied to the benchmark, and were compared among one another. For the specific problem, no significant differences between the models were noticed, especially at higher grid resolutions, and all were able to reproduce the flow, as prescribed by the further apart (especially ELB). The poor performance of ELB in the benchmark is surprising. However, the method is said to have an edge over others, especially in high Re flows which exceed the range of the presented study [8]. ELB also does not require any parameter tuning, unlike MRT, which makes its use easier for the end user. Similarly, RLB also did not appear to have an edge over other methods, apart from great stability, but like ELB it involves no parameter tuning, and it is (in these authors' humble opinion) even simpler to implement. Also surprising is poor stability of TRT at Re bm = 100 that caused the need to use a different in this case. The instability developed from the inlet boundary condition (data not shown), which suggests its poor compatibility with TRT. Among the collision models, there is no clear-cut winner in this application. The flow around a confined cylinder was studied and compared to diverse data from the literature, mainly data for an unconfined cylinder. For the purpose of these computations, the MRT collision model was used. Re crit for the cylinder with the blockage ratio 0.125 was determined to be 55.0 ± 0.5 , which agrees with the earlier prediction from the literature. The results from the laminar and transition vortex shedding modes agree well with the literature. A correlation for St in 55 ≤ Re ≤ 300 was suggested. adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.