Comparison of Three-Dimensional Numerical Methods for Modeling of Strut Effect on the Performance of a Vertical Axis Wind Turbine

Abstract

This paper compares three different numerical models to evaluate their accuracy for predicting the performance of an H-rotor vertical-axis wind turbine (VAWT) considering the influence of struts. The strut of VAWTs is one factor that makes the flow feature around the turbine more complex and thus influences the rotor performance. The focus of this study is placed on analyzing how accurately three different numerical approaches are able to reproduce the force distribution and the resulting power, taking the strut effect into account. For the 12 kW straight-bladed VAWT, the blade force is simulated at three tip speed ratios by the full computational fluid dynamics (CFD) model based on the Reynolds-averaged Navier–Stokes (RANS) equations, the actuator line model (ALM), and the vortex model. The results show that all the models do not indicate a significant influence of the struts in the total force over one revolution at low tip speed ratio. However, at middle and high tip speed ratio, the RANS model reproduces the significant decrease of the total tangential force that is caused due to the strut. Additionally, the RANS and vortex models present a clear influence of the struts in the force distribution along the blade at all three tip speed ratios investigated. The prediction by the ALM does not show such distinctive features of the strut impact. The RANS model is superior to the other two models for predicting the power coefficient considering the strut effect, especially at high tip speed ratio.

1. Introduction

Renewable energy has attracted more and more attention in recent years, and wind power is one of the energy resources that can be a promising potential alternative to conventional electricity sources. The developed nations account for most of the world’s carbon dioxide emissions, and climate scientists warn that we must keep global warming below 1.5 degrees. There is no doubt that the capacity of clean energy, including wind power, will keep growing to cover global electricity demand in the future.

There are two main types of wind turbine configurations, which are the horizontal-axis wind turbine (HAWT) and the vertical-axis wind turbine (VAWT). Nowadays, HAWTs are taking the lead in the wind energy industry because of their higher energy efficiency. VAWTs have demonstrated an ability to fulfill certain energy generation requirements that cannot be achieved by HAWTs, for instance, under conditions such as high wind turbulence or wind fluctuations [1]. Another advantage of VAWTs is that the yaw control system is not needed due to their ability to take omni-directional wind. Additionally, a floating offshore VAWT does not need to be very tall to capture offshore wind. This leads to a lower support structure, which can mitigate the cost for operation and maintenance [2].

The aerodynamics of VAWTs are characterized by complex and unsteady behavior due to their inherent nature. For example, the relative wind speed seen by the blade is varied during rotation. One of the important phenomena for VAWTs is the dynamic stall, which can occur particularly at low tip speed ratios, where the blade suffers a large amplitude in the variation of angle of attack. During the dynamic stall, the vortex is formed and released from the leading edge, and the lift is reduced dramatically as soon as the vortex separates from the trailing edge. Moreover, the wake generated from the blades, struts, tower, and other components in the upwind region disturbs the flow in the downwind region. This reduces the torque downwind, which is lower compared to upwind.

Understanding the flow features of VAWTs is essential to predict and maximize the power efficiency. One of the numerical approaches employed widely today is the actuator line model (ALM). This model was developed by Sørensen and Shen [3,4] and was implemented by Mikkelsen and Troldborg [5,6]. In the ALM theory, the actual geometry of the rotor blade is replaced with lines carrying body forces corresponding to the loading of the blades in order to solve the Navier–Stokes equations. Thus, it is not necessary to build a very refined mesh of the rotor geometry for solving the boundary layer around surfaces [7]. Porté-Agel et al. [8] and Lu and Porté-Agel et al. [9] examined the characteristics of the wind turbine wakes in a wind farm considering the effect of the atmospheric boundary layer. The wakes of VAWTs with various parameters for the ALM, such as solidities, were examined by Hezaveh et al. [10].

Another computational tool commonly used for calculating the aerodynamics of VAWTs is the vortex model. In the vortex method, the Navier–Stokes equations are formulated in terms of vorticity to track the movement of vortices. Larsen [11] proposed a two-dimensional vortex model for VAWTs, and it was applied by Holme [12] and Wilson [13]. Strickland et al. [14] first used the three-dimensional vortex method for Darrieus turbine flows, and Cardona [15] modified the model developed by Strickland to reflect the phenomenon of flow curvature, which improved the prediction of the blade forces. Scheurich et al. [16] used Brown’s vorticity transport model to simulate the aerodynamic performance and the wake structure for the Darrieus-type VAWT. Zanon et al. [17] used the two-dimensional vortex panel method with the model for flow separation and vortex shedding to simulate the evolution of the wake from the VAWT.

