Aerodynamic Performance and Wake Flow of Crosswind Kite Power Systems

This paper presents some results from a computational fluid dynamics (CFD) model of a multi-megawatt crosswind kite spinning on a circular path in a straight downwind configuration. The unsteady Reynolds averaged Navier-Stokes equations closed by the k−ω SST turbulence model are solved in the three-dimensional space using ANSYS Fluent. The flow behaviour is examined at the rotation plane, and the overall (or global) induction factor is obtained by getting the weighted average of induction factors on multiple annuli over the swept area. The wake flow behaviour is also discussed in some details using velocity and pressure contour plots. In addition to the CFD model, an analytical model for calculating the average flow velocity and radii of the annular wake downstream of the kite is developed. The model is formulated based on the widely-used Jensen’s model which was developed for conventional wind turbines, and thus has a simple form. Expressions for the dimensionless wake flow velocity and wake radii are obtained by assuming self-similarity of flow velocity and linear wake expansion. Comparisons are made between numerical results from the analytical model and those from the CFD simulation. The level of agreement was found to be reasonably good. Such computational and analytical models are indispensable for kite farm layout design and optimization, where aerodynamic interactions between kites should be considered.


Introduction
It is known that winds become generally stronger and steadier at higher altitudes, and this is why reaching higher altitudes for harnessing wind power is very appealing. Wind power density is proportional to wind speed cubed; thus, if the wind speed doubles (2×), the wind power density octuples (8×), assuming the airflow density to remain almost constant. In addition, steadier winds at higher altitudes enhance the capacity factor. Interestingly, all these advantages may be gained using airborne wind energy (AWE) technologies. AWE concerns accessing and harnessing high-altitude (i.e., from several hundreds of meters to several kilometers) wind energy using tethered flying kites or aerostatic airborne devices (e.g., balloons and auto-gyros). Among different types of AWE technologies, crosswind kite power systems (CKPSs) look favourable to most AWE developers. Unlike a static kite which is only subjected to the incoming wind, CKPSs exploit much higher apparent wind speeds, and consequently larger aerodynamic forces, by flying at high speed in the transverse direction with respect to the incoming wind [1]. In addition to gaining larger aerodynamic forces, the crosswind principle allows the kite to access wind power from a much larger effective capture area while flying closed-loop patterns.
Compare this to a static kite which only harvests from a region of the sky corresponding to its projected cross-section (and/or the total rotor area of the turbine(s) it carries); see Figure 1 which shows schematic drawings of ground-based (or lift mode or pumping mode) and on-board (or drag mode) power generation via CKPSs. For more details about various AWE technologies, particularly CKPSs, please refer to refs. [2][3][4][5]. Schematic drawings of (a,b) ground-based (or lift mode or pumping mode), and (c) on-board (or drag mode), power generation via CKPSs. Subfigure (a) shows the reel-out or power generation phase, and subfigure (b) shows the reel-in or power consumption phase (adapted from [6]).
Most studies on CKPSs neglect the wind-kite aerodynamic interactions. In other words, the flow retardation or "induction" effects caused by a flying kite are neglected based on the notion that the kite flies over a significantly large area compared to the kite planform area. However, several researchers have attempted to include such effects in the aerodynamic modelling of CKPSs. Some examples are the works of [6][7][8][9][10][11][12][13].
In addition to wind-kite aerodynamic interactions, kite-kite and wake-kite interactions for CKPSs have been, for the most part, overlooked. Unlike the extensive research that exists on the layout design and optimization of wind farms for conventional wind turbines, very few studies exist on the subject for CKPSs. The two most relevant studies were performed by Fagiano [14] and Faggiani and Schmehl [15] who considered solely the mechanical interference between kites, i.e., kites were spaced such that the tethers did not become entangled with each other. Fagiano [14] developed an algorithm for maximizing the average net power per square km of land, generated in a farm composed of 4-kite units. To avoid aerodynamic interference between kites, it was assumed that the polyhedra limiting the tether-kite system flight regions do not intersect and that the maximum flight elevation of the downwind kites is lower than the minimum elevation of the upwind ones. Also, in the work by Faggiani and Schmehl [15], wake interaction effects were assumed to be negligible for the following two reasons: (i) the kite operating in pumping mode (i.e., ground-based generation) sweeps a very large airspace, and (ii) the kites can be flown at different altitudes and maneuvered such that perfect alignment with the wind is prevented.
However, it is not known yet whether avoiding only the mechanical interference between kites and staggering them spatially, as discussed above, will be sufficient for optimal layout design of kite farms. When designing and optimizing a farm of conventional wind turbines, wake flow from individual turbines and the aerodynamic interference between them are the primary factors to be considered. Thus, examining the aerodynamic features of individual kites, such as their wake flow, as well as the aerodynamic interference between neighbouring kites seems essential. Such studies are rare in the literature. The most relevant studies are refs. [12,[16][17][18][19] where the wake flow developed downstream of a single CKPS or multiple CKPSs was studied computationally.
In the works presented in [16,18,19], an in-house large-eddy simulation (LES) framework was used to examine the wake flow of a single CKPS or multiple CKPSs. In their simulations, they modelled the kite using the actuator line/segment technique, where the aerodynamic forces acting on the kite were calculated and added to the momentum equations through a forcing term. They concluded that the wind-kite aerodynamic interactions (in other words, induction factor), even if limited, cannot be neglected for CKPSs. Moreover, based on their numerical results, they suggested that the velocity deficit in the wake of a kite system with the optimal power output is half of the deficit for a conventional wind turbine. On the other hand, in refs. [12,17] an unsteady Reynolds-averaged Navier-Stokes (URANS) flow solver was coupled to a two-equation turbulence model to investigate wake flow features of a CKPS. In their simulations, the geometry of the kite was modelled fully, and thus, aerodynamic interactions between the kite and flow were computed directly. Their results confirmed the wake flow development to large distances from the kite with significant streamwise flow velocity deficit for a multi-megawatt kite system. Some CFD results reported in the present paper have previously been presented at the AWEC (Airborne Wind Energy Conference) 2019 [12].
Today, several kinematic (also called engineering or explicit) wake models exist for conventional wind turbines; some early models are those presented in [20][21][22][23]; some recent models are those developed in [24,25]-for a comprehensive list, please refer to ref. [26]. Although some studies, such as that by Andersen et al. [27] cast doubts on the robustness of the early wake models, other studies like the one by Kaldellis et al. [26] show, on the contrary, that such models provide "acceptable representation of the wake behaviour". The early wake models are, in fact, very popular in the wind energy community, mostly due to their fairly simple forms; for more details on the analytical wake models, the reader is referred to refs. [26,[28][29][30]. Unlike conventional wind turbines, there has been very little effort to develop analytical wake models for CKPSs. Only recently, Kaufman-Martin et al. [31] developed a model based on the entrainment hypothesis for annular wakes [32]. In general, a system of three ordinary differential equations should be solved simultaneously to find the variation of the wake radii and flow velocity. They also obtained closed-form equations by assuming that the annular wake does not drift radially. A good agreement was found between their results and the CFD results presented in [16]. This paper presents some computational fluid dynamic (CFD) results from an ongoing research on the aerodynamic modelling of a multi-megawatt CKPS. These numerical results primarily intend to highlight some features of the wake flow of the CKPS in the near and far wake regions. For example, it will be shown that the initial expansion region within which the static pressure recovers and reaches the atmospheric pressure extends up to one gyration radius. Also, a maximum of 10% flow velocity deficit can be observed at a point inside the wake located at a distance of 10 times the gyration radius from the kite. These CFD results are particularly important considering the fact that experimental studies and field tests on the wake flow of crosswind kite systems are not publicly available. In addition to the CFD model, this paper also presents the development of a simple analytical wake model. This includes closed form equations for obtaining the wake radii and flow velocity. This model is developed based on the widely used Jensen's [20] model developed for circular wakes generated by conventional wind turbines. To verify the results from the analytical model, they are compared with those from the CFD model.
The rest of the paper is organized as follows: in Section 2, first, the theoretical framework used for the CFD simulation is briefly discussed. The set-up for the CFD simulation, including but not limited to geometric and meshing details, boundary conditions, solution type, and convergence criteria are described next. Then, some numerical results on the induction factor and the wake flow, such as velocity components at various streamwise locations are presented and discussed. In Section 3, the derivation of the analytical wake model is given, and in Section 4, the numerical results from the analytical wake model are compared with those from the CFD model. Finally, Section 5 provides a summary of the main contributions and findings of the present paper.

Theoretical Background
Turbulent flow is characterized by random fluctuations created by the presence of numerous eddies. The increase in turbulence has important effects on the fluid flow parameters, such as density, pressure, and velocity. In general, there are three different methods for turbulent flow calculations: (1) Reynolds-averaged Navier Stokes (RANS) equations 'closed' by a turbulence model, (2) large-eddy simulation (LES), and (3) direct numerical simulation (DNS). In the present paper, three-dimensional (3-D), incompressible, unsteady, Reynolds-averaged Navier-Stokes (URANS) equations are solved computationally [33,34]: where u i (i = 1, 2, 3) are the mean (or time-averaged) flow velocity components; ρ and p are the density and time-averaged pressure, respectively; ν is the kinematic viscosity; also, S ij = (1/2)(∂u i /∂x j + ∂u j /∂x i ) are the mean rate of strain tensor (alternatively, µS ij are the Newtonian or laminar stresses; µ being the dynamic viscosity), and τ ij = −ρu i u j are called Reynolds or turbulent stresses; u i being fluctuating velocity components; x i and t are physical coordinates and time, respectively; finally, the overbar indicates the mean value of a quantity, which has been dropped from u i and p, for convenience.
Here, Menter's [35] shear-stress transport (SST) turbulence model, commonly referred to as the k − ω SST turbulence model, is adopted. This means that two additional partial differential equations should be solved simultaneously with Equations (1) and (2): where the first equation accounts for the turbulent mixing energy, i.e., the k-equation, and the second governs the specific turbulent dissipation rate, i.e., the ω-equation [36].
In the above equations, ρ k and ρ ω are the turbulent Prandtl numbers for k and ω and are calculated using blending functions; µ t is the turbulent viscosity; G k is the generation of turbulent kinetic energy due to mean velocity gradients. Also, G ω represents the generation of ω, which is calculated as in the k − ω model of Wilcox [37] (or the standard k − ω model); Γ k and Γ ω represent the effective diffusivity of k and ω, respectively, and D ω is the cross-diffusion term. Lastly, S k and S ω are the defined source terms given by the user. For more detailed explanation on the terms, please refer to the user guide of ANSYS Fluent [36].
The SST model utilizes the modified form of the original k − ω model of Wilcox [37]-by accounting for the effect of the transport of the principal turbulent shear stress τ ij -in the inner region of the boundary layer, and it switches to the standard k − ε model away from the wall boundaries (i.e., the outer region) and in free shear flows. The k − ω SST turbulence model is particularly superior in the prediction of adverse pressure gradient flows often encountered over the blades of wind turbines. To date, many researchers have adopted the k − ω SST turbulence model for CFD simulations of the flow around wind turbines: refs. [38][39][40][41][42][43] for horizontal-axis wind turbines, and refs. [44,45] for vertical-axis wind turbines, just to name a few. Figure 2 shows the CFD solution domain. The distance between the inlet and the outlet, i.e., the length of the simulation domain, is 20 km, and the width and the height of the domain are identical and equal to 9 km. The kite geometry is fully modelled, which is essentially a rectangular wing of span b = 53.94 m, chord length c = 3.72 m (aspect ratio AR = 14.5), and Clark Y airfoil sections. The Clark Y airfoil section was chosen in this study due to its popularity in general aviation and in the field of aerodynamics. Clark Y's relatively high lift-to-drag ratio, gentle stall, and nearly flat bottom make it very attractive in terms of aerodynamics and construction.

CFD Simulation Set-Up
The inlet of the simulation domain is placed 4.5 km upstream of the rotor plane. The reasonably large size of the computational domain, especially the large distance between the rotor plane and the outlet, prevents possible interactions between the wake flow and boundaries and thus ensures the reliability of the flow simulation results. Small downstream distances sometimes lead to solution divergence and sometimes to an overestimate of the power output coefficient. The flow blockage ratio in the solution domain, i.e., the solid area divided by the flow area, is calculated as 0.052%, that is negligible. A uniform spatio-temporal flow velocity of U ∞ = 8.33 m/s with a turbulence intensity of 1% is imposed at the inlet-the velocity inlet condition. The turbulence intensity is commonly defined as the ratio between the root mean square of turbulent velocity fluctuations to the mean flow velocity [33]. The rest of the boundary conditions are: zero gauge pressure at the outlet (i.e., pressure outlet), symmetry over the top, bottom and sides of the domain, and no-slip over all the faces of the kite; see Figure 2a). The SIMPLE algorithm is used for pressure-velocity coupling. Second order upwind scheme is utilized for spatial discretization of the momentum, turbulent kinetic energy, and specific dissipation rate equations. First order implicit scheme is used for time discretization of the transient term.
The nominal power output of the kite is approximately 5.5 MW, which is calculated based on Loyd's [1] equations for a lift mode CKPS flying at the average altitude of 355 m where the freestream density is ρ ∞ = 1.1752 kg/m 3 . The kite is flying continuously on a circular path of radius R = 123.3 m with the angular velocity Ω = 0.738 rad/s which corresponds to a linear velocity V c = 91.03 m/s. Therefore, the kite sweeps an annular area with the inner diameter D i = 2R − b = 192.66 m and outer diameter D o = 2R + b = 300.54 m. It is assumed that the actual kite is tethered and is subjected to a 12.5 m/s wind normal to the swept area. This corresponds to a wind speed in wind power class 7 (following the National Renewable Energy Laboratory's classification) at 10 m and assuming α = 1/7 in the power law velocity profile. The kite is reeling-out at one-third of the wind speed (to reach maximum power output according to Loyd's formulation), so that the relative velocity (i.e., v rel = (1 − 1/3) × 12.5 = 8.33 m/s) becomes exactly the same as the inlet flow velocity considered in the CFD simulations. It is further assumed that the kite would fly such that the aerodynamic efficiency, defined as χ = C L (C L /C D ) 2 [6], reaches the optimal value of χ = 1.23 × (1.23/0.1074) 2 = 161.3 and that there is no induction. The optimal value of χ may be calculated using tables for the lift and drag coefficients of Clark Y airfoil section versus angle of attack (with the Reynolds number of Re ≈ 20 × 10 6 ), which have been corrected for the effects of a finite aspect ratio; see ( [46], Chapter 5). It should be noted that C D = 0.1074 also includes the effective drag coefficient due to the tether [47]. ANSYS workbench meshing tool [48] is used for mesh generation, and the simulation domain is divided into stationary and rotating domains to allow for the implementation of the sliding mesh method (SMM). The SMM is a special case of dynamic mesh motion, where nodes move rigidly, provided that cell zones are connected through non-conformal interfaces. The kite is placed in the middle of the rotating domain, which rotates counterclockwise; see Figure 2b). Wake flow analysis requires a smooth, uniform mesh growth rate; therefore, the patch independent method is used to create the transitional mesh between the rotating and stationary domains. This method provides additional control over the maximum element size and leads to a smooth transition of the mesh. The maximum element size of the rotating domain is at the interface. A sphere of influence (of radius of 640 m) is created in the downstream region of the kite to effectively capture the physics of the wake flow; see Figure 2c,d. The advanced size function is also enabled with refinement options set to proximity and curvature -this helps to adequately capture and mesh small edges of an airfoil. The global and local mesh controls options are utilized to dictate the mesh refinements regions, faces and edges. In the solution domain, in total, there are 6.9 million cells and 3.8 million nodes. The maximum mesh skewness is 0.83.
The upper and lower sides of the airfoil are split in half to apply different mesh densities with more concentrations towards the leading and trailing edges. A bias factor of 2 is applied towards the leading edge while it is 10 towards the trailing edge. 22 inflation layers are created around the kite geometry to reliably simulate flow in the boundary layer and to accurately capture flow separation; see Figure 2e). The first layer height is defined by the y + requirements (i.e., y + ∼ 1) with a steady growth rate of 1.2. The airfoil mesh is swept across the kite surface to ensure better mesh quality and uniform inflation layers. The total number of hexahederal elements that are swept across the kite surface is around 3 millions.
The CFD simulation is performed using ANSYS Fluent, which employs a finite volume method for solving the Navier-Stokes equations, on Compute Canada clusters (accessed 12 cores and 192 GB of RAM). For the present study, the transient flow analysis is adopted. As mentioned earlier, the sliding mesh method is employed to model the transient flow effects between the stationary and rotating domains. The rotational velocity of the rotating domain is set at the angular velocity of the kite. The solution time step was considered to be equivalent to ∆θ = 1 deg of rotation (i.e., ∆t = ∆θ/Ω = 0.0236 s, CFL number = 0.43). This time step was judged to be sufficiently small considering best practices adopted for CFD simulation of similar systems, such as H-type Darrieus [44], and Savonius [49] wind turbines. The simulation was ran up to 51 cycles (each cycle took almost 24 h CPU time) to ensure the wake flow is fully developed up to large distances (i.e., 12R) from the rotor. Two criteria were established to ensure the solution is converged and reached post-transient conditions. The first criterion is that the residual values should fall below 10 −5 . The second criterion is met by monitoring the average value of the axial flow velocity over multiple concentric annuli across the wake width at several points from the rotor (e.g., at x = 5R and x = 8R). The solution is assumed to have converged when for any two successive revolutions of the kite, the average difference between the velocities at equivalent time steps within the period is less than 1%.
The methodology used in the CFD simulation of the kite, which was discussed above, is verified by considering computational results from two other simulations performed for an airfoil (Appendix A) and for a wing (Appendix B). This may be an acceptable verification approach given the lack of experimental data on the aerodynamics of crosswind kites. Although some flow features in the kite simulation, such as wake rotation, is absent from the airfoil/wing simulations, still, these two simulations and the kite simulation share some common flow features, such as boundary layer development, tip vortices generation, and possible flow separation. A similar simulation setup and mesh setting as those mentioned above have been used in the airfoil/wing simulations, and, wherever possible, a similar geometry and flow Reynolds number have been adopted; please see Appendices A and B for more details.

