Open AccessArticle

Highly Stable Lattice Boltzmann Method with a 2-D Actuator Line Model for Vertical Axis Wind Turbines

Luca Cacciali

Martin O. L. Hansen

^1,*

and

Krzysztof Rogowski

DTU Wind and Energy Systems, Fredensborgvej 399, 4000 Roskilde, Denmark

Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology, 00-665 Warsaw, Poland

Author to whom correspondence should be addressed.

Energies 2024, 17(19), 4847; https://doi.org/10.3390/en17194847

Submission received: 6 September 2024 / Revised: 19 September 2024 / Accepted: 23 September 2024 / Published: 27 September 2024

(This article belongs to the Topic Advances in Wind Energy Technology)

Download

Browse Figures

Figure 1
D2Q9 velocity set for a 2-D fluid-flow problem. "> Figure 2
R500 case study [<a href="#B31-energies-17-04847" class="html-bibr">31</a>]. "> Figure 3
Computational domain. "> Figure 4
Sensitivity of the relative wind speed Vrel to the increasing mesh resolution Δh for simulations with ε equal to 0.5c. "> Figure 5
Numerical model and mesh with boundary conditions. "> Figure 6
Effect of the corrections to the normal force coefficient Cn* for ε = 0.5c (left-hand side) and ε = c (right-hand side). BASELINE = uncorrected curve; FC = flow curvature; DS = dynamic stall; T = Smagorinsky turbulence model. "> Figure 7
Effect of the corrections to the tangential force coefficient Ct* for ε = 0.5c (left-hand side) and ε = c (right-hand side). "> Figure 8
Effect of the corrections to the angle of attack α for ε = 0.5c (left-hand side) and ε = c (right-hand side). "> Figure 9
Blade normal force coefficient Cn* vs. azimuthal angle ϑ (left-hand side) and tangential force coefficient Ct* (right-hand side) at a high tip-speed ratio (λ = 4.5). Comparison between PIV experimental data (EXP), Lattice Boltzmann Method (LBM), Scale-Adaptive Simulation (SAS), and Vortex Method (VM). "> Figure 10
Blade normal force coefficient Cn* vs. azimuthal angle ϑ (left-hand side) and tangential force coefficient Ct* (right-hand side) at a low tip-speed ratio (λ = 2.0). Comparison between PIV experimental data (EXP), Lattice Boltzmann Method (LBM), Scale-Adaptive Simulation (SAS), and Vortex Method (VM). "> Figure 11
Axial induction factor aϑ vs. azimuthal angle ϑ (left-hand side) and angle of attack α vs. azimuthal angle ϑ (right-hand side) for tip-speed ratios λ of 4.5 and 2.0. "> Figure 12
Normalized streamwise velocity ux,n at different stations x/D downstream of the turbine for λ = 4.5 (left-hand side) and λ = 2.0 (right-hand side). "> Figure 13
Contours of normalized out-of-plane vorticity for the combined field ωzc/u0 for λ = 4.5 (top) and λ = 2.0 (bottom). "> Figure 13 Cont.
Contours of normalized out-of-plane vorticity for the combined field ωzc/u0 for λ = 4.5 (top) and λ = 2.0 (bottom). ">

Versions Notes

Abstract

A 2-D Lattice Boltzmann Method, designed to ensure stability at high Reynolds numbers, is combined with an Actuator Line Model to compute the loads on a two-bladed vertical axis wind turbine. Tests on the kernel size at a high mesh resolution reveal that a size equal to half of the full chord length yields the most accurate results. The aerodynamic load solution is validated against a fully resolved Scale-Adaptive Simulation (SAS) output, demonstrating high correlation, and enabling an assessment of near wake and downstream effects. The model’s adaptability to various rotor operating conditions is confirmed through tests at high and low tip-speed ratios. Additionally, a Biot–Savart-based Vortex Model (VM) is employed for further comparison, showing good agreement with the Lattice Boltzmann output. The results indicate that the Highly Stable Lattice Boltzmann Method integrated with the Actuator Line Model enhances the accuracy of flow field resolution and effectively captures complex aerodynamic phenomena, making it a valuable tool for simulating vertical axis wind turbines.

Keywords:

lattice Boltzmann method; LBM; VAWT; SAS; CFD; vortex method

1. Introduction

Different numerical methods of differing sophistication exist to calculate the performance of a Vertical Axis Wind Turbine (VAWT), as thoroughly detailed in [1]. These methods range from momentum-based methods such as the Double Multiple Streamtube model (DMS) [2,3], Vortex Method (VM) [1], Actuator Line Method (ALM) [4], and Computational Fluid Dynamics (CFD) with fully resolved boundary layer on the blade surfaces [5,6,7,8]. Integrating ALM with CFD reduces the computational cost of a simulation compared to a fully resolved CFD simulation. However, the accuracy of the approximation relies significantly on the choice of the kernel size parameter, which is crucial for achieving a proper distribution of forces around the actuator line [9]. Sørensen and Shen [4] pioneered the ALM for a Horizontal Axis Wind Turbine (HAWT), avoiding resolving the boundary layer. Shamsoddin and Porté-Agel [10] investigated the wake of a Vertical Axis Wind Turbine (VAWT) in the boundary layer using large eddy simulation (LES) combined with an ALM. They examined variations in the power coefficient for a 1-MW VAWT, identifying the optimal solidity and tip-speed ratio. Their wake analysis revealed maximum velocity deficits and turbulence intensities at specific downstream distances, with significant vertical asymmetries in the turbulence intensity and momentum flux fields. Their research also demonstrated the effectiveness of coupling ALM with LES for advanced applications, such as analyzing both near and far wakes [11].

Bachant et al. [12] applied the ALM model to VAWTs. They observed that RANS ALM tends to overpredict power coefficients at high tip-speed ratios due to unmodeled factors such as drag from spokes and limitations in aerodynamic data. Martinez-Ojeda et al. [13], aiming to optimize VAWT wind farms, proposed a hybrid wake actuator–cylinder model within an OpenFOAM v2012 solver. This model replaces turbines with volumetric forces, reducing the computational cost required for detailed meshes around the blades. Using the steady-state actuator–cylinder model [14], this two-dimensional approach requires minimal computational time and offers higher fidelity than simpler models by resolving the viscous wake through RANS simulation. Mendoza and Goude [15] investigated the wake produced by a VAWT and its interaction with terrain surface roughness using an ALM implemented in OpenFOAM. Their simulations, validated against experimental data, focused on different terrains and similar inflow velocities. The results demonstrated the model’s ability to reproduce the wake pattern, capturing key flow parameters, vorticity concentrations, and turbulence production.

From the perspective of integrating ALM with a fluid flow solver, the lattice Boltzmann method (LBM) offers an alternative to traditional CFD, standing out as a mesoscale approach among particle-based solvers for addressing complex fluid dynamics problems. The lattice Boltzmann equation (LBE) translates microscopic particle interactions into observable macroscopic fluid behavior based on the kinetic theory of gases. The LBE describes the evolution of a probability distribution function f(x,c,t) at time t over a discrete lattice in space between x and x + dx, with mesoscale particle speed between c and c + dc. Due to the direct link between f and the macroscopic fluid density and velocity, the value of f is determined directly at the lattice nodes. Through collision and streaming, f develops, resulting in macroscopic fluid properties such as velocity, density, and pressure [9].

The LBM offers several advantages over traditional CFD, including an easier implementation than conventional methods. Because the most computationally intensive operations are confined to individual nodes, LBM is well-suited for high-performance computing on parallel architectures [16], reducing the run time. LBM is useful in simulating mass-conserving flows [16,17,18] and moving boundaries [17]. Another advantage of LBM is that it avoids the need for upwind schemes, typically used in advection-dominated CFD methods. Instead, LBM is based on mesoscopic fluid distribution, which implies that particles are streamed along lattice directions to simulate fluid flows, relying on algorithms for collision (relaxation) and streaming.

A disadvantage of LBM is its higher memory requirement than CFD [19]. The numerous variables propagating to and from each node make LBM a memory-intensive method. Because LBM is time-dependent, it is more suitable for simulating transient rather than steady flows [20].

