The lattice Boltzmann method for wind farm simulations: a review

High-fidelity numerical simulations are currently a key method to gain insight into flows in very large wind farms. However, these simulations are extremely costly in terms of computational resources, and the progress in computational efficiency has been outrun by the growing size of wind farms and the need for simulations. The adoption of massively parallel hardware, namely graphics processing units (GPUs), by the wind energy community has begun, but the numerical structure of the Navier–Stokes equations hinders the efficient use of such hardware. This is one of the reasons that the lattice Boltzmann method has gained increasing attention in recent years. By construction, this method for simulating the Navier–Stokes equations is perfectly parallelizable and well suited for massively parallel hardware. However, as with every new method, the foundation of the method is not widely known by the wind energy community and often met with doubt. This review paper collects the various methods necessary for a potential GPU-resident wind farm flow solver based on the lattice Boltzmann method. Furthermore, it discusses various aspects of the application of the lattice Boltzmann method to wind farm flows and related flow regimes. It also identifies gaps in the current literature and aims to direct future research on the lattice Boltzmann method for wind farm flows.

1 Introduction

High-fidelity wind farm simulations are an essential tool for both research and industrial applications in the field of wind energy. The continuous increase in wind turbine size, combined with the expanding scale of wind farms, has led to more complex and computationally intensive simulations. In industrial contexts, it has been empirically demonstrated that the ability to complete simulations overnight is a critical benchmark for acceptance (Asmuth et al., 2023; Löhner, 2019). However, the computational demands of large eddy simulation (LES) for wind energy with conventional methods are significantly too high for industrial application and limit scientific progress. To address this challenge, researchers are exploring novel approaches to conducting wind simulations.

One promising approach involves making use of the parallel processing capabilities of graphics processing units (GPUs), which can significantly reduce both simulation time and associated costs. The ExaWind project (Sprague et al., 2020), featuring AMR-Wind and Nalu-Wind, attempts to port traditional approaches to GPUs. The proprietary solver GRASP by Whiffle (Schalkwijk et al., 2012), based on the Dutch Atmospheric Large-Eddy Simulation, is emerging as a commercial application of LESs, and FastEddy (Sauer and Muñoz-Esparza, 2020) uses a numerical discretization specifically tailored to be more suitable for GPUs.

A fundamentally different approach to LES is based on the lattice Boltzmann method (LBM). Rooted in kinetic theory, the LBM offers extensive potential for parallelization, enabling efficient and rapid simulations with significant applicability to wind energy. Due to these advantages, the method has attracted growing attention in the wind energy community. However, its relative novelty has left several aspects insufficiently investigated or not yet fully understood by researchers. As a result, its application in wind energy has remained limited to date, constrained by both its novelty and historical limitations.

The objectives of this review paper are twofold. First, we aim to collect the numerous methods necessary to implement an LES model based on the LBM for wind energy, provide a brief introduction to the underlying theory, and identify gaps in the current literature. Second, we want to demonstrate to the research community that the LBM is a proven, accurate, and efficient method for the simulation of wind farm flows and familiarize the community with some of the basics of the LBM.

We limit the scope of our review to methods relevant to LESs of wind farms in the atmospheric boundary layer (ABL) and mostly to horizontal axis wind turbines (HAWTs). For any solver to be suitable for wind farm LES, we consider the following criteria to be necessary:

• Low Mach number. ABL flows occur at Mach numbers well below the incompressible threshold (Ma<0.1). Therefore, the solver does not need to simulate compressible flows.

• Stability at high Reynolds number. ABL flows are characterized by high Reynolds numbers (Re>10⁷; Stull, 1988), thus requiring low numerical diffusivity and high stability of the numerical method of the solver.

• Advanced subgrid models. The accuracy of LESs for ABL flows depends heavily on the parametrization of subgrid-scale (SGS) turbulent fluxes (Gadde et al., 2021).

• Wall modelling. Given that ABL simulations are inherently wall bounded, turbulence scales near the surface (within the lowest 10 % of the ABL; Kawai and Larsson, 2012) becomes progressively smaller. Since resolving the smallest scales near the wall is computationally prohibitive, wall models are essential for accurately predicting flow behaviour at the closest nodes from the surface.

• Turbine modelling. Fully resolving rotor dynamics is computationally expensive, especially for large wind farms. Actuator models like the actuator line model (ALM) (Sørensen and Shen, 2002), actuator disc model (ADM) (Sørensen and Myken, 1992), and the more recent actuator sector model (ASM) (Mohammadi et al., 2024) offer computationally efficient alternatives by modelling the effect of the turbine via body forces. Hence, the solver must have the capability to integrate body forces and must be stable for large jumps in force over relatively small volumes.

• Complex terrain and forested areas. Onshore wind farms are often situated in remote, forested, or hilly regions, where terrain complexity significantly impacts wind conditions. Effective simulation requires incorporating these effects, as discussed in Elgendi et al. (2023) for complex terrains and Arnqvist (2015) for forest environments.

• Thermal stratification. Atmospheric stability significantly influences critical flow parameters such as wind shear and turbulence intensity. The extensive literature on this topic (Sathe et al., 2013; van den Berg, 2008; St. Martin et al., 2016; Wharton and Lundquist, 2012) underscores the importance of incorporating thermal stratification for accurate simulations.

The remainder of this paper is structured as follows: Section 2 presents the theoretical background and is divided into three subsections: Sect. 2.1 introduces the isothermal lattice Boltzmann method and reviews methods purely related to LBM, Sect. 2.2 discusses specific aspects of wind turbine modelling in the LBM framework, and Sect. 2.3 outlines approaches for incorporating thermal stratification. In Sect. 3 we review results from various LBM applications and address several common questions regarding the method. Finally, Conclusions are presented in Sect. 4.

2 The theory behind the lattice Boltzmann method

The following section gives an overview of the theory behind the LBM and specific aspects related to wind energy.

2.1 The isothermal lattice Boltzmann method

In contrast to conventional numerical methods used in wind energy, which solve the discretized Navier–Stokes equations (NSEs) directly, the LBM is rooted in kinetic theory. While continuum mechanics, from which the NSEs are derived, describes macroscopic scales and individual particles that are considered at the microscopic scale, kinetic theory addresses the intermediate, commonly referred to as the mesoscale. Through Chapman–Enskog analysis, the NSE can be derived from the Boltzmann equation in the limit of dense gases; the Boltzmann equation is therefore the more general equation. However, we want to stress that, although the LBM can be derived from the Boltzmann equation, the LBM solves the NSE, not the general Boltzmann equation. Historically, the LBM evolved from lattice gas cellular automata (LGCA, McNamara and Zanetti, 1988). From the LGCA, the LBM has inherited its excellent parallelizability, encapsulated by the principle “all nonlocality is linear, all nonlinearity is local” (Geier et al., 2015a). This subsection will give a short introduction to kinetic theory, explain the derivation of the LBM, and describe the different models necessary and suitable for a wind farm flow solver.

2.1.1 The particle distribution function

The Boltzmann equation describes the evolution in time and space of the particle distribution function (PDF) $f(\mathit{x},t,\mathit{\xi })$. The PDF describes the likelihood of measuring a particle with microscopic velocity $\mathit{\xi }=[\mathit{\xi },\mathit{\upsilon },\mathit{\zeta }{]}^T$ at location x and time t. The ensuing paragraph follows the description by Krüger et al. (2017). The moments M and central moments C of the PDF are defined as

$$\begin{array}{}\text{(1)}& {\displaystyle }{M}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}=\int {\mathit{\xi }}^{\mathit{\alpha }}{\mathit{\upsilon }}^{\mathit{\beta }}{\mathit{\zeta }}^{\mathit{\gamma }}f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi }\text{(2)}& {\displaystyle }{C}_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}=\int (\mathit{\xi }-u{)}^{\mathit{\alpha }}(\mathit{\upsilon }-v{)}^{\mathit{\beta }}(\mathit{\zeta }-w{)}^{\mathit{\gamma }}f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi },\end{array}$$

where $\mathit{u}=[u,v,w{]}^T$ is the macroscopic velocity. Note that this notation follows Geier et al. (2015b), where the indices denote the order of the moment in that velocity component. The order of the moment is then the sum of the indices, i.e. M₁₀₀ is a first-order moment, C₁₁₀ is a second-order central moment, and so on. The moments can be related to a variety of macroscopic quantities. The density ρ and the macroscopic velocity can be recovered from the PDF by taking the zeroth and first-order moments, respectively:

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\rho }& {\displaystyle }=\int f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi }=M_{\mathrm{000}}\text{(4)}& {\displaystyle }\mathit{\rho }\mathit{u}& {\displaystyle }=\int \mathit{\xi }f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi }={\left[M_{\mathrm{100}},M_{\mathrm{010}},M_{\mathrm{001}}\right]}^T.\end{array}$$

The pressure p has to be related to the density by an equation of state, typically the ideal gas law

$$\begin{array}{}\text{(5)}& p=\mathit{\rho }RT,\end{array}$$

where R is the specific gas constant of the fluid and T the fluid's temperature. It is possible to extend the LBM to non-ideal fluids, for example, using the Shan–Chen or the free energy models (Shan and Chen, 1993; Swift et al., 1995), or the approach discussed in Sect. 4.3.3 of Krüger et al. (2017). However, we will not discuss these models in this review due to their limited applicability in atmospheric boundary layer flows. Second-order moments also have macroscopic meaning. For example, the total energy E and the inner energy e can be recovered via the trace of the second-order moments and central moments, respectively:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\rho }E& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(M_{\mathrm{200}}+M_{\mathrm{020}}+M_{\mathrm{002}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\int {\left|\mathit{\xi }\right|}^{\mathrm{2}}f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi },\text{(7)}& {\displaystyle }\mathit{\rho }e& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(C_{\mathrm{200}}+C_{\mathrm{020}}+C_{\mathrm{002}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\int {\left|\mathit{\xi }-\mathit{u}\right|}^{\mathrm{2}}f{\mathrm{d}}^{\mathrm{3}}\mathit{\xi },\end{array}$$

and the central moments of second order and the stress tensor σ are related by

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\sigma }}_{x^{\mathit{\alpha }}{y}^{\mathit{\beta }}{z}^{\mathit{\gamma }}}=-C_{\mathit{\alpha }\mathit{\beta }\mathit{\gamma }}\forall \mathit{\alpha }+\mathit{\beta }+\mathit{\gamma }=\mathrm{2}.\end{array}$$

The molecules in a gas tend towards an equilibrium state, described by the Maxwell distribution that only depends on temperature, velocity, and density:

$$\begin{array}{}\text{(9)}& {\displaystyle }{f}^{\mathrm{eq}}(\mathit{\xi },\mathit{u},\mathit{\rho })=\mathit{\rho }{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }RT}}\right)}^{\mathrm{3}/\mathrm{2}}\exp \left({\displaystyle \frac{-{\left|\mathit{\xi }-\mathit{u}\right|}^{\mathrm{2}}}{\mathrm{2}RT}}\right).\end{array}$$

It can be shown that the equilibrium distribution corresponds to a fluid without viscous stresses. Therefore, the viscous stresses must be contained in the non-equilibrium part of the PDF $f^{\mathrm{neq}}=f-f^{\mathrm{eq}}$. A common way to analyse the Boltzmann equation is the Chapman–Enskog analysis. It is based on a perturbation expansion of the PDF around f^eq:

$$\begin{array}{}\text{(10)}& f=f^{\mathrm{eq}}+\mathit{ϵ}{f}^{\left(\mathrm{1}\right)}+{\mathit{ϵ}}^{\mathrm{2}}{f}^{\left(\mathrm{2}\right)}+\mathrm{\dots },\end{array}$$

with the parameter ϵ being on the order of the Knudsen number, the ratio of the length of mean free path of the particles to a macroscopic length scale.

2.1.2 The lattice Boltzmann equation

The Boltzmann equation is the transport equation of the PDF. Taking the total derivative of the PDF f yields

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{D}f}{\mathrm{D}t}}& {\displaystyle }={\displaystyle \frac{\partial f}{\partial t}}+{\displaystyle \frac{\partial f}{\partial x_i}}{\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}+{\displaystyle \frac{\partial f}{\partial {\mathit{\xi }}_i}}{\displaystyle \frac{\mathrm{d}{\mathit{\xi }}_i}{\mathrm{d}t}}\text{(12)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial f}{\partial t}}+{\displaystyle \frac{\partial f}{\partial x_i}}{\mathit{\xi }}_i+{\displaystyle \frac{\partial f}{\partial {\mathit{\xi }}_i}}{\displaystyle \frac{F_i}{\mathit{\rho }}}=\mathrm{\Omega },\end{array}$$

