Computational aerodynamics for leading-edge inflatable kite design

Soft-wing kites are morphing, bridled, tensile lifting surfaces used for wind-assisted ship propulsion and airborne wind energy applications. Their swept-back planform, pronounced anhedral, and unconventional leading-edge geometry induce complex aerodynamic behaviour that challenges conventional modelling approaches. For leading-edge inflatable (LEI) kites, pressure-side separation induced by the inflated tubular leading edge renders classical inviscid methods insufficient, thereby necessitating sectional input from higher-fidelity approaches.

This study presents and applies a computationally efficient aerodynamic framework to an LEI kite by coupling a vortex step method (VSM) with RANS (Reynolds-averaged Navier–Stokes)-derived airfoil polars validated against wind tunnel measurements. The RANS simulations were used to train a machine learning surrogate model to facilitate parametric design studies. The application of machine learning to LEI kite aerodynamics is novel, and the surrogate model achieves R²> 0.98 across the considered parameter space.

Three-dimensional load predictions for the TU Delft V3 LEI kite were evaluated against wind tunnel data and reference three-dimensional RANS simulations. Within the operational incidence range $\mathit{\alpha }\in [-\mathrm{1},\mathrm{10}]$°, the predicted lift and drag agree with measurements to within 9 % and 13 %, respectively. Across this range, the framework reproduces the measured aerodynamic trends more consistently than the reference 3D RANS results, while reducing the computational cost by several orders of magnitude.

A rigid-body stability analysis indicated static stability in roll, pitch, and yaw, as well as aerodynamic damping in all three axes within the quasi-steady model.

Parametric analyses revealed inherent trade-offs between aerodynamic efficiency and stability, motivating the adoption of multi-objective optimization strategies.

The validated framework provides high predictive accuracy at low computational cost and forms a foundation for rapid aerodynamic analysis, stability assessment, design optimization, and aero-structural coupling in the conceptual and preliminary design phases.

1 Introduction

Soft-wing kites have expanded beyond their original use in marine sports to renewable energy applications, most notably airborne wind energy (AWE) and wind-assisted ship propulsion (WASP). A soft kite is a tensile structure consisting of a lightweight membrane wing supported by a bridle line system that transfers the aerodynamic loads and serves as the actuation mechanism (Oehler and Schmehl, 2019). These membrane wings fall into three types: ram-air wings, with two stacked membrane layers connected by membrane ribs and inflated by the incoming flow; leading-edge inflatable (LEI) wings, with a pre-inflated tubular frame supporting a single membrane layer; and single-skin wings, with a single membrane layer. Soft-wing kites are morphing aerodynamic control surfaces that deform through aeroelastic coupling (Breukels et al., 2013) and in response to actuation-induced changes in the bridle line system geometry (Hummel et al., 2019; Elfert et al., 2024).

Because the bridle line system transfers the distributed aerodynamic load of the membrane wing directly into tension, soft-wing kites achieve high lift-to-weight ratios in crosswind operation. For example, a 14 m² LEI kite sustained a controlled tether force of 3 kN during crosswind figure-eight manoeuvres, corresponding to a load factor of approximately 20 relative to the 15 kg airborne system mass and 60 relative to the 5 kg wing alone (Van der Vlugt et al., 2019).

Design approaches common for sports kites rely heavily on experimental testing to determine the aerodynamic performance. Likewise, in AWE and WASP, control strategies rely on empirical measurements of control responses to inform gain tuning. When scaling kites to sizes required for commercially meaningful wind energy harvesting, such reliance on experiments becomes prohibitively costly, restricting iterative development cycles and increasing risks.

Instead, numerical simulation tools can be used effectively to assess kite designs, control responses, and stability derivatives, and to optimize flight paths in a virtual environment, without reliance on extensive prototype testing. To realize the full potential of these tools, fluid–structure interaction (FSI) phenomena must be accounted for by modelling the two-way coupling between aerodynamic loads and structural deformation (Candade et al., 2020). This FSI modelling has typically combined aerodynamic models with high-fidelity structural solvers through iterative coupling (Bosch et al., 2014; Lolies et al., 2019; Folkersma et al., 2020; Thedens, 2022). In these partitioned frameworks, the aerodynamic and structural solvers exchange loads and displacements in successive iterations until convergence. These iterative procedures make computational speed a critical factor, particularly in the context of design and optimization.

Simulating the aerodynamics of soft-wing kites is challenging owing to their complex geometry, which are characterized by a swept-back planform, low aspect ratio, and pronounced anhedral. A recent review synthesized the state of the art in AWE aerodynamics, spanning prototype architectures, and numerical and experimental methods from potential-flow models to high-fidelity Reynolds-averaged Navier–Stokes (RANS) and large-eddy simulation (LES) approaches (Castro-Fernández et al., 2026). Within this context, 3D RANS computational fluid dynamics (CFD) simulations of LEI kites have been conducted to explore variations in Reynolds number, angle of attack, and sideslip angle (Deaves, 2015; Viré et al., 2020, 2022). These studies, performed with k-ω shear stress transport (SST) turbulence and γ-Re_θ transition models, revealed characteristic separation zones and 3D crossflow effects typical of curved, low-aspect-ratio wings. Despite providing valuable verification data, the high computational cost of these simulations makes them unsuitable for iterative optimization.

Lower-fidelity methods have been developed, including lifting-line (LL) and vortex-lattice formulations adapted for curved planforms and low aspect ratios (Gaunaa et al., 2011; Leloup et al., 2013; Branlard et al., 2022; Duport, 2018), as well as non-linear vortex-lattice and panel-based formulations capable of capturing moderate flow separation (Kappel, 2012; Leuthold, 2015). Panel and vortex-lattice methods fail for LEI kites because they cannot accurately capture flow separation and subsequent reattachment, which dominate the pressure side of the wing throughout most of its flight envelope (Folkersma et al., 2019).

A proposed multiple-wake extension (Leuthold, 2015) partially mitigates this limitation but ultimately suffers from the same underlying problem: when the reattachment line intersects a separation line, the model cannot robustly resolve the local flow topology. Accurately maintaining flow tangency and impermeability in these regions requires prohibitively fine body and wake discretization, leading to numerical instability and computational costs that undermine the method’s intended efficiency.

Taken together, these considerations motivate the use of LL-based formulations as a computationally efficient basis for aerodynamic simulation, as they retain the dominant 3D induction effects while avoiding the prohibitive cost of high-fidelity CFD.

The concept underlying these formulations dates back to Prandtl (1918), who modelled a finite wing as a bound vortex with a trailing vortex sheet and computed the sectional lift from the Kutta–Joukowski theorem (Kutta, 1902), evaluated with the vector sum of the freestream and wake-induced velocities. Two later developments are central to the present work. First, Pistolesi (1932) showed that the flow-tangency condition is best enforced at the three-quarter-chord point, an insight that Weissinger (1947) embedded in a discrete horseshoe-vortex formulation with improved robustness for swept wings (Dorn and DeYoung, 1947). Second, beginning with Tani (1934) and formalized by Piszkin and Levinsky (1976), the strict no-through-flow condition of classical potential-flow theory was replaced by a circulation-matching approach, in which the sectional circulation is iterated to consistency with prescribed 2D airfoil polars, yielding a hybrid viscous–inviscid formulation that admits non-linear sectional aerodynamics, including stall. Barnes (1997) combined both elements in the vortex step method (VSM), which Damiani et al. (2019) introduced to the AWE field, and Cayon et al. (2023) extended it to aero-structurally coupled kite simulations.

The predictive capability of any such formulation, however, hinges on implementation details and on the availability of representative airfoil polars for each spanwise wing section.

The present study focuses on kites with LEI wings, although the general modelling approach can be applied to other types of kite. For conventional airfoils, aerodynamic polars are commonly generated using tools such as XFOIL (Drela, 2013), which couples an inviscid panel method to an integral boundary layer formulation in a viscous–inviscid interaction framework. For LEI kites, however, the inflated leading-edge tube induces persistent pressure-side separation and recirculation across the full range of angles of attack (Folkersma, 2022), thereby violating the thin, predominantly attached boundary layer assumptions underlying such methods. Consequently, inviscid solvers are unlikely to be suitable for predicting the 2D flow around such geometries, requiring higher-fidelity viscous simulations, typically based on the RANS equations (Folkersma et al., 2019). However, this considerably increases computational demand, especially in the context of iterative aero-structural coupling, where the VSM discretizes the wing into multiple deforming spanwise sections, each requiring its own polar data.

To reduce computational effort, Breukels (2011) conducted 2D RANS CFD simulations for a parameterized family of LEI airfoils and derived a correlation model for airfoil polars. In this formulation, C_l is represented by a third-degree polynomial in the angle of attack α, whereas C_d and C_m are represented by second-degree polynomials. The polynomial coefficients are expressed as functions of the airfoil thickness ratio t and camber parameter κ, facilitating the efficient interpolation of aerodynamic polars within the considered range of airfoil geometries. While substantially reducing computational cost, the approach introduces several limitations, including discontinuities arising from empirical extrapolation beyond the simulated angle-of-attack range. Watchorn (2023) performed new 2D RANS CFD simulations and reported systematic deviations with Breukels for comparable airfoils.

A fast model alone is not sufficient; numerical models require rigorous validation to earn industry trust and deliver meaningful impact. Although 3D RANS CFD data can benchmark lower-fidelity methods in incidence and sideslip (Duport et al., 2019), wind tunnel testing is preferred. To remove structural and aeroelastic effects, a rigidized subscale model of the TU Delft LEI V3 kite was tested with a six-component load balance and stereoscopic particle image velocimetry to quantify forces and flow fields (Poland et al., 2026a, b). In a separate campaign on LEI airfoils, a stethoscope, infrared cameras, zigzag tape, pressure strips, and a pitot-tube rake were used to measure loads, assess forced transition, and characterize the flow (van Lith et al., 2026).

This study investigates the aerodynamic modelling of soft-wing kites using the V3 kite as the reference system, leveraging its openly available datasets and further establishing it as a community benchmark. The aerodynamic behaviour was simulated using an extended VSM. To provide reliable predictions at low computational cost, machine learning (ML) models were trained using 2D RANS CFD simulations. The substitution of potential-flow-based surrogate models with an RANS-trained ML model was evaluated against flat-plate data and 2D wind tunnel measurements (van Lith et al., 2026). The influence of the source of 2D input polars on 3D performance was assessed by comparison with wind tunnel experiments (Poland et al., 2026b) and 3D RANS CFD simulations (Viré et al., 2022). The study further examines VSM convergence characteristics, establishes recommended best practices, and derives rigid-body stability derivatives. Finally, kite design considerations are discussed by analysing the effects of variations in airfoil shape, planform configuration, and tow-point position.

The remainder of this paper is organized as follows. Section 2 describes the VSM, Sect. 3 introduces the 2D RANS CFD setup, and Sect. 4 introduces the deployed ML model generation and evaluation for 2D LEI airfoils. The results are discussed in Sect. 5, and design considerations are provided in Sect. 6.

2 Vortex step method

The implemented VSM largely follows Cayon et al. (2023), which builds on the formulations of Ranneberg (2015), Damiani et al. (2019), and Jonkman (2021). The implemented method assumes steady flow around the lifting surface, such that the solution is restricted to the spatial domain and the starting vortex, presumed far downstream, can be neglected. For circulation matching, the flow is decomposed into a near-field region in the vicinity of each spanwise section and a surrounding far-field region dominated by 3D induction effects. The effect of the flow in the near-field region is characterized by non-linear viscous airfoil polars that govern sectional aerodynamic lift, drag, and moment, while the far-field flow induced by the bound and trailing vortices is treated as inviscid, incompressible, and irrotational within the vortex filament framework. These regions are coupled at a single spanwise location per section, namely the control point, where the Kutta–Joukowski theorem relates sectional circulation to the integrated lifting force. The VSM computes the aerodynamic properties of a given wing through a number of steps, illustrated in Fig. 1 and described in order in Sect. 2.1–2.5.

Figure 1 Flowchart of the VSM, integrating the table-of-contents structure. The three stacked layers represent preprocessing (top layer), iterative circulation matching (middle layer, Sect. 2.4), and postprocessing (bottom layer). After geometry discretization and initialization, the spanwise circulation distribution Γ is solved iteratively until convergence, after which the sectional and global aerodynamic coefficients are evaluated.

2.1 Wing geometry discretization

To ensure methodological consistency, the geometry of the TU Delft V3 kite was adopted from previous CFD studies (Viré et al., 2020, 2022) and wind tunnel experiments with a rigid scale model (Poland et al., 2026a, b).