The LBM has been used for wake studies of Horizontal Axis Wind Turbines by the University of Uppsala [21,22] applying the ELBE code [23]. Asmuth et al. [22] achieved similar wake statistics by applying laminar and turbulent inflow to their ALM assisted by the LBM solver. LBM was also used for calculating the flow over complex terrain [24] applying the PALABOS code [25]. Ribeiro and Muscari [26] simulated wind farms with an LBM-ALM combined with a sliding mesh technique, validating their results with Navier–Stokes solutions. The sliding mesh technique allows for targeted refinement around the blades and tip vortices, improving computational efficiency and providing significant advantages over traditional methods.

Finally, Rullaud et al. [27] developed a solver for VAWTs based on the LBM-ALM coupling. Their model was validated with experimental data and outputs of other simulation methods. The LBM solver of Rullaud’s work was developed with an MRT collision operator [28] and a forcing term modeled according to Kupershtokh’s equations [29]. Although showing promising results compared to experiments, the MRT collision operator requires several numerical attempts to set shear and bulk viscosity factors and stabilize the factors related to higher-order moments to comply with applications at high Reynolds numbers. This operation is not straightforward, and is more time-consuming in fine grids [19].

The Simplified and Highly Stable Lattice Boltzmann Method (SHS-LBM) [30] is proposed to achieve more accurate and efficient VAWT performance simulations while enhancing the numerical approach’s robustness. Although the computational time in SHS-LBM is higher than Conventional LBM [30], no additional parameters should be tuned in SHS-LBM [9].

The ALM is used to compute the force field generated from the lift and drag forces, but instead of directly solving the Navier–Stokes equations, an LBM is employed to calculate the flow past the VAWT in time.

The two-bladed reference turbine illustrated in [31,32] is numerically investigated using the integrated ALM-LBM for validation. Static airfoil data are interpolated using a table lookup procedure and adjusted for dynamic stall and flow curvature. The analysis is further simplified by using the polar database of the symmetric airfoil NACA 0018, which removes the additional variables associated with camber. The behavior of the reference VAWT is explored under stable (high tip-speed ratio) and unstable (low tip-speed ratio) operating conditions, validating the outputs with those achieved from a traditional ALM-CFD and a VM.

With these simulations, this paper aims to demonstrate the attractiveness of the LBM as a reliable alternative for applying the Actuator Line Model to vertical axis wind turbines. Based on the SHS-LBM technique, this model differs significantly from Rullaud’s [27]. It is developed according to the equations of Chen and Shu [30], ensuring higher accuracy and stability, particularly at high Reynolds numbers, at which no special treatment is needed.

The output of the ALM-LBM is validated with that of a 2-D Scale-Adaptive Simulation (SAS), demonstrating excellent agreement, and is subsequently compared with the solution of a VM. SAS is a turbulence modeling technique in CFD that adapts to resolved flow scales, allowing more accurate simulations of unsteady turbulent flows than traditional methods, without the high computational cost of full Large Eddy Simulation (LES).

Section 2 outlines the method, starting from the model description. Section 2.1 introduces the SHS-LBM and the ALM in Section 2.2. Section 2.3 illustrates the details of ALM corrections, including flow curvature correction (Section 2.3.1), dynamic stall (Section 2.3.2), and the Smagorinsky turbulence model (Section 2.3.3). Section 3 presents the numerical simulations, starting with sensitivity analysis on kernel size and mesh resolution in Section 3.1. Section 3.2 provides validation using a 2-D SAS and VM. Simulations at high and low tip-speed ratios are covered in Section 3.3 and Section 3.4, respectively. Lastly, a wake analysis is presented, examining how vorticity is shed downstream and how wake asymmetry diffuses.

2. Method

2.1. Simplified and Highly Stable Lattice Boltzmann Method

The SHS-LBM integrated with the ALM was implemented in a customizable MATLAB R2023a code, enabling its extension to 3-D modeling. The mesoscale particles in the LBM are obliged to follow defined directions according to Mohamad [33]. In 2-D fluid-flow simulations, the D2Q9 velocity set is used (Table 1, Figure 1), where the streaming velocity is characterized by nine velocity vectors denoted by the index k. The spatial discretization is given in a cartesian grid as (i,j), where x = (x_i,y_j) = (i∙Δx, j∙Δy).

The velocity set is built according to the Gauss–Hermite quadrature rule and must ensure continuity and momentum conservation. Lattice isotropy is a requirement to be satisfied by the weights w_i up to fifth-order integration via Hermite polynomials [19].

The discrete lattice Boltzmann equation (LBE) states the relationship between the distribution function at the new location and time step and the distribution function at the previous location and time step, with τ denoting the relaxation time. The LBE (Equation (1)) is not solved explicitly in the SHS-LBM, but predictor and corrector steps are developed [30].

f_k (x + c_k × Δt,t + Δt) = f_k (x,t) − 1/τ × [f_k (x,t) − f_k ^eq (x,t)]

(1)

Since the time step in the SHS-LBM must be equal to the lattice spacing (Δt = Δx, Δy), a unitary lattice particle speed (c = Δx/Δt) implies the speed of sound and fluid pressure to yield c_s = c/√3 = 1/√3 and p = ρ/c_s² = ρ/3, respectively.

The proportionality between τ and the kinematic shear viscosity ν and the time step is thus specified in Equation (2). To comply with the physical-scale Reynolds number Re, ν is calculated directly from Equation (3), with D denoting the physical turbine diameter (i.e., 1.0 m in this work) and u₀ denoting the lattice-scale streamwise inflow velocity of 0.1. An advantage of this method is that the simulation remains stable at high Re, also satisfying the necessary stability condition, i.e., τ > ½. [19].

τ = ν/(c_s² × Δt) + ½

(2)

ν = D × u₀/Re

(3)

Since zeroth and the first moments of the non-equilibrium distribution f_k ^neq are zero [30], the macroscopic variables ρ(x,t) and u(x,t) are recovered only from the moments of the distribution f_k, i.e., conservation of mass and momentum. However, intermediate density ρ* and momentum J* are first calculated in the predictor step (Equations (4) and (5)), avoiding any external forcing term.

ρ*(x,t) = Σ_k f_k ^eq (x − c_k × Δt, t − Δt)

(4)

J*(x,t) = Σ_kc_k × f_k ^eq (x − c_k × Δt, t − Δt)

(5)

The equilibrium distribution f_k ^eq in Equation (6) is second-order accurate, with w_k denoting the weights of the velocity set and c_k denoting discrete velocities. Density ρ and fluid velocity u are those corresponding to the previous streaming indices (x − c_k × Δt, t − Δt).

f_k ^eq (x,t) = w_k × ρ × [1 + (u × c_k)/c_s² + ½ × (u × c_k)²/c_s⁴ − ½ × (u × u)/c_s²]

(6)

The boundary conditions are applied to the intermediate variables. Hence, the macroscopic variables ρ(x,t) and u(x,t) are achieved from the corrector step disclosed in Equations (7)–(9). The macroscopic density equals the intermediate density (Equation (7)), and the macroscopic momentum J(x,t) is corrected according to Equation (8), including the external forcing term f_x generated from the ALM.

ρ(x,t) = ρ*

(7)

J(x,t) = J* + (τ − 1) × Σ_k c_k × f_k ^eq (x + c_k × Δt, t) - (τ − 1) × ρ(x,t − Δt) × u(x,t − Δt) + f_x × Δt

(8)

Finally, the macroscopic velocity u(x,t) is recovered from the momentum divided by ρ(x,t) in Equation (9). Next, the boundary conditions are applied to macroscopic variables.

u(x,t) = J/ρ

(9)

2.2. Actuator Line Model

The ALM relies on the 2-D blade element model, developed by interpolating lift and drag coefficients (C_l and C_d) from the static polar database at each angle of attack α. Notably, α is determined by solving the velocity triangle on the blade’s reference frame. The normal and tangential velocity components are denoted by u_n and u_t in Equations (10)–(12), the radian velocity denoted by ω, the rotor radius denoted by R, and the azimuthal angle denoted by ϑ. Since no pitch mechanism is introduced, the angle of attack α coincides with the flow angle ϕ.

u_n = u_x × sin(ϑ) − u_y × cos(ϑ)

(10)

u_t = u_x × cos(ϑ) + u_y × sin(ϑ) + ω × R

(11)

α = tan⁻¹ (u_n/u_t)