where Ω is the collision operator. In the Boltzmann equation, this collision operator has to account for all possible collisions between particles and is therefore rather cumbersome. Since the actual operator is of no interest for the LBM, it is omitted here, and approximations of it will be outlined in Sect. 2.1.3. On the left-hand side, the expanded form consists of three terms. The first term is a time derivative, the second term is a convective term, and the third term describes an acceleration. Accelerations and forces are related by Newton's second law, thus forces can be directly represented in the Boltzmann equation. In comparison to the Navier–Stokes equations, no diffusion or non-linearity is found on the left-hand side. All of this behaviour is contained in the collision operator. This separation of convection and non-linearity is the source of many of the advantages the LBM has over traditional computational fluid dynamics (CFD) approaches. It can be shown via the aforementioned Chapman–Enskog analysis that the Boltzmann equation recovers the Euler equations with f=f^eq and with a first-order approximation of the Navier–Stokes equations – see Krüger et al. (2017, p. 26) for further details.

The lattice Boltzmann equation can be obtained from discretizing the Boltzmann equation. For details of the typical discretization on a Cartesian grid, the reader is referred to Krüger et al. (2017, p. 71–98); for brevity's sake only some key results from the discretization will be discussed here. Other discretizations exist, such as for curvilinear grids, presented in Reyes Barraza and Deiterding (2020), and super-sonic speeds, for example, Saadat et al. (2020). However, we have deemed them not relevant here, as we are interested in atmospheric boundary flows. Typically, the lattice Boltzmann method is understood to be discretized on a cubic grid of nodes, with spacing Δx and a constant time step of Δt. The cubic grid is related to the fact that the Boltzmann equation has to be discretized not only in space and time but also velocity space. For wall-bounded flows, such as the ABL, 27 discrete velocities are recommended by Kang and Hassan (2013). This is called a D3Q27 velocity set. Another common choice is the D3Q19 set. A detailed explanation of the construction of velocity sets can be found in Shan et al. (2006). From discretizing the PDF at the discrete velocities ${\mathit{c}}_{ijk}=[i,j,k{]}^T\mathrm{\Delta }x/\mathrm{\Delta }t$, we obtain the populations $f_{ijk}\left(\mathit{x}\right)=f(\mathit{x},{\mathit{c}}_{ijk})$. The D3Q27 velocity set consists of all possible combinations of $i,j,k\in -\mathrm{1},\mathrm{0},\mathrm{1}$. Applying the method of characteristics to the Boltzmann equation, we can obtain the discretized Boltzmann equation:

$$\begin{array}{}\text{(13)}& f_{ijk}(\mathit{x}+{\mathit{c}}_{ijk}\mathrm{\Delta }t,t+\mathrm{\Delta }t)=f_{ijk}(\mathit{x},t)+\mathrm{\Delta }t{\mathrm{\Omega }}_{ijk}(\mathit{x},t).\end{array}$$

Again, the equation is rather simple in comparison to discretized Navier–Stokes equations. It is important to note that Eq. (13) is often split into two operations: a collision step, where the collision operator is computed, and then a propagation step, where the post-collision populations are distributed to the neighbouring nodes. This is equivalent to a split into a potentially non-linear and a non-local step. We want to emphasize that, after a small change of variables, this approximation is already second-order accurate, as can by shown by integration along characteristics or Strang-splitting (Dellar, 2013).

There are a few more aspects of the discretization we need to discuss. First, to satisfy the isotropy requirements necessary for recovering NSE (Krüger et al., 2017), the lattice speed of sound is commonly fixed to $c_{\mathrm{s}}=\frac{\mathrm{1}}{\sqrt{\mathrm{3}}}\frac{\mathrm{\Delta }x}{\mathrm{\Delta }t}$. While this value is standard in LBM, it is not unique: for specific applications (possibly involving modified equilibrium distributions), it can be adjusted but doing so necessitates a full redesign of the discrete velocity set and its corresponding weights. As such cases lie beyond the scope of this review, they are not further considered here. Since the standard LBM solves the weakly compressible Navier–Stokes equations, typically restricted to flows with Mach number $\mathit{Ma}=\frac{V_{\mathrm{0}}}{c_{\mathrm{s}}}<\mathrm{0.1}$, where V₀ is a reference velocity, the time step is restricted to

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }t<{\displaystyle \frac{\mathrm{\Delta }x}{\sqrt{\mathrm{3}}{V}_{\mathrm{0}}}}\mathit{Ma}.\end{array}$$

It should be noted that it is common practice in the LBM to set Ma to an arbitrary value, namely as high as possible while limiting compressibility effects, since it increases time step size and computational speed. Second, we have left out the force term here. Different methods of force-integration exist, which are discussed in Krüger et al. (2017, p. 231ff), but the most common approach is developed by Guo et al. (2002). The final term of Eq. (13) that will need to be discussed is the collision operator.

2.1.3 Collision operators

As mentioned previously, the Boltzmann collision operator is not of interest for the LBM; rather, the collision operator has to be chosen so that the LBM solves the intended equations, in our case the Navier–Stokes equations.

The Bhatnagar–Gross–Krook (BGK) operator. The first collision operator proposed was the Bhatnagar–Gross–Krook operator, introduced in Bhatnagar et al. (1954) to analyse the Boltzmann equation. It takes the following form:

$$\begin{array}{}\text{(15)}& {\mathrm{\Omega }}_{ijk}^{\mathrm{BGK}}={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left(f_{ijk}^{\mathrm{eq}}-f_{ijk}\right).\end{array}$$

It can be viewed as a relaxation of the populations towards the equilibrium with the relaxation time τ. The relaxation time is related to the viscosity ν by

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathit{\nu }=c_{\mathrm{s}}^{\mathrm{2}}\left(\mathit{\tau }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t\right).\end{array}$$

While it is numerically very efficient, it is only suitable for low Reynolds and Mach numbers. With the BGK operator as a starting point, three different approaches to improve its stability and accuracy at higher Reynolds numbers have been developed. One approach is to apply different types of regularizations, while only using one relaxation time. The other approach consists of methods that increase the number of tunable parameters and transform the distributions to moment or cumulant space. For a timeline of the development of these approaches, and a theoretical and numerical comparison, the reader is referred to Coreixas et al. (2019, 2020). The third branch of collision operators, the so-called entropic methods, extend the applicability of the LBM by using different velocity sets and enforcing the H-Theorem from kinetic theory, defined in Boltzmann (1872). This approach can even be extended to simulate supersonic thermohydrodynamic flows – see, for example, Frapolli et al. (2016). However, in the flow regime of interest here, the two former methods appear much simpler and faster (Coreixas et al., 2019). We will therefore limit our review to regularized and multi-relaxation time LBMs. We also omit other ways of solving the LBM, such as the simplified and highly stable LBM by Chen et al. (2017).

The regularized LBM. The possibility of improving the accuracy of the first-order approximation was first shown by Latt and Chopard (2006). This method is based on the reconstruction of non-equilibrium parts of the populations through the discrete analogue of Eq. (8). In the literature, it is often referred to as projection regularization as it can be interpreted as a projection of the non-equilibrium parts of the populations onto the Hermite basis vectors. The projection regularization approach was extended to a recursive regularization to increase stability at high Reynolds and Mach numbers in Malaspinas (2015). The recursive regularization is accomplished through a recursive relationship between the moments of the equilibrium population and the non-equilibrium parts of the PDF. An extension to include thermal and high-order LBMs was presented in Coreixas et al. (2017). By adding a finite-difference reconstruction of the non-equilibrium populations, the regularized approach was extended for explicit and implicit LES in the hybrid recursive regularized (HRR) LBM by Jacob et al. (2018). It introduces a tunable parameter for hyperviscosities, which can be set to add dissipation, similar to a classical explicit subgrid-scale model. A stability analysis of the athermal models can be found in Wissocq et al. (2020). It shows the drastic increase in stability obtainable by using the regularization approaches without introducing any new parameters. One disadvantage of this method is the need to compute additional finite differences. Thus, it relies on non-local quantities which goes against the locality paradigm of the LBM and is clearly disadvantageous in terms of computational performance.

Multi-relaxation times (MRT) LBM. The other large and actively developed group of collision operators originate from the idea of transforming the populations to some moment-like space and applying different relaxation times for these moments. It is possible to choose different relaxation rates because only the relaxation rate of the second-order moments is related to a physical property, namely the viscosity. The first such method proposed is the multi-relaxation time method presented in d'Humières (1994). A careful analysis of the 2D case in Lallemand and Luo (2003) showed how the different relaxation rates are connected and gave some insight into the choice of the higher-order relaxation rates. Following this work, relaxation in different moment spaces were proposed, such as the cascaded LBM by Geier et al. (2006) and the central moment LBM in Dubois et al. (2015). The approach of the cascaded LBM was further developed to the cumulant LBM, introduced in Geier et al. (2015b). It exhibits improved Galilean invariance, although the mechanism used is not limited to the cumulant LBM as the authors stress in Geier et al. (2020), is stable for very high Reynolds number flows, and can be modified for fourth-order accuracy in diffusion (Geier et al., 2017a, b). By adding a finite-difference term, even fourth-order accuracy in convection can be accomplished (Geier and Pasquali, 2018). Furthermore, it is suitable for underresolved simulations, as shown in Geier et al. (2020). The method consists of transforming the populations into cumulant space. Since cumulants are statistically independent, the authors argue this is the only moment space where independent relaxation rates are valid. In its original form, it only relies on the local quantities; however, the transformation into cumulant space is cumbersome.

2.1.4 Boundary conditions

One of the major difficulties in the transition from Navier–Stokes solvers to the LBM is the question of boundary conditions. In the LBM, a boundary condition describes a method of determining unknown populations at the domain boundary. Boundary conditions can be divided into two groups: wet-node approaches and link-wise approaches. An illustration of both methods can be seen in Fig. 1. Wet-node approaches assume that the boundary lies on the boundary node and the macroscopic quantities are used to reconstruct all populations at the node. In link-wise approaches, the boundary is considered to lie between a solid and a fluid node.

Figure 1 Schematic illustration of the two methods for imposing boundary conditions. The left side shows the wet-node approach, with unknown populations indicated by dashed arrows. The right side presents the link-wise approach, using Miller indices to denote the opposite y direction.

Subsequently, only the populations that would stream from the solid domain into the fluid domain have to be determined. Depending on the geometry of the boundary and the conditions imposed on it (velocity gradients, velocities, and/or density), one has to solve an underdetermined or overdetermined system. A plethora of closures for these systems of equations can be found in the literature.

It is easy to see that periodic boundary conditions are analogous to Navier–Stokes solvers and are fairly simple to implement.

Another rather straightforward boundary condition is that of a solid wall. Several approaches exist, but the most common is the bounce-back (BB) approach. As the name suggests, the populations bounce back from the wall to the node they originated from. It can be shown that this method is second-order accurate if the wall is located half the cell width away from the original node. Different interpolation methods exist for walls located at arbitrary distances from the original node. By adding a source of momentum, this approach can be extended to moving walls. Further details can be found in Krüger et al. (2017, p. 175–189).

Full slip boundaries can be implemented similarly to solid walls, but instead of bouncing back, the populations are bounced forward, as described in Krüger et al. (2017, p. 206–208). However, this process ensures symmetry only for straight boundaries.

The biggest difficulty for LBM solvers lies in the open boundaries. Neither the prescription of velocities nor pressures is straightforward, and this is still an active area of research. By prescribing velocity and density, it is possible to find the equilibrium populations on the boundaries. However, since the stresses at equilibrium are zero, this creates an inconsistency if the inflow is not uniform. Another way of prescribing the velocity on the boundary is by using the bounce-back method as described above and by prescribing a wall velocity normal to the wall. See Latt et al. (2008) for a detailed discussion of velocity boundary conditions and Krüger et al. (2017, p. 200) for a discussion of all bounce-back boundary conditions. In classical Navier–Stokes solvers, the outlet condition typically describes a pressure, used as a reference pressure in the simulation. In the weakly compressible LBM, there is no need to prescribe a reference pressure, since the pressure is calculated from the density. So instead, a reference density is needed, which is 1 in the non-dimensionalized form. Thus, the outlet condition has to fulfil two criteria: advection of momentum out of the domain while not exciting spurious oscillations.

The anti-bounce-back approach is similar to the bounce-back approach for solid walls, with the difference that the bounced populations have a flipped sign. To compute the bounced populations, a wall-normal velocity has to be prescribed, which has to be extrapolated from inside the domain. The approach is described in Krüger et al. (2017, p. 200). Another approach is proposed in Appendix F of Geier et al. (2015b), which relies on copying missing populations on the face of the domain from nodes inside the domain. The reflection of pressure waves, which are rather common in the LBM, is a common issue at open boundaries (see, e.g., Krüger et al., 2017, p. 519). Various damping methods have been found that filter these waves – see, for example, Xu and Sagaut (2013) or Appendix F of Geier et al. (2015b). Another option is to introduce a sponge layer: a zone where artificial damping is applied to push the flow towards a steady-state solution. By dissipating kinetic energy in this region, the sponge layer helps to avoid non-physical reflections at the downstream outflow boundary. This technique is also being used in Navier–Stokes solvers, as discussed in Colonius (2004).

2.1.5 Refinement

Within one domain there often exist areas with different requirements in cell size. This is usually the case near walls and, in the case of wind farm simulations, the area of the turbines and the wakes. For traditional finite-difference or finite-volume methods, it is not a problem to change the cell size continuously, within certain bounds. This is not the case with the LBM. Since the ratio of cell width to time step determines the lattice speed of sound, changing the cell width leads to either a change in the speed of sound or the time step. A change in speed of sound is undesirable since it leads to refraction of acoustic waves, for example, as discussed in Schönherr et al. (2011). On the other hand, changing the time step is very disadvantageous for parallelization and would require interpolations of populations. To avoid these problems, the domain can be partitioned into refinement zones. The location of the refined nodes relative to the coarse nodes differs between authors. Either some of the fine nodes coincide with the coarse nodes or all of the fine nodes are placed between coarse nodes. These approaches are called vertex centred and cell centred, respectively. The concept of refinement was first introduced into the LBM by Filippova and Hänel (1998a). The authors use coinciding nodes and compute the populations for the finer grid by using a second-order interpolation of populations of the coarse nodes. While this approach offers flexibility in the refinement ratio, it is not accurate or stable enough for high Reynolds flows. Based on the cascaded LBM, Geier et al. (2009) describe the so-called bubble functions, which make use of the information about derivatives of the velocity that is stored locally in the central moments. These bubble functions are second-order accurate interpolations of the momentum. Thus, in theory, an arbitrarily spaced and oriented fine grid could be placed within the coarser grid. However, the authors limit their examination to a ratio of 2. The bubble functions are adopted to the cumulant LBM and extended to 3D in Kutscher et al. (2019). A different approach, which is based on the hybrid recursive regularized LBM, is presented in Feng et al. (2020). Due to the enforcement of a global entropy equation, correction terms are calculated during the streaming step.

2.1.6 Large eddy simulation

Resolving all flow scales of the ABL is prohibitively expensive. Therefore, simulations of the ABL are limited to Reynolds-averaged Navier–Stokes (RANS) or LES. The LBM can be used for both types of simulation, but LES is more common and of higher relevance for ABL flows, so we will limit our discussion to LES. There exist multiple ways of accounting for the contributions of the subgrid-scale turbulence. The path towards an LES formulation specific to the LBM is described in Sagaut (2010). However, a simpler approach is to employ an eddy viscosity model and to recompute the relaxation frequency with the effective viscosity,

$$\begin{array}{}\text{(17)}& {\mathit{\tau }}_{\mathrm{eff}}={\displaystyle \frac{\mathit{\nu }+{\mathit{\nu }}_{\mathrm{sgs}}}{c^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}.\end{array}$$

Malaspinas and Sagaut (2012) prove that, in the case of the weakly compressible LBM, the filtered Navier–Stokes equations are recovered. Therefore, all the models already developed for Navier–Stokes can be reused in LBM-LES. Note, however, that this is only true for the weakly compressible case. For higher-order discretizations, the effective relaxation time is much more involved and cannot easily be expressed in viscosities. A variety of SGS models have been used in the LBM, such as the classic Smagorinsky model and its refined versions, the wall-adapted local eddy viscosity model (WALE) and the shear-improved Smagorinsky (SISM) (Smagorinsky, 1963; Nicoud and Ducros, 1999; Lévêque et al., 2007); newer models such as the Vreman model (Vreman, 2004), the coherent structure model (Kobayashi, 2005), the QR model (Verstappen, 2011), and the anisotropic minimum dissipation (AMD) model (Rozema et al., 2015) have also been applied occasionally. Turbulence models that only rely on the entries of the strain-rate tensor can be implemented particularly efficiently, as the strain-rate tensor is available locally from the second-order central moments (see Eq. 8). This is the case for the Smagorinsky and QR models. More complex subgrid-scale models such as the Lagrangian-averaged scale dependent model (Bou-Zeid et al., 2005) or models requiring the solution of additional transport equations for turbulence kinetic energy, such as the Deardorff model (Deardorff, 1980) have not been applied, likely due to their high computational cost and additional memory requirements.

Many higher-order collision operators introduce tunable parameters, which can be used for implicit LESs. For example, in Geier et al. (2020), the authors show that the cumulant LBM intrinsically accounts for the high wave number contributions. Gehrke and Rung (2022b) introduce a model that dynamically modulates a Smagorinsky-style eddy viscosity based on third-order cumulants.

2.1.7 Wall modelling

Due to the very high Reynolds numbers in ABLs, the flow close to the ground cannot be sufficiently resolved for LESs. Especially between the first node and the true location of the wall, it is necessary to model the influence of the wall. For this purpose, a number of different relations for the velocity near the wall, called wall functions, have been found, such as Reichardt's, Spalding's, and Musker's laws (Reichardt, 1951; Spalding, 1961; Musker, 1979) or Monin–Obkhov similarity theory (Stull, 1988). When approaching the wall, eddies become smaller, thus making it difficult for subgrid-scale models to correctly predict turbulent viscosity near the wall. Typically, turbulence models overpredict eddy viscosity near the wall. While wall modelling already poses a problem in classical solvers, those difficulties are exacerbated in LBM solvers, due to the difficulties of prescribing shear stresses directly. Numerous wall methods have been proposed in recent years, based on a variety of approaches. In this study, we classify these methods into four distinct categories, analysing the advantages and limitations of each. We present a schematic representation of all approaches in Fig. 2.

Direct population reconstruction. This approach relies on the direct reconstruction of the unknown population at the first wall-adjacent node. This principle was employed for the first implementation of a wall model in LBM-LES by Malaspinas and Sagaut (2014). The authors propose reconstructing the populations at the node closest to the wall as the sum of equilibrium and non-equilibrium parts. The equilibrium is computed from density and velocity at the node obtained with a wall function, and non-equilibrium parts from estimating the shear stress and employing Eq. (16). One disadvantage of the model is that it requires an estimation of the turbulent viscosity near the wall for the reconstruction of the non-equilibrium parts of the populations. Since this model is based on a wet-node boundary approach, it was initially limited to straight boundaries. The approach was extended in Wilhelm et al. (2018) to curved boundaries. First, macroscopic quantities are interpolated to an arbitrary reference point normal to the wall, then the velocity at the node is computed via the wall model, and derivatives are computed from one-sided finite differences. In a final step, the equilibrium and non-equilibrium parts of the populations are computed. Note that this model was originally developed in the context of solving RANS equations with the LBM; however, the method is equally applicable to LES. Further refinements in the computation of macroscopic quantities at the wall are presented in Degrigny et al. (2021).

Another approach to model curved boundaries via reconstruction was proposed by Haussmann et al. (2020), where an interpolated bounce-back step is applied first, followed by a velocity correction at the boundary nodes. This correction is performed by setting the equilibrium part according to the velocity, computed with the wall function while keeping the non-equilibrium parts unchanged.

Wall function bounce (WFB). WFB, introduced by Han et al. (2021a, b), is a variation of the classical bounce forward method, equivalent to a free slip boundary condition. This approach applies a drag term F_τ to the diagonally reflected populations, effectively decelerating the flow in the vicinity of the wall. This drag term is computed from the wall shear stress given by a wall function. Since the wall shear stress is imposed directly, no turbulent viscosity is needed. However, it is unclear whether the model can be extended to a D3Q27 lattice or to curved boundaries.

Immersed virtual wall (IVW). Kuwata and Suga (2021) introduce yet another approach for wall modelling, named the specular reflection bounce-back method. It is very similar to the WFB approach; however, the drag force is applied to all populations, not just those coming from the wall. Because the specular reflection is difficult to extend to curved boundaries, the authors propose a virtual wall layer approach, where an arbitrary number of virtual fluid layers is introduced below the wall. A drag force F_wm is imposed at the node neighbouring the wall to replace the viscous force F_ν with the force due to the wall shear stress F_τ:

$$\begin{array}{}\text{(18)}& F_{\mathrm{wm}}=F_{\mathit{\tau }}-F_{\mathit{\nu }}.\end{array}$$

A bounce-back scheme is applied at the lowest level of virtual fluid nodes. Thus, the model can easily be applied at curved boundaries and any lattice, and does not require the computation of a turbulent viscosity. However, additional nodes need to be included in the domain. Furthermore, by applying the force to all populations, the shear near the wall is not correctly calculated.

Wall velocity approach. The final class of approaches considered relies on determining a fictitious wall velocity u_w, which is then used to modify the BB scheme. Pasquali et al. (2020) leverage the additional information about the stress tensor provided by the cumulants of the populations to compute a skin friction coefficient and from that u_w. This approach is local to the first node off the wall and is applicable to curved boundaries. However, it does rely on a relation of skin friction coefficient and bulk velocity. As noted in the study, it should be applied exclusively with the cumulant collision operator, as less accurate methods fail to provide the necessary precision for an accurate stress tensor approximation. A viscosity is needed to relate cumulants and stresses.

The approach introduced by Asmuth et al. (2021), called the inverse momentum exchange method (iMEM), computes a wall velocity so that the momentum exchanged between the wall and the node closest to the wall matches the wall shear stress computed from a wall function. This method is independent of the collision model, and no turbulent viscosity is required. However, Asmuth et al. present only a formulation for straight boundaries; an extension to curved boundaries is not yet available, although the authors note that it is possible.

Figure 2 Schematic representation of the different wall-modelling approaches in LES-LBM. Symbols in red indicate the required outputs for model implementation, while black denotes known information. (a) Population reconstruction approaches, (b) wall function bounce, (c) immersed virtual wall method, (d) wall velocity approaches.

2.1.8 Complex geometries

The main difficulty of applying the LBM to complex boundaries is the Cartesian grid, on which most LBMs are based. A variety of strategies exist to either circumvent or address this issue. A simple method is a stepwise approximation of the geometry and the use of simple bounce-back methods. This introduces an error in the approximation of the geometry, while keeping the boundary condition nominally second-order accurate. By refining the region around the boundary, the approximation error can be reduced.

If a better approximation of the geometry is necessary, interpolated bounce-back (IBB) methods are another option. IBB is based on the fact that populations travel exactly Δx during one time step. Therefore, a fictitious population is created by interpolation that will travel Δx and end up at the boundary node. Methods differ on the interpolation scheme. The original IBB proposed by Bouzidi et al. (2001) uses a linear or quadratic interpolation. The resulting scheme is independent of the collision operator and stable; however, it is not local if the boundary is closer than $\frac{\mathrm{\Delta }x}{\mathrm{2}}$. It is also not mass conserving, as shown in Krüger et al. (2017, p. 447). Furthermore, the distance at which the fictitious population has to be interpolated has to be known. In the case of a stationary boundary, this only has to be done once, but if the boundary moves, this computation has to be repeated at every time step. An immersed solid wall approach is presented in Feng et al. (2019b) that computes the equilibrium and non-equilibrium populations based on interpolation of macroscopic quantities. It requires multiple interpolations of macroscopic quantities and is therefore likely to be computationally expensive. However, an extension to include wall models appears straightforward, although not explicitly mentioned in the study.

Another class of boundary conditions is based on extrapolation methods. There exists a variety of methods with different advantages and disadvantages. The interested reader is pointed to Krüger et al. (2017, p. 455–463).

A different approach to implementing complex geometries is via the immersed boundary method (IBM). The idea of IBM, introduced in Peskin (2002), is to model the geometry via Lagrangian marker points that exert a force onto the fluid. The geometry is therefore completely decoupled from the lattice. This makes it suitable for complex and moving geometries, for example, in particle-laden flows. Since the IBM models the boundary via forces, the same methods can be used for Navier–Stokes solvers and LBM solvers.

IBMs can be split into two groups: explicit methods and direct forcing. In explicit methods, the marker points exert a force proportional to the virtual deformation of the boundary due to the fluid. This introduces a proportionality constant, which is effectively a spring constant that has to be determined. High values of the spring constant lead to higher accuracy of the geometry but can introduce instability, while the opposite is true for low values – see Krüger et al. (2017, p. 474f). The direct forcing approaches construct a force, such that the velocity post-collision has the correct value. A variety of ways to compute this force exist. The original variant by Feng and Michaelides (2005) computes the force explicitly by computing the force of the fluid on the node and then imposing a balancing force. Another approach, called the implicit velocity correction-based IB-LBM, computes the force of all marker points at the same time and is introduced in Wu and Shu (2009). This requires the inversion of a matrix of size N², where N is the number of marker points. However, this matrix depends only on the geometry of the boundary. Thus, if the boundary is rigid, it only has to be done once. The previously mentioned approaches all aim at imposing a no-slip condition at the boundary. However, if IBM is to be used to model a complex terrain in high Reynolds flows, a wall model has to be included as well. We could not find an example of IB-LBM used with wall model in the literature. The problem has been addressed with classical Navier–Stokes solvers, though, for example, in Chester et al. (2007).

2.2 Simulating wind turbines and wind farms in the LBM

To simulate wind turbines and wind farms in a flow solver, the effect of the turbine on the flow must be accounted for. Full-rotor simulations that resolve the geometry of the rotor require either the mesh to move with the rotor or to remesh when the rotor has moved. Both approaches have been applied in the LBM, and we will discuss these examples in Sect. 3.7.

However, often the exact flow at the blades is not of interest, and the rotor and tower can be modelled via body forces using actuator line (Sørensen and Shen, 2002) or actuator disc methods (Sørensen and Myken, 1992). Here, we only want to highlight aspects of these methods in relation to the LBM. More information on these methods can be found in, for example, Mikkelsen (2003), Troldborg (2009), and Martínez-Tossas et al. (2015). Since both methods rely solely on imposing forces on the flow, they can easily be introduced to the LBM.

Generally, the actuator line is considered more accurate. However, it is usually also more costly because it limits the maximal size of the time step, since the tip of the actuator line should not move more than one cell per time step (Troldborg, 2009). With $\mathit{\lambda }=\frac{\mathit{\omega }R}{V_{\mathrm{0}}}$ as the tip-speed ratio and ω as the rotational speed of the turbine, the time step is limited by

$$\begin{array}{}\text{(19)}& \mathrm{\Delta }t<{\displaystyle \frac{\mathrm{\Delta }x}{V_{\mathrm{0}}\mathit{\lambda }}}.\end{array}$$

With λ≈5–10 for modern turbines, this limit is significantly lower than the Courant–Friedrichs–Lewy (CFL) condition, namely $\mathrm{\Delta }t<\frac{\mathrm{\Delta }x}{V_{\mathrm{0}}}$, which is a typical stability criterion in implicit incompressible Navier–Stokes solvers. However, the time step size of the LBM is very small, as described in Sect. 2.1.2. Thus, for $\mathit{\lambda }<\frac{\sqrt{\mathrm{3}}}{\mathit{Ma}}\approx \mathrm{17}$, the time step limit for ALM is more relaxed than the limit for the LBM itself, assuming Ma=0.1. This makes the application of the ALM very attractive for the LBM. Indeed, to the best of our knowledge, only the ALM has been used in combination with the LBM so far, as we discuss further in Sect. 3.7.

To summarize, wind turbines can be simulated in the same way in the LBM as in Navier–Stokes-based solvers. However, one of the main drawbacks of the ALM in implicit Navier–Stokes solvers is not present in the LBM.

2.3 Thermal stratification

Since the LBM is based on kinetic theory and, as was shown before, the total energy E and internal energy e can be computed from the distribution functions, it would be fair to assume that the LBM can be extended to include temperature as well. While the full Navier–Stokes–Fourier equations can be derived from the Boltzmann equation, a discretization of such equations would require the use of more than 27 discrete velocities in three dimensions and the interaction between nodes that are not direct neighbours, as shown, for example, in Qian (1993). These so-called multi-speed methods suffer from high computational cost and instability (McNamara et al., 1995). However, in the ABL, it is sufficient to solve the Navier–Stokes equations and an advection-diffusion equation (ADE) for temperature and couple them via the Boussinesq approximation (Stull, 1988). Alternatively, Feng et al. (2021b) also develop an approach to incorporate the anelastic approximation into an LBM model. Transport equations for other scalars, such as humidity or turbulent kinetic energy, are sometimes considered in atmospheric flows as well (Deardorff, 1980; Vollmer et al., 2017). A variety of approaches have been developed to solve coupled Navier–Stokes and advection-diffusion equations with the LBM. The so-called hybrid methods discretize the advection-diffusion equation directly with finite-difference or finite-volume methods, while the double distribution function (DDF) approach solves the ADE using a second set of populations and an altered LBM scheme. The approaches can also be combined if multiple ADEs are being solved, in particular to reduce the memory cost of using the DDF approach.

A recent review of both approaches can be found in Sharma et al. (2020), therefore we include here only studies with particular relevance to our application or more recent studies.

The derivation of the hybrid methods is rather straightforward and was first introduced in Filippova and Hänel (1998b), and improved again in Filippova and Hänel (2000). Further progress is made in Lallemand and Luo (2003) by utilizing an MRT collision operator and a thorough analysis of instabilities in the thermal LBM. Three-dimensional simulations for convective flows were first conducted in Mezrhab et al. (2004) and showed good results for low Rayleigh and Reynolds numbers. An implementation of the hybrid model for a multi-GPU was proposed in Obrecht et al. (2012) and tested in Obrecht et al. (2013). Obrecht et al. utilize a Euler-forward time integration and a finite-difference scheme based on all direct neighbours, analogous to the discretization of the LBM. The authors stress, based on Lallemand and Luo (2003), that the finite-difference operator needs to have the same symmetry as the LBM stencil operator to improve stability. However, no further explanation is provided. A model combining the hybrid recursive regularized LBM with a finite-volume solver for scalars is presented in Feng et al. (2019a, 2021a). The finite-volume solver employs a Euler-forward discretization for the time integration and the MUSCL scheme (Van Leer, 1977) and central differences for the advective and diffusive terms, respectively. The authors use this model to simulate two active scalars, namely humidity and potential temperature. A similar approach with second-order Runge–Kutta time integration is used in Feng et al. (2021b) to solve the anelastic Navier–Stokes equations. The model is shown to yield accurate results in 2D test cases, but no 3D cases are shown.

In the DDF approach, the ADE is solved via the LBM. The temperature θ, or any other scalar, is the zeroth moment of the populations g_ijk,

$$\begin{array}{}\text{(20)}& \mathit{\theta }=\sum g_{ijk}.\end{array}$$

First-order moments of equilibrium are related to the advective flux,

$$\begin{array}{}\text{(21)}& \mathit{\theta }\mathit{u}=\sum {\mathit{c}}_{ijk}{g}_{ijk}^{\mathrm{eq}}.\end{array}$$

Note that in contrast to the LBM for the NSE, the first-order moments of the populations and the first-order moments of the equilibria are not equal. A similar relation of relaxation time τ to diffusivity α as Eq. (16) exists:

$$\begin{array}{}\text{(22)}& \mathit{\alpha }=c_{\mathrm{s}}^{\mathrm{2}}\left(\mathit{\tau }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\Delta }t\right).\end{array}$$

The application of the LBM to solve other than the Navier–Stokes equations and specifically ADE began shortly after the LBM was first developed. One early exploration of the LBM to solve ADE can be found in Wolf-Gladrow (1995). The author shows that an LBM scheme very similar to the one used for fluid flow is able to approximate the ADE. For a more extensive explanation of the LBM for ADE, the reader is referred to Krüger et al. (2017, p. 297–321), which also covers boundary conditions. Different collision operators can be used, and many of the operators proposed for the momentum equations have been adopted for ADE. An overview of proposed collision operators can be found in Gruszczyński and Łaniewski-Wołłk (2022), highlighting the fact that only few models have been proposed for three dimensions. Furthermore, the overview misses the model proposed in Yang et al. (2016), which is based on the factorized cascaded LBM. The same model is used by Adekanye et al. (2022) to simulate an active scalar.

In a wall-modelled LES, the temperature model also needs to employ an SGS and a wall model. A variety of models that compute an eddy diffusivity exist – see, for example, Stoll et al. (2020) and Gadde et al. (2021). Making use of the effective diffusivity approach, analogous to the effective viscosity approach and Eq. (22), enables the use of such models in the DDF approach. In finite-difference solvers, the application of a wall model is straightforward by applying Monin–Obukhov theory to compute the heat flux for a given stability, as done in Feng et al. (2019a, 2021a). In contrast to the Navier–Stokes model, it is possible to prescribe a heat flux with a bounce-back approach the same way as a wall velocity is specified. Another approach for a wall model is proposed by Kuwata and Suga (2021). The authors propose to use the same immersed virtual wall approach as described in Sect. 2.1.7 to model the wall heat flux as a source term on the first node.

In summary, the LBM can and has been augmented by a separate solver for temperature, either based on finite differences or an altered LBM scheme. Nevertheless, not many models have been proposed, especially in three dimensions. The development of the LBM for thermal simulations lags behind that of the isothermal LBM, especially for wall-modelled LES.

Figure 3 Overview of recent key advancements in research on simulating atmospheric flows using LBM-LES, presented as a timeline.

2.4 Summary of the available methods in the lattice Boltzmann method

The application of the LBM in wind energy is a relatively recent research field, with many aspects still actively being explored. To a large extent, this is due to the fact that a number of crucial advancements in the fundamentals of the LBM for wind energy applications only happened during the past decade. An overview of these milestones is collected in a timeline shown in Fig. 3. Today, the LBM is equipped with several accurate and sufficiently robust collision operators to handle high Reynolds number flows as needed for atmospheric boundary simulations or wind turbine wake flows. The same subgrid models that are used in Navier–Stokes-based solvers can be used in the LBM, thus enabling LES of atmospheric boundary layer flows. A number of wall-modelling approaches and refinement methods exist to enable large computational domains while reducing computational cost. Wind turbines can be simulated either as full-rotor simulations or using actuator line models, therefore also making it possible to simulate wind farms with a large number of rotors. Methods for extending solvers to include thermal stratification, forest canopies, and complex terrain exist, enabling the simulation of realistic conditions wind turbines are exposed to. Thus, the methods necessary for the LBM to become a comprehensive, reliable, and fast LES tool capable of meeting the increasing simulation requirements of wind energy research and industry exist.

3 Discussion

In the following section, we discuss different applications of the LBM in relation to wind farm flows. We gather studies to demonstrate how various aspects of the LBM have been examined, beginning from its capability to simulate highly resolved turbulent flows to atmospheric boundary layers and wind farm flows. At the same time, we identify areas where more research is needed. We structure our discussion around a set of questions we have asked ourselves or have been asked during the last years as we have been conducting research on the LBM.

3.1 Is the LBM suitable for turbulent flow simulations?

All flow phenomena relevant for wind farm simulations, be it ABL flows or wind turbine wakes, are highly turbulent. Thus, a numerical scheme used for these applications needs to be suitable for very high Reynolds numbers. Interestingly, this capability is still frequently questioned, particularly from people less familiar with the method, including the wind energy LES community. (It should be noted that we can only refer to anecdotal evidence for this scepticism. Examples are discussions about these capabilities at conferences or requests from reviewers to add fundamental validation cases to application-oriented papers, even though various previous publications have covered these topics extensively.) Indeed, the initial BGK collision model, and even later MRT approaches, do suffer from numerical instabilities at high Reynolds numbers. However, these issues have been remedied by modern collision models, as outlined in Sect. 2.1.3. In the following, we discuss these capabilities and provide numerical evidence to support that claim. We do so by highlighting several studies that perform direct numerical simulations (DNSs) and slightly underresolved simulations of turbulent wall-bounded flows.

Validations of LBM solvers against benchmark DNSs of various flow conditions have been successfully carried out since the early 1990s, affirming the LBM's suitability for bulk flow studies.

Jahanshaloo et al. (2013) provide a comprehensive review of studies validating LBM using fully or pseudo-spectral DNS. The pioneering work by Benzi and Succi (1990) examined turbulent flow using the BGK operator within a 2D square domain with periodic boundary conditions, evaluating LBM’s efficiency relative to spectral methods. At low Reynolds numbers, the study finds good agreement in both enstrophy and energy time series, as well as in the energy spectra. The authors also noted that, while both methods have comparable computational costs, LBM demonstrates improved scaling at higher grid resolutions.

Lammers et al. (2006) conducted resolved simulations of fully developed channel flows at a friction-based Reynolds number of ${\mathit{Re}}_{\mathit{\tau }}=\frac{u_{\mathit{\tau }}H}{\mathit{\nu }}=\mathrm{180}$, where $u_{\mathit{\tau }}=\sqrt{\frac{{\mathit{\tau }}_{\mathrm{w}}}{\mathit{\rho }}}$ is the friction velocity based on the wall shear stress τ_w and H is the channel half-height. They used the BGK method on a D3Q19 lattice, comparing the results with those from pseudo-spectral solvers. The results show remarkable agreement, concluding that LBM is as reliable as Chebyshev pseudo-spectral codes for DNS of turbulent flows, with the authors stating: “This removes any doubt that LBMs have inferior performance in resolved DNS.” Furthermore, they emphasized the LBM's numerical efficiency, achieving a 5× speed-up over the pseudo-spectral solver.

A comparison of various collision operators (BGK with D3Q19 and D3Q27, MRT, and cascaded LBM) was conducted by Freitas et al. (2011). The study finds that the BGK model showed better agreement with reference results in simulations of turbulent channel flows at Re_τ=200. However, this finding contrasts with results from their lid-driven cavity test, where the BGK models exhibits instabilities at higher Reynolds numbers, while MRT and cascaded models maintain good agreement. They concluded that universal stability assessments across methods remain challenging. This highlights the fact that the MRT operator is not able to improve the stability of the LBM in general.

Another comparison of collision operators was performed in Nathen et al. (2018), this time of the BGK, MRT, and regularized (Latt and Chopard, 2006) operators. In studying dissipation of the Taylor–Green vortex at varying Reynolds numbers and resolutions, they find that the BGK operator is accurate if the simulation is sufficiently resolved, while the MRT model is accurate in underresolved simulations, yet becomes unstable at too high resolutions. The regularized model is always stable but very dissipative.

The work of Gehrke et al. (2017) compared the performance of BGK, MRT, and cumulant collision operators in DNSs of turbulent channel flows at moderate Reynolds numbers (Re_τ=180). This study finds that DNSs with all three collision operators produce excellent results and that the cumulant operator is stable even at coarser grid resolutions, where BGK and MRT models exhibit non-physical oscillations. Gehrke et al. (2017) emphasized the resolution-dependent damping nature of the cumulant model, noting that it acts similarly to an inherent subgrid-scale model.

In conclusion, LBM has been used to simulate turbulence for 35 years and has been shown to be as accurate as traditional approaches based on discretizing the NSE, while being numerically significantly more efficient. However, it is only with the recent introduction of more advanced collision operators, especially the regularized and cumulant collision operators, that stability and accuracy at lower resolutions have become feasible in LBM.

3.2 Is the LBM applicable to LESs? Which lattice and collision operators are most suitable? How are subgrid scales accounted for?

A large number of LES studies have been conducted with LBM, utilizing explicit or implicit methodologies to model the subgrid scales, and employing various lattices and collision operators. In implicit LES, it is assumed that the effects of subgrid-scale motions are modelled accurately enough by the collision operator without additional explicit models. Again, we only highlight a few studies of particular importance to wind farm LESs.

A review of early studies of LES-LBM can be found in Jahanshaloo et al. (2013). The authors gathered a large collection of studies with different collision operators and subgrid-scale models. They find that, in general, boundary treatment is of high importance for accurate results.

In Kang and Hassan (2013), the authors examined the influence of the velocity discretization on the results of simulations of square ducts and round pipes. They find that a D3Q19 lattice is sufficient if the axes of the lattice align with the main axes of the walls. However, the lack of isotropy in the lattice negatively impacts the results if the duct is rotated 45° or in the round pipes. Similar problems related to the isotropy of the lattice have been reported in Asmuth et al. (2020b) when simulating wind turbine wakes.

The MRT D3Q19 model was applied, together with the Smagorinsky SGS model, to simulate a channel flow at Re_τ=180 by Wu et al. (2011). They furthermore simulated a passive scalar identified as temperature with an MRT D3Q7 model. The results match the reference DNS data reasonably well considering the differences in the numerical approaches.

Wang et al. (2014) also simulated a channel at Re_τ=180, using a BGK D3Q19 model with the Smagorinsky model for LES and DNS. They find an excellent agreement with the reference data.

Gehrke and Rung (2022a) presented results from simulations of a periodic hill with the cumulant LBM at different resolutions. They show that the well-conditioned parameterized cumulant LBM is capable of performing implicit LESs for highly resolved LESs by limiting the size of higher-order cumulants using so-called limiters (Geier and Pasquali, 2018). However, at very coarse resolutions, the damping of higher-order cumulants does not suffice and warrants the use of an additional SGS model. In a subsequent study (Gehrke and Rung, 2022b), the same authors developed a cumulant-based SGS model along a resolution-sensitive regularization. The authors simulated turbulent channel flows at ${\mathit{Re}}_{\mathit{\tau }}=\mathrm{180},\mathrm{550},\mathrm{2000}$ with 12, 24, and 48 nodes per channel height H using the cumulant LBM with a Smagorinsky SGS model and a wall function fitted to data from DNS. The authors find excellent agreement in mean velocities and shear stress even at low resolutions of 12 nodes per half-height, with the exception of the highest Reynolds number, where the mean velocities in the bulk are consistently underpredicted in the low-resolution case. Turbulence energy production is accurately predicted throughout the channel with the exception of the first or first two nodes. An additional simulation at Re_τ=5200 also agrees excellently with reference data obtained from DNS.

Spinelli et al. (2023) compared different collision operators and turbulence models. The MRT and HRR models were tested together with the Vreman and WALE SGS models, and the cumulant model was applied without any explicit SGS model. They simulated the flow past a cylinder at Re=3900, based on the cylinder diameter. The cumulant operator consistently outperforms other collision operators significantly and matches reference data very well. The authors reported that the MRT model is unstable with a D3Q27 lattice, and they therefore employed the MRT model with a D3Q19 lattice.

Overall, we find that LBM-LES has been applied successfully to wall-bounded flows, and even the most simplistic collision operators can produce accurate results at high resolutions. However, more advanced collision operators, especially the cumulant LBM, can produce highly accurate results also for more complex cases. If walls are aligned with the lattice directions, good agreement can be found with D3Q19; however, if the walls are misaligned or strong rotational features are present, D3Q27 is generally more suitable. At high resolutions, the implicit damping from the cumulant and HRR method suffice to achieve good results; however, at a coarse resolution, explicit SGS models are required. However, the interaction of the damping from the cumulant limiter or the hybrid finite-difference term of the HRR method with explicit SGS models has not yet been examined in detail. Some authors, for example, Asmuth et al. (2021), choose to deactivate the implicit damping when using an explicit model, while Spinelli et al. (2023) forwent the use of an explicit SGS model when using the cumulant limiter. In a direct comparison, the cumulant collision operator outperforms all other collision operators in terms of accuracy.

3.3 How do the different wall models for LBM perform in LESs?

Even with a suitable collision operator and an SGS model, wall models are required to reduce the computational cost of LES, be it of full-rotor simulations or the ABL. As discussed in Sect. 2.1.7, a number of models have been proposed in the literature, but we have not yet discussed how they perform.

The wall model approach introduced in Malaspinas and Sagaut (2014) was tested in the same paper in simulations of turbulent channel flows. The flow was modelled using a D3Q19 MRT LBM with a Smagorinsky SGS model, and the Musker's law is applied. Simulations at ${\mathit{Re}}_{\mathit{\tau }}=\mathrm{950},\mathrm{2000}$ show very good agreement in mean velocity with reference data. At Re_τ=20 000, the agreement is not as good but still satisfactory given the used resolution. Feng et al. (2021a) validated the generalized reconstruction method presented in Wilhelm et al. (2018) for simulations of atmospheric boundary layers using ProLB and the HRR model by simulating a neutral pressure-driven boundary layer. The mean velocities compare well to the log law and reference data in the vicinity of the wall. However, the momentum fluxes are consistently underpredicted compared to a range of reference cases. We will discuss further results presented in that publication later on.

Haussmann et al. (2019) compared different velocity boundary conditions and wall functions using a D3Q19 BGK model with a Smagorinsky turbulence model and a similar wall-modelling approach to Malaspinas and Sagaut (2014) but employed a three-layer wall function. They find that the reconstruction of the populations based on Guo's extrapolation scheme is significantly more accurate than the simple equilibrium reconstruction in simulations of turbulent channel flows at Re_τ=1000. At Re_τ=2000 and Re_τ=5200, results deteriorate at low resolutions, and at least 20 nodes per channel half-height are required. A systematic examination of wall function, SGS model, and collision operator is conducted by Spinelli et al. (2024), using the wall model approach from Haussmann et al. (2019). A turbulent channel flow at Re_τ=1000 was simulated with three different resolutions, namely 10, 20, and 40 nodes per H. The wall functions used comprise the Musker's, Reichardt's and power laws. We only want to report the findings of the collision operator here, where again the cumulant LBM without an explicit SGS model consistently outperforms the HRR and MRT models also tested.

The partial-slip-velocity-based model by Pasquali et al. (2020) was validated in the same publication, with simulations of turbulent channel flows at ${\mathit{Re}}_{\mathit{\tau }}=\mathrm{950},\mathrm{2000},\mathrm{16}\mathrm{000}$ and with resolutions of 10 and 20 nodes per H. They employed the cumulant collision operator without an explicit SGS model. They compared both approaches for computing the wall shear stress and find that computing the wall shear stress from Musker's law results in good agreement to reference data for Re_τ=950 and Re_τ=2000, while the cumulant-based wall model approach yields quite large deviations. At the highest Reynolds number (16 000), the cumulant-based approach yields better agreement at higher resolutions.

In Kuwata and Suga (2021), the authors compared both wall model approaches proposed in the same article to DNS of a turbulent channel flow at Re_τ=5200. They employed a D3Q27 MRT model, the Smagorinsky SGS model, and Musker's law. Both approaches yield good results for first- and second-order statistics; however, the immersed virtual wall approach consistently underpredicts the mean velocity. Due to the limitations of the specular reflection approach to straight walls, only the immersed virtual wall approach is further investigated. At low resolutions of 10 nodes per channel half-height, the accuracy of mean velocity prediction is reduced, while the shear stress is still in very good agreement with the reference solution. The model also performs well at lower (Re_τ=500) Reynolds numbers, whereas at higher (Re_τ=10 000), the underprediction of mean velocity increases. Finally, the effect of the thickness of the virtual wall was assessed, and it is observed that a lower thickness yields better agreement with the reference data. Xue et al. (2023) also implemented the specular reflection method by Kuwata and Suga (2021) and applied a D3Q19 MRT together with the Smagorinsky SGS model to simulate turbulent channel flows at ${\mathit{Re}}_{\mathit{\tau }}=\mathrm{1000},\mathrm{2000},\mathrm{5200}$, with 10, 20, and 30 nodes per H. (Note that in their publication, Xue et al. denoted the channel half-height δ and the channel height H.) The wall shear stress was computed from Reichardt's law. Notably, they employed a synthetic turbulence generation method at the inflow. Far enough downstream of the inlet, the results agree well with DNS reference data in mean velocities and stresses, even at low resolutions.

The WFB approach was evaluated in Han et al. (2021b) by simulating a turbulent channel flow at Re_τ=640 and Re_τ=2003 with a D3Q19 MRT and the Smagorinsky SGS model. The wall shear stress was computed from Spalding's law. Results show an improvement over not using a wall model; however, there is a consistent overprediction of mean velocity at the lower Reynolds number. At the higher Reynolds number, the error reduces. The wall model is further validated in Han et al. (2020) against measurements of the flow past a rectangular block, with a Reynolds number based on the block width of 40 800. SGS stresses are computed with the WALE model. The comparison to reference simulations and experimental data shows good agreement in mean velocities and turbulence kinetic energy (TKE); however, the authors report numerical oscillations attributed to the collision operator. In Han et al. (2021a), the effect of the wall model was examined in the same setup as the previous studies. They find that the application of the wall model reduces grid requirements and is able to produce results of the same quality as a simple bounce-back scheme at significantly coarser resolutions, reducing computational cost.

Finally, Asmuth et al. (2021) validated the wall model approach introduced in the same study in an isothermal pressure-driven boundary layer. The cumulant LBM, in conjunction with the AMD SGS model and Monin-Obukhov similarity theory, were used, and the flow statistics were compared to results from a pseudo-spectral solver. Different methods to determine the wall shear stress were compared in terms of their effect on the log-layer mismatch. The authors show that the elevated Schumann–Götzbach model of Maronga et al. (2020) yields the best results. In a subsequent grid sensitivity study, excellent agreement of mean quantities and stresses is found, and even higher-order statistics show decent agreement. The iMEM approach was further examined in Gehrke and Rung (2022b), which we already partially discussed in Sect. 3.2. Recall that the authors simulated turbulent channel flows at a range of Reynolds numbers and resolutions, and found excellent agreement in the bulk. Only at the highest Reynolds number do they find consistently lower velocities than the reference. Furthermore, the mean velocity and turbulence energy production at the first node off the wall is consistently too low for the two higher Reynolds numbers.

To conclude this section, we find that many wall model approaches have shown good agreement in the simple case of a turbulent channel flow. However, a consistent and direct comparison of the different wall-modelling approaches has not yet been conducted. Additionally, the effect on wall models from wall function, collision operators, and subgrid-scale models has not been explored in detail either. Furthermore, most wall models were only ever applied at Reynolds numbers, significantly lower than what is found in the atmospheric boundary layer, with the exception of iMEM model and the reconstruction method. Finally, not all wall models can be extended to complex geometries, and comparisons of wall models in complex geometries are missing entirely.

3.4 How are complex geometries handled in LBM in the context of turbulent flows?

Dealing with complex geometries, such as complex terrain surrounding wind farms, poses a number of challenges. Since the LBM is based on a Cartesian grid, it does not have the same flexibility to adapt the grid to the terrain as a finite-volume method. However, as we discussed in Sect. 2.1.4, interpolated boundary conditions can be used to approximate the shape of a geometry more accurately than a simple staircase approximation. Another challenge is the use of wall models. Some of the approaches presented in Sect. 2.1.7 can inherently be used to model to interpolated boundaries, while others can be adapted; however, this is not possible for all approaches. Finally, the data layout of a solver has to be suitable for the efficient representation of complex geometries.

An early example of complex geometries represented in a GPU-resident solver can be found in Bernaschi et al. (2010). The solver employs the indirect addressing scheme that can also be found in modern solvers such as VirtualFluids, allowing for an efficient representation of arbitrary geometries in a fashion suitable for GPUs. Jin et al. (2015) used a highly resolved, D3Q19 BGK-based simulation to study the effect of roughness on channel flows. The authors highlight the efficiency of the method, consequently allowing an even higher resolution. A more recent publication that considers the flow over a hill can be found in Schubiger et al. (2020). The authors assessed the ability of LBM to predict the flow over a well-studied benchmark using the open-source, CPU-resident solver Palabos, utilizing only interpolated bounce-back boundaries without a wall model. The authors find reasonable agreement between the LBM, reference simulations conducted with RANS and DES, and the experiments, despite the simplistic modelling setup. A comparison of the computation time shows that the LBM solver is about five times faster than the DES solver. The study highlights the continued need for an accepted wall-modelling approach for complex geometries. As already mentioned in Sect. 3.2, Gehrke and Rung (2022a) simulate the flow over a periodic hill using cumulant LBM with interpolated bounce-back boundaries but without turbulence or wall model at different resolutions. At the highest resolution, the results match the reference solutions very well and using different Mach numbers also has little influence on the quality of the solution. At lower resolutions, the quality of the solution deteriorates.

The extended reconstruction method for curved boundaries was used in conjunction with the HRR model to simulate flow around an NACA0012 airfoil in Degrigny et al. (2021), and excellent agreement to reference data is observed. In addition to the turbulent channel flow, Haussmann et al. (2019) also conducted simulations of a Coriolis mass flow meter using the wall model, bounce-back, and interpolated bounce-back. Note that unless the interpolated bounce-back is applied, all walls were approximated with a staircase approximation. Nevertheless, the wall model yields generally good agreement with the measured pressure drop across the mass flow meter.

In Kuwata and Suga (2021), the authors apply the IVW model described in Sect. 2.1.7 to the periodic hill test case in addition to the channel flow discussed in Sect. 3.3. Musker's law was used to compute the shear stress, and a D3Q27 MRT combined with the shear-improved Smagorinsky model was used for the bulk flow. The results of the simulations of the periodic hill mostly agree well with the reference solutions; however, the recirculation areas exhibit significant differences in the mean velocities and Reynolds stresses. Overall, the results are of similar quality as those reported by Gehrke and Rung (2022a) with the same resolution, who did not use a wall model.

Overall, we find that the representation of complex geometries is well established; however, simulation of real-world complex terrain remains sparse. Furthermore, only two studies use a wall model in combination with complex geometries. This clearly presents a gap in the current research and is likely related to the fact that many wall-modelling approaches have only been formulated for straight boundaries.

3.5 Is the LBM suitable to simulate large domains and the atmospheric boundary layer?

Despite the lack of wall models in complex geometries, a rapidly growing number of studies of wind in the atmospheric boundary layer have been presented in recent years. Many of these studies focus on urban flows, where wall models are typically not as important since drag terms often dominate the flow (Lenz et al., 2019). However, these studies demonstrate the suitability of the LBM for wind energy research and are therefore also included here.

The earliest application of LBM to urban flows was reported in Fan et al. (2004), which is also the first implementation of LBM on a cluster of GPUs. The BGK operator without any turbulence model was used. The methodology is very simplistic, yet many more recent studies use essentially the same methodology, albeit with the inclusion of a turbulence model. The authors focussed on the computational performance and report almost 50 MNUPS on a cluster of 32 NVidia GeForce FX 5900 Ultra GPUs, each with 128 MB of memory. We mention this study to highlight how far the methodology and hardware have come in the last 20 years. It also demonstrates the early embrace of GPUs by the LBM community.

Modern applications of LBM to urban flows began with Onodera et al. (2013), where simulations of a 10 km×10 km domain with a grid spacing of 1 m of the urban area of Tokyo are presented. A D3Q19 BGK model was used, walls are modelled with the bounce-back scheme, and the coherent structure SGS model was used. Scaling tests up to 1000 GPUs were performed, and good scalability is obtained. Additionally, a mesh with $\mathrm{10}\mathrm{080}\times \mathrm{10}\mathrm{240}\times \mathrm{512}=\mathrm{5.28}\times {\mathrm{10}}^{\mathrm{10}}$ nodes was simulated on 4032 GPUs. It marks one of the largest LESs to this day. The same setup was used by Ahmad et al. (2017) and Inagaki et al. (2017), presenting results on wind gust index, and turbulence statistics and structures, respectively. Further improvements of the model with respect to scaling across a very large number of GPUs are presented in Onodera et al. (2018) and Onodera and Idomura (2018), using adaptive mesh refinement and techniques to reduce internode communication. These studies utilized the cumulant collision operator. They report excellent computational performance and decent scalability across a very large number of GPUs, highlighting the LBM's suitability for very large simulations.

Watanabe et al. (2020) discussed simulations of plant canopies based on central moment collision operators. The forest canopy was represented by a drag force model. Comparisons with reference simulations confirm that the approach is capable of reproducing the canopy flow very accurately.

The study by King et al. (2017) employed a D3Q19 BGK model to simulate the flow through a building facade. The comparison with experiments and a traditional CFD code is within the expected variation, while the LBM approach significantly reduces the computation time.

In Lenz et al. (2019), the cumulant LBM solver VirtualFluids was used to study the feasibility of real-time simulations of urban flows by comparing to measurements made in the Basel UrBan Boundary Layer Experiment. Boundaries are represented, with the bounce-back method and mesh refinement used to reduce the computational cost. Good agreement with the measurements is reported, and quality criteria for this simulation are reached.

In contrast to the previous studies, Jacob and Sagaut (2018) utilized the HRR collision model implemented in ProLB to simulate an urban flow. No subgrid-scale model was used, and the study mentions the use of a wall model but does not give any details regarding its implementation. The results were compared to field and wind tunnel measurements in several points and show a fair agreement. Further studies of urban and micrometeorological flows are presented in Jacob et al. (2021).

Buffa et al. (2021) also used the HRR approach with the wall model from Malaspinas and Sagaut (2014) and a synthetic eddy model for the inlet to study the wind loads on a single high-rise building. The authors find very good agreement with experimental data, given an adequate setup.

A study examining forest canopy modelling in the LBM was conducted by Shao et al. (2022). The study employs a BGK D3Q19 model with a Smagorinsky model. The forest is modelled with a drag term. A simulation of a boundary layer with a forest on the ground was compared to reference simulation, and reasonable agreement is found in the mean velocities, while noticeable discrepancies in the fluctuations are present.

A major step towards enhancing the applicability of LBM-LES for wind energy problems is to establish a reliable coupling between LBM-LES and mesoscale atmospheric simulations. This challenge is not unique to LBM; in fact, coupling mesoscale models with NS-LES has been an active research topic for many years (Haupt et al., 2019). Although reviewing these techniques lies outside the scope of this work, it is worth noting that such attempts remain rare for LBM-LES. This is likely because the use of LBM for atmospheric flows is still relatively recent. However, Onodera et al. (2021) have explored such a coupling by interfacing WRF with LBM-LES. Their approach introduces an extended nudging and data-assimilation zone around the LBM-LES domain, where WRF data are interpolated in space and time, and additional turbulence is generated via roughness blocks. At the boundaries, the authors reconstruct the equilibrium populations directly from the mesoscale fields, while the non-equilibrium components are derived from the turbulence developed in the nudging region. This suggests that employing LBM-LES in a mesoscale-coupling framework does not introduce fundamentally new difficulties compared to conventional NS-LES/mesoscale coupling but rather presents the same challenges. Applying the LBM itself for mesoscale simulations is hindered by its limitation to cubic cells.

Overall, we find a large increase in the number of studies over recent years. Most studies employ the D3Q19 MRT model with a Smagorinsky subgrid-scale model. However, groups with a longer history tend to more sophisticated collision operators and the D3Q27 models. The HRR model is also well established but only implemented in a single solver. Furthermore, despite claims from various authors, the actual computational efficiency of HRR is not reported in any of the papers reviewed here. Different approaches to thermal coupling are used throughout the literature, with DDF and hybrid approaches being used about equally; neither show great advantages or disadvantages over the other. Most studies reviewed do not employ a wall model, and only two wall models have been validated in simulations of an atmospheric boundary layer. Furthermore, none of the reviewed studies utilize a precursor method and instead rely on uniform velocity or periodic boundary conditions. Since the resolution and domain size are very similar to wind farm simulations, LBM is shown to be suitable for simulations of wind farms.

3.6 How do the two thermal LBM approaches perform, and has the LBM been applied to stratified ABLs?

An isothermal description of the atmospheric boundary layer is often sufficiently accurate for simulating individual wind turbines or small farms; however, the role of thermal stratification on various aspects of wind farm flows is gaining much attention. Therefore, a potential LBM solver should also have the capability of simulating thermally stratified flows. As discussed in Sect. 2.3, two approaches are commonly used to simulate thermally stratified flow: the DDF and the hybrid approach. As the review by Sharma et al. (2020) notes, the vast majority of thermal LBM applications are limited to 2D flows.

The study by Wu et al. (2011) simulated the transport passive scalar in a turbulent channel flow using a DDF approach with the MRT collision operator, as mentioned previously. They compared their results with DNS simulations and show a reasonable agreement.

A similar methodology was used by Ren et al. (2018); however, the scalar was coupled via the Boussinesq approximation, and the Vreman SGS model was used. Simulations of a turbulent channel flow at Re_τ=180 and varying Rayleigh numbers were compared with reference data obtained from DNS. The authors find good agreement, even in second-order statistics.

A more advanced collision model for the second populations, namely a cascaded model on a D3Q15 lattice, was employed in Hajabdollahi and Premnath (2018) to study natural convection at Rayleigh numbers up to 10⁵ by means of DNS. They find very good agreement with reference data. Guo et al. (2023) used a DDF approach to conduct a DNS of a stably stratified flow. They employed the BGK collision operator for both distributions, with a D3Q7 lattice for the temperature, a D3Q19 for the momentum, and apply the Boussinesq approximation to couple momentum and temperature. The results agree very well with DNS results obtained with a pseudo-spectral solver. A comparison of the computational performance reveals that the pseudo-spectral solver is significantly more efficient, due to looser requirements in mesh generation and larger time steps. Kuwata and Suga (2021) present the only simulation of a DDF approach that also includes a wall model. However, the scalar is passive. The results of a turbulent channel flow agree well with the reference data.

Watanabe et al. (2021) simulated the transport of a passive scalar through a plant canopy with the DDF approach by augmenting the model presented in Watanabe et al. (2020) with a second distribution function modelled with a single relaxation time operator. Simulation results were compared to that of a Navier–Stokes-based solver, and excellent agreement is found in the concentration of the scalar. An extension of the model used in Wang et al. (2020a) to simulating stratified flows over a ridge is presented in Wang et al. (2020b). The temperature field was simulated using a D3Q7 MRT model. Results exhibit the expected physical behaviour, and compared well to laboratory and numerical results. However, significant oscillations are present in the flow field.

In Onodera et al. (2021), the cumulant method was augmented with a finite-difference solver for temperature and concentration. The advection-diffusion equations were discretized using second-order finite differences in space and explicit Euler-forward in time. The evolution of a plume of concentration in an urban environment is simulated. The simulation was driven by large-scale meteorological data obtained from the Weather Research and Forecasting model. Results were compared with results from a measurement campaign, and decent agreement is found.

In Feng et al. (2019a), ProLB was extended with a hybrid approach to solve temperature and humidity equations with a finite-volume approach as described in Sect. 2.3. Excellent agreement to reference data is reported. The previously mentioned work from Feng et al. (2021a) gives an extensive description of the models used in ProLB and presents validation studies. Under stable and convective conditions, results agree well with the reference data despite lower resolution. Various other benchmark cases with more complex conditions with fair agreement are also reported. Jacob et al. (2021) applied the method described in Feng et al. (2021a) to a range of benchmark cases of urban and micro-meteorological flows, and report decent to good agreement throughout. The authors also report high computational efficiency, although metrics are not provided.

To summarize, we find that the double distribution function and hybrid approach have been used to model temperature fields and have been applied to model stratified atmospheric boundary layers. While Sect. 2.3 showed that there is much more theoretical discussion of the DDF model, the hybrid approach is preferred in application to stratified boundary layer flows.

3.7 How has the LBM been applied to LESs of wind turbines and wind farms?

So far, we have discussed many individual aspects of LES related to wind energy. However, combining all of the individual methods and models, such as a turbine model, wall model, and SGS model, in a single solver is not always straightforward. Furthermore, some models depend, for example, on specific collision operators and are maybe incompatible with others. Nevertheless, we found that a number of publications employ the LBM to study either single wind turbines or wind farm flows. Interestingly, wind turbines were first simulated as full-rotor simulations, and actuator models were only introduced later on.

The first application of LBM to wind turbine simulations was conducted in Perot et al. (2012). The authors simulated a fully resolved NREL Phase VI model with the commercial LBM solver powerFLOW, compared the results to experimental data, and found good agreement of blade forces and wake quantities. Xu (2016) conducted simulations of up to three NREL Phase VI model turbines using the same solver. The results were compared to experimental data, and qualitatively good agreement was found; however, a quantitative comparison was not reported. Deiterding and Wood (2016a, b) present results obtained from conducting rotor-resolved simulations of the NREL Phase VI model turbine and a Vestas V27, therefore reporting the first simulation of an, albeit small, full-scale turbine. They employed the AMROC solver, which includes adaptive mesh refinement and a third-order scheme for the equilibrium distribution. They reported good agreement to the experimental data in all cases.

The first simulation of a wind turbine using an actuator line method was presented in Rullaud et al. (2018) using a 2D LBM solver based on the MRT collision operator and validating against measurements of two vertical axis wind turbines. A first validation study of a 3D LBM simulation coupled with the ALM was presented in Asmuth et al. (2019), employing the cumulant LBM implemented in the GPU-resident LBM solver elbe without an explicit SGS model. The authors studied the 5MW turbine in uniform inflow conditions. Parameter studies of the smearing width and the Mach number are presented and compared to results obtained with the Navier–Stokes solver EllipSys3D. The results of the velocity deficit in the wake show excellent agreement. Asmuth et al. (2020b) used the same methodology, albeit with a Smagorinsky SGS model with uniform and turbulent inflow, providing a more thorough examination of the wake. Again, excellent agreement to the reference solution obtained with EllipSys3D was found. Compressibility in LBM-ALM simulations was studied by Asmuth et al. (2020a). The effect on the wake was found to be negligible, although changes in density were observed in the vicinity of the actuator line. The methodology was compared to measurements and other solvers in Asmuth et al. (2022). Two benchmark cases were defined based on measurements from the DanAero experiments, featuring two 2.3 MW turbines. The LBM simulations feature the iMEM wall model. The authors conclude that the models were able to capture integral quantities such as torque and thrust well, even in waked inflow, and highlight the difficulty of generating inflow conditions for LES based on met-mast measurements. The LBM model performs similar to the other solvers presented in the benchmark. The same methodology, but implemented in the GPU-resident LBM solver VirtualFluids, was further validated in Korb et al. (2023b), where we presented validation against the SWiFT benchmark, including turbine response and wake measurements. The setup utilized a precursor to generate realistic inflow, the induction correction model proposed by Meyer Forsting et al. (2019) and the iMEM wall model. We find very good agreement in the turbine response and good agreement with the wake measurements, on par with the results obtained by traditional methods in the original benchmark by Doubrawa et al. (2020). In Asmuth et al. (2023), we examined the requirements for the industrial adoption of LES in terms of capabilities and computational performance, and compare the results with a simulation of 64 DTU 10 MW turbines (see Bak et al., 2012 for the definition) using the same methodology as described in Korb et al. (2023b).

Schottenhamml et al. (2022) presented an implementation of the ALM in the LBM code waLBerla and compared it to simulations conducted with SOWFA. The study also employs the cumulant LBM in combination with the Smagorinsky model. The actuator line model was modified to use two alternative formulations for the Gaussian smearing kernel. Simulations of the DTU 10 MW reference turbine were compared. Excellent agreement of blade forces was reported, while the velocity deficit in the wake differs slightly. The study features scaling experiments up to 40 GPUs and finds that waLBerla scales well to multiple GPUs across a number of nodes, provided the workload is large enough. The approach was further developed by Schottenhamml et al. (2024), where the algorithm was detailed, and validation via comparisons with the NewMexico experiment was presented. Furthermore, the excellent scaling of waLBerla to a total of 120 GPUs is demonstrated, thus proving its suitability to simulate even very large wind farms, although the scaling tests only feature a single turbine. A further implementation of an ALM in a cumulant LBM solver was presented in Watanabe and Hu (2024). The authors presented results of a weak scaling test, scaling from four to 256 GPUs, assigning each GPU the same number of fluid nodes and a single turbine. They demonstrated excellent scaling of the fluid solver. However, the calculations related to the ALM were not parallelized and therefore do not scale well. Furthermore, a simulation of a single Vestas V80 turbine with a diameter of 80 m was presented. The results compare fairly well to results from a finite-difference Navier–Stokes solver. Finally, results from simulating eight Vestas V80 turbines in a row are presented, and the importance of a rotor-speed controller is demonstrated. The same model was used in Watanabe and Hu (2025) to simulate the Horns Rev I wind farm comprising 80 Vestas V80 turbines under 51 inflow directions. Each simulation features 1.65 billion nodes and was conducted on 32 GPUs. In Watanabe et al. (2026), the same authors extended their model to use the interpolated bounce-back approach to model nacelle and tower of the turbine, and compare their results to results from the NTNU blind test wind tunnel study (Krogstad and Eriksen, 2013) with good agreement to the reference data.

An actuator line embedded in a rotating mesh was presented by Ribeiro and Muscari (2023). Simulations of the NREL Phase VI rotor and the NREL 5 MW turbine with powerFLOW are presented and achieve excellent agreement with reference simulations. However, the presented case is somewhat unusual in that the domain includes a very fine mesh around the blade, with resolutions much higher than is typical for ALM simulations. The authors also conclude that the approach is feasible for highly resolved simulations but might not merit its application in the typical scenario of a coarse mesh. The same sliding mesh approach was used in Ribeiro et al. (2025), comparing blade-resolved simulations and highly resolved actuator line simulations to measurements obtained from a rotor in a water channel. The study finds excellent agreement between the blade-resolved simulations and the actuator lines. Furthermore, asymmetric rotors and the effect of asymmetry on tip vortex breakdown were examined.

In a study by Qiu et al. (2025), the authors implemented a modified actuator line method, where the blades, hub, and tower of the turbine were modelled as rigid bodies with constrained joints and torque, and pitch and yaw controllers. Furthermore, they employed the entropic LBM. No data on the computational performance were provided. The resulting blade loads feature some discrepancies to reference data near the tip. However, the wake deficit agrees reasonably well. Wang et al. (2025) implemented the actuator line method in ProLB, and validations of single wind turbine and wind farm flows under different thermal stability are presented. Validation of the wake of a single turbine in neutral, convective, and stable boundary layers yields good agreement with reference data from another LES. Simulations of a model wind farm in a wind tunnel comprising 30 turbines also agree closely to measurement data. The same wind farm was also simulated in convective and stable boundary layers, and the effect of stability on the wind turbine wakes was investigated. This marks the first study simulating wind turbines and wind farms in thermally stratified boundary layers using the LBM. However, the chosen setups do not include a capping inversion as is present in the atmospheric boundary layer, thus not representing realistic atmospheric conditions.

The high computational efficiency of LBM makes it a prime candidate for use in combination with data-driven and machine learning methods. Two studies have utilized the LBM in combination with machine learning approaches. In Korb et al. (2021), we used the methodology developed in Asmuth et al. (2020a) to create a training environment for a deep reinforcement learning approach to optimize a wind farm controller. While we ultimately failed to improve the wind farm performance, we showcased how LBM enables new data-driven methods. In a second study (Asmuth and Korb, 2022), a similar setup was used to generate training data for a convolutional neural network to predict the average velocity and turbulence intensity in the wake of a single turbine. We find excellent agreement with the test data, while the computational time to infer a single case was negligible.

To give an overview of the collision operators employed in wind energy, we conducted a systematic review. The methodology and selected papers are listed in the Appendix A. Although not exhaustive, this sample is considered sufficiently representative to provide an overview of the current use of collision operators in wind energy LBM research. The results are shown in Fig. 4. The BGK model is widely used for studies of model turbines and vertical axis wind turbines (VAWTs), i.e. flows with a lower Reynolds number. However, for simulations of full-scale horizontal axis wind turbines, the cumulant operator is by far most commonly used. No studies of full-scale VAWT were found.

Figure 4 Distribution of collision operators employed in wind-energy-focused LBM literature. Results are differentiated by type of turbine simulated and simulations of full-scale or lab-scale conditions.

To conclude, a variety of different approaches for using the LBM for wind farms have been presented in the literature, with the most thoroughly studied approach being the cumulant LBM in combination with an actuator line. However, many of these studies are conducted in idealized conditions. Recently, studies in more realistic conditions have been presented, including wall models and realistic inflow conditions, either through synthetic turbulence or precursor–successor setups. One study has also simulated wind farms that consider thermal stratification, albeit not fully representing the complexity of the atmospheric boundary layer due to the lack of a capping inversion. Furthermore, several studies simulating a large number of turbines have been presented and have proven the feasibility of LBM for simulating large wind farms. Even under more realistic conditions, the method remained computationally very efficient and is thus a viable candidate for the industrial application of LES.

Figure 5 Evolution of computational performance per single GPU over time, distinguished by GPU architecture generation. Unfilled markers indicate the use of a D3Q19 stencil, while filled markers represent the use of D3Q27.

3.8 Which computational performance can be expected from LBM solvers, and are they really faster than Navier–Stokes-based solvers?

As highlighted multiple times throughout this paper, the exceptional parallelization capabilities of the LBM makes it particularly well suited for GPU implementation. The collision operator of the LBM has a low arithmetic intensity, meaning that computational performance is mostly limited by memory bandwidth (Harris, 2005; Lehmann et al., 2022). Collision operators that increase arithmetic intensity can therefore be implemented without reducing computational cost. On the other hand, choosing a smaller lattice increases the computational performance since less data need to be loaded per grid point. Since the first application of LBM on a GPU cluster in 2004, substantial advancements in hardware have led to a significant increase in computational power. Million node/lattice updates per second (MNUPS/MLUPS) are a common measure of computational efficiency in the LBM community. MNUPS can also be converted to calculate the ratio of simulated time-to-wall time ($T_{\mathrm{sim}}/T_{\mathrm{wall}}$) by a simple relation:

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{T_{\mathrm{sim}}}{T_{\mathrm{wall}}}}={\displaystyle \frac{\mathrm{\Delta }tP}{N}}={\displaystyle \frac{\mathrm{\Delta }x\mathit{Ma}P}{\sqrt{\mathrm{3}}{V}_{\mathrm{0}}N}},\end{array}$$

