Surrogate based aeroelastic design optimization of tip extensions on a modern 10MW wind turbine

Advanced aeroelastically optimized tip extensions are among rotor innovation concepts which could contribute to higher performance and lower cost of wind turbines. A novel design optimization framework for wind turbine blade tip extensions, based on surrogate aeroelastic modeling is presented. An academic wind turbine is modelled in an aeroelastic code equipped with a near wake aerodynamic module, and tip extensions with complex shapes are parametrized using 11 design variables. The design space is explored via full aeroelastic simulations in extreme turbulence and a surrogate model is fitted 5 to the data. Direct optimization is performed based on the surrogate model, seeking to maximize the power of the retrofitted turbine within the ultimate load constraints. The presented optimized design achieves a load neutral gain of up to 6% in annual energy production. Its performance is further evaluated in detail by means of the near wake model used for the generation of the surrogate model, and compared with a higher fidelity aerodynamic module comprising a hybrid filament-particle-mesh vortex method with a lifting-line implementation. A good agreement between the solvers is obtained at low turbulence levels, 10 while differences in predicted power and flap-wise blade root bending moment grow with increasing turbulence intensity.

The power performance and ultimate loads of every design are evaluated in a single load case, comprising an IEC specific 55 DLC1.3 (IEC, 2005) simulation at 8m/s. This case is considered representative for determining the average power performance in below rated power operation and the range of peak loading, since the turbine operation ranges from low power production to full rated power within the simulation time.

Tip extension parametrization
The definition of the tip extension design variables and their design space is probably the most important step in the described 60 optimization process. The variables have been chosen in a way which enables a general blade stretching design capability.
Their range is a result of many prior parametric studies, and it is limited to ensure the validity of the aerodynamic modeling.
The 11 chosen variables and their extend in the design space are shown in Table 1 Fig. 3. The main controller parameters changing with the tip addition account for the response in below rated operation, by scaling the rpm-generator torque quadratic gain (K opt), and varying the fine pitch setting (pitch opt).

Pre/post-processing
For every design evaluation loop, the HAWC2 case files are pre-processed, executed and post-processed on a single CPU. The top-level process is shown in Fig. 4. In the pre-processing MATLAB script, the baseline HAWC2 input files are modified in order to generate each tip extension design case. In the HAWC2 model, 10 additional structural and aerodynamic sections are added on the new part of the blade beyond 100% and the sections between 97.5%-100% are modified. The rpm-generator torque 85 quadratic controller gain and the fine-pitch setting are also modified. A case folder with all the HAWC2 input is assembled from all necessary files.
In the post-processing MATLAB script, the output time series files of HAWC2 are processed, and performance statistics are extracted. For the optimization, the mean generator power and the ultimate blade root flapwise bending moment are extracted.