(12)

Lift and drag forces are determined from the interpolated C_l and C_d in Equations (13) and (14), introducing c and V_rel, as the blade chord and the relative wind speed, respectively.

l = ½ × ρ × c × V_rel² × C_l

(13)

d = ½ × ρ × c × V_rel² × C_d

(14)

Normal and tangential forces (f_n and f_t) are specified as a function of lift l and drag d in Equations (15) and (16). Hence, the blade loads in Cartesian coordinates f_x_,b and f_y_,b are developed according to Equations (17) and (18).

f_n = l × cos(α) + d × sin(α)

(15)

f_t = l × sin(α) − d × cos(α)

(16)

f_x_,b = f_n × sin(ϑ) − f_t × cos(ϑ)

(17)

f_y_,b = - f_n × cos(ϑ) − f_t × sin(ϑ)

(18)

Introducing the spacing r = (x − x_p_,b)² + (y − y_p_,b)² with x_p_,b and y_p_,b denoting the actual blade position in the global reference frame, and denoting the blade loads f_x_,b and f_y_,b for each blade j, the total force distribution in axial and crosswise directions (f_x, f_y) is derived from Equations (19) and (20), dictating a Gaussian-shaped force distribution in the cartesian space, with a kernel radius ε shaping the flow field around the airfoil.

f_x = Σ_b f_x_,b /(π × ε²) × (−(r/ε)²)

(19)

f_y = Σ_b f_y_,b /(π × ε²) × (−(r/ε)²)

(20)

According to Rullaud et al. [27], alternative definitions can better represent physical phenomena, such as non-isotropic projections that simulate force distribution around an airfoil, as discussed by Churchfield et al. [34]. However, the cylindrical force distribution along the entire blade span provides a more accurate representation of the geometric discontinuity at the blade tip [35]. Consequently, this model allows for a more precise resolution of the vortices shed from the blade tips [36].

Martinez-Tossas et al. [37] suggested that ε should be proportional to the blade chord c, specifically 0.2c, to minimize the deviation between the solutions of the actuator line and potential flow for a flat plate or Joukowsky airfoil. Following this recommendation, Melani et al. [38] used the same kernel size as Martinez-Tossas in their simulations, with a grid resolution of dh ≤ ε/2. This approach is consistent with Troldborg [39], who advised that ε should be at least twice the size of the local grid cell to reduce spurious oscillations in the flow field.

2.3. ALM Corrections

Flow curvature and dynamic stall corrections are used on the airfoil data, and a subgrid-scale turbulence model is added.

2.3.1. Flow Curvature

The effect of the curvature causes a variation in the apparent camber and incidence angle of the rotating blades [40]. Virtual incidence significantly impacts blade aerodynamics by shifting the lift coefficient (C_l) curve to the left, especially at a high blade chord to turbine radius ratio (c/R). A conformal mapping technique [41] was derived by modeling a flat plate in a circular trajectory, allowing correcting the angle of attack to account for the flow curvature. The same correction was developed in an ALM-CFD model [42]. Using uncorrected airfoil data, the geometric angle of attack is corrected as in Equation (21), with V_rel, ω, c, and x_0r denoting the relative wind speed, rotational speed, chord, and blade normalized mounting point, respectively.

α = α + ω × c/V_rel × (x_0r + 1/4)

(21)

At high rotor speeds, V_rel ≈ ω × R leads the flow correction model to shift the C_l curve leftward by a factor of ¼ × (c/R). The blade is mounted at the mid-chord, so its normalized distance from c/4 is accounted for by setting the mounting point x_0r to 0.

2.3.2. Dynamic Stall

The modeling of dynamic stall [43] encompasses boundary layer separation and vortex shedding while interacting with the airfoil. As shown in [44], inviscid lift coefficient (C_l_,inv), fully separated lift coefficient (C_l_,fs), and separation point (f_s^st) are derived from known polar data, covering both positive and negative angles of attack. The dynamic stall for the lift coefficient is governed by Equation (22), where f_s denotes the separation function defined between [0, 1]. The lift coefficient C_l is a linear combination of the inviscid flow with no separation (extrapolated from the static airfoil data in the linear region) and fully separated flow for a flat plate with a sharp leading edge [45]. Thus, f_s = 0 corresponds to a fully separated flow, and f_s = 1 to a completely inviscid flow. A single value f_s = f_s^st gives exactly the static lift coefficient for a given angle of attack, and the dynamic stall model of Equation (23) always forces f_s to go towards f_s^st.

C_l = f_s × C_l_,inv(α) + (1 − f_s) × C_l_,fs(α)

(22)

The ordinary differential equation in a discrete form, Equation (23), is solved for the separation function f_s at the new time step (t + Δt), in which τ_ds denotes the constant time delay influencing blade loads due to cyclical variations in the angle of attack (Equation (24)).

When the boundary layer begins to separate from the trailing edge and gradually moves toward the leading edge at increasing angles of attack, Øye’s model predicts this time delay.

f_s(t + Δt) = f_s^st + (f_s(t) − f_s^st) × e^−Δt/τds

(23)

τ_ds = 4 × c/V_rel

(24)

2.3.3. Smagorinsky Turbulence Model

The Smagorinsky model [46] is applied to high Re flows, where the inertial region and energy dissipation are exhibited. The largest turbulent structures are responsible for the energy dissipation rate resolved in large eddy simulation (LES), leading to the computation of the eddy viscosity (ν_t). Although the finest scales responsible for dissipation are not captured [47], the model allows recovering some aspect of the anisotropic behavior of turbulence by accounting for the strain rate tensor. The subgrid-scale stresses are proportional to the local strain rate of the resolved flow (

\bar{S}

_ij) [48], with i and j denoting the Cartesian directions. In Einstein’s notation,

{\bar{S}}_{i j} = \frac{1}{2} \times (\frac{\partial u_{j}}{\partial x_{j}} + \frac{\partial u_{i}}{\partial x_{i}})

(25)

The magnitude of the large-scale strain rate tensor is given by

| \bar{S} | = (2 \times Σ_{i j} {\bar{S}}_{i j} {\bar{S}}_{i j})^{\frac{1}{2}}

(26)

Equation (27) yields the eddy viscosity ν_t, with C = 0.1 denoting the Smagorinsky constant and Δh denoting the mesh resolution. Hence, the total kinematic viscosity is given by the sum of bulk kinematic viscosity ν₀ and turbulent viscosity ν_t (Equation (28)), yielding the relaxation factor for the LBM (Equation (29)). Applications with the Smagorinsky turbulence model and SHS-LBM were investigated by Jiang et al. [49].

ν_{t} = (C \times Δ h)^{2} \times | \bar{S} |

(27)

ν_tot = ν₀ + ν_t

(28)

τ_tot = ν_tot/(c_s² × Δt) + ½

(29)

3. Simulations

The two-bladed H-Darrieus rotor tested at the Open Jet Facility (OJF) at TU Delft [31,32] is investigated with 2-D simulations. The post-processing of Particle Image Velocimetry (PIV) data provided both 2-D and 3-D flow fields [31], as well as azimuthal blade force profiles [32] for the turbine. The PIV analysis established a benchmark for validating and comparing numerical models. The rotor was operated at tip-speed ratios (λ) of 4.5 and 2.0, corresponding to average diameter-based Reynolds numbers of 160,000 and 80,000, respectively. Velocity fields were captured at eight azimuthal positions within the midplane of the blade. While velocity fields were directly derived from PIV, the loads were determined using an integral approach. The primary objective of evaluating these tip-speed ratios was to discern the impact of dynamic stall, which is absent at the higher λ but predominant at the lower λ. Given the complexity of numerically modeling the dynamic stall phenomenon, this benchmark for a VAWT experiencing dynamic stall is valuable.

The rotor characteristics are presented in Figure 2 and Table 2, and the LBM operating conditions are detailed in Table 3, illustrating stable operation at λ = 4.5 and unstable operation at λ = 2.0, marked by dynamic stall extending over a significant portion of one rotor revolution. The rotor has a blade aspect ratio of 16.67, i.e., the ratio of blade length to chord, and a blade pitch angle of zero degrees.