where Δt is still the size of the time step, P is the computational performance, and N is the number of grid points. The second expression is obtained by employing Eq. (14). Figure 5 illustrates the evolution of computational performance, expressed in MNUPS per GPU over time and with a distinction between GPU architecture generation. The data displayed are gathered from a compilation of LBM studies mainly focused on atmospheric flow simulations on GPUs, and which we have already discussed or presented earlier in this paper (see Table B1 for a list of the included studies). We clearly observe a steady increase in computational performance and that the current state of the art performance is 2000–3000 MNUPS. Therefore, given Eq. (23) and assuming Δx=10 m, $V_{\mathrm{0}}=\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, $N=\mathrm{50}\times {\mathrm{10}}^{\mathrm{6}}$, P=2000 MNUPS, and an artificial Ma=0.1, $\frac{T_{\mathrm{sim}}}{T_{\mathrm{wall}}}=\mathrm{2.3}$, thus simulating more than twice as fast as real time – which is in line with values reported, for example, in Asmuth et al. (2023). Furthermore, we find that the D3Q27 lattice has now replaced D3Q19 as the preferred choice. We have also computed a best-fit logarithmic regression, indicating that, on average, computational performance has doubled approximately every 2.75 years. While Fig. 5 depicts the computational performance per GPU, it is also important to address the efficiency of multi-GPU configurations. In early studies (Fan et al., 2004; Onodera et al., 2013), multi-GPU implementations suffered from significant performance losses. However, advancements in parallelization algorithms have substantially mitigated these inefficiencies, rendering the associated losses almost negligible in modern implementations (Onodera et al., 2018; Schottenhamml et al., 2022). Thus, LBM is especially suited for very large simulations. Furthermore, there has been a very rapid increase in the available memory per GPU, thus enabling fairly large simulations to be conducted on a single node, as noted, for example, in Asmuth et al. (2023), where a single data centre GPU was used to simulate a wind farm consisting of 64 turbines.

Direct comparisons to Navier–Stokes-based solvers in terms of computational performance remain sparse to this day. Furthermore, it is difficult to establish how such a comparison should be carried out since a traditional finite-volume code has very different requirements than an LBM solver. The same quality of solution typically requires more cells in an LBM solver due to the fact that all cells have a fixed aspect ratio and stretching is not possible. This difficulty is exacerbated by the fact that many LBM solvers run on GPUs, requiring few compute nodes, whereas most Navier–Stokes solvers run on massively parallel CPU systems.

A direct comparison of a CPU-based, explicit finite-difference solver for the weakly compressible Navier–Stokes equations and a CPU-based LBM solver can be found in Wichmann et al. (2021). For the solution of a transient problem, both codes require similar computational time. In terms of time to solution, Asmuth et al. (2020b) compared results from an LBM simulation to those from EllipSys3D, a CPU-based finite-volume solver and found that the wall time of the LBM simulation run on a single-consumer GPU is only a third of the FV solver, executed on a compute cluster. In terms of processor time, a reduction on the order of 10³ is found. Schottenhamml et al. (2022) compared a CPU and GPU solver to SOWFA and report reductions in wall time of a factor of 73 and 473, respectively. However, wall time is strongly dependent on the number of processors used. A more robust measure is the energy to solution, which is total electric energy required for the solver to arrive at the solution. In Korb et al. (2023b), we compared to many different solvers in terms of energy to solution and find that our solver reduces energy to solution by 2 to 3 orders of magnitude. Nevertheless, it still requires an order of magnitude more energy than a RANS solver.

So far, every comparison has shown that the LBM on GPU is significantly more efficient than a CPU-based Navier–Stokes solver. Comparisons to GPU-based Navier–Stokes solvers have not yet been performed and are difficult to do based solely on values reported in literature. For example, Sauer and Muñoz-Esparza (2020) report that FastEddy achieves roughly 90 time steps per seconds for a grid of approximately 1.25×10⁵ nodes on a single NVidia V100 GPU. Converting those numbers to the familiar metric in LBM suggests a performance of around 123 MNUPS, which is lower than all LBM solvers gathered here on similar hardware. However, the physics considered in the simulation are also significantly more complex. A fair comparison would have to be done on the same, or at least very similar, hardware and simulate a clearly defined benchmark case. Furthermore, the exact metric on which to evaluate the benchmark will have to be established, since the LBM has other limitations on the spatial and temporal discretization as finite-volume or finite-difference solvers. Therefore, solvers must be compared in terms of metrics, taking into account the trade-off between computational effort and accuracy, such as the energy to solution vs accuracy metric presented in Doubrawa et al. (2020).

Another practical concern for simulating large domains is the available memory. Naïve implementations of the LBM require two sets of population per node, one each for pre- and post-collision populations. However, advanced streaming algorithms, such as the AA pattern (Bailey et al., 2009) or the esoteric schemes (Geier and Schönherr, 2017; Lehmann, 2022), can avoid the need for a second set of populations. Another common approach is employing single-precision numbers or even custom representations (Lehmann et al., 2022). Memory requirements are seldom noted in publications, but there are several examples in the literature. In Asmuth et al. (2023), we mention that a grid of 187 million nodes fits in 40 GB of memory, obtaining a memory footprint of 0.21 GB per million nodes for a solver using esoteric twist and single-precision floating point numbers. Similarly, Watanabe and Hu (2024) report a memory footprint of less than 0.36 GB per million nodes. For comparison, Sauer and Muñoz-Esparza (2020) note that a grid of 32 million nodes fits onto a single GPU with 32 GB memory, while a grid of 72 million does not, resulting in a memory footprint of more than 0.44 GB per million nodes. Again, the comparison is insufficient as Sauer and Muñoz-Esparza (2020) considered significantly more complex physics. Nevertheless, it indicates that the memory requirements for Navier–Stokes and LBM solvers are of a similar order of magnitude.

4 Conclusions

This review has gathered the methods necessary to implement a solver for wind farm flows in the atmospheric boundary layer based on the lattice Boltzmann method. We have described various collision operators for the bulk flow, necessary boundary conditions, subgrid scale models, and wall-modelling approaches, and how wind turbines can be represented in the flow. We presented a short overview of methods to extend the isothermal LBM to thermally stratified flows. Finally, we gathered applications of the previously described methods, highlighting relevant studies of wall-bounded flows, both DNS and LES, as well as complex geometries and thermally stratified flows, and reviewed simulations of single wind turbines, wind farms, and atmospheric boundary layer flows.

We find that a variety of methods suitable for high Reynolds number flows exist and that the LBM has been successfully applied to such flows many times. Subgrid-scale models from the Navier–Stokes world can be readily applied in the LBM, and the LBM is well suited for simulating turbulence, with both DNS and LES. Wall modelling posed a challenge for a long time, but a variety of approaches have been proposed in recent years and successfully tested, even in isothermal atmospheric boundary layers. While complex geometries are relatively simple to simulate with the LBM, only two approaches to incorporate wall models have been extended to curved boundaries, and applications of these methods in the context of atmospheric boundary layers remain sparse. The actuator line method is well suited to the LBM due to the small time step of the LBM, and has been used in various studies of single wind turbines and wind farms of up to 80 turbines. However, we could not find any study applying an actuator disc method. Approaches to simulate thermal stratification exist, but we could find only a small number of studies actually using these approaches in 3D flows. One exception is a number of studies conducted with the solver ProLB, which has been used to simulate thermally stratified atmospheric boundary layers in a number of studies, including one study of a wind farm. The numerical performance of the LBM has steadily improved with improved hardware, and the LBM outperforms Navier–Stokes-based solvers consistently due to its parallel nature.

In summary, we find that all of the methods required for wind farm simulations with LBM-LES have been developed and validated, and that LBM-LES is well suited for large-scale simulations, especially using massively parallel hardware. Nevertheless, no simulation combining all of the discussed features – that is, a wind farm in a thermally stratified atmospheric boundary layer – has yet been conducted.

Appendix A: Collision operators used for wind energy applications

To identify a large body of papers using LBM for wind energy, we conducted a Google Scholar search using two sets of keywords: “LBM” and “wind energy” for the first set, and “LBM” and “wind turbine” for the second. The papers listed in Table A1 include all articles up to and including the 10th results page. Studies that did not appear in this search but which are cited in this paper were also included. To define the scope of wind energy research and apply basic quality control, several categories of papers were excluded:

• Studies focusing on urban flows,

• Articles addressing only the acoustic properties of wind turbines,

• Theses and manuscripts,

• Atmospheric studies that do not explicitly target wind energy applications,

• Research on related but distinct technologies (e.g. helicopters, propellers, tidal turbines),

• Studies conducted at Reynolds numbers <100 000,

• Studies in two dimensions.

Asmuth et al. (2020a)Asmuth et al. (2021)Asmuth and Korb (2022)Asmuth et al. (2023)Asmuth et al. (2022)Asmuth et al. (2020b)Asmuth et al. (2019)Cui et al. (2024)Deiterding and Wood (2016a)Deiterding and Wood (2016b)Dighe et al. (2020)Feng et al. (2019a)Feng et al. (2021a)Khan (2018)Korb et al. (2023a)Korb et al. (2023b)Korb et al. (2021)Laloglu and Alpman (2023)Li et al. (2021)Perot et al. (2012)Qiu et al. (2025)Ribeiro et al. (2025)Ribeiro and Muscari (2023)Schubiger et al. (2020)Schottenhamml et al. (2024)Schottenhamml et al. (2022)Shao et al. (2022)Venkatraman et al. (2023a)Venkatraman et al. (2023b)Venkatraman et al. (2022)Wang et al. (2025)Wang et al. (2024)Watanabe and Hu (2024)Watanabe et al. (2020)Watanabe et al. (2024)Watanabe and Hu (2025)Watanabe et al. (2025)Wood and Deiterding (2015)Wu et al. (2021a)Wu et al. (2021b)Xu et al. (2024)Xu (2016)Xue et al. (2021)Ye et al. (2020)Yuan et al. (2021)Zhiqiang et al. (2018)

Table A1 Articles used to create Fig. 4.

Appendix B: Computational performance

Fan et al. (2004)Bernaschi et al. (2010)Obrecht et al. (2013)Obrecht et al. (2015)Onodera et al. (2018)Onodera and Idomura (2018)Lenz et al. (2019)Asmuth et al. (2020b)Asmuth et al. (2021)Onodera et al. (2021)Schottenhamml et al. (2022)Shao et al. (2022)Asmuth et al. (2023)Watanabe and Hu (2024)Schottenhamml et al. (2024)Watanabe and Hu (2025)

Table B1 Publications listing achieved computational performance on NVidia GPUs.

^a When available, the performance of a single GPU in single precision was reported. However, some studies only provided the performance of GPU clusters. In such cases, we estimated single-GPU performance by dividing the cluster's performance by the number of GPUs, assuming no parallelization losses. Although such losses have been significantly reduced in recent years, they are not entirely negligible, making some of the reported figures conservative estimates. ^b Some inconsistencies have been observed in the reported performance of their simulations. However, the performance values remain within a credible range.