Figure 2 Photograph of the TU Delft V3 kite taken in 2012, with wing panels superimposed for illustration purposes. The quarter-chord and three-quarter-chord points for each panel are indicated. The wing reference frame is defined with the x axis aligned with the mid-span chord line, the y axis oriented towards the left tip (when looking from the front) parallel to the LE tube at mid-span, and the z axis oriented in the direction of the bridle lines.

Representative airfoils were extracted as chordwise slices through the wing (see Fig. A1 in Appendix A). The discrete spanwise locations lie on the spanwise reference line defined by connecting the quarter-chord points, obtained at each station as

$$\begin{array}{}\text{(1)}& {\mathit{x}}_{c/\mathrm{4}}={\mathit{x}}_{\mathrm{LE}}+\mathrm{0.25}({\mathit{x}}_{\mathrm{TE}}-{\mathit{x}}_{\mathrm{LE}}).\end{array}$$

At each location, the section plane contains the global x direction and is oriented in such a way that its normal is perpendicular to the local surface normal at ${\mathit{x}}_{c/\mathrm{4}}$.

A formulation with planes aligned with the global x axis was selected rather than one oriented locally normal to the quarter-chord line. Although both formulations are equivalent in the linear lift regime, the latter introduces additional geometric non-linearities once C_l(α) departs from linear behaviour, which may bias the predicted load distribution (Fritz, 2024).

The wing is discretized into n_p panels, each constituting a fundamental element of the solver. The number of panels is not required to coincide with the number of chordwise slices used to define the geometry. For instance, when adopting a uniform panel distribution, both geometric properties (leading-edge (LE) and trailing-edge (TE) coordinates) and aerodynamic properties (2D polars) defined at the chordwise slices may be interpolated to an arbitrary value of n_p.

Because the panels are rectangular, they do not form a geometrically continuous surface representation. Instead, the panel system serves to define the bound and trailing vortex filaments, rather than to reproduce the exact wing geometry.

2.2 Vortex filament system

Following Kelvin’s circulation theorem and Helmholtz’s vortex theorems, all vortex elements must form a closed system; vortices cannot spontaneously appear or vanish within the fluid (Batchelor, 1967). Accordingly, the bound vortex filament, placed at the quarter-chord line of each panel, is connected to a trailing vortex filament that extends in the chordwise direction to the TE (see Fig. 3). At the TE, these are attached to wake filaments, which extend to infinity in the freestream direction u_∞, in accordance with the frozen wake assumption (Ranneberg, 2015). Aligning the wake with the undisturbed free stream maintains consistency with the linearized LL formulation (Phillips and Snyder, 2000). Allowing the wake direction to depend on the locally induced velocity would introduce geometric non-linearity and require a free-wake treatment; and for strongly swept or highly anhedral configurations, spanwise-varying wake directions may additionally lead to filament convergence or intersection in a straight, semi-infinite filament model.

At infinity, the wake filaments are closed by an additional filament, thereby completing the vortex system. This closing filament is not explicitly modelled, as its aerodynamic influence at infinite distance is assumed to be negligible. The resulting vortex structure is known as a horseshoe vortex. Helmholtz’s theorems further require that the circulation strength remains constant along the length of each vortex filament (Batchelor, 1967); therefore, each panel has a single circulation value.

Figure 3 Illustration of a single panel, showing the incoming freestream velocity u_∞ in blue and the circulation Γ around each filament in red. The local panel reference frame – defined by the axes x_p, y_p, and z_p – is indicated. The positions of the LE, TE, quarter-chord point $c/\mathrm{4}$, and three-quarter-chord point $c\mathrm{3}/\mathrm{4}$ are also shown. Additionally, the panel chord c and panel width w are displayed.

The induced velocity vector u_ind at an evaluation point x due to a finite, 3D-bound vortex filament or chord-aligned trailing vortex filament is computed using the Biot–Savart law (Damiani et al., 2019):

$$\begin{array}{}\text{(2)}& {\mathit{u}}_{\mathrm{ind}}={\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{{\mathit{r}}_{\mathrm{1}}\times {\mathit{r}}_{\mathrm{2}}}{\Vert {\mathit{r}}_{\mathrm{1}}\times {\mathit{r}}_{\mathrm{2}}{\Vert }^{\mathrm{2}}}}\left({\mathit{r}}_{\mathrm{0}}\cdot \left({\displaystyle \frac{{\mathit{r}}_{\mathrm{1}}}{\Vert {\mathit{r}}_{\mathrm{1}}\Vert }}-{\displaystyle \frac{{\mathit{r}}_{\mathrm{2}}}{\Vert {\mathit{r}}_{\mathrm{2}}\Vert }}\right)\right),\end{array}$$

where Γ is the circulation strength, r₁ and r₂ are vectors from the filament endpoints to x, and r₀ is the vector along the filament segment.

For the wake filaments, which are modelled as semi-infinite vortex filaments, the induced velocity at x is given by Ranneberg (2015):

$$\begin{array}{}\text{(3)}& {\mathit{u}}_{\mathrm{ind},\mathrm{wake}}={\displaystyle \frac{\mathrm{\Gamma }}{\mathrm{4}\mathit{\pi }}}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{r}}_{\mathrm{1}}\cdot {\widehat{\mathit{u}}}_{\mathrm{\infty }}}{\Vert {\mathit{r}}_{\mathrm{1}}\Vert }}\right){\displaystyle \frac{{\mathit{r}}_{\mathrm{1}}\times {\widehat{\mathit{u}}}_{\mathrm{\infty }}}{\Vert {\mathit{r}}_{\mathrm{1}}\times {\widehat{\mathit{u}}}_{\mathrm{\infty }}{\Vert }^{\mathrm{2}}}},\end{array}$$

where r₁ is the vector from the filament origin to the evaluation point, and ${\widehat{\mathit{u}}}_{\mathrm{\infty }}={\mathit{u}}_{\mathrm{\infty }}/\Vert {\mathit{u}}_{\mathrm{\infty }}\Vert$ is the wake-direction unit vector.

A vortex core model was employed to avoid singular behaviour by reducing the induced velocity at points located within a specified core radius, ϵ. For bound vortex filaments, the core radius is defined as a fixed fraction of the filament length; here, 5 % was adopted (Garrel, 2003; Damiani et al., 2019). For trailing and semi-infinite wake filaments, a viscous core model is adopted (Bhagwat and Leishman, 2002), with the core radius given by

$$\begin{array}{}\text{(4)}& \mathit{ϵ}=\mathrm{2}\sqrt{{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{0}}\mathit{\nu }\mathrm{\ell }}{\Vert {\mathit{u}}_{\mathrm{\infty }}\Vert }}},\end{array}$$

where α₀= 1.25643 is the Oseen parameter, ν the kinematic viscosity, and ℓ the relevant filament length.

For bound and finite trailing filaments, core regularization is applied for $d_{⟂}<\mathit{ϵ}$ by evaluating the Biot–Savart velocity at the radially projected point x_proj on the core boundary and scaling linearly:

$$\begin{array}{}\text{(5)}& {\mathit{u}}_{\mathrm{ind},\mathrm{core}}\left(\mathit{x}\right)={\displaystyle \frac{d_{⟂}}{\mathit{ϵ}}}{\mathit{u}}_{\mathrm{ind}}\left({\mathit{x}}_{\mathrm{proj}}\right),d_{⟂}<\mathit{ϵ}.\end{array}$$

The perpendicular distance to the filament centreline is

$$\begin{array}{}\text{(6)}& d_{⟂}={\displaystyle \frac{∥{\mathit{r}}_{\mathrm{1}}\times {\mathit{r}}_{\mathrm{0}}∥}{∥{\mathit{r}}_{\mathrm{0}}∥}},\end{array}$$

and x_proj is constructed by keeping the axial projection of x onto the filament and setting its radial distance to ϵ.

For semi-infinite wake filaments, ϵ represents a viscous core of the far wake and is used only to desingularize the Biot–Savart kernel near the centreline. The corresponding perpendicular distance is

$$\begin{array}{}\text{(7)}& d_{⟂,\mathrm{\infty }}=∥{\mathit{r}}_{\mathrm{1}}\times {\widehat{\mathit{u}}}_{\mathrm{\infty }}∥,\end{array}$$

and for $d_{⟂,\mathrm{\infty }}<\mathit{ϵ}$, the induced velocity is evaluated at the projected point x_proj without additional $d_{⟂,\mathrm{\infty }}/\mathit{ϵ}$ scaling:

$$\begin{array}{}\text{(8)}& {\mathit{u}}_{\mathrm{ind},\mathrm{core}}\left(\mathit{x}\right)={\mathit{u}}_{\mathrm{ind}}\left({\mathit{x}}_{\mathrm{proj}}\right),d_{⟂,\mathrm{\infty }}<\mathit{ϵ}.\end{array}$$

2.3 Aerodynamic influence coefficient matrix

To reduce computational cost during the circulation-matching iterations, the aerodynamic influence coefficients were precomputed once. Let $\mathbf{A}\in {\mathbb{R}}^{\mathrm{3}{n}_{\mathrm{p}}\times n_{\mathrm{p}}}$ denote the block aerodynamic influence coefficient (AIC) matrix that maps the spanwise circulation vector $\mathbf{\Gamma }\in {\mathbb{R}}^{n_{\mathrm{p}}}$ to the induced velocity at each control point. For panel $i=\mathrm{1},\mathrm{\dots },n_{\mathrm{p}}$, the induced velocity is given by

$$\begin{array}{}\text{(9)}& {\mathit{u}}_{\mathrm{ind},i}=\sum _{j=\mathrm{1}}^{n_{\mathrm{p}}}{\mathbf{A}}_{ij}{\mathrm{\Gamma }}_j,\end{array}$$

where ${\mathbf{A}}_{ij}\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{1}}$ denotes the (i,j) block of A and represents the induced velocity at control point i due to a unit circulation on panel j.

To ensure compatibility with non-linear airfoil polars, the near-field 2D self-induced contribution of panel i is subtracted at its own control point to avoid double-counting induction effects already embedded in the sectional polar data (Piszkin and Levinsky, 1976; Ranneberg, 2015; Damiani et al., 2019; Cayon et al., 2023). Define the per-unit-circulation 2D self-induced influence vector at control point i as

$$\begin{array}{}\text{(10)}& {\mathit{a}}_{\mathrm{ind},\mathrm{2}\mathrm{D},i}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{{\widehat{\mathit{y}}}_{\mathrm{p},i}\times \left({\mathit{x}}_i-{\mathit{x}}_{\mathrm{0},i}\right)}{\Vert {\mathit{x}}_i-{\mathit{x}}_{\mathrm{0},i}{\Vert }^{\mathrm{2}}}},\end{array}$$

where ${\widehat{\mathit{y}}}_{\mathrm{p},i}$ is the local spanwise unit vector, and ${\mathit{x}}_i-{\mathit{x}}_{\mathrm{0},i}$ is the displacement from the bound-vortex centre of panel i to its control point. For a given circulation Γ_i, the corresponding 2D self-induced velocity is

$$\begin{array}{}\text{(11)}& {\mathit{u}}_{\mathrm{ind},\mathrm{2}\mathrm{D},i}={\mathrm{\Gamma }}_i{\mathit{a}}_{\mathrm{ind},\mathrm{2}\mathrm{D},i}.\end{array}$$

The effective induced velocity used for circulation matching is then

$$\begin{array}{}\text{(12)}& {\mathit{u}}_{\mathrm{eff},i}={\mathit{u}}_{\mathrm{ind},i}-{\mathit{u}}_{\mathrm{ind},\mathrm{2}\mathrm{D},i}.\end{array}$$

This subtraction is applied only to the self-induced contribution of panel i and not to the cross-induced contributions from other panels.

2.4 Iterative solving for circulation

With the AIC matrix assembled, the spanwise circulation distribution was determined iteratively. At iteration k, the induced velocity vectors at all control points were computed collectively as

$$\begin{array}{}\text{(13)}& \mathrm{vec}\left({\mathbf{U}}_{\mathrm{ind}}^{\left(k\right)}\right)=\mathbf{A}{\mathbf{\Gamma }}^{\left(k\right)},\end{array}$$

where vec(⋅) denotes the stacking of the three velocity components over all control points, and ${\mathbf{U}}_{\mathrm{ind}}^{\left(k\right)}\in {\mathbb{R}}^{n_{\mathrm{p}}\times \mathrm{3}}$ is the corresponding reshaped matrix whose ith row contains ${\mathit{u}}_{\mathrm{ind},i}^{\left(k\right)}$.

The effective relative velocity at the control point of panel i was then obtained by adding the local freestream velocity u_∞,i to the induced velocity:

$$\begin{array}{}\text{(14)}& {\mathit{u}}_{\mathrm{rel},i}^{\left(k\right)}={\mathit{u}}_{\mathrm{\infty },i}+{\mathit{u}}_{\mathrm{ind},i}^{\left(k\right)}.\end{array}$$

