A dual-lattice hydrodynamic-thermal MRT-LBM model implemented on GPU for DNS calculations of turbulent thermal flows

The lattice Boltzmann method (LBM) is a cellular automata method developed as a more stable and computationally effective version of the Lattice Gas Automata (LGA) method. The LGA method, as described in Frisch et al. (1986), tracks discrete particles through a collection of nodes with static connections. When two particles are in the same node at the same time a collision step is carried out and the design of this collision operator impacts the dynamic behavior of the “gas.” In contrast, the LBM does not track individual particles and instead uses ensemble averages in a control volume as quantities of interest. Because of this, the method is often referred to as a mesoscopic method (Chen and Doolen, 1998; Li et al., 2016), where every control volume holds the velocity information of the particles contained in it, i.e. a discretized version of the local velocity distribution. Recently, the method has gained increased interest due to its favorable parallelizability properties, especially for implementations on GPUs. As demonstrated in Delbosc (2015) the expected reduction in calculation time is approximately two orders of magnitude for a GPU implementation as compared to a CPU implementation. Usage of GPUs thus enables fast DNS and LES simulations of complex fluid phenomena and is especially well suited for problems with complex geometries and evenly distributed local Reynolds numbers, such as porous media flows (Gao et al., 2015).

LBM models capable of accurately resolving turbulent structures have been reported by Suga et al. (2015), Geier et al. (2006), Premnath and Banerjee (2011) and Abdul-Kadhim et al. (2019) among others, and as shown in Delbosc et al. (2014), Delbosc (2015) and Li et al. (2016) hydrodynamic models can be combined with a convective-diffusive LBM model to model thermal fluids. Therefore, a combination of these models is assumed to yield accurate models of resolved turbulent thermal fluid phenomena. In this article, a GPU implementation of an LBM model which is stable for low viscosities and thermal diffusivities is reported in detail. The model is tested against three cases for which validation data is available. Excellent agreement is obtained.

The evolution of a gas is governed by the Boltzmann equation and it is from this equation that the LBM is derived according to:

Here, f is the distribution function which, in general, is a function of space (x_i), time (t), particle velocity (c_j) and Ω is a collision operator. If the velocity distribution is discretized, the set of values in it can be written as a vector f_i. The vector c_j is then a discrete set of velocity vectors called e_ij. Usually the procedure is separated into two steps, a collision step and a streaming step. With a discretized version of the distribution function we have the following set of equations:

Here e_ij is, as mentioned earlier, called the discretized velocity vector which depends on how the spatial discretization of the distribution is done and notice that equation (2) is the collision step and equation (3) is the streaming step. Also notice that f_i is the discretized distribution function after being relaxed while $f˜i$ is the distribution before relaxation. The control volume average density (ρ) and velocity (u_i) can be obtained from the distribution function by the relations (Chen and Doolen, 1998):

These quantities are always preserved during the collision step which follows naturally from continuity and momentum preservation. A natural assumption for the velocity distribution of a gas within an isolated control volume is that it should tend toward the Boltzmann distribution as time goes to infinity. This assumption was first used by Bhatnagar et al. (1954). Because discrete time steps are being taken, a relaxation parameter (τ) is assigned which decides how quickly the control volume approaches the equilibrium state. The collision operator can then be written as:

Here, τ is a relaxation parameter being related to the diffusivity of the gas, i.e. the viscosity for the hydrodynamic equations and thermal diffusivity for the diffusion equation. The term $fieq$ is the equilibrium distribution, the form of which decides the set of equations that are solved. These models are referred to as single relaxation time (SRT) models.

Moments are constructed from the discretized velocity vector and they have a physical interpretation such as momentum or density. In general, the moment $μxnymzo$ can be written as (Premnath and Banerjee, 2009):

where the summation convention is not used for the index j. When the moment operator acts on a distribution f_i the desired physical quantity is calculated. For example, $fieix$ gives the x-direction momentum.

The LBM can be used to solve many different coupled PDEs depending on the choice of equilibrium distribution. Some more exotic examples include the shallow water equations (Tubbs and Tsai, 2019) and economic simulations of income distribution (Cerdá et al., 2013). In this work, the forms of the equilibrium distribution will be chosen so that the Navier–Stokes equations for the hydrodynamic modeling and the advection-diffusion equation for the thermal modeling can be solved.

