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.
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
Cn* and
Ct* plotted in
Section 3.2 and
Section 3.3 represent the normal and tangential force coefficients equated as
fn/(½
ρRu02) and
ft/(½
ρRu02), respectively, to align with experimental data of Tescione et al. [
31]. In these definitions,
fn and
ft denote the blade loads,
u0 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
R2 [
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,
R2 is calculated by measuring the deviation between the reference points
Fref and the numerical outputs
F to the deviation between
Fref and the average value of the observed data
Favg. It is expressed as:
The SAS solution was first averaged over the last ten revolutions to smooth the final curve, from which Favg was also calculated. This approach allowed the calculation of more accurate deviations Fref(ϑ) − 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].
ux1 = u0 and Uy1 = v0 were imposed at the inlet boundary. A second-order extrapolation is set for mass density such as ρ1 = 4/3 × ρ2 − 1/3 × ρ3 while pressure yields simply p1 = ρ1 × cs2. 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.,
ε = 1
c, 0.5
c, 0.25
c) 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 Vrel 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.5
c 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
3 and 1.79 × 10
−5 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
−6 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
−8 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
Cl and
Cd. 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 Cn* and Ct* 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 Cl and αv result in poor aerodynamic performance, often leading to counter-torque. Power extraction is absent where the tangential force Ct* is negative, which occurs when the virtual incidence αv inverts while the lift coefficient Cl remains positive.
At high rotor speeds, the effect of dynamic stall on the angle of attack is limited. It increases Cn* at the cusp upwind and reduces it slightly downwind. This trend also applies to Ct*, 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
Cl at negative angles of attack, as shown in [
40]. A positive torque is generated when the lift coefficient
Cl 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 Cn* and Ct* 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
Cn* of a kernel size
ε = 0.5c and the reference (SAS) is high (
R2 = 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
Cn* 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
Ct* 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
Ct* (
R2 = 0.937) and the reference.
However,
Ct* 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 Cn* and Ct* 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
Cn* and
Ct* 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
Cn* and 0.920 for
Ct*. 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 Cl(α) and drag coefficient Cd(α) 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 (ux/u0 and uy/u0), while the vorticity is normalized using the chord and free-stream velocity (ωzc/u0).
Despite being normalized to
u0, 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.