For each panel i, the sectional lift coefficient $C_{\mathrm{l},i}^{\left(k\right)}$ was obtained from the 2D polar using the effective angle of attack ${\mathit{\alpha }}_{\mathrm{eff},i}^{\left(k\right)}$:

$$\begin{array}{}\text{(15)}& {\mathit{\alpha }}_{\mathrm{eff},i}^{\left(k\right)}=\arctan \left({\displaystyle \frac{{\mathit{u}}_{\mathrm{rel},i}^{\left(k\right)}\cdot {\widehat{\mathit{z}}}_{\mathrm{p},i}}{{\mathit{u}}_{\mathrm{rel},i}^{\left(k\right)}\cdot {\widehat{\mathit{x}}}_{\mathrm{p},i}}}\right),\end{array}$$

where ${\widehat{\mathit{x}}}_{\mathrm{p},i}$, ${\widehat{\mathit{y}}}_{\mathrm{p},i}$, and ${\widehat{\mathit{z}}}_{\mathrm{p},i}$ denote the local chordwise (LE-to-TE), spanwise, and panel-normal (upward) unit vectors, respectively.

The circulation at each panel was updated using the Kutta–Joukowski theorem, normalized by the local reference velocity magnitude:

$$\begin{array}{}\text{(16)}& {\mathrm{\Gamma }}_{\mathrm{new},i}^{\left(k\right)}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\Vert {\mathit{u}}_{\mathrm{rel},i}^{\left(k\right)}\times {\widehat{\mathit{y}}}_{\mathrm{p},i}{\Vert }^{\mathrm{2}}}{\Vert {\mathit{u}}_{\mathrm{\infty },i}\times {\widehat{\mathit{y}}}_{\mathrm{p},i}\Vert }}{c}_i{C}_{\mathrm{l},i}^{\left(k\right)},\end{array}$$

where c_i denotes the local chord length.

To ensure numerical stability, under-relaxation was applied globally to the circulation array:

$$\begin{array}{}\text{(17)}& {\mathbf{\Gamma }}^{(k+\mathrm{1})}=(\mathrm{1}-f_{\mathrm{relax}}){\mathbf{\Gamma }}^{\left(k\right)}+f_{\mathrm{relax}}{\mathbf{\Gamma }}_{\mathrm{new}}^{\left(k\right)}.\end{array}$$

The convergence metric δ_Γ was defined as the normalized maximum change in the global circulation distribution between successive iterations:

$$\begin{array}{}\text{(18)}& {\mathit{\delta }}_{\mathrm{\Gamma }}={\displaystyle \frac{max\left|{\mathbf{\Gamma }}^{(k+\mathrm{1})}-{\mathbf{\Gamma }}^{\left(k\right)}\right|}{max\left|{\mathbf{\Gamma }}^{(k+\mathrm{1})}\right|}}.\end{array}$$

2.5 Converting to global properties

Once convergence was reached, the local aerodynamic force on each panel was computed as

$$\begin{array}{}\text{(19)}& {\mathit{F}}_{\mathrm{i}}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{u}_{\mathrm{rel},i}^{\mathrm{2}}\left(C_{\mathrm{l},i}{\widehat{\mathit{n}}}_i+C_{\mathrm{d},i}{\widehat{\mathit{t}}}_i\right)c_i{w}_i,\end{array}$$

where ρ is the fluid density, w_i is the panel width, and C_l,i and C_d,i are the sectional lift and drag coefficients, respectively, evaluated at the three-quarter-chord control point. In other words, the local angle of attack evaluated at the three-quarter chord was used to obtain the aerodynamic coefficients from the 2D input data. The magnitude u_rel,i corresponds to the in-plane component of the relative velocity vector u_rel,i:

$$\begin{array}{}\text{(20)}& u_{\mathrm{rel},i}=∥{\mathit{u}}_{\mathrm{rel},i}\times {\widehat{\mathit{y}}}_{\mathrm{p},i}∥.\end{array}$$

The force direction was set by an in-plane effective-flow direction following Ranneberg (2015), Jonkman (2021), and Damiani et al. (2019). This differs from the approach suggested by Li et al. (2022) and implemented by Cayon et al. (2023); the choice is motivated by comparing simulations of two geometries with wind tunnel measurements, detailed in Appendix B.

For each panel $i=\mathrm{1},\mathrm{\dots },n_{\mathrm{p}}$, the in-plane drag-direction unit vector is defined as

$$\begin{array}{}\text{(21)}& {\widehat{\mathit{t}}}_i={\displaystyle \frac{\cos \left({\mathit{\alpha }}_{\mathrm{eff},i}\right){\widehat{\mathit{x}}}_{\mathrm{p},i}+\sin \left({\mathit{\alpha }}_{\mathrm{eff},i}\right){\widehat{\mathit{z}}}_{\mathrm{p},i}}{∥\cos \left({\mathit{\alpha }}_{\mathrm{eff},i}\right){\widehat{\mathit{x}}}_{\mathrm{p},i}+\sin \left({\mathit{\alpha }}_{\mathrm{eff},i}\right){\widehat{\mathit{z}}}_{\mathrm{p},i}∥}}.\end{array}$$

The lift was oriented in the normal direction ${\widehat{\mathit{n}}}_i$, which was constructed as

$$\begin{array}{}\text{(22)}& {\widehat{\mathit{n}}}_i={\displaystyle \frac{{\widehat{\mathit{t}}}_i\times {\widehat{\mathit{y}}}_{\mathrm{p},i}}{\Vert {\widehat{\mathit{t}}}_i\times {\widehat{\mathit{y}}}_{\mathrm{p},i}\Vert }}.\end{array}$$

The local aerodynamic moment on panel i was defined as

$$\begin{array}{}\text{(23)}& {\mathit{M}}_i=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{u}_{\mathrm{rel},i}^{\mathrm{2}}{C}_{\mathrm{m},i}{\widehat{\mathit{y}}}_{\mathrm{p},i}{w}_i{c}_i^{\mathrm{2}},\end{array}$$

where C_m,i is the sectional moment coefficient; this moment represents the pitching moment about the quarter-chord point of the panel, consistent with the definition of C_m,i.

The contributions of all panels were summed and projected onto both the global coordinate directions ${\widehat{\mathit{e}}}_x$, ${\widehat{\mathit{e}}}_y$ and ${\widehat{\mathit{e}}}_z$, and the body-fixed aerodynamic directions: drag ${\widehat{\mathit{e}}}_{\mathrm{D}}$, aligned with u_∞; lift ${\widehat{\mathit{e}}}_{\mathrm{L}}$, perpendicular to drag; and side ${\widehat{\mathit{e}}}_{\mathrm{S}}$, orthogonal to both drag and lift. With $U_{\mathrm{\infty }}=\Vert {\mathit{u}}_{\mathrm{\infty }}\Vert$, the resulting non-dimensional force and moment coefficients are

$$\begin{array}{ll}{\displaystyle }{C}_j& {\displaystyle }={\displaystyle \frac{\mathrm{2}}{\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{2}}S}}\sum _{i=\mathrm{1}}^{n_{\mathrm{p}}}{\mathit{F}}_i\cdot {\widehat{\mathit{e}}}_j,\\ \text{(24)}& {\displaystyle }& {\displaystyle }\text{for }{\widehat{\mathit{e}}}_j\in \mathit{\{}{\widehat{\mathit{e}}}_x,{\widehat{\mathit{e}}}_y,{\widehat{\mathit{e}}}_z,{\widehat{\mathit{e}}}_{\mathrm{D}},{\widehat{\mathit{e}}}_{\mathrm{L}},{\widehat{\mathit{e}}}_{\mathrm{S}}\mathit{\}},{\displaystyle }{C}_{\mathrm{M},j}& {\displaystyle }={\displaystyle \frac{\mathrm{2}}{\mathit{\rho }{U}_{\mathrm{\infty }}^{\mathrm{2}}S{c}_{\mathrm{ref}}}}\sum _{i=\mathrm{1}}^{n_{\mathrm{p}}}\left[{\mathit{M}}_i+\left({\mathit{r}}_{\mathrm{c}/\mathrm{4},i}-{\mathit{r}}_{\mathrm{0}}\right)\times {\mathit{F}}_i\right]\cdot {\widehat{\mathit{e}}}_j,\\ \text{(25)}& {\displaystyle }& {\displaystyle }\text{for }{\widehat{\mathit{e}}}_j\in \mathit{\{}{\widehat{\mathit{e}}}_x,{\widehat{\mathit{e}}}_y,{\widehat{\mathit{e}}}_z\mathit{\}},\end{array}$$

where S is the projected planform area, c_ref the reference chord length, n_p the number of panels, ${\mathit{r}}_{c/\mathrm{4},i}$ the quarter-chord location about which the moment is defined for each panel, r₀ the chosen global reference point, and ${\widehat{\mathit{e}}}_j$ the projection direction. The centre of pressure r_cp was determined by

$$\begin{array}{}\text{(26)}& {\mathit{r}}_{\mathrm{cp}}={\mathit{r}}_{\mathrm{0}}+{\displaystyle \frac{\mathit{F}\times {\mathit{M}}_{\mathrm{0}}}{\Vert \mathit{F}{\Vert }^{\mathrm{2}}}},\end{array}$$

where F is the combined aerodynamic force vector of the wing, M₀ is the combined moment vector about the reference point r₀, and $\Vert \mathit{F}\Vert$ denotes the Euclidean norm of F. The resulting vector r_cp lies along the line of action of the aerodynamic force and represents the unique point where this line intersects the wing surface. This intersection defines the effective location at which the resultant force can be considered to act, and it provides a critical indicator of aerodynamic stability, especially in relation to the trim state of the kite.

3 Airfoil aerodynamics

This section outlines the modelling framework used to generate aerodynamic airfoil data. It describes the profile parameterization, mesh generation, and simulation setup employed for the CFD analyses, followed by the convergence monitoring procedure.

3.1 Profile geometry parameterization

The parameterization presented here extends the earlier approach by Breukels (2011), who introduced two fundamental parameters: the circular LE tube diameter t and the camber height κ, shown in Fig. 4. Subsequently, informed by the experimental findings of Den Boer (1980), who observed a forward shift of the maximum camber position with increasing angle of attack, the chordwise location of maximum camber η was incorporated as an essential design parameter (Watchorn, 2023). Another step was to reduce the trailing-edge height from $t/\mathrm{2}$ to zero, in accordance with the conventional airfoil definition (Kappel, 2012). A TE reflex angle δ is introduced, as this is often used to enhance the static pitching moment C_m and reduce the risk of front stall.

The LEI kite profile is constructed by defining a circular LE tube, smoothly extending it towards the maximum camber location using a Bézier curve named the front spline, and subsequently connecting it to the TE via another Bézier curve called the rear spline (see Fig. 4). Tangent continuity is enforced at both the junction between the curves and at the TE, ensuring aerodynamic smoothness. In this way, the profile shape can be characterized through six design parameters, summarized in Table 1.

Figure 4 Overlay of the LEI kite profile parameterizations using Bézier curves.

Table 1 Overview of the profile shape parameters used in the Bézier-curve-based parameterization.

Each Bézier curve is defined by four control points [$P_{\mathrm{0}},P_{\mathrm{1}},P_{\mathrm{2}},P_{\mathrm{3}}$], with the parametric coordinate $t_{\mathrm{b}}\in [\mathrm{0},\mathrm{1}]$ controlling the interpolation:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\left[\begin{array}{c}x\\ y\end{array}\right]& ={\left(\mathrm{1}-t_{\mathrm{b}}\right)}^{\mathrm{3}}{P}_{\mathrm{0}}+\mathrm{3}{\left(\mathrm{1}-t_{\mathrm{b}}\right)}^{\mathrm{2}}{t}_{\mathrm{b}}{P}_{\mathrm{1}}\\ & +\mathrm{3}\left(\mathrm{1}-t_{\mathrm{b}}\right)t_{\mathrm{b}}^{\mathrm{2}}{P}_{\mathrm{2}}+t_{\mathrm{b}}^{\mathrm{3}}{P}_{\mathrm{3}}.\end{array}\end{array}$$

The front spline connects the seam to the maximum camber point. The seam, i.e. the stitched junction at which the canopy is attached to the LE tube, is located on the tube circumference at the seam angle. Because a tensioned membrane cannot sustain a bending moment, the canopy must depart the tube surface tangentially; tangent continuity is therefore enforced at the seam. The seam angle is determined by a dynamic search algorithm that increases the angle stepwise until the apex of the resulting front spline coincides with the prescribed maximum camber point, thereby selecting the smallest seam angle for which a tangent departure is compatible with the prescribed camber. This generalizes earlier fixed seam-angle definitions (Watchorn, 2023), which can produce invalid geometries for weakly cambered profiles. On the physical kite, the seam additionally acts as the roughness element that trips the boundary layer, supporting the fully turbulent modelling assumption discussed in Appendix C.