Induction Factor
According to the axial momentum theory [39,50,51], flow undergoes a velocity reduction and a simultaneous pressure increase when approaching an energy extracting device, such as a wind turbine or a crosswind kite. The magnitude of flow velocity reduction at the turbine/kite is commonly expressed as aU ∞ . In general, the axial (or streamwise) induction factor a may be defined as the ratio between the axial flow velocity induced by an energy extracting/generating device and the freestream velocity, U ∞ . In other words, a may be defined as a = 1 − (U/U ∞ ), where U is the flow velocity component normal to the plane of rotation at any point on that plane. The acceptable range of the mean value of a for a wind turbine or a crosswind kite is 0 < a < 1/2 while for an energy generating device, such as a propeller, a < 0. The induction factor can be linked to the power coefficient of the device, which for a wind turbine-like device can reach at best 16/27 ≈ 59.3% (known as the Betz-Joukowsky limit [52]) which occurs when a = 1/3. When the flow passes through an energy extracting device, the static pressure undergoes a sudden reduction while the flow velocity remains almost unchanged. Figure 3 shows the contour plots of the instantaneous axial induction factor at the kite location while the kite flies counterclockwise on a circular path. The kite is shown as a white rectangle in the figure. More precisely, the white region is the pressure side of the kite, which is subjected to the incoming wind flow. In the figure, positive values (i.e., warm colours) indicate flow retardation, while negative values (i.e., cool colours) indicate flow acceleration. Induction values higher than 1 may indicate a reverse flow which can occur close to the kite tips. At the tips, the axial velocity 'induced' by tip vortices is typically larger than that at in-board points. From Figure 3, flow retardation in the wake of the kite and acceleration in the regions close to the front and sides of the kite as well as in the inner (or core) region of the annulus swept by the kite are evident. The inner region may be defined as the area between the axis of rotation and the innermost tip of the kite. It is also evident that a greater flow retardation occurs close to the kite trailing edge, which gradually diminishes at larger distances from the kite on the circular path. This appears reasonable since during the time that the kite is advancing in a cycle, the regions already swept by the kite will be exposed to fresh undisturbed flow and thus be recovering from the effects of the presence of the kite.  Following the 'reduced axial velocity method' [53,54], which was originally developed for conventional wind turbines, the flow area is divided into multiple concentric rings (or annuli) over which the axial flow velocity is averaged. From the above equation for a, the average induction factor over each ring is obtained, where a r,i being the average induction over the ith ring; a r,i may be considered to be effectively acting at a radial distance corresponding to the center of the ring. Thus, one can obtain a radial distribution of the average induction factor, i.e., a r (r), which for the kite in the present study is shown in Figure 4. On the other hand, the global induction factor, a, is obtained by averaging the radial induction values: a = ∑ n i=1 w i a r,i , where w i = A r,i /A s is the ratio between the ith ring area and the total area swept by the kite. This global induction factor is an indicator of the overall aerodynamic performance of the kite. For the kite whose radial distribution of the average induction factor is shown in Figure 4, the global induction factor is calculated as a = 0.127 which is in good agreement (≈0.6-6% relative error) with the values obtained analytically [11,47] for the same kite. It should be noted that in ref. [11] the arithmetic mean of induction values over the rings covering the swept area was reported as the CFD-calculated induction factor. It is also interesting to note that a obtained computationally using a URANS flow solver closed by the k − ε turbulence model is only different by less than 2% from that obtained here by employing the k − ω SST turbulence model; refer to [11] for more details.  Figure 5 shows the distribution of the instantaneous static pressure coefficient, C p , and the dimensionless axial (or X-component) flow velocity, U/U ∞ , in the YZ-plane (6R × 6R in size) at different dimensionless axial distances X/R = 0.25, 0.50, 0.75, and 1, from the plane of rotation. The static pressure coefficient is defined as C p = (p − p ∞ )/(1/2ρ ∞ U 2 ∞ ), where p and p ∞ are the local and freestream static pressures, respectively. Typically, when the flow passes through an energy extracting device, the static pressure undergoes a sudden drop. As the wake flow moves downstream, it expands, where the static pressure increases and reaches the ambient pressure while the flow velocity inside the wake decreases. This marks the end of the initial expansion region, the length of which for conventional wind turbines is estimated to be about one diameter [28]. Nevertheless, some researchers, such as Frandsen et al. [23] and Bastankhah and Porte-Angel [24] argued that it is difficult to exactly determine the length of this region.

