Low-and High-Fidelity Aerodynamic Simulations of Box Wing Kites for Airborne Wind Energy Applications

: High aerodynamic efﬁciency is a key design driver for airborne wind energy systems as it strongly affects the achievable energy output. Conventional ﬁxed-wing systems generally use aerofoils with a high thickness-to-chord ratio to achieve high efﬁciency and wing loading. The box wing concept suits thinner aerofoils as the load distribution can be changed with a lower wing span and structural reinforcements between the upper and lower wings. This paper presents an open-source toolchain for reliable aerodynamic simulations of parameterized box wing conﬁgurations, automating the design, meshing, and simulation setup processes. The aerodynamic tools include the steady 3D panel method solver APAME and the CFD-solver OpenFOAM, which use a steady Reynolds-Averaged Navier–Stokes approach with k-ω SST turbulence model. The ﬁnite-volume mesh for the CFD-solver is generated automatically with Pointwise using eight physical design parameters, ﬁve aerofoil proﬁles and mesh reﬁnement speciﬁcations. The panel method provided accurate and fast results in the linear lift region. For higher angles of attack, CFD simulations with high-to medium-quality meshes were required to obtain good agreement with measured lift and drag coefﬁcients. The CFD simulations showed that the upper wing stall lagged behind the lower wing, increasing the stall angle of attack compared to conventional ﬁxed-wing kites. In addition, the wing tip boundary layer separation was delayed compared to the wing root for the straight rectangular box wing. Choosing the design point and operational envelope wisely can enhance the aerodynamic performance of airborne wind energy kites, which are generally operated at a large angle of attack to maximise the wing loading and tether force, and through that, the power output of the system.


Introduction
In a recent outlook on world energy and carbon emissions, the International Energy Agency (IEA) stressed that offshore and onshore wind capacity is not only expected but also needs to expand substantially [1]. However, the agency also conceded that the limited availability of areas suitable for wind energy extraction would restrict this growth. Airborne wind energy (AWE) can reach higher wind resources using kites, needs fewer construction materials, and can potentially produce power levels similar to conventional wind turbines [2]. Kites for airborne wind energy systems (AWES) can be subdivided into three categories: fixed-wing kites, soft-wing kites, and hybrid kites [3]. The present paper focuses on a fixed-wing kite design.
Several concepts with fixed-wing kites have been pursued so far. The most common is the monoplane configuration. Prominent examples are the Makani M600, designed for a nominal power of 600 kW [4], and the Ampyx Power AP3, designed for a power of 150 kW [5]. A second configuration used for AWE is the box wing. The concept is similar to the biplane, connecting the tips of the two stacked wings with an aerofoil to form a closed loop. While the box wing was widely adopted for meteorological and man-lifting kites at the beginning of aviation due to its inherent structural stability and load-bearing  [15]. Both the biplane and box wing have a height-to-span ratio of 0.2. Adapted from Gall and Smith [10].
There have been several analyses on box and tandem wing designs in the last decade, mostly for transport aircraft. However, the aerodynamic analyses are nonetheless relevant to this study. In an effort to minimise induced drag, Khan [16] conducted a parametric study utilising vortex lattice algorithms and Euler inviscid simulations. The study revealed that the height-to-span ratio was the most critical design parameter, with an increase in this value resulting in a significant reduction of induced drag. Gagnon and Zingg [17] also utilised inviscid simulations to vary parameters, such as twist, sweep, and centre of gravity, for a transonic box wing regional jet configuration. After analysis, it was determined that the box wing configuration exhibits superior performance in terms of pressure drag and stability compared to a traditional wing design. However, no structural analysis or viscous drag study was conducted in their paper. Using a vortex panel method with corrections based on a designated set of drag polars, Andrews and Perez [18] conducted a parametric study to assess the performance of conventional and box wing aircraft. By altering parameters, including wing area, height-to-span ratio, and stagger, they could determine that the box wing design has the potential to exhibit increased aerodynamic efficiency compared to the conventional design. Finally, Bauer et al. [19] conducted an analysis using computational fluid dynamics (CFD) with a Reynolds-Averaged Navier-Stokes (RANS) solver and the k-ω turbulence model to evaluate the performance of a small-scale monoplane and biplane kite equipped with wind turbines. They adjusted geometrical parameters, such as aspect ratio, area, and aerofoils, and found that the biplane had a significantly higher lift coefficient than the monoplane.
Developing an automated process for generating a finite-volume CFD mesh from a parameterised box wing geometry is the main aim of this work, which is a follow-up of the graduation project [20] of the second author of the present paper. Generally, this step is the most time-consuming aspect of a CFD analysis. A fully parametrised geometry allows one to easily change the geometry or use an optimiser to obtain optimal design configurations. A second objective is to integrate this automatic mesh generation process into a reliable simulation framework for the prediction of the aerodynamic performance of the wing. Reducing the engineering time required for the aerodynamic design process is a major focus of this work. This challenge applies not only to the AWE sector but also to various other disciplines within and beyond the aerospace industry. Time deduction is realised by automating and identifying tasks that do not necessarily require user input, including selecting a favourable geometry, creating a mesh and pre-processing for CFD simulations, supervising the simulation, and post-processing. The framework is available in open source [21].
Regarding AWES, the automated framework is verified and employed to examine the aerodynamic efficiency of recent concepts, including the box wing, which proposes aerodynamic benefits when contrasted with traditional fixed-wing designs. Another possible benefit could be in terms of sustainable material choices. Carbon-fibre-reinforced polymers (CFRP) are used extensively to stiffen fixed-wing kites to withstand the high loading present during the pumping cycles. Carbon is a synthetic fibre, which requires a large amount of energy to manufacture [22]. The box wing can have a shorter wing span, and the structure also has geometrical self-stiffening properties. These two characteristics can decrease the need for CFRP materials.
The remaining paper is structured as follows. In Section 2, the proposed methodology is outlined, including the validation of the computational approach. In Section 3, the results for the investigated box wing are presented and discussed in detail, and in Section 4, conclusions are drawn.