In the Bhatnagar–Gross–Krook (BGK) model, the vector $fieq$ represents the equilibrium state given by the Boltzmann distribution and because exponentials are computationally expensive to evaluate, Taylor expansion version is used (Qian et al., 1992; Chen and Doolen, 1998):

In this case, the summation convention is applied over the index j but not i. Also, w_i are weights and details regarding how to derive these quantities can be found in Kruger et al. (2017). This version of $fieq$ is only valid when the lattice velocity c is equal to 1, i.e. natural units are used. The BGK-SRT collision model has drawbacks as compared to the multiple relaxation time (MRT) and cascaded types because it has stability issues and is consistently outperformed in all areas (D’Humières et al., 2002). Although the process by which the relaxation is carried out differs between the MRT and SRT models the equilibrium function is usually the same as equation (8).

The thermal equilibrium distribution can be written in the following form (Li et al., 2016; Yoshida and Nagaoka, 2010; Delbosc et al., 2014):

In equation (9), T is the temperature of the lattice point defined analogously to the density ρ, see equation (4) as:

u_j is the local velocity in coordinate direction j, this velocity can either be specified or obtained from a hydrodynamic solver. The variable ϵ depends on the choice of lattice and for the case of a D3Q7 lattice as used in this work the value is ϵ = 4 while for a D3Q6 lattice ϵ = 2.

The BGK-LBM and MRT-LBM collision operator, see equations (2) and (3), can in its general form be written as:

Here, M_ij is the moment matrix, $Skj^$ is a diagonal matrix called the relaxation coefficient matrix and F_i is some general forcing term. For the special case of a BGK model the relaxation coefficient matrix is the identity matrix times a constant. F_i is a general source term which, in general, can be a function of f_i depending on the forcing scheme. Changing the values of the relaxation matrix affects the numerical stability. This general form holds true for any MRT-LBM model but the choice of M_ij, $Ski^$ and $fjeq$ will affect the set of equations being solved and the numerical stability.

Two types of boundary conditions are applied in this work, periodic boundary conditions and bounce-back boundary conditions. The periodic boundary condition connects the domain along a specified axis as illustrated in Figure 1. The application of the boundary condition needs to be made with proper consideration for the structures present in the physical system that is being modeled. This is because the condition acts as a box filter for the large-scale structures and may remove some physically important processes if not set up in a proper way (Forslund et al., 2021). For the validation cases presented, the same domains as the reference cases are used.

The bounce-back boundary condition is a second-order accurate no-slip wall condition (Mohamad, 2011; Bouzidi et al., 2001; Ginzburg and D’Humières, 2003) which inverts all the velocities going into the wall which is then streamed back during the next time-step. The following operation is applied to the wall element at each time-step:

For the implementation, a discretized velocity vector e_ij is chosen in such a way that an inversion of the vector along the element axis i corresponds to the reflection operation specified by equation (12). An effect of the bounce-back boundary condition is that the wall distance for the elements closest to the wall is half an element side-length (Mohamad, 2011).

In this work, the construction of the collision operator Ω(f_i) will be designed in such a way that two different types of partial differential equations are solved for. These are the hydrodynamical equation which describes an incompressible fluid, sometimes called the Navier–Stokes equations:

where u_i is the velocity, p is the modified pressure which is equal to the physical pressure times density, a_i is a source term and ν is the kinematic viscosity. For the thermal part the collision operator is adapted so that the convection-diffusion equation is solved for:

Here, ϕ can be any variable which are convected by the flow, ϕ_s is a source term and α is the diffusion coefficient. For the case of a thermal field, ϕ = T and α is referred to as the thermal diffusion.

In this section, the implementation of the thermal and hydrodynamic model is described in detail.