Solving the Navier–Stokes equations directly by computational fluid dynamics (CFD) simulations is the common practice to analyze the aerodynamics for a wide range of engineering applications, and VAWTs are no exception. Ferreira [18] compared the CFD model with different turbulence models for a VAWT flow. Lanzafame et al. [19] developed a two-dimensional CFD model for the Darrieus turbine and underlined the important role of modeling for turbulence transition. Balduzzi et al. [20] carried out sensitivity analysis to identify the effective CFD setting, including the turbulence models and domain dimensions for an H-type Darrieus turbine. Many studies employ 2D models due to the high computational cost of reproducing the three-dimensionality of the flow around VAWTs, but 2D simulations could give an overestimated prediction, as shown by Howell et al. [21] and McLaren et al. [22]. The reason for this overestimation is that flow phenomena such as blade tip vortices and flow divergence are neglected in 2D simulations.

There are several advantages in the CFD approach, which has the capability of full modeling for the boundary layer flow on the blade surface with no approximation of the turbine geometry. The CFD method can handle any blade shape, while the ALM and the vortex method only work for airfoil profiles with known force coefficients. In addition, the lift and drag coefficients used in these models do not strictly represent the values that VAWTs actually produce during rotation, due to phenomena such as flow curvature or dynamic stall.

While the struts must be strong and stiff enough to prevent excessive deflections and carry the weight of blades and themselves [23], it is known that the struts cause the drag loss and consequently reduce the thrust of VAWTs [21,24]. For the reasons mentioned above, CFD simulations solving the wall boundary are suitable to capture the detailed flow characteristics of VAWTs. There are some studies, such as those by De Marco et al. [25], Elkhoury et al. [26], and Marsh et al. [27], that conducted CFD simulations with strut components considered in the models, and their results clearly indicated the influences of struts on the rotor performance.

This study compares the accuracy and limitations of three different 3D numerical approaches for predicting the aerodynamic characteristics of the VAWT with the strut effect taken into account. These three methods are the CFD model based on the Reynolds-averaged Navier–Stokes (RANS) equations, the ALM, and the vortex model. The presence of struts are not negligible and thus need to be modeled properly when the aerodynamic forces need to be predicted with high accuracy. The VAWT studied is a three-bladed H-rotor that has 6.5 m diameter and 6 m hub height. When the VAWT has struts mounted with an angle like this turbine, the struts generate lift forces and give a more prominent contribution to the rotor thrust. The blade forces of the turbine are simulated for two cases when the struts are and are not included in the models. The differences of the predicted forces are discussed for three tip speed ratios in order to investigate how appropriately each model reproduces the influence of struts.

2. Turbine Geometry and Test Conditions

Table 1. Parameters of the VAWT.

Property	Value
Rated power	12 kW
Number of blades	3
Diameter	6.5 m
Hub height	5.75 m
Blade length	5 m
Blade profile	NACA 0021
Chord length at left	25 cm
Blade pitch angle	2${}^{\circ }$
Strut profile	NACA 0025

lists the main parameters for the reference VAWT model. This turbine was designed by the Division of Electricity, Uppsala University. It is located at Marsta in Sweden, standing in open-air conditions over a flat surface. The radius of the rotor is $R=3.24$ m, and the hub has a height of 5.75 m. The rotor diameter is denoted by $D=2R$. The turbine consists of three blades, and they have the cross-section of a NACA 0021 profile. The length of the blades is 5 m. The chord length is $c=25$ cm. The blades are tapered linearly on both sides, the chord length being 1 m from the blade tip and 15 cm at the tip. As shown in Figure 1, two inclined struts are attached to the 27% position of the blade length at both ends. The cross-section of the struts is designed based on a NACA 0025 profile with modification for the blunt trailing edge. The chord length varies linearly from 32 cm at the root up to 20 cm at the attachment.

Figure 2 shows the top view of the turbine. The angular velocity of rotation and the azimuth angle are denoted by $\omega$ and $\theta$, respectively. $U_{\infty }$ represents the freestream velocity at hub height. The tip speed ratio is defined as $\lambda =R\omega /U_{\infty }$.

Table 2. List of the tip speed ratios λ studied in this paper.

λ	U∞ (m/s)	ω (rpm)
2.6	6.6	49.9
3.4	6.4	64.8
4.2	5.3	65.0

lists the three tip speed ratios studied in this paper. The normal and tangential forces are represented by $F_n$ and $F_t$, whose positive values are defined to be in the outward and rotational direction.

3. Numerical Methods

3.1. Rans Model

The open source code OpenFOAM [28], based on the finite volume method, is used to conduct incompressible flow simulations. Unsteady three-dimensional RANS simulations are performed with the realizable $k-\epsilon$ turbulence model [29]. DES or LES would also be appropriate to reproduce the flow field around the VAWT, and indeed, there are many studies using these models to analyze the aerodynamics of the VAWT. The purpose of this study is to calculate the time-averaged force; therefore, the RANS model will be adequate, considering its accuracy and computational cost.

