Investigation and Optimisation of High-Lift Airfoils for Airborne Wind Energy Systems at High Reynolds Numbers

: The potential of airfoil optimisation for the speciﬁc requirements of airborne wind energy (AWE) systems is investigated. Experimental and numerical investigations were conducted at high Reynolds numbers for the S1223 airfoil and an optimised airfoil with thin slat. The optimised geometry was generated using the NSGA-II optimisation algorithm in conjunction with 2D-RANS simulations. The results showed that simultaneous optimisation of the slat and airfoil is the most promising approach. Furthermore, the choice of turbulence model was found to be crucial, requiring appropriate transition modeling to reproduce experimental data. The k - ω - SST - γ - Re θ model proved to be most suitable for the geometries investigated. Wind tunnel experiments were conducted with high aspect ratio model airfoils, using a novel structural design, relying mostly on 3D-printed airfoil segments. The optimised airfoil and slat geometry showed signiﬁcantly improved maximum lift and a shift of the maximum power factor to higher angles of attack, indicating good potential for use in AWE systems, especially at higher Reynolds numbers. The combined numerical and experimental approach proved to be very successful and the overall process a promising starting point for future optimisation and investigation of airfoils for AWE systems.


Introduction
Airborne wind energy (AWE) is a relatively novel technology in renewable energies, aiming at the utilisation of stronger and more constant winds in altitudes not accessible with conventional wind turbines. A common feature of all AWE systems is the use of some form of tethered aircraft to reach altitudes between 200 and 1000 m [1]. The key advantages of AWE systems include the utilisation of additional wind resources at sites where wind energy harvesting is not feasible with conventional wind turbines. A significantly higher capacity factor can be achieved due to more constant winds in higher altitudes resulting in better suitability for base-load generation. Furthermore, the materials needed for a single system can be reduced by up to 90% compared to a modern wind turbine and consequently a substantial reduction in levelised cost of energy can be achieved [2]. Typical challenges when designing and operating AWE systems include the added complexity and dependence on a robust and reliable control system [3] as well as high-performance materials. Furthermore, existing regulations have to be revised to accommodate the operation of high-altitude tethered aircraft for power generation and more research is needed on the social acceptance of AWE systems [4].
Although all AWE systems share the general setup of a tethered aircraft of one form or another, a vast number of fundamentally different concepts exist, which are developed by a number of companies and universities. An excellent overview and classification of existing approaches can be found in [5,6].
The focus of the present work is the development and investigation of airfoils for ground-gen (pumping-kite generator) systems, operating in crosswind flight mode and Wind 2023, 3 using rigid wings, which have been developed by a number of companies, e.g. EnerKite, TwingTec or Ampyx Power. Figure 1 illustrates the operation of such a system, which consists of three distinct phases. During the traction (reel-out) phase, the kite flies in repeating patterns while reeling out the tethers from a winch connected to a generator and thus generating electrical energy. Figure-eight and circular trajectories are most common here, both resulting in the same performance under the same boundary conditions [7]. Systems relying on multiple tethers usually operate using figure-eight patterns to avoid twisting. At the end of the traction phase, the kite has to be retracted. To keep the energy consumed during the retraction phase to a minimum, two main approaches exist [8]. While soft kites are usually brought to an elevation of close to 90 • to retract without crosswind motion, rigid wings can be pitched down to low or negative angles of attack to more rapidly retract. The third phase is take-off and landing of the kite. A comprehensive summary of different approaches for rigid wings and soft kites can be found in [8],while an excellent overview of the general operation of different AWE systems is given in [9]. Besides control aspects, which were not in the scope of the presented work, the fundamental challenge when designing a rigid wing for power generation is covering the three distinct operating phases with partially contradicting aerodynamic requirements for optimal power yield. During take-off and landing, a high glide ratio C l C d is required (C l : lift coefficient; C d : drag coefficient), while during retraction a low glide ratio is beneficial. The deciding factor for a high yield during the power generation (reel-out phase) is the power factor Cl 3 C 2 d [11]. When considering total energy yield, additionally, a high-lift coefficient C l is required to generate high traction forces, while keeping the wing area and therefore the weight of the system low.
When designing wings for power generation, these aerodynamic criteria already greatly affect the choice of 2D airfoil geometries. The rather unique requirements for high lift and power factor at high Reynolds numbers, to allow for scalability of the AWE system to larger wings and therefore higher yields per single system, somewhat limit the pool of suitable airfoils. Although there are several airfoils with high lift available, e.g. S1223 [12] or FX 63-137 [13], they are usually not optimised for high Reynolds numbers required for large-scale AWE wings and power factors at high angles of attack. Attaining all the aforementioned requirements, including those for take-off, landing and retraction, can potentially be achieved via the extension to high-lift configurations, i.e., the addition of a slat.
The S1223 was picked as a starting point, since it already meets some of the aforementioned requirements quite well. It has also already been investigated for use on unmanned aerial vehicles [14]. A high-lift coefficient of C l = 2.2 is achieved at relatively low drag, and preliminary low-fidelity simulations using XFOIL [15] showed similar potential for higher Reynolds numbers. While featuring moderate stall characteristics, the concave pressure recovery with aft loading additionally results in a high (negative) pitching moment [12,16]. This is generally considered favourable when designing rigid wings for AWE, since it allows the tethered aircraft to be pitched down to negative angles of attack during retraction.
This paper presents an initial assessment of the aerodynamic optimisation potential of airfoil-slat configurations for airborne wind energy. Simulation and experimental data will be presented for an optimised geometry. Furthermore, experimental and simulation data for the S1223 airfoil at high Reynolds numbers will be presented for baseline comparison.