For a 3D gas lattice automata there are many possible choices of lattice, these include the D3Q7, D3Q13 and D3Q19 among others (Mohamad, 2011; Geier et al., 2006). In this notation, a DnQm lattice would have n spatial dimensions and m number of distributions. As discussed in Geier et al. (2006), popular velocity sets such as those mentioned in the preceding sentence are insufficient to adjust all the relevant moments independently, for example, the μ_xyz moment which is important for the short wavelength flow structures by its impact on the shear rate advection. Because properly resolving turbulent structures are of interest for the hydrodynamic part, a D3Q27 lattice set is used, the basis vectors are defined by:

In these equations $⌊⌋$ is the flooring operator and % is the modulo operator, the subscript index 0, 1 and 2 is the coordinate direction which are also sometimes referred to as e_ix, e_iy and e_iz. Written out explicitly e_ij becomes:

This ordering differs as compared to most works and is useful because the wall bounce-back condition, in this formulation, is equivalent to a reversal of the streaming operator. The D3Q27-distribution is visualized in Figure 2, where red are length 1, blue are length $2$ and green are $3$⁠; Each link corresponds to a different weight, there are four distinct weights these are ω₀ = 8/27 (grey), ω₁ = 2/27 (red), ω₂ = 1/54 (blue) and ω₃ = 1/216 (green). Written out the weight vector w_i becomes:

In this work, the same moment basis as that presented in Premnath and Banerjee (2011) is used. In terms of the vector set the following orthogonal matrix basis is applied:

In this equation, the same index indicates element-wise multiplication instead of the more common summation convention. Although a pair-wise orthogonal moment basis is chosen for this model it may be noted that the moment basis does not need to be orthogonal to obtain good numerical stability as shown by Fei et al. (2019). The choice of relaxation factors in the matrix is key for the stability of the computations. As noted by the solver implemented in Suga et al. (2015) more relaxation is required for the energy terms. Ideally it is preferred that all relaxation terms except for those related to moment, density and viscosity should be exactly 2, as to minimize the impact of the moment relaxation on the solution. Guided by this and a limited parameter sweep, the following values for the relaxation matrix are used:

The relaxation factors for the preserved moments, i.e. physical quantities that does not change dues to conservation laws (density and momentum), are set to zero. Moment 4–8 relates to the viscosity tensor and is changed to adapt the viscosity. The rest of the relaxation factors are set to a value which promotes stability. The relation between the viscosity and relaxation factors for the moments s_ν can be derived via a Chapman–Enskog analysis. This has been done before in, for example, Premnath et al. (2009), and the following relation between the relaxation and viscosity holds true for the D3Q27 lattice:

The forcing scheme used for the body force is first-order and is applied to the postcollision distribution function, hence the forcing term is described by the expression (Suga et al., 2015):

where summation is not carried out across the j index and a_i is the fluid acceleration in coordinate direction j.

The model implemented for the thermal D3Q7-model is described in Yoshida and Nagaoka (2010) and the basis vectors are defined by:

A visualization of the distributions can be seen to the right in Figure 2. Each link corresponds to a different weight, there are two distinct weights which are ω₀ = 1/4 (grey in the middle), ω₁ = 1/8 (red). Written out the weight vector w_i becomes:

The choice of moment basis results in the following moment matrix:

In this equation, $(eij)2$ is the square of the norm of the discretized velocity distribution for link i. The relaxation matrix $Skj^$ for this thermal model can contain additional off-diagonal terms dependent on whether the diffusion tensor is anisotropic or not. Because only fully isotropic media will be considered in this work, the following form of the relaxation matrix is obtained:

Here, s₀ = 0 because the value of it does not affect the numerical solution at all and it can be omitted. The factors s₄, s₅ and s₆ do not affect the leading-order approximation but impact the error terms. For the CUDA implementation, these terms are kept as compile-time constants that can be adjusted for accuracy and stability. For the D3Q7 lattice, a Chapman–Enskog analysis results in a different relation between the thermal diffusivity α and the corresponding relaxation factors s_α related to it according to:

The numerical speed of sound c_s differs between the D3Q27 and the D3Q7 model but as noted by Hajabdollahi and Premnath (2018) the models can still be used in conjunction. The scheme used for the thermal source is first-order and is applied to the postcollision distribution function, the source term described by the expression:

where S is the source term.