The PIMPLE algorithm, which is a combination of the PISO and SIMPLE algorithm, is applied to couple the pressure-velocity equations. Every time step size is adapted to the maximum Courant number specified, and this value is set to be below 0.9 here. The average of the time step increment is $2\times {10}^{-6}$ s when the tip speed ratio is 3.4, and it is equivalent to the time of rotating by azimuth angle of $0.{0008}^{\circ }$. This number can be considered small enough according to the sensitivity analysis conducted by Rezaeiha et al. [30], who state that the azimuthal increment should be lower $0.1^{\circ }$.

The schematic of the computational domain and the boundary conditions are shown in Figure 3. The domain is divided into a rotational inner part and a stationary outer part. The rotational region, which is represented as the circular area in the figure, has 1.5D diameter. The rotor is located in the center of the domain. Both domain dimensions in the cross-stream and the vertical directions are 20.1D. The inlet velocity is represented by the log law. The log profile $U_{log}$ of a function of the height z is defined as

$$U_{log}=\frac{U^{*}}{K}ln\left(\frac{z-z_0}{z_0}\right)$$

(1)

where

$$U^{*}=\frac{K{U}_{\infty }}{ln\left(\frac{z_{hub}+z_0}{z_0}\right)}$$

(2)

is the frictional velocity, $K=0.41$ is the Kármán constant, $z_0=0.025$ m is the roughness length, and $z_{hub}$ is the hub height. At the inlet, the turbulence kinetic energy k is determined so that the turbulence intensity $I=\sqrt{2k/3}/U_{\infty }$ is approximately equal to 0.1%. The pressure at the outlet is set to be zero. The slip condition is assumed at boundaries in the cross-stream and the vertical directions, but the wall boundary conditions are applied at the bottom boundary $z=0$ to represent the ground.

Figure 4 shows the discretized mesh in the wake area, the rotor region, and the area around the blade surface at the hub height. A circular line in the middle picture represents the sliding interface boundary between the stationary and rotating parts. The height of the first grid point on the blade and strut surfaces is 2 mm, and the flows of the boundaries for these surfaces are modeled by the wall function. The mesh resolution satisfies the requirement for the wall functions, that is, $y^{+}<300$. The sensitivity of the mesh resolution and the domain size to the simulation results has been investigated previously [31]. The numbers of mesh cells are 19 and 33 million for the case containing only blades and the case including both blades and struts, respectively. Except for the strut surface region, the meshes between the two cases have the same resolution for the entire region.

First, the initial field is computed by solving the steady-state flow. Then, the simulations are run for more than 21 revolutions to develop the wake by a few rotor-diameter distances using a coarse mesh. The last revolution is solved using the finer mesh explained above.

The first-order upwind scheme is used to discretize the convective terms, and then the second-order upwind total variation diminishing (TVD) scheme is used for simulating the last revolution. The central difference scheme is applied for discretization of the diffusive terms. The bounded first-order implicit scheme is used for the time differencing.

The parallel computations are conducted on the Tetralith cluster, which is provided by the National Supercomputer Center at Linköping University. The calculations were performed on 256 processors using 8 nodes of the Intel Xeon Gold 6130 (produced by Intel^®, Mountain View, CA, USA), 2.1 GHz clock frequency computer. The computational domain is divided into 256 subdomains, and one processor is assigned to each subdomain. When the tip speed ratio is 2.6, it takes around 225 and 140 h to calculate one revolution for the cases with and without struts, respectively.

3.2. Actuator Line Model (ALM)

The actuator line model (ALM) is an unsteady three-dimensional model employed to investigate the flow fields around turbine rotors. It is based on the classic blade element model coupled with a solver for the governing equations of mass and momentum conservation (Navier–Stokes). This model splits the blades, struts (in the case of VAWT studies), and tower into line elements. These elements are modeled as a two-dimensional airfoil by using lift and drag coefficients $C_L$ and $C_D$, respectively. The dynamic stall model considered on the force coefficients was the Leishman–Beddoes model [32] with the modifications of Sheng et al. [33] and Dyachuk [34]. For this study, the ALM has been implemented by employing the library turbinesFoam developed by Bachant et al. [35,36,37]. Their filtered version of the Navier–Stokes equations and the incompressible LES approach were considered for simulating turbulence effects. This model has been previously validated and used for studying both HAWTs and VAWTs in a wide range of operational conditions [38,39].

First, the flow velocity is sampled to provide the relative velocity and angle of attack as a vector for each element. The local inflow velocity is calculated by averaging the velocity at axisymmetric points surrounding the quarter of the chord at a radial distance. In this code, the velocity is sampled at 20 points and a radial distance of $2ϵ$, where $ϵ$ is the parameter defined below. Later, the dynamic stall model calculates the unsteady blade force coefficients, imparting them back into the solver as body forces. The lift and drag forces per spanwise length unit, $f_l$ and $f_D$, are calculated by