Design Process
The airfoil and slat configuration developed here is to be used for an AWE wing that is as light and mechanically simple as possible. Therefore, in contrast to most conventional high-lift configurations, the developed slat shall be fixed, i.e., not retractable, without any moving parts and therefore effective over the whole flight envelope. Configurations of this kind can already be seen in prototype wings, e.g. by TwingTec. Fixed passive slats have been investigated for use on conventional wind turbines, although they are often applied to thick root airfoils here [17]. Other studies mainly focus on an increase in lift and airfoils not particularly suitable for use in AWE [18,19]. The focus of this study is the further improvement of a highly cambered airfoil using a fixed slat, with regard to the aforementioned aerodynamic requirements of AWE wings. Additionally, the scope will be limited to a "thin slat" with a uniform thickness distribution, to keep the manufacturing process as simple as possible.
The starting point for the optimisation is the well-known S1223 airfoil [12]. A thin slat was added and multiple optimisations were run at a design Reynolds number of Re = 1.5 × 10 6 , with high lift and power factor being the optimisation goals. The optimisation runs included geometrical variation of the slat and airfoil separately, as well as slat and parts of the airfoil simultaneously, to find promising candidates for further investigation. The process constitutes an initial attempt at the optimisation of airfoils for AWE and will be further adapted and refined in future work.
Geometry modeling and parametrisation, as well as optimisation were performed using CAESES, while the simulations in the optimisation loop were conducted using Ansys Fluent (see Section 2.3).

Parametrisation
The airfoil was described using the thickness distribution and camber line, which were modeled using B-splines. Each B-spline consists of fixed start and end points with four additional control points. The slat camber line was described by a B-spline with three control points, each of which was free to move. The thickness was fixed to 0.2% airfoil chord c and an elliptical leading edge with an aspect ratio of 0.75 was added. An illustration of the parametrisation is shown in Figure 2. The parameters were chosen to provide a reasonable balance between flexibility during the generation of new geometries and the number of evaluations, i.e., simulations, necessary to sufficiently cover the parameter space. Adding additional degrees of freedom to the slat has shown not to provide any significant improvements, so following some attempts with an added fourth point, it was kept at three control points.

Optimisation Process
Parametrisation, optimisation and data exchange were all handled inside CAESES. A schematic overview of the general process is given in Figure 3. Numerous optimisations were run to fine tune the process, including different settings, constraints and optimisation parameters.