The intermediate control points of the front spline are positioned along the tangent lines emanating from the seam and the maximum camber point, with the LE curvature parameter ϕ specifying their relative distance along these tangents. In this way, ϕ shapes the curvature distribution: larger values concentrate curvature near the centre of the curve, whereas smaller values shift it towards the endpoints. As surface curvature sets the local pressure gradient, ϕ effectively redistributes the suction along the front spline. Within the physically consistent range, however, this redistribution only slightly alters the pressure distribution and has a negligible effect on the integral coefficients (Masure, 2025), as the dominant flow features of LEI airfoils (the suction peak on the LE tube and the pressure-side recirculation) are governed primarily by t, κ, and η. For this reason, ϕ is held fixed during dataset generation (Sect. 4).

The rear spline connects the maximum camber location to the TE. The intermediate control points are governed by the camber tension parameter λ, which adjusts profile smoothness, and by the TE reflex angle δ. To enhance model robustness and reduce the number of free parameters, the third control point is positioned at 20 % of the distance from the second rear control point to the TE, measured from the second rear control point.

3.2 Meshing

To ensure a smooth computational mesh, the LEI kite profile geometry requires specific adaptations. In particular, the sharp corner at the junction between the lower canopy and the LE tube is smoothed using a fillet. In addition, a canopy thickness t_canopy is introduced to represent a realistic finite thickness, to distinguish top from bottom airflow, and to enable a rounded TE. The optimal fillet size and t_canopy were determined through a sensitivity study. Masure (2025) evaluated seven leading-edge fillet sizes on two reference profiles, spanning enclosed cross-sectional areas of 0.13–45.52 cm² (conventional) and 4.22–137.12 cm² (highly cambered) for c= 1 m. The smallest numerically robust configurations, corresponding to enclosed areas of 0.98 and 4.22 cm², informed a fixed fillet construction rule that scales with the tube diameter and seam angle, and is applied uniformly to all geometries in the dataset.

Figure 5 Fully structured mesh (575 × 201) for an LEI kite profile with shape parameters t= 0.08, η= 0.15, κ= 0.08, δ= 0°, λ= 0.2, and ϕ= 0.65. The top panel shows the airfoil, and the lower panels depict, from left to right, the LE tube, the O-grid, and the TE.

A structured O-grid mesh with approximately 1.15 × 10⁵ cells, shown in Fig. 5, was generated for each profile using a hyperbolic-Laplace solver in Pointwise (Pointwise, 2025), automated via a Python script.

An O-grid topology was selected primarily for its insensitivity to angle-of-attack variations, which reduces the effort required for case setup, as only one mesh needs to be generated per profile. As the chord is normalized to unit length, the 2D mesh is extruded by a single cell over one chord length in the spanwise direction, as required by OpenFOAM's 3D formulation. In the wall-normal direction, the 201 cell layers are generated by hyperbolic extrusion with a fixed growth ratio of 1.1 following Watchorn (2023), starting from a first-layer height of 7.8 × 10⁻⁶ c, chosen to ensure $y_{\max }^{+}<$ 1 over the entire airfoil at Re = 5 ×10⁶. This places the circular far-field boundary at approximately 1.4 × 10⁴ chord lengths from the profile, sufficiently far for the freestream boundary conditions not to influence the computed force coefficients (Masure, 2025). Relative to the next-finer level, the selected production resolution of 575 × 201 cells changed the lift, drag, and moment coefficients by less than 2 % across the operational angles of attack, which is taken as the remaining discretization error.

For all cases, smoothing parameters for hyperbolic extrusion in Pointwise were applied uniformly to ensure mesh quality and numerical stability. The explicit smoothing parameter was set to 8 for transverse smoothing, while the implicit parameter was set to 16.0. Kinsey–Barth smoothing, with a value of 5, was used in the wall-normal direction to address regions of strong concavity. The volume smoothing parameter was set to its default value of 0.5.

Node spacing was uniform near the LE and TE, while a hyperbolic tangent distribution was applied along the canopy surface to minimize abrupt grid transitions. Mesh independence was assessed by Masure (2025), who independently varied the tangential (275–775 nodes) and wall-normal (161–221 layers) resolution of the structured O-grid for representative conventional and highly cambered profiles.

The non-dimensional wall distance y⁺ is defined as

$$\begin{array}{}\text{(28)}& y^{+}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{\infty }}y{u}_{\mathit{\tau }}}{{\mathit{\mu }}_{\mathrm{\infty }}}},\end{array}$$

where ρ_∞ denotes freestream density, μ_∞ the dynamic viscosity, and u_τ the friction velocity.

3.3 Simulation setup

By assuming steady-state and incompressible flow, the steady RANS equations were solved using OpenFOAMv2006, with pressure-velocity coupling managed internally by the semi-implicit method for the pressure-linked equations (SIMPLE) algorithm (OpenFOAM, 2006). For an overview of the numerical schemes, the reader is referred to Masure (2025). The revised k-ω SST model by Menter (1994) was employed due to its robustness in capturing both near-wall and outer-region flow physics, supported by prior implementations in aerodynamic studies of LEI kites (Folkersma et al., 2019; Viré et al., 2020, 2022; Watchorn, 2023). All simulations were performed under the assumption of fully turbulent flow, owing to the measured turbulence intensities, the surface roughness present on both suction and pressure sides, and the high Reynolds numbers; a more elaborate discussion of this modelling choice is provided in Appendix C.

The inlet velocity scalar u_∞, set to 1 m s⁻¹, was projected in the direction of the angle of attack α to define the inlet velocity vector

$$\begin{array}{}\text{(29)}& {\mathit{u}}_i=u_{\mathrm{\infty }}\left[\begin{array}{c}\cos \mathit{\alpha }\\ \sin \mathit{\alpha }\\ \mathrm{0}\end{array}\right].\end{array}$$

The turbulent kinetic energy k_i is defined as

$$\begin{array}{}\text{(30)}& k_i={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\left(u_{\mathrm{\infty }}I\right)}^{\mathrm{2}},\end{array}$$

where I represents the freestream turbulence intensity. The specific dissipation rate ω_i at the inlet was defined as

$$\begin{array}{}\text{(31)}& {\mathit{\omega }}_i={\displaystyle \frac{k_i}{{\mathit{\nu }}_{\mathrm{\infty }}}}{\left({\displaystyle \frac{{\mathit{\nu }}_t}{{\mathit{\nu }}_{\mathrm{\infty }}}}\right)}^{-\mathrm{1}},\end{array}$$

where ${\mathit{\nu }}_t/{\mathit{\nu }}_{\mathrm{\infty }}$ represents the eddy viscosity ratio. Values of I= 2 % and ${\mathit{\nu }}_t/{\mathit{\nu }}_{\mathrm{\infty }}=$ 10 were selected based on the sensitivity studies of Folkersma et al. (2019) and Viré et al. (2020).

Boundary conditions for velocity u, kinematic pressure p_k, turbulent kinetic energy k, specific dissipation rate ω, and eddy viscosity ν_t were applied as summarized in Table 2. Far-field boundaries employed non-reflecting types; the airfoil surface was defined as a stationary, impermeable no-slip wall with appropriate wall functions for turbulence; and spanwise boundaries were set to empty.

Table 2 Boundary conditions and units for velocity, pressure, and turbulence variables.

3.4 Convergence monitoring

Convergence was monitored through both equation residuals and aerodynamic force coefficients. Specifically, residuals associated with momentum, pressure, and turbulence variables u_x, u_y, p, k, and ω were required to fall below 10⁻⁶ (Ferziger et al., 2020), while the force coefficients C_l, C_d, and C_m were expected to converge to within 10⁻⁴.

In OpenFOAM, each outer iteration reports two residuals per variable: the initial residual, evaluated on the newly assembled equation before the inner linear solve and thus measuring how well the previous iterate satisfies it; and the final residual, which only reflects the inner linear-solver settings. Outer-loop convergence is therefore monitored through the initial residual, on which the automatic stopping criterion acts. Because this criterion considers residuals only, a simulation could terminate before the force coefficients satisfy the 10⁻⁴ criterion. To prevent such early termination and ensure consistency across all cases, the stopping threshold was set slightly below the nominal residual criterion, at 8 × 10⁻⁷ for all flow variables. Each simulation was limited to a maximum of 5000 iterations. Cases not reaching the required thresholds within this limit were classified as non-converged; these typically occurred at high or negative angles of attack, where steady RANS solutions become unreliable owing to increased unsteadiness, e.g. flow separation.

For improved numerical stability, the velocity relaxation factor was set to 0.8, which dampened residual oscillations and improved convergence consistency. Representative histories for a typical LEI profile at α= 6° are shown in Fig. 6, where final flow residuals fell below 10⁻⁶ and force residuals just below 10⁻⁴.

Figure 6 Flow residuals at α= 6° for a profile with t= 0.0782, η= 0.178, κ= 0.0953, δ= 0.351, λ= 0, ϕ= 0.1.

4 Machine learning

This section presents the development and validation of a machine learning (ML) surrogate model trained on the aerodynamic database generated in Sect. 3. It describes the generation and filtering of the input data; the training and tuning of multiple regression algorithms; and the evaluation of their predictive accuracy for lift, drag, and moment coefficients.

4.1 Defining input profiles

Training data for the ML model were generated by varying the profile parameters, as shown in Table 3, over a range that encompasses physically viable profile designs. As varying ϕ within the marginal physically consistent range was found to have minimal aerodynamic influence, it was fixed at ϕ= 0.65, allowing greater resolution of the remaining parameters. The angle-of-attack range and sampling density were chosen to optimize accuracy in the lift polar and to maintain computational efficiency, with fewer points in the linear region (2° $\le \mathit{\alpha }\le$ 10°) and negative α, and a finer step of 2° near stall. With ϕ held fixed, the sampled input space is six-dimensional (five varied-shape parameters and the angle of attack), analysed separately for three Reynolds numbers: Re = 1 × 10⁶, representative of small kites or static flight; Re = 5 × 10⁶, corresponding to typical in-flight operation; and Re = 2 × 10⁷, characteristic of larger kites and/or higher inflow velocities (Cayon et al., 2025; Masure, 2025).

Sampling Re at discrete levels is consistent with the fully turbulent modelling assumption (Appendix C). With transition fixed by the seam roughness elements and the elevated freestream turbulence encountered in operation, the regime changes associated with free transition are excluded. The force coefficients then vary weakly, smoothly, and monotonically with Re between the sampled levels, in line with turbulent skin-friction scaling, as quantified in Appendix D.

This is corroborated numerically by the 3D RANS results of Viré et al. (2022), which showed a negligible effect on the integral forces over a factor-two change in Re. It was also corroborated experimentally by the wind tunnel measurements presented in the companion paper (Poland et al., 2026b), in which the force coefficients varied only gradually and without abrupt regime changes across Reynolds numbers from 1.3 × 10⁵ to 6.1 × 10⁵. Only the lowest measured Reynolds number indicated an earlier stall onset, consistent with transitional effects at Reynolds numbers well below the fully turbulent regime considered here.

Table 3 Design parameter ranges of the generated data for Re = 1 × 10⁶, 5 × 10⁶, and 2 × 10⁷.

In preprocessing, the profiles were filtered out if the camber tension control point λ was positioned too far aft, i.e. $\mathit{\eta }+\mathit{\lambda }\ge$ 0.8, or if the camber height κ was excessively large, i.e. κ> 2.5t. After filtering, 3464 configurations remained per angle of attack, resulting in 45 032 simulations per Reynolds number.

4.2 Input data generation

Case generation was fully automated: for each profile-parameter combination in the design space, the scripts executed geometry construction, mesh generation, and case setup for the required angles of attack and Reynolds numbers without manual intervention. The resulting RANS CFD simulations were performed on the HPC-12 cluster at Delft University of Technology, with each case decomposed over four cores, requiring a total of 1.4 × 10⁵ core hours, corresponding to an average of approximately 1 core hour per simulation. A limited number of cases did not converge within 5000 iterations, predominantly at high angles of attack beyond stall, and were subsequently removed from the data.

Additionally, a z-score filtering approach (Dekking et al., 2005),

$$\begin{array}{}\text{(32)}& z={\displaystyle \frac{x-\mathit{\mu }}{\mathit{\sigma }}},\end{array}$$

where μ represents the mean and σ the standard deviation, was used to identify and remove statistical outliers. Specifically, for each α and force coefficients, any data point with an absolute z score greater than 3 standard deviations from the mean was classified as an outlier and removed, illustrated in Fig. 7.