The outputs C_n* and C_t* plotted in Section 3.2 and Section 3.3 represent the normal and tangential force coefficients equated as f_n/(½ρRu₀²) and f_t/(½ρRu₀²), respectively, to align with experimental data of Tescione et al. [31]. In these definitions, f_n and f_t denote the blade loads, u₀ denotes the free stream velocity in lattice units [l.u.], and R denotes the turbine radius.

Sensitivity on kernel size, mesh size, and final validation are pursued with the coefficient of determination R² [50] assessing the numerical solution’s accuracy relative to observed data or analytical solution. SAS is used as a reference for validation instead of experimental data in this work. Hence, R² is calculated by measuring the deviation between the reference points F_ref and the numerical outputs F to the deviation between F_ref and the average value of the observed data F_avg. It is expressed as:

R² = 1 − Σ_ϑ (F_ref(ϑ) − F(ϑ))/Σ_ϑ (F_ref(ϑ) − F_avg)

(30)

The SAS solution was first averaged over the last ten revolutions to smooth the final curve, from which F_avg was also calculated. This approach allowed the calculation of more accurate deviations F_ref(ϑ) − F(ϑ).

3.1. Kernel Size and Mesh Resolution

The chosen blockage ratio of 5% complies with the CFD literature [51,52], so a domain with a 20D cross-stream width was deemed acceptable (Figure 3). The distance from the inlet boundary to the rotor shaft was set to 10D to prevent overestimating turbine performance. The same distance was set downstream of the rotor to minimize the influence of the domain outlet on the turbine performance, according to Rezaeiha et al. [53].

u_x₁ = u₀ and U_y₁ = v₀ were imposed at the inlet boundary. A second-order extrapolation is set for mass density such as ρ₁ = 4/3 × ρ₂ − 1/3 × ρ₃ while pressure yields simply p₁ = ρ₁ × c_s². The second-order extrapolation is executed on the velocities at the other boundaries. The extrapolation at the boundaries at each rotor side avoids the flow acceleration that occurs with symmetry boundary conditions.

To allow the wake to develop completely within the computational domain [53], data sampling was performed after 40 rotor revolutions. Flow curvature correction and dynamic stall corrections were activated for the sensitivity analysis.

The actuator line’s Gaussian distribution ideally mirrors the chordwise pressure distribution. Although achieving high precision is often impractical [39], the sensitivity of the solution to the kernel size ε was investigated proportionally to the chord length (i.e., ε = 1c, 0.5c, 0.25c) using coarse and fine meshes. Multiple tests assessed the solution’s independence from the mesh size as the mesh resolution Δh increased from 0.040 to 0.005.

The relative wind speed denoted as V_rel was used as a control parameter in multiple tests assessing the solution’s independence from the mesh size. Larger ε yield smoother solutions while low ε result in noticeable oscillations and spurious solutions. These oscillations occur mainly at a low mesh resolution, insufficient to resolve each ε accurately.

Figure 4 shows only the sensitivity of the loads to the kernel size at the finest mesh resolution which allowed the highest correlation with SAS data, as depicted in Section 3.2. The solutions for mesh resolutions of 0.005 and 0.010 appear nearly identical, except for a slight deviation downwind. Therefore, the mesh resolution of 0.005 was considered accurate to deal with the next simulations with sufficient computational efficiency. The simulation executed with ε = 0.5c generated the smoothest solution, providing proof of Δh independence.

3.2. Validation with a 2-D CFD Approach and VM

SAS is a modern turbulence modeling approach used in CFD. It bridges the gap between traditional Reynolds-Averaged Navier–Stokes (RANS) models and the computationally expensive Large Eddy Simulation (LES) and hybrid RANS/LES models. It adjusts its approach based on the resolved scales of flow, improving the simulation of unsteady turbulent flows like dynamic stall phenomena, without the high cost of full LES. The SAS approach offers a powerful tool for simulating unsteady turbulent flows more accurately than traditional URANS models. Its ability to adapt to the resolved scales of the flow makes it particularly effective for complex phenomena such as dynamic stall in vertical axis wind turbines [54,55].

Unlike conventional URANS, the SAS method incorporates an additional turbulence scale, known as the von Kármán length scale, which enables the model to identify and resolve turbulent structures as they form. While URANS is limited to capturing only large-scale unsteady motions, the SST-SAS model dynamically adapts to the resolved flow scales, allowing it to capture a broader spectrum of turbulence, particularly in regions with flow separation. This capability allows the SAS model to behave similarly to LES in areas dominated by unsteady flow, while still retaining the stable flow characteristics of RANS. As a result, SAS provides a more refined turbulence resolution in complex flow conditions, enhancing the accuracy of simulations without significantly increasing computational costs. The SAS approach, proposed by Menter and Egorov [56] and implemented in ANSYS Fluent, is used in classical ω-based turbulence models such as the well-known two-equation k-ω SST model or the relatively new Transition SST model. A challenge with these formulations is that the transport equation for ω (or ε) resolves the dissipative scales rather than the larger scales of turbulence. To address this, Menter and Egorov extended the ω equation by building on a more systematic approach, originally developed by Rotta, where an exact transport equation for turbulent kinetic energy (k) is multiplied by the turbulence length scale (L). Rotta’s equation focuses on capturing the large turbulent scales, making it a solid foundation for developing models on a term-by-term basis [56,57].

In this study, the same 2-D turbine rotor model was prepared (Figure 5), omitting struts and the rotating tower. Rogowski et al. [58] and Rezaeiha et al. [53] showed that under such conditions, the 2-D CFD approach yields reasonable results for aerodynamic loads, while the effects of tip losses are less significant.

In the case of the URANS approach with SAS turbulence modeling, a rectangular computational domain with dimensions length to width of 35D by 20D was utilized. The rotor axis was positioned at a distance of 10D from the vertical front edge of the domain, matching the setup of the LBM model.

Within the rectangular computational domain, a circular core with a diameter of 1.5D was used around the rotor to allow its rotation during the simulation, employing the sliding mesh technique. The following boundary conditions were applied: velocity inlet (at the domain inlet), pressure outlet (at the outlet), and symmetry (on the tunnel side walls). The interface condition ensured data exchange between the stationary rectangular domain and the rotating core. CFD 2-D simulations were performed under the same flow conditions as in the LBM approach. At the domain inlet, turbulence intensity and scale were assumed to be 1% and 0.001 m, respectively, and these parameters were not further investigated. They were established based on previous analyses [59]. It is also known that the turbulence intensity in the TU Delft tunnel was low (0.24%) [32]. Assuming the rotor is at a significant distance from the domain inlet, the kinetic energy of turbulence will decrease significantly [60], and as a result, the turbulence intensity experienced by the rotor blades will be comparable to real wind tunnel conditions. The simulations used incompressible Navier–Stokes equations assuming air density and viscosity of 1.225 kg/m³ and 1.79 × 10⁻⁵ kg/(m∙s), respectively.

The numerical grid around the rotor blades was based on the findings of Rogowski et al. [61] and Rogowski [59]. It consisted of a structured grid near the blade edges and an unstructured mesh with quadrilateral elements in the remaining area. The total number of nodes on the airfoil edges was 830. The thickness of the first layer of the mesh near the airfoil edges was 3.6 × 10⁻⁶ m, ensuring that the wall y⁺ parameter did not exceed 1.0. The growth rate of the grid normal to the profile edge was set to 1.1, according to [62]. The number of structured mesh layers is 40.

A mesh sensitivity test was presented in the papers [59,61]. The computational domain along with the numerical grid is presented in Figure 3. The 2-D CFD simulations were conducted for 25 full rotor revolutions, assuming a time step size equivalent to an azimuth increment Δθ of 0.5 degrees. Rezaeiha et al. [53] conducted several analyses on this investigated wind turbine for Δθ between 0.5 and 0.05 degrees. Their findings indicated that within this range of azimuth increments, the difference in the rotor power coefficient is negligible for a tip-speed ratio of 4.5. Simulating 25 full rotor revolutions is sufficient to ensure that after 10 full revolutions, the effect of the initial conditions is neglected [59].