Algorithm
The NSGA-II genetic algorithm [20] as implemented in CAESES was chosen, as it is well suited for constrained multi-objective problems and generally achieves good convergence for nonlinear problems. A study using NSGA-II in conjunction with XFOIL can be found in [21]. Furthermore, it enables parallel evaluation of up to the population size cases (depending on the computational resources available) and combines design-space exploration and optimisation. Mutation and crossover probabilities were set to 0.01 and 0.9, respectively. Depending on the number of active optimisation parameters, population sizes of up to 128 were run for at least 10 generations each. Since the slat obviously was not to collide with the airfoil but still needed a large enough parameter range to allow for diverse shapes and positions, many generated variants in one generation were immediately dismissed by the prescribed constraints; hence, larger initial populations proved to be advantageous.
Some of the most promising candidates of NSGA-II optimisations were then subjected to a tangent-search algorithm to seek further optimisation potential in the immediate vicinity.

Objectives
The objective for the NSGA-II algorithm implemented in CAESES was set to maximise the power factor Cl 3 C 2 d and the lift coefficient C l at a design Reynolds number of Re = 1.5 × 10 6 . The original plan included MSES for the simulation of multiple operating points up to maximum lift C l,max , but had to be altered to accommodate the significantly higher computational cost of 2D-RANS simulations (see Section 2.3). The airfoils presented here were evaluated at an angle of attack of α = 15 • , which was deemed to be a reasonable operating point to start with, since it marks the onset of stall for the S1223, but should be well within the operating range of an airfoil and slat configuration. Although optimisation at only one angle of attack is not considered optimal and to be seen as more of an initial assessment, since, for example, the actual C l,max of the configuration cannot be assessed explicitly during the optimisation, quite promising results could be achieved. In future iterations, with the experience gained during the process and by using an open source solver to avoid licence availability being a limiting factor, three angles of attack per variant will be evaluated to better cover the whole range of operating points.

Design Parameters
With all non-fixed control-points of airfoil and slat having two degrees of freedom each, a maximum of 22 parameters were available. Since such a high number of parameters results in an unfeasible amount of evaluations to cover the parameter space, several variants of reduced parameter numbers were assessed. Initial runs included the addition of a slat to a predetermined airfoil (6 parameters) and several attempts to optimise just the airfoil with a reduced parameter space. Since the performance of the starting point (S1223) was already considered to be quite good, modifying just parts of the airfoil was deemed a valid approach.
The airfoil and slat geometry described later on originated from an optimisation run with ten parameters; all three slat control-points (x and y coordinates), the third airfoil camber line control point (x and y coordinates), the fourth airfoil camber line control point (y coordinate) and the overall thickness factor.

Constraints
For structural reasons, the airfoil was constrained to have a maximum thickness of more than 15% chord. The slat thickness was kept at 0.2% chord. Furthermore, a number of geometrical constraints were imposed on the relative position of control points, in order to exclude nonsensical geometries, e.g. collisions of slat and airfoil or concave slats.

Numerical Setup for Optimisation Evaluations
The computationally efficient viscous-inviscid low-fidelity solver MSES [22] should originally be used, to allow for a maximum number of variants to be evaluated. Subsequently, the findings were to be validated using 2D RANS simulations. However, the thin slat turned out to be a major challenge for the low-fidelity code. The sharp leading edge of the slat and its orientation in relation to the free stream velocity were identified as the main issue, especially during the mesh generation. Even though it was possible to simulate some thin slat and airfoil configurations, nearly all new configurations required a significant amount of manual input and customisation. Hence, the process was found not to be feasible for the simulation of large numbers of automatically generated geometries of this particular type. Therefore, the decision was made to instead employ RANS simulations right away. Steady-state 2D simulations were set up for each evaluation, using the k-ω-SST turbulence model with the γ-Re θ transition model [23]. The choice of turbulence model will be discussed in Section 3.2.
Simulations were run on unstructured hexahedral meshes generated with Ansys meshing, consisting of about 160k cells, including multiple refinement areas and prism layers (15 at the slat, 25 at the airfoil) to ensure a dimensionless wall distance of y + < 1 over the whole surface. The solver used during the optimisations was Ansys fluent using the SIMPLE algorithm. Convective fluxes were modeled with a second order upwind scheme. Although the use of commercial CFD software enabled a relatively quick and easy setup and integration of the simulations into the optimisation framework (see [24]), it was not considered optimal, mainly because the number of parallel evaluations was limited by licence availability. For future works, a tool-chain only depending on the open-source packages gmsh and OpenFOAM for the simulations is to be implemented.

