Unsteady effects on a pitching airfoil at conditions relevant for large vertical axis wind turbines

Interest in offshore vertical axis wind turbines (VAWTs) has created a need to study VAWTs at much higher Reynolds numbers than they have previously been studied at. VAWTs are characterised by unsteady aerodynamics, and high Reynolds numbers have the potential to alter the blade aerodynamics significantly. Here, results are reported on an airfoil that is pitched sinusoidally around zero angle of attack at Rec = 1 × 106, at a reduced frequency k = 0.15 and amplitudes of 5° < α < 20°. Since the static stall angle is not exceeded, no stall effects occur. Nevertheless, lift, drag and moment coefficients show noticeable hysteresis. Despite the high amplitudes of oscillation, lift and moment coefficients show reasonable agreement with unsteady aerodynamic theories by Theodorsen and Motta et al. with regard to the width of the hysteresis loops, but have noticeably steeper slopes.


Introduction
Vertical axis wind turbines (VAWTs) have several advantages compared to conventional horizontal axis wind turbines (HAWTs) [1]. They are agnostic to wind direction, have a simpler mechanical design, and have the potential to interact beneficially with each other when arranged in farms [2]. However, the aerodynamics of VAWTs are more complex than those of HAWTs, and cannot be approximated by steady one-dimensional models for optimisation and performance predictions. The non-dimensional parameters that characterise the aerodynamics of a VAWT are the Reynolds number, the tip speed ratio λ, and the solidity σ: Here, ν is the kinematic viscosity, U ∞ is the freestream velocity, D is the turbine diameter, ω is its angular velocity, c is the blade chord length, and n is the number of blades. The Mach number is assumed to be M a ≤ 0.3, so that compressibility effects are negligible. Most VAWTs deployed to date are significantly smaller than commercial HAWTs. Modern offshore HAWTs have diameters up to 220 m and produce up to 12 MW [3]. In contrast, the largest VAWT deployed to date was 64 m in diameter [4], while most are much smaller. Nevertheless, the literature suggests that VAWTs have certain advantages over HAWTs particularly with regard to offshore applications [5]. Offshore VAWTs would be considerably larger and experience higher wind speeds than previously studied VAWTs, leading to significantly higher Reynolds numbers. While most existing VAWTs operate at Re D ≤ 3 × 10 6 , offshore VAWTs might operate at Re D ≥ 10 7 . This creates the need for detailed study of Reynolds number effects on the aerodynamic behaviour of VAWTs. In a recent experimental study, Miller et al. (2018) showed Reynolds number effects on the power coefficient up to Re D = 4 × 10 6 [6].
Studies of VAWT aerodynamics at the blade level can provide insight into such Reynolds number effects. The parameters that characterise an individual turbine blade are the Reynolds number based on the chord length and mean effective velocity, the reduced frequency k, the angle of attack α, the velocity ratio u rel , and the aspect ratio AR: Here, U r is the relative velocity seen by the blade, f is the pitching frequency of the airfoil, and s is the airfoil span. Since U r changes throughout a rotation, Re c and k are here defined as averages based on the mean relative velocity U r following [7]. This blade-based parameterisation can be related to the turbine-based parameterisation as follows: The reduced frequency, k, which is a measure of the unsteadiness of the flow over the blade, is therefore solely a function of the turbine geometry. The angle of attack, α, and the velocity ratio, u rel , are functions of λ and the azimuthal position of the blade θ. The exact relationship depends on the deceleration of the flow ahead of the turbine. In BEM-type models, this effect is captured by the local induction factor. There is no analytical expression for the velocity of the flow at the turbine, and thus no precise analytical expression for α and u rel . However, they are often approximated by the following equations, in which the induction factor is assumed to be zero: Thus, the amplitudes of the angle of attack oscillationα and of the velocity ratio u rel are solely a function of λ. Unsteady effects are known to occur in both attached and separated flow conditions. The most prevalent unsteady effect is dynamic stall, which occurs when the angle of attack rapidly exceeds the stall angle. This sudden increase in angle of attack leads to the formation of a leading edge vortex which then convects downstream, causing a temporary increase in the lift force before it drastically collapses [8]. This phenomenon affects not only the performance, but also the structural integrity of the turbine.
Even if the blades do not exceed the stall angle at any point during their trajectory, unsteady phenomena still occur. Theodorsen developed analytical expressions for the time-dependent lift and moment coefficients of a flat plate pitching and heaving at small amplitude in a potential flow [9]. They account for added mass and shedding of vorticity into the wake. For the purpose of this study, pure pitching motion will be considered. Following [10], the equations for the lift and quarter-chord moment coefficients of an airfoil pitching sinusoidally around its half-chord are as follows: C L,th = π (2 C(k) sin(ωt) + k (C(k) + 1) cos(ωt)) (2) Here, C(k) is the lift deficiency factor. The equations given are exact only for a flat plate oscillating at small angles of attack in an inviscid flow. Nevertheless, they have shown reasonable  [12]. Motta et al. studied the effect of airfoil thickness on the unsteady lift and moment coefficients of a pitching airfoil and find that below an inversion point in the reduced frequency, an increase in thickness leads to an increase in hysteresis, whereas above the inversion point the trend is reversed [13]. The inversion point is itself a function of airfoil thickness. They modified the Theodorsen model to include an empirical correction for arbitrary airfoil thickness which modifies equations 2 and 3 as follows: Expressions for the P L(c), P M (c), and C(k) are given in [13]. Most experimental studies of unsteady lift and moments on airfoils in attached flow are limited to amplitudes of oscillation belowα = 10 • . Both numerical simulations and analytical models can be applied to capture the aerodynamic behaviour and performance of VAWTs. However, they often require tabulated airfoil data as inputs. Large Eddy Simulations (LES) can provide detailed accounts of the flow field, but are computationally intensive, especially at high Reynolds numbers [14]. To reduce the computational power required, LES is often combined with a blade level model that requires a lookup table. The most commonly applied analytical model is the double multiple streamtube model, a BEM-type model that uses momentum theory and an aerodynamic force analysis to iteratively solve for spatially-resolved induction factors, with significantly less computational effort [15]. The tabulated airfoil data are generally taken from steady airfoil experiments in which Re c , the airfoil section and its angle of attack are matched, e.g. [16]. Such data do not account for unsteady effects on the blades due to the oscillating angle of attack and inflow velocity. Instead, corrections are often applied to account for dynamic stall effects at high amplitudes of oscillation [17].
Under what conditions a VAWT experiences dynamic stall as opposed to attached flow unsteady effects is the subject of ongoing investigations. The static stall angle α ss depends on Re c , the blade profile, and its aspect ratio AR, while the amplitude of oscillationα depends on λ and the induction factor a. The occurrence of stall effects therefore likely depends on all of these variables. They are presumably less likely to occur at high α ss and lowα, i.e. at high Re c , on thick blades with low AR, at high λ, and high a. Ferreira et al. state that dynamic stall generally occurs for λ ≤ 5 [18], whereas Buchner et al. find that at 7 × 10 4 ≤ Re c ≤ 14 × 10 4 no leading edge separation occurs for λ ≥ 3.5 [7]. With regard to the design of offshore VAWT, the nature of unsteady effects at high Re c in particular requires further study.
The present study captures the unsteady effects on a VAWT blade at high Reynolds number. Experiments were conducted in a pressurised wind tunnel, enabling high Reynolds numbers without the need for high velocities or large models. A pitching airfoil with a NACA 0021 profile was studied at Re c = 1 × 10 6 and k = 0.15. The freestream velocity was held constant over the cycle. The airfoil was pitched at varying amplitudes around a mean angle of α 0 = 0 • . While the angle of attack variation of a VAWT blade is given by equation 1, it was approximated as sinusoidal following [19]. In all cases, the amplitude of oscillation was below the static stall angle of the airfoil at this Reynolds number, and corresponded to tip speed ratios of λ ≥ 2.9. The attached flow hysteresis in lift and moment coefficients were compared to the models by Theodorsen [9] and Motta et al. [13].