Wake Flow
As discussed in [28] and also confirmed by the plots shown in Figure 5, the regions immediately downstream the rotation plane are featured by non-uniform pressure and axial flow velocity deficit (i.e., C p < 0, U/U ∞ < 1). As seen, by increasing the axial distance from X/R = 0.25 to 1, C p vanishes increasingly over the flow area. It is evident that the static pressure gradually increases by X/R and eventually reaches the ambient static pressure at X/R = 1 or soon after. The static pressure increase is accompanied by a flow velocity reduction as may be confirmed from Figure 5e-h. Thus, one may take X = R as the end of the initial expansion region for the present crosswind kite. Nevertheless, it is too early to generalize this observation for other kite systems. Additional CFD simulations for kites of different scales, and more favourably, experimental studies are required to reach a general conclusion.   Figure 6 shows the distribution of the instantaneous, dimensionless axial flow velocity, U/U ∞ , in the YZ-plane at X/R = 0, 2, 4, 6, 8, and 10. These plots were obtained after the kite completed 50 cycles of rotation. At X/R = 0, as seen from Figure 6a, the flow velocity distribution in the YZ-plane is asymmetric, obviously due to the presence of the kite. Moreover, as also discussed in Section 2.3.1, flow acceleration (i.e., U/U ∞ > 1) can be observed in the front and sides of the kite. As seen from Figure 6b-f, as the wake moves further downstream, the axial flow velocity becomes nearly axisymmetric. Also, the annular wake flow becomes thicker as X/R is increased, which means that the wake expands towards both inner and outer undisturbed flows. Figures 7 and 8 show the distribution of the instantaneous dimensionless flow velocity in the vertical (or Y-axis) and lateral (or Z-axis) directions, represented by V/U ∞ and W/U ∞ , respectively. At X/R = 0, as seen from Figure 7a, in the inner flow region and for Y > 0, V/U ∞ has a predominantly negative sign while it is primarily positive for Y < 0, indicating that the flow is convergent in that region. On the other hand, in the outer region, the flow is generally divergent (or expanding), where the kite pushes the flow up in the positive half-plane and thus V/U ∞ > 0 while it pushes down the flow in the negative half-plane and thus V/U ∞ < 0. In Figure 7b,c, arrows (which represent the in-plane flow velocity vector) confirm the flow expansion but also show that the expansion is diminished as X/R is increased. Moreover, it is observed that the wake rotates in the clockwise direction as it moves downstream. This makes perfect sense since the kite spins in the counterclockwise direction. As a result, a flow asymmetry appears in the inner flow region and area swept by the kite, where V/U ∞ > 0 and V/U ∞ < 0 on most parts of the left and right semicircles (of the radius of approximately R), respectively. Similar observations can be made from plots shown in Figure 8, which are not repeated here for the sake of brevity.