For the 2-D CFD approach, a pressure-based solver was utilized. The turbulence model employed was the four-equation Transition SST model with the scale-resolving SAS model. For such a small blade-based Reynolds number of 160,000, the transition model is more adequate than the classical two-equation turbulence models, such as the k-ω or k-ε models [61,63]. The SIMPLE scheme was used for pressure–velocity coupling. The bounded central differencing scheme was applied for the spatial discretization of the governing equations. The second-order scheme was also applied for the remaining turbulence model equations. The second-order transient scheme was additionally used in the transient formulation. Under-relaxation factors were set to 0.3 for pressure, 1 for density and body forces, and 0.8 for other quantities. The convergence criteria level was set to 10⁻⁸ for all iterated values. The maximum number of iterations per time step was set to 20.

The 2-D VM was developed by determining the bound circulation from the lift coefficient using the Kutta–Joukowski theorem, enabling the calculation of induced velocities through a simplified form of the Biot–Savart law. To satisfy the conservation of total circulation through Kelvin’s theorem, new vortices are introduced at each time step and subsequently convected into the wake. This model was developed following the approach outlined by Bangga et al. [1].

3.3. High Tip-Speed Ratio

The baseline achieved from static airfoil data is compared with the output incorporating the ALM corrections at high rotor speeds (λ = 4.5) in Figure 6, Figure 7 and Figure 8. Flow curvature (FC) and dynamic stall (DS) modify the angle of attack from the uncorrected curve (BASELINE), impacting the static C_l and C_d. Additionally, the Smagorinsky turbulence model (T) is enabled to account for the effect of eddy kinematic viscosity. The solutions for ε = 0.03 and ε = 0.06 are shown on the left and right sides, respectively.

The flow curvature causes C_n* and C_t* to increase upwind and decrease downwind from the baseline, but the shift in blade loads is more pronounced for the normal force. Low values of C_l and α_v result in poor aerodynamic performance, often leading to counter-torque. Power extraction is absent where the tangential force C_t* is negative, which occurs when the virtual incidence α_v inverts while the lift coefficient C_l remains positive.

At high rotor speeds, the effect of dynamic stall on the angle of attack is limited. It increases C_n* at the cusp upwind and reduces it slightly downwind. This trend also applies to C_t*, suggesting higher power extraction than the baseline upwind and lower power extraction downwind.

In regions where |α_v| < |α_g|, the rotor’s performance is suboptimal due to very low C_l at negative angles of attack, as shown in [40]. A positive torque is generated when the lift coefficient C_l and the virtual incidence angle α_v have the same sign. The flow curvature correction induces higher angles of attack upwind, where α_v is positive, and lower angles of attack downwind, where α_v is negative. The virtual incidence increases the maximum static angle of attack by several degrees, yet it prevents the occurrence of deep stall even at high rotor speeds (see Figure 8).

The Smagorinsky model introduces a change in the angle of attack, which is more pronounced upwind, leading to a decrease in C_n* and C_t* in the upwind region and an increase in the downwind region for ε = 0.5c. However, for ε = c, the model is almost ineffective, as it barely alters the angle of attack. This suggests that the Smagorinsky model depends on kernel size, whereas other corrections remain independent. Moreover, turbulence is inherently integrated into the LBM model, as τ relaxes the LBE at high Reynolds numbers. Additionally, while the Smagorinsky model is effective for 3-D LES, it may lead to inaccuracies in 2-D simulations. This is because it assumes isotropic turbulence and tends to over-dissipate energy, resulting in poor capture of 2-D turbulent structures. Since solutions independent of kernel size and mesh resolution are desirable, the Smagorinsky model is excluded from the next simulations and charts.

The correlation between the predicted normal force coefficient C_n* of a kernel size ε = 0.5c and the reference (SAS) is high (R² = 0.989), with greater deviations observed downwind. None of the numerical models account for the interference of the rotating shaft, which partly explains the divergence of C_n* in the downwind path compared to the experimental data (PIV) (Figure 9, left-hand side). According to Castelein et al. [32], the uncertainty in normal force measurements is relatively low, whereas it is higher for tangential force measurements. A negative tangential force coefficient C_t* is observed from leeward to nearly upwind, and from windward to downwind, consistent with the findings of Castelein et al. A high correlation is presented also by the LBM C_t* (R² = 0.937) and the reference.

However, C_t* does not accurately predict the experimental data (Figure 9, right-hand side). This discrepancy is mainly observed downwind, where the collision between the blades and the shaft’s wake causes a sudden decrease in tangential loads due to turbulence, as observed experimentally.

The LBM output is thus compared with the VM output. In this case, the correlation of C_n* and C_t* is also convincing (R² = 0.979 and R² = 0.923, respectively). The downwind deviation is expected, as the VM is less accurate in regions with higher turbulence, such as the wake. The potential flow-based VM does not account for turbulence and viscosity, leading to incomplete modeling of fluid dynamics effects like flow separation and wake interaction with blades. In contrast, these effects are captured by the LBM.

3.4. Low Tip-Speed Ratio

The SAS output at the lowest tip-speed ratio in Figure 10 shows a highly unsteady flow with marked oscillations, indicating a permanent dynamic stall condition. Conversely, the C_n* and C_t* simulated with the VM closely capture most of the separated flow polar data, resulting in a more coherent pattern.

A test executed with the LBM for Re = 80,000 and λ = 2.0 results in a solution almost superimposable on the VM output, with R² of 0.994 for C_n* and 0.920 for C_t*. In this case, the most suitable kernel size is chosen equal to the blade chord for the best approximation, even though the solution at ε = 0.5c shows a slight deviation. This observation results from the significant reduction in induction in both the LBM and VM at λ = 2.0 compared to λ = 4.5 (Figure 11, left-hand side), due to the effects of unsteady flow. The increase in the angle of attack (Figure 11, right-hand side) illustrates a decline in aerodynamic performance, leading to a decrease in tangential force and, consequently, a reduction in the expected mechanical torque.

Simulating a rotor at a low tip-speed ratio is challenging due to continuous boundary layer separation, a complex three-dimensional phenomenon. An accurate simulation of turbulent flow patterns necessitates a fully resolved flow around the airfoil, which can be computationally intensive. Assuming that the ALM effectively captures the trend of blade loads, refining the mesh size and kernel size is a viable strategy to enhance the detail and accuracy of the simulation.

As previously mentioned, the SAS approach shows certain oscillations in aerodynamic loads compared to ALM. This is not due to the mesh resolution used in the CFD calculations, but rather to the turbulence model γ − Re_ϑ. The turbine blade, operating under such flow conditions (with a chord-based Reynolds number below 200,000), experiences the formation of laminar separation bubbles on its surfaces. The transition turbulence model allows for a reasonable approximation of these phenomena.

It is also important to emphasize that although the ALM approach uses airfoil polars measured for the airfoil in low Reynolds number regimes and the effects of laminar separation bubbles are visible in the characteristics of lift coefficient C_l(α) and drag coefficient C_d(α) in the ALM approach, these forces are treated as point forces and introduced into the flow as a certain distribution. In contrast, the SAS approach allows for obtaining similar average aerodynamic characteristics of the airfoil as in the case of measurements, but it considers these forces as instantaneous distributions of the pressure coefficient, considering laminar separation bubbles.

In the case of such a low λ, laminar separation bubbles are not the only phenomena occurring during the flow and influencing load oscillations. Important aspects include, in the upstream part of the rotor, dynamic stall, which in the SAS approach depends on the turbulence model itself, while in the ALM approach, it depends on the applied dynamic model or its absence. Moreover, in the downstream part of the rotor, blade–wake interaction plays a significant role.

3.5. Wake Analysis

Contours of the streamwise and crosswise flow velocities, along with the instantaneous vorticity in the vertical plane, are developed for the rotor at various downstream x/D stations for both tip-speed ratios. The flow velocities are normalized with the free-stream flow speed (u_x/u₀ and u_y/u₀), while the vorticity is normalized using the chord and free-stream velocity (ω_zc/u₀).

Despite being normalized to u₀, the streamwise velocity exceeds 1 at the rotor sides. This is attributed to the blockage, which accelerates the flow bypassing the rotor. This effect is somehow reduced by choosing an extrapolation boundary condition instead of a symmetry boundary condition for these domain boundaries [9]. However, the impact of this difference is limited, especially considering the maximum prescribed blockage ratio of 5%.

The wake profile valid for λ = 4.5 illustrated on the left-hand side of Figure 12 exhibits moderate asymmetry, with an off-center velocity minimum and a steep gradient. As the wake convects downstream, the velocity profile becomes smoother, with enhanced wake expansion and a slight deflection towards the leeward side due to the asymmetric load distribution.