$$f_L=\frac{1}{2}\rho c{C}_L{\left|V_{\mathrm{rel}}\right|}^2$$

(3)

and

$$f_D=\frac{1}{2}\rho c{C}_D{\left|V_{\mathrm{rel}}\right|}^2$$

(4)

respectively, where $\rho$ is the flow density.

The coefficients $C_L$ and $C_D$ are obtained by a linear interpolation from a table for a specific $\alpha$, and then the body forces are resolved by using the interpolated coefficients together with the blade element method. The forces in the lines of elements are added into the momentum conservation equation as a source of body force per unit of density. Note that the Reynolds number varies during revolution, and the ALM as well as the vortex method includes this Reynolds number dependence in the coefficients.

In this paper, the results obtained using two different coefficient are presented, namely those from Sheldahl and Klimas [40] and from the XFOIL program [41]. These two versions of the ALM are denoted as ALM-SK and ALM-XF, respectively.

3.2.1. Calculated Force Distribution

The forces in the solver are smoothly distributed over mesh cells to avoid numerical instabilities caused by high gradients. The source term of the forces is projected onto the element location from its maximum value using a three-dimensional Gaussian kernel.

The implemented smoothing function ${\eta }_F$, which is multiplied by the computed local force on the actuator line element and dispensed on a cell with a distance $|\overrightarrow{d}|$ from the quarter chord location of the actuator line element, is expressed as

$${\eta }_F=\frac{1}{{ϵ}^3{\pi }^{3/2}}exp\left[-{\left(\frac{|\overrightarrow{d}|}{ϵ}\right)}^2\right]$$

(5)

with $ϵ$ representing the smoothing width parameter. This parameter is determined from the maximum among three values: a quarter of the chord length, the mesh size, and the momentum thickness due to drag forces, that is,

$$ϵ=max\left[\frac{c}{4},4\sqrt[3]{V_{\mathrm{cell}}},\frac{c{C}_D}{2}\right]$$

(6)

where $V_{\mathrm{cell}}$ is the volume of the cell.

3.2.2. End Effects

The lift force distribution must drop to zero at the blade tips for the finite span, as a consequence from Helmholtz’s second vortex theorem. From Prandtl’s lifting theory, the geometric angle of attack $\alpha$ can be expressed as a function of the adimensional location $\theta$ [42]:

$$\alpha \left(\theta \right)=\frac{2S}{\pi c\left(\theta \right)}\sum _{n=1}^N{A}_{n}sin\theta +\sum _{n=1}^{N}n{A}_n\frac{sinn\theta }{sin\theta }+{\alpha }_{L=0}\left(\theta \right)$$

(7)

where S is the total span length, $c\left(\theta \right)$ is the chord length as function of the spanwise location, and N represents the sampling number of elements. To solve the Fourier coefficients $A_n$, this equation is expressed in the matrix form

$$\left[{\alpha }_m\right]-{\alpha }_{L=0}=\left[D_{mn}\right]\left[A_n\right]$$

(8)

where

$$D_{mn}=\left[\frac{2b}{\pi c_m}sinn{\theta }_m+n\frac{sinn{\theta }_m}{sin{\theta }_m}\right].$$

(9)

With the calculated $A_n$, the circulation distribution can be obtained as

$$\mathrm{\Gamma }\left(\theta \right)=2S{V}_{\infty }\sum _{n=1}^N{A}_{n}sinn\theta$$

(10)

which gives the correction factor of the distributed lift coefficient

$$C_L\left(\theta \right)=\frac{-\mathrm{\Gamma }\left(\theta \right)}{\frac{1}{2}c{V}_{\infty }}$$

(11)

through the Kutta–Joukowski theorem. $C_L\left(\theta \right)$ varies between 0 and 1, the former representing the value at the tips location.

3.2.3. Computational Settings

The number of total cells is 17 million for all cases. The cell sizes in the rotor region and in the wake are approximately 0.11 m and 0.21 m, respectively. The entire domain size, the distances of the rotational center to each boundary, and the wind condition at inlet (such as the log law profile and the turbulence level) are the same as those in the RANS model. The numbers of discretization points along the blade and strut are both 30.

The time step size is approximately 0.0048 s, which corresponds to azimuth angle of $1.9^{\circ }$ per step. The computation is run on the same cluster computers as what are used in the RANS model. For the cases with and without struts, around 40 and 50 revolutions are calculated, and it takes around 17 and 14 h to complete the simulations by using 128 processors, respectively.

3.3. Vortex Method