Methodology
An outline of the toolchain is given by the flowchart in Figure 2. Three main groups of processes were defined, which can operate independently of each other provided the required inputs are pre-defined (illustrated in Figure 2 by colour and number).

1.
Generating a geometry and solving using the 3D panel method (yellow).
Even though many tasks were automated in this toolchain, several tasks still required user interference. These tasks can be divided into three user categories:

2.
Mesh quality and simulation performance monitoring (for a detailed description of quality judgement, one may refer to [20]).

3.
Post-processing of CFD simulation results. The flowchart depicts the three main groups connected by two decision blocks of User Category 2. The first process group entails an interaction between MATLAB (a programming and numeric computing platform) and APAME (an inviscid 3D panel method algorithm). The details on how the panel method is implemented in APAME can be found in [23]. First, a parametrization of the geometry is performed. This allows for possible optimisations of the geometry based on the chosen parameters. The tool also enables manual geometry modifications, enabling the user to generate a new mesh quickly based on any given change. Based on the literature (see Section 1) and a parametric study on the influence of parameters on the lift coefficient, the eight most important physical parameters were identified, visible in Figure 3 and Table 1, and five different aerofoil sections were chosen. Aerofoils were defined at the root and tip of the upper and lower wing and for the winglet. This set of parameters translated to a wireframe geometry converted into nodes and geometrical panels for APAME. The (panel) mesh was built in MATLAB. The end result consists of both the reference geometry and its corresponding aerodynamic performance. It should be noted that the reference aerodynamic design referred to in the flowchart is actually the reference geometry that meets the user's specific requirements for aerodynamic performance. Once the user performance requirements were satisfied, the reference aerodynamic design was translated to the second main process.  For this work, the performance requirements of Process 1 were based on the equivalent wing criteria. Equivalent means a wing similar to another AWE kite. For example, an equivalent wing can be defined with respect to the MegAWES AWE reference kite [24] to compare the performance of a monoplane and box wing design in AWE applications. The following five conditions are proposed:

1.
Equal (or higher) tether force, which can be assumed to be proportional to the total lift.