Figure 7 Aerodynamic coefficients versus angle of attack for the Re = 5 × 10⁶ dataset, illustrating filtering of cambered profiles, camber-limited profiles, and removed data.

Profiles with a maximum camber smaller than tube radius, i.e. $\mathit{\kappa }<\frac{t}{\mathrm{2}}$, are defined as camber-limited profiles and represent 0.23 % of the training dataset. For these cases, the number of parameters affecting the front and rear splines was reduced to t and δ, as the removed parameters had little to no effect. The camber-limited profiles exhibit distinct aerodynamic characteristics, producing lower lift and stalling earlier at approximately α= 10°. Consequently, they exhibit relatively high drag coefficients post-stall and increased pitching moments, indicating a stronger nose-up tendency in pre-stall flight conditions.

Building a predictive model within this design space requires the definition of both training and test datasets. A separate test dataset was constructed at different angles of attack from those used in the training set to verify the predictions at angles on which the model was not trained, shown in Table 4. After removing non-converged cases and statistical outliers, as summarized in Table 4, the final dataset comprises 42 183 training samples and 5016 test samples, corresponding to a test–train split of 10.6 %.

The test dataset exhibited a higher overall convergence rate, which was attributed to its narrower angle-of-attack bounds, making convergence less prone to failure near the limits, i.e. at negative α or under stall conditions.

Table 4 Data filtering steps and resulting sizes for the Re = 5 × 10⁶ dataset.

4.3 Model setup

A separate ML model was trained for each Re to predict C_l, C_d, and C_m from any given set of the design parameters and the angle of attack, resulting in a model with a 7D input space (six shape parameters and α) and a 3D output space. Because α enters as a continuous input feature, the model can be evaluated at any requested angle of attack, not only at the values contained in the training data.

By contrast, Re was deliberately not included as a model input: with only three sampled levels, a tree-based regressor would represent the Re-dependence as piecewise-constant steps between the training values. Instead, a separate model was trained per Re, and for intermediate Reynolds numbers the sectional polars are obtained by interpolating between the bracketing models in log Re, which yields a smooth and monotonic Re-dependence. The resulting surrogate model is applicable within the trained range Re $\in [\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}},\mathrm{2}\times {\mathrm{10}}^{\mathrm{7}}]$ under the fully turbulent assumption (Appendix C); transitional flow at substantially lower chord-based Reynolds numbers falls outside its scope.

The quantities C_l, C_d, and C_m were selected for the output because they represent key aerodynamic metrics and allow for the determination of the chordwise centre of pressure x_cp, which is a crucial stability indicator. The centre of pressure is determined through

$$\begin{array}{}\text{(33)}& x_{\mathrm{cp}}=x_{\mathrm{ref}}-{\displaystyle \frac{C_{\mathrm{m}}c}{\sqrt{C_{\mathrm{l}}^{\mathrm{2}}+C_{\mathrm{d}}^{\mathrm{2}}}\sin \left({\tan }^{-\mathrm{1}}\left(\frac{C_{\mathrm{l}}}{C_{\mathrm{d}}}\right)+\mathit{\alpha }\right)}},\end{array}$$

which yielded values within 1 % of those obtained by pressure integration (Masure, 2025).

Model development began with standardizing the filtered training datasets, after which several regression algorithms were evaluated using their default settings (Pedregosa et al., 2011). The selected models were chosen for their strong baseline performance reported by Pandala (2021), and linear regression was included for reference. The full ML workflow is described in Appendix D and illustrated in Fig. D1.

After an initial model run with default settings, hyperparameters were systematically fine-tuned using grid-search cross-validation. In this context, validation refers to the internal process of developing the regression model across different hyperparameter combinations using subsets of the training data before evaluating the final model's performance on an independent test dataset. The coefficient of determination R² (Dekking et al., 2005), defined as

$$\begin{array}{}\text{(34)}& R^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}},\end{array}$$

was selected as the scoring metric, where y_i denotes the true value for each data point, ${\widehat{y}}_i$ the corresponding model prediction, $\stackrel{\mathrm{‾}}{y}$ the mean value of the observed data, and n the number of points in the dataset. The hyperparameter search utilized a k-fold cross-validation approach, dividing the training set into k= 10 subsets or folds, with 90 % of the data used for training and 10 % for validation in each iteration. This was feasible because the dataset was sufficiently large to provide ample training data for improving model accuracy, measured here by R², while still reserving a large enough test set to ensure reliability, as quantified by a low standard deviation.

Tree-based regressors yield piecewise-constant responses, so although the model returns a value for any requested α, the raw prediction varies in fine steps rather than continuously. Where a continuously differentiable polar is required, e.g. for derivative-based analyses or within iterative solvers, the model output is sampled at the training angles of attack. These samples are then fitted with a piecewise cubic Hermite interpolating polynomial (PCHIP) spline, whose shape-preserving, monotone formulation guarantees C¹ continuity without introducing spurious oscillations. For a representative profile, raw predictions queried at arbitrary α deviate from the PCHIP-fitted polar by less than 0.03 in C_l (RMS 0.006), well below the regression error itself, confirming that the spline is a smoothness refinement rather than a required postprocessing step.

4.4 Model selection

After tuning was finalized, the performance of each model was assessed on the test dataset using three metrics: R², the percentage root mean square error

$$\begin{array}{}\text{(35)}& {\mathrm{RMSE}}_{\mathit{\%}}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n}{\sum }_{i=\mathrm{1}}^n{\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}}{\langle \left|y_{\mathrm{test}}\right|\rangle }}\times \mathrm{100},\end{array}$$

and the percentage standard deviation

$$\begin{array}{}\text{(36)}& {\mathit{\sigma }}_{\mathit{\%}}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{n-\mathrm{1}}{\sum }_{i=\mathrm{1}}^n{\left(r_i-\overline{r}\right)}^{\mathrm{2}}}}{\langle \left|y_{\mathrm{test}}\right|\rangle }}\times \mathrm{100},\end{array}$$

where n is the number of test samples, y_i and ${\widehat{y}}_i$ are the true and predicted values, $\langle \left|y_{\mathrm{test}}\right|\rangle$ is the mean absolute value of the test targets, and r_i and $\overline{r}$ are the residuals and their mean, respectively.

Table 5 shows the averaged results per model, from which extremely randomized trees (ERTs) were selected as they outperform the other algorithms, achieving the highest R², and lowest RMSE_% and σ_%.

Table 5 ML model verification using R², RMSE %, and σ_% for the Re = 5 × 10⁶ model.

4.5 Model evaluation

As shown in Table 6, the ERT exhibited high predictive accuracy with R²> 0.98 across all Re. The standard deviation σ_% decreased with increasing Re, indicating greater prediction consistency and improved alignment with CFD predictions at higher Re.

Table 6 Regression model performance metrics per Re.

To further evaluate the model's performance, three distinct and representative profiles were selected from the test dataset, detailed in Table 7. The third profile is notable as it represents a camber-limited configuration. The prediction metrics for these profiles are summarized in Table 7 and illustrated in Fig. 8 for C_l, C_d, and C_m.

In the plots, the ML model is evaluated at the angles of attack present in the training data, and the PCHIP spline introduced in Sect. 4.3 provides a C¹-continuous polar through these predictions, from which the coefficients are interpolated at the test angles of attack for comparison with the true test data.

Overall, profile 2, with t= 0.086 close to the LE tube thickness of the mid-span profile of the V3 kite, exhibited the highest R² values and the lowest RMSE_% values, indicating the greatest predictive certainty. No consistent trend was observed in a specific force coefficient; for example, R² and RMSE_% were highest either for C_d or for C_m. In general, a lower R² corresponded to a higher RMSE, with the RMSE_% showing larger absolute changes when mispredictions occurred. For instance, in the drag polar of profile 2, two clear mispredictions increased the RMSE_% by a factor of 2 compared with C_l, where only a single misprediction was present. The drag polars of profiles 1 and 3 revealed a similar pattern: in profile 1, one large misprediction resulted in an RMSE_% approximately four times greater.

Table 7 Prediction performance metrics for profiles 1–3.

Figure 8 Verification of the ML model using lift, drag, and moment polar for the Re = 5 × 10⁶ dataset.

5 Results

This section reports the principal findings of the study, beginning with flow field observations and subsequently presenting a quantitative comparison of LEI airfoil polars. In addition, the VSM was examined with respect to its convergence characteristics and sensitivity to airfoil-level inputs; rigid-body stability derivatives were derived; and a concise synthesis of best practices and limitations was provided.

5.1 Flow field observations

The flow topology at α= 2° and α= 8° is shown in Fig. 9, indicating the pressure-coloured streamlines around the V3 kite mid-span airfoil and the corresponding surface distributions of C_p and C_f. From the pressure-side streamlines, it is evident that a recirculation region is present, splitting into two zones. The first is a corner eddy close to the tube; the second, larger zone is more sensitive to inflow angle variations, decreasing in size as α increases. The rotation direction of the pressure-side recirculation zone can be inferred from the friction coefficient; a negative friction coefficient indicates counterclockwise rotation.

At α= 8°, the C_p indicates a comparatively larger negative pressure peak at the suction side and shows a distinct upward pressure peak at the centre of the recirculation zone. Where the recirculation zone stretches to the airfoil's surface, the friction coefficient C_f displays a peak in the opposite (downwards) direction. Increasing α from 2 to 8° leads to a more concentrated and more negative peak. Furthermore, double values are plotted for both C_p and C_f in the region where, along a vertical line, one encounters three pressure-side surfaces, a distinct feature of an LEI airfoil.

Figure 9 Flow fields at α= 2° and α=8° for the V3 kite at Re = 1 × 10⁶. The figures present the surface pressure coefficient C_p and skin-friction coefficient C_f as functions of the non-dimensional chordwise coordinate $x/c$.

5.2 Simulating airfoil aerodynamics

To assess the feasibility of inviscid methods for generating aerodynamic polars, NeuralFoil, a physics-informed surrogate model trained on tens of millions of XFOIL simulations and analytical solutions, was utilized (Drela, 2013; Sharpe, 2023, 2024).

Although XFOIL, on which NeuralFoil was trained, is not a purely inviscid method but a viscous–inviscid interaction solver, its integral boundary layer formulation assumes a thin shear layer coupled to an outer potential flow. Once separation becomes extensive and the flow no longer resembles a thin boundary layer attached to the surface, the mathematical foundation of the model is violated, and the predicted lift, drag, and pressure distributions can no longer be considered physically reliable (Drela, 2013). As LEI airfoils exhibit persistent pressure-side separation and recirculation that contravene these assumptions, the inclusion of NeuralFoil permits a structured assessment of the limitations of such viscous–inviscid approaches when applied to geometries dominated by extensive separated flow.

In this study, the largest model size, xxxlarge, was used; it achieves maximum errors of less than 0.3 % in the lift, drag, and moment coefficients relative to XFOIL (Sharpe, 2024). The boundary layer transition was set to $x/c=$ 0 to mimic the fully turbulent assumption used in the deployed 2D RANS CFD simulations.

To put the LEI airfoil aerodynamic trends into perspective, the lift and drag wind tunnel measurements of a flat plate were used (Wick, 1954). Because moment measurements were not available, the moment polar was computed following classical flat-plate airfoil theory:

$$\begin{array}{}\text{(37)}& C_{\mathrm{m}}=-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\sin \mathrm{2}\mathit{\alpha },\end{array}$$

where α is in radians.

Figure 10 plots the lift, drag, and moment coefficients as a function of α for the mid-span airfoil of the V3 kite (shown in Appendix A), comparing classical flat-plate theory, flat-plate wind tunnel measurements, Breukels' regression model, NeuralFoil, 2D RANS CFD results, and wind tunnel measurements (van Lith et al., 2026).

Figure 10 The figure shows the airfoil C_l, C_d, and C_m as functions of α, obtained from different sources evaluated at Re =10⁶. The NeuralFoil and CFD results use the ML profile-1 outline with minimal edge fillet, shown in Appendix A, to mimic the airfoil wind tunnel profile that had no edge fillet. The Breukels results used a larger edge fillet, corresponding to the CFD profile-1 outline in Appendix A. The numerical models, except for the flat-plate model, assume fully turbulent flow. This was mimicked in the wind tunnel measurements using zigzag strips, applied on the suction side at 3 % and on the pressure side at 45° of the tube, which were confirmed to be of sufficient height to induce transition (van Lith et al., 2026).

The predictions generated by Breukels’ regression model were considered invalid outside the range $\Vert \mathit{\alpha }\Vert \le$ 20°. At these higher angles of attack, the extrapolated polars exhibited unphysical behaviour, including lift and drag discontinuities near $\mathit{\alpha }=\pm$20° and negative drag coefficients for $\mathit{\alpha }<-$20°, which were attributed to insufficient smoothing between the empirical expressions used to extrapolate the model beyond the CFD-simulated range. Furthermore, assuming symmetry and antisymmetry about α= 0° was considered inappropriate due to the airfoil's inherent geometric asymmetry.