At x/D = 1, velocity gradients are sharper due to the blade’s interaction with the flow. The presence of positive and negative crosswise velocity gradients on either side indicates the formation of trailing vorticity, which is then shed into the wake. The skewness of the wake is more pronounced at a tip-speed ratio of λ = 2.0 depicted on the right-hand side of Figure 12, where the velocity profiles exhibit strong oscillations due to the rotor’s unsteady operating regime. As the asymmetry diffuses downstream, it suggests ongoing wake self-induction.

Figure 13 shows a contour plot of out-of-plane vorticity at λ = 4.5 (on the top) and λ = 2.0 (on the bottom). At the highest tip-speed ratio, the blades undergo periodic changes in the angle of attack and relative velocity vector, leading to a cyclical variation in the bound circulation around the airfoil.

This variation results in an equivalent amount of circulation being shed downstream, forming pairs of vortex sheets that are convected into the wake from the trailing edges of the airfoils. Since the rate of change of circulation is not constant and varies in sign as it crosses the midplane axis both upwind and downwind, the resulting vortex sheets also change accordingly [31]. The deformation of the wake begins roughly one rotor diameter downstream, initiating the roll-up process that generates large-scale vortices and causes the wake to skew slightly toward the leeward side. This roll-up is driven by peaks in vorticity resulting from the interaction of overlapping blade wakes.

After about three rotor radii, the individual blade wakes converge into two continuous vortex sheets, which remain distinct and evenly spaced as they extend downstream. At the lowest tip-speed ratio, large, unstable vortex structures are shed from the trailing edges of the blades. These vortices are more pronounced and irregular due to the blades’ diminished control over the airflow. As these chaotic structures merge, turbulence increases, leading to greater energy dissipation. The blades experience a higher angle of attack, leading to flow separation on the suction side, generating additional vortices and degrading the rotor’s performance.

4. Conclusions

An actuator line model integrated with the LBM enables accurate resolution of the flow field generated by a two-bladed vertical axis wind turbine. The outputs of the SHS-LBM and SAS show a high degree of correlation, validating the model. Only minor deviations are observed when the LBM result is compared to the 2-D VM solution.

The SHS-LBM, incorporating flow curvature and dynamic stall corrections, effectively adapts to stable and unstable operating conditions. Conversely, the Smagorinsky turbulence model introduces inaccuracies, indicating its limited applicability in 2-D flows. Turbulence is inherently accounted for in the SHS-LBM, as the relaxation time is determined from a very low kinematic viscosity.

The choice of the kernel size proves to be a critical phase before grid refinement, with a size equal to half the chord length proving effective for accurate flow field resolution, particularly at λ = 4.5. However, the optimal value for ε remains uncertain [37], as it may vary across different tip-speed ratios, reflecting a limitation of the method. Despite this uncertainty, the LBM-ALM output shows a high correlation with the URANS-SAS numerical solution, and good agreement is also noted with the 2-D VM at high tip-speed ratios. During unsteady rotor operation, reduced induction leads to a solution approximating that of the VM.

Contour plots of streamwise and crosswise flow velocities reveal increased wake skewness, illustrating the formation and merging of vortices into continuous vortex sheets a few rotor diameters downstream, due to overlapping blade wakes. This indicates a stable spatial wake pattern, with instability developing further downstream, ultimately leading to a turbulent wake as root vortices destabilize early in the flow.

The model demonstrates a strong correlation between numerical solutions and the SAS outputs. Considering the reduced computational time and excellent parallelization capabilities, the lattice Boltzmann solver emerges as a leading tool for simulating vertical-axis wind turbines, offering accuracy and efficiency comparable to traditional CFD methods.

Author Contributions

Conceptualization, L.C. and M.O.L.H.; methodology, L.C. and M.O.L.H.; software, L.C. and M.O.L.H.; validation, L.C., M.O.L.H. and K.R.; formal analysis, L.C., M.O.L.H. and K.R.; investigation, L.C., M.O.L.H. and K.R.; resources, L.C., M.O.L.H. and K.R.; data curation, L.C., M.O.L.H. and K.R.; writing—original draft preparation, L.C.; writing—review and editing, L.C., M.O.L.H. and K.R.; visualization, L.C.; supervision, M.O.L.H.; project administration, M.O.L.H. All authors have read and agreed to the published version of the manuscript.

Funding

The part concerning 2-D LBM received no external funding. The part concerning 2-D CFD was supported by POB Energy of Warsaw University of Technology within the Excellence Initiative: Research University (IDUB) programme (Grant No. 1820/355/Z01/POB7/2021), and through the infrastructure of the Interdisciplinary Centre for Mathematical and Computational Modelling (ICM) at the University of Warsaw under computational allocations no. GB83-33.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

VAWT—Vertical Axis Wind Turbine; SHS-LBM—Simplified and Highly Stable Lattice Boltzmann Method; MRT—Multi-Relaxation Time; ALM—Actuator Line Model; CFD—Computational Fluid Dynamics; RANS—Reynolds-Averaged Navier–Stokes; URANS—Unsteady Reynolds-Averaged Navier–Stokes; SAS—Scale-Adaptive Simulation; VM—Vortex Method; FC—Flow Curvature; DS—Dynamic Stall.

References