A Simple Analytical Wake Model
In Section 2.4, the length of the initial wake expansion region for the present study crosswind kite system was obtained computationally to be about R. Nevertheless, for simplicity and inspired by the work of [20,21,23], it is assumed that the initial expansion behind the kite occurs immediately; in other words, the initial expansion region has a negligible length (Assumption 1). From the momentum theory, the axial flow velocity where the initial expansion region ends, is obtained from [39,50,51]: (7) Figure 9 shows schematically the annular wake (dotted line) downstream of a crosswind kite in the straight downwind configuration, where the kite flies essentially on a plane (referred to as rotation plane) normal to the incoming wind which is considered to be steady and uniformly distributed (Assumption 2). The fact that a crosswind kite sweeps an annulus in the sky permits considering an annular wake flow. It should be noted that this is the wake developed past the initial expansion region. Considering Assumption 1, one may then conclude that the wake essentially starts developing immediately after the rotation plane, as shown in Figure 9. It is assumed that the wake flow is entirely axial (i.e., with no rotational motion) and axisymmetric and that it expands linearly as a function of the axial coordinate X (Assumption 3). Motivated by the assumption made in [20,21] for wind turbines, it is further assumed that the wake has a starting width equal to the width of the annulus swept by the kite (Assumption 4); that is, D w,i (X = 0) = D i and D w,o (X = 0) = D o , where D w,i and D w,o are the inner and outer diameters of the wake, respectively; see Section 2.2 for the definition of D i and D o . The wake radius expansion is assumed to occur at two different rates towards the outer and inner flow regions, which are represented by κ o and κ i , respectively (Assumption 5). The rate of wake expansion, the so-called "entrainment" or "decay" constant, for on-shore conventional wind turbines is typically taken between 0.05 and 0.1 [20,23,24]. A commonly-used semi-empirical equation for estimating the entrainment constant for wind turbines is κ = 0.5/ ln(h/z 0 ), where h is the hub height, and z 0 is the surface roughness length [55]. the inlet AD is located where the initial expansion ends (marked by X = 0), which is immediately after the rotation plane (see Assumption 1). The outlet, on the other hand, is placed upstream of the critical point, i.e., X < X cr . The critical point (marked by X cr ) is defined as the point where the wake inner diameter vanishes. The side boundaries AB and CD are parallel to the freestream. At the inlet, as seen from the figure, the velocity of the flow going into CV1 along radial range 0 ≤ r < (D i /2) is uniform and of magnitude U i . The flow velocity along (D i /2) < r < (D o /2) is U 0 (see Equation (7)), and that is U ∞ for (D o /2) < r < (D w,o /2).
The flow velocity inside the wake at any point along the X-axis is considered uniform. This is equivalent to the so-called "top-hat" distribution of the velocity deficit in the wake, which is used for simplicity, in contrast to a more accurate but also more complex Gaussian or bell-shaped distribution (see, e.g., ref. [24] which considers the Gaussian velocity profile for conventional wind turbines). Moreover, the uniform distribution automatically satisfies "self-similarity" of wake flow velocity often assumed in the analytical wake models for conventional wind turbines. Thus, at the outlet BC, the velocity of the flow leaving CV1 along the radial range of (D w,i /2) < r < (D w,o /2) is uniform and is denoted by U w . Moreover, the flow velocity within the inner region 0 ≤ r < (D w,i /2) is assumed uniform and is denoted by U c . No inflow or outflow occurs across AB and CD since they are, in fact, streamlines which are parallel to the freestream.
axis of rotation/symmetry A B C D r 1 Figure 10. A schematic drawing showing the control volume CV1 (hatch pattern) whose boundaries (dashed line) are specified by ABCDA. CV1 encloses a small part of the freestream, the inner flow and the wake flow (dotted line) up to X < X cr . The inlet is specified by AD (X = 0), the outlet by BC (0 < X < X cr ), and the side boundaries by AB and CD. The thick line at X = 0 represents the area swept by the kite; also, r represents the radial distance from the axis of rotation/symmetry.
Considering the control volume CV1, the equation for the conservation of mass may be written as: Using the assumption that the wake expands linearly with the axial distance X from the rotation plane, we can write where 2κ o and 2κ i are the rates at which the wake diameter expands towards the freestream and the axis of rotation, respectively. From Equations (8) and (9), one can obtain which by defining Using the control volume CV2 shown in Figure 11, the dimensionless wake velocity for X > X cr (or alternatively ξ > ξ cr ) is obtained as and the dimensionless counterpart of Equation (9) can be written as axis of rotation/symmetry A B C D Figure 11. A schematic drawing showing the control volume CV2 (hatch pattern) whose boundaries (dashed line) are specified by ABCDA. CV2 enclose a part of the freestream, the inner flow and the wake flow (dotted line) up to X > X cr .
A simplifying assumption is to consider the flow velocity in the inner region to remain constant and equal to the freestream velocity, i.e., υ i = υ c = 1; a similar assumption was made in [31]. From CFD results presented in Section 2.4, it is seen that υ c may actually reach values above 1, which means that flow may accelerate inside the inner region; nevertheless, this increase of flow velocity is minimal and still the assumption of υ c = 1 may be reasonable. Moreover, from Equation (7), υ 0 = 1 − 2a; thus, Equations (11) and (12) may be simplified to and respectively.
Equations (14) and (15) provide simple mathematical equations to obtain the flow velocity in the wake as a function of ξ; moreover, with Equation (13) one can calculate the wake size as it expands downstream. Nevertheless, one should note the limitations of the model, the first being the unrealistic assumption of top-hat distribution of the velocity deficit and the second being the fact that the model solely considers the conservation of mass. In fact, one can mathematically show that the present model violates the conservation of momentum law in the wake, and the violation grows with ξ [32]. See Appendix C which provides some expressions to quantify the deviation from the conservation of momentum law.