The vortex method is based on using the vorticity as the discretization variable, instead of the velocity. The current method falls into the category of free vortex methods, where the vortex elements propagate with the flow, making the method Lagrangian. The velocity field can be constructed from the vorticity field using Biot Savart’s law. Here, one assumes an infinitely large simulation domain where the vorticity is zero in all regions outside of the wake of the turbine, and vorticity only exists in the turbine wake. The advantage is that the method automatically places vorticity in the regions of interest, and one does not have to create calculation particles in the other regions. The disadvantage with the current implementation is that the method only works for constant inflows, and windshear and turbulence are therefore not included.

The vortex method implemented here is a vortex filament method. The vortices are modeled as a lattice of straight lines, connecting at the grid points. To propagate the vortices, vortex filaments are created between each grid point in the lattice, and the flow velocity is calculated at the grid points, which are propagated with the flow velocity. See Figure 5 for illustration. To avoid the singular behavior at the center of the vortex, the vorticity at each point along the filament was distributed evenly over a sphere with radius given by $1.2max(R\mathrm{\Omega }\Delta t,d_{neighbor})$, where R is the turbine radius, $\mathrm{\Omega }$ is the rotational speed, $\Delta t$ is the time-step, and $d_{neighbor}$ is the maximum distance to any of the neighbor points in the vortex lattice.

The advantage of the filament method is that it automatically fulfills the divergence-free condition of the vorticity field, which removes the issues with vortex stretching and divergence relaxation of the particle methods. The disadvantage is that there is no clear path to implement diffusion or vortex merging within the filament approximation. This makes the model only suitable to model the near-wake region where the flow is relatively ordered, but the method has issues in the far-wake region where the wake decay causes the flow to become turbulent, which causes extensive stretching of the filaments.

The vortex method can easily be combined with panel methods to solve the potential flow problem. This is useful for the modeling of the blades. Here, a simplified approach where the blades are modeled as flat plates is used. The blades are then modeled using vortex rectangular vortex rings, where each vortex ring is comprised of four vortex filaments. The no-penetration boundary condition can be fulfilled by setting the relative normal velocity to zero at the center of each vortex ring. The circulation around the blades can be calculated by fulfilling the Kutta condition that the flow should be parallel to the flat plate at the trailing edge. This is accomplished by using one extra control point behind the trailing edge, with the same condition that the relative velocity normal to the blade should be zero here as well. To account for viscous effects and stall, the method was adjusted according to the following procedure: First, the relationship between the circulation and the angle of attack was determined with static 2D simulations of a single blade. Second, the following procedure was used for each time step to calculate lift and drag forces:

1.

The circulation around each blade was calculated according the method described above. The blade was divided into multiple segments along its span and one circulation was calculated for each segment.

2.

The correlation between circulation and angle of attack previously determined was used to obtain a corresponding angle of attack for each blade segment.

3.

An external force model was used to calculate the lift and drag forces for the blade. This model is the same as the one used for the ALM-SK simulations and uses static airfoil data combined with a dynamic stall model, as described in Dyachuk and Goude [43].

4.

Calculate the new circulation from Kutta Joukowski’s lift formula $f_L=-\rho V_{rel}\mathrm{\Gamma }$, where $f_L$ is the lift force per unit span, $\rho$ is the density, $V_{rel}$ is the relative velocity (which is calculated from the vorticity field), and $\mathrm{\Gamma }$ is the new circulation.

5.

The circulation that is to be added to the flow can now be determined as the change in circulation from the previous step.

As a simplification, to account for the fact that the blade circulation after this correction procedure no longer fulfills the Kutta condition, it was chosen to reduce the blades to lifting lines during the propagation step for the vortex filaments. This approximation is validated in Dyachuk et al. [44] and works because it is the calculation of the angle of attack that is most important for the model to achieve correct results. Validations of the vortex method can be found in Mendoza and Goude [39].

The vortex method solved all blades and all struts at the same time for the potential flow calculation. The position of the panels were adjusted to ensure that the connection between the struts and the blades are always at the edge of a panel for both the blades and the struts. This gave a slightly different number of blade segments for the two cases: the case without struts were simulated using 60 blade segments, while the case with struts used 57 blade segments and 30 strut segments.

Each simulation was run for 25 revolutions with 120 time steps per revolution, which corresponds to azimuth angle of $3^{\circ }$ per step. The calculations are performed on 12 CPU cores of the Intel Xeon E5645 processors (produced by Intel^®, CA, USA) at 2.4 GHz clock speed and the GPU of the GeForce 480 RTX. The computational time is approximately 67 and 12 h for the cases with and without struts, respectively.

4. Results and Discussion