Surrogate based optimization setup
The SBO framework is setup based on the MATLAB code package MATSuMoTo (Müller, 2013(Müller, , 2014, which is the MATLAB Surrogate Model Toolbox for deterministic computationally expensive black-box global optimization problems with continuous, integer, or mixed-integer variables that are formulated as minimization problems. The SBO framework determines the 95 design variable sets and send them to the pre-processor to execute the HAWC2 cases, in parallel CPU processing. The general SBO algorithm works as follows: -Generate initial design sets 5 https://doi.org/10.5194/wes-2020-108 Preprint. Discussion started: 15 October 2020 c Author(s) 2020. CC BY 4.0 License. The objective function is a very important part of this study, since it determines, which direction in the design space the SBO takes by evaluating new design variable sets. The objective function is defined as a weighted sum of the mean generator power and the ultimate blade root flapwise bending moment. Since we do not pursue any purely load alleviation driven designs but load-neutral power increase designs, the objective function is based only on the maximization of power when the loads are neutral or negative compared to the baseline. When loads are higher than 2% (an empirical limit accounting for model 110 uncertainty), the objective has a 90% weight on loads and 10% on power. A smooth Gaussian filter is used for the transition between neutral and higher loads (Fig. 5).

Surrogate modelling
For generating the initial sample set, MATLAB's Latin Hypercube design is used, with the maximin option and 20 iterations.
The minimum sample size is used which is 3*d+1, where d is the number of design variables, in our case 11. In the initial set, 115 a reduced cubic polynomial regression model is fitted. The choice of the surrogate model is decided based on prior studies of accuracy, comparing it with quadratic regression polynomials and radial basis functions.

Optimization
Using the fitted surrogate model on the initial set, a global optimization approach is followed, utilizing MATLAB's genetic algorithm with default settings. The best performing design point is chosen for a HAWC2 evaluation, together with points 120 created by randomly perturbing the best point found so far. Also a set of points that is uniformly selected from the whole variable domain is generated (using again a Latin Hypercube design) and the score is calculated over both sets of points.
Hence, it is possible to improve the global fit of the surrogate model and new areas of the variable domain where the global optimum may be located can be detected. Based on the available number of CPUs, 20 iterations are chosen, resulting in a total number of 174 HAWC2 evaluations, including the initial sample set.

Results
The progress of the optimization and the results for the whole set of evaluated design samples is discussed here. The characteristics of the best converged design are also discussed in detail.

Optimization results
The progress plot showing the best value of the objective function during the evaluation of each sample is shown in Fig 6. 130 It is seen that the objective function value is improved considerably from the starting samples and practically converges after  The best design in terms of the minimum value of the objective function comprises a tip extension with a length close to the limit of the defined length (7%), with all 11 design variables listed in Table 2. The design is shown to be a slender, backwards swept and highly upwind prebent tip shape, with fore positioning of the shear center, lower rotor speed and higher pitch settings. In Fig. 8, the 3D geometry of the baseline and optimized blade tip is compared.
The blade centerline of the optimized design is compared to the baseline in Fig. 9, where the sweep and prebend offsets

145
The performance of the optimized design is evaluated in terms of AEP in its IEC wind class I, and the lowest average wind speed class III. The 'clean' power curve is defined by steady uniform wind speed inflow from cut-in to cut-out with 1m/s steps, with no wind shear or turbulence. The higher fidelity aerodynamic module MIRAS is also used to run the same cases for comparison. The results are shown in Table 3, with the power curves plotted in Fig. 12. We see that MIRAS over-predicts the AEP for the baseline with around 1% deviation, and under-predicts the increased AEP due to the tip by around 1-2%.

150
The performance of the optimized design is also evaluated in the DLC1.3 (ETM) case which is used in the optimization, performed with the NW method against two different fidelity models. The BEM model implemented in HAWC2 is used as the lower fidelity method and the lifting line aerodynamic module implemented in MIRAS, which is employed as the higher fidelity solver. To ensure a meaningful comparison between the solvers, simulations of the extreme turbulent case have been carried out as follows: Firstly, to run a free turbulent simulation in MIRAS, the velocity defined Mann turbulent box used in 155 the optimization procedure has been transformed into a particle cloud by computing the curl of the velocity field. This cloud is slowly released one diameter upstream the turbine and it develops as it convects downstream towards the rotor plane. The released turbulent particles interact freely with the turbine wake. Vortex simulations with and without the turbine are performed.
In the simulation without a turbine, the local velocities are extracted every time step at the rotor plane position, in a 64x64 mesh with a cell size of approximately 6.25m. These velocities differ from the initially defined turbulent box due to the downstream 160 development of the flow in MIRAS. Such velocities are used to generate a new turbulent field which will be loaded in the HAWC2-NW and HAWC2-BEM simulations. Such turbulent field will mimic the turbulence seen by the turbine in MIRAS,  although the presence of the turbine and its wake will modify the turbulent field and such phenomena can not be accounted for.
In order to reduce uncertainties related to the turbine control in HAWC2, both the azimuthal rotor position as well as the pitch angle for each one of the blades are forced to be the same as the ones computed in the MIRAS simulation (with turbine), as 165 shown in Figures 13. In order to have a smooth start of the HAWC2-BEM and HAWC2-NW simulations, the first 60 seconds of the rotational speed signal from MIRAS have been modified using a hyperbolic tangent function. Differences between the simulations are therefore mainly related to the wake and flow modelling. The visible rotor speed differences (left plot of Figure   13) appear between the baseline rotor and the rotor with optimized tip, and the different fidelity levels clearly operate at the same rotor speeds from 100 seconds simulated time. The pitch angles also agree well between fidelity levels. The main offset 170 between the baseline and extended rotors is the minimum pitch angle, that is reduced to -1 degree by the optimization routine, see Table 2.
General statistics of the aerodynamic power and the flap-wise root bending moment are presented in Table 4. Regarding the power, it seems like the BEM method is slightly closer to the LL predictions in both mean and standard deviations for both rotor designs. However, regarding the MxBR it is for all quantities except of the min predicted value of the optimized blade 175 that the NW model is closer to the LL calculations. This is specially remarkable when looking at the standard deviation, where     In order to study the influence of the turbulence level in the results, the extreme turbulence level of the DLC 1.3 (40%) has been down-scaled to obtain inflow fields with a range of turbulence intensities from 0 to 40%. Statistics of the difference in the standard deviation of the signal are small for the optimized tip at low TI and larger for the baseline, as the TI increases the differences increase and align for both rotors. A similar picture is observed in the behaviour of the minimum root bending moment values, although generally differences between the rotors increase with the TI. In terms of the mean values, differences respect the LL predictions grow with increasing TI, with the engineering models predicting a higher moment at low TI and a lower one at high TI.  A detail analysis is carried out for the 40% turbulent case. Figure 16 shows the time signal of the aerodynamic power for both blade designs with the three different fidelity models. Generally there is a good agreement between the three solutions.
However, there is a small power offset, where the LL predictions are most of the time slightly larger than the BEM and NW predictions.
The time variation of the predicted flap-wise root bending moment is shown in Figure 17. Opposite to what is observed 195 for the power signal, there is no clear offset between the LL and the NW/BEM model predictions. It is visible that the BEM calculations experience larger high frequency variations compared to the higher fidelity models, exhibiting a larger standard deviation as it has previously been shown.
13 https://doi.org/10.5194/wes-2020-108 Preprint. Discussion started: 15 October 2020 c Author(s) 2020. CC BY 4.0 License. The radial distributions of mean and standard deviation of AOA for the baseline blade (solid lines) and extended blade (dashed lines) are shown in Figure 18. It can be seen that the mean AOA on the tip increases by roughly 4 degrees at the start 200 of the tip extension, which is mainly due to the twist distribution, see Figure 10. The MIRAS and the NW results both show decreased mean AOA inboard of the tip and an increased AOA on the tip itself. This is likely due to a load redistribution due to two factors: the offset of the trailed vorticity at the very blade tip and the velocity induced by the curved bound vorticity on the tip extension, Li (2018). In the standard deviation of the AOA (right plot of Figure 18) it can be seen that all codes predict the standard deviation to increase on the tip extension compared to the outboard part of the baseline blade. This is partly because 205 the sweep angle reduces the fraction of the relative velocity that is in the planes of the aerodynamic sections, while it does not affect the wind speed changes due to turbulence. Thus the same turbulence will lead to larger variations in AOA on the swept extended tip part of the blade than on the straight outboard part of the original blade. It can also be seen in all fidelity levels that the AOA is varying less on the section between 80 and 95 meter radius on the extended blade than on the original blade,