Comparison between CFD and Analytical Results
We adopt the method proposed in [27] to obtain the wake radii and the average wake velocity from the CFD results presented in Section 2. The wake boundary (or width) in the radial direction at any streamwise coordinate ξ is considered as the point where the flow velocity reaches 99% of the freestream velocity-Kaufman-Martin et al. [31] considered the surface where the flow velocity reached the freestream velocity as the wake boundary. The average wake velocity was also obtained by integrating the area under the velocity profile between the wake inner and outer radii. Figure 12 shows the comparison between CFD results (circle markers) and those obtained via the analytical model for four different sets of entrainment constants κ i and show, respectively, the variation of υ w and δ w (δ w is used to refer to both δ w,o and δ w,i ) as a function of ξ. The first three analytical sets are basically different combinations of 0.1 and 0.05, which are typical entrainment constants used for conventional wind turbines onshore; see Section 3. The last set of κ i and κ o was found by fitting lines to the CFD wake inner and outer radii curves.
As seen from Figure 12a, the CFD wake velocity recovers quite fast within ξ ≤ 2 with υ w increasing above 0.9, while it almost plateaus at higher values of ξ. The analytical model with the given entrainment constants generally underestimates the value of υ w . The best agreement is achieved with κ i = κ o = 0.1 although the curves with other values of κ i and κ o may lead to a better agreement beyond ξ = 10 based on their trends. As seen from Figure 12b, according to CFD results, the wake expands almost linearly with ξ but with two different rates within ξ 2 and ξ 2 wake regions. In addition, the rate of wake expansion is quite different towards the inner and outer flow regions. Considering the quality of agreement over both υ w and δ w , particularly for larger values of ξ, all analytical results, except those obtained via κ i = κ o = 0.05, appear to reasonably match the CFD results. See Table 1 for a quantitative comparison between analytical and CFD results. In the table, the mean and standard deviation of the relative difference between the two sets of results for wake velocity and wake radii have been provided.  Figure 13 shows the wake velocity profile, more precisely, the radial distribution of the time-averaged wake velocity, obtained from the CFD simulation (solid line) and that obtained via the analytical model (dashed line) at different axial distances from the rotation plane for (a) κ i = κ o = 0.1, and (b) κ i = 0.091, κ o = 0.058. As seen, the CFD velocity profile looks like a Gaussian profile while the analytical one is a top-hat profile as was assumed in the beginning. From CFD profiles, it is evident that the flow velocity in the inner flow region is slightly higher than the freestream velocity (i.e., υ w > 1.0) for ξ ≤ 2 while it gradually decreases by increasing ξ but remains very close to the freestream velocity. Table 1. Comparison between analytical and CFD results, where the mean and standard deviation of the relative difference between the two sets of results are given.