Preliminary Geometry
Out of several thousand variants created over a number of optimisation runs, the most promising candidates were evaluated for a range of angles of attack and compared. One of the best performing was then chosen as a preliminary geometry and subjected to more detailed numerical and experimental investigation. Figure 4 shows the newly generated geometry, designated EW01, in comparison to the initial geometry of S1223.

Numerical Setup
All simulations outside the optimisation were conducted using OpenFOAM. The results presented are from 2D, steady-state RANS simulations, using the k-ω-SST turbulence model with the γ-Re θ transition model [23], as implemented in OpenFOAM as kOmegaSSTLM. Convective terms were approximated using a second order upwind scheme, the solver uses the SIMPLE algorithm. The turbulence intensity at the inlet was set to Tu = 0.01, while the pressure was fixed at the outlet. The boundary layer at the slat and airfoil were fully resolved with y + < 1 over the whole surface in all simulations.

Mesh
The airfoils were placed at the centre of a circular domain with a radius r = 20c. The unstructured hexahedral mesh was generated using gmsh. The results obtained were found to be equivalent to those on similar Ansys generated meshes, with the added benefit of reduced meshing time and using open-source software. Refinement areas have been placed around the airfoil and slat (see Figure 5) and prism layers were added to the slat (15 layers) and airfoil (25 layers) to ensure y + < 1 and a full coverage of the boundary layers. A mesh independency study with the force coefficients as indicators resulted in the use of meshes with about 172k cells (airfoil and slat) and about 123k cells (airfoil without slat). Cell sizes are shown in Table 1. Further refinement did not result in relevant changes of the force coefficients over the whole operational range.

Turbulence Model
Since the flow field around high-lift configurations is rather complex, including strong adverse pressure gradients, complex wakes, separated areas and potential re-attachment, the right choice of turbulence model is crucial. As a decent general purpose model, k-ω-SST [25] seemed a reasonable choice here. However, the use of the standard formulation early on showed rather severe problems for higher angles of attack which had to be addressed. Compared to experimental as well as lower fidelity simulations, even though the data matches quite well for low angles of attack, a very early onset of stall could be observed for configurations with and without slat. Figure 6 shows lift and drag data for S1223 and EW01 obtained with the same setup but different turbulent models as well as experimental data from [12] for comparison.  Figure 6. Influence of transition model for S1223 (blue) and EW01 (green).
Different approaches for the mitigation of possible problems when using k-ω-SST for the simulation of airfoils exist in the literature. Matyushenko et al. [26] observed an overprediction of C l,max and suggest a modification of k-ω-SST model coefficients. Although significant improvements could be achieved regarding the lift prediction, the modification resulted in a severe underestimation of drag and hence is not suitable for the intended prediction of the power coefficient.
A more promising approach is to include a transition model to more accurately simulate the complex boundary layers typical for airfoils and high-lift configurations. Khan et al. [27] successfully used the k-k L -ω model to improve the predictive accuracy for a symmetrical airfoil. This was also tested for slat and airfoil configurations investigated in this paper, but resulted in a drastic over-prediction of the maximum lift (see Figure 6).
The best match with experimental data for both configurations was achieved using the k-ω-SST-γ-Re θ model [23], which was tested extensively for flat plates, airfoils and high-lift configurations [28] and symmetrical airfoils [29]. The influence of transition modeling for thin slat and airfoil configurations will be discussed in more detail in Section 5.2.
To obtain more insight into the post-stall behaviour of the EW01 configuration, additional IDDES simulations were conducted. However, since only the standard k-ω-SST model was available as the underlying RANS model, the prediction of the stall point without an appropriate transition model still suffered from the aforementioned issues, so the simulations will not be presented here.

