Airfoil optimisation for vertical-axis wind turbines with variable pitch

To advance the design of a multimegawatt vertical-axis wind turbine (VAWT), application-specific airfoils need to be developed. In this research, airfoils are tailored for a VAWT with variable pitch. A genetic algorithm is used to optimise the airfoil shape considering a balance between the aerodynamic and structural performance of airfoils. At rotor scale, the aerodynamic objective aims to create the required optimal loading while minimising losses. The structural objective focusses on maximising the bending stiffness. Three airfoils from the Pareto front are selected and analysed using the actuator cylinder model and a prescribed-wake vortex code. The optimal pitch schedule is determined, and the loadings and power performance are studied for different tip-speed ratios and solidities. The comparison of the optimised airfoils with similar airfoils from the first generation shows a significant improvement in performance, and this proves the necessity to properly select the airfoil shape.


INTRODUCTION
While horizontal-axis wind turbines (HAWTs) are widely studied and have proven their capabilities, the alternative vertical-axis wind turbine (VAWT) concept is not yet common. To achieve competitiveness compared with the standard HAWT, improving the rotor performance is a key in order to make VAWTs commercial. Some researchers stated that VAWTs could achieve maximum power coefficients in line with current HAWTs 1,2 and even have the potential to decrease the cost of energy in urban applications 3 or for large-scale offshore wind farms 4,5 . To unlock this potential, research is performed on various aspects related to rotor aerodynamics such as the effect of rotor configuration (blade solidity, 6 double-rotor interaction, 7 three-dimensional losses 8,9 ), the effect of blade profile (airfoil selection 1,10,11 and roughness sensitivity 8,12,13 ), and the effect of using load control (trailing edge flaps, [14][15][16] circulation control, 17 and blade pitching [18][19][20] on the turbine efficiency.
This research focusses on the understanding of the desired airfoil characteristics and the design of optimised airfoil shapes for multimegawatt VAWTs with individual blade pitching.

Background
In the early years of the VAWT history, researchers mainly opted for symmetric airfoils since it was believed important that the behaviour was similar in the upwind and downwind part of the rotor and because abundant information was available for symmetric airfoils. Pioneering research in the airfoil design for VAWTs has been performed at Sandia National Laboratories at Tokai University and at Delft University of Technology. At Sandia National Laboratories, Klimas 1 investigated the desired characteristics of VAWT airfoils. Klimas stated that modest maximum lift coefficients and sharp stall are important for stall regulation together with a low and wide drag bucket for unsteady effects. Around the 1980s, the University of Tokai 10 sets a list of conditions for high aerodynamic efficiency, and they derived from simple streamtube models that a large lift gradient was important. Recently, airfoil optimisation for VAWTs has received new interest. At Delft University, Claessens 12 designed airfoils for a small-scale rotor, but contrary to earlier research, a smooth stall behaviour was preferred for noise reasons. Ferreira et al 11,13 generated new airfoils using a dual-objective genetic algorithm including an aerodynamic and structural optimisation function, which is validated experimentally in the DeepWind project. 21 Furthermore, researchers such as Kirke, 8 Claessens, 12 and Ferreira 11 all agreed that the sensitivity to roughness of the airfoils should be small.
Because VAWT rotor blades are subjected to a cyclic angle of attack variation, the turbine efficiency can be significantly improved by optimising the angle of attack at each point of the azimuthal cycle. Already in the 1980s, variable pitch has been investigated. Theoretical studies about the optimal and prescribed pitch-angle cycles as a function of tip-speed ratio and blade solidity have been studied by various researchers such as Natuvetty and Gunkel, 22 Zervos et al, 23 and Vandenberghe and Dick 24 or more recently by Kirke, 8 Kosaku,25 or Houf. 19 In addition to increasing the turbine efficiency, Kirke 26 stated that variable blade pitch can also assist self-starting and reduce the vibration of VAWTs.

Research objective
Previous airfoil optimisations for VAWTs exclude blade pitching; this research focusses on integrating individual blade pitch into the airfoil design.
The goal of this research is threefold: 1. Identify the desired rotor loading and determine the corresponding pitch cycle.
2. Optimise airfoils for VAWTs taking into account blade pitching and early transition because of surface roughness.
3. Analyse the performance of a VAWT with fixed and variable blade pitching.

Simulation of the 2D VAWT
The VAWT simulations are performed using the two different codes. The first model is a vortex model developed by Ferreira. 27 The second model is the actuator cylinder (AC) model developed by Madsen. 28 The vortex code U2DiVA is a two-dimensional vortex code following Prandtl's lifting-line theory in which the circulation distribution is determined. PW2DiVA is a variant of it using the prescribed wakes derived from the unsteady two-dimensional panel/vortex model U2DiVA 27 for different solidities and thrust coefficients. While the free-wake model generates higher fidelity results, it is computationally more expensive.
The prescribed-wake approach, on the other hand, is relatively quick in computation and only takes a few seconds to converge (compared with a couple of hours for the free-wake model). The forces of the blades are determined using the lift and drag polars of an airfoil. With the airfoil forces known, the bound vortex is determined with the Kutta-Joukowski theorem, and the induced velocity at a certain location is given by the Biot-Savart law. The vortex is consequently shed into the prescribed wake.
The AC is a 2D flow model extending the actuator disk concept. An actuation surface is introduced that coincides with the swept area of the rotor, which is for a VAWT, a cylindrical surface. The reaction of the blade forces is applied on the flow as body forces. The solution builds on the 2D Euler equations and the equation of continuity. The induced velocities are prescribed by the volume forces and induced forces and consist of a linear and non-linear solution. The so-called Mod-Lin solution uses only the linear version of the AC model and a correction to account for the non-linear part, instead of solving the computationally expensive non-linear solution. This method requires less resources and still delivers accurate results. 29 The computation time is in the same order as for the prescribed-wake vortex code.
All models predict similar VAWT performance. For a VAWT equipped with two blades, with a solidity of 0.06, running at a tip-speed ratio of 3, the results are given in Figure 1. For this case, the thrust coefficient is 0.51. In Figure 2, a highly loaded case with a thrust coefficient of 0.85 is presented for a two-bladed VAWT with a solidity of 0.1 and tip-speed ratio of 4. The polars used in these calculations are the inviscid polars C l = 1.1 · 2 and C d = 0. The representation of the angles and the sign convention of the load coefficients are presented in Figure 3. These examples show that the codes are in good agreement in the upwind part of the rotor. The differences between the three models increase with rotor loading, by increasing either tip-speed ratio or rotor solidity. Despite this, the difference in power coefficient predicted by the AC, U2DiVA, and PW2DiVA remains small. These findings are in agreement with previous studies of Ferreira et al. 30 Note that U2DiVA does take into account the dynamic effects while the other codes do not. The discontinuities in the vortex-model results (mainly downstream) are due to the interaction of the blades with the wake.

Rotor loading and pitch cycle
According to Madsen, 31 the optimal rotor loading is determined by optimising the pressure-jump distribution along the actuator surface, in which the actuator surface coincides with the swept area of the turbine. The pressure jump is created by the load distribution over the actuator surface.   For an HAWT, Glauert 32 has proven from the derivation of the Betz limit that uniform loading on the actuator disk is optimal. For an HAWT, this means that the circulation on the disk is uniform and that the shed wake is composed of trailing vortices mainly at the edge of the actuator disk. The strength of the wake is defined by the total bound circulation on the blade. The optimal shed vorticity is given by a constant value of C l cV rel , according to the Kutta-Joukowski theorem. c represents the chord length, V rel is the relative velocity, and C l is the lift coefficient.
Similarly, Wilson 33 and Madsen 29 have shown that for a VAWT with infinitely high tip-speed ratio, the maximum power coefficient is found for a rotor loading approaching uniform distribution of normal loading in the upwind half of the rotor and a uniform distribution of normal loading  in the downwind half of the rotor. The difference in magnitude between the upwind and downwind distributions defines the non-conservative force field that is responsible for the energy exchange and the generation of the wake. In 2D, this load case results in the same velocity (and wake) field as the 2D actuator line, which represents the 2D actuator disk, and should therefore have the same optimum. Because the AC is not limited to two points of releasing vorticity as the 2D actuator disk, the solution of the 2D AC might result in a power coefficient above the Betz limit, as proposed by Madsen. 29 Assuming infinite tip-speed ratio, the loading distribution discussed in the previous paragraph (uniform distribution of normal loading upwind and a uniform distribution of normal loading downstream, which differ by a constant) can be translated to a distribution of bound circulation that tends to be 0 both upwind and downwind but with a jump at the poles of the cylinder, corresponding to the jump of loading at the poles. Kelvin's and Helmoltz' theorems imply that this change of the bound circulation (although tending to be 0) will generate a finite amount of vorticity at the poles, creating the wake.
We can now approximate the case of finite number of point vortices at a finite tip-speed ratio. We can also imply that these point vortices The optimal performance of a VAWT is thus prescribed by the variation of the circulation at the edges of the actuator. This means that for large solidities, a smaller d(C l V rel )| upwind-downwind is required than for smaller solidities. To demonstrate this, Figure 4A is provided in which three optimal circulation distributions normalised with the chord are presented.
To achieve the optimal distribution of circulation, an infinite amount of combinations is available. By redistributing the loads over the upwind and downwind part of the rotor, different circulation distributions can be achieved leading to similar power outputs. To illustrate this, Figure 4B is provided in which three different circulation distributions are presented with the same power output. The possibility to obtain multiple distributions with the same power conversion was also demonstrated by Madsen 31 and Ferreira et al. 34 They stated that when having the same average loading, the power conversion will not change beyond what can be attributed to the accuracy of the model.
Ideally, all extracted energy from the flow should be converted into mechanical power, ie, the torque T causing the turbine to rotate. However, this is limited by the effect of drag, as can be understood by Equation 1. Assuming that the chord and velocity field are fixed (relative velocity, , which is a fair assumption for a fixed tip-speed-ratio case, the optimal circulation distribution is the one corresponding to the one giving the optimal jump in circulation and at the same time minimises drag throughout the rotation.
Variable pitch can be used to alter the angle of attack to achieve the desired circulation. Because the induced velocity field is not truly constant, the pitch cycle is adapted iteratively to approach as closely as possible the target circulation distribution using an optimisation scheme.
Constraints could be added to the pitch rate and maximum pitch deflection to avoid solutions that could be infeasible in a real application.
The required pitch cycle and the corresponding angle of attack, velocity, and loading distribution for a VAWT with solidity of 0.1 and tip-speed ratio of 4 are presented in Figure 5. This schedule is calculated using the AC model, but similar results were found with the prescribed-wake vortex model (PW2DiVA). The results are in agreement with previous findings from, eg, the second-order harmonic pitch from Vandenberghe and Dick 24 and the variable pitch law from Zang et al. 20 The case with pitch yields around 4.5% more power.
So far, inviscid polars were used. However, one should realise that using viscid airfoil polars instead might not always allow to operate at constant circulation as prescribed by the optimal loading. Especially at low tip-speed ratios or solidities, a too large lift coefficient is required to fulfil the optimal circulation condition. As an example of this, the optimal pitch for a rotor with solidity of 0.08 and tip-speed ratio of 3 is presented for the inviscid and viscid polars of the NACA0018 in Figure 6. Between the azimuthal angles of 90 • and 280 • , the angle of attack is constant yielding maximum lift coefficient, and the circulation is below the optimal. One could expect that the lower the solidity and tip-speed ratio become, the more the optimal angle of attack distribution will converge to a uniform angle of attack distribution.

Airfoil optimiser
To perform the airfoil optimisation, a multiobjective optimisation scheme is built. It uses the build-in MATLAB genetic algorithm function gamultiobj to find a balance between the aerodynamic and structural needs of an airfoil, since the structural and aerodynamic optima are not Further, the airfoils are discretised using the class-shape transformation (CST) surface parameterisation method, 36 and constraints are defined to result with realistic airfoil shapes. The optimisation scheme starts with an initial generation, which is further developed through an evolutionary process from one generation to another for a predefined number of generations. The solution of the optimisation algorithm creates a Pareto front.

Aerodynamic objective function
In the work of Ferreira et al, 11 a relation is derived to transform the aerodynamic objective at rotor scale to an objective on airfoil level of which the main steps will be repeated.
On rotor scale, the objective is to approach the optimal loading. As explained before, the optimal rotor loading is given by a unique distribution where c · |C l | = const.
The final expression of the aerodynamic objective function is given in Equation 3 and corresponds to what is found by Ferreira. 11 The aerodynamic objective at rotor scale was to extract as much energy and approach optimal loading. On airfoil scale, this is translated to maximize the lift gradient over drag in the angle of attack range encountered by the rotor. The term b( ) represents a weighting function defining the occurrence of a certain angle of attack over one rotation. The domain in which the integral is evaluated is defined by a range of angle of attack Δ over which the turbine operates and a minimum angle of . can be varied from the minimum angle of attack for which the lift polar is given till the maximum angle of attack minus the angle of attack range and can be understood as a fixed pitch. Δ is a variable, since pitch allows to change the angle of attack to the desired one.
To compute the aerodynamic properties of the airfoils, the optimiser is connected to XFOIL. The lift-curve slope and drag coefficient are calculated as a function of the angle of attack at a fixed Reynolds number. Because wind-turbine blades often suffer from roughness due to dust, insect accumulation, or surface damage, it is important to account for possible early turbulent transition in the airfoil optimisation. In the optimisation, this is taken into account by considering the airfoil sensitivity to promote turbulent transition. This is realised by adapting the limit amplification factor N crit of the e N transition model to 9.00 (normal sensitivity) and 0.01 (high sensitivity). 37 This second case will cause transition to occur closer to the leading edge.

Structural objective function
The structural objective function calculates the bending stiffness in flapwise direction of the airfoil per wall thickness, defined in reference to the centroid of the airfoil, 11 as shown in Equation 5. This parameter is defined as dominant.
Wind-turbine blades are normally stiffened inside by, eg, I-beams or a box structure. However, a full structural optimisation in terms of internal structure and material is outside the scope of this work.

Airfoil discretisation
To perform an optimisation, it is desirable to limit the number of design variables. By parameterising the airfoil geometry, the model complexity can be reduced. A typical method to represent an airfoil is the CST surface parametrisation method. 36 In the CST method, an airfoil shape is represented by an analytical expression consisting of three parts, ie, the class function C, the shape function S, and the trailing-edge function T . The expression is given in Equation 6. The class function defines the general geometric shape, while the shape function defines the specific shape of the object using a sum of Bernstein polynomials. The last function controls the trailing-edge thickness. 36 The CST method builds on the coefficients defining these functions. The amount of coefficients prescribing the shape function depends on the order N of the Bernstein polynomial used for the upper or lower surface. The trailing-edge function introduces one more coefficient.
By representing a variety of existing airfoils by the CST parameterisation, it is found that the average error between the original airfoil coordinates and the parameterised coordinates approaches an asymptotic value from an order of the Bernstein polynomial of 6, as proven by Figure 7. Therefore, an order of 6 per surface should suffice to capture different airfoil shapes in this optimisation.

Constraints
To obtain realistic airfoils, constraints should be implemented into the optimisation algorithm. Infeasible airfoils should be excluded. On the basis of the author's experience, the following constraints are set to the optimiser: • The Bernstein polynomial coefficients should remain between −1 and 1. This constraint is set to limit the maximum local thickness of the airfoil.
• The upper and lower surface of the airfoil should not intersect.
• The airfoil thickness should gradually decrease in the last 20% of the airfoil chord. This constraint prevents infeasible trailing-edge shapes.
• The angle formed at the trailing edge should be larger than 8 • . This value is based on previous experience.
• The trailing-edge thickness should be within 0.001 and 0.05 of the chord length. The minimum value is set because of manufacturing reasons while the maximum is set due to limitations of XFOIL to model thick trailing edges.
• XFOIL should be able to calculate (converge) the aerodynamic polars of the airfoils in clean and rough state for at least an angle of attack range from −15 • to 15 • .

Validity objective functions
The objective functions are analysed by studying over 2000 existing airfoils. The airfoil database consists of a broad range of airfoils designed for wind energy and aircraft applications. The airfoil polars are calculated using the XFOIL at a fixed Reynolds number of 5E6, and the performance of the VAWTs is determined using both the AC and the vortex code with prescribed wake; however, only the results calculated with the AC will be presented. Figure 8A shows the relation between the moment of inertia I xx and the maximum obtained power coefficient C p | max (for a tip-speed ratio of 4 and a solidity between 0.05 and 0.13). In general, the larger the moment of inertia of the airfoil, the lower the power coefficient is, especially for of the data can be related back to the discrete intervals used for the solidity and the limited code accuracy. Also, the aerodynamic objective function assumes an invariant velocity field and constant angle of attack distribution while this is not fully true, though an acceptable assumption.

Airfoil-operating conditions considering the pitch cycle
One of the input parameters in the aerodynamic objective function is the distribution of the angles of attack (b[ ]). By allowing variable pitch, the angle of attack can be altered to any desired angle of attack yielding maximum efficiency. A rotor operating at an angle of attack distribution corresponding to uniform circulation will yield optimal power. To include this in the airfoil optimisation, the angle of attack distribution function b( ) should be modified compared with the original one derived by Ferreira et al. 11 Figure 9A provides the angle of attack over one rotation with a fixed and a variable pitch schedule optimised for maximum power. The second case shows a more evenly distributed angle of attack between the upwind and downwind part of the rotor. Both distributions are plotted

Final input parameters
The final optimisation routine is run for a specific set of input parameters. The final input parameters used to run the optimiser to generate airfoils for the VAWT including pitch are summarised in Table 1. The weight function derived earlier is used. Further, the polars are calculated at a Reynolds number of 6.4 · 10 6 , which is around the average Reynolds number encountered by a multimegawatt turbine and is still accurately modelled by XFOIL. The clean and rough airfoil state are weighted with a weight factor w of 0.5.
The population size is set in relation to the number of control points, 38 namely, to 10 times the number of control points. 39 Because the airfoils are prescribed by the CST method with 13 variables, the population size is set to 130. By studying the evolution of the Pareto front, it is found that the genetic algorithm seems to have converged, in the region of interest, after 30 generations. To make sure that the genetic algorithm does not tend to converge to one particular airfoil family, different initial populations are used to run the optimisation code. The different initial generations include several NACA four-series and six-series airfoils, wind-energy-application airfoils from the Delft University (DU), Wortmann (FX), and Althaus (AH) airfoils. Despite the large initial source population, it cannot be guaranteed that the initial variety of the population covers all relevant possibilities.

Optimised airfoils
From the optimisation routine, thousands of airfoils are generated and analysed. The results are combined, and the Pareto front shown in Figure 10A is obtained. Each data point represents a different airfoil.
The Pareto front presents a strictly decreasing trend between the aerodynamic objective function and the structural objective function. The thinnest airfoil on the Pareto front has a thickness around 15%c, saying that thin airfoils not necessarily perform better aerodynamically. Thicker airfoils lead to better performances up to the point where the increased lift slope is cancelled out by the losses because of viscous effects. From the final Pareto front, three airfoils are selected with a thickness of 20%c, 25%c, and 30%c (see Figure 10B). They are named as DU17VAWT200,  DU17VAWT250, and DU17VAWT300, respectively. The coordinates of the airfoils are provided in the Appendix (Table A1). Because the weight function in the aerodynamic objective function is close to symmetric, the optimised airfoils are nearly symmetric too.

Polars
The lift and drag polars calculated with XFOIL of the three selected airfoils are presented in Figure 11A-D in clean and rough state, respectively.
The lift slope increases with the relative airfoil thickness, as has been confirmed earlier for Joukowski airfoils by Katz and Plotkin 40 and Hoerner and Borst 41 (for airfoils up till a certain thickness). Also, the minimum drag increases with airfoil thickness. Higher transition sensitivity, realised by decreasing the N-factor, causes the lift gradient to decrease and drag coefficient to increase.

Performance
The performance of the optimised airfoils is evaluated using the AC model and the prescribed-wake vortex model (PW2DiVA) with fixed and variable pitch. The power coefficient is calculated for a two-bladed VAWT application with different blade solidities and tip-speed ratios, and the results are presented in Figure 12. These results are calculated using the AC model and use the polars with high sensitivity to turbulence (N crit = 0.01). Uncoloured areas indicate the convergence limit of the model.
For the 25% thick airfoil (DU17VAWT250), Figure 13A is provided. In this figure, the difference in power coefficient for fixed and variable pitch normalised with respect to the fixed pitch case is presented. Depending on the tip-speed ratio and solidity gains above 30% is achieved.
For highly loaded cases, the reachable gains are limited while they increase with decreasing tip-speed ratio and solidity. In Figure 13B, the angle of attack range over which the turbine operates with variable pitch is presented. This range decreases with increasing solidity and tip-speed ratio. This follows the expectations as explained in Section 2.2. The optimal rotor loading is given by a constant d(C l cV rel )∕d , and thus, decreasing velocity or chord length requires a larger lift coefficient and thus a larger angle of attack. Note that the maximum angle of attack range corresponds to the angle of attack range between minimum and maximum lift coefficient. For the other airfoils, similar observations can be made.
Although maximum gains up to 30% can be achieved, the overall power coefficient gain that can be reached is limited from 2.17% to 7.96% depending on the airfoil thickness. The maximum achievable power performance of the three optimised airfoils using the polars with N crit = 9.00 with both fixed and variable pitch is very similar for the three airfoils. However, using the polars with N crit = 0.01, a clear degradation in performance can be identified for the thicker airfoils. The gains from using the variable pitch instead of fixed pitch for the N crit = 9.00 case are  significantly smaller then in case of N crit = 0.01, and therefore, variable pitch narrows the loss because of turbulence sensitivity. This information is presented in Figure 14 and Table 2.

Optimised airfoils vs four-digit NACA
To discuss the performance of the optimised airfoils, a comparison will be made to airfoils with a similar thickness that were part of the initial generation of the optimisation, ie, four-digit symmetric NACA airfoils. The comparison consists of three parts: (1) maximum power coefficient, (2) optimal blade solidity, and (3) pitch schedule. Figure 15 and Table 3 present a comparison of the maximum achievable power coefficient for the three optimised airfoils (DU17VAWT200, DU17VAWT250, and DU17VAWT300) against NACA airfoils with similar thickness (NACA0020, NACA0025, and NACA0030). The polars are calculated at Re = 6.4E5 with N crit = 9.00 and N crit = 0.01. The performance is calculated for a two-bladed VAWT with the AC model. The power coefficient represents the maximum achievable power coefficient when considering variable pitch for a tip-speed ratio interval from 2.5 to 6.0 and a solidity between 0.05 and 0.13. The optimised airfoils clearly outperform the reference NACA airfoils, especially in N crit = 0.01 state. This was expected since in the airfoil optimisation procedure, the emphasis was set on polars with higher sensitivity to turbulent transition.
Because the chord is adjusted to reach optimal d(C l cV rel )∕d , the gain in power can be mainly attributed to the decrease in drag.
For every tip-speed ratio, one particular blade solidity yields maximum power output. In Figure 16A, the optimal solidity as a function of tip-speed ratio is presented for the optimised DU17VAWT250 airfoil and the NACA0025 airfoil, both for N crit = 9.00 and N crit = 0.01, and with variable pitch. The figure indicates that a smaller blade chord is required for the optimised airfoils. Minimising the chord could be beneficial in terms of material cost and loading in parked position. For the other airfoils, similar observations are made. In Figure 16B, the pitch schedules for the same airfoils with N crit = 9.00 and N crit = 0.01 are presented. These pitch schedules correspond to the case in which the tip-speed ratio is 3.5 and the blade solidity is as provided in Figure 16A. Although for the different cases, the fixed pitch and optimal angle of attack distribution are different, the pitch schedule is rather similar. However, the optimised airfoils require smaller changes in pitch. FIGURE 15 Maximum achievable power coefficient for the three optimised airfoils (DU17VAWT200, DU17VAWT250, and DU17VAWT300) and four-digit NACA airfoils (NACA0020, NACA0025, and NACA0030) for N crit = 9 and N crit = 0.01 with variable pitch [Colour figure can be viewed at wileyonlinelibrary.com]

Limitations
The validity of the optimisation results and performance results is limited by the assumptions made and the codes used.
Firstly, in the derivation of the aerodynamic objective function, an assumption is made that the velocity field is constant. This is not completely true. When changing the rotor loading, the induced velocity field will be affected. However, for a fixed tip-speed ratio, the induced velocity is only a minor portion of the relative velocity at the rotor blades, and therefore, this assumption is expected not to affect the results significantly.
Secondly, the airfoil optimiser builds on the XFOIL calculation of the lift and drag polars in clean and rough state. In general, the lift gradient is overestimated while the drag is underestimated. The drag effect of thick trailing edges is underestimated as well. 42,43 XFOIL still misses some accuracy in the stall prediction. 44 Although small variations on these parameters exist, this has only a small impact on the VAWT simulations.
Further, note that there is not accounted for phenomena such as flow curvature (cambering effect of the airfoil) or dynamic stall. Also, the fact that VAWTs operate at a range of Reynolds numbers throughout the operation is not taken into account.
Thirdly, the power calculations are limited by the capabilities of the models. Both codes make assumptions while computing the induced velocities and lack in accuracy at high rotor loadings. To make sure singularities and computational instabilities are filtered out, the AC model and the prescribed-wake vortex code are always used in parallel. The power coefficients are also relatively high due to the neglecting of struts and tower or any unsteady effects such as dynamic stall or blade-vortex interaction. Also, the modelling is two-dimensional, causing the 3D effect such as tip losses to be neglected. Furthermore, the pitching moment of the blades is neglected, as well as the additional velocity component due to pitching.
Ideally, the optimal pitch schedule should be determined freely, without any constraint on the required loading. However, this causes the code converges to its numerical limits. Therefore, the variable pitch schedule is determined iteratively using an optimisation scheme assuming that the circulation upwind and downwind should be constant.
The optimisation routine performed in this work already considers the aerodynamic and structural needs of airfoils to operate in a VAWT including blade pitching; however, additional requirements could be added. These could include internal blade structure, aeroelastic behaviour, production costs, circulation control, unsteady effects, or turbulent inflow. However, this is out of the scope of this research and irrelevant when analysing the effect of variable blade pitching to airfoil optimisation.

CONCLUSION
In this work, an airfoil optimisation routine is set up to design blade shapes for a multimegawatt VAWT with variable pitch. The aerodynamic objective at rotor scale is to extract as much energy and approach optimal loading. On airfoil scale, this is translated to maximise the lift gradient over drag in the angle of attack range encountered by the rotor. The angle of attack distribution for a variable pitch turbine is more evenly distributed compared with the fixed pitch turbine. The optimiser is connected to XFOIL to obtain the polars and takes into account the sensitivity to turbulence. The structural objective function focusses on maximising the bending stiffness in flapwise direction.
Using this airfoil optimisation routine, airfoils are generated that form a Pareto front. The minimum airfoil thickness of the Pareto front is around 15%c, saying that thin airfoils not necessarily perform better aerodynamically. Also, the optimised airfoils are nearly symmetric, which is a direct result of the angle of attack distribution to be close to symmetric as well. Three airfoils from the Pareto front with a maximum thickness from 20% till 30% are selected. They are referred to as DU17VAWT200, DU17VAWT250, and DU17VAWT300. The airfoils are analysed using the two VAWT simulation codes, namely, the AC model and a prescribed-wake vortex code. The maximum achievable power coefficient is rather constant for the three airfoils in clean state while thicker airfoils suffer more from roughness causing a degradation in performance. Introducing variable pitch allows to increase the achievable power coefficient and narrows the loss because of early turbulent transition.
Comparing the performance of the optimised airfoils with reference airfoils from the initial generation with similar thickness showed the improved efficiency of the newly designed airfoils. The optimised airfoils outperform the reference airfoils in terms of power coefficient, while they allow to decrease the optimal blade solidity for a given tip-speed ratio and pitch-angle range.