The main dependent variable of the LBM is the fluid particle density distribution function, f, whose transport is governed by the lattice Boltzmann equation

where, α is the directional index of the discrete velocity vector e_α , x is the position of the fluid particle, t is the time, dt is the time step, and τ is the non-dimensional relaxation time. f_α and are the discretized local density distribution function and its equilibrium, respectively. Equation  (1) is the linear governing equation of the standard LBM. The left-hand side of Eq.  (1) represents the streaming motion of fluid particles, whereas the right-hand side is the BGK collision term.^{[1, 19]}

is a function of the local fluid density ρ , velocity u, and sound speed c_s (or temperature):

which is derived from the Maxwellian distribution function. The weighting functions w_α are determined by satisfying necessary conditions of isotropy of the fourth-order tensor of velocities and Galilean invariance.^[1]

The macroscopic (continuum) fluid properties such as density ρ , velocity u, and pressure p are functions of f_α according to

The relaxation time τ is a function of the fluid density ρ , the kinematic viscosity ν , the sound speed c_s, and the time step dt, according to

The algorithm of the DMD used here is from Ref.  [13]. A brief description is given in the following part. First the data, which are produced by the LBM simulation and separated in time by a constant step △ t_s, is represented in the form of a snapshot sequence and given by a matrix It is assumed that the linear map A that maps a given v_j onto the following one v_{j+ 1} is

Therefore eigenvalues of A give the growth rate and frequency of each mode. Snapshots become linearly independent with sufficient number. This can be written in the vector form as

where r is the residual. By using singular value decomposition (SVD), S can be solved with minimizing residual r. So the Ritz values of S are used to approximate those of A, and the Ritz vectors of S are used to construct the modes as shown in Ref.  [13].

In this work, flows around a cylinder especially at Reynolds number, Re= 50, were carefully studied, while other simulations for Re ∈ [40, 160] were also carried out for validation. Therein, Re is calculated with the inlet velocity U, diameter of the cylinder D, and the viscosity ν according to Re = U × D/ν . A rectangle computation domain with size L_x × L_y = 61D × 50D was used. The cylinder was placed in the domain with its center located at x = 16D and y = 25D. The wall of the cylinder was considered to be non-slip and treated with the bounce-back scheme. A slip boundary condition was imposed at the top and bottom boundaries. A constant velocity U was prescribed at the inlet (x = 0) boundary. A full-developed boundary was applied at the outlet boundary (x = L_x). The GPU-based LBM solver used here is the open source code, Sailfish.^[20]

As mentioned before, a GPU works with lower frequency than a CPU, but has access to fast memory and contains many more cores on a single device. In the present study, simulations were performed on from single to four Tesla K20M cards. There are 13 processors on every GPU card, and each processor has 196 cores. In other words, there are 2496 cores on each Tesla K20M card.

The performance of the code was tested with difference mesh size, block size, and number of GPU cards before productive simulations, and is presented in Fig.  1. In particular, the performance versus the number of lattices per cylinder diameter, N_D = 40, 60, and 80 are given in Fig.  1(a). The clock time used for every 1000 time steps is given for comparison. The lower clock time used the better numerical performance it has. It can be seen that with the increase of N_D, more and more clock times are needed. On the other hand, with the increase of the number of GPU cards, lower clock time is needed. N_D = 80 lattices are chosen in the following simulations.

Fig.  1.

	Figure Option ViewDownloadNew Window
	Fig.  1. Performance of the code. (a) Performance versus mesh resolution, (b) performance versus block size, and (c) scalability study.

Every LBM node in the code is processed in a separate GPU thread, and these threads are grouped into one dimensional blocks for data exchange.^[20] The block size is thus another tunable parameter that has an effect on the performance, as shown in Fig.  1(b). For a simple case, the larger number block size would result in better performance; for complex cases, such a relationship is dependent on the cases. In the present study, the best performance can be found when the block size is 64.

The performance of the code on the GPU is also compared with that on the CPU with another open source code of LBM, Palabos, ^[21] for a mesh with N_D = 80. The later run on a supercomputer Blue Waters, ^[22] which has 22640 Cray XE6 nodes (each containing 2 AMD Interlagos processors and 16 cores per processor) and 4224 Cray XK7 nodes (each containing one AMD Interlagos processor). As shown in Fig.  1(c), up to 10192 GPU cores (or 52 processors, 4 Tesla K20M cards) and 8192 CPU cores (or 512 processors) were used, respectively. Both codes have good scalability on corresponding computing platforms. Moreover, the performance of 10192 GPU cores is similar to that of about 800 CPU cores. Considering the difference of clock speeds (2.3  GHz for CPU and 706  MHz for GPU), such performance is acceptable. Anyway, one can see that the performance of several GPU cards is comparable with that of a traditional cluster or even supercomputer.

First of all, the resulting drag coefficient, C_D(= F_D/(0.5× ρ U²)), and Strouhal number, St(= fD/U), for different Re numbers in the range of [40, 160] are compared with those from references. It can be seen from Fig.  2(a) that the obtained results of C_D agree well with those from Zdravkovich’ s simulation.^[23] For St, the dash line in Fig.  2(b) is from the three-term fit by Willason and Brown, ^[24] where St is given by

for Re ∈ [50, 200]. Again, the obtained St via the LBM agrees well with the regression formula and the data from Ref.  [25] via the high order method. Especially for cases Re = 40 and Re = 50, which are around the bifurcation point Re_c = 46.6, the current solver can fully capture the different flow regime.