Experiment Setup
To complement the simulations and confirm the findings, experiments were conducted for the S1223 airfoil and the optimised EW01 geometry. The S1223 served as a starting point for the optimisation and was therefore measured for baseline comparison. Furthermore, up-to-date published experimental data for S1223 are only available for low Reynolds numbers, e.g. in [12]. The Reynolds numbers in the experiments conducted here ranged from Re = 3 × 10 5 up to Re = 7 × 10 5 .

Facilities and Measurement Techniques
The experiments were conducted in the large wind tunnel (GroWiKa) at TU Berlin. The closed-circuit tunnel features a measurement section of 2 m by √ 2 m, in which velocities of up to 70 m s can be achieved with a free-stream turbulence intensity lower than 0.5%. The free-stream velocity is measured via the differential pressure at the nozzle.
The model airfoils were mounted vertically in the measurement section with splitter plates at both ends. The gap between the model and the splitter plates was covered by a 3D-printed skirt to ensure a minimal gap and pressure equalisation. The setup is shown in Figure 7.
Force and pressure measurement data were acquired using a National Instruments DAQ-system with a sampling frequency of f s = 1 kHz and 10,000 samples per measurement point.

Pressure Measurements
To measure the airfoil drag, pressure measurements were conducted downstream of the model using a wake-rake with 58 pressure tubes (see Figure 8). The pressure tubes were distributed over a total length of 500 mm, with distances between the tubes varying from 4 mm to 16 mm (higher resolution in the middle). The tubes were 100 mm long with inner and outer diameters of d in = 1 mm and d out = 1.6 mm. The wake-rake was positioned 0.5c behind the airfoil trailing edge at α = 0 • and at the mid span of the model. Two Prandtl probes were positioned in the same plane in the free stream above the outermost tubes of the rake. To determine the difference in velocity and static pressure at these positions, an additional pressure sensor was installed between the static pressure taps of the probes. The static pressure of the left (suction side) Prandtl probe was also taken as reference pressure for the wake-rake differential sensors. The position of 0.5c was chosen as a compromise between spatial resolution and the resolution range of the pressure sensors, in order to minimise the differences in velocity and static pressures at the outermost positions of the rake.
The drag coefficient was calculated using Equation (1) [30], with P and P ∞ being the total pressure in the wake and free stream, p and p ∞ the respective static pressures, and using trapezoidal numerical integration.
HDO-Series pressure sensors by Sensortronics with a nominal load of ±1000 Pa and a maximum error of < 0.5% full scale were used. The measurement range of the pressure sensors set the upper limit of Reynolds numbers achievable with the setup.

Model Airfoils
The dimensioning of the wing models results from two conflicting design objectives. On the one hand, especially for high cambered airfoils, a significant influence of threedimensional effects due to corner vortices on the suction side of the outer regions of the model had to be expected. To minimise this influence, a high aspect ratio is required, which in turn results in models, especially with slim high-lift airfoils, being pushed to their structural limits rather quickly and thereby limiting the achievable Reynolds numbers. Since model construction using GFRP/ CFRP or machined metal is rather time-consuming and very costly, a novel structural design was developed to drastically reduce costs and resources needed, while still providing sufficient mechanical stiffness for investigations at high Reynolds numbers.
The principal construction is shown in Figure 9. The model mainly consists of 119 mm spanning, 3D-printed plastic (PET-G) segments, manufactured using the FDM method. The segments were mounted on a stainless steel profile with a square cross section. In between the segments, 1 mm thick laser-cut stainless steel sheets were placed to provide additional stiffness. Airfoil segments and steel sheets were fixed in place by two metal bolts and glued together. For airfoil and slat configurations, a thin CFRP slat was glued to the steel sheets. The models consist of ten airfoil segments with a total chord length of c = 300 mm, resulting in a span of s = 1200 mm and an aspect ratio of AR = 4. After assembly, the models were sanded, primed and painted to provide a smooth surface.
Of the two models investigated in this paper, S1223 was constructed using the described new process, while EW01 was built as a more traditional, full CFRP model. Both models showed virtually no deformation at the highest Reynolds numbers investigated.