When simulating fluid dynamical systems where the density varies with temperature to the extent that the liquid cannot be treated as isothermal, the Boussinesq approximation can be used, given that the density difference is sufficiently small. A model for handling non-isothermal fluids is necessary if natural convection influence the flow. The approximation assumes that the expansion related to compressibility is negligible, therefore the density differences impact on the dynamical behavior can be approximated as a difference in the gravitational acceleration acting on the fluid elements, i.e.:

where a_i is the acceleration vector acting on the fluid element, g_i is the gravitational acceleration, β is a dimensional coupling constant and T is the temperature. For numerical stability, the variation of T should be kept around zero to avoid large pressure gradients in the simulation domain.

There are some considerations that need to be made, regarding memory alignment and access, to enable optimal performance of the code. The computational domain consists of an equirectangular lattice of cubic elements where each element holds the data presented in Table 1.

In Table 1, f_i is the double-buffered D3Q27 and D3Q7 distribution for the hydrodynamic and thermal LBM solvers, respectively, u_i is the velocity vector, ρ is the density, T is the temperature and the wall character is a flag that is used to mark whether the current cell is a wall cell or not. The solver is implemented so that 64-bit floating point accuracy can be used instead of 32-bit, but this leads to significantly longer calculation times on most GPU architectures since they are not optimized for double precision. Because the Mach number is low, it is not important to save the velocities in each time-step, therefore a saving interval is set to some specified number of time-steps to save memory bandwidth in the main loop.

As indicated in Delbosc et al. (2014), Delbosc (2015) memory alignment, when the memory is transferred to the kernels, is important for optimal performance. Therefore, the memory is arranged as indicated in Figure 3. The memory arrangement, called coalesced memory, is necessary for the transfers to be aligned and therefore enabling the usage of the entire memory bus. This, however, introduces a problem because the streaming step is necessarily noncoalesced for this memory arrangement. It was concluded in Delbosc et al. (2014) that reading in values from the surrounding nodes before executing the collision kernel is faster than saving to the surrounding nodes after the collision calculations. This is due to uncoalesced reads being faster than uncoalesced writes. Therefore, the streaming step is carried out at the beginning of the combined collision-streaming kernel.

The overall program structure consists of an initialization stage and three main loops (see Figure 4). The inner loop updates the distributions using the MRT-LBM algorithm. After n − 1 operations a special version of the kernel is called that saves monitored variable data to GPU memory. The data is then sampled and transferred to the host memory. Once the outer sweep is carried out m times, data is saved to the hard drive and the program either terminates or carries out another outer sweep. This loop structure is designed to avoid the progressively slower bottlenecks of host RAM memory and hard drive transfer speeds by minimizing the amount of transfer calls. With a proper choice of the intervals n and m, the impact of saving data on simulation speed is negligible.

The stream and collision kernel is divided into four steps, the streaming step which reads in the surrounding distributions, the fluid collision step which applies equation (11) to the D3Q27 model, the thermal collision step which applies the same equation to the D3Q7 model and finally the saving step which saves all relevant variables to the GPU memory.

Because the simulation memory is stored on the GPU it is possible to use the memory to render the flow fields directly by using graphical libraries such as Open Graphics Library (OpenGL) or DirectX. The solver presented in this article can be coupled to an OpenGL live renderer by defining the memory/variables of interest for rendering as a glBufferData() object. This gives the ability to both CUDA and OpenGL to access the memory of the variables of interest such as density, velocity and temperature. Video examples of live rendering are provided as appended additional material to the article.

The code is open source and available for download from the repository https://gitlab.com/c8383/thermal-mrt-lbm

The solver is validated by comparison to three thermal DNS test cases, the plane channel flow case as presented in Kawamura et al. (1998), a porous media test case as presented in Chu et al. (2019) and a Rayleigh–Bénard convection case compared to the data presented in Wörner (1994). The two-way coupling between the thermal and hydrodynamic lattice is only applied to the Rayleigh-Bénard case.