Figure 6 shows the instantaneous velocity and vorticity fields in the vertical plane intersecting the center of the rotor for $\lambda =3.4$, which are calculated by the RANS model, the ALM, and the vortex model for both cases with and without struts. To indicate the swept area, the plots include an illustration of the blades and struts located at $\theta ={90}^{\circ }$ and ${270}^{\circ }$ that are not actually in these positions in the simulation. They show the magnitude of the flow velocity normalized with $U_{\infty }$ and the y component of the non-dimensionalized vorticity, which is calculated by $(\nabla \times U)c/U_{\infty }$. Figure 7 shows the velocity and vorticity fields in the horizontal plane at $z/D=0.97$, which is the middle between the upper strut attachment point and the hub height. The velocity and vorticity illustrated are the values in the y and z directions, respectively. Each of the three blades in Figure 6 and Figure 7 is located at the same azimuth angle for all plots, with one of the blades passing at $\theta ={60}^{\circ }$.

It can been seen from all plots in Figure 6 that the vortices are shed from the blade tips and travel downstream. The pictures from the RANS model represent the influence of the strut clearly in both regions within the rotor and the wake where the velocity is dramatically dropped and the strength of vortices is high over a large area. The vortex model also displays large differences between two cases. The vortex model represents well the complex vortical structure caused by the struts both within the rotor area and in the wide range of the wake, while the other models almost only reproduce clear vortices released from the blade tips. These pictures do not show as noticeable of an influence of the struts for both the ALM-SK and ALM-XF models as for the RANS and vortex models, but the reduction of the velocity in the wake can be slightly observed in the ALM.

Additionally, the pictures indicate that the wind shear effect is reproduced in the RANS model and the ALM. As a result, the vortices shed from the upper tips convect faster than those from the lower. Note that the vortex model does not consider the wind shear.

The above findings are also confirmed from Figure 7, that is, the RANS and vortex models show more pronounced differences between the two cases than the ALMs. It is observed from the pictures of the RANS model that the vortices at the edge of the wake travel downstream, keeping their strength when the struts are not included. When the struts are present, the strength of the vortices decays rapidly, but instead, they spread out to a wide area in the wake. The behavior of this vortex diffusion for the case with struts is identified not only at this height $z/D=0.97$, but also at almost all blade heights that are not presented in the paper. The pictures from the RANS and vortex models clearly show vortices generated from the struts as well as those from the blades, which can affect the flow around the blades rotating in the downwind side. Both the ALM-SK and ALM-XF produce strong vortices that convect downstream, almost keeping their pattern. However, there is no remarkable differences caused due to the struts. The fact that the wake predicted by the ALM starts to break up further downstream compared to the RANS model coincides with the results observed by Troldborg et al. [45]. They explain, as the reason for this wake behavior, that the fully resolved model is able to simulate fluctuations of the blade force where the stall occurs, and this makes the wake more unstable. There is another study by Kalvig et al. [46], who simulated the wake velocity for the HAWT with different numerical models. They present the differences observed in the predicted wake velocity profiles, as well as the blade force. The results predicted by the fully resolved RANS method match better to the measurements than those by the actuator disk and line methods.

Figure 8 shows the instantaneous velocity profile in x direction along the vertical line, $0\le z/D\le 1.8$, at $y=0$ for three downstream locations, $x/D=1$, 2, and 5. The horizontal dotted lines indicate the height of the upper and lower blade tips. The strut attachment points are located around the height of $z/D=0.7$ and 1. Figure 9 shows the velocity profile in x direction along the horizontal line, $-1\le y/D\le 1$, at $z/D=0.97$. The presented data are those for the tip speed ratio of 3.4.

The RANS model shows significant differences due to the struts in the area close to the rotor, as seen in the plots of both the horizontal and vertical velocity profiles. The vertical velocity profiles at $x/D=1$ and 2 show two sharp decreases right at the strut attachment height. The plots from the two ALMs are almost identical between the cases with and without struts in the horizontal velocity profiles, but the vertical velocity curves show larger velocity deficit in the wake when the struts are included. It is also observed that the RANS model and the ALM predict close wake width. In addition, some profiles, for example the plot of the ALM-SK at $x/D=1$, show the wake asymmetry. In other words, the maximum velocity deficit is located not at $y=0$, but slightly at the positive side of y. The velocity profiles of the vortex model represent intensive fluctuations when the struts are included.

For all three models, the normal and tangential forces presented in this paper are all defined from the quarter chord position to make the tangential force the force that generates the torque. To verify that the torque can be calculated using the tangential force defined above, Figure 10 is plotted, where the torque and the force multiplied by the radius during one revolution are compared. The curves presented are those for the case with struts at $\lambda =3.4$ simulated using the RANS model. The values shown are of one blade and two struts, normalized by $C_T=T/0.5\rho rA{U}_{\infty }^2$, where T is the torque and A is the swept area of the rotor. In the OpenFOAM code, the torque is computed by multiplying the force with the distance between the rotor axis and the position of each cell where the force is generated. Both curves have close values overall. There are some small differences, and it is considered that these differences are caused due to the assumption of constant radius over all the surfaces. However, as these differences are minor, the performance of the turbine can be evaluated by checking the tangential force.