Results
Experimental and simulation data for S1223 and EW01 will be presented and compared in the following section. The lift measurements presented were obtained via force measurements, while the drag data originate from the pressure measurements in the wake. Experimental drag data are therefore only shown for angles of attack lower than that of C l,max .

Baseline Airfoil (S1223)
Experimental and numerical lift and drag coefficients over the angle of attack α are shown in Figure 10 for Re = 3 × 10 5 and Re = 5 × 10 5 , as well as simulation data for Re = 1.5 × 10 6 , which was the design Reynolds number for the optimised slat and airfoil geometry. Additionally, the lift coefficients measured by Selig and Guglielmo [12] at Re = 2 × 10 5 are shown for reference. The experimental lift polars are relatively close to the data obtained by Selig and Guglielmo. The slightly higher overall lift and difference in slope might not necessarily be an effect of the slightly higher Reynolds number, but may be explained by suspected stronger three-dimensional effects at higher angles of attack in the original experiments [12]. This seems to be a reasonable assumption, since the airfoil aspect ratio was only 2.8, compared to 4 in the experiments presented here. In the range measured, the lift shows no strong Reynolds number dependency. However, the Reynolds numbers measured were relatively close to each other. The overall drag is slightly increased for the higher Reynolds number.
The match between simulation and experimental data in general is considered to be quite good. The prediction of the overall trend over the whole range, the pressure side stall between α = − 3 • and α = − 4 • , as well as the position of C l,max were accurately predicted. The simulations show a higher dependency on the Reynolds number than found in the experiments. This might be a result of slightly different prediction of the transition onset, possibly caused by the non-zero roughness of the real-life model. The drag coefficients match the experimental data as well as could be expected using RANS simulations, even though a slightly different slope and lower drag at small angles of attack were predicted. Overall, the drag polars do not differ as much from each other, compared to the experimental data.
In the range measured and simulated, S1223 does not exhibit excessive Reynolds number dependency regarding lift. The lift for the highest Reynolds number simulated (Re = 1.5 × 10 6 ) represents the best match to the measured data at lower Reynolds numbers, so no further increase in maximum lift is to be expected. The simulations fail to capture the increase in drag with increasing Reynolds number.