Plane channel flow is a standard case for validating CFD codes and is also of academic interest due to the flow properties dependence on the nature of near-wall turbulence. Because of this, there are many sources of reliable DNS data for both the turbulent and thermal statistics of this case, some examples are: Alcántara-Ávila et al. (2021) and Kawamura et al. (1998). For the thermal plane channel flow case there are walls at the top and bottom where a constant temperature and no-slip condition are applied. On the other walls periodic boundary conditions are used, see Figure 5. The temperature source is varied so that the average temperature of the domain is equal to 100 (natural units). The flow is driven by a constant pressure gradient F_x(t) which is set to a value so that the time-average Re_τ is constant. The mesh consists of cubic elements of unit size in a lattice of size nX × nY × nZ = 192 × 288 × 576.

There are three Reynolds numbers that are commonly used to characterize the dynamics of the flow for this case, these are:

where Re is based on the mean channel stream-wise velocity $U¯$⁠, the height of the channel L and kinematic viscosity ν. Re₀ uses the channel central velocity U₀ and channel half height δ = L/2. Re_τ is called the friction Reynolds number and is based on the friction velocity which is defined as a function of the wall shear as:

where the last approximation can be made because ρ is kept at approximately 1 across the flow field. The wall shear stress τ_w is defined as:

where z is the wall normal coordinate and z = 0 is at the wall. From these quantities the wall shear distance z⁺ and wall shear velocity u⁺ can be defined as:

The dimensionless number which characterizes the thermal dynamics is the Prandtl number, which is given as:

where α is the thermal diffusivity. The numerical values used for the simulation are:

Natural units are used for all values and are therefore dimensionless. The initial velocity distribution is a completely uniform x-direction velocity profile with added numerical noise. The initial temperature profile is a uniform profile of T = 100. The flow field is allowed to develop for 2,500,000 time-steps which is equivalent to approximately 80 wash-outs before data sampling begins. The data is sampled over the time-span of 2,500,000 time-steps with samples being taken every 40th time-step. Complete field data of standard deviations and mean values are exported every 100,000 time-steps, or approximately 3.5 wash-outs.

The thermal porous media flow is validated against one of the cases presented in Chu et al. (2019). The geometry consists of cylinders with a quadratic cross section placed in a staggered arrangement, see Figure 6. The flow and thermal dynamics in the porous media are characterized by two dimensionless numbers, these are Re and the Prandtl number Pr, Re is here defined as:

where U_Darcy is the stream-wise velocity averaged over both the solid and liquid phases, d_c is the square rod side-length and intrarod distance as illustrated in Figure 6. α is the thermal diffusivity. The flow is driven by a varying force F_x(t) which enforces a constant Darcy velocity, and a source term for the temperature keeps a constant average temperature $T¯$ of 1 across the entire domain. The mesh consists of cubic elements of unit size in a lattice of size nX × nY × nZ = 256 × 256 × 256. The numerical values used for the simulation are:

The initial velocity is zero for all components in the whole domain and the average temperature is set to 1 for all fluid elements. The flow is allowed to initialize for 6,000,000 time-steps, or approximately 460 wash-outs before sampling begins. Values are sampled every 40th time-step. Complete field data of standard deviations and mean values are exported every 400,000 time-steps, or approximately 30 wash-outs.

For Rayleigh–Bénard convection there are two relevant dimensionless numbers (Wörner, 1994; Yilmaz, 2020), these are usually chosen to Pr, see equation (39), and the Rayleigh number Ra defined as:

and using natural units g = 1, ΔT = 1 Ra can be written as:

The geometry and boundary condition of the case is presented in Figure 7 (Wörner, 1994). The flow is driven by the acceleration g and temperature difference between the top wall T_w = 0 and bottom wall T_w = 1. The domain is periodic in the x and y directions with the side length to height ratio being 8. The mesh consists of cubic elements of unit size in a lattice of size nX × nY × nZ = 512 × 512 × 64. The values used for the variables in the simulation are:

These simulation values yield Pr = 0.71 and, as a result, Ra = 630 × 10³. The velocity field is initially set to zero everywhere with numerical noise added, the thermal field is initialized as T = 0.5 for all fluid elements. The flow is allowed to stabilize for 2,000,000 time-steps, or approximately 100 eddy turnover cycles, before sampling of data starts. Values are sampled every 40th time-step. Complete field data of standard deviations and mean values are exported every 200,000 time-steps.

Here the results for the plane channel flow, porous media and the Rayleigh–Bénard convection are presented.

The validation case for the plane channel flow is carried out for Re_τ ≈ 180. A volume render of the instantaneous temperature field from the LBM simulations clarifies that the turbulence influences the temperature field near the wall, see Figure 8. Next the hydrodynamical statistics of the LBM case is validated against the DNS data presented in Kawamura et al. (1998). The hydrodynamical statistics of the LBM case (blue) and Kawamura et al. (dashed red) are presented in Figure 9. The components of the kinetic energy budget are the production (P), the dissipation (ϵ_k), the turbulent convection (T), the viscous diffusion (V) and the pressure transport (p_T) which are marked in the figure. For further details on the kinetic energy budget terms and their significance the reader is referred to Abe et al. (2001). Overall there is excellent agreement between the LBM simulations and Kawamura et al. data. Regarding U⁺, which is defined in equation (38), the main difference is at small y+ values. Here, U⁺ is larger for the LBM simulations because the LBM case near wall resolution is lower and the wall resolution is at a distance of ≈ y + = 1. At the wall U⁺ is, however, zero for both cases, which is, naturally, not captured with the logarithmic scale applied.

Regarding the Reynolds stresses $ui′uj′$ there are only minor differences for $u′2$ and $v′2$ for limited ranges of y⁺. There are also minor differences for the energy budget terms. These smaller discrepancies are judged to be attributed to normal numerical errors. The comparison of the thermal statistics, see Figure 10, is likewise in excellent agreement with the same deviation as for U⁺ at low y⁺ due to the LBM wall resolution.

For the porous medium studied there are complex interactions between the thermal field and flow field structures. This is particularly apparent in, for example, the horseshoe vortices folding across the cylinders, see the volume plot of an instantaneous thermal field in Figure 11. By averaging results from the LBM simulations, flow and thermal data over time can be compared to the pore scale data presented in Chu et al. (2019). This comparison yields that the average trends of the flows are the same, also in this case, indicating excellent agreement between the cases, see Figure 12.

An instantaneous thermal field for the LBM simulations of Rayleigh–Bénard convection can be seen in Figure 13. The temperature field takes the form of a typical Rayleigh–Bénard convection pattern consisting of plumes where the fluid sinks or rises, hence so-called Bénard cells are developed. A comparison to Wörner (1994) of the thermal statistics is presented in Figure 14. The comparison shows excellent agreement between the cases.

For GPU implementations, the goal is to have a maximum utilization of the memory bandwidth (Delbosc, 2015; Delbosc et al., 2014). For this algorithm the theoretical maximum Lattice Updates per Second (LUPS) can be derived as:

where BW is the bandwidth of the graphics card, f_n is the amount of distributions being streamed and T_s is the size in bytes of the type being used for the distribution. For the RTX8000 cards used in this work, the theoretical maximum becomes $672e9/(2×34×4)≈2470MLUPS$⁠. The performance is tested on a benchmark test case consisting of a cubic lattice of elements n³, the case is initialized and then the code is executed for 2,000 time-steps. The longest execution timing, the $n3=3843$ case is approximately 60 s which means that thermal throttling effects are negligible. A summary of the values given from performance measurements can be seen in Table 2. The maximum performance is around 80% of the theoretical max capacity, this is the same utilization degree as presented in Delbosc et al. (2014) and Delbosc (2015). The performance begins to scale linearly when the amount of elements reaches approximately 1 M which is also on par with the results presented in Delbosc et al. (2014) and Delbosc (2015).

An MRT-LBM model with hydrodynamic (D3Q27) and thermal (D3Q7) lattices capable of simulating thermal turbulent flows has been implemented in the CUDA framework. The source code for the solver is openly available at the repository https://gitlab.com/c8383/thermal-mrt-lbm. The model is successfully validated with existing DNS data from three different turbulent thermal flow cases. Furthermore, the performance is on par with earlier published results from the area.

This work was made possible by funding from the Swedish Research Council (Grant No. 2017-04390).

In the interest of transparency, data sharing and reproducibility, the author(s) of this article have made the data underlying their research openly available. It can be accessed by following the link here: https://gitlab.com/c8383/thermal-mrt-lbm