2.
Identical aerofoil profiles are utilised in the upper and lower wings' lifting areas.

3.
The wingspan remains either equivalent or smaller.

4.
The reference surface is either the same or smaller.

5.
The mass is either the same or smaller.
Translating these criteria into an objective function and constraints can lead to an optimisation which finds the equivalent wing for a given monoplane design. Even though Process 1 is implemented in the framework, running the optimisation is not part of this work but will be performed in a future study.
The second process group entails an interaction between MATLAB and the commercial mesh generation software Pointwise [25]. In MATLAB, all the commands for generating the 3D CFD mesh were exported to the native Pointwise scripting language Glyph. Pointwise could then execute, resulting in a mesh which can be used by OpenFOAM® (a free and open-source tool for computational fluid dynamics) [26]. Producing a surface mesh of excellent quality is essential as it is employed in generating the volume mesh, using the boundary layer extrusion from said surface mesh. Therefore, the possible conflicting regions had to be identified. For this design, the four wingtips were labelled as hazardous due to the 90º bending. The inner regions of corners are especially critical because as cells move away from the surface, they have a tendency to compress rather than to expand. This can lead to boundary layer mesh collision issues. This problem was solved using geometric relations to set the last boundary layer size, given a certain bent region surface mesh size and a growth ratio, as the initial dimension of the straight wing segment's elements. By employing this technique, greater uniformity was achieved in the boundary layer's outer edge, resulting in a smoother transition to isotropic cells. Since the geometric relationships were formulated for a 2D planar bent, an empirical correction factor was chosen to accommodate for both 3D effects and aerofoil thickness. Taking into account these two characteristics in the layer propagation would have significantly complicated the forward propagation procedure. Appendix A.2. in [20] shows the mathematical relations together with graphical illustrations.
Pointwise software operates in a bottom-up fashion, signifying that the meshing process advances gradually. First, the given wireframe was translated into nodes, connected by curves, which form a surface when grouped. Next, the curves were discretized to create the 1D mesh (Pointwise: Nodes), while the surfaces were used to produce the 2D mesh (Pointwise: Domains). Subsequently, the far-field zone, as well as the locally refined areas in the 3D mesh, were specified. Finally, the 3D mesh (Pointwise: Block) was created by incorporating the 1D and 2D meshes along with the specific meshing parameters. Once the mesh satisfies the quality requirements specified by the user, the CFD mesh is ready to be simulated in the OpenFOAM environment. The reason for using a commercial meshing tool is that open-source meshing algorithms produce low-quality meshes. Building a volumetric mesh generation tool from scratch would be too time-consuming for this research project.
The third process group entails an interaction between MATLAB and OpenFOAM. The user must provide the desired flow conditions given a good CFD mesh. Once the mesh has been parallelised into the available cores, the OpenFOAM solver (v2006) can be started. To achieve this, the "decomposeParDict" configuration is employed. To minimise the number of CPU boundaries, the Scotch decomposition algorithm is utilised. It is preferred due to its automated nature that does not require any user inputs [27]. The simulations use the SimpleFOAM solver (Semi-Implicit Method for Pressure-Linked Equations) [28] using steadystate RANS equations [28,29] with a k-ω SST turbulence model. The specific OpenFOAM v2006 implementation of the turbulence model is based on [30]. Fully turbulent boundary layers are assumed as it is likely that the surface will be rough, tripping the flow close to the leading edge at this Reynolds number. Table 2 shows the solver control settings and Table 3 the initial and boundary condition values. There is no established standard for the farfield dimensions, but Athadkar and Desai [31] found that using a far-field which extended 10 chords upstream to 15 chords downstream produced satisfactory results for their 2D aerofoil simulation. To prevent the boundary conditions from affecting the flow over the wing, the far-field box is extended from 15 chords upstream to 40 chords downstream and 15 chords in all other directions. This domain ensures that the wing's flow remains unaffected by the surrounding environment while staying at a reasonable computation cost. A more detailed description of the OpenFOAM setup and parameters can be found in [20]. Once the simulation results are ready, they are post-processed using ParaView.