Optimised Airfoil with Slat (EW01)
As pointed out earlier, accurate simulations without an appropriate transition model turned out to be challenging at best and nearly impossible for the slat and airfoil geometries investigated here. The cause of the inability of the standard k-ω-SST model to predict the stall points of configurations investigated in this work deserves some attention. The optimised EW01 will serve as an example here.
The main issue is the predicted premature onset of stall and hence underestimation of C l,max when using plain k-ω-SST (see Figure 6). When comparing lift polars from simulations with and without the γ-Re θ transition model (dark-green and light-green graphs in Figure 6), the first major deviation occurs between α = − 3 • and α = − 4 • . Here, the sudden increase in lift and decrease in drag is not captured by the standard k-ω-SST simulation. This "dent" in the polars was generally confirmed by the experiments conducted (see Figure 11). Figure 12 shows the turbulence intensity for simulations with and without γ-Re θ transition model (see Figure 6) at critical angles of attack. For the γ-Re θ model, the approximate point of transition of the boundary layer on the main airfoil is marked with a dashed line. The sudden increase in lift and decrease in drag between α = 5 • and α = 6 • seems to be caused by a downstream jump of the transition point between those angles of attack. Consequently, the size of the area with small scale separation on the rear of the airfoil decreases, starting at α = 6 • . This shift could not be predicted by the pure k-ω-SST model and therefore leads to an increasingly separated flow on the rear part of the airfoil for increasing α. While at low angles of attack the predictions are roughly in agreement, the thicker and more developed turbulent boundary layer here does not withstand the adverse pressure gradient as well as the transitioned one and hence leads to vast separated areas for α ≥ 9 • .
The earlier onset of transition for α < 6 • on the main airfoil may additionally be triggered by the large separated area at the slat pressure side for low angles of attack. Note that the size of this area is also slightly overestimated by the simulations without the γ-Re θ transition model. Figure 11 shows experimental and numerical lift and drag coefficients over the angle of attack α for Re = 5 × 10 5 and Re = 7 × 10 5 , as well as simulation data for Re = 1.5 × 10 6 (design Reynolds number for EW01).
The experimental polars are quite close to each other for the Reynolds numbers measured, with the drag polars being nearly identical. The lift for higher angles of attack is marginally higher with increasing Reynolds number. The aforementioned "dent" between α = 5 • and α = 6 • can be observed for both polars, although slightly less pronounced than in the simulations.  The simulation data, again, show a higher degree of Reynolds number dependency than suggested by the experiments. The maximum lift C l,max and α stall increase significantly with higher Reynolds numbers, while the drag polars, again, only show slight differences for small α. The effect of the transitioning boundary layer mentioned above seems to be more pronounced at lower Reynolds numbers, which could be expected. Furthermore, the stall at negative angles of attack is predicted to occur more rapidly than suggested by the experiments. Apart from the obvious limitations of RANS simulations, a highly probable source of deviations between experimental and simulation data is the influence of the slat mounting points, which could not be accounted for in 2D simulations.
The velocity field at Re = 1.5 × 10 6 (black graph in Figure 11) for selected angles of attack is illustrated in Figure 13. The upper end of α until which the airfoil is in full pressure side stall is marked by α = − 3 • , which is a desirable state for efficient retraction of an AWE wing. The flow then rapidly attaches and exhibits a large separated area between slat and airfoil, as well as high velocities near the main airfoil leading edge (see α = 0 • ); both steadily decrease with increasing α. The effects of the shift of the transition point between α = 5 • and α = 6 • described earlier can clearly be observed by the more attached flow at the rear part of the airfoil suction side from α = 6 • onwards. The separated area at the slat pressure side is decreasing further, but is still observable at the angle of attack the optimisations were run at (α = 15 • ). The flow stays mostly attached until α = 21 • and then starts to separate, starting from the trailing edge. Since only steady state RANS simulations were run and hence simulation data after stall should be taken with a grain of salt, it can be assumed that the consequent drop in lift happens much more gradually than might be expected based on the concave shape of the suction side of the main airfoil. The experimental data (see Figure 11) support this assumption. Figure 13. Velocity-field of EW01 at Re = 1.5 × 10 6 for selected α.

Discussion
Compared to S1223, being the optimisation starting point, the maximum C l at the design Reynolds number of Re = 1.5 × 10 6 could be increased significantly (by about 0.6), which is an improvement that one should expect when adding a slat to an airfoil. Figure 14 shows a comparison of power factors, which was the second optimisation objective, calculated from the simulation data. Even though the maximum power factor could only be increased for higher Reynolds number, the main benefit here is the shift of high power factors to larger angles of attack and therefore regions of higher lift. For AWE wings, this potentially means a higher absolute traction force at the point of maximum and therefore increased yield. Furthermore, the higher lift allows to compensate for additional drag forces, e.g. caused by three-dimensional effects or tether drag, when considering the power factor for the whole system. Additionally, the glide ratio C l C d drops to lower values significantly earlier when decreasing the angle of attack (starting at α = 5 • ), potentially reducing the amount of energy required during the retraction phase. EW01, Re = 500k Re = 700k Re = 1.5m S1223, Re = 300k Re = 500k Re = 1.5m Figure 14. Power factor for S1223 and EW01.