Figure 11 shows the comparison of the normal force over one revolution obtained using the three models for two cases with and without struts. The plots in the top, middle, and bottom rows represent the forces for the tip speed ratios $\lambda$ of 2.6, 3.4, and 4.2, respectively. The plots also include the forces measured by Dyachuk et al. [47] for this VAWT. The measurement was conducted using load cells for detecting tension and compression, which were installed to one of the blades between the struts and the hub. The data bins measured when the wind speed can be considered steady were extracted to be used for calculation. The data presented here are those obtained by averaging around 8 revolutions. A more detailed description of the experimental method is explained in reference [47]. To assess the reliability of the numerical models, the differences between the measured and simulated normal forces are presented in Appendix A.

The results show that two curves with and without struts cases are quite similar for each of the models at all tip speed ratios, which indicates that any of the models represent no significant impact of struts. When compared with measurement, the amplitude is overpredicted by all three models in the upwind side, especially for $\lambda$ of 2.6 and 3.4. The prediction agrees well in the upwind region for $\lambda$ of 4.2, but the steep increase up to 225${}^{\circ }$ is not reproduced by any of the models. There is a drop of the measured force at around 270${}^{\circ }$ for $\lambda$ of 3.4 and 4.2, which is also not present in the prediction at all.

Figure 12 shows the comparison of the tangential force over one revolution obtained using all the models. It can be seen that the RANS model clearly reproduces the strut influence, since the amplitude of the force is decreased in the downwind side at middle and high tip speed ratio, that is, $\lambda$ of 3.4 and 4.2. It decreases less at $\lambda$ of 2.6 because the interaction between the blade and the wake from the struts is more critical for higher tip speed ratio. In addition, the negative peak at around $\theta ={20}^{\circ }$ obtained by the RANS model has larger amplitude when the strut is considered for all tip speed ratios. This means that the strut affects even in the upwind side. It can also be found that more significant differences among all the models can be observed in the tangential force than in the normal force. These discrepancies are more pronounced for smaller tip speed ratios. This can be considered to happen due to the dynamic stall, which is likely to occur at low tip speed ratio. When $\lambda$ is 2.6, the force predicted by the RANS model fluctuates at $\theta ={213}^{\circ }$ for both cases with and without struts. A further investigation of the flow fields in the RANS simulation has revealed that the flow is completely separated and the vortex is shed from the blade surface at this azimuth angle, which indicates that the RANS model is able to capture the dynamic stall. The ALM-SK also represents the large drop of the tangential force at $\theta ={223}^{\circ }$, while the amplitude of the force predicted by the ALM-XF remains high even in the downwind side.

Figure 13 shows the comparison of the normal force distribution $F_n^{*}$ predicted by the three models for all the tip speed ratios when the blade is located at $\theta ={90}^{\circ }$ and ${270}^{\circ }$. $F_n^{*}$ is the forces per unit length and is calculated by $F_n^{*}=f_n/l$, where $f_n$ is the normal force of each segment, and l is the segment length. Note that in the RANS model, the blade is split into 25 segments, and the force of each segment is recorded during simulations. The values of $F_n^{*}$ at $\theta ={90}^{\circ }$ and ${270}^{\circ }$ are negative and positive, respectively, so these values are plotted on the left and right sides in one plot. The plots in the left, middle, and right columns illustrate the forces for the tip speed ratios $\lambda$ of 2.6, 3.4, and 4.2, and the plots in the top and bottom rows show the forces simulated without and with struts. The strut attachment points are at the normalized height of 0.27 and 0.73.

When the struts are not considered, the amplitude of the forces predicted by the three models is close at $\lambda$ of 3.4 and 4.2. The shapes of the plotted lines are similar, being flat at the center and becoming small toward both ends. However, there are some differences seen at $\lambda =2.6$. The ALM-XF overpredicts the force at both azimuth angles, as the flow predicted is not properly stalled. The ALM-XFs uses the blade force data calculated based on solving the boundary layer equations, and these data may not be reliable for deep stall. The ALM-SK and the vortex models use the same force model and are therefore close each other.

Although there is no significant influence of struts in the total forces observed in Figure 11, Figure 13 indicates that the struts strongly affect the force distribution at all tip speed ratios. This influence is especially significant in the downwind side, $\theta ={270}^{\circ }$. It is obvious from the RANS and vortex models at all three tip speed ratios that the curves have irregular profile shapes around the strut attachment points, which means that the models are able to reflect the influence from the struts. The influence of struts cannot be seen so noticeably in the ALM at all tip speed ratios compared to the other models. It is because the struts are considered as independent wings with their own tip vortex at the end, while the other models consider them as a coupled system.

