1 Introduction

The lattice Boltzmann method (LBM) ^[1-3] has been proved to be an effective computational tool for simulating and modeling many complex fluid problems such as multiphase flows ^[4-5], porous media ^[6], and turbulent flows ^[7-12]. In contrast to the conventional numerical solution of macroscopic equation, i.e., the Navier-Stokes equation (NSE), the LBM represents the macroscopic averaged properties by the evolution of the statistical distribution of microscopic particles in term of the discrete kinetic theory. The main advantages of using LBM include easy implementation of boundary conditions, short codes, and natural parallelism. Thus, the LBM has gained more and more attention recently ^[13].

Although the simplest LBM that used a single relaxation time (SRT) parameter, based on the Bhatnagar Gross Krook (BGK) approximation, has successfully simulated various hydrodynamics problems, it has some weaknesses, especially, the low numerical instability which leads to an unstable pressure oscillation throughout the entire computational domain. Therefore, the multiple relaxation time (MRT) lattice Boltzmann equation (LBE) has been proposed for improving the numerical stability ^[14-17].

Generally, the LBM is used as a direct numerical simulation (DNS) method for low Reynolds number (Re < 10 000) flows because the required computer memory is limited ^[18]. With the MRT-LBE, cavity flows with high Reynolds numbers (up to 1×10⁶) have been simulated ^[17]. With a further higher Reynolds number, even the MRT-LBE method will become unstable, or the computation cost will become too high. An improvement is to add the large eddy simulation (LES) technique to the MRT-LBE model ^[7] so that the small scale turbulence will be mimicked using a suitable model, e.g., the Smagorinsky eddy viscosity model, but the large scale eddy can be directly simulated.

Because of the significant improvement of numerical stability of MRT-LBE and the use of LES to model small scale eddies but directly simulate large scale eddies, these two approaches have been combined together to simulate turbulent flows. Many good results show that the above approach (hereafter called the MRT-LBE-LES model) is able to provide more accurate predictions of turbulent flows ^{[8-12, 19]}.

However, the computation efficiency of an MRT-LBE-LES model for turbulent flows is still too low to solve any practical flow problems unless a parallel computation can be used to accelerate the simulation pace. Even though the parallel programs of LBM based on the message passing interface (MPI) standard which requires an expensive computer clusters with hundred cores are available ^[12], the initial investment is high, or the access is limited. Recently, the graphic process unit (GPU) has been developed into an inexpensive and novel massively parallel computing device with tremendous computational power and high memory bandwidth ^[20]. In particular, after NVIDIA introduced a general purpose parallel computing platform and programming model called CUDA in November 2006, the parallel programming on GPU becomes relatively easy for software developers familiar with a standard programming languages such as C ^[21]. Thus, a GPU with CUDA has been installed in a personal computer to accelerate many applications, such as computational finance, molecular dynamics, seismic processing, and computational fluid dynamics ^[22].

Since an LBM updates a node only with the nearest neighbor nodes, this method can take advantage of the massively parallel capability provided by GPU and CUDA. For example, the MRT-LBE models ^[23-24] have been accelerated by using a NVIDIA 8800 Ultra graphics card with CUDA for the simulation of the square array cylinders and the moving sphere in a pipe, respectively; the two-dimensional (2D) lattice Bhatnagar Gross Krook (LBGK) codes' ^[25] computation efficiency has been improved by using a NVIDIA GTX 280 CUDA-enabled graphics card. In this study, we develop a GPU-based MRT-LBE-LES model for simulating turbulent flows. In order to evaluate the computation efficiency of GPU-based MRT-LBE-LES program and the ability of MRT-LBE-LES model to simulate turbulent flows, 2D lid driven cavity flows for various Reynolds numbers, up to 10 7, are chosen as test cases. Application of this GPU-based MRT-LBE-LES model for a three-dimensional (3D) weak turbulent flow (with Re=1×10⁴) is also available ^[19].

This study is organized as follows. In Section 2, the LES with the Smagorinsky eddy viscosity model is extended to the MRT-LBE. In Section 3, the implementation details of the MRT-LBE-LES model with CUDA-GPU are introduced. In Section 4, the simulation results of 2D lid driven cavity flows with the help of CUDA-GPU parallel computing are presented, and conclusions are summarized in Section 5.