Concluding Remarks
In this paper, using the CFD model, some aspects of the aerodynamics in the vicinity and downstream of a multi-megawatt CKPS were examined. URANS equations closed by the k − ω SST turbulence model were solved numerically inside a large sim-ulation domain. The radial distribution of the average axial induction factor as well as the global induction factor were obtained by following the reduced axial velocity method. Next, the variation of the static pressure and the axial flow velocity in small axial distances from the rotation plane were obtained. For the kite considered in this study, the initial expansion region was found to extend up to one gyration radius from the rotation plane, which is comparable to that of conventional wind turbines. The wake rotation caused an asymmetry in the vertical and lateral wake flow velocities while the axial flow velocity became increasingly axisymmetric as the axial distance from the rotation plane was increased. Although the flow velocity deficit in the wake reduced as the axial distance was increased, flow velocity deficits as large as 10% could be observed at large distances such as 10 times the gyration radius from the rotation plane. This means potentially a maximum of 27% reduction in the available wind power. CFD models can provide much details about the flow dynamics around kite power systems. However, it is not practical to use such models for, for example, conducting preliminary kite farm studies (where many simulations are normally required) due to their prohibitive CPU time.
In addition, in the present paper, the viability of developing a simple analytical wake model for CKPSs was proved. The present model, and those which will be developed subsequently, will pave the way for a better estimation of energy capture in kite farms and their layout design and optimization. These are crucial for estimating the levelized cost of energy for CKPSs, which directly relates to commercialization prospects of such systems. The model was built based on a similar model, known as Jensen's model in the literature, for conventional wind turbines. Comparisons between numerical results obtained from the proposed analytical model and those from the CFD simulation showed that the analytical model may be able to provide accurate information regarding the average wake flow velocity and its radii. The simple form of the analytical model permits fast computations, which in turn make the model ideal for conducting preliminary as well as optimization studies. Nevertheless, the level of accuracy of the analytical model appeared to be largely dependent on the magnitude of entrainment constants. Thus, future studies should highlight the correlation between the design parameters of a kite system and the entrainment constants. In the absence of publicly available experimental data, more computational studies, which preferably consider the wind gradient due to atmospheric boundary layer, would be necessary for finding a suitable range and/or closed-form equation for the entrainment constants to be used for CKPSs. Moreover, more rigorous analytical wake models such as those satisfying also the conservation of momentum and those considering a Gaussian-like velocity profile are desirable.