Summary and Conclusions
The potential of airfoil optimisation to meet the specific requirements of AWE systems was investigated and joint experimental and numerical investigations of the S1223 airfoil and an optimised airfoil and slat configuration were conducted.
Starting from S1223, airfoil and slat configurations were generated, mainly using the NSGA-II optimisation algorithm, where the slat was constrained to a "thin slat" with a uniform thickness distribution. The use of lower fidelity methods proved not to be feasible for this specific set of requirements, so the evaluations had to be conducted using 2D RANS simulations. Several variants of optimising slat and airfoil at a time or both simultaneously were investigated. The simultaneous optimisation of slat and airfoil proved to be the most promising approach and one of the best candidates was picked for further investigation.
During preliminary numerical investigations, it was found that even at high Reynolds numbers, simulations using the k-ω-SST turbulence model provided unsatisfactory results (Section 3.2). Results best in line with experimental data could be achieved using k-ω-SST with transition modeling using the γ-Re θ model. Using this approach, it was found that for thin slat and airfoil configurations, even at high Reynolds numbers, the effects of transitioning boundary layers play a significant role in the accurate prediction of separation and hence should always be included (see Figure 12).
More detailed 2D RANS simulations for a range of operating points were conducted at three Reynolds numbers using OpenFOAM for the S1233 airfoil and the airfoil and slat configuration designated EW01.
Experimental investigations for both airfoils were conducted in the large wind tunnel (GroWiKa) at TU Berlin. To minimise three-dimensional effects, a new vertical test rig was developed, allowing for an aspect ration of AR = 4 and Reynolds numbers up to Re = 1.2 × 10 6 (depending on the airfoil). Lift forces were acquired using force measurements while drag was determined using pressure measurements in the wake behind the airfoils.
A novel design for model airfoils with sufficient stiffness for high Reynolds number investigations and greatly reduced manufacturing cost was developed, relying mostly on 3D-printed airfoil segments. The airfoils manufactured using this method proved to be competitive with respect to full CFRP models.
Experimental and simulation data were compared. For the S1223, investigated for reference, the measured lift matches experimental data available from the literature reasonably well and is suspected to be slightly more accurate due to the setup allowing for a higher aspect ratio of the airfoil. Experiments were conducted for higher Reynolds numbers than available until now (Re = 3 × 10 5 and Re = 5 × 10 5 ), as well as simulations for those Reynolds numbers and Re = 1.5 × 10 6 . The simulation data for C l matches the experiments reasonably well, although a slightly higher Reynolds number dependency is predicted at lower Reynolds numbers. Numerical drag data, although reasonably close to the experimental data, struggle to predict the slight Reynolds number dependency observed in the measurements.
The optimised EW01 configuration was investigated thoroughly using experiments and simulations. It shows increased maximum lift and a shift of the maximum power factor to higher angles of attack, as well as an increased power factor at Re = 1.5 × 10 6 , compared to the S1223. Again, an increased Reynolds number dependency is predicted by the simulations, compared to the experiments, likely resulting from the influence of the slat mounting points not being considered in the simulations.
The optimised airfoil and slat geometry shows a significantly improved C l,max , as well as improved power factor at high Reynolds numbers and high angles of attack and therefore good potential for use in AWE systems. The methods employed for optimisation produced promising results, although the process will be further improved in future work. Using a "thin slat" proved to be a viable approach, not only to reducing manufacturing effort, but also regarding aerodynamic performance. When simulating such configurations using CFD, proper modeling of boundary layer transition should always be taken into account. Overall, the combined numerical and experimental approach to quickly and accurately simulate optimisation candidates and if needed conduct experiments for validation using relatively inexpensive models proved to be very promising for future optimisation and investigation of airfoils for AWE systems. Funding: This research was funded by the German Federal Ministry for Economic Affairs and Climate Action, project "EnerWing_xM", grant number 0324355E.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors would like to thank Friendship Systems AG for providing CAESES licenses.

Conflicts of Interest:
The authors declare no conflict of interest.