In terms of lift, the NeuralFoil predictions align well with wind tunnel measurements and 2D CFD predictions up to approximately α= 13°. Beyond this range, however, NeuralFoil overpredicts lift and delays stall. For drag, between $\mathit{\alpha }=-$5 and 10°, the NeuralFoil and Breukels' predictions exhibit similar trends to those of CFD simulations and wind tunnel measurements but are off in magnitude. At α> 20°, the C_d predictions of NeuralFoil and the CFD approach the flat-plate wind tunnel measurement in slope and magnitude.

Among all numerical approaches, 2D CFD simulations exhibit the closest agreement with wind tunnel measurements. Within the linear regime, from approximately $\mathit{\alpha }=-$1° to α= 10°, the average percentage errors are 4 % and 16 % for lift and drag, respectively. The largest discrepancies occur outside the linear regime, in the negative-lift regime and in the prediction of the stall angle, where 2D CFD indicates a slightly delayed stall onset.

The largest discrepancies among the model predictions arise in the moment coefficient, for which wind tunnel measurements were unavailable for both the LEI airfoil and the flat plate. Breukels' model yields the greatest deviation from all other approaches. From $\mathit{\alpha }=-$3 to 3° or so, the CFD results closely follow the flat-plate prediction for C_m. Beginning around α= 3°, the CFD-predicted C_m diverges from the flat-plate trend and instead flattens, maintaining an almost constant negative value from α= 5° to α= 12°. Beyond this point, a sharp decrease in C_m is observed, coinciding with the onset of stall. NeuralFoil is unable to capture this behaviour, and also deviates strongly with negative angles of attack.

5.3 VSM convergence

The VSM permits arbitrary spanwise panel distributions; for instance, the area between existing strut locations may be divided into equal parts, arranged with a uniform distribution or spaced according to a cosine-based scheme – commonly used in LL models to concentrate panels near the tips (Phillips and Snyder, 2000). However, for the V3 kite, the cosine distribution proved ineffective and failed to converge when the number of panels n_p exceeded 50, whereas a uniform panel distribution yielded the most robust and reliable convergence as n_p increased.

Within this paper, the simulations were considered converged when the metric δ_Γ< 1 × 10⁻⁶, and a relaxation factor of f_relax= 0.01 was employed to ensure numerical stability.

Figure 11 demonstrates the convergence of aerodynamic coefficients by showing the percentage error relative to the n_p= 150 case for α= 2, 8, and 16°. The n_p= 150 case was selected as the endpoint because the VSM often failed to converge for larger values of n_p. The results were computed using airfoil data obtained from CFD simulations fitted with PCHIP, and the used V3 kite airfoils are shown in Fig. A2 in Appendix A.

The convergence analysis demonstrated that convergence deteriorated with increasing α, with markedly less smooth trends at α= 16°, consistent with VSM's post-stall behaviour, discussed in Sect. 5.4 and 5.6. Because the solution converged with increasing n_p, the largest panel count – n_p= 150 – was adopted throughout the remainder of the paper.

Figure 11 Convergence analysis of the VSM for C_L, C_D, and C_M,y at α= 2, 8, and 16°. The error relative to the reference case with n_p= 150 is shown.

5.4 VSM sensitivity to airfoil aerodynamic input

To place the VSM results with different airfoil inputs in context, wind tunnel measurements and accompanying confidence interval (CI) (Poland et al., 2026b) and 3D RANS CFD simulations (Viré et al., 2022) are plotted in Fig. 12, all at Re = 5 × 10⁵ except for the 3D CFD simulations, which are at Re = 1 × 10⁶. The 3D CFD simulations demonstrated that both the inclusion of struts and the difference in Reynolds number between Re = 1.0 × 10⁶ (with struts) and Re = 5.0 × 10⁵ (without struts) had a negligible effect on the integral forces (Viré et al., 2022). The plotted 3D CFD results differ from those presented in Viré et al. (2022), as they were corrected for a geometric offset of Δα= 1.02°; further details are provided in the appendix of the companion paper (Poland et al., 2026b).

The VSM results, generated with 2D RANS CFD PCHIP-fitted data, differ slightly from the geometry used in the wind tunnel and 3D CFD, as struts are omitted and differences between the profiles' pressure sides emerge near the tips, an artefact of relying on the parameterization scheme and an automatic meshing script (see Fig. A2 in Appendix A). The VSM results based on NeuralFoil, meanwhile, use the geometries with struts. The pitching moment C_M,y was defined with respect to the tow point, defined as the location where the bridle lines converge beneath the wing, corresponding to the attachment point of the tether or kite control unit in AWE systems. The tow point of the V3 kite is located at x= 1.16 m, y= 0 m, and z= −11 m, relative to the mid-span LE point. The results are presented in Fig. 12.

Figure 12 The C_L, C_D, and C_M,y are plotted over an α, for β= 0°, for varying VSM simulations, wind tunnel measurements, and CFD simulations (Viré et al., 2020; Poland et al., 2026b).

The V3 kite operates at α≈ 1° during reel-in and α≈ 8° during reel-out (Cayon et al., 2025), placing the nominal operating conditions largely within the near-linear portion of the polar. Accordingly, comparisons are assessed over $\mathit{\alpha }\in [-\mathrm{1},\mathrm{10}]$° to bracket the operational range and its associated variability.

Within this band, all results shown in Fig. 12 exhibited a consistent linear increase in lift, except for the wind tunnel measurements, which showed a change in slope at around α= 9°. At negative α, the VSM results based on 2D CFD data aligned most closely with the measurements. At higher α, where LL models are known to be inaccurate, the VSM 2D CFD predictions surprisingly agreed more closely with the measurements than the 3D CFD data. In drag, the sensitivity of the VSM predictions to the airfoil-level inputs was also evident: the VSM 2D CFD case – the only one without struts – yielded the lowest drag. Considering both lift and drag, the VSM with 2D CFD input reproduces the measurements more closely than the 3D CFD predictions. Averaged over $\mathit{\alpha }\in [-\mathrm{1},\mathrm{10}]$°, the respective mean absolute errors are 9 % for lift and 13 % for drag.

The sensitivity of the VSM predictions is most evident in the pitching-moment coefficient C_M,y, for which no 3D CFD data were available. The VSM results deviated from the measured trend at negative α and became less negative at higher angles. Within the mentioned α-band, all simulations predicted a continually decreasing C_M,y, whereas the measurements showed a nearly constant and near-zero C_M,y.

Overall, it is concluded that the VSM predicted the most consistent trends when using 2D CFD-derived aerodynamic input, primarily because the associated lift and drag polars capture stall behaviour more accurately.

5.5 Rigid-body stability derivatives

Rigid-body stability derivatives were obtained by perturbing the apparent velocity angles and the body rates about a trimmed baseline and evaluating the resulting changes in aerodynamic force and moment coefficients. Angle derivatives with respect to the angle of attack α and sideslip β were reported per radian, whereas rate derivatives with respect to the body rates p, q, and r (roll, pitch, and yaw) were non-dimensionalized to hat-rates using

$$\begin{array}{}\text{(38)}& \widehat{p}={\displaystyle \frac{pb}{\mathrm{2}{U}_{\mathrm{\infty }}}},\widehat{q}={\displaystyle \frac{q{c}_{\mathrm{MAC}}}{\mathrm{2}{U}_{\mathrm{\infty }}}},\widehat{r}={\displaystyle \frac{rb}{\mathrm{2}{U}_{\mathrm{\infty }}}},\end{array}$$

where $c_{\mathrm{MAC}}=S/b=$ 2.346 m represents the mean aerodynamic chord with projected surface area S= 19.413 m² and projected span b= 8.274 m. Sensitivity $\partial C/\partial \mathit{\xi }$ was then approximated using a central finite-difference scheme:

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{\partial C}{\partial \mathit{\xi }}}\approx {\displaystyle \frac{C\left(\mathit{\xi }+{\mathrm{\Delta }}_{\mathit{\xi }}\right)-C\left(\mathit{\xi }-{\mathrm{\Delta }}_{\mathit{\xi }}\right)}{\mathrm{2}{\mathrm{\Delta }}_{\mathit{\xi }}}},\end{array}$$

with perturbation steps $\mathrm{\Delta }\mathit{\alpha }=\mathrm{\Delta }\mathit{\beta }=$ 0.005 rad and $\mathrm{\Delta }p=\mathrm{\Delta }q=\mathrm{\Delta }r=$ 0.1 rad s⁻¹. The reported derivatives were confirmed to lie close to the central-difference plateau: varying the angular step over 0.03–1.1° and the rate step over 0.005–0.1 rad s⁻¹ altered every derivative by less than 0.3 % and 1 %, respectively, with the exception of the two yaw-moment derivatives $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$ and $\partial C_{\mathrm{M},z}/\partial \widehat{p}$, which varied by up to 0.4 % and 3 %.

The perturbation steps were applied around the baseline trim angle, the angle at which $C_{\mathrm{M},y}=$ 0 and $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }<$ 0, further discussed in Sect. 6. The trim angle was found at α_trim= 7.469°, with $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }=-$2.708 rad⁻¹, indicating local pitch stability. Moments were computed about the tow-point, i.e. the junction where the bridle lines connect beneath the wing, as this represents the physical point about which the kite rotates.

The resulting derivatives, reported in the aircraft reference frame in Table 8, represent instantaneous rigid-body stability characteristics; they do not capture unsteady aerodynamic effects or the influence of deformation on the dynamic response. As deformation was not included, control effects are likewise absent, because the kite wing relies on deformation for actuation. Despite this limitation, the aerodynamically derived stability derivatives still provide insight into the kite’s intrinsic stability characteristics, particularly given the scarcity of comparable data in the literature and the reference-model status of the V3 kite.

Table 8 Computed rigid-body stability derivatives at the trim angle α_trim= 7.469°. All quantities are expressed in the aircraft reference frame, where the x axis is defined positive from the TE to the LE, the y axis is positive towards the right wing, and the z axis is positive downwards, forming a right-handed coordinate system consistent with standard aerodynamic force and moment sign conventions.

The signs of the static derivatives – $\partial C_{\mathrm{M},x}/\partial \mathit{\beta }<\mathrm{0}$, $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }<\mathrm{0}$, and $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }>\mathrm{0}$ – indicated static stability in roll, pitch, and yaw, respectively (Etkin and Reid, 1996). The negative diagonal rate derivatives – $\partial C_{\mathrm{M},x}/\partial \widehat{p}$, $\partial C_{\mathrm{M},y}/\partial \widehat{q}$ and $\partial C_{\mathrm{M},z}/\partial \widehat{r}$ – indicate aerodynamic damping in roll, pitch, and yaw, with yaw damping an order of magnitude weaker, consistent with the absence of a vertical stabilizing surface. Aerodynamic damping alone does not guarantee dynamic stability; resolving the dynamic modes requires a full modal analysis, including inertia and unsteady effects.

5.6 Best practices and limitations

VSM, like other LL methods, is sensitive to implementation details, and inconsistent force evaluation can lead to non-physical performance improvements when these models are embedded in optimization loops.

For this reason, the present work has prioritized internal consistency and documented limitations over incremental extensions whose theoretical implications are not fully validated.

Examples of excluded incremental extensions include accelerating convergence by initializing the solver with a previously converged circulation distribution. This strategy was omitted because the solution exhibited sensitivity to the initial condition, particularly in the vicinity of stall.

Extensions proposed for swept and sideslip configurations were explored but not adopted. The first concerned a sweep-consistent kinematic projection, in which the sectional kinematics were evaluated in the plane normal to the local spanwise direction when coupling the 3D induced velocity to 2D polars, rather than in the global projection plane used in the present VSM formulation; this modification produced only modest changes in the predicted loads. The second involved a modified near-wake treatment, in which the trailing vorticity was initialized locally perpendicular to the lifting line before aligning with the far-wake direction (Goates and Hunsaker, 2021), instead of being directly aligned with the free stream. This modification reduced the predicted lift of the V3 kite by nearly 40 % and was therefore not retained.

At high angles of attack, convergence issues may arise when the lift curve becomes non-bijective, such that multiple values of α correspond to the same C_l, particularly in post-stall regimes (Anderson et al., 1980; Chattot, 2004; Mukherjee and Gopalarathnam, 2006; Hunsaker, 2007; Gallay and Laurendeau, 2015; Dias, 2016; Simonet et al., 2024). Similar discontinuities were observed for the V3 kite in post-stall. To address potential spanwise discontinuities, a simple smoothing algorithm was implemented to damp sharp local gradients in the circulation field during iteration. However, no consistent results were obtained, as the smoothing effectiveness varied with both the number of panels and the geometry. It was therefore decided to omit the proposed algorithm. The model of Simonet et al. (2024) was also tested, but it likewise did not yield consistent results and was therefore excluded.