Fig.  2.

	Figure Option ViewDownloadNew Window
	Fig.  2. Comparison of predicted drag coefficient and Strouhal number with the results in references: (a) CD versus Re, (b) St versus Re.

In this subsection, we focus on the flow at Re = 50. The time history of lift force coefficient (C_L) on the cylinder is shown in Fig.  3. As in Refs.  [14] and [16], different time scales can be identified in different time intervals. In Fig.  3(a), there is an exponential growing of C_L in time interval [400, 800]. Within this interval, amplification, A, can be studied by means of Stuart– Landau amplitude equation^[14]

where σ _r in the first part of the right-hand side of Eq.  (7) is a linear growth rate, and equal to 0.012 in the present study as shown in Fig.  3(b). Besides, in the interval [1200, 1500] the flow ultimately approaches the saturated limit cycle. In the following part, the DMD analysis is carried out for these two different intervals, where the snapshots are from the LBM simulation with sampling interval Δ t_s = 0.3333  s.

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. Time history of lift coefficient (a) CL and (b) | CL| .

The Ritz values and energy spectrum of a data set within time interval [1200, 1500] are shown in Figs.  4(a) and 4(c). The Ritz values of the three most energetic ones, corresponding to modes 0 to 2, and even higher order transient modes lies on the unit circle in Fig.  4(a), indicating that the corresponding modes are neither growth nor decay. They are partitioned equally with angle θ relative to St in the way St = arctan(θ )/2π Δ t_s. Such a phenomenon was also found in Refs.  [14] and [16]. It can be interpolated that the flow’ s periodic orbit exhibits harmonic-like behavior. The results from the DMD are also compared with the Koopman eigenvalues of the NS equations from Ref.  [14] in Fig.  4(b). This is accomplished by the logarithmic mapping of the eigenvalues according to μ = log(λ )/Δ t_s. The values of μ for modes 0 to 2 and their complex conjugate are at axis y = 0, and the same as those from the Koopman analysis. Moreover, the ones for modes 3 to 5 can be observed in Fig.  4(b), while they cannot be seen clearly due to a visualization reason in Fig.  4(a). Horizontally, they are distributed in a similar way to the corresponding Koopman eigenvalues in Fig.  4(b).

Fig.  4.

	Figure Option ViewDownloadNew Window
	Fig.  4. Ritz values computed from the DMD on the LBM results for flow in time interval [1200, 1500]. (a) Ritz values, (b) Ritz values versus the Koopman eigenvalues, (c) energy spectrum.

All modes marked in Fig.  4 are shown in Fig.  5 in vorticity. Modes 0 to 5 correspond to the mean flow, global mode, second harmonic mode, shift mode, transient global mode, and transient second harmonic mode respectively in the Koopman analysis. Vortex structures in modes 0, 1, 3, and 4 are symmetric, which correspond to large-scale convection of fluid structures in the wake, while modes 2 and 5 are anti-symmetric and correspond to the alternative shedding of the Von Karman vortex street.

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. Several main modes for flow in time interval [1200, 1500]. Modes 0– 5 correspond to panels (a)– (f).

The DMD results of the exponential growing regime in interval [400, 800] are shown in Figs.  6 and 7. It can be seen from Fig.  6(a) that all the Ritz values lie outside the unit cycle. So in Fig.  6(b), all the corresponding values are positive. This is different from the results in Ref.  [14], in which Bagheri collected snapshots in the whole transit range including the limit cycle as well as the nonlinear growing part and stated that the results are sensitive to the snapshots in their paper and their collection is not reasonable. It can be interpolated, with the data set we selected here, that an almost linear dynamics is observed and so the growth rate should be positive. We also noticed that even if there are double the number of snapshots, the DMD results are the same. It can be seen in Fig.  6(b), that the growth rate of the most energetic growing mode, mode 1, is 0.012. This value is the same as the σ _r of the lift coefficient within the exponential interval. The growth rate of mode 2 is twice that of mode 1. Besides in Fig.  6(c), the energy spectrum is not as regular as in Fig.  4(c). This is due to the higher order non-linearity, which needs a cluster of complex modes to describe it.

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. Ritz values computed from the DMD on the LBM results for flow in time interval [400, 800]. Panel  (a) for Ritz values, (b) Ritz values versus the Koopman eigenvalues, (c) energy spectrum.

Mode 0, in Fig.  7(a), has a longer tail than the equivalent mode in the limit cycle in Fig.  5(a). This is because mode 0 is of the base flow, where the perturbation grows, rather than the mean flow of a fully developed flow past a cylinder. Mode 1 has a symmetrical vortex structure as shown in Fig.  7(b), like the equivalent mode in the limit cycle in Fig.  5. It is similar to the results of the linear stability analysis in Refs.  [26] and [27], in which the leading eigen-mode has the same symmetrical vortex structure. Since mode 1 has the same value of growth rate as that of C_L, this means within this interval the growing of C_L amplitude is dominated by large-scale convection of fluid structures in the wake. Mode 2 has an anti-symmetrical vortex structure as shown in Fig.  7(c). The amplitude of the vorticity in the vortex center in this mode keeps growing downstream, this indicates that the alternative shedding of the Von Karman vortex street starts from a convective instability.

Fig.  7.

	Figure Option ViewDownloadNew Window
	Fig.  7. Several main modes for flow in time interval [400, 800]. Panels (a)– (c) correspond to modes 0– 2.