Methodology
Experimental data were obtained in the High Reynolds number Test Facility (HRTF) at Princeton University, a closed-loop wind tunnel that can achieve pressures up to 230 atm and The airfoil is equipped with 32 pressure taps, which provide the pressure distribution around the surface of the airfoil. A force balance records lift, drag and moments on the airfoil. Generally, the force balance provided higher quality measurements than the pressure transducers at the low pressures and velocities used in this study. Thus, all data presented here were calculated using the force balance measurements. The sampling frequency was 20 kHz in all cases. Phaseaveraged quantities were computed from all available cycles. The number of cycles for each run, as well as the number of data points per cycle are given in table 1. The accuracies of the pressure, density, and viscosity measurements are given in [6], and correspond to a systematic error of ±0.97% in Re c and ±0.4% in k. Slight variations in Reynolds number up to ±7.3% occurred throughout an experimental run due to slight changes in temperature and pressure inside the tunnel. The force balance used was a JR3 multi-axis load cell with a manufacturer-provided range of 200 N with an accuracy of ±0.5 N for the forces, and a range of ±25 Nm with an accuracy of 0.0625 Nm for the moments. Because the experimental setup was designed for much higher Reynolds numbers, the forces recorded by the force balance were small compared to its range. However, zeroes were subtracted to account for any systematic offsets, and phase averaging was applied to reduce random fluctuation, so that the true error is likely much smaller than the manufacturer-provided accuracies suggest. While some scatter  in the data is expected due to variations in the flow between cycles, the predominant source of random error in this setup is electrical noise in the measuring equipment. Electrical noise is not expected to affect the phase-averaged quantities as long as they are fully converged. All data were filtered before phase averaging, and subsequently smoothed to achieve better convergence. However, the measurements of C d and C m atα = 5 • appear to not be fully converged. The exact uncertainty in the experimental setup is currently unknown and will be the subject of future investigation. The drag due to the supporting rod and the end plates was not subtracted, so that the drag coefficients reported should be considered to have an unknown offset. However, hysteresis behavior should be qualitatively correct.
The airfoil was oscillated sinusoidally around a mean angle of attack, α 0 , using a stepper motor. The reduced frequency was held constant across all experimental runs at k = 0.15. This corresponds to a VAWT chord-to-diameter ratio of approximately 1:7. Equation 1 suggests that both the angle of attack and the relative velocity of a VAWT blade vary quite drastically at low tip speed ratios. For example, at λ = 2, the angle of attack amplitude isα = 30 • , and the amplitude of velocity oscillation isÛ r = 0.447 U r . However, these calculations do not account for variations in induction factor. A non-zero induction factor reducesα corresponding to any given λ, or conversely reduces the λ corresponding to any givenα. In a comparison of VAWT models, Ferreira et al. found induction factors to be as high as a = 0.25 on the front side and a = 0.60 on the back side of the turbine [21]. As a result, the angle of attack was shown not to exceed α = 10 • on the front side and α = −5 • on the back side. Table 2 provides a general idea of which λ correspond to theα used in this study by considering two exemplary values of a. In reality, a is a function of azimuthal position β. The expected oscillations in wind speed, which according to equation 1 range fromÛ r = 0.087 U r forα = 5 • toÛ r = 0.326 U r forα = 20 • , were neglected because the current experimental setup does not allow for periodic variations in tunnel wind speed.
The experimental data is compared to the analytical predictions of Theodorsen [9] and Motta et al. [13]. Despite the use of endplates, the low AR of the airfoil lead to a shallower lift curve slope ∂C l ∂α and a higher stall angle α ss than is expected for a 2D airfoil. Furthermore, the finite cross-section of the wind tunnel and associated blockage effects lead to a steeper ∂C l ∂α and a higher C l,max . Existing corrections were developed for steady flows, and their validity for unsteady flows is unclear. Furthermore, no established blockage corrections for wind tunnels with round crosssections are known to the authors. Therefore, no corrections were applied to the data presented here. The analytical predictions of Theodorsen and Motta et al. are based on thin airfoil theory and therefore assume a static lift curve slope of ∂C l ∂α = 2π. In order to compare these 2D theories to the experimental data, an aspect ratio correction described by [22] was applied in reverse. However, the use of the geometric aspect ratio AR g = 1.5 in the correction did not project the potential flow lift curve C l = 2πα onto the experimental static lift curve. This is likely due to the presence of the endplates and blockage effects. Therefore, in order to provide an accurate  Figure 1. Lift coefficient for a NACA 0021 airfoil sinusoidally pitching around a zero mean angle of attack at reduced frequency k = 0.15 and Re c = 10 6 . Note the differing axis scales.
comparison between the 2D theories and the experimental data, the effective aspect ratio used for the correction was empirically selected such that the correction projected C l = 2πα onto the experimentally obtained steady lift curve. Thus, the aspect ratio used for the correction was AR e = 2.8, and the correction should be regarded as purely empirical.

Results
Due to the high Reynolds number and the relatively low aspect ratio of the wing model, the static stall angle α ss = 21 • is higher than in most VAWT studies reported in the literature. In steady conditions, at Re c = 1 × 10 6 , the maximum lift coefficient, attained at α ss , was C L,max = 1.31. The airfoil began to stall gradually from the trailing edge starting at around 20 • . Thus, even  Figure 2. Drag coefficient for a NACA 0021 airfoil sinusoidally pitching around a zero mean angle of attack at reduced frequency k = 0.15 and Re c = 10 6 . Drag due to supporting rod was not subtracted, so values are purely qualitative at the highest amplitude of oscillation tested here,α = 20 • , the airfoil does not enter the stall regime at any point during the oscillation. Figure 1 shows the lift coefficient, C l , as a function of angle of attack for the various amplitudes tested. No stall effects are visible, as expected. For all tested amplitudes, the lift coefficient follows a clockwise hysteresis loop. Despite the high amplitudes of oscillation no amplitudedependent changes in the shape of the hysteresis loop are visible. Furthermore, the width of the hysteresis loop appears to agree with the theoretical predictions. Motta et al. predicted that for k larger than the phase inversion point at k i = 0.144, added mass effects will dominate over the circulatory contribution, causing a clockwise hysteresis loop [13]. At the inversion point, no hysteresis is predicted to occur. The present study was conducted just above the phase inversion point at k = 0.15. Thus, the clockwise hysteresis loops observed are consistent with the theory. Above the inversion point, an increase in airfoil thickness is predicted to reduce the width of the hysteresis loop. However, figure 1 indicates that this effect is small, if at all present. The experimental data show a steeper loop, leading to a maximum C l as high as the static C l at those angles. Both Theodorsen and Motta et al. predict a decrease in C l compared to the static condition. This decrease is due the lift deficiency factor C(k), which accounts for the shedding of vorticity into the wake. The experimental data suggests that this effect might be overestimated for the conditions studied here.
The drag coefficient is shown as a function of angle of attack in figure 2. The drag measurements were obtained using the force balance, and the drag due to the supporting rod was not subtracted. Since the supporting rod is cylindrical, we expect no additional drag hysteresis due to its presence. As such, the absolute drag values will not represent the drag on the airfoil, but have an unknown offset. However, we expect the characteristics of the drag hysteresis loop to be representative of the airfoil. As with the lift coefficient, the drag coefficient follows a hysteresis loop that widens asα increases, where the drag is higher as the angle moves away from zero and lower as it returns to zero. Both the qualitative shape and the direction of the drag curve are in agreement with those found by Leishman and Beddoes at angles below the stall angle [11]. In fact, they found that the pressure drag became negative on the downstroke. The models by Theodorsen and Motta et al. do not predict drag.
The moment coefficient around the quarter chord point, C m, c