Figure 14 shows the comparison of the tangential force distribution $F_t^{*}$, which is calculated as $F_t^{*}=f_t/l$, predicted by the three models when the blade is located at $\theta ={90}^{\circ }$ and ${270}^{\circ }$, respectively. These figures represent the results for all the three tip speed ratios and for the cases without and with struts, as presented in Figure 13.

Figure 14 confirms the findings mentioned in Figure 13 that the influence of struts are simulated properly by the RANS and vortex models. It is seen in Figure 14 that the strut affects the force distribution significantly, even in the upwind side. It is also observed that the pattern of the force distribution caused by the strut looks complicated when the blade moves in the downwind side. The fluctuation of the force at each height is large, especially in the curves simulated by the RANS model at all three tip speed ratios. The ALM does not produce a pronounced difference caused by the strut in the tangential force distribution or in the normal force distribution. At $\lambda =2.6$, the tangential force predicted by the ALM-XF is much higher than those by other models, especially in the downwind side.

As stated earlier, the ALM predicts the influence of struts less than the RANS and vortex models do. The ALM assumes that struts are individual blades, which each have a tip vortex assigned to them. The model will push the force down to zero at the tip of the strut due to the tip vortex corrections used, as it otherwise can give unfeasible forces at the blade tips. The main issue appears because the ALM does not fulfill the no-penetration boundary condition on the blades, which one should do to get a model that can resolve the tip effect by itself without the need for tip correction models. Without this modeling, the ALM will not capture the effect that the bound circulation can propagate from the struts into the blades, similar to what happens when a winglet is added to a blade. As it is the propagation of bound circulation from the struts to the blades that is causing the jump in forces on the blade at the attachment points of the struts, this jump will not be modeled by the ALM.

Figure 15 shows the power coefficient $C_p$ calculated by all the models for both the cases without and with struts. $C_p$ is defined as $\omega T_{ave}/0.5\rho A{U}_{\infty }^3$, where $T_{ave}$ is the average torque of three blades during one revolution. This result shows that the RANS model predicts the change of the $C_p$ due to the strut most significantly among all the models. The difference of $C_p$ is largest at $\lambda =4.2$ for all the three models, and these decreases of $C_p$ by considering struts are 24%, 13%, 14%, and 9% for the RANS model, the ALM-SK, the ALM-XF, and the vortex model, respectively. Only for the vortex model, $C_p$ at $\lambda =2.6$ is a little higher when considering struts than when not. Additionally, it is noted that $C_p$ values predicted by these models are not close to each other at any of the tip speed ratios. At $\lambda =4.2$, the RANS model considering struts gives the lowest $C_p$ among all. The RANS model represents the decrease of the tangential force caused by struts in both the upwind and downwind sides, which results in the low values of $C_p$. Due to the low tangential force calculated in the upwind side, $C_p$ is underestimated by the ALM-SK at $\lambda =2.6$. The ALM-XF show little difference of $C_p$ among three tip speed ratios compared to the other models. Although the vortex model reproduces the clear strut influence in the force distribution, the difference of $C_p$ by considering struts is smallest among all the models at all three tip speed ratios.

5. Conclusions

This paper compares the performance of three different numerical methods for predicting the aerodynamic force and the power of the VAWT with consideration of the strut influence. Since the strut of VAWTs can affect the flow fields and consequently the rotor power, this study aims to investigate how precisely the numerical models are able to reproduce the influence of struts on the turbine aerodynamics. The 12 kW H-rotor VAWT has been studied for the tip speed ratios of 2.6, 3.4, and 4.2, and the blade force is simulated using the RANS model, the ALM, and the vortex model. The results at low tip speed ratio show that all three methods do not indicate significant differences of the total forces during one revolution that are caused by the presence of struts. However, the RANS model represents clear influence of the strut in the tangential force at middle and high tip speed ratio, especially when the blade moves in the downwind side. For the force distribution along the blade height, both the RANS and vortex models show the significant difference when the strut is considered at all three tip speed ratios. The ALM predicts the slight reduction of the velocity in the wake due to the struts but does not reproduce the influence of struts on the force distribution as clearly as the other models do. Of particular importance is that the RANS model is able to simulate the large influence of the strut in the prediction of the power coefficient. The influence is significant, especially at high tip speed ratio, and the decrease of the power coefficient due to the strut is predicted as 24% by the RANS model, while those of the other models are less than 14%. The vortex model shows larger differences due to struts in the force distribution than the ALM, but the difference of the power coefficient by considering struts is slightly higher for the ALM than the vortex model. These results show the necessity of properly selecting the numerical model, depending on how high an accuracy is required for prediction of the performance of the turbine.