Bangga, G.; Dessoky, A.; Wu, Z.; Rogowski, K.; Hansen, M.O.L. Accuracy and Consistency of CFD and Engineering Models for Simulating Vertical Axis Wind Turbine Loads. Energy 2020, 206, 118087. [Google Scholar] [CrossRef]
Strickland, J.H. Advanced Energy Projects Dept. In The Darrieus Turbine: A Performance Prediction Model Using Multiple Streamtubes; SAND, Sandia Laboratories: Albuquerque, NM, USA, 1979. [Google Scholar]
Paraschivoiu, I. Double-Multiple Streamtube Model for Studying Vertical-Axis Wind Turbines. AIAA J. Propuls. Power 1988, 4, 370–378. [Google Scholar] [CrossRef]
Sørensen, J.N.; Shen, W.Z. Numerical Modelling of Wind Turbine Wakes. J. Fluids Eng. 2002, 124, 393–399. [Google Scholar] [CrossRef]
Hansen, M.O.L.; Sørensen, D.N. CFD Model for Vertical Axis Wind Turbine; European Wind Energy Conference and Exhibition: Copenhagen, Denmark, 2001; pp. 485–487. [Google Scholar]
Rogowski, K. Numerical Studies on Two Turbulence Models and a Laminar Model for Aerodynamics of a Vertical-Axis Wind Turbine. J. Mech. Sci. Technol. 2018, 32, 2079–2088. [Google Scholar] [CrossRef]
Rogowski, K.; Hansen, M.O.L.; Galih, B. Performance Analysis of a H-Darrieus Wind Turbine for a Series of 4-digit NACA Airfoils. Energies 2020, 13, 3196. [Google Scholar] [CrossRef]
Belabes, B.; Paraschivoiu, M. Numerical Study of the Effect of Turbulence Intensity on VAWT Performance. Energy 2021, 233, 121139. [Google Scholar] [CrossRef]
Cacciali, L. Lattice Boltzmann Method for VAWT Aerodynamics. Wind Energy. Master’s Thesis, Technical University of Denmark, DTU, Lyngby, Denmark, 2024. [Google Scholar]
Shamsoddin, S.; Porté-Agel, F. A Large-Eddy Simulation Study of Vertical Axis Wind Turbine Wakes in the Atmospheric Boundary Layer. Energies 2016, 9, 366. [Google Scholar] [CrossRef]
Troldborg, N.; Sørensen, J.N.; Mikkelsen, R. Numerical Simulations of Wake Characteristics of a Wind Turbine in Uniform Inflow. Wind Energy 2010, 13, 86–99. [Google Scholar] [CrossRef]
Bachant, P.; Goude, A.; Wosnik, M. Actuator Line Modeling of Vertical-Axis Turbines. arXiv 2018, arXiv:1605.01449. [Google Scholar]
Martinez-Ojeda, E.; Ordaz, F.J.S.; Sen, M. Vertical-Axis Wind-Turbine Computations Using a 2D Hybrid Wake Actuator-Cylinder Model. Wind Energy Sci. 2021, 6, 1061–1077. [Google Scholar] [CrossRef]
Madsen, H. The Actuator Cylinder: A Flow Model for Vertical-Axis Wind Turbines. Ph.D. Thesis, Aalborg University Centre, Aalborg, Denmark, 1982. [Google Scholar]
Mendoza, V.; Goude, A. Wake Flow Simulation of a Vertical Axis Wind Turbine Under the Influence of Wind Shear. J. Phys. Conf. Ser. 2017, 854, 012031. [Google Scholar] [CrossRef]
Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Oxford University Press: Oxford, UK, 2001. [Google Scholar]
Nourgaliev, R.R.; Dinh, T.N.; Theofanous, T.G.; Joseph, D. The Lattice Boltzmann Equation Method: Theoretical Interpretation, Numerics and Implications. Int. J. Multiph. Flow 2003, 29, 117–169. [Google Scholar] [CrossRef]
Dünweg, B.; Ladd, A.J.C. Advances in Polymer Science; Springer: Berlin/Heidelberg, Germany, 2008; pp. 1–78. [Google Scholar]
Krüger, T.; Kusumaatmaja, H.; Kuxmin, A.; Shardt, O.; Silva, G.; Viggen, E.M. The Lattice Boltzmann Method. In Principles and Practice; Springer: Cham, Switzerland, 2017. [Google Scholar]
Geller, S.; Krafczyk, M.; Tölke, J.; Turek, S.; Hron, J. Benchmark computations based on Lattice-Boltzmann, Finite Element and Finite Volume Methods for laminar Flows. Comput. Fluids 2006, 35, 888–897. [Google Scholar] [CrossRef]
Asmuth, H.; Olivares-Espinosa, H.; Nilson, K.; Ivanell, S. The Actuator Line Model in Lattice Boltzmann Frameworks: Numerical Sensitivity and Computational Performance. J. Phys. Conf. Ser. 2019, 1256, 012022. [Google Scholar] [CrossRef]
Asmuth, H.; Olivares-Espinosa, H.; Ivanell, S. Actuator Line Simulations of Wind Turbine Wakes Using the Lattice Boltzmann Method. Wind Energy Sci. 2020, 5, 623–645. [Google Scholar] [CrossRef]
Glessmer, M.S.; Janssen, C.F. Using an Interactive Lattice Boltzmann Solver in Fluid Mechanics Instruction. Computation 2017, 5, 35. [Google Scholar] [CrossRef]
Schubiger, A.; Barber, S.; Nordborg, H. Evaluation of the Lattice Boltzmann Method for Wind Modelling in Complex Terrain. Wind Energy Sci. 2020, 5, 1507–1519. [Google Scholar] [CrossRef]
Latt, J.; Malaspinas, O.; Kontaxakis, D.; Parmigiani, A.; Lagrava, D.; Brogi, F.; Belgacem, M.B.; Thorimbert, Y.; Leclaire, S.; Marson, F.; et al. Palabos: Parallel Lattice Boltzmann Solver. Comput. Math. Appl. 2021, 81, 334–350. [Google Scholar] [CrossRef]
Ribeiro, A.F.P.; Muscari, C. Sliding Mesh Simulations of a Wind Turbine Rotor with Actuator Line Lattice-Boltzmann Method. Wind Energy 2023. Early View. [Google Scholar] [CrossRef]
Rullaud, S.; Blondel, F.; Cathelain, M. Actuator-Line Model in a Lattice Boltzmann Framework for Wind Turbine Simulations. J. Phys. Conf. Ser. 2018, 1037, 022023. [Google Scholar] [CrossRef]
d’Humières, D.; Ginzburg, I.; Krafczyk, M.; Lallemand, P. Multiple-Relaxation-Time Lattice Boltzmann Models in Three Dimensions. Phil. Trans. A 2002, 360, 437–451. [Google Scholar] [CrossRef] [PubMed]
Kupershtokh, A.L.; Medvedev, D.A.; Karpov, D.I. On Equations of State in a Lattice Boltzmann Method. Comput. Math. Appl. 2009, 58, 965–974. [Google Scholar] [CrossRef]
Chen, Z.; Shu, C. Simplified and Highly Stable Lattice Boltzmann Method—Theory and Applications. In Advances in Computational Fluid Dynamics; World Scientific Publishing Co. Pte. Ltd.: Singapore, 2021; Volume 5. [Google Scholar]
Tescione, G.; Ragni, D.; He, C.; Simão Ferreira, C.J.; Van Bussel, G.J.W. Near Wake Flow Analysis of a Vertical Axis Wind Turbine by Stereoscopic Particle Image Velocimetry. Renew. Energy 2014, 70, 47–61. [Google Scholar] [CrossRef]
Castelein, D.; Ragni, D.; Tescione, G.; Simão Ferreira, C.J.; Gaunaa, M. Creating a Benchmark of Vertical Axis Wind Turbines in Dynamic Stall for Validating Numerical Models. In Proceedings of the 33rd Wind Energy Symposium, AIAA Scitech, Kissimmee, FL, USA, 5–9 January 2015. [Google Scholar]
Mohamad, A.A. Lattice Boltzmann Method—Fundamental and Engineering Applications with Computer Codes, 2nd ed.; Springer: Cham, Switzerland, 2019. [Google Scholar]
Churchfield, M.; Schreck, S.; Martinez-Tossas, L.; Meneveau, C.; Spalart, P. An Advanced Actuator Line Method for Wind Energy Applications and Beyond. In Proceedings of the 35th Wind Energy Symposium, Grapevine, TX, USA, 9–13 January 2017. [Google Scholar]
Jha, P.K.; Schmitz, S. Actuator Curve Embedding an Advanced Actuator Line Model. J. Fluid Mech. 2018, 834, R2. [Google Scholar] [CrossRef]
Mohamed, O.S.; Melani, P.F.; Balduzzi, F.; Ferrara, G.; Bianchini, A. An Insight on the Key Factors Influencing the Accuracy of the Actuator Line Method for Use in Vertical-Axis Turbines: Limitations and Open Challenges. Energy Convers. Manag. 2022, 270, 116249. [Google Scholar] [CrossRef]
Martínez-Tossas, L.A.; Churchfield, M.J.; Meneveau, C. Optimal Smoothing Length Scale for Actuator Line Models of Wind Turbine Blades Based on Gaussian Body Force Distribution. Wind Energy 2017, 20, 1083–1096. [Google Scholar] [CrossRef]
Melani, P.F.; Balduzzi, F.; Ferrara, G.; Bianchini, A. Tailoring the Actuator Line Theory to the Simulation of Vertical-Axis Wind Turbines. Energy Convers. Manag. 2021, 243, 114422. [Google Scholar] [CrossRef]
Troldborg, N. Actuator Line Modeling of Wind Turbine Wakes. Ph.D. Thesis, Technical University of Denmark, Lyngby, Denmark, 2008. [Google Scholar]
Migliore, P.G.; Wolfe, W.P.; Fanucci, J.B. Flow Curvature Effects on Darrieus Turbine Blade Aerodynamics. J. Energy 1980, 4, 49–55. [Google Scholar] [CrossRef]
Goude, A. Fluid Mechanics of Vertical Axis Turbines: Simulations and Model Development. Doctoral Thesis, Uppsala University, Uppsala, Sweden, 2012. [Google Scholar]
Dyachuk, E.; Goude, A.; Bernhoff, H. Dynamic Stall Modeling for the Conditions of Vertical Axis Wind Turbines. AIAA J. 2014, 52, 72–81. [Google Scholar] [CrossRef]
Øye, S. Dynamic Stall Simulated as Time Lag of Separation. In Proceedings of the 4th IEA Symposium on the Aerodynamics of Wind Turbines, Rome, Italy, 20–21 November 1990. [Google Scholar]
Branlard, E. Wind Turbine Aerodynamics and Vorticity-Based Methods. In Fundamentals and Recent Applications. Research Topics in Wind Energy; Springer: Cham, Switzerland, 2017; Volume 7. [Google Scholar]
Hansen, M.O.L. Aerodynamics of Wind Turbines, 3rd ed.; Routledge: London, UK, 2015. [Google Scholar]
Smagorinsky, J. General Circulation Experiments with the Primitive Equations. Mon. Weather. Rev. 1963, 91, 99–164. [Google Scholar] [CrossRef]
Bailly, C.; Comte-Bellot, G. Turbulence; Springer: Cham, Switzerland, 2015. [Google Scholar]
Versteeg, H.K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics. In The Finite Volume Method, 2nd ed.; Pearson: London, UK, 2007. [Google Scholar]
Jiang, L.; Gu, X.; Wu, J. A Simplified Lattice Boltzmann Method for Turbulent Flow Simulation. Adv. Appl. Math. Mech. 2021, 14, 1040–1058. [Google Scholar] [CrossRef]
Balduzzi, F.; Bianchini, A.; Ferrara, G.; Ferrari, L. Dimensionless Numbers for the Assessment of Mesh and Timestep Requirements in CFD Simulations of Darrieus Wind Turbines. Energy 2016, 97, 246–261. [Google Scholar] [CrossRef]
Perzon, S. On Blockage Effects in Wind Tunnels—A CFD Study. In Proceedings of the SAE 2001 World Congress, Detroit, MI, USA, 5–8 March 2001. [Google Scholar]
Barlow, J.B.; Rae, W.H.; Pope, A. Low Speed Wind Tunnel Testing, 3rd ed.; John Wiley and Sons, Inc.: Hoboken, NJ, USA, 1999. [Google Scholar]
Rezaeiha, A.; Kalkman, I.; Blocken, B. CFD Simulation of a Vertical Axis Wind Turbine Operating at a Moderate Tip Speed Ratio: Guidelines for Minimum Domain Size and Azimuthal Increment. Renew. Energy 2017, 107, 373–385. [Google Scholar] [CrossRef]
Rogowski, K.; Hansen, M.O.L.; Maroński, R.; Lichota, P. Scale Adaptive Simulation Model for the Darrieus Wind Turbine. J. Phys. Conf. Ser. 2016, 753, 022050. [Google Scholar] [CrossRef]
Rezaeiha, A.; Montazeri, H.; Blocken, B. CFD Analysis of Dynamic Stall on Vertical Axis Wind Turbines Using Scale Adaptive Simulation (SAS): Comparison against URANS and hybrid RANS/LES. Energy Convers. Manag. 2019, 196, 1282–1298. [Google Scholar] [CrossRef]
Menter, F.R.; Egorov, Y. The Scale-Adaptive Simulation Method for Unsteady Turbulent Flow Predictions. Part 1: Theory and Model Description. Flow Turbul. Combust. 2010, 85, 113–138. [Google Scholar] [CrossRef]
Menter, F.; Hüppe, A.; Matyushenko, A.; Kolmogorov, D. An Overview of Hybrid RANS–LES Models Developed for Industrial CFD. Appl. Sci. 2021, 11, 2459. [Google Scholar] [CrossRef]
Rogowski, K.; Michna, J.; Ferreira, C.S. Numerical analysis of Aerodynamic Performance of a Fixed-Pitch Vertical Axis Wind Turbine Rotor. Adv. Sci. Technol. Res. J. 2024, 18, 97–109. [Google Scholar] [CrossRef] [PubMed]
Rogowski, K. CFD Computation of the H-Darrieus Wind Turbine—The Impact of the Rotating Shaft on the Rotor Performance. Energies 2019, 12, 2506. [Google Scholar] [CrossRef]
Carbó Molina, A.; De Troyer, T.; Massai, T.; Vergaerde, A.; Runacres, M.C.; Bartoli, G. Effect of Turbulence on the Performance of VAWTs: An Experimental Study in Two Different Wind Tunnels. J. Wind. Eng. Ind. Aerodyn. 2019, 193, 103969. [Google Scholar] [CrossRef]
Rogowski, K.; Królak, G.; Bangga, G. Numerical Study on the Aerodynamic Characteristics of the NACA 0018 Airfoil at Low Reynolds Number for Darrieus Wind Turbines Using the Transition SST Model. Processes 2021, 9, 477. [Google Scholar] [CrossRef]
ANSYS Inc. ANSYS Fluent Theory Guide Version 22; ANSYS Inc.: Canonsburg, PA, USA, 2022. [Google Scholar]
Michna, J.; Rogowski, K. Numerical Study of the Effect of the Reynolds Number and the Turbulence Intensity on the Performance of the NACA 0018 Airfoil at the Low Reynolds Number Regime. Processes 2022, 10, 1004. [Google Scholar] [CrossRef]

