Design of the LRP airfoil series using 2D CFD

This paper describes the design and wind tunnel testing of a high-Reynolds number, high lift airfoil series designed for wind turbines. The airfoils were designed using direct gradient- based numerical multi-point optimization based on a Bezier parameterization of the shape, coupled to the 2D Navier-Stokes flow solver EllipSys2D. The resulting airfoils, the LRP2-30 and LRP2-36, achieve both higher operational lift coefficients and higher lift to drag ratios compared to the equivalent FFA-W3 airfoils.


Introduction
One of the main challenges to the continuous up-scaling of wind turbines is to maintain low blade mass while achieving high power capture and the necessary stiffness of the blade. An important parameter to achieve high stiffness is the cross-sectional thickness along the blade. Modern slender blade designs therefore move towards using airfoils of higher relative thickness also on the outer part of the blade. This trend poses considerable challenges to the aerodynamic design, since thick airfoils in general are less efficient than thinner ones.
Airfoil design using direct methods for shape optimization coupled to solvers such as XFOIL is a fairly well established method, and has at DTU Wind Energy been used to design a number of airfoil series now in use on modern wind turbines [1,2]. It is, however, well-known that the accuracy of panel codes decreases for airfoils with thick trailing edges and relative thicknesses greater than 30% to 36%, making use of such codes more challenging for design of thick airfoils. Although it can also be challenging to predict the correct performance of thick airfoils with Computational Fluid Dynamics (CFD), CFD is more accurate than panel codes, particularly close to stall and for thick airfoils.
In this paper a new optimization framework for airfoils is presented, which is based on the open-source framework OpenMDAO that has been interfaced to a range of simulation codes. For this work, the framework has been coupled to both the panel code XFOIL and the Navier-Stokes solver EllipSys2D and a robust and efficient optimization strategy using a combination of XFOIL and 2D CFD was developed. Since objectives and constraints typically consist of a range of evaluations both at different angles of attack, but also fully turbulent and transitional simulations, a nested parallel approach is used in order to speed up the optimization runs.
The tool was used to design an airfoil series, the Light Rotor Project (LRP) airfoil series, aimed for operation at high Reynolds numbers and designed based on a series of aerodynamic and structural constraints derived through a close collaboration between DTU Wind Energy and Vestas. The key objectives were to obtain airfoils that could operate with high design lift coefficients, achieving high performance for rough surface conditions while not compromising performance under clean conditions, while ensuring high structural stiffness through a series of geometric constraints. The resulting airfoils, the LRP2-30 and LRP2-36 airfoils, are presented in this paper. Numerical predictions are presented over a range of Reynolds numbers along with wind tunnel measurements carried out in the Stuttgart Laminar Wind Tunnel at Re = 3 × 10 6 . Finally, the new airfoils are compared to measurements carried out in the same wind tunnel on the FFA-W3-301 and FFA-W3-360 airfoils.

Numerical Optimization Setup Using 2D CFD
The optimization framework developed for this work is written in Python and uses OpenMDAO [3] as the basic framework for executing and connecting the different components in the optimization. Figure 1 shows a flowchart where the major components are included: The Driver, which in this case is the optimization algorithm, the airfoil geometry builder, which generates a parameterized airfoil shape, the mesh generator, which takes an airfoil shape as input and generates a CFD ready mesh, the flow solver which computes the flow over the airfoil and resulting forces acting on it, and finally the objective function which computes a measure of goodness based on the flow solution. In the following a brief overview of each of the components is provided.

Driver
Through interfaces provided by OpenMDAO a wide range of optimizers are available in the optimization tool. However, the driver used in this work, was developed in-house and is based on a gradient-based Sequential Linear Programming (SLP) optimizer, developed at DTU Wind Energy, also used in the codes AirfoilOpt and HAWTOPT [4].

Airfoil Shape Parameterization
To effectively design an airfoil a suitable parameterization of the airfoil shape is necessary, which with as few parameters as possible can produce smooth generalized shapes. In the present work, Bezier curves are used since this class of curves are characterized by a high degree of controllability, while still maintaining gentle curvature changes. Figure 2 (a) shows a fit of two 5th order Bezier curves to the NACA-0030 airfoil, where the curves meet at the leading edge. The fit is achieved using a least squares minimization of the difference between the original and fitted coordinates. To ensure C1-continuity across the curve boundary at the leading edge, the starting points of the two Bezier curves must coincide, and the first two control points of the two curves must be co-linear, which has been constrained to a line in the x=0 plane. Subsequent manipulation of the airfoil shape can now be achieved by moving the position of the control points. In Figure 2 (b) the last interior control point on the pressure side is moved to y=0, with a subsequent fairly large change in shape noticeable even at x/c=0.15. An alternative to Bezier curves are B-spline curves, which offer more local control of the shape, which for general shape fitting can be viewed as an advantage. It does, however, also result in a more "lively" curve, which could require more geometric constraints to avoid e.g. sudden curvature changes.

Mesh Generation
A central part of a robust CFD-based optimization strategy is the meshing, which needs to be very robust in order to produce valid meshes for all design points. The meshes for the optimization were generated using the hyperbolic mesh generator HypGrid [5] using 256 cells in the chordwise direction and 128 cells in the normal direction and a first cell height of 5 × 10 − 6, which results in a y+ of less than 2 for a Reynolds number of 9×10 6 . For the subsequent evaluation of the airfoils, meshes with 768×256 cells were used. Figure 3 shows a typical mesh generated during the optimization.  Figure 3. 2D CFD mesh produced using the automated meshing scripts.

Airfoil Flow Simulation
The CFD simulations were carried out using EllipSys2D [6,7,8] using an automated setup, which meshes the airfoil and executes a series of simulations over a range of angles of attack.
The simulations were all steady state simulations since for low to moderate angles of attack, steady state is sufficiently accurate. The airfoil could be evaluated assuming either a fully developed turbulent boundary layer on the airfoil surface with the k − ω SST model by [9], as well as simulations where the boundary layers were allowed to transition freely from laminar to turbulent flow using either the γ − Re θ correlation based transition model of Menter et al. [10,11] or the Drela Giles transition model. For the optimizations the correlation based model was used due to its robustness, while the subsequent evaluations were carried out using the Drela Giles model. To evaluate whether the simulations were sufficiently mesh independent for an optimization context, the flow over the airfoil was evaluated using a finer mesh with 768×256 cells around the airfoil. Figure 4 shows airfoil polars computed using the 768×256 and 256×128 meshes for fully turbulent and transitional flow, respectively. Evidently, both the fully turbulent and free transition simulations are in good agreement for the two mesh resolutions, with a sufficiently small but consistent difference between solutions on the two meshes. Although the most accurate solution possible would be desirable for the airfoil optimization, the coarser mesh was chosen, since the overall performance of the airfoil is captured very accurately using this mesh.

Problem Formulation
The objective function was defined as where n is the number of flow cases for which the airfoil was evaluated using the CFD solver, w = [0, 1] is a weight assigned to each case, clcd norm is a normalization factor, and Cl(i) and Cd(i) are the lift and drag coefficients for each of the cases. The weights w were adjusted in the individual optimizations to obtain a desired property, and clcd norm was set according to realistic targets for each angle of attack evaluated in the objective function.
The initial angles of attack that were evaluated were α init = [2, 7, 12], but to allow the optimizer to locate the best operational range of the airfoil the variable ∆α was also added as a design variable, which controlled the offset of the angles of attack at which the flow solution was computed: Besides the basic bounds on the design variables, different types of constraints were imposed in the airfoil optimizations on both geometric and flow-based parameters. The primary geometric constraints were placed on parameters thickness as function of chord, camber line as function of chord, camber line angle as function of chord, curvature of the airfoil shape, while flow based constraints were enforced on e.g. transition point locations. To minimize the difference between the clean and rough C l−max , the transition point close to C l−max was forced towards the leading edge. The geometric constraints were chosen such that the structural properties of the airfoil would resemble those of the FFA-W3 airfoil, ensuring ample room for e.g. the spar caps in the internal structure. For a detailed study of the trade-offs between aerodynamic and structural objectives in the design of airfoils, refer to the work by Bak et al. [2]. The flow based constraints were chosen such that the designs would tend towards an efficient airfoil under clean conditions, while not compromising the performance for rough conditions, thus minimising roughness sensitivity.
The design variables in the problem were the x and y coordinates of the control points of the two Bezier curves describing the airfoil shape to allow for maximum flexibility for the optimizer. The y-coordinates of the trailing edge points were also free, allowing the optimizer to generate flatback or Gurney flap-like shapes.
The procedure followed in this work had two steps, which involved a preliminary optimization using XFOIL with an initial shape fitted to a NACA00XX airfoil, followed by an optimization on the resulting airfoil using EllipSys2D, with the same design variables, constraints and objective function. With six flow cases per objective evaluation, a total of 21 design variables and a total of 20 optimization iterations, a total of 2640 CFD simulations were carried out for each optimization (including gradient evaluations). The six cases used to evaluate the objective were computed concurrently on a total of 48 cores, and with an average simulation time of 60-70 seconds, an optimization completed in approximately eight hours.    Tables 1 and 2 show a summary of the performance of the two LRP airfoils compared to the FFA-W3 and Risø C2 airfoils. The operational point was set to be five degrees below the point of maximum lift, α oper = α C l−max -5 deg.. Compared to the FFA-W3 airfoils the LRP airfoils are designed to operate at a higher lift level, while still maintaining good L/D performance. As such the LRP airfoils achieve higher L/D ratios that the FFA-W3 airfoils both at turbulent flow and free transition conditions. The LRP airfoils have very similar performance to the Risø C2 airfoils, although the C2 airfoils achieve a higher L/D ratio. This is due to the fact that the structural constraints imposed on the LRP airfoils result in a thicker pressure side towards the trailing edge, which has a direct impact on the drag.     Table 1. Performance summary for the LRP2-30 airfoil compared to the Risø C2-30 and FFA-W3-301 airfoils computed at Re = 9 × 10 6 using EllipSys2D.

Wind Tunnel Measurements
have a turbulence intensity of 0.1%. Note that drag wake rake measurements are only available in the attached flow regions, since measurements outside this region are not sufficiently accurate due to flow unsteadiness. While these are lower Reynolds numbers than the airfoils were designed for, such measurements are important for validating the predicted performance of the airfoils.  Table 2. Performance summary for the LRP2-36 airfoil compared to the Risø C2-36 and FFA-W3-360 airfoils computed at Re = 9 × 10 6 using EllipSys2D. Figures 9 and 10 show a comparison of the measured and computed lift and drag polars of the LRP airfoils. The agreement for clean conditions is quite good with only a small overestimation of C l−max , although for C l > 0.8, the measured drag coefficient is consistently higher than predicted. This was identified as being partly due to a small trailing edge separation, which was not seen in the transitional simulations but only in the fully turbulent ones at the lowest Reynolds number of 3.0 × 10 6 . It is at present unclear what is the explanation for the overall higher drag level in the experimental data compared to computations. The measurements with zig-zag tape on the leading edge at x/c=[0.05, 0.1] on the suction and pressure side, respectively, show a poorer performance than the fully turbulent CFD simulations. A likely explanation for this could be that the tape due to its height not only trips the boundary layer, but also introduces an additional disturbance of the boundary layer that is not modelled in the simulations. Figures 11 and 12 show a comparison between the LRP2-30 and LRP2-36 airfoils with the equivalent FFA-W3-301 and FFA-W3-360 airfoils, also tested in the Stuttgart Laminar Wind Tunnel. As also shown in Tables 1 and 12 the LRP airfoils achieve an overall higher lift level for both rough and clean conditions, while maintaining comparable or improved L/D ratios.

Conclusion
In this paper a new optimization tool based on 2D CFD was presented. The tool was used to optimize two airfoils, a 30% and 36% airfoil, the LRP2-30 and LRP2-36, respectively. The predicted performance of the new airfoils showed that they achieved high operational lift coefficients of 1.57 and 1.32, respectively, while maintaining high L/D ratios. Their performance  was improved compared to the FFA-W3 airfoils in terms of the maximum lift and lift to drag ratio. Results from wind tunnel tests in the Stuttgart Laminar Wind Tunnel on the LRP2-30 and the LRP2-36 airfoils at Reynolds numbers below the design Re, also showed an improved performance with a higher lift level while maintaining a good L/D ratio.

Acknowledgments
The present work was funded by the project Light Rotor, EUDP 64010-0107. Wind tunnel measurements were co-funded by the Light Rotor project and Vestas. Computation were made possible by the use of GORM PC-cluster at the DTU central computing facility at DTU Risø Campus.