2 MRT-LBE with turbulence modeling 2.1 MRT-LBE

The MRT-LBE describes the fluid flow by using two processes: collision and stream processes of particle distribution function (PDF), f, on each lattice node, r_i, in which the collision step happens on the moment space M, and the stream step is executed in the velocity space V. Although the details are available in Ref.^[15], a summary is repeated here to make this a completed study.

The evolution equation of MRT lattice Boltzmann model that can also be called the generalized LBE or the moment method, for Q velocity directions on the D-dimensional lattice, can be written as

(1)

where the bold symbols denote the Q-dimensional column vectors, i.e., f=(f₀, f₁, ···, f_Q-1)^T, m=(m₀, m₁, ···, m_Q-1)^T, and the superscript T stands for the transposal operator. In Eq.(1), the left-hand side describes the particle stream from the lattice node r_i to the nearest location r_i+e_αδt at the particle velocities e_α along the α-direction with the time interval δt in the velocity space V. The right-hand side of Eq.(1) represents the particle collision step, i.e., the relaxation process, executed in the moment space M, in which the moment vector m(r_i, t) returns to the local equilibrium value m^(eq) (r_i, t) at a rate determined by the diagonal relaxation matrix S, and then, the transformation matrix M (with size Q×Q) linearly transforms the PDF vector f to its moment vector m, that is, m=M· r and r=M^-1 · m.

Here, the MRT-LBE model with nine (i.e., α=0, 1, ···, 8) discrete particle velocities in the 2D lattice (D2Q9 MRT-LBE model in Fig. 1) is selected, and these directions are defined as e_α=0=(0, 0), e_α=1-4=c(cos ((α-1)π/2), sin ((α-1)π/2)), and , sin ((α-1)π/2+π /4)), where c=δx/δt stands for the magnitude of the particle velocity, and δx refers to the dimensionless lattice length. For simplicity, δx and δt are set equal to 1, that is, c=δx=δt=1. In particular, the corresponding transformation matrix M is defined as

(2)

For the D2Q9 MRT-LBE model, the nine components {m_α(r_i, t)|α=0, 1, ···, 8} of the moment vector m= (ρ, e, ε, j_x, q_x, j_y, q_y, p_xx, p_xy)^T are set as the below order: m₀=ρ is the mass density, m₁=e stands for the kinetic energy, m₂=ε is related to the kinetic energy square, m₃=j_x and m₅=j_y are the two components of the momentum, m₄=q_x and m₆=q_y stand for the two components of the energy flux, and m₇=p_xx and m₈=p_xy are the stress tensors. The nine components {m_α ^(eq) (r_i, t)|α=0, 1, ···, 8} of the equilibrium moment vectors m^(eq)=(ρ, e^(eq), ε^(eq), j_x, q_x^(eq), j_y, q_y^(eq), p_xx^(eq), p_xy^(eq))^T consist of two parts, i.e., the conserved moment that includes mass density m₀^(eq)=ρ and two compone-nts of the momentum m₃^(eq)=j_x, m₅^(eq)=j_y, and the non-conserved moments that are shown as

(3a)

(3b)

(3c)

Thus, the equilibrium PDF f_α^(eq) (r_i, t) which corresponds to the equilibrium moments m^(eq) is defined as

(4)

where ω_α is the weight parameters, ω₀=4/9, ω_1-4=1/9, ω_5-8=1/36, and stands for the lattice sound speed. The macroscopic velocity u and the density ρ are functions of f_α (r_i, t), i.e., .

In Eq.(1), the diagonal relaxation matrix S is given by

(5)

where s_i (i=0, 1, ···, 8) stands for the relaxation rate corresponding to nine velocity directions. Since the incompressible NSE could be deduced from the MRT-LBE by using the multi-scale technique Chapman-Enskog expansion, the viscosity υ, which is obtained from this expansion under the condition s₇=s₈=1/τ where τ is the dimensionless relaxation time corresponding to the viscosity, is shown as

(6)

2.2 LES with D2Q9 MRT lattice Boltzmann model

In the LES, any physical variable of the government equations, e.g., f, is decomposed into a resolved component and an unresolved component f' by using a filter function so that the filtered equations for the large eddies (resolved component) could be directly computed, and the effects of the smaller eddies (unresolved component) are replaced by using a selected model to represent it, e.g., the Smagorinsky eddy viscosity, υ_t, model ^{[7, 26]}. With the use of LES in the D2Q9 MRT-LBE, the total viscosity υ is expressed as the sum of the molecular viscosity υ₀ and the eddy viscosity υ_t, i.e., υ=υ₀+υ_t^[7]. Thus, the corresponding total relaxation time τ can also be separated into two components: τ=τ₀+τ_t=3(υ₀+υ_t)+0.5, where τ₀=3υ₀+0.5 and τ_t=3υ_t are the dimensionless relaxation time. In the Smagorinsky model ^[26-27], the eddy viscosity υ_t is defined as

(7)

where C_sm is the Smagorinsky constant, Δ stands for the filtered size, which is set as the lattice length, i.e., Δ=δx=1 because of the use of uniform lattice, S^* is the magnitude of the shear strain rate tensor S_ij, and can be computed directly from the second-order moments of PDF, Q_ij. It has been proved ^[10] that

(8)

where s_xx is the relaxation rate for the second-order moments, p_xx and p_xy, that is, s_xx=s7=s8=1/τ. Thus, the eddy viscosity can be estimated as

(9)

where the details on how to calculate Q^* were given in Ref.^[10].

By combining υ=υ₀+υ_t, τ=τ₀+τ_t, and Eq.(9), the eddy viscosity can be obtained as follows:

(10)

2.3 Boundary condition

For 2D cavity flows, the no-slip boundary condition, i.e., the half-way bounce back boundary condition ^[12], is applied on all the boundaries. Because of the use of this boundary condition, the physical boundary is specified between the fluid nodes and the inner fictional ghost nodes. Thus, the collision step is only applied on the fluid nodes. Because of this setup, the unknown PDFs of the fluid nodes that are immediately inside of the boundary can be estimated as follows ^[12] :

(11)

where f_α is the PDF of e_α=-e_α, and U_lid stands for the velocity on the boundary, e.g., the lid speed orzero for a stationary boundary.

3 Parallel algorithm (GPU-based MRT-LBE-LES model) 3.1 Introduction of CUDA

CUDA^TM is a parallel computing platform and programming model that enables the dramatic increase in computing performance by using the power of a GPU. After the release of the CUDA, the parallel computational technique based on GPU has been successfully used in various fields ^[21-22]. The CUDA is designed to support various programming languages or application programming interfaces, e.g., FORTRAN, C, MATLAB, and OPENCL. Thus, the researchers familiar with the standard programming languages can easily code the CUDA programs for solving complex computational problems. The software environment of CUDA consists of the CUDA Toolkit and the CUDA Software Development Kit (SDK). The available operation systems include Windows XP, Vista, 7, Linux, Unix, and Mac ^[21-22].

As shown in Fig. 2, a CUDA program consists of a host part that runs on the central processing unit (CPU) side in serial and a device part that runs on the GPU side in parallel. The device part code, called kernels, can be invoked by the serial code on the CPU side. The primary tasks of serial codes include (ⅰ) initial clarifications, (ⅱ) allocation and de-allocation memory spaces, and (ⅲ) data transmissions between the host memory and the device memory. When the data are transferred to device memory, many threads (up to thousands) are launched by the kernel functions to exploit the data parallelism. In order for an easy programming, these threads with independent one-dimensional (1D), 2D, or 3D thread index are divided into many equally-shaped thread blocks according to the data sizes, and these thread blocks are organized into 1D or 2D grid of thread blocks, as shown in Fig. 2. Because of GPU hardware limitations, the number of threads in one thread block and the number of thread blocks in a grid should be properly chosen.

3.2 Implementation details

The solving of an MRT-LBE model consists of three major calculation parts: (ⅰ) initialization, (ⅱ) collision, stream, and boundary processing, and (ⅲ) convergence check, see the flowchart in Fig. 3(a). The main task of the first part is to initialize all the variables, e.g., the PDFs and the corresponding moments according to the velocity and Reynolds number. Since this is a one-time process and only spends little time, it is unworthy of the concurrent execution. The second part (i.e., collision, stream, and boundary functions) is repeatedly calculated many times for updating all the PDFs. This part takes about 98% of the total computing time in the CPU-only code, and therefore, this part is rewritten as three kernels executed on the GPU device in parallel. The third part is to check convergence. For laminar flows, the convergence is to check whether the results reach the steady state. For turbulent flows, it is used for determining whether the lattice size or the Smagorinsky constant is enough. The reason is that if the lattice size or the Smagorinsky constant is not right, the program will blow up because the relaxation time is close to 0.5. The number of executions of this part is much smaller than that for the second part, and does not have much advantage to be parallelized. Thus, it is run in serial in CPU, see the flowchart in Fig. 3(b).

A 2D lid driven cavity flow deals the flow field within a square enclosure that is driven by a selected uniform velocity on the top lid from left to right. All the other three walls are fixed. Cavity flows are typically selected to verify a numerical scheme, and thus, it is used here again to demonstrate the efficiency of GPU parallel computing. The cavity enclosure can be fitted by using a uniform lattice with lattice units GX and GY in the x-and y-directions, respectively. In other words, there are (GX+1) and (GY+1) nodes in the x-and y-directions, respectively. Although GX=GY in this special case, different symbols are used to map the lattice to the thread. In order to have efficient parallel computation, two ghost layers of nodes (solid circles in Fig. 4(a)) are added around the real boundary. The inner ghost layer is to incorporate the half-way bounce back boundary condition into the LBM, and the outer layer is used for improving the efficiency of using parallel computing in the GPU kernel with the minimum ‘if statement’. More details will be given later.

In the CPU-only program, the PDFs on the D2Q9 MRT-LBE model are often set as a 3D array, such as f[x][y][k], where x and y stand for a node coordinate in the Cartesian system, i.e., 0 ≤ x ≤ GX and 0 ≤ y ≤ GY, and k indicates the nine different directions of particles moving. For the GPU-based program, the 3D array is rewritten as nine 1D arrays on the host memory, i.e., f₀ [P], f₁ [P], ···, f₈ [P], where P=(GX+1)×y+x represents the node index because of the change of a 2D array into a linear addressed memory stored in a 1D array by the row-major convention. Thus, the transfer from the 2D node array to the 1D array starts from the bottom row and moves upward, as shown in Fig. 4(a). For example, Node (3, 3) will change to Node ((GX+1)×3+3). For avoiding the read-write conflict on the GPU side, two sets of these nine linear 1D arrays, i.e., f_{0_old}[P], f_{1_old}[P], ···, f_{8_old}[P] and f_{0_new}[P], f_{1_new}[P], ···, f_{8_new}[P], are allocated on device memory, one set for reading data, and another for writing data.

After the transfer of data from the host to the device memory, one thread which has the general index T=((GX+1)×y_global+x_global) is mapped to a lattice node (GX+1)×y+x, and these threads are equally divided into many thread blocks in one grid, as shown in Fig. 4(b), where x_global and y_global are the global coordinates of one thread, x_t and y_t stand for the local coordinates of one thread in one thread block, x_b and y_b are the coordinates of one block in the grid, and x_t and y_t stand for the numbers of threads in the x-and y-directions in one block, respectively. Thus, the global coordinates of a thread, x_global and y_global, are defined as x_global=x_b×x_t+x_t and y_global=y_b×y_t+y_t, and each thread block has x_t×y_t threads which is chosen by the node assignment.

In our GPU-based program, three kernels are repeatedly invoked for particles: collision, stream, and boundary condition processing, respectively.

Case 1 Collision kernel

This kernel is responsible for the particle collision. Since the particle collisions happen only on the fluid nodes, the computation is only carried out on fluid nodes. As a consequence, the execution of collision kernel is done simultaneously in GPU instead of using a `for loop' in the CPU-only program. Herein, since each thread is mapped to a fluid node, only one ‘if statement’ is used in the collision kernel for choosing the threads that correspond to the fluid nodes. The eddy viscosity is computed by Eq.(10) before the calculations of the post-collision moments of PDFs in the above process.

Case 2 Stream kernel

This kernel is responsible for particle streaming. The particles on all fluid nodes will propagate to the adjacent nodes after the particle collisions on each fluid node, and the PDFs on these fluid nodes that are next to the boundary nodes are still unknown, e.g., the fluid nodes on the x=2 column, as shown in the left part of Fig. 4(a), where the solid arrows stand for the post-stream PDFs in the arrow directions. For example, only three fluid nodes, i.e., (2, GY-3), (3, GY-3), and (3, GY-2), will give particles with the north, northwest, and west directions to the corner fluid node (2, GY-2), respectively. Thus, on the (2, GY-2) node, there are three solid arrows in the north, northwest, and west directions, respectively. No arrows in the other directions on this node show that the PDFs with the other directions are unknown. It is the same to the hatched box on the left side of Fig. 4(a), e.g., the boundary node (1, GY-1) only receives particles from the fluid node (2, GY-2), and so one solid arrow is on it. The stream process can also be viewed as one set of nine procedures, i.e., the PDFs with the same direction on the all fluid nodes are given to the neighbor nodes, totally eight directions (the PDFs with the e₀-direction are at rest). In this kernel function, the nine procedures can be simultaneously performed by the corresponding threads on GPU because only the fluid nodes give particles to the adjacent nodes. In addition, one thread mapped with the selected fluid node has the thread index T=(GX+1)×(y_global)+x_global. Therefore, the thread indexes of the nearest eight neighbor nodes can be searched by the equation N=(GX+1)×(y_global+δt×e_α)+(x_global+ δt×e_α), where α=1, 2, ···, 8. Note that δt=1 and the two components of discrete particle velocity e_α are either 0 or 1. Thus, δt and e_α are only used to find the thread indexes of neighbor nodes.

Case 3 Boundary kernel

The kernel is responsible for boundary conditions, i.e., after the stream kernel, all the unknown PDFs on the fluid nodes are updated using Eq.(11). In the CPU-only program, this process can be easily performed even without ‘if statement’ by returning the PDFs with e_α, which has been updated on the boundary in the stream step, back to the unknown PDFs with e_α on the fluid nodes. Unfortunately, in the GPU-based program, the implementation of the above procedure needs too many ‘if statements’ for choosing the specified threads that correspond to the lattice nodes. This leads to a significant loss of the parallel efficiency. In order to avoid this, another outer ghost layer is introduced (see Fig. 4(a)). The approach proposed in this study is that the PDFs which have the direction of the dotted and solid arrows on the hatched box of the right wall (the GX-1 column) give back to the adjacent left fluid nodes for updating all the unknown PDFs according to Eq.(11)(see the right part of Fig. 4(a)). For example, on the node (GX-1, 2), the solid arrow with the southeast direction is for returning to the solid arrow with the northwest direction on the node (GX-2, 3), the solid arrow with the east direction is for returning to the solid arrow with the west direction on node (GX-2, 2), and the dotted arrow with the northeast direction is for returning to the dotted arrow with the southwest direction on node (GX-2, 1), where the solid and dotted arrows on the right wall are the PDFs that have and have not been updated in the stream step, respectively. In this process, even though the PDFs with the directions of the dotted arrows are assigned to the opposite PDFs of the neighbor nodes (the dotted line arrows on the GX-2 column), it has no effect on the results due to the fact that no changes happen on the PDFs of the fluid nodes, and the macroscopic variables are computed by the PDFs of the fluid nodes only. This has the advantage that the half-way bounce back boundary condition could be easily rewritten as the GPU kernel, and the parallel efficiency of GPU program can be improved significantly by using a minimum ‘if statement’. The additional lattice nodes only request a small extra amount of memory and a negligible small computing time for data copying.

With all that implemented, the massive parallel computing can be used more smoothly, and the efficiency of the MRT-LBE-LES program is excellent for numerical simulations of the 2D lid driven flow. With the capability of using the same double precision in the GPU and the CPU programs, they have the same accuracy.

4 Results and discussion 4.1 2D lid driven laminar cavity flow

Here, the numerical experiments on 2D cavity flows with Re=100, 400, 2 000, 7 500, 8 400 (i.e., the laminar flow) are performed for testing the GPU-based program of the D2Q9 MRT-LBE-LES model with C_sm=0. The dimensionless uniform top driven velocity in the lattice Boltzmann system, U_lid, is set to be equal to 0.1, and the Reynolds number is defined as Re=U_lidL_ref/υ, where υ is the dimensionless kinetic viscosity of the fluid, and L_ref=GX-3=GY-3 is the dimensionless characteristic length used in the lattice Boltzmann system, where the extra 3 is from the half-way bounce back boundary condition and the extra layer for a better parallel efficiency to process the half-way bounce back boundary condition in the GPU kernel. Therefore, the total lattice units in each line are GX.

A uniform lattice unit of 259 (i.e., GX=GY=259, L_ref=GX-3=256) is selected in the x-and y-directions, respectively. The nine components in the diagonal collision matrix are given by s₀=s₃=s₅=0, s₄=s₆=1.2, s₁=s₄-0.1, s₂=s₁-0.1, and s₇=s₈=1/τ. The steady situation of a laminar flow can be claimed if the difference of normalized kinetic energy, δ, between two consecutive checking times is less than a small value, e.g., 10^-12,

(12)

For laminar flows, the simulation results (except for Re=8 400 because no data can compare with it) at the steady situation (i.e., velocity profiles at the center of cavity, see Fig. 5) agree well with the reference ^{[18, 28-30]}. The velocity profiles change from a smooth curve at low Re to a nearly straight line with sharply bends near the boundaries at high Re. The locations of the primary and secondary vortices for different Reynolds numbers also agree well with the references ^{[18, 28-30]} (see Table 1).

In this study, the maximum Reynolds number of 2D lid driven cavity flows which can be computed with the steady solution is 8 400. To obtain the steady solution for the cavity flow with a higher Reynolds number (i.e., Re >8 400), the number of lattice size is increased to 3 075² for allowing the cavity flow to have steady solutions ^[31]. However, we observe that the solutions are not convergent but oscillating with time. This suggests that the 2D lid driven cavity flow field is unstable when Re is larger than 8 400, and the time-averaged results should be used to judge whether the solution converges to a steady state.

4.2 2D lid driven cavity turbulent flow

As Re> 8 400, the computed instantaneous velocity of the lid driven flow will oscillate despite the selection of large lattice units. This indicates the start of changing flow pattern from laminar to turbulent flows. For a fixed number of lattice unit, e.g., 259² lattice units, the D2Q9 MRT-LBE-LES model with C_sm=0 will continue to produce results, although unstable, until Re> 29 796. At that time, the fluctuation is too large, and the program crashes because the relaxation time τ is too close to 0.5. This suggests that keeping increasing the lattice units can make the relaxation time sufficiently larger than 0.5, and thus, it is possible to use the MRT-LBE model for turbulent flows. However, the higher the Reynolds number (i.e., the larger the lattice number), the more the computational cost, and even with the help from a GPU, the computation cost will be too high to do the task on a personal computer for large (more than 3 300²) lattice sizes. That is the reason to implement the LES for simulating high Reynolds number flows with a reasonably selected coarse lattice.

To evaluate the efficiency and capability of the GPU-based MRT-LBE-LES program in simulating turbulent flows, the 2D lid driven flows with five high Reynolds numbers, i.e., 5×10⁴, 1×10⁵, 5×10⁵, 1×10⁶, and 1×10⁷ are checked. The dimensionless driven velocity at the lid U_lid is still 0.1. Since the instantaneous velocity field is fluctuating with time, it is impossible to find a stable flow field unless the time-averaged flow is used. The selection of duration for taking average may depend on the following factors: the number of lattice units, the Smagorinsky constant, the initial running time before taking time average, and the duration of average. Here, we find that the initial running time of 500T (where T=L_ref/U_lid) and averaged time of 150T are sufficient to warrant a stable time-averaged 2D lid driven flow. A longer initial running time and longer average time make little difference in the time-averaged velocity profile. Thus, the initial running time and duration for average are set as 500T and 150T, respectively, for all cases. This is one of the major differences compared with the early studies ^{[7, 17]}.

In addition, possible values for the Smagorinsky constant, from 0.05 to 0.3, are tested to address its effects on the time-averaged results. Different lattice units are also checked to examine its effect. For checking the computation efficiency in different thread numbers, the lattice nodes are arranged as a multiple of 16, i.e., 128², 256², 512², ···, 3 072² lattice nodes. The selection of the finest lattice units, i.e., 3 072² is limited by our graphic card memory.

The results suggest that the Smagorinsky constant C_sm has great effect on the time-averaged velocity, if the lattice units are the same for the fixed Re (see Fig. 6). For instance, the results with C_sm=0.07 are close to those given by Chai et al. ^[17], but smaller than the results with C_sm=0.14 on the 127², 255², and 511² lattice units. It may be because that as C_sm increases, the eddy viscosity also increases, and thus, the momentum can be further transferred into the cavity. It appears that for the same C_sm, there is a minimum requirement for the lattice units. For example, with C_sm=0.07 and Re=5×10⁴ and 1×10⁵, the time-averaged velocity profiles show that the selected minimum lattice units, i.e., 127², is acceptable, but it may be too coarse to provide comparable results for higher Reynolds numbers. For Re=5×10⁵ and 1×10⁶, the results on the 511² lattice units are better to match with the results of Chai et al. ^[17] when compared with that using 255² lattice units, and the lattice units of 127² is not acceptable at all.

After the Smagorinsky constant and the number of lattice are determined, the model computed streamlines are constructed (see Fig. 7). For Re=5×10⁴ and Re=1×10⁵, it can be seen that a primary vortex dominates the entire flow field with several secondary vortices appearing in the three corners, and the shapes and locations of these secondary vortices will change vastly with time. Although only one primary vortex dominates the entire flow field, the center for the primary core is relatively stable for Re=5×10⁴, but the center starts to rotate for Re=1×10⁵, and these secondary vortices fluctuate more. When Re increases to 5×10⁵, it appears that a secondary vortex may be released from the boundary and moves into the area mainly occupied by the primary vortex. When Re=1×10⁶, two primary vortices rotate around the cavity center, and more and more secondary vortices are produced in corners. Generally, as Re increases, the flow field becomes more turbulent and more oscillatory.

For further evaluating the eddy viscosity term in the D2Q9 MRT-LBE-LES model, the simulations of 2D cavity flows at Re=1×10⁷ are performed with 255², 511², and 1023² lattice units. Because there are no other numerical or/and experimental data at this high Reynolds number, no comparison of velocity profiles can be made. The numerical results show that the time-averaged velocity is also greatly affected by the Smagorinsky constant C_sm for the same lattice size at the fixed Re (see Fig. 8). It is assumed that for Re=1×10⁷, C_sm should be slightly reduced, and thus C_sm=0.065 is selected for the three different lattice sizes. It is also observed that the instantaneous streamline contour plots are very different after the 500T initial running.

As the lattice units increase but no change on C_sm, the primary vortex is split and gradually replaced with many secondary vortices (see Figs. 9 and 10(a)). This may be due to the use of the same C_sm for different lattice numbers. Since C_smΔ is corresponding to the mixing length ^[26], the maximum length is limited to Δ (for the traditional mixing length hypothesis, the length could be much larger, e.g., half of the water depth). For this reason, the Smagorinsky eddy viscosity model only addresses small scale, local eddy viscosity, and cannot describe the variation of eddy viscosity globally. For a fixed Reynolds number, if the lattice number is small, the physical lattice length Δ represents must be large. For this reason, C_smΔ should increase as the lattice number decreases. If the local turbulent eddy viscosity becomes excessively large, it will affect the velocity and produce unrealistic velocity profiles. As a consequence, the selection of lattice size will affect the simulation results for a fixed C_sm. This also explains why the time-averaged velocity increases with C_sm when the lattice size remains unchanged (see Fig. 6).

4.3 Performances of MRT-LBE-LES with GPU

In order to accelerate the 2D MRT-LBE-LES program, a personal computer (AMD Phenom Ⅱ 1100T CPU) with 16 GB memory and one CUDA-enable graphic card (GTX 580) is built for this purpose. The GTX 580 has 512 CUDA cores and 192.4 GB/s memory bandwidth and can deliver over 0.2 TeraFLOPs in double precision operations. The software environment for the computer consists of Windows XP professional 64bit, CUDA Toolkit version 4.1, and GPU driver 286.19 ^[21].

When taking time-average to show the mean velocity field, as the typical approach for turbulent flows, the time-averaged flow pattern is expected (see Fig. 10(b)). In contrast with the instantaneous streamlines plot, these small vortices, except those at the three corners, are smoothed out by taking time-average over a sufficiently long period, e.g., 150T. Moreover, the total kinetic energy k defined as , for this high Re(1×10⁷) is much stronger than that for the cavity flow with Re=1×10⁵ (see Fig. 11), where the notation stands for the time averaged over 150T, and u' and v' are the velocity fluctuations in the x-and y-directions, respectively.

The computational efficiency of lattice Boltzmann simulations is often measured according to how many million lattice nodes could be updated per second (MLUPS) ^[23-25]. Since the calculation of the eddy viscosity has negligible effect on the MLUPS, the 2D lid driven flow at Re=7 500 is taken as the example (see Table 2). In general, the calculation efficiency of GPU-based program increases with the number of threads per block, but Table 2 indicates that the efficiency is about the same when more than 32 threads per block are used. The best performance is achieved with 128 threads per block for different grid sizes. The comparison of MLUPS between GPU-based (128 threads) and CPU-only program on the 512² lattice shows a ratio of 338/4.317, about 78 times improvements.

5 Conclusions

In this study, the implementation of D2Q9 MRT-LBE-LES model based on the CUDA technology is presented, and the numerical codes are validated by the 2D lid driven flows over a large range of Reynolds numbers.

For the 2D lid driven flows at low Reynolds numbers, i.e., laminar flows, the results agree well with the numerical results from the others. The maximum Reynolds number of 2D cavity flows which can be simulated for the steady solution by using the MRT-LBE-LES model is Re=8 400. When the flows are changed from laminar flows to turbulent flows, i.e., Re=5×10⁴, 1×10⁵, 5×10⁵, 1×10⁶, and 1×10⁷, the time-averaged variables should be used because the instantaneous variables are always fluctuating. The modeling results show that MRT-LBE-LES model can well simulate the 2D lid driven flows at high Reynolds numbers, up to 1×10⁷.

As the change of C_sm and lattice number may alter the local eddy viscosity, the entire flow field will be affected. As the increase of C_sm or the decrease of physical lattice size roughly means the increase of local eddy viscosity (if the shear strain rate tensor remains the same), excess momentum transfer from the top lid to the water below increases the mean velocity, and thus, larger and probably unrealistically large velocity may result. This implies that there is a minimum requirement on the lattice number, which depends on the Reynolds number, for successful simulations of turbulent flows. More studies are necessary to clearly address this issue.

Since the double precision floating point is used in the CPU-only and GPU-based program, there is the same computational accuracy. An increase of 76 times in computation speed when compared with CPU-only codes shows that the GPU massively parallel computation method is excellent. This advantage could be further enhanced if the limitation caused by the GPU hardware, e.g., 0.2 TeraFLOPs in double precision FLOPs and 1.5GB memory for the GTX 580 graphic card, is relaxed. Simulations with finer resolutions (e.g., 4 096²) lid driven flows or 3D turbulent flows with Re> 1×10⁴ will become possible. Other alternatives, e.g., use of multiple GPUs, or a better GPU, e.g., Tesla C2050, or a new GPU-architecture, Kepler, which delivers 3 to 4 times better double precision FLOPs performance, when compared with that for a GTX 580, is possible.

Since using a CUDA-enabled GPU and a regular desktop PC is inexpensive, smaller, and less power-demanding than a cluster with multi-CPU, the LBM with CUDA-GPU technique can be more frequently applied for simulating fluid dynamic phenomena.

Acknowledgements The authors would like to thank Prof. Lishi LUO for providing the initial MRT-LBE CPU code. This work is supported by College of William and Mary, Virginia Institute of Marine Science for the study environment.