Figure 1. D2Q9 velocity set for a 2-D fluid-flow problem.

Figure 2. R500 case study [31].

Figure 3. Computational domain.

Figure 4. Sensitivity of the relative wind speed V_rel to the increasing mesh resolution Δh for simulations with ε equal to 0.5c.

Figure 5. Numerical model and mesh with boundary conditions.

Figure 6. Effect of the corrections to the normal force coefficient C_n* for ε = 0.5c (left-hand side) and ε = c (right-hand side). BASELINE = uncorrected curve; FC = flow curvature; DS = dynamic stall; T = Smagorinsky turbulence model.

Figure 7. Effect of the corrections to the tangential force coefficient C_t* for ε = 0.5c (left-hand side) and ε = c (right-hand side).

Figure 8. Effect of the corrections to the angle of attack α for ε = 0.5c (left-hand side) and ε = c (right-hand side).

Figure 9. Blade normal force coefficient C_n* vs. azimuthal angle ϑ (left-hand side) and tangential force coefficient C_t* (right-hand side) at a high tip-speed ratio (λ = 4.5). Comparison between PIV experimental data (EXP), Lattice Boltzmann Method (LBM), Scale-Adaptive Simulation (SAS), and Vortex Method (VM).

Figure 10. Blade normal force coefficient C_n* vs. azimuthal angle ϑ (left-hand side) and tangential force coefficient C_t* (right-hand side) at a low tip-speed ratio (λ = 2.0). Comparison between PIV experimental data (EXP), Lattice Boltzmann Method (LBM), Scale-Adaptive Simulation (SAS), and Vortex Method (VM).

Figure 11. Axial induction factor a_ϑ vs. azimuthal angle ϑ (left-hand side) and angle of attack α vs. azimuthal angle ϑ (right-hand side) for tip-speed ratios λ of 4.5 and 2.0.

Figure 12. Normalized streamwise velocity u_x_,n at different stations x/D downstream of the turbine for λ = 4.5 (left-hand side) and λ = 2.0 (right-hand side).

Figure 13. Contours of normalized out-of-plane vorticity for the combined field ω_zc/u₀ for λ = 4.5 (top) and λ = 2.0 (bottom).

Table 1. Velocity set in explicit form.

Velocities c_k	Length \|c_k\|	Weight w_k
(0,0)	0	4/9
(±1,0), (0,±1)	1	1/9
(±1,±1)	√2	1/36

Table 2. Rotor characteristics.

Number of Blades	2
Radius [m]	0.50
Blade chord [m]	0.06
Airfoil	NACA 0018

Table 3. Operating conditions.

λ	u₀ [l.u.]	Re
4.5	0.1	160,000
2.0	0.1	80,000

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Cacciali, L.; Hansen, M.O.L.; Rogowski, K. Highly Stable Lattice Boltzmann Method with a 2-D Actuator Line Model for Vertical Axis Wind Turbines. Energies 2024, 17, 4847. https://doi.org/10.3390/en17194847

AMA Style

Cacciali L, Hansen MOL, Rogowski K. Highly Stable Lattice Boltzmann Method with a 2-D Actuator Line Model for Vertical Axis Wind Turbines. Energies. 2024; 17(19):4847. https://doi.org/10.3390/en17194847

Chicago/Turabian Style

Cacciali, Luca, Martin O. L. Hansen, and Krzysztof Rogowski. 2024. "Highly Stable Lattice Boltzmann Method with a 2-D Actuator Line Model for Vertical Axis Wind Turbines" Energies 17, no. 19: 4847. https://doi.org/10.3390/en17194847

APA Style

Cacciali, L., Hansen, M. O. L., & Rogowski, K. (2024). Highly Stable Lattice Boltzmann Method with a 2-D Actuator Line Model for Vertical Axis Wind Turbines. Energies, 17(19), 4847. https://doi.org/10.3390/en17194847

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu