Evaluation of the lattice Boltzmann method for wind modelling in complex terrain

The worldwide expansion of wind energy is making the choice of potential wind farm locations more and more difficult. This results in an increased number of wind farms being located in complex terrain, which is characterised by flow separation, turbulence and high shear. Accurate modelling of these flow features is key for wind resource assessment in the planning phase, as the exact positioning of the wind turbines has a large effect on their energy production and lifetime. Wind modelling for wind resource assessments is usually carried out with the linear model Wind Atlas Analysis and Application Program (WAsP), unless the terrain is complex, in which case Reynolds-averaged Navier–Stokes (RANS) solvers such as WindSim and Ansys Fluent are usually applied. Recent research has shown the potential advantages of large-eddy simulation (LES) for modelling the atmospheric boundary layer and thermal effects; however, LES is far too computationally expensive to be applied outside the research environment. Another promising approach is the lattice Boltzmann method (LBM), a computational fluid technique based on the Boltzmann transport equation. It is generally used to study complex phenomena such as turbulence, because it describes motion at the mesoscopic level in contrast to the macroscopic level of conventional computational fluid dynamics (CFD) approaches, which solve the Navier–Stokes (N–S) equations. Other advantages of the LBM include its efficiency; near-ideal scalability on high-performance computers (HPCs); and ability to easily automate the geometry, the mesh generation and the post-processing. However, the LBM has been applied very little to wind modelling in complex terrain for wind energy applications, mainly due to the lack of availability of easy-to-use tools as well as the lack of experience with this technique. In this paper, the capabilities of the LBM to model wind flow around complex terrain are investigated using the Palabos framework and data from a measurement campaign from the Bolund Hill experiment in Denmark. Detached-eddy simulation (DES) and LES in Ansys Fluent are used as a numerical comparison. The results show that there is in general a good agreement between simulation and experimental data, and the LBM performs better than RANS and DES. Some deviations can be observed near the ground, close to the top of the cliff and on the lee side of the hill. The computational costs of the three techniques are compared, and it has been shown that the LBM can perform up to 5 times faster than DES, even though the set-up was not optimised in this initial study. It can be summarised that the LBM has a very high potential for modelling wind flow over complex terrain accurately and at relatively low costs, compared to solving N–S equations conventionally. Further studies on other sites are ongoing.

1 Introduction

In order to assess the wind resource for both the planning and the assessment of wind farms, measurements and simulations of the prevailing wind conditions are required. Simulations are especially crucial in the observation of flows over complex terrain due to the large impact of steep inclines on the flow conditions. If the terrain shows only weak topographic changes or low hills, linear models can be used to make fast and sufficiently accurate yield forecasts (Berg and Kelly, 2019). The extremely low computational effort and ease of use make such models the current industry standard. Due to their simplified formulation, however, such models fail in complex terrain and the predictions can be unreliable (Bowen and Mortensen, 1996). For complex flows, non-linear methods that solve the Navier–Stokes (N–S) equations are better suited. The successful use of Reynolds-averaged Navier–Stokes (RANS) models has been demonstrated in many studies (e.g. Ferreira et al., 1995; Maurizi et al., 1998; Kim et al., 2000; Castro et al., 2003), and they are being used increasingly in the industry. This is reflected by the recent development of wind-energy-specific tools using RANS-based computational fluid dynamics (CFD), including WAsP CFD (Bechmann, 2012) and WindSim (Dhunny et al., 2016). The RANS equations govern the transport of the averaged flow quantities, with the whole range of the scales of turbulence being modelled using turbulence closure schemes. The RANS-based modelling approach therefore greatly reduces the required computational effort and resources compared to fully resolved methods and is widely adopted for practical engineering applications. A more detailed modelling of turbulence is possible using large-eddy simulation (LES). LES lies between direct numerical simulation (DNS) and turbulence closure schemes. The idea of this method is to compute the mean flow and the large vortices exactly. The small-scale structures are not simulated, but their influence on the rest of the flow field is parameterised by a heuristic model. However, the computational effort and the demands on the computational mesh increase drastically compared to RANS simulations, due to the need to resolve the small and important dynamic eddies in the boundary layer. Recent studies of the Bolund Hill blind test also show that it is still a great challenge to achieve sufficiently accurate predictions using LES (Bechmann et al., 2011; Diebold et al., 2013; Ma and Liu, 2017; DeLeon et al., 2018). This is because, to accurately resolve the small-scale turbulent structures near walls at high Reynolds numbers, an extremely fine grid resolution is required.