Discrepancies in drag, most pronounced when comparing numerical results to wind tunnel measurements (Poland et al., 2026b), may arise from several sources. One notable factor is the absence of a spanwise viscous drag contribution in the baseline model. A correction for this effect, proposed by Gaunaa et al. (2024), accounts for spanwise flow by decomposing the local velocity vector into spanwise and perpendicular components. The correction was omitted from the current VSM implementation due to its empirical nature and limited impact, as it increased drag by approximately 1 %.

To aid convergence, smooth sectional polars are recommended. Computational cost can be reduced by using fewer panels or by switching from the open-source Python version to the open-source Julia version, which is about 30 to 50 times faster, and runs in approximately 3 ms with 36 panels and in 25 ms with 150 panels on a 12th Gen Intel^® Core™ i7-1265U CPU. The 3D RANS CFD simulations required 19 394 s of CPU time (Viré et al., 2022), indicating a speed-up on the order of 10⁶, although the processor specifications for the 3D RANS CFD computations were not documented.

6 Aerodynamic design considerations

The aerodynamic design of kites presents distinct challenges relative to conventional aircraft, stemming fundamentally from their reliance on structural tension for flight integrity and from differences in the centre of rotation. Kites are uniquely susceptible to two contrasting stall phenomena: front stall, which generally occurs at low, negative angles of attack; and back stall, which manifests at high, positive angles of attack.

A front stall is often an irrecoverable failure mode. Once the kite initiates a nose-down motion, it is driven into negative angles of attack, where it produces minimal or negative lift. As the aerodynamic loading vanishes, line tension drops and the tensile structure collapses, reducing the kite to an uncontrollable mass of fabric.

In contrast, a back stall usually allows for recovery.

During the pitching-back motion, the kite typically retains part of its line tension and thus much of its shape. Control authority can often be restored by depowering, which effectively reduces the angle of attack and stops the stall. When the pitching-back motion is asymmetric, or when depowering is applied too late or insufficiently, the kite often rotates naturally in such a way that its heaviest component – the LE – moves downwards first. As the kite accelerates downwards, the magnitude of the apparent wind velocity increases, which reduces the effective angle of attack. Once the angle of attack re-enters the attached-flow regime, the separated shear layer can reattach, thereby restoring aerodynamic loading. The resulting renewed line tension re-establishes control authority before ground impact, often with sufficient time to steer the kite and keep it airborne.

In-flight, bridled soft-wing kites naturally align themselves until the net aerodynamic force and the tether attached to the tow point become parallel, at which point the moment arm vanishes and $C_{\mathrm{M},y}=\mathrm{0}$ (Cayon et al., 2026). This equilibrium defines the trim angle of attack α_t, which remains nearly constant because the tow point offers negligible rotational resistance. For a fixed bridle geometry, this equilibrium occurs at a single inflow angle, as corroborated experimentally for two kites, including the V3 kite (Cayon et al., 2025). The angle changes primarily during the transition phases, because powering and depowering alter the front–rear bridle length ratio, which shifts the effective tow-point location and thereby modifies the pitching-moment balance that sets α_t.

For AWE systems, the design objective is to maximize the system's $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ (Loyd, 1980) at the trim angle while maintaining sufficient stability margins to prevent front stall, which requires careful coordination between the wing's aerodynamic characteristics and the bridle configuration.

This section therefore presents parametric studies to investigate the effects of airfoil shape in Sect. 6.1, wing planform geometry in Sect. 6.2, and tow-point location in Sect. 6.3. Gravitational contributions are neglected in this analysis, which is considered a reasonable assumption given the low weight of soft kites relative to aerodynamic forces (Cayon et al., 2026, 2025). Although aero-structural coupling and unsteady effects are not included, the results provide valuable insight into aerodynamic trends and design trade-offs, and should be interpreted in the context of these limitations.

6.1 Effect of airfoil shape variation

To quantify the sensitivity of kite aerodynamic loads to airfoil geometry, the ML model was deployed to systematically vary the geometric parameters κ, η, t, δ, and λ. The resulting airfoil polars were subsequently used as input for the VSM to evaluate the kite aerodynamic loads. All the airfoil slices, shown in Appendix A, were uniformly varied by the same percentage relative to the V3 kite reference, while ensuring that all parameter values, shown in Table 3, remained within the bounds used for ML model training. As the V3 kite had a reflex angle of δ= 0°, percentage-based variation was not feasible; therefore, each profile was modified instead by an equal absolute value. The resulting variations are plotted in Fig. 13, with the exception of λ, which, within the parameter bounds used for training the ML model, had the least visible profile change.

Figure 13 Effect of changing the profile shape on the aerodynamic performance and pitching moment of the airfoil and wing.

The parameter study revealed that the effect on aerodynamic forces is more pronounced than the effect on the pitching moment. It furthermore showed that trends in 2D and 3D aerodynamic performance are not always aligned. For example, while the 2D ratio $C_{\mathrm{l}}^{\mathrm{3}}/C_{\mathrm{d}}^{\mathrm{2}}$ exhibits a clear peak near α= 6°, the $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ gradually increases beyond this angle. Likewise, variations in t affected C_m differently depending on α; the 2D case exhibited a near constant trend at a 10° angle of attack, whereas the corresponding impact on the 3D C_M showed a more uniform, decreasing trend across the investigated range.

For AWE, the aerodynamic design objective is to maximize the system ratio $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$, evaluated at the corresponding trim angle. Compared to the V3 kite, the results indicate that an optimal design would have smaller values of κ, larger values of η, smaller values of t, and a small δ, although for the latter, the opposite holds for α> 8°.

6.2 Effect of planform variations

The effect of varying the wing planform geometry was analysed by simulating three distinct planforms. A cubic Bézier curve was fitted to the quarter-chord line of the V3 kite wing using four control points (P₁, P₂, P₃, and P₄) that preserved the endpoint positions at tip and root, as illustrated in Fig. 14. Two shape parameters, namely the distances from P₁ to P₂ and from P₃ to P₄, were employed to capture arc fullness and tip curvature. These distances were initially obtained by non-linear least-squares minimization of the total point-to-curve distance between measured quarter-chord points. The anhedral angle ϕ_an, also shown in Fig. 14, was defined as the angle between the line from P₁ to P₂ and the line from P₁ to P₄.

Figure 14 Effect of changing the V3 kite's aspect ratio and ϕ_an, from its base properties, while maintaining the same flat surface area.

When reconstructing the V3 using cubic Bézier curves, an aspect ratio of AR= 5 was found, with AR defined as

$$\begin{array}{}\text{(40)}& \mathrm{AR}={\displaystyle \frac{b^{\mathrm{2}}}{S}},\end{array}$$

where b= 11 m denotes the flat wingspan, not identical to the projected wingspan and measured as the quarter-chord arc length in the (y,z) plane when laid out flat.

A second planform with the same non-dimensional arc shape, but with a 25 % larger aspect ratio, was generated by keeping the shape parameters fixed, and scaling span and chord. A third planform was produced by increasing the anhedral angle by 25 % to ϕ_an= 41.73°. The flattened surface areas of the reference model and enlarged AR planform were nearly identical, where a reduction of less than 0.5 % was observed for the planform with enlarged ϕ_an.

The aerodynamic effects are presented for both α and β sweeps in Fig. 15. The trim angles, computed about the tow point, were α_t= 6.89° for the baseline configuration, α_t= 6.23° for the AR-scaled configuration, and α_t= 7.20° for the ϕ_an-scaled configuration. As the maximum difference of 0.97° was not considered substantial in view of the applied shape changes, and to isolate the effect of the planform change, the β sweeps of all configurations were evaluated at the baseline trim angle of 6.89°.

Figure 15 Effect of changing the V3 kite, from its baseline, is shown by plotting the force and moment coefficients over an α and β sweep, the latter evaluated at the baseline trim angle α_t= 6.89°.

Increasing the aspect ratio AR, in line with classical aircraft aerodynamic theory, resulted in higher $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ values. This was accompanied by a more positive yaw-moment derivative $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$. From this, one could conclude that a kite with a larger AR would generate a greater restoring yaw moment in response to an offset in β, thereby enhancing yaw static stability.

Increasing the anhedral angle ϕ_an reduced $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$.

In addition, the lateral–directional coefficients exhibited stronger coupling: the side-force coefficient C_S increased in magnitude for a given β, indicative of enhanced turning capability.

The results demonstrated a fundamental trade-off in kite design. Increasing AR yielded a higher lift-to-drag ratio, whereas increasing the ϕ_an lowered $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ but enhanced the turning ability.

6.3 Effect of tow-point location

The sensitivity of the kite to variations in tow-point position was analysed by displacing its location in the x and z directions, as illustrated in Fig. 16, which also includes a V3 kite mid-span slice showing the LEI airfoil at the design distance, spanning precisely the figure x limits. For each tow-point configuration, the trim angle α_t and its corresponding $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ ratio were determined using the VSM. Since the VSM tends to lose accuracy in post-stall conditions, the analysis was restricted to the range $\mathit{\alpha }\in [-\mathrm{3},\mathrm{13}\mathit{°}]$.

The figures include a grey-shaded area indicating where $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }>\mathrm{0}$, i.e. where the kite is not statically stable in pitch and no physically meaningful trim angle can be found. Additionally, the centre of gravity (CG) and the range of centre of pressure (CP) positions corresponding to the obtained trim angles were calculated and plotted.

Figure 16 Effect of tow-point displacement on the V3 kite aerodynamic characteristics. The left panel shows the mid-span airfoil, the CP range, and the CG. All plots show isolines of varying parameters. The left plot shows the trim angle α_t as a function of tow-point location in the x and z directions. The other plots show the isolines for $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$, $\partial C_{\mathrm{M},x}/\partial \mathit{\beta }$, $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }$, and $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$, with statically unstable regions marked by grey shading.

The results indicate that the trim angle α_t is more sensitive to variations in x than in z. The same trend is observed for $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$, $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }$, and $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$. The largest trim angles α_t, corresponding to the largest values of $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$, occur when the tow point is positioned aft of the CP. However, positioning the tow point excessively close to the TE should be avoided, as it may lead to static pitch instability. More restrictive is the yawing-moment derivative $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$, which shows that parts of the region yielding the highest $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$ coincide with $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }<\mathrm{0}$, implying static instability in yaw.

The $\partial C_{\mathrm{M},x}/\partial \mathit{\beta }$ is the only parameter that is more sensitive to variations in z than in x, and it indicates static roll stability at the design tow-point location. The results show that raising the tow point can eliminate the restoring roll moment. This implies the existence of a boundary in tow-point height at which the effective roll moment arm about the wing changes sign, rendering the kite statically unstable in roll. As the CG lies above this boundary, placing the tow point at the CG location would lead to static roll instability. This behaviour is consistent with expectations for an anhedral wing in conventional aircraft (Etkin and Reid, 1996), and it underscores a key distinction: kites that rotate about a sufficiently low tow point require anhedral for roll stability, whereas aircraft that rotate about the CG require dihedral.

The tow-point variation analysis shows that to avoid static pitch and yaw instability, the tow point must be placed sufficiently far forward, which becomes more restrictive at higher z values. The stability gains of placing the tow point forward come at the cost of a reduced trim angle and, thereby, a lower aerodynamic performance.

7 Conclusions

A computationally efficient aerodynamic framework for soft-wing kite design was developed by coupling a vortex step method (VSM) with CFD-informed airfoil polars obtained from more than 10⁵ RANS CFD simulations. The resulting novel RANS-trained machine learning (ML) surrogate model predicts lift, drag, and moment coefficients for parameterized LEI airfoils within the trained design space; three models were trained at Re = 1 × 10⁶, 5 × 10⁶, and 2 × 10⁷, achieving R²> 0.98 on held-out test data, with intermediate Reynolds numbers covered by interpolation in log Re.

Within the operational, near-linear incidence band $\mathit{\alpha }\in [-\mathrm{1},\mathrm{10}]$°, selected to bracket the reported reel-in and reel-out angles of attack of the V3 kite, the underlying 2D CFD data agreed with wind tunnel measurements with mean errors of 4 % in lift and 16 % in drag. When these polars were used as VSM inputs, 3D lift and drag predictions matched wind tunnel results within this band on average by 9 % and 13 %, respectively, reproducing measured lift and drag trends more consistently than the available 3D RANS CFD results. The latter required 19 394 s of CPU time, whereas the present framework evaluated the n_p= 150 case in approximately 25 ms on the stated CPU, corresponding to an estimated speed-up of 𝒪(10⁶) for the compared configurations, subject to differences in hardware, parallelization, and included preprocessing. This computational efficiency facilitates optimization, aero-structural coupling, and real-time control applications.