Results and Discussion
Before generating aerodynamic data related to any innovative box wing design, the required tools had to be validated. Both APAME and OpenFOAM were used to investigate the box wing performance. This validation indirectly validated both aerodynamic solvers' mesh resolution and quality. The validation data came from experimental results given by Gall and Smith [10]. Table 4 lists the geometrical parameters of the box wing examined, while Table 5 specifies the flight conditions. The examined box wing has a constant chord that equals the mean aerodynamic chord (MAC), and with the exception of the winglet, the aerofoil profile remains constant across all sections of the wing. Although Gall and Smith [10] did not provide the cant radius, the value of 15% of the winglet length was obtained from the research conducted by Thedens [32]. Furthermore, it was not explicitly stated whether the upper wing's stagger was positive (with the upper wing positioned backward) or negative (with the upper wing positioned forward). It should be noted that the typical definition of stagger is usually the opposite sign. To save time and avoid having to conduct a CFD mesh resolution analysis for every potential stagger configuration, a stagger value of 0 was assumed for the kite in the volumetric mesh of OpenFOAM. Table 4. Geometrical parameters of the box wing for validation [20].  Figure 4 shows the lift coefficient versus the angle of attack given by APAME, for a different total number of panels. As the panel method is limited to inviscid solutions, only the linear region of the lift curve was considered. In addition, due to the numerical stability issues of APAME, the panels in the chordwise direction were forced to be kept constant at 59 panels. It can be seen that the lift coefficient remained nearly the same, regardless of the number of panels used in the discretisation. This might mean that the grid convergence happens at a lower resolution than the minimum value examined or that the design performance is insensitive to spanwise panel resolution. As the computational cost of the lowest resolution was deemed acceptable, it was decided to proceed with this number of panels in the reference design creation rather than exploring even lower resolutions.   Table 4 [20]. Figure 5 depicts the variation on the lift curve depending on the value and sign of the stagger with respect to the experimental data, given the lowest grid resolution. Three stagger values were evaluated, zero, negative, and positive mean aerodynamic chord. On average, the relative error with respect to the experimental results did not exceed 5% on any of the mesh resolutions and stagger configurations. Figure 6 shows a parametric study on both the height-to-span ratio and the stagger-to-chord ratio. Focusing on the staggerto-chord ratio suggested a more gentle lift slope can be obtained from a decreased stagger to negative infinity. APAME demonstrated satisfactory conformity with the experimental findings in the linear range of the lift coefficient, rendering it a suitable choice for this study.  Table 4 [20]. For the CFD analysis, the stagger remained at zero. Even though the lowest average relative error (3.5%) was achieved by a negative stagger (-MAC) in the previously mentioned panel method results, the errors were close enough (4.1%) to leave the stagger in the CFD kite design at zero. A modification would have required another mesh refinement study on the CFD mesh as described below, which is computationally expensive. To evaluate the accuracy of the solution from OpenFOAM, numerical results for aerodynamic coefficients (lift and drag) were analysed at three distinct resolutions. Strictly speaking, a grid convergence analysis is performed on one region, where the grid is resolved in three resolutions. In this work, it was chosen to have three different regions where grid refinements were applied. The computer specifications and runtime of both CFD meshing and simulations are found in Tables 6 and 7. A complete grid convergence analysis was performed using Richardson extrapolation [33]. The primary objective of this technique was to obtain solutions for varying grid levels and extrapolate them to ∆x = 0, which would correspond to the exact values of the coefficients under investigation. Three different grid resolutions were selected, as illustrated in Figure 7: To validate the obtained results, the CFD simulations were all compared against the experimental results. The three meshes were generated using Pointwise using the cell types shown in Figure 8 and the refinements mentioned before in Figure 7. A more detailed explanation of grid parameters can be found in [20]. Figure 9 illustrates the deviation of the CFD simulations from the experimental findings. While the lift coefficient trend agreed well with the experimental results, the maximum lift coefficient (C L max ) was over-predicted.

