Vortex model of the aerodynamic wake of airborne wind energy systems

. Understanding and modeling the aerodynamic wake of airborne wind energy systems (AWESs) is crucial for estimating the performance and deﬁning the design of such systems, as tight trajectories increase induced velocities and thus decrease the available power, while unnecessarily large trajectories increase power losses due to the gravitational potential energy exchange. The aerodynamic wake of crosswind AWESs ﬂying circular trajectories is studied here with vortex methods. The velocities induced at the AWES from a generic helicoidal vortex ﬁlament, trailed by a position on the AWES wing, are modeled with an expression for the near vortex ﬁlament and one for the far vortex ﬁlament. The near vortex ﬁlament is modeled as the ﬁrst half rotation of the helicoidal ﬁlament, with its axial component being neglected. The induced drag due to the near wake, built up from near vortex ﬁlaments, is found to be similar to the induced drag the AWES would have in forward ﬂight. The far wake is modeled as two semi-inﬁnite vortex ring cascades with opposite intensity. An approximate solution for the axial induced velocity at the AWES is given as a function of the radial (known) and axial (unknown) position of the vortex rings. An explicit and an implicit closure model are introduced to link the axial position of the vortex rings with the other quantities of the model. The aerodynamic model, using the implicit closure model for the far wake, is validated with the lifting-line free-vortex wake method implemented in QBlade. The model is suitable to be used in time-marching aero-servo-elastic simulations and in design and optimization studies.


Introduction
Airborne wind energy (AWE) refers to the field of wind energy in which tethered airborne systems are used to harvest wind power at high altitudes. Airborne wind energy systems (AWESs) are typically classified based on their flight operations, which can be crosswind, tether-aligned and rotational as described by Vermillion et al. (2021). Electric power is generated with onboard wind turbines and transferred to the ground through the tether (Fly-Gen) or generated directly on the ground by a moving or fixed ground station (Ground-Gen). This work focuses on crosswind AWESs, and the results are applicable to both Ground-Gen and Fly-Gen systems which are characterized by a single wing. Loyd (1980) first derived the power equation for crosswind AWESs, based on given lift and drag coefficients of the AWES. To properly evaluate system performance, spe-cial care should be given to the estimation of the drag coefficient, as it includes contributions from the tether drag (see Trevisi et al., 2020a, for more details) and from the induced drag. The latter is the result of the induced velocities generated by the lifting surface and its wake on the lifting surface itself. For wind-harvesting devices, the induced velocities effectively reduce the incoming wind and thus the available power. Generating high aerodynamic lift could be counterproductive for power generation, as induced velocities grow consequently. For flying wings, the induced velocities effectively rotate the inflow field of an induced angle. As lift is, by its definition, perpendicular to the local velocity, it is rotated by this induced angle, generating an aerodynamic force component parallel to the motion, namely the induced drag. High aerodynamic lift then generates high induced drag, which needs to be compensated by the propulsive thrust.
F. Trevisi et al.: Vortex model of the aerodynamic wake of AWESs As AWESs are wind-harvesting devices, the evaluation of the induced velocity is of extreme importance to properly estimate the power production. The effect of the induced velocity can be taken into account by reducing the incoming wind, as typically done for wind turbines, or by including them in the induced drag estimation, as typically done for flying wings.
Learning from conventional wind energy, the induced velocities for AWE are estimated in the literature with momentum-based and vortex-based methods. High-fidelity computational fluid dynamics (CFD) methods are typically used to study the wake characteristics and validate lowerfidelity models.
Initially, the applicability of momentum methods for AWE was doubted (Loyd, 1980;Archer, 2013;Costello et al., 2015), as the definition of swept area is not settled but by intuition way larger than the AWES wing. Recently, however, momentum methods have been developed for AWE, with De Lellis et al. (2018) and Kheiri et al. (2018Kheiri et al. ( , 2019 independently generalizing momentum theory to compute the induced velocity. This is derived by equating the aerodynamic lift of the AWES with the thrust applied to the annulus swept by the wing, with the momentum formulation. The power equations for Fly-Gen and Ground-Gen AWESs are then derived by reducing the incoming wind of the induced velocities. They find that Fly-Gen AWESs have higher power generation potential compared to Ground-Gen. Akberali et al. (2021) extend the derivation for AWESs with an elevation and a side-slip angle. Kaufman-Martin et al. (2021), Kheiri et al. (2022) and Karakouzian et al. (2022) study, with the aim of understanding the interaction of AWESs in a wind farm, the downstream AWES wake shape and characteristics with momentum and mass conservation considerations, finding a good agreement with CFD results. Gaunaa et al. (2020) point out that using momentum theory to evaluate the induction of an AWES, which is described by 3D polars, is not physically consistent. Momentum theory is indeed used in wind turbine aerodynamics to compute the local velocity triangle of the airfoil (2D polars) in the wind turbine blade. If the AWES is described by 3D polars (i.e., the drag coefficient includes the induced drag coefficient that the AWES would have in a forward flight) in the momentum formulation, then a part of the wake would be counted twice. If, instead, momentum theory is used to evaluate the velocity triangle of an airfoil (2D polars) in the AWES wing, then root and tip corrections would be needed to take into account that the rotor is not a disc built up of an infinite number of wings but one single wing. The root and tip corrections for AWESs would however differ largely from conventional wind turbines corrections, as these are developed for blades extending almost to the rotation axis and need to be reworked. Gaunaa et al. (2020) then build a vortex-based engineering model, which is physically consistent, and it is tuned with CFD simulations. Another vortex method, based on the vortex tube model, is developed by Leuthold et al. (2019) for an AWE system composed of more wings flying in the same annulus to compute the induction at the AWESs mid-span. This model takes as inputs the system thrust coefficient, relative radius and reel-out factor.
As the field counts a limited number of prototypes with a limited number of flying hours, the benchmark of engineering models is typically done with higher-fidelity codes, instead of experiments. For aerodynamic models, CFD studies then represent the reference. Haas and Meyers (2017) describe the wake characteristics for a given aircraft in a Fly-Gen and Ground-Gen circular path with a large-eddy simulation (LES) setup. Aerodynamic forces, applied with an actuator line technique, are computed to impose an induction of 1/3 at the kite location. Haas et al. (2019) further develop the same LES framework by including an optimal control problem for Ground-Gen AWESs in non-turbulent and turbulent sheared inflow conditions. Mehr et al. (2020) investigate the aerodynamic interaction of the onboard wind turbines with the main wing in a crosswind circular maneuver with a viscous vortex particle method. The main wing wake is however removed one span downstream from the trailing edge, not resolving the helicoidal wake. Branlard et al. (2022) extend the AeroDyn module of OpenFAST through the open-source lifting-line vortex code OLAF to support arbitrary collections of wings, rotors and towers. Complex geometries, such as AWESs, can then be handled by OLAF.
The present work builds on the considerations from Gaunaa et al. (2020) to develop a vortex-based engineering model, which is physically consistent. The use of vortex methods, described by Branlard (2017), allows for modeling the induction at the AWES due to its own wake. The goal of this work is to analytically study the wake produced by AWESs flying circular trajectories by modeling how geometrical and aerodynamic quantities influence the induced velocities. This model can be used to find how the induced velocities vary along the wing span and to refine the power equations.
An analytical model of the helicoidal wake will inform about the validity bounds of the assumption of a straight wake employed in the engineering models (Bauer et al., 2018;Trevisi et al., 2020b, among others) and in computational studies focusing on the aerodynamic characterization of the AWES (Viré et al., 2020;Castro-Fernández et al., 2021, among others).
Modeling the induced velocities along the wing span will instead improve the engineering aerodynamic models for load estimation (Damiani et al., 2019) and the flight dynamic models (Trevisi et al., 2021a).
A refined power equation, including the contribution from the helicoidal wake, will inform the designer in an intuitive way about the effects of changing geometrical and aerodynamic quantities on the overall performance. This could be used in low-fidelity system design and optimization studies (Bauer et al., 2018;Trevisi et al., 2021b, among others). The power equation refinement will be performed in an upcoming work to limit this paper's scope and length. This paper is organized as follows: in Sect. 2, the modeling of the induction produced by a helicoidal vortex filament is introduced. The helicoidal vortex filament can be described as a near vortex filament, modeled by a half revolution of the filament in the rotor plane, and a far vortex filament, modeled by an infinite series of vortex rings. In Sect. 3, the near-wake model and the relative induced drag coefficient are derived. In Sect. 4, the far-wake model, which depends on the radial (known) and axial (unknown) position of the vortex rings, is derived. An explicit and an implicit model are proposed to find the axial position. In Sect. 5, the model is validated with QBlade, which uses a lifting-line free-vortex wake method. In Sect. 6, the main conclusions are discussed.

Modeling of a helicoidal vortex filament
In this section, the modeling framework of the induced velocities produced by a non-expanding helicoidal vortex filament is introduced. As the induced velocities at a given point are found by integrating the contributions of the vortex filaments composing the aerodynamic wake, they are first studied in this section. The formulation introduced in this section will be used to derive a near-and a far-wake model in the next sections.
AWESs flying crosswind typically follow a circular or a figure-of-eight pattern. The present work analyzes the wake produced by circular trajectories, which are contained in a plane perpendicular to the incoming wind, as in Fig. 1. It is here assumed that the AWES has a constant rotational speed, and its wings are on the rotational plane. This condition is also considered by the other studies concerning AWES wakes (e.g., De Lellis et al., 2018;Kheiri et al., 2018;Gaunaa et al., 2020;Kaufman-Martin et al., 2021). Moreover, Trevisi et al. (2022a) show that these are the optimal trajectories for a Fly-Gen maximizing thrust power with constant inflow. Optimal trajectories for a Fly-Gen maximizing electric power including wind shear can be understood as a perturbation of this solution. Additionally, one can study the flight mechanics of the AWES for this condition (Trevisi et al., 2021a).
Referring to Fig. 1, the AWES moves along a circular trajectory in the plane (e 1 , e 2 ) with radius R 0 . The versor e 3 points upwind, and therefore the wind velocity is v w = −v w e 3 . For Ground-Gen AWESs, the relative wind velocity at the AWES is is the reelout velocity of the tether, assumed to be only along the axial direction. For Fly-Gen AWESs, the relative wind velocity at the AWES coincides with the wind speed v r = v w , as the tether is fixed at the ground station. The vortices trailed by the AWES are transported downwind by the wind and have a helicoidal shape. Following the assumptions of Prandtl and Goldstein in the derivation of the tip-loss correction for wind turbines, no distortion and no expansion of the wake are assumed (Branlard, 2017).
The geometry of the helicoidal vortex filament is shown in Fig. 2a. The radius of the filament is R f , and the pitch is h f . The induction is evaluated at a generic point p j , with radius R j . With these definitions, the radial difference R between the evaluation point and the filament is and it is normalized with the evaluation radius as (2) Note that η = 0 when R = 0, η = 1 when R j → +∞ and η → −∞ when R j → 0. The normalized position of the vortex filament with respect to the mid-span turning radius R 0 (see Fig. 1) is indicated as where y f indicates the position of the vortex filament along the wing span direction. The filament radius R f can then be expressed as The helix can be modeled as where R 0 R j = 1−η 1−η 0 (Eq. 4), θ ∈ [0, ∞[ is the angular parameter of the helix and it is assumed that the helix pitch of any vortex filament h f trailed by the wing is equal to the helix pitch h 0 of the vortex filament at the wing center. The normalized torsional parameter λ 0 (Branlard, 2017) is the ratio  between the projected circumference length 2π R 0 and the helix pitch h 0 of the vortex filament: The induced velocities w j,f produced by the filament of vorticity at a point p j is found using the Biot-Savart law: where and r = p j −l. As the vortex filaments are trailed by the wing and the evaluation point is typically on the aircraft itself, p j is modeled as where θ j and φ j define the 3D position of the evaluation point and are assumed to be small (see Fig. 2a). r can then be expressed as When looking for the axial induced velocity w j,f = w j,f · e 3 , the third component of the numerator of Eq. (7) becomes where the subscript ∼ indicates a periodicity in θ . The norm of r squared, removing the second-order terms in θ j and φ j and using the trigonometric relation cos(θ − θ j ) ≈ cos(θ ) + sin(θ )θ j , can be expressed as where a distinction between the periodic D ∼ and the nonperiodic part D / is made. Finally, the Biot-Savart law (Eq. 7) can be written considering Eqs. (11) and (12) as where the term outside the integral models the induced velocity produced by a semi-infinite straight vortex filament pointing in the −e 1 direction ( Fig. 2b) on the point p j (θ j = 0, φ j = 0) (Eq. 9). The integral, denoted with ϒ, can be understood as a shape factor multiplying it.

By splitting the integration interval as
, the integral can be rewritten as an infinite summation of integrals. By properly rewriting the non-periodic terms, Eq. (13) becomes In Eq. (14) a distinction between near and far vortex filament is made, noting that the near vortex filament (i.e., the first integral) models the induction produced by the first half revolution of the filament, whereas the far vortex is modeled through the infinite summation of integrals.
The integrals modeling the near and the far vortex filament have terms proportional to θ in the denominator. These terms physically represent the contribution to |r| from the projection of r along e 3 , and they have a maximum value at θ = ±π , which is also where the projection of r in the (e 1 , e 2 ) plane is largest. This suggests that the overall contribution of the terms proportional to θ to the integral should be limited.
By neglecting the terms proportional to θ, the model of the helicoidal vortex filament in Fig. 2a reduces into the model of Fig. 2b (i.e., half a vortex ring plus a semi-infinite vortex ring cascade). Note that, with the split of the integration interval introduced in Eq. (14), the axial distance of the vortex rings corresponds to the helix pitch and the integration bounds of the far-wake integrals are symmetric.
The axial velocity at point p j , induced by the idealized helicoidal vortex filament, is then approximated as which in its extended form is In Fig. 3, the shape factor ϒ of the helicoidal vortex filament in Fig. 2a (solid, Eq. 13) and of the vortex ring system of Fig. 2b (dashed, Eq. 15), computed for typical values of 1 1−η 0 1 λ 0 , is shown. A detailed analysis of λ 0 for the aerodynamic wake of AWESs is given in Sect. 4. Values of η in the figure are shown up to η = −5, which corresponds to an evaluation radius of R j = 1 6 R f . The velocities induced by the helicoidal vortex filament are well approximated by the vortex ring system. When looking at the velocities induced at the AWES by the vorticity trailed by AWES itself, values of η are expected to be close to 0, which is when the two models are the most similar. These results justify the decision of neglecting the terms proportional to θ in the modeling.
The integrals involved in Eq. (16) now have a closed-form solution, which will be used to derive a near-and a far-wake model in the next sections. This expression is derived by assuming that the angular distance of evaluation point p j from the filament origin (i.e., θ j and φ j ) is small, which is the case when the velocities induced by the vortex filaments trailed by the AWES are evaluated on the AWES itself, and by approximating the helicoidal vortex filament in Fig. 2a with the vortex ring system of Fig. 2b to compute the induction, which is a good idealization when the normalized torsional parameters λ 0 of the helicoidal vortex filament are large.

Near-wake model
In Sect. 2, the axial induced velocity produced by a helical vortex filament is modeled as half a vortex ring plus a semiinfinite vortex ring cascade, as in Fig. 2b. In this section, a near-wake model is developed based on this idealization. The velocities induced at the AWES by the near vortex filaments trailed by the AWES itself are studied here so that values of η (Eq. 2) are expected and consequently assumed to be small. This section is organized as follows: in Sect. 3.1, the closed-form solution of the velocity induced by the near vortex filament and its linearized version are derived. In Sect. 3.2, the circulation of an elliptic wing flying in a circular trajectory is introduced to find the related trailed vorticity. In Sect. 3.3, the velocities induced by the trailed vortices on the lifting line of the elliptic wing are found. This allows for defining an induced drag coefficient modeling the near wake.

Velocity induced by the near vortex filament
The axial velocity at point p j induced by the near vortex filament of intensity d is (Eq. 16) By changing the integration variable to θ * = θ − θ j and performing the integration, the shape factor of the near vortex filament ϒ n has a closed-form solution: where K θ * 2 m is the incomplete elliptic integral of the first kind and E θ * 2 m is the incomplete elliptic integral of the second kind with m = 4 (η−1) η 2 . As θ j is assumed to be small throughout all this work, Eq. (18) can be linearized with respect to θ j as where K π 2 m and E π 2 m are the complete elliptic integral of the first and second kind respectively. Equation (19) can be linearized with respect to η, as its values for AWESs are expected to be small. Therefore, ϒ n becomes (20) Figure 4 shows the trend of ϒ n obtained with Eq. (18) (solid lines) and with the linearized version (Eq. 20) (dashed lines). The shape factor of the near vortex filament ϒ n can be interpreted as a corrective factor to the solution of the straight vortex filament due to the filament curvature. If the evaluation point is in the same radial direction of the vortex filament origin (i.e., for θ j = 0), the following regions of η can be analyzed: η → 1 is obtained for R j R f and ϒ n → 0, meaning that no velocities are induced; η > 0 is obtained for R j > R f and ϒ n < 1, meaning that fewer velocities are induced compared to the straightfilament case; and ϒ n → 1, meaning that the filament curvature is negligible and the solution coincides with the straightfilament case; η < 0 is obtained for R j < R f and ϒ n > 1, meaning that more velocities are induced compared to the straightfilament case; η → −∞ is obtained for R j R f and ϒ n → π , as the evaluation point is close to the center of the half circle.
If the evaluation point is not in the same radial direction of the vortex filament origin (i.e., for θ j = 0), the solution tends to ϒ n (θ j = 0) when η → 1 or η → −∞ because the evaluation point is far from the filament, so the effect of θ j vanishes. For θ j > 0, the evaluation point is downstream of the filament origin (see Fig. 2a) and ϒ n → 2 when η → 0, meaning that twice the velocities induced by the semi-infinite straight filament (equivalent to one infinite straight filament) are found. For θ j < 0, the evaluation point is upstream of the filament origin and ϒ n → 0 when η → 0, meaning that no velocities are induced.

Circulation trailed by an elliptic wing
Straight elliptical wings in forward flight have minimum induced drag and therefore are here chosen to analyze the near wake. To express the ratio between the wing coordinate y j = R j − R 0 and the turning radius R 0 , the parameter η j is introduced as for which at the outer wing tip y j = b/2 takes the value κ 0 : which is named the inverse turning ratio in this work. The position of the trailed vorticity y f = R f −R 0 is defined through the angular position α as with η 0 defined in Eq.
(3). The aerodynamic lift per unit length at the location of the trailed vorticity on the wing can be computed with the Kutta-Joukowski relation as where ρ is the air density, (y f ) is the bound circulation, and the velocity u is composed by the rigid-body velocity u 0 R 0 +y f R 0 = u 0 (1 − η 0 ) and the relative wind velocity v r : which can be approximated with the rigid-body velocity u(y f ) ≈ u 0 (1−η 0 ) for high wing speed ratios λ = u 0 /v r and a small η 0 . For the AWES to be in equilibrium, the total forces and moments needs to be null. Assuming that the roll moment L is only due to the main wing lift 1 , L is The bound circulation can be taken to set the lift distribution to be symmetric along the wing span: where is the symmetric circulation the wing would have in a straight motion and η 0 is assumed to be small. Note that loads are typically distributed along the wing span by means of ailerons, which create discontinuities in the lift distribution. The lift distribution is modeled here to be a continuous function for simplicity. The symmetric circulation the elliptic wing would have in forward flight is (Anderson, 2017) = 0 sin(α) = 2b where C L is the wing lift coefficient, AR is the wing aspect ratio and w is the induced velocity the same wing would have in forward flight. The trailed vorticity, according to Helmholtz' law, is the derivative of the bound vorticity with respect to the wing span coordinate: 3.3 Induced drag coefficient due to the near wake for an elliptic wing The axial induced velocity w n j at a given point p j along the lifting line (θ j = 0) due to the near wake is found by integrating the effects of the trailed vorticity γ along the wing span as The induced velocities on the lifting line can be found numerically and fitted between y j = −b/2 and y j = b/2. The induced velocity at a point on the lifting line given by the near wake is In Fig. 5, a comparison of the analytic approximation (Eq. 31) with the solutions obtained numerically is shown. The induction computed with straight vortex filaments or with the near vortex filaments are almost equivalent. This suggests that the development of lifting-line methods assuming straight trailed vortex filaments for time simulators (e.g., Damiani et al., 2019) or for flight dynamics analyses (e.g., Trevisi et al., 2021a) are a good approximation of the near wake for elliptic wings and for the inverse turning ratio κ 0 considered in the example or lower values. Recall that assuming straight trailed vortex filaments is equivalent to approximating ϒ n (η, θ j = 0) in Eq. (30) with ϒ n = 1. ϒ n goes to 1 for values of η going to 0, as shown in Fig. 4. Larger values of η in modulus are obtained for larger values of κ 0 . For an increasing κ 0 , the approximation of straight trailed vortex filaments is then expected to perform less well. Moreover, Figure 5. Induced velocities computed analytically w n (Eq. 31) and numerically with curved wake (Eq. 20 is used) and no-roll circulation w ∪ ( ), with straight wake and no-roll circulation w ( ), and with straight wake and symmetric circulation w ( ).
for wing types with more trailed vorticity near the wing tips (i.e., where values of η are larger in modulus), the approximation of straight trailed vortex filaments is also expected to perform less well.
The induced change in the angle of attack for the local profile is Finally, the wing drag coefficient can be found by integrating the projection of the local lift coefficient along the velocity direction as where κ 0 is assumed to be small. The induced drag of an elliptical wing, due to the near wake, is then found to be similar to the induced drag of the same wing in forward flight. The glide ratio related to the near-wake drag coefficient is 4 Far-wake model In Sect. 2, the axial induced velocity produced by a helical vortex filament is modeled as half a vortex ring plus a semiinfinite vortex ring cascade, as in Fig. 2b. In Sect. 3, a nearwake model is built based on this idealization, while in this section a far-wake model is developed. This section is organized as follows: in Sect. 4.1, an approximation of the velocities induced by two semi-infinite vortex ring cascades with opposite intensity is given as a function of their radial and axial distance. In Sect. 4.2, a drag coefficient modeling the far wake for an elliptic wing is found. In Sect. 4.3 an explicit closure model and in Sect. 4.4 an implicit closure model to find the vortex rings axial distance are introduced.

Velocity induced by the far vortex filaments
In the far wake, the wake is assumed to be already rolled up into two separate vortices, one for each of the wing tips. The velocities induced by the far wake are assumed to be constant over the wing span; therefore they are evaluated at the wing center of R j = R 0 . The velocity induced by the cascade of vortex rings of intensity on the wing center is where φ j in Eq. (16) is neglected. The integral has a closedform solution, so ϒ f z,k takes the form with m = 4(η 0 −1) The outer and inner trailed vortexes are assumed to be located at R f = R 0 ± y v such that R = ∓y v and η 0 = ∓η v = ∓ y v R 0 with an intensity of = ± 0 respectively. The velocity at the wing center induced by the two vortex ring cascades is The summations can be solved numerically for different values of λ 0 and η v , and its solution can be fitted as a function of these two parameters. The approximated solution of the summation is Figure 6 shows the comparison between the solution obtained numerically and the approximation given in Eq. (38).

Induced drag coefficient due to the far wake for an elliptic wing
For an elliptical wing, the rolled-up trailed vortices are located approximately at the center of trailed vorticity y v = π 8 b (Gaunaa et al., 2020) for the outer and inner wing respectively such that their non-dimensional radial coordinate is η v = πb 8R 0 = π 4 κ 0 . The approximate solution for the induced velocity, using the fitted solution of Eq. (38) into Eq. (37) and considering 0 = 2b u 0 C L π AR , is and the relative induced drag coefficient due to the far wake is The glide ratio related to the far-wake drag coefficient is As one of the final goals of this work is to model the system glide ratio G to refine a power equation by properly including the contribution from the wake, the correct drag coefficient shall include contributions from the near and far wake. The system glide ratio is therefore where C D0 is the summation of the aircraft viscous and pressure drag C D,v and the equivalent tether drag C D,t = C ⊥ D t L t 4A , where C ⊥ is the drag coefficient of the tether section, D t and L t are the tether diameter and length and A is the main wing area (Trevisi et al., 2020a).
The system glide ratio can be then seen as the sum of G 0 , G n and G f connected in parallel (as resistors connected in parallel in an electric circuit) as Equation (43) models the system glide ratio of an AWES with an elliptic wing under the assumption of no wake expansion for a generic normalized torsional parameter λ 0 . A closure model is then necessary to link λ 0 with the other aerodynamic and geometric parameters of the model. Two closure models are proposed in the following two sections.

Explicit closure model for the normalized torsional parameter
The first closure model assumes that the axial velocity of the vortex filaments is equal to the axial velocity at the AWES, similar to what assumed by de Vaal et al. (2014) in the derivation of a free-vortex wake model for wind turbines based on vortex ring elements. The axial velocity of the vortex filaments is assumed to be v r (1 − a z ), where a z is the axial induction at the AWES wing center. The helix pitch h 0 can be approximated with the distance covered by the vortex filaments moving downwind in the revolution period. The revolution period is the ratio between the circumference length and the AWES mid-span tangential velocity 2πR 0 u 0 : where λ = u 0 v r is the wing speed ratio. Considering the definition of the normalized torsional parameter λ 0 = 2πR 0 h 0 (Eq. 6), λ 0 can be linked to the wing speed ratio λ as .
The axial induction is the sum of the near-wake and the farwake induced velocities, normalized with the relative wind speed: By inserting the axial induction a z into Eq. (45), the wing speed ratio λ is which can be rewritten as This equation links λ (wing speed ratio) to λ 0 (normalized torsional parameter of the helicoidal wake) and to the induced change in the angle of attack produced by the near wake C L π AR and by the far wake 1 4π C L πAR κ π/2 0 λ 3/2 0 . For AWES flying crosswind in steady state the propulsive lift is balanced by the drag such that the force balance along the AWES longitudinal direction is null. The glide ratio G is then equal to the wing speed ratio λ (Fig. 7a). By setting equal Eqs. (43) and (48), G 0 is found to be equal to λ 0 . This means that the normalized torsional parameter of the helicoidal wake λ 0 is equal to the glide ratio related to the constant drag term G 0 . A graphical representation of this equality is given in Fig. 7b. λ 0 is the ratio between the AWES longitudinal velocity and the wind velocity reduced by the induced velocities at the AWES (Eq. 45). If the induced velocities due to the near and far wake are subtracted from the incoming wind, the induced drag related to the near and far wake should not be considered in the drag evaluation. Therefore, the only remaining drag term is related to the drag at zero lift. For the AWES to be in equilibrium along the longitudinal direction, λ 0 needs to be equal to G 0 . Note that the two sketches in Fig. 7 show the same physics. In Fig. 7a the effect of the induced velocities is taken into account by including induced drag terms, as typically done for flying wings. In Fig. 7b the effect of the induced velocities is taken into account by reducing the incoming wind, as typically done for wind turbines. The resultant longitudinal velocity u 0 is equal for the two sketches.
The system glide ratio (Eq. 43) can then be formulated by considering λ 0 = G 0 , leading to a new definition: The system glide ratio is now also dependent on the flight trajectory radius through the inverse turning ratio κ 0 . For a large turning radius (i.e., small κ 0 ) the contribution from the far wake vanishes and Eq. (49) coincides with the glide ratio expression for the same wing in straight motion. Figure 8 shows the ratio between the far-and near-wake induced drag (i.e., the expression 1 4π κ π/2 0 G 3/2 0 ). AWES designs tend to have high values of G 0 , obtained by minimizing C D0 and operating at high lift coefficient C L , and to have a large κ 0 (i.e., small turning radius R 0 ). Small turning radii minimize the vertical height of the trajectory and therefore reduce the potential energy exchange over the loop, reducing power losses (Trevisi et al., 2022a). The induction of the far wake however penalizes solutions with high G 0 and high κ 0 , highlighting the need of a trade-off between different disciplines in the design. These considerations, among others, justify the development of the multidisciplinary design, analysis and optimization framework T-GliDe (Trevisi et al., 2022b).
In Sect. 5, a comparison between the results produced by this closure model and by a lifting-line free-vortex wake method, implemented in QBlade, is performed. Unfortunately, the comparison does not show a satisfactory agreement between the two methods for high-loading cases. Therefore, an implicit closure model is proposed in the next section.

Implicit closure model for the normalized torsional parameter
The second proposed closure model for the normalized torsional parameter estimates the wake velocity by considering the axial and the radial induced velocities at the AWES. It is assumed that the velocity of the vortex filaments in the far wake is equal in modulus to their velocity when they are trailed (i.e., at the AWES). The velocity at the AWES is composed by the axial velocity v r (1−a z ) and by the radial velocity v r a r . Note that the radial velocity is only due to velocities induced from the far wake. The helix pitch h 0 can then be found by considering that the axial velocity of the far vortex filaments is v r (1 − a z ) 2 + a 2 r : The normalized torsional parameter can then be written as where a z is given in Eq. (46) and a r for an elliptic wing can be approximated as (see Appendix A for the derivation) The normalized torsional parameter (Eq. 51) can be found iteratively by considering the axial (Eq. 46) and the radial (Eq. 52) induction and the definition of the glide ratio (Eq. 43), which is equal to the wing speed ratio G = λ when the AWES is in equilibrium.
The estimation of radial induced velocities is important not only to compute the axial velocity of the far vortex filaments but also to estimate the AWES aerodynamic forces and moments. In the case of a straight wing with the wing span aligned with the radial direction (as in this paper), the radial induced velocities do not generate any aerodynamic force on the wing, as they are aligned with the span. On the contrary, the radial induced velocities generate aerodynamic forces and moments in the case of a straight wing with the span not aligned with the radial direction (i.e., the AWES is yawed or rolled) or in the case of a wing characterized by a dihedral and/or sweep angle. Moreover, the radial induced velocities also influence the angle of attack of the vertical stabilizer.

Comparison with a lifting-line free-vortex wake method
To validate the newly introduced aerodynamic model, a comparison with higher-fidelity code is needed. QBlade is chosen for the comparison, as it uses a thoroughly validated liftingline free-vortex wake method (Marten et al., 2015) for the aerodynamic modeling of conventional wind turbines. The wake is modeled explicitly with Lagrangian vortex elements, resulting in a fully resolved velocity distribution around the AWES.

Lifting-line free-vortex wake setup
As QBlade is being used to validate the new model, the free-vortex wake model is idealized accordingly. The airfoils are described by idealized polars with constant drag coefficient C D0 and a constant slope of the lift coefficient with respect to the angle of attack of C l,α = 2π . The AWES is modeled as a horizontal axis wind turbine with one blade and no tower. The hub radius of the wind turbine is set to be equal to the radial position of the inner wing tip. The elliptical wing is then defined in the blade definition. Large structural rigidity is given to avoid aeroelastic phenomena. Gravity is set to zero to look for steady results. The wing is discretized with 20 aerodynamic panels. For a given setup (C D0 , κ 0 , AR), different lift coefficients can be obtained by pitching the wing. The wing is free to move on the circular trajectory so that its velocity is determined by the force equilibrium along the motion direction. Once quantities reach a constant value as a function of time, the wing speed ratio, which is equal to the glide ratio, is estimated from the rotational speed. The wing lift coefficient is obtained by knowing the thrust force produced and the rotational speed. In Fig. 9, the wake shown in the graphic interface of QBlade is reported for a case with AR = 20, κ 0 = 0.15, C D0 = 0.05 and 0 • of pitch, which leads to C L = 0.55. A plane, defined by the AWES position and wake center, is used to visualize the axial velocity.

Comparison of the normalized torsional parameter
The first comparison concerns the downstream position of the vortices. The helix pitch can be found with h 0 = 2πR 0 λ 0 (Eq. 6), where λ 0 is the normalized torsional parameter. In this paper, an explicit (Sect. 4.3) and an implicit (Sect. 4.4) closure model for the normalized torsional parameter are proposed. Figures 10 and 11 show the downstream vortex position, obtained by taking the curl of the velocity field, in the plane containing the AWES (see Fig. 9 for the plane illustration) for a setup with AR = 20, κ 0 = 0.15 and C D0 = 0.05. Figure 10 shows a low-loading case with C L = 0.55, while Fig. 11 shows a high-loading case with C L = 1.3. The blue vertical  lines highlight the downstream position of the vortices predicted by the explicit model, while the black line is that of the implicit model. The dotted red lines highlight the radial position of the vortices assumed by the analytical model. For the low-loading case, both models accurately predict the vortex position. For the high-loading case, the implicit model greatly outperforms the explicit model, showing that the radial induced velocity at the AWES significantly contributes to the velocity of the far vortex filaments. In the high-loading case, the wake expands slightly due to the radial induced velocities. The wake expansion is neglected in the analytical modeling.

Comparison of axial induced velocities
In this section, the axial induced velocities computed with the analytical model just introduced and the implicit closure model are compared with the results of QBlade and two literature models related to the far wake. The high-loading case introduced in the previous section (AR = 20, κ 0 = 0.15, C D0 = 0.05 and C L = 1.3) is here used to show the results.
The bound circulation is retrieved from QBlade, and the trailed vorticity γ can be obtained by differentiating the bound vorticity with respect to the spanwise coordinate. The induced velocities due to the near wake w n can then be found by integrating Eq. (30) numerically. The induced velocities due to the far wake w f are computed with Eq. (39). The first model from the literature is the far-wake model by Gaunaa et al. (2020): where C L is taken from QBlade and G is found iteratively. This model is coupled to the near-wake model introduced in this paper, as the near-wake model from Gaunaa et al. (2020) is tuned with CFD results and does not agree with the proposed near-wake model. The second model from the literature, used to estimate the induction a f m. , is the momentum-based model by Kheiri et al. (2018): where C D,k = C D0 + C 2 L πAR and C L is taken from QBlade. Also this model is coupled to the near-wake model introduced in this paper such that a f m. can be seen as the induction applied to the 3D polars of the AWES. Figure 12 shows the induction at the lifting line computed with QBlade, with the present model and with two models from the literature. The momentum-based model overestimates the induced velocities because a part of the wake is counted twice if the AWES is described by 3D polars, as explained by Gaunaa et al. (2020). The proposed model and the far wake from Gaunaa et al. (2020) coupled to the proposed near-wake model accurately predict the induction at the lifting line.

Comparison of radial induced velocities
In this section, the radial induced velocities computed with the analytical model (Eq. 52) are compared with the results of QBlade. Figure 13 shows the radial velocities for the high-loading case (AR = 20, κ 0 = 0.15, C D0 = 0.05 and C L = 1.3). A good agreement between the QBlade output and the analytic solution is found.

Comparison of axial induced velocities at the AWES tail
In this section, the axial induced velocities in the plane normal to the wing span passing by the center of the wing predicted with the proposed model and with a literature analytical model are compared with the results from QBlade. The induced velocity in the symmetry plane of the AWES due to the trailed vorticity in the near wake can be found as a function of θ j with where γ is taken from QBlade. To find the velocity field, the velocities induced by the bound vorticity and from the far wake are considered too. Phillips et al. (2002) study the velocities induced at the tail from a wing in forward flight. For an elliptical wing with no sweep, the induced velocities can be reformulated as As the velocities induced by the near wake are similar to the velocities induced by the wing in forward flight, the contribution from the far wake is added to this expression. In Fig. 14, the induced change in the angle of attack found in QBlade is compared with the two solutions for the high-loading case (AR = 20, κ 0 = 0.15, C D0 = 0.05 and C L = 1.3). The proposed model is in good agreement with  QBlade, and the model by Phillips et al. (2002) can still capture the main trend. These models can then be used to size the horizontal stabilizer.

Comparison of glide ratios
The estimation of the system glide ratio is crucial for power production estimation. Figure 15 shows the glide ratio as a function of the lift coefficient for a case with AR = 20, κ 0 = 0.15 and C D0 = 0.05 computed considering straight wakes (G 0 G n ), the explicit model (G 0 G n G f expl , Eq. 49), the implicit model (G 0 G n G f impl , Sect. 4.4) and the far-wake model from Gaunaa et al. (2020) coupled to the proposed near wake (G 0 G n G f G. ). The squares are the solutions obtained with QBlade. The explicit closure model results in a conservative glide ratio estimation for high loading, as it predicts the helix pitch to be smaller than what is found in QBlade. Using the implicit closure model and the model from Gaunaa et al. (2020) for the far wake results in accurate prediction of the results of QBlade.
In Appendix B, more validation cases are reported. Generally, the implicit closure model results in an underestimation of the QBlade glide ratio estimation, while the model from Gaunaa et al. (2020) in an overestimation. The two methods are in accord and have the same computational cost.

Conclusions
In this work, a detailed analysis of the airborne wind energy system (AWES) aerodynamic wake is carried out with vortex methods and models to find the induced velocities at the AWES.
Under the assumption of steady crosswind operations and a non-expanding wake, the expression for the velocities induced in the neighborhood of the AWES from a helicoidal vortex filament, trailed by a position on the AWES wing, is divided into an expression modeling the near vortex filament and one for the far vortex filament.
The near vortex filament is modeled as the first half rotation of the helicoidal filament, where the axial component of the filament is neglected. The velocity induced by the near vortex filament is expressed in terms of elliptic integrals, and it is linearized to a more intuitive expression. The induced drag coefficient modeling the near wake, built up from the contributions of the near vortex filaments, is found to be similar to the drag coefficient the same wing would have in forward flight. This suggests that models assuming straight trailed vortex filaments are a good approximation of the near wake only.
The far wake is modeled as two semi-infinite vortex ring cascades with opposite intensity. The related induced velocities depend on the radial position of the rings, which is known, and on the axial distance of the rings, which is unknown. Two closure models are proposed to link the axial distance of the rings with the other physical quantities of the model. An explicit model is derived by assuming that the vortex filaments move downstream with the axial velocity at the AWES center. An implicit model is derived by assuming that the far-wake vortices move downstream with a velocity equal to the velocity, in modulus, at the AWES center. To find the modulus of the velocity at the wing center, the radial velocity induced by the far wake is also derived.
To validate the newly introduced model, a comparison with a lifting-line free-vortex wake method, implemented in QBlade, is performed. A good agreement between the implicit model and the free-vortex wake results is found. The proposed aerodynamic model will be used to refine the power equations of Ground-Gen and Fly-Gen AWESs in an upcoming work. The model is suitable for time-marching aero-servo-elastic simulations of AWESs. Indeed, the near wake can be modeled with engineering vortex-based aerodynamic models, developed to take into account the curvature of the trailed vortex filaments, and the far-wake model can take as inputs the average values over the loop. Moreover, the model can be used in design and optimization studies. intensity of = ± 0 respectively. The velocity at the wing center induced by the two vortex ring cascades is The summation can be solved numerically for different values of λ 0 and η v , and its solution can be fitted as a function of these two parameters. The approximated solution of the summation is ∞ k=1 ϒ r,k (η 0 = −η v , λ 0 ) + ∞ k=1 ϒ r,k (η 0 = η v , λ 0 ) ≈ π 12 η π/2 v λ 1.1 0 . (A7) Figure A1 shows the comparison between the solution obtained numerically and the approximation given in Eq. (A7).
Appendix B: Comparison of glide ratios Figure B1. Glide ratio as a function of the lift coefficient for a case with AR = 20, κ 0 = 0.15 and C D0 = 0.1. Figure B2. Glide ratio as a function of the lift coefficient for a case with AR = 20, κ 0 = 0.2 and C D0 = 0.05. Figure B3. Glide ratio as a function of the lift coefficient for a case with AR = 20, κ 0 = 0.2 and C D0 = 0.1. Figure B4. Glide ratio as a function of the lift coefficient for a case with AR = 10, κ 0 = 0.15 and C D0 = 0.05. Figure B5. Glide ratio as a function of the lift coefficient for a case with AR = 10, κ 0 = 0.15 and C D0 = 0.1. Figure B6. Glide ratio as a function of the lift coefficient for a case with AR = 10, κ 0 = 0.2 and C D0 = 0.05. Figure B7. Glide ratio as a function of the lift coefficient for a case with AR = 10, κ 0 = 0.2 and C D0 = 0.1.

Appendix C: Nomenclature
Latin symbols A Wing area a r Radial induction a z Axial induction AR Wing aspect ratio b Wing span C D System drag coefficient C D0 System drag coefficient at zero lift C Di Induced drag coefficient C L Wing lift coefficient G C L /C D : glide ratio Latin symbols G 0 C L /C D0 : glide ratio related to the drag coefficient at zero lift G f C L /C f Di : glide ratio related to the farwake drag coefficient G n C L /C Intensity of the rolled-up trailed vortices κ 0 b/(2R 0 ): inverse turning ratio λ u 0 /v r : wing speed ratio λ 0 Normalized torsional parameter of the helicoidal wake φ j Angular position of the evaluation point p j around e 1 starting from e 2 ρ Air density θ Angular parameter of the helix θ j Angular position of the evaluation point p j around e 3 starting from e 2 ϒ n Shape factor modeling the near-wake axial induction ϒ f z,k Shape factor modeling the axial induction produced by the kth vortex ring ϒ r,k Shape factor modeling the radial induction produced by the kth vortex ring  (Trevisi et al., 2023) and opened with MATLAB or other open-source programming languages (e.g., Python, Octave).
Author contributions. FT conceptualized and developed the research methods, produced the results and wrote the draft version of the paper. AC and CEDR supervised the research, supported the methodology conceptualization and reviewed the paper.
Competing interests. At least one of the (co-)authors is a member of the editorial board of Wind Energy Science. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Review statement. This paper was edited by Roland Schmehl and reviewed by Emmanuel Branlard and one anonymous referee.