The detached-eddy-simulation (DES) method is a combination of LES and RANS. With this method, the flow is mostly calculated by LES, but the flow and vortices in wall regions are modelled by RANS. This method promises a strong reduction in the computational effort and the mesh requirements compared to LES. In addition, boundary layer modelling using RANS models makes it possible to use surface roughness models (Bechmann and Sørensen, 2010).

An alternative to solving the N–S equations with great potential is the lattice Boltzmann method (LBM). The LBM has become more and more popular in recent years and is being continuously developed further. The LBM has also been used successfully for initial studies in the field of wind energy. Most of these studies focus on the simulation of flows around wind turbines and wind farms or analyse the wake behaviour of turbines (e.g. Deiterding and Wood, 2016; Asmuth et al., 2019). Studies have shown that the LBM is a valid alternative to conventional CFD methods and has many advantages (Wang et al., 2018). The main advantage of the method is its almost ideal scalability. This makes the application on high-performance computers (HPCs) attractive, but LBM codes based on graphics processing units (GPUs) have also been implemented recently (Schönherr et al., 2011; Onodera and Idomura, 2018). This makes it possible to perform computationally intensive LESs on a simple desktop in a reasonable time (Asmuth et al., 2019). However, the LBM has not yet been assessed a great deal for the calculation of wind fields in complex terrain for wind energy applications.

The goal of this present paper is therefore to evaluate the capabilities of the LBM for wind modelling in complex terrain. Ansys Fluent is used as a reference for comparisons, using both a RANS and a DES approach. The paper starts with a brief introduction of the theories behind the LBM and the conventional N–S-based CFD calculations in Sect. 2, then introduces the simulation method applied in Sect. 3, discusses the results in Sect. 4, and finishes with the conclusions in Sect. 5.

2 Lattice Boltzmann method theory

2.1 Numerical method and governing equations