Parameter
Regarding the drag coefficient, the CFD results corresponded to the experimental values at low angles of attack (α ≤ 5), but at higher angles, the CFD results were lower. Two factors, the stagger and the cant radius, could be accountable for this discrepancy. The CFD simulations utilised zero stagger because of project-associated computational power availability and unknown sign, while the experimental box wing had non-zero stagger. This variation primarily affects the drag coefficient, but at high angles of attack, it could also impact the lift coefficient. Stagger affects the angle of incidence for each wing, and its impact is based on the level of flow interaction, which is subject to vertical separation. This vertical separation is mainly responsible for the efficiency factor having an effect on lift-induced drag [34]. Additionally, it modifies the wing area exposed to flow at non-zero angles of attack, leading to alterations in the total form drag. Another factor that may have contributed to the differences between CFD simulations and experiments is the cant radius. The experiment description did not specify the radius, and the value of 15% of the winglet could have led to variations in lift and drag coefficients.  The Richardson extrapolated value showed relative errors to the experiment in C L , which remained below 7% in most cases, except for C L max . Nevertheless, notable variations in the drag coefficient occurred at α = 12.3 • and 14.5 • . To predict lift accurately, medium-or high-resolution meshes were satisfactory, but for post-stall analysis, a high-resolution mesh would be preferable. The subsequent findings originate from the high-resolution mesh. Figures 10 and 11 show two 3D illustrations of the flow field around the box wing design at a 21 • angle of attack. This is a post-stall condition, where highly separated flow regions can be expected. Both figures show how the flow at the root section was fully separated on both wings, generating large regions with strong vorticity and reversed flow components. These visualisations can help better understand the phenomena occurring at high angles of attack.  For aerodynamic analyses, three points were selected on the lift coefficient curve based on their significance. The first point is located at α = 12.3 • , which is the end of the linear region, where the lift performance bends from the linear curve towards a maximum lift coefficient because of the separation of the boundary layer. The second point is positioned at α = 18.5 • and represents the location on the curve where the maximum lift coefficient is observed. This is followed by a sharp decline in lift caused by the pronounced separation of the boundary layer. The final point, located at α = 21 • , corresponds to the wing in post-stall, where significant non-linear effects are anticipated.
The velocity and pressure coefficient contours of the three angles of attack are then shown in more detail for three different sections along the wingspan in Figures 12-14, namely the root (symmetry plane), mid (middle in between root and tip), and tip section (the region where the wing-winglet connection begins). Figure 12 shows the velocity and pressure coefficient contour plots for the root section. At lower angles of attack, the flow remains attached, with high suction peaks close to the wing's leading edge, while the pressure coefficient approaches 0 towards the trailing edge. At the point of maximum lift, boundary layer separation is evident in the form of a substantial separation bubble that stretches from the upper wing leading edge to the trailing edge. Additionally, the root section of the lower wing experiences severe separation, resulting in a mismatch between the pressure at the trailing edge and the free stream pressure. As the lift coefficient has not dropped yet, one can assume the root stall did not greatly influence the lift yet. The wing's drag increased significantly, which can be related to these separated flows. During post-stall, the flow does not reattach on the upper wing.   [20]. Figure 13 shows the velocity and pressure coefficient contour plots for the midsection. Similar behaviour can be observed for the root section; however, the stall seems delayed. The upper wing did not stall until after the maximum lift, as shown in Figure 13c [20]. Figure 14 shows the velocity and pressure coefficient contour plots for the tip section. Generally, a smaller suction peak than the other wing sections was present. If Figure 14b is closely examined, localised separation and reattachment at the upper wing top side were already observed at lower angles of attack. Increasing the angle of attack to the stall angle or post-stall seemed to enlarge the separation bubble on the leading edge of the upper wing and increased the boundary layer separation of the lower wing. The upper wing seemed to have attached flow beyond the wing stall angle. The pressure between the free-stream and the upper-wing aerofoil was not matched, and the wake deflected upwards.  Overall, one can observe that the root stalled earlier than the tip. This is a common characteristic of straight (rectangular) wings as opposed to swept-back wings, where the tip tends to stall before the root. The 3D effects induced by the winglet positively affect the aerodynamic performance. A broad pattern concerning the relative separation of the upper and lower wings can be identified, where the separation of the upper wing lags behind that of the lower wing. This phenomenon is predictable, as when separation begins on the lower wing, the top surface flow is less forcefully redirected downwards than it would be under attached flow conditions. Therefore, the upper wing does not have to accelerate the flow as much to satisfy the Kutta condition. In principle, the upper wing will experience a lower effective angle of attack than the lower wing, reducing the lift force and delaying stall.
Zooming out to the whole box wing, the stall mechanism that seemed to be present is classified as thin-airfoil stall (as defined by Gault [36]). However, the areas in the vicinity of the upper wing's tip exhibit an anomalous stall classification, characterised by the occurrence of leading-edge stall and trailing-edge stall.
At various positions downstream of the box wing, the vorticity in the normal plane is illustrated in Figure 15 for maximum lift. In general, the vortex distribution exhibited a comparable pattern to conventional wings equipped with winglets, shedding counterclockwise vortices at the right wingtip and clockwise vortices at the left wingtip. However, the winglet on the box wing stretches to the upper wing, resulting in the fusion of the tip vortices into one that spans from the bottom to the upper wing. Strong vortices in the opposite direction of the wing tips on the same side can be observed near the box wing's root section due to the separated flow. The vortices' magnitude is proportionate to the separated regions, but their intensity is inferior to the wing structure's natural vortices. Consequently, the natural vortices dominated the vortex direction during the recombination, causing the large vortex regions to dissipate rapidly. Additionally, small regions of opposite rotating vortices were present near the wingtip of the upper wing at y = b/2 − R. The geometrical shape of the upper wing top surface explains these vortices as the flow has difficulties reattaching to the concave shape.