Acknowledgments:
The supports for this research project by New Leaf Management and Concordia University are acknowledged. The authors are also grateful to Compute Canada for providing computational resources.

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

Appendix A. 2-D Airfoil CFD Simulation
Here, we report on the simulations with a NACA 0012 airfoil section held in a flow of Re = 6 × 10 6 . The rationale for choosing NACA 0012 instead of Clark Y (which is the airfoil section of the kite) was the availability of experimental results for the NACA airfoil. The angle of attack was varied from α = −17.5 to 17.5 deg with a step of 1 deg for small to moderate values of α and with a step of 0.5 deg in the near-stall region. The lift and drag coefficients of the airfoil, represented by c l and c d , respectively, were obtained and compared with experimental results found in [56,57]. Figure A1 presents the numerical results on (a) the variation of c l as a function of α, and (b) the variation of c d as a function of α. As seen, the present CFD results and experimental ones are in very good agreement with each other, even in the stall region. Simulations were performed using a URANS flow solver (Ansys Fluent) closed by the k − ω SST (y + < 1) turbulence model. The spatial discretization was set to second order upwind scheme. However, for higher angles of attack, the solver was changed to densitybased implicit solver for the sake of accuracy and the flow discretization was set to second order upwind with pseudo transient option switched on. A C-shape domain of radius 100 m has been created upstream of the airfoil while the downstream domain is a rectangle of the length of 150 m and width of 200 m; see Figure A2. The boundary conditions are velocity inlet over the upstream and side boundaries, pressure outlet over the downstream boundary, and wall over the airfoil profile. Figure A2 also shows the structured mesh in different regions. The mesh quality parameters were set according to the criteria described in the user guide of ANSYS. The number of nodes and elements in the solution domain are 613,363 and 612,150, respectively, with the maximum 'skewness' of 0.135 (which is well below 0.9 that is the maximum limit recommended by ANSYS) and 'orthogonal quality' of 0.585 (which is well higher than 0.1 that is the minimum recommended value). There are 260 elements around the airfoil which is split into segments to impose higher mesh density around the leading and trailing edges.