Interest in the LBM has been growing in the past decades as an efficient method for computing various fluid flows, ranging from low-Reynolds-number flows to highly turbulent flows (e.g. Chen and Doolen, 1998; Filippova et al., 2001). The first LBM models struggled with high-Reynolds-number flows due to numerical instabilities. To solve this problem, various adaptations such as regularised finite difference (Latt and Chopard, 2006), multiple relaxation time (d'Humieres, 2002) or entropic methods (Ansumali and Karlin, 2000) have been developed.

The LBM has the following advantages over N–S: (1) it is a linear equation with only local non-linearity, making it more stable and perfectly scalable; (2) the dissipation is introduced locally by the collision term and does not depend on the lattice; and (3) the relaxation time includes both the regular viscous effects and its higher-order modifications. A description of the LBM can be found, for example in Chen and Doolen (1998). The governing equations describe the evolution of the probability density of finding a set of particles with a given microscopic velocity at a given location:

$$\begin{array}{}\text{(1)}& f_i\left(\mathit{x}+{\mathit{c}}_i\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)=f_i(\mathit{x},t)+{\mathrm{\Omega }}_i(\mathit{x},t)\end{array}$$

for $\mathrm{0}\le i<q$, where c_i represents a discrete set of q velocities, f_i(x, t) is the discrete single-particle distribution function corresponding to c_i and Ω_i is an operator representing the internal collisions of pairs of particles. Macroscopic values such as density ρ and the flow velocity u can be deduced from the set of probability density functions f_i(x, t), such as

$$\begin{array}{}\text{(2)}& \mathit{\rho }=\sum _{i=\mathrm{0}}^{q-\mathrm{1}}{f}_i,\mathit{\rho }\mathit{u}=\sum _{i=\mathrm{0}}^{q-\mathrm{1}}{f}_i{\mathit{c}}_i.\end{array}$$

The set of allowed velocities in the LBM is restricted by conservation of mass and momentum and by rotational symmetry (isotropy). Some of the most popular choices for the set of velocities are D2Q9 and D3Q27 lattices, referring to nine velocities in 2D and 27 velocities in 3D. For both of these lattices, the speed of sound in lattice units is given by $c_{\mathrm{s}}=\mathrm{1}/\sqrt{\mathrm{3}}$ (Succi, 2001). The collision operator Ω_i is typically modelled with the Bhatnagar–Gross–Krook (BGK) approximation (Bhatnagar et al., 1954). It assumes that the fluid locally relaxes to equilibrium over a characteristic, non-dimensional timescale τ. The relaxation time τ determines how fast the fluid approaches equilibrium and is thus directly dependent on the viscosity of the fluid. The corresponding equilibrium probability density function $f_i^{\left(\mathrm{eq}\right)}$ is defined as

$$\begin{array}{}\text{(3)}& {\mathrm{\Omega }}_i=-{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[f_i\left(\mathit{x},t-f_i^{\left(\mathrm{eq}\right)}(\mathit{x},t)\right)\right].\end{array}$$

The equilibrium distribution function $f_i^{\left(\mathrm{eq}\right)}$ is a local function that only depends on density and velocity in the isothermal case. It can be computed thanks to a second-order development of the Maxwell–Boltzmann equilibrium function (Qian, 1992):

$$\begin{array}{}\text{(4)}& f_i^{\left(\mathrm{eq}\right)}=w_i\mathit{\rho }\left[\mathrm{1}+{\displaystyle \frac{{\mathit{c}}_i\cdot \mathit{u}}{c_{\mathrm{s}}^{\mathrm{2}}}}+{\left({\displaystyle \frac{{\mathit{c}}_i\cdot \mathit{u}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}}\right)}^{\mathrm{2}}-{\displaystyle \frac{{\mathit{u}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{s}}^{\mathrm{2}}}}\right],\end{array}$$

where w_i refers to the Gaussian weights of the lattice. A Chapman–Enskog expansion, based on the assumption that f_i is given by the sum of the equilibrium distribution plus a small perturbation $f_i^{\mathrm{1}}$,

$$\begin{array}{}\text{(5)}& f_i=f_i^{\left(\mathrm{eq}\right)}+f_i^{\left(\mathrm{1}\right)},\mathrm{with}{f}_i^{\left(\mathrm{1}\right)}\ll f_i^{\left(\mathrm{eq}\right)}\end{array}$$

can be applied to Eq. (1) in order to recover the exact N–S equation for quasi-incompressible flows in the limit of long wavelengths (Chapman et al., 1990). The lattice pressure is thus given by $p=c_{\mathrm{s}}^{\mathrm{2}}\mathit{\rho }$, and the lattice viscosity is linked to the BGK relaxation parameter through

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

The numerical scheme is divided into two steps:

• a collision step where the BGK model is applied,

  $$\begin{array}{}\text{(7)}& f_i\left(\mathit{x},t+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)=f_i(\mathit{x},t)+{\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[f_i^{\left(\mathrm{eq}\right)}(\mathit{x},t)-f_i(\mathit{x},t)\right];\end{array}$$

• a streaming step,

  $$\begin{array}{}\text{(8)}& f_i\left(\mathit{x}+{\mathit{c}}_i,t+\mathrm{1}\right)=f_i\left(\mathit{x},t+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right).\end{array}$$

In the collision step, particle populations interact and change their velocity directions according to scattering rules. This operation is completely local which makes the LBM well suited for parallelism. The streaming step consists of an advection of each discrete population to the neighbour node located in the direction of the corresponding discrete velocity. Since a boundary node has fewer neighbours than an internal node, some populations are missing at the boundary after each iteration. These populations need to be reconstructed, which is the purpose of the implementation of boundary conditions in the LBM.

2.2 Turbulence modelling

Turbulence leads to the appearance of eddies with a wide range of length scales and timescales, which interact with each other in a dynamically complex way. Given the importance of the avoidance or promotion of turbulence in engineering applications, it is no surprise that a substantial amount of research effort is dedicated to the development of numerical methods to capture the important effects due to turbulence. The methods can be grouped into the following four categories:

• turbulence models for Reynolds-averaged Navier–Stokes (RANS) equations,

• large-eddy simulation (LES),

• detached-eddy simulation (DES),

• direct numerical simulation (DNS).

In this work, LES was applied for the LBM simulations. LES is an intermediate form of turbulence calculation which simulates the behaviour of the larger eddies. The method involves spacial filtering, which passes the larger eddies and rejects the smaller eddies. The effects on the resolved flow (mean flow plus large eddies) due to the smallest, unresolved eddies are included by means of a so-called sub-grid-scale model. It is assumed that the sub-grid scales have the effect of a viscosity correction, which is proportional to the norm of the strain-rate tensor at the level of the filtered scales: $\mathit{\nu }={\mathit{\nu }}_{\mathrm{0}}+{\mathit{\nu }}_{\mathrm{T}}$. ν_T is defined as

$$\begin{array}{}\text{(9)}& {\mathit{\nu }}_{\mathrm{T}}=\left(C\mathrm{\Delta }{)}^{\mathrm{2}}\right|S|,\end{array}$$

where C is the Smagorinsky constant, Δ is the grid size and the tensor norm of the strain rate is defined as $\left|S\right|=\sqrt{S:S}$. The value of the Smagorinsky constant depends on the physics of the problem and usually varies between 0.1 and 0.2 far from boundaries (Davidson, 2015).

3 Simulations

3.1 The Bolund Hill experiment

The Bolund field campaign took place from December 2008 to February 2009 on Bolund Hill in Denmark. Bolund Hill is a 130 m long (east–west axis), 75 m wide (north–south axis) and 11.7 m high hill, situated near the Risø Campus of the Technical University of Denmark. Details of the experiment are described in Bechmann et al. (2011). The campaign showed dominant wind directions from the west and south-west. Thus the wind has an extensive upwind fetch over the sea before encountering land, leading to a “steady” flow on the windward side of the hill. The wind first encounters a 10 m vertical cliff, after which air flows back down to sea level on the east side of the hill. A recirculation zone and a flow separation are expected due to this abrupt change in slope. During the campaign, 35 anemometers were deployed over the hill. The location of the measurement devices can be seen in Fig. 1. Instrumentation included 23 sonic anemometers, 12 cup anemometers and 2 lidars. At each measurement location, the three components of the wind velocity vector and their variances were recorded for four different dominant wind directions, three westerly winds originating from the sea (268, 254 and 242^∘) and one easterly wind originating from the land (95^∘). The mean wind speed during the measurements was around 10 m s⁻¹, leading to a Reynolds number of $Re=uh/\mathit{\nu }\approx {\mathrm{10}}^{\mathrm{7}}$ with the free-stream velocity u=10 m s⁻¹, the hill height h and the kinematic viscosity ν. The measured values are 10 min averages of measurements sampled at 20 Hz for sonic anemometers.

Figure 1 A contour map of Bolund Hill with meteorological masts denoted from M0 to M9 (Bechmann et al., 2011).

3.2 Simulation set-up

3.2.1 Boundary conditions

Palabos

The LBM flow solver used in this work was the Palabos open-source library (Latt et al., 2009). The Palabos library is a framework for general-purpose CFD with a kernel based on the LBM. The use of C$++$ code makes it easy for experienced programmers to install and run the code on any machine. It is thus possible for experienced modellers to set up fluid flow simulations with relative ease and to extend the open-source library with new methods and models, which is of paramount importance for the implementation of new boundary conditions.

To calculate the wind fields with Palabos in this work a 525 m long (east–west axis), 250 m wide (north–south axis) and 40 m high domain with a uniform grid resolution of $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=\mathrm{0.5}$ m was used, leading to a total cell count of 46 million. Palabos has the capability to include multiple grid sizes (octree grid structure) to refine the grid near the hill and coarsen it in the far field; however these techniques were not in the scope of this first study.

There are no turbulence closure models or surface roughness models implemented in the Palabos library yet; therefore the water surfaces were prescribed as free-slip bounce-back nodes and the ground surfaces were modelled using regularised bounce-back nodes (Malaspinas et al., 2011; Izham et al., 2011). The bounce-back scheme in this first study was chosen due to its simple implementation and robustness. There are more sophisticated models, like the immersed boundary method, which may provide better accuracy than the staircase approximation of bounce-back nodes, which will be investigated in further studies.

The inlet profile was described according to the Bolund Hill blind test specification for the westerly wind case. The logarithmic velocity profile is defined as

$$\begin{array}{}\text{(10)}& u\left(z_{\mathrm{agl}}\right)={\displaystyle \frac{u_{*\mathrm{0}}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{z_{\mathrm{agl}}}{z_{\mathrm{0}}}}\right),\end{array}$$

with κ=0.4, the friction velocity $u_{*\mathrm{0}}=\mathrm{0.4}$, the elevation above ground level $z_{\mathrm{agl}}=z-gl$ (gl=0.75 m) and the roughness length z₀=0.0003 m. Additionally, a time-varying fluctuation in the wind speed, corresponding to the given turbulence intensity value, was superposed. The logarithmic wind profile was updated every second during the simulation. The atmospheric boundary layer (ABL) was considered neutral, and thermal effects are therefore neglected. Both the sides and the top of the domain were modelled as free-slip walls (zero normal velocity). The outlet was set to a constant pressure. Each simulation was run for 600 s with a time step Δt=2.89 ms, leading to around 10 advection times. The relaxation time τ and the viscosity ν were chosen to respect Eq. (6) and to stabilise the solution. The Smagorinsky constant was set to 0.14.

Fluent

Ansys Fluent contains the broad, physical-modelling capabilities needed to model flow, turbulence, heat transfer and reactions for industrial applications, ranging from airflow over an aircraft wing to combustion in a furnace, from bubble columns to oil platforms, from blood flow to semiconductor manufacturing and from clean room design to wastewater treatment plants. For the Fluent simulations in this work the mesh was created with the new improved Fluent meshing tool; additionally the domain was extended to 830 m × 450 m × 60 m, and two mesh refinement zones near the hill were implemented. The mesh resolution ranged from 0.5 m near the hill up to 15 m in the far field, resulting in a total cell count of 10 million. A roughness length of z₀=0.3 mm was prescribed for the water surface and a roughness length of z₀=15 mm for the ground surfaces. The RANS simulation, using the k–ω SST turbulence model, was used to initialise the flow and turbulence quantities for the DES. Each simulation was run for 600 s with a time step Δt of 50 ms, leading to around seven advection times for the DES Fluent simulations. The inlet velocity was described as discussed before. According to the blind test, the turbulent kinetic energy (TKE) at the inlet was set to 0.928 m² s⁻². For the DES model a synthetic turbulence generator scheme was used to generate a time-dependent inlet condition. It uses a Fourier-based synthetic turbulence generator. This method is inexpensive in terms of computational time compared with the other existing methods while achieving high-quality turbulence fluctuations (ANSYS, 2019). A free-slip boundary condition is used for all the flow variables at the top and the side boundaries. The Smagorinsky constant was set to 0.14.

4 Results and discussion

4.1 Flow comparisons

The calculated velocity magnitude fields at a vertical plane through the position of met mast M3 for each measurement technique are shown in Figs. 2 and 3. The LBM results are shown in Fig. 2, in terms of the averaged velocity magnitude over the simulation time (Fig. 2a) and the instantaneous velocity magnitude at time t=600 s (Fig. 2b). Figure 3a and b show the averaged velocity magnitude over the simulation time for RANS simulations and DESs, and Fig. 3c shows the instantaneous velocity magnitude at time t=600 s for DESs. It is interesting to note the separation region as the wind flows over the sharp edge of the hill, as well as the highly separated flow at its rear side.

Figure 2 Velocity field over the hill along the B line (LBM results): (a) averaged velocity field, (b) instantaneous velocity field at t=600 s.

Figure 3 Velocity field over the hill along the B line (Fluent results): (a) averaged velocity field (RANS), (b) averaged velocity field (DES), (c) instantaneous velocity field (DES).

4.2 Performance comparisons

For a quantitative comparison, the same methodology is used as described by Bechmann et al. (2011) for the wind flow along the 270^∘ axis (Case 1) as shown in Fig. 1. This involved investigating the difference between measurements and simulations after the mast M0 by comparing and quantifying the changes in the wind field as both changes in speed (so-called “speed-up”) and in direction (so-called “turning”). Speed-up is defined as

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }{S}_{\mathrm{m}}={\displaystyle \frac{\langle \stackrel{\mathrm{‾}}{s}{\rangle }_{z_{\mathrm{agl}}}-\langle \stackrel{\mathrm{‾}}{s_{\mathrm{0}}}{\rangle }_{z_{\mathrm{agl}}}}{\langle \stackrel{\mathrm{‾}}{s_{\mathrm{0}}}{\rangle }_{z_{\mathrm{agl}}}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{s}$ is the mean wind speed at the sensor location and $\stackrel{\mathrm{‾}}{s_{\mathrm{0}}}$ is the mean wind speed at the mast M0. Turning is defined as the difference between the wind direction at the measurement point and that at M0. The comparison is made in Fig. 5 for two different elevations, 2 and 5 m a.g.l. (above ground level), and for the four masts along the B line (M7, M6, M3 and M8). The simulation results for the speed-up (Fig. 4) shows good agreement with experimental data for all simulation techniques at 5 m a.g.l., with all deviations lower than 7.1 % and the average speed-up error for each simulation technique shown in Table 1. The average speed-up error is defined as

$$\begin{array}{}\text{(12)}& R_{\mathrm{s}}=\mathrm{100}\left(\mathrm{\Delta }{S}_{\mathrm{s}}-\mathrm{\Delta }{S}_{\mathrm{m}}\right),\end{array}$$

where S_m is the measured speed-up and S_s is the simulated speed-up defined by Eq. (11). Table 1 also allows the three simulation techniques to be compared to each other. The 2 m a.g.l. results show higher deviations in general, with the average speed-up error for each simulation technique shown in Table 1. The highest discrepancy can be seen at M6, which is probably due to the separation bubble observed in the velocity fields in Fig. 2a. The experiment showed reduction in wind speed at M6, whereas the simulations all show an increase in wind speed. This leads to the conclusion that the actual separation bubble is larger than the simulated one. This could be due to an error in the 3D-model capture of the overhang of the hill noted in previous studies (Lange et al., 2017). Furthermore, all the simulation techniques underpredicted the negative speed-up in the highly separated region of M8 compared to the experiment. The reason for this is probably due to the well-known difficulty of correctly simulating the separation point in CFD. As this effect is particularly pronounced at a height of 2 m a.g.l., it may be due to the fact that the lower measuring points lie within the boundary layer and the used models were not able to capture the near-wall flow entirely correctly, perhaps due to the assumptions regarding surface roughness.

As shown in Table 1, the most accurate overall prediction was the LBM simulation, with an average error over all the measurement positions of 8.0 %. The RANS and DES mean errors are 10.0 % and 17.3 %, respectively. All three methods showed significantly more accurate results at 5 m a.g.l. than at 2 m a.g.l.

Figure 4 Speed-up along Bolund Hill. Wind direction is from 270^∘.

Figure 5 Turning along Bolund Hill. Wind direction is from 270^∘.

Figure 6 Turbulence kinetic energy (TKE) along Bolund Hill. Wind direction is from 270^∘.

Table 1 Average speed-up error.

For the turning of the wind, a similar behaviour can be observed. The results match the experimental data very well at 5 m a.g.l., with all deviations lower than 3.0 % and the average turning error for each simulation technique shown in Table 2. As for the speed-up, the deviations in turning are higher at 2 m a.g.l., with the average turning error for each simulation technique shown in Table 2. The highest discrepancy can be seen at M8. Met mast M8 is located at the lee side in the recirculation zone of the hill. All the simulation results struggle to capture the flow accurately in terms of the turning. This could be due to the inaccuracy in predicting the exact separation location on the rear of the hill, as mentioned above.

Figure 7 Scatter plot of wind speeds, measurement against simulation results.

Table 2 Average turning error.

Further analysis using the entire set of measurement data is shown in Fig. 7, in which a comparison between the simulation and experimental data for all three simulation methods is shown. Overall there is a good agreement between the measurements and simulated results. M2 and M6, both right after the edge of the cliff, show the biggest mismatch due to the detached flow after the edge of the hill, as discussed above. The next two figures show the ratio of simulated wind speeds to measured wind speeds as a function of elevation (Fig. 8) and measurement location (Fig. 9). The biggest deviation between the data can again be seen at lower heights and at masts M2, M6 and M8. Between the simulation methods, the LBM shows the highest averaged deviation of the ratios. The DES and RANS model both perform better in this comparison. This may be due to both these models using the k–ω SST turbulence model and incorporating the surface roughness to calculate the near-wall turbulence. The reason for the DES model performing worse than the RANS model is unclear at this point and requires further investigation.

Figure 8 Ratio of simulation results to experimental wind speeds as a function of elevation. The dotted grey line represents the average value.

Figure 9 Ratio of simulation results to experimental wind speeds as a function of measurement location. The dotted grey line represents the average value.

Finally Fig. 6 shows the simulated turbulence kinetic energy (TKE) compared to the measurements for met masts M7, M6, M3 and M8. Overall there is not a particularly good agreement between the measured and simulated data. All the simulations show similar but lower values for TKE. Especially at M6, at the cliff, the deviation is the highest. This discrepancy is being further investigated, but several authors have reported difficulties in simulating a horizontally homogeneous ABL flow in at least the upstream part of computational domains (Blocken et al., 2007; Zhang, 1994).

Figure 10 Comparison of computational time per CPU core and million cells.

4.3 Cost comparisons

In this section, the performance of the simulation techniques is compared in terms of the computational costs. This has been done because the overall cost of a simulation is an important factor for modellers, who need to choose the most suitable model for a given wind energy project. The results of this work have been used in order to develop a new method for helping wind modellers choose the most cost-effective model for a given project. This was done by first defining various parameters for predicting the skill and cost scores before carrying out the simulations as well as for calculating skill and cost scores after carrying out the simulations. Weightings were then defined for these parameters, and values assigned to them for a range of tools, including the ones applied in the present work, using a template containing predefined limits in a blind test. This allowed for a graph of predicted skill score against cost score to be produced, enabling modellers to choose the most cost-effective model without having to carry out the simulations beforehand. More details can be found in Barber (2020).

Figure 10 and Table 3 summarise the computational costs for the three different techniques applied in this paper. It can clearly be seen that the LBM performed 5 times faster than the DES and only slightly slower than the steady RANS simulation. This is due to its explicit formulation and exact advection operator. Furthermore, each of the collision and streaming processes is independent at each lattice, which makes the method so suitable for parallelisation. This advantage extends also to other types of high-performance hardware like graphics processing units (GPUs). Some studies of GPU-based LBM solvers show promising results in this field (Asmuth et al., 2019; Schönherr et al., 2011; Onodera and Idomura, 2018). The performance of this LBM simulation could be increased by adapting the code to use different grid sizes, depending on the flow, therefore reducing the overall cell count drastically. Incorporating the same grid refinement zones as used in the Fluent simulation, while maintaining the extended domain zone, the cell count for the resulting grid would decrease by a factor of 5 to 7. Work on this is ongoing.

Table 3 Computational time. All simulation were run on 80 cores (Intel Xeon E5-2630v4 2.2 GHz).

^* Computational time per CPU core and million cells.

5 Conclusion

In this study, LES using the LBM framework Palabos was implemented to calculate the wind field over the complex terrain of Bolund Hill. Advantages of the LBM include its efficiency; near-ideal scalability on high-performance computers (HPCs); and capabilities to easily automate the geometry, the mesh generation and the post-processing.

The results were compared to RANS simulations and DESs using Ansys Fluent and field measurements. In general there was a good agreement between simulation and experimental data. The average wind speed-up error compared to measurements was 8.0 % for the LBM, 17.3 % for DES and 10.0 % for RANS. The average wind turning error compared to measurements was 2.7^∘ for the LBM, 2.0^∘ for DES and 0.2^∘ for RANS. Some deviations could be observed near the ground, close to the top of the cliff (M2) and on the lee side of the hill (M8). Larger deviations could be observed for the TKE calculation, especially at met mast M6, which is positioned right after the edge of the cliff. This corresponds to previous work, which shows difficulties in correctly resolving the TKE.

The computational costs of these three models were compared, and it has been shown that the LBM, even in this set-up of the simulation which has not yet been fully optimised, can perform 5 times faster than DES and lead to slightly more accurate results.

It can be summarised that the LBM may be applicable to modelling wind flow over complex terrain accurately at relatively low costs if the challenges raised in this work are addressed. Further studies on other sites are ongoing.

Appendix A: Non-dimensioning procedure

The non-dimensioning procedure used in this study is carried out according to the similarity theory. It consists of two steps. First a physical system is converted into a dimensionless system, independent of the original physical scales but also independent of simulation parameters. In a second step, the dimensionless system is converted into a discrete simulation. Thus the dimensionless level (D) links the physical system (P) with the discrete lattice Boltzmann system (LB). The solutions to the incompressible Navier–Stokes equations for example depend only on the Reynolds number (Re). Thus, the three systems are defined to have the same Reynolds number. The transition from P to D is made through the choice of a characteristic length scale l₀ and timescale t₀, and the transition from D to LB is made through the choice of a discrete space step Δx and time step Δt (Latt, 2008).