Conclusions
The main goal of this research was to develop an automated system that enables aerodynamic analysis of box wing configurations with both low and high levels of fidelity with the application to airborne wind energy (AWE) in particular. The framework can be used for two main purposes. The first is to evaluate several box wing designs to find ones equivalent to a conventional fixed-wing monoplane kite through low-fidelity optimisations. The second is to evaluate a given box wing design's performance based on geometrical parameters using CFD and an automatically generated volumetric mesh. The latter is described in detail in this paper.
First, the box wing configuration was parameterized utilising eight different physical parameters in addition to five distinct aerofoils. As a result, a broad spectrum of box wing kites can be designed. Using an optimiser, one can automate the design process by specifying certain performance criteria. Several criteria are proposed to develop an equivalent wing to a monoplane AWE kite. However, finding an equivalent wing is part of future studies.
Second, the APAME and OpenFOAM mesh generation was automated. The initial aerodynamic analysis software, APAME, utilises a steady panel method and is considered a low-fidelity solver. A relatively low-resolution mesh proved suitable for quickly generating satisfying lift coefficients on the linear slope. The high-fidelity aerodynamic solver, Open-FOAM, uses steady-state RANS with a high-quality mesh. Satisfying lift and drag coefficients were achieved in all regions with mesh refinements in the proximity of the wing and the wake, up to 10 MAC downstream. Nonetheless, for studying post-stall, a high-resolution mesh incorporating a greater refinement region in the wake is highly recommended.
Finally, the CFD computations showcase significant aspects of the straight rectangular box wing aerodynamics, notably the upper wing's stall being delayed relative to the lower wing. This happened due to the decrease in the local angle of attack of the upper wing, caused by the separated flow on the top surface of the lower wing. As a result, the stall angle was increased with respect to a monoplane using similar aerofoils. The stall characteristics are important in AWE applications, as the kite will operate most of its time in high-lift conditions to increase power production.
This work's high degree of automation enables box wing designs to be applied in various scenarios for AWE systems in future studies. The emphasis on automation led to simplifying the aerodynamic analyses, only performing a CFD simulation at three angles of attack. To enhance the analysis of box wing designs and compare them with conventional AWE kites, future work should focus on designing a detailed equivalent wing, incorporating a comprehensive structural model that can accurately predict mass and resistance to aerodynamic loads. An equivalent wing would open the door to a new reference kite that other research groups can use to analyse and accelerate the research into box wing design for AWE applications.