Rigid-body analysis yielded a novel, numerically derived set of static and rate (damping) stability derivatives for the TU Delft V3 reference kite. The configuration was found to be statically stable in roll, pitch, and yaw for the examined trim condition, and the computed rate derivatives indicate aerodynamic damping within the present quasi-steady aerodynamic model. A dedicated modal analysis, including unsteady aerodynamics, is required to characterize the associated dynamic modes.

Parametric analysis revealed that configurations maximizing aerodynamic efficiency, expressed as $C_{\mathrm{L}}^{\mathrm{3}}/C_{\mathrm{D}}^{\mathrm{2}}$, do not coincide with those improving static stability margins (e.g. via $\partial C_{\mathrm{M},y}/\partial \mathit{\alpha }$ and $\partial C_{\mathrm{M},z}/\partial \mathit{\beta }$), indicating a multi-objective design problem. Within the present quasi-steady, rigid-geometry analysis, the roll-stability mechanism reverses: anhedral is stabilizing for kites rotating around the tow point, whereas dihedral is stabilizing for aircraft rotating about the centre of gravity.

The framework provides a validated and computationally efficient basis for conceptual and preliminary kite design. Future work should integrate aero-structural coupling to resolve wing deformation and its feedback on aerodynamic loading, for example, by updating spanwise sectional polars in response to the local deformed airfoil geometry and operating point. In addition, optimization studies would benefit from incorporating variations in bridle geometry and tow-point positioning to quantify their coupled influence on aerodynamic performance and stability.

Appendix A: Profile distribution

Geometric differences exist between the various simulation approaches. To illustrate these differences, Fig. A1 shows a spanwise slice taken at 25 % of the chord length of the V3 kite at α= 6°, plotted for half of the span.

Figure A1 Slice distribution along a representative spanwise section of the V3 kite at α= 6°, evaluated at 25 % chord. The figure is selected for illustrative purposes and does not correspond to the exact physical position.

Figure A2 Profiles used for simulations and their differences. A 2D CFD was also run on the dashed blue airfoil of profile 1 to better match the profile used in the 2D wind tunnel study, which had no edge fillet (van Lith et al., 2026).

As shown in Fig. A2, a total of 19 chordwise slices are defined, ranging from “1” at the mid-span to “19” at the tip. The three simulation techniques differ in the airfoil geometry across these profiles, illustrated in Fig. A2.

The 3D RANS CFD analysis conducted by Viré et al. (2020) employed a smoothed geometry, which incorporated substantial edge fillets downstream of the LE. The 2D RANS CFD simulations, undertaken to facilitate direct comparison with the 3D RANS CFD results, similarly included large edge fillets but did not account for the struts. The ML model, by comparison, was trained on 2D CFD datasets generated from geometries designed to approximate the physical configuration as closely as possible. Consequently, it cannot predict geometries with such large edge fillets; hence, only fillets of significantly reduced dimension are present relative to the other simulated geometries.

As the tip region was approached, the 2D slice intersected the LE tube at an increasingly shallow angle; that is, less perpendicular and more parallel to the tube axis. This progression increased the disparity between the respective geometries, most apparent for profiles 18 and 19 in Fig. A2. As the trained ML model focuses on the LEI airfoil shape, it does not capture these changes.

Appendix B: Evaluating force direction

Figure B1 Effect of varying the force-direction evaluation point, shown using VSM simulations performed at the Reynolds numbers of their respective wind tunnel (WT) measurements, which are included for comparison. Plotted are C_L, C_D, and $C_{\mathrm{L}}/C_{\mathrm{D}}$ against α for the TU Delft V3 LEI kite and the ram-air kite (Poland et al., 2026b; Belloc, 2015).

Li et al. (2022) revisited the chordwise placement of the control point in LL-type formulations by embedding unsteady 2D thin airfoil theory into a generalized 3D framework. Building on the classical thin airfoil solution, they showed that the quasi-steady lift magnitude is formally linked to the effective angle of attack evaluated at the three-quarter-chord point, consistent with the Kutta condition and with the location where the no-through-flow boundary condition is enforced in Weissinger-type discretizations. In contrast, the direction of the aerodynamic force is governed by the local flow angle at the quarter-chord point, where the bound vortex resides.

Li et al. (2022) validated their correction against blade-resolved RANS CFD simulations of wind turbine rotors with finite thickness and viscous effects, and demonstrated improved agreement for coned configurations. For a 15° coned rotor, the theoretically consistent two-point approach reduced power prediction errors by 8 %–10 % relative to simplified one-point implementations, resulting in noticeably closer agreement with CFD.

To assess the applicability of the correction, derived from unsteady thin airfoil theory, to highly swept, anhedral kite wings, VSM simulations were performed for the TU Delft V3 LEI kite and a ram-air kite using two alternative implementations: (i) a standard three-quarter-chord evaluation for both force magnitude and direction ($\mathrm{3}/\mathrm{4}c$), and (ii) Li's suggested approach in which the magnitude is evaluated at three-quarter-chord while the direction is evaluated at quarter-chord ($\mathrm{1}/\mathrm{4}c$). The validation dataset comprised wind tunnel measurements for both geometries, supplemented by 3D RANS CFD for the V3 LEI kite.

Figure B1 shows that, although both implementations predict C_L accurately and are in close agreement with the measurements and CFD, the $\mathrm{1}/\mathrm{4}c$ method systematically underpredicts C_D relative to the reference data. As a result, the predicted $C_{\mathrm{L}}/C_{\mathrm{D}}$ is consistently overestimated. Both kites exhibit similar trends, with the discrepancy being most pronounced for the V3 LEI kite. In contrast to the wind turbine results of Li et al. (2022), the standard $\mathrm{3}/\mathrm{4}c$ implementation yields better overall agreement with the experimental and numerical reference data, particularly for the key metric $C_{\mathrm{L}}/C_{\mathrm{D}}$.

The reason why the correction degrades, rather than improves, the drag prediction accuracy is not yet understood. Consequently, it cannot be concluded that the method proposed by Li et al. (2022) is fundamentally incorrect. Instead, based on the present validation results, it is inferred that, for the investigated geometries, the standard $\mathrm{3}/\mathrm{4}c$ implementation provides the closest agreement with the measurements and is therefore adopted throughout this work.

Li et al. (2022) additionally proposed a single-evaluation-point pitch-rate correction for unsteady cases. However, evaluating this correction requires validation data with known pitch histories, which are currently unavailable for these configurations; the present comparison is therefore limited to steady conditions.

Appendix C: Assuming fully turbulent flow

In 2019, Folkersma et al. (2019) showed with 2D RANS CFD that explicitly modelling boundary layer transition affected the aerodynamic predictions for Re < 2 × 10⁷ under a low freestream turbulence intensity of approximately 2 %. This finding motivated the subsequent use of natural-transition models in 3D RANS CFD studies of the V3 kite (Viré et al., 2020, 2022). Recent field measurements of the Kitepower B.V. V9 kite system, however, reported turbulence intensities of between 5 % and 30 %, with a median of about 8 % (Cayon et al., 2025; Kitepower, 2025), indicating a substantially more turbulent atmospheric environment than that assumed by Folkersma et al. (2019). Elevated freestream turbulence intensities are known to promote earlier boundary layer transition, both in terms of lower critical Reynolds numbers (Fransson et al., 2005; Zhang et al., 2020) and reduced critical roughness heights (Jain et al., 2024). Taken together, these observations suggested that the operational flow regime of the industrial kite system was closer to fully turbulent conditions than to the low-disturbance environment assessed in the original 2D study.

In addition to the atmospheric turbulence, the full-scale kite geometry introduced significant surface roughness elements not considered by Folkersma et al. (2019), which might trigger transition. On the suction side, the canopy was attached to the leading-edge (LE) tube by a zigzag-patterned stitching seam located on top of a forward step running along the span, while the pressure side featured a closing seam where the canopy wrapped around the LE tube. The obstacle height associated with the LE seam, representing the combined effect of the stitching and the step, was measured as k= 0.6 mm (Masure, 2025). In a 2D LEI wind tunnel experiment, a zigzag roughness height of 0.3 mm was found to already be sufficient to trip the boundary layer at Re = 5 × 10⁵ (van Lith et al., 2026). This is consistent with classical roughness-induced transition theory. The onset of transition by a roughness element of height k is commonly characterized as

$$\begin{array}{}\text{(C1)}& k={\displaystyle \frac{R{e}_{k,\mathrm{crit}}\mathit{\nu }}{U_k}},\end{array}$$

where ν is the kinematic viscosity and U_k the local velocity at the roughness height (Braslow and Knox, 1958). One finds that the roughness scales inversely with the velocity. In other words, the faster the flow is, and thereby the higher the Reynolds number, the smaller the roughness height needs to be. Given that 0.3 mm was found to be sufficient at Re = 5 × 10⁵, 0.6 mm was deemed sufficient at Re = 3 × 10⁶ – the approximately averaged Reynolds number encountered during the V3 flight (Cayon et al., 2025).

For the pressure side, van Lith et al. (2026) reported that the pressure-side flow close to the LE tube should be characterized by a different Reynolds number. Using the tube diameter as a characteristic length, an effective Reynolds number was found that is roughly 1 order of magnitude smaller than on the suction side. Masure (2025) found the closing seam height to be approximately 10 times larger than the suction-side zigzag seam, the direct opposite. In terms of the roughness Reynolds number, these two effects approximately compensated each other, suggesting that the pressure-side seam also operated above the critical roughness threshold.

Although detailed boundary layer measurements on the pressure side were not available, and thus a definitive experimental confirmation could not be provided, the combined evidence of (i) substantially higher atmospheric turbulence levels than in previous low-disturbance CFD studies, (ii) supercritical stitching seams on at least the suction side and likely also on the pressure side, and (iii) the trend towards higher Reynolds numbers with continued industrial upscaling provided the basis for assuming fully turbulent boundary layers in the presented 2D CFD simulations. This modelling choice was considered more representative of the operational flow conditions than not modelling the suction- and pressure-side seams, and letting the flow transition naturally over a smooth surface.

Appendix D: Machine learning model

Building the machine learning (ML) model involved several steps, illustrated in Fig. D1. In the blue boxes on the top-left, the data creation steps are shown. This involved defining a parameter grid, running 2D RANS CFD simulations, and splitting the data into training and test sets. Both datasets were preprocessed, and outliers were removed as explained in Sect. 4.

Figure D1 ML model setup flowchart, from data preprocessing to model evaluation.

In the purple block, the model setup is illustrated. Building the ML model began with standardizing the filtered training datasets using a standard scaler to ensure consistent weighting of all input parameters. Multiple regression algorithms were then evaluated using their default settings (Pedregosa et al., 2011), with the selected models chosen for their strong baseline performance reported by Pandala (2021). Linear regression (LR) was also included for reference. The hyperparameters of the selected algorithms were subsequently optimized using a grid-search procedure with k-fold cross-validation, using the coefficient of determination R² as the scoring metric. The dataset size was sufficient to allow reliable training–validation splits while still retaining a large and representative test set.

After hyperparameter tuning, the final models were evaluated on the independent test set. As the tree-based model response is piecewise constant in α, the outputs were additionally fitted with a PCHIP spline where a C¹-continuous polar was required.

Table D1 Overview of ML models, hyperparameters, and R² values.

To quantify the Reynolds-number sensitivity across the three per-Re models, the polars of the three representative test profiles of Table 7 were evaluated with each model. Table D2 reports C_l and C_d at α= 6°, representative of the reel-out operating condition of the V3 kite flight in 2019 (Cayon et al., 2025). Across the factor of 20 in Re, C_l increases by 1 %–8 % while C_d decreases by 14 %–26 %. The drag reduction between Re = 1 × 10⁶ and 5 × 10⁶ corresponds to a fitted power-law exponent of approximately −0.145, in close agreement with the turbulent skin-friction scaling C_f∼ Re${}^{-\mathrm{1}/\mathrm{7}}$, although the decay flattens towards the highest Reynolds number. Across the operational band $\mathit{\alpha }\in [-\mathrm{1},\mathrm{10}]$°, the three independently trained models remain closely spaced and consistently ordered; where the ordering is not strictly monotonic, the differences are within the regression uncertainty (Table 6). The largest sensitivity occurs at stall onset, where the stall angle increases with Re, consistent with expected boundary layer behaviour. The gradual, physically consistent variation across the sampled levels supports obtaining polars at intermediate Reynolds numbers via interpolation in log Re (Sect. 4).

Table D2 Reynolds-number sensitivity of the per-Re ML models: predicted C_l and C_d at α= 6° for the test profiles of Table 7.