Appendix B. 3-D Wing CFD Simulation
We present the wing simulation in this section. We adopt similar mesh settings as those used in the 2-D simulation (see Appendix A) for a rectangular wing planform with a Clark-Y airfoil section and aspect ratio of 6. Following experimental results presented in [58], the angle of attack was varied from α = 0.1 to 22.8 deg, and the Reynolds number was varied form Re = 3.1 × 10 6 at the zero-lift to Re = 2.59 × 10 6 at the maximum lift coefficient. The lift and drag coefficients of the wing, represented by C L and C D , respectively, were obtained and compared with the experimental results, as shown in Figure A3. The figure shows (a) the variation of C L as a function of α, and (b) the variation of C D as a function of α. As seen, the present CFD results and experimental ones are in generally good agreement. Simulations were performed using a RANS flow solver (Ansys Fluent) closed by the k − ω SST (y + < 1) turbulence model. The rationale for using the steady-state solver was to reduce the computational time. The spatial discretization used for momentum equation was second order upwind. As seen from Figure A4, the wing is placed in a cuboid solution domain. The distance between the upstream boundary and the leading edge is 60 m and that between the trailing edge and the downstream boundary is 150 m. The distance between the mid-span and the upper/lower wall is 60 m and that between the mid-span and the side walls is 90 m. The domain consists of 2 regions. There is an inner region which is close to the wing, where hexahedral elements are swept over the wing with 260 elements over the edge of the airfoil. This region is connected to the outer region through a non-conformal interface. The edge of the airfoil is split into segments to impose biased element sizes on the leading edge and trailing edge. Inflation layers are created around the airfoil to accurately simulate flow in the boundary layer and flow separation.
The first layer height is defined by the y + requirements with a growth ratio of 1.2. The total number of hexahedral elements swept across the wing surface is around 2.9 million; see Figure A4. The wake region is further refined to simulate separated flow from the wing surface. The boundary conditions are velocity inlet over the upstream and the sidewalls, pressure outlet over the downstream boundary, and wall over the wing surface. The minimum orthogonal quality is 0.193, and the maximum skewness for the mesh is 0.805.
It should be noted that at high values of α in the post-stall region where flow separation is significant, one should normally use a transient solution to ensure an accurate calculation of lift and drag coefficients; however, because of the adequate mesh refinement in the wake region behind the wing, the results from the present CFD study and experimental ones show fairly good agreement with each other as evidenced by the plots in Figure A3.

Appendix C. Deviation from the Conservation of Momentum
Considering CV2 ( Figure 11) and applying the conservation of mass, the difference between axial momentum out and in, which, in fact, represents the deviation from the conservation of momentum, may be written as Using Equations (7) and (15) and assuming υ i = 1, the expression for the momentum change is simplified to (A2) Note that since we are considering the wake at large distances from the rotation plane, Equation (15) has been used.
From Equation (A2), one can easily confirm that ∆M < 0, meaning that a theoretical axial momentum loss/deficit appears as a consequence of satisfying only the conservation of mass for the wake flow. Also, it can be shown that d|∆M|/dξ > 0, meaning that the deviation from the conservation of momentum increases as ξ is increased, and in the limit of infinitely large axial distances from the rotation plane, i.e., ξ → ∞, ∆M → −8a 2 (δ 2 o − δ 2 i ). To get a sense of the magnitude of the momentum deficit, we can normalize ∆M with respect to the thrust (or axial load) over the swept area [50]: Thus, one can write From Equation (A4), in the limit of very large ξ, ∆M/T → −2a/(1 − a). This means that if the energy harvesting device is lightly loaded, i.e., a 1, which is often the case for small-to medium-scale crosswind kites, then the deviation from the conservation of momentum is negligible and the proposed wake model may be reliable. On the other hand, for conventional wind turbines (or large-scale kite systems), which could have induction factors as high as a = 1/3, the deviation from the conservation of momentum could be as large as the thrust acting on the turbine (or kite), and the simple wake model may not be very reliable, particularly for large distances from the rotation plane.