Aero-servo simulations of an airborne wind energy system using geometry-resolved computational fluid dynamics

Airborne wind energy (AWE) is an innovative and promising technology for harnessing wind energy, often achieved through the use of tethered aircraft flying in crosswind patterns. A comprehensive understanding of the unsteady interactions between the wind and the aircraft is required for developing efficient, reliable, and safe AWE systems. High-fidelity simulation tools are essential for accurately predicting these interactions. To provide meaningful insights into crosswind flight maneuvers, they should incorporate the coupled nature of aerodynamics, dynamics, and control systems. Moreover, local aerodynamic phenomena, such as flow separation, play a significant role in the overall performance of the system and must be represented accurately. Capturing these phenomena requires resolving the complete geometry of the aircraft. Therefore, this work presents a geometry-resolved computational fluid dynamics framework of an AWE system, encompassing all lifting surfaces and integrating movable control surfaces, referred to as the virtual wind environment (VWE). Unlike existing models that only consider linear combinations of individual aerodynamic effects, the VWE addresses the challenge of combining the relevant aerodynamic interactions specific to crosswind flight motion. This VWE is coupled to the dynamics and control framework of an AWE system, enabling geometry-resolved aero-servo simulations. We demonstrate the novel simulation framework by tracking a pre-optimized one-loop power cycle in the VWE coupled to model predictive control, achieving 96 % of the reference power.

1 Introduction

Airborne wind energy (AWE) is an innovative and promising technology to harness wind energy and convert it into electricity, for example using tethered aircraft flying in crosswind patterns. A key advantage of AWE systems is their potential to operate at higher altitudes than conventional wind turbines, where the wind is stronger and more consistent (Diehl, 2014). Additionally, AWE systems require significantly less material for the same power generation, as the need for the tower and blade material near the axis of rotation is eliminated. Two main operation modes for energy conversion are currently pursued within emerging companies and academia: onboard generation (fly gen) and on-ground generation (ground gen). Furthermore, both soft kites and fixed wings are utilized as tethered aircraft (Cherubini et al., 2015; Vermillion et al., 2021). This work focuses on the modeling and simulation of fixed-wing ground-gen systems with a fixed ground station, which operate using a so-called pumping cycle. This cycle consists of a reel-out phase, during which the aircraft pulls on the tether and energy is extracted, and a reel-in phase, during which the tether is rewound, consuming a portion of the energy converted during the reel-out phase.

AWE systems experience a dynamic interaction between the aircraft, the atmosphere, and the controller. The presence of non-ideal wind conditions can induce unsteady aerodynamic phenomena, and its prediction remains an open challenge (Vermillion et al., 2021). Reliable control systems are crucial for steering the aircraft and operating the system safely. The fixed-wing system employs control surfaces to steer the aircraft, similar to conventional aircraft. Accurate prediction of aerodynamic forces and moments while deflecting the control surfaces and their impact on the system dynamics is crucial for designing safe control systems. Moreover, precisely tracking the intended flight path is important for the system's performance. The varying flow velocity encountered during crosswind flight maneuvers also influences local aerodynamic phenomena, such as flow separation. Because these phenomena significantly affect the overall performance of the system, they must be accurately represented in simulations, necessitating resolving the complete geometry of the aircraft in the simulations.

To study local aerodynamic phenomena, several studies have investigated various designs of fixed-wing AWE systems using computational fluid dynamics (CFD). Eijkelhof et al. (2023) developed an aerodynamic toolchain for the design analysis of AWE systems with a box-shaped wing using steady Reynolds-averaged Navier–Stokes (RANS) simulations. Vimalakanthan et al. (2018) conducted RANS simulations of a double-fuselage aircraft that included control surfaces. However, both studies focus on a steady horizontal flight, neglecting crosswind motion and setting the control surfaces in a fixed position. Kheiri et al. (2022) examined the wake flow of an AWE system using unsteady RANS for both the aerodynamics and the wind simulation, assuming circular flight motion and considering only the main wing, while omitting control surfaces. Castro-Fernández et al. (2021) incorporated more complex motion in their aerodynamic simulations by prescribing crosswind flight motion to a panel representation of the wing, without control surfaces, using the vortex lattice method (VLM). In the work of Fasel et al. (2019), a fully coupled aero-servo-elastic framework was used with a 3D panel method representation of the wing, focusing on optimization studies for morphing wings. Similarly, the work of Haas et al. (2022) and Crismer et al. (2024) presented an actuator line representation of the wing within a large-eddy simulation (LES) framework coupled to 3-degrees-of-freedom (DOF) and 6-DOF aircraft representations, respectively. However, in Haas et al. (2022), only the main wing is considered, and Crismer et al. (2024) uses an analytical model for control surface deflections. In both cases, the local aerodynamics are calculated from pre-computed steady airfoil polars, thus neglecting unsteady effects. In conclusion, none of the previously conducted studies on AWE have combined the aircraft's motion in geometry-resolved CFD simulations coupled to the system dynamics and control. Current state-of-the-art AWE simulations often rely on tracking pre-defined flight trajectories using modeling and optimal control tools, such as AWEbox (De Schutter et al., 2023). These tools typically employ aerodynamic models based on pre-computed coefficients or fast analytical approaches, such as the stability-derivative-based model proposed by Malz et al. (2019). While these aerodynamic models are computationally efficient, they exhibit notable limitations. Specifically, they fail to account for various aerodynamic effects, including unsteady phenomena and interactions, such as the influence of the rotational speed on the effectiveness of the ailerons, which can only be captured through flight testing or high-fidelity CFD. Such effects can lead to violations of the constraints defined during trajectory optimization, potentially resulting in degraded performance or even structural failure. Given the high cost of flight testing, accurate CFD tools are essential to predict and mitigate these effects during the design and planning stages.

This work presents a geometry-resolved CFD framework coupled to a controller, which is described in the next paragraph. The CFD framework integrates both the motion of the AWE system and the movement of control surfaces, including ailerons, elevators, and rudders. We refer to this comprehensive CFD framework as the virtual wind environment (VWE), and it uses the Chimera/overset technique, previously applied to simulate control surface deflections for general aircraft (Capsada and Heinrich, 2018). This technique offers flexibility by enabling complex grid configurations with multiple moving components through the decoupling of the background grid from the grids of the moving components. In this way, the VWE allows the combination of all individual aerodynamic contributions to be captured (such as the combination of the angle of attack, side-slip angle, angular rates, and control surface deflections), in contrast to analytical models that can only make predictions of individual contributions derived from a limited set of data points. Additionally, it offers improved accuracy in simulating localized flow phenomena, such as flow separation, which lower-fidelity techniques like the vortex-lattice method (VLM) cannot capture adequately.

The VWE is then coupled with the AWE system dynamics and the model predictive control (MPC) capability from AWEbox, enabling geometry-resolved aero-servo simulations. This coupling approach is illustrated in Fig. 1. The rigid-body motion and control surface deflection rates, computed by AWEbox, drive the motion of the aircraft component grids within the VWE. As these grids move, the VWE updates the flow field and determines the resulting forces and moments, which are then fed back into AWEbox. This coupling is demonstrated by tracking the pre-optimized one-loop power cycle in the VWE using an MPC-based controller. Additionally, the forces and moments derived from the VWE are compared to those predicted by the analytical aerodynamic model (AAM) embedded in AWEbox, highlighting the differences.

Figure 1 The aero-servo coupling approach.

The structure of this paper is outlined as follows. Section 2 introduces the reference aircraft studied. Section 3 describes the VWE. Section 4 covers the AWE system dynamics and control capabilities from AWEbox, explaining the trajectory optimization process and defining the controller used in this study. The coupling between the VWE and AWEbox is detailed through a stepwise approach in Sect. 5. Section 6 presents and discusses the results, and the conclusions are drawn in Sect. 7.

2 Reference aircraft

In this work, we consider a representation of an existing academic reference AWE aircraft, MegAWES (Eijkelhof and Schmehl, 2022). This aircraft has a wing area of 150 m², a wing span of 42.47 m, a root chord of 4.46 m, and a mass of 6885.2 kg. It is designed to have an electrical power output of up to 3 MW at 22 m s⁻¹ wind speed. This aircraft consists of a wing, two ailerons, an elevator (all-moving horizontal tail), two rudders (all-moving vertical tail), and two fuselages. In this work, the focus is on the lifting surfaces, so the fuselages are omitted from the aerodynamic models.

Note that the geometry used here differs slightly from that of Eijkelhof and Schmehl (2022), specifically regarding the aileron location, which extends from 62.0 % to 95.3 % of the half span. Additionally, the aileron gap is increased to 0.4 m to facilitate overset connectivity (as explained further below) with an allowable grid size for this simulation. There is also a rudder offset of 0.5 m from the elevator leading edge to prevent overlap between these components and to enable overset connectivity.

The center of gravity (CG) is located at $[-\mathrm{1.67},\mathrm{0},-\mathrm{0.229}]$ m in the geometry axis system, $[x,y,z{]}_{\text{G}}$, located at the leading edge of the main wing, as shown in Fig. 2.

Figure 2 (a) The MegAWES reference aircraft (Eijkelhof and Schmehl, 2022) and its representation in the VWE viewed from (b) the top and (c) the side. The values between the red brackets have been modified in this work, as explained in the text.

3 Virtual wind environment

A VWE is constructed using a geometry-resolved CFD framework in ANSYS Fluent, incorporating 6-DOF rigid-body motion and moving control surfaces. This makes the simulation suitable for analyzing the aerodynamics related to complex maneuvers and power cycles for airborne wind energy systems. This section first outlines the models and numerical settings employed. Subsequently, a detailed description of the aircraft grid, wind flow domain, and boundary conditions is provided. Finally, we explain the overset technique, which is used to connect the grids for various aircraft components and the background wind domain.

3.1 Flow model and numerical settings

The flow physics are modeled using the incompressible unsteady RANS equations with the k–ω SST model. Wall functions are used to model the boundary layer near the walls. Pressure–velocity coupling is achieved using a coupled scheme. Spatial and temporal discretization is implemented using a first-order upwind scheme for the convective terms in the momentum equations and a first-order implicit scheme with a time step of 5 ms, respectively.

3.2 Aircraft component grids

For each lifting surface component, an individual structured grid with C topology is constructed (Fig. 3). The grid domain of the main wing extends with a radius of 5 times the root chord in front of the wing, and the wake zone extends to 10 times the chord. The chord is divided into 132 cells, with refinements near the leading edge, near the trailing edge, and at the aileron location (at 75 % chord) to enable overset connectivity with the aileron. This overset connectivity allows for interpolation of the flow variables between the separate grids, as explained in Sect. 3.4. The wing domain is divided into 64 cells in the radial direction with a growth rate of 1.15. The size of the first grid cells from the wall corresponds to a y⁺ value of approximately 50, which is within the valid range for wall functions between 30 and 500. While this approach reduces the accuracy in prediction flow separation, it provides a reasonable overall impression of the flow around the aircraft and limits the computational expense of the simulation. In the spanwise direction, the wing is discretized into 202 divisions, with 70 divisions at the location of each aileron. The grid of the wing is based on the grid refinement study that was performed in Pynaert et al. (2023).

The C-grid radius and wake zone length of the elevator mesh are set to 1.75 and 3.5 times the elevator chord, respectively. For the rudder, these values are 1 and 2 times the rudder chord, respectively. The chords of both the rudder and the elevator are divided into 152 cells, with 36 cells in the radial direction and a growth rate of approximately 1.2. For these components, the size of the first grid cells from the wall yields a y⁺ value of around 100. In the spanwise direction, both the elevator and the rudder are discretized into 20 divisions. All aircraft component grids together consist of 6.3 million cells.

Figure 3 Cross section of the grids illustrating the C topology for (a) the wing, (b) the ailerons, (c) the elevator, and (d) the rudders (Pynaert et al., 2024).

3.3 Wind flow domain and boundary conditions

A rectangular grid is constructed to simulate a wind flow domain of $\mathrm{600}\times \mathrm{600}\times \mathrm{600}$ m with a uniform cell size of 3 m (Fig. 4a). This rectangular grid comprises a total of 8 million cells. This work presents a proof of concept for a single power cycle simulation without accommodating wake development, making this domain size sufficient for the current purpose. However, a larger background size would be necessary for conducting wake studies, as demonstrated in Haas et al. (2022) and Crismer et al. (2024). The decomposition between the background grid and the aircraft grid, combined with the overset method (see Sect. 3.4), facilitates the simulation of the 6-DOF motion of the aircraft within a large computational domain while ensuring adequate refinement near the wall for assessing the aircraft's local aerodynamics.

The following boundary conditions are applied to the wind flow and aircraft domain boundaries to simulate an AWE system operating within the atmospheric boundary layer (ABL), as shown in Fig. 4. To model the ABL, a logarithmic velocity profile is used, expressed as

$$\begin{array}{}\text{(1)}& v_{\text{w}}\left(z\right)={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{z+z_{\mathrm{0}}}{z_{\mathrm{0}}}}\right).\end{array}$$

In this equation, $u_{*}=\mathrm{0.3829}$ m s⁻¹ is the friction velocity, representing the reference wind velocity scale; κ=0.42 is the von Kármán constant; z denotes the height; and z₀=0.0002 m is the ground surface roughness height, characteristic of offshore conditions (Wieringa, 1992). The combination of u_*, z₀, and κ corresponds to a reference wind speed of u_ref=12 m s⁻¹ at the reference height of z_ref=100 m.

The specific dissipation rate (ω) profile is defined by Eq. (2), as proposed by Yang et al. (2009), to minimize the inconsistency between inlet profiles and the rough-wall formulation in the k–ω SST model:

$$\begin{array}{}\text{(2)}& \mathit{\omega }={\displaystyle \frac{u_{*}}{\mathit{\kappa }\sqrt{C_{\mathit{\mu }}}}}{\displaystyle \frac{\mathrm{1}}{z+z_{\mathrm{0}}}}.\end{array}$$

Here, C_μ=0.09 is a turbulence parameter. Although Yang et al. (2009) proposed a formulation for the turbulent kinetic energy (k) profile, it is not applied here due to limitations in defining it for the atmospheric boundary layer (ABL) conditions used. Instead, a default value of k=1 m² s⁻² is applied at the inlet. This value remains above k=0.7 m² s⁻² across the domain and reaches a maximum of k=1.75 m² s⁻² at ground level.

This logarithmic profile approximates the atmospheric surface layer and demonstrates the simulation's capability to incorporate specific wind profiles within the domain, which can be readily replaced with alternative profiles if needed. The entire domain is initialized using these inlet conditions. A uniform pressure of 1 atm is imposed at the outlet, and symmetry conditions are applied to the sides and top of the wind flow domain. The bottom boundary of the domain, representing the ground, is set as a stationary, no-slip wall with a roughness height z₀ to align with the imposed logarithmic wind velocity at the inlet. A moving no-slip wall condition is applied to the surfaces of the aircraft (wing, ailerons, elevator, and rudders).

Figure 4 Complete flow domain with boundary conditions at (a) the wind flow domain and (b) the aircraft component domains (Pynaert et al., 2024).

3.4 Overset technique

The background grid (wind flow domain) and the various aircraft component grids are coupled using overset boundary conditions at the boundaries of the aircraft component domains. This method enables the simulation of the aircraft's rigid-body motion, including deflected control surfaces, without deforming or re-generating the mesh. Figure 5 illustrates the connectivity between the background grid and the wing, elevator, and rudder grids in both the body-fixed x–z plane and the x–y plane.

In the overset technique, specific cell types (donor and receptor) are assigned to cells at the overset boundary. The flow solution from donor cells is interpolated and transferred to the receptor cells of other components, while no interpolation is used for the other cells. To assign these cell types, the grid priority method is applied, giving the highest priority to the control surfaces – aileron, rudder, and elevator – followed by the main wing, with the background grid having the lowest priority. For components of equal priority, a boundary distance-based priority method is used. This ensures that the overset cells are positioned as far as possible from moving boundaries, which promotes solver convergence.

Figure 5 Overset cell types in the wing component grid in the (a) body-fixed x–z plane at y=1.3 m and in the (c) body-fixed x–y plane at z=0 m. Overset cell types in the aileron grid in the (b) body-fixed x–z plane at y=15 m. Overset cell types in the elevator component grid in the (d) body-fixed x–y plane at z=1 m. Overset cell types in the rudder component grid in the (e) body-fixed x–y plane at z=1 m. Calculated cells are highlighted in green, donor cells in red, and receptor cells in blue (Pynaert et al., 2024).

4 AWE system dynamics and control

To simulate realistic trajectories aimed at maximizing power output, the AWEbox toolbox (De Schutter et al., 2023) is employed. AWEbox provides capabilities for the modeling and optimal control of both single- and multi-aircraft AWE systems and is built on CasADi, a nonlinear optimization framework using the algorithm differentiation tool Autodiff and inter-point optimizer IPOPT. This section outlines the specific capabilities from AWEbox that are utilized to build the aero-servo coupling. These capabilities include the formulation of the AWE system dynamics and an AAM of the aircraft. Additionally, the periodic optimal control problem (POCP) formulation is employed to generate a reference trajectory for the MegAWES aircraft. The final objective is to fly this trajectory within the VWE, utilizing the MPC toolbox for effective flight path tracking.

4.1 AWE system dynamics

This work considers the AWE system using 6-DOF aircraft dynamics and assumes a straight tether with mass and drag. The dynamics are represented using two reference frames: a body-fixed reference frame, $[x,y,z{]}_{\text{B}}$, located at the CG of the aircraft, and an inertial frame, $[x,y,z{]}_{\text{I}}$, positioned at the ground station. In the body-fixed frame, the x axis points to the rear of the aircraft, the z axis points upward, and the y axis extends towards the right wing, forming a right-handed coordinate system. In the inertial frame, the x axis aligns with the wind direction, the z axis points upward, and the y axis completes the right-handed system. These coordinate systems are illustrated in Fig. 6. This paper uses the following convention: lowercase italic letters represent scalars, bold lowercase letters denote vectors, and bold uppercase letters indicate matrices.

The state variables $\mathit{x}=\left(\begin{array}{cccccccc}\mathit{q}& \dot{\mathit{q}}& \mathbf{R}& \mathit{\omega }& \mathit{\delta }& l& \dot{l}& \ddot{l}\end{array}\right)$ of the system include the aircraft's position q and the velocity $\dot{\mathit{q}}$ in the inertial frame; the direct cosine matrix (DCM) R, representing the orientation of the aircraft; the angular velocity ω in the body-fixed frame; the aileron deflection δ_a; the rudder deflection δ_r; and the elevator deflection δ_e. The control surface deflections are grouped in δ. Additionally, the states include the tether length l, its reel-in/reel-out speed $\dot{l}$, and its acceleration $\ddot{l}$. The control inputs $\mathit{u}=\left(\dot{\mathit{\delta }}\textcolor{red}{}\right)$ to the system are the deflection rates of the control surfaces $\dot{{\mathit{\delta }}_{\text{a}}}$, $\dot{{\mathit{\delta }}_{\text{r}}}$, and $\dot{{\mathit{\delta }}_{\text{e}}}$, collected in $\dot{\mathit{\delta }}$, as well as the tether reel-in/reel-out jerk .

The relevant system parameters p include the wind speed v_w(z), defined by u_ref, z_ref, and z₀; the aircraft mass m_W; and the aircraft's inertia tensor J, defined in the body-reference frame. Figure 6 provides a visualization of the key system states, controls, and parameters. A comprehensive overview of the AWE system parameters and constraints can be found in Appendix B.

Figure 6 Visualization of the key states, controls, and parameters of the AWE system dynamics (Pynaert et al., 2024).

The system dynamics model employed in this study is based on the formulation presented in Gros and Diehl (2013), derived using Lagrangian mechanics. The resulting translation dynamics for a single aircraft is expressed as

$$\begin{array}{}\text{(3)}& (m_{\text{W}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{m}_{\text{T}})\ddot{\mathit{q}}+\mathit{\lambda }\mathit{q}={\mathit{f}}_{\text{e,I}}-(m_{\text{W}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_{\text{T}})g{\mathit{l}}_z.\end{array}$$

In this equation, ${\mathit{l}}_z={\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]}_{\text{I}}^{\mathrm{T}}$, and λ is the algebraic Lagrange multiplier associated with the constraint c. The total external force f_e,I, expressed in the inertial frame I, acting on the system comprises the aerodynamic force of the aircraft, f_I, and the tether drag force, f_T,I. The calculation of the tether drag is performed by dividing the tether into five segments, applying a multi-segment drag model as described in De Schutter et al. (2023). This model utilizes a constant tether drag coefficient, C_D,T, set to 1.2. The tether mass, m_T, is defined by

$$\begin{array}{}\text{(4)}& m_{\text{T}}={\mathit{\rho }}_{\text{T}}l{\displaystyle \frac{\mathit{\pi }{D}_{\text{T}}^{\mathrm{2}}}{\mathrm{4}}}.\end{array}$$

In this equation, D_T and ρ_T are the tether diameter and density, respectively, whose values are given in Appendix B.

The rotational dynamics is given by

$$\begin{array}{}\text{(5)}& \mathbf{J}\dot{\mathit{\omega }}={\mathit{m}}_{\text{e,B}}-\mathit{\omega }\times \mathbf{J}\mathit{\omega }.\end{array}$$

The total external moment m_e,B, expressed in the body-fixed frame B, is equal to the aircraft's aerodynamic moment, m_B, without contribution from the tether, as the tether is attached to the CG. In the next section, we discuss the aircraft's aerodynamic forces f_I and moments m_B in depth.

The aircraft is constrained to ensure that the distance between the aircraft's CG and the origin matches the tether length, enforcing a straight tether:

$$\begin{array}{}\text{(6)}& c={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{q}}^{\mathrm{T}}\mathit{q}-l^{\mathrm{2}})=\mathrm{0}.\end{array}$$

In Malz et al. (2019), it was found that the straight-tether assumption is adequate for estimating power generation in a small-scale airborne wind energy (AWE) system (specifically, the AP2 developed by the former Ampyx Power). In contrast, the study of Heydarnia et al. (2025), based on the MegAWES aircraft, concluded that the straight-tether assumption can lead to the overestimation of harvested power by up to 33 %. Future work will focus on incorporating tether sag into the system dynamics model. The DCM R contains the unit vectors of the body-fixed frame in the inertial frame. This non-minimal coordinate representation requires an orthonormality constraint:

$$\begin{array}{}\text{(7)}& c_{\text{R}}={\mathbf{P}}_{\text{ut}}({\mathbf{R}}^{\mathrm{T}}\mathbf{R}-\mathbf{I})=\mathrm{0}.\end{array}$$

In this equation, I is the identity matrix, and the operator P_ut is used to select the six upper triangular elements of a matrix (De Schutter et al., 2023).

Both the dynamics equations and the constraints c, c_R need to be enforced. An order reduction technique is applied to obtain an index-1 differential-algebraic equation by differentiating c twice with respect to time. Consistency conditions $(c,\dot{c},c_{\text{R}})=\mathrm{0}$ must be enforced at an arbitrary time point in the trajectory. The system's kinematics are integrated in time using an explicit Euler scheme, which is consistent with the time integration of mesh movement in ANSYS Fluent. For a more detailed explanation of the AWE system dynamics and kinematics in AWEbox, the reader is referred to De Schutter et al. (2023).

4.2 Analytical aerodynamic model

To complete the dynamic model in AWEbox, the AAM proposed in Malz et al. (2019) has been adapted for the MegAWES aircraft. This model is expressed by Eq. (9), where the superscript a refers to the AAM. Note that this model uses a different axis system, referred to as the aerodynamic axis system $[x,y,z{]}_{\text{A}}$, with the x axis pointing forward and the z axis pointing downward. In this equation, ρ(z) represents the air density and is modeled according to the international standard atmosphere (Archer, 2013), S is the wing surface area, and v_a is the apparent wind velocity, which is a function of both the wind velocity v_w(z) (as described in Eq. 1) and the aircraft velocity $\dot{\mathit{q}}$:

$$\begin{array}{}\text{(8)}& {\mathit{v}}_{\text{a}}=\left[v_{\text{w}}\right(z),\mathrm{0},\mathrm{0}{]}^{\mathrm{T}}-\dot{\mathit{q}}.\end{array}$$

The transformation matrix T transforms a vector in the aerodynamic frame to the body-fixed frame. The force coefficients C_x, C_y, and C_z are functions of the angle of attack α, sideslip angle β, roll rate p, pitch rate q, yaw rate r, and control surface deflections δ_a,e,r, as described by Eq. (11). The roll moment coefficient C_l, the pitch moment coefficient C_m, and the yaw moment coefficient C_n are computed similarly. The stability derivatives C_i,j represent the contributions of the quantity $j=\mathit{\{}\mathit{\alpha },\mathit{\beta },p,q,r,{\mathit{\delta }}_{\text{a,e,r}}\mathit{\}}$ to the forces in the i direction and the moments about the i axis. These derivatives are computed with the aid of the VWE using simple flight maneuvers (see Appendix A). The stability derivatives are then fitted to a second-order polynomial function of α. The polynomial coefficients, c₂, c₁, and c₀, are summarized in Tables A1 and A2.

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\mathit{f}}_{\text{I}}^{\text{a}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\left(z\right)\left|\right|{\mathit{v}}_{\text{a}}|{|}^{\mathrm{2}}S\mathbf{RT}\left[\begin{array}{c}{C}_x\\ C_y\\ C_z\end{array}\right],\\ & {\mathit{m}}_{\text{B}}^{\text{a}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }\left(z\right)\left|\right|{\mathit{v}}_{\text{a}}|{|}^{\mathrm{2}}S\mathbf{T}\left[\begin{array}{c}b{C}_l\\ c{C}_{\mathrm{m}}\\ b{C}_n\end{array}\right]\end{array}\text{(10)}& {\displaystyle }\mathbf{T}=\left[\begin{array}{ccc}-\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right]\text{(11)}& {\displaystyle }\begin{array}{rl}\left[\begin{array}{c}{C}_x\\ C_y\\ C_z\end{array}\right]& =\left[\begin{array}{c}{C}_{x,\mathrm{0}}\\ C_{y,\mathrm{0}}\\ C_{z,\mathrm{0}}\end{array}\right]+\left[\begin{array}{c}{C}_{x,\mathit{\beta }}\\ C_{y,\mathit{\beta }}\\ C_{z,\mathit{\beta }}\end{array}\right]\mathit{\beta }\\ & +\left[\begin{array}{ccc}{C}_{x,p}& C_{x,q}& C_{x,r}\\ C_{y,p}& C_{y,q}& C_{y,r}\\ C_{z,p}& C_{x,q}& C_{x,r}\end{array}\right]\left[\begin{array}{c}bp\\ cq\\ br\end{array}\right]{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\left|\right|{\mathit{v}}_{\text{a}}\left|\right|}}\\ & +\left[\begin{array}{c}{C}_{x,{\mathit{\delta }}_{\text{a}}}\\ C_{y,{\mathit{\delta }}_{\text{a}}}\\ C_{z,{\mathit{\delta }}_{\text{a}}}\end{array}\right]{\mathit{\delta }}_{\text{a}}+\left[\begin{array}{c}{C}_{x,{\mathit{\delta }}_{\text{e}}}\\ C_{y,{\mathit{\delta }}_{\text{e}}}\\ C_{z,{\mathit{\delta }}_{\text{e}}}\end{array}\right]{\mathit{\delta }}_{\text{e}}+\left[\begin{array}{c}{C}_{x,{\mathit{\delta }}_{\text{r}}}\\ C_{y,{\mathit{\delta }}_{\text{r}}}\\ C_{z,{\mathit{\delta }}_{\text{r}}}\end{array}\right]{\mathit{\delta }}_{\text{r}}\end{array}\end{array}$$

While this model includes the primary aerodynamic effects required to simulate 6-DOF aircraft maneuvers, it relies on the following assumptions. The stability derivatives are calculated using CFD simulations at a constant flight speed of 80 m s⁻¹, and consequently, the Reynolds number is assumed constant during the computation of these derivatives. However, both the flight speed and the Reynolds number vary during flight. Additionally, the model assumes a linear relationship between β, p, q, r, δ_a,e,r, and their effects on the force and moment coefficients. Furthermore, this model is quasi-steady and, therefore, independent of time, while in reality, unsteady aerodynamic effects occur. These limitations are addressed in the VWE, and the resulting differences are discussed in the Results section.

Figure 7 Illustration of the aerodynamic properties relevant to the AAM.

4.3 Flight path generation

The reference flight path (x^r(t), u^r(t)) considered in this work is generated by solving a POCP with a time period T, which is treated as an optimization variable. The POCP is formulated by Eqs. (13)–(16) (De Schutter et al., 2023). The objective function combines the power and penalties on the reference control actuation u^r(t) to prevent actuator fatigue, on the sideslip angle β(t) to avoid side forces, and on angular accelerations $\dot{\mathit{\omega }}\left(t\right)$ to prevent maneuvers that are too aggressive. These variables are collected in $\widehat{w}\left(t\right)$ and weighted by the matrix W. The optimization variables include the reference system states x^r(t), the reference control inputs u^r(t), the reference algebraic Lagrange multiplier λ^r(t), and the time period T. The power output of the system is determined by

$$\begin{array}{}\text{(12)}& P\left(t\right)=F_{\text{T}}\left(t\right)\dot{l}\left(t\right)=-\mathit{\lambda }\left(t\right)l\left(t\right)\dot{l}\left(t\right).\end{array}$$

Here, F_T represents the tension of the tether. The AWE system dynamics (including the AAM) and kinematics, represented by F, and the inequality constraints for path generation, represented by h_g, must be satisfied at every time step. Equation (15) bundles the following constraints. First, constraints are applied to ensure that the flight envelope (angle of attack and sideslip angle) is not violated. Furthermore, constraints ensure that the maximum tether force is not exceeded. Finally, aircraft orientation constraints prevent collision between the aircraft and the tether. Additionally, bounds are imposed on flight altitude, tether length, speed, acceleration, aircraft angular velocity, control surface deflections and their rates, and the time period T. Finally, the reference initial state x^r(0) must be equal to the reference final state x^r(T) to enforce the periodicity of the trajectory (Eq. 16).

$$\begin{array}{}\text{(13)}& {\displaystyle }\underset{{\mathit{x}}^{\text{r}}\left(t\right),{\mathit{u}}^{\text{r}}\left(t\right),{\mathit{\lambda }}^{\text{r}}\left(t\right),T}{min}{\displaystyle \frac{\mathrm{1}}{T}}\underset{\mathrm{0}}{\overset{T}{\int }}(-P(t)+\widehat{w}(t{)}^{\mathrm{T}}\mathbf{W}\widehat{w}\left(t\right))\mathrm{d}t\text{(14)}& {\displaystyle }\text{s.t.}\mathbf{F}\left({\dot{\mathit{x}}}^{\text{r}}\right(t),{\mathit{x}}^{\text{r}}(t),{\mathit{u}}^{\text{r}}(t),{\mathit{\lambda }}^{\text{r}}(t),\mathbf{p})=\mathrm{0},\forall t\in [\mathrm{0},T]\text{(15)}& {\displaystyle }{\mathit{h}}_{\text{g}}\left({\dot{\mathit{x}}}^{\text{r}}\right(t),{\mathit{x}}^{\text{r}}(t),{\mathit{u}}^{\text{r}}(t),{\mathit{\lambda }}^{\text{r}}(t),\mathbf{p})\le \mathrm{0},\forall t\in [\mathrm{0},T]\text{(16)}& {\displaystyle }{\mathit{x}}^{\text{r}}\left(\mathrm{0}\right)-{\mathit{x}}^{\text{r}}\left(T\right)=\mathrm{0}\end{array}$$

4.4 Flight path tracking

The final capability used from the AWEbox toolbox is the MPC, which is used to track the reference flight trajectory. This controller, called at the current time ${\widehat{t}}_{\mathrm{0}}$, solves an optimal control problem during the simulation to steer the aircraft toward the reference flight path in an optimal manner. The optimal control formulation, with a moving time horizon T_h, is given by Eqs. (17)–(20) (Gros et al., 2013). The objective is to minimize the difference between the system states x(t) and controls u(t) over the upcoming time horizon and the optimal reference states x^r(t) and controls u^r(t), which are determined using the method described in the previous section. The weighting matrices Q_c, R_c, and P_c are used to track the states, the controls, and the terminal cost, respectively. For this problem, equal weightings are assigned to each state and control variable. The optimization variables are the system states x(t), control inputs u(t), and the algebraic Lagrange multiplier λ(t). Similar to the optimal control problem for flight path generation, the flight dynamics F (including the AAM) and the path tracking constraints h_t must hold over the MPC time horizon. The tracking constraints h_t are more relaxed than the generation constraints h_g to allow for more controllability, and they are summarized in Appendix B. The initial state $\mathit{x}\left({\widehat{t}}_{\mathrm{0}}\right)$ is set equal to the current state estimate ${\widehat{\mathbf{x}}}_{\mathrm{0}}$.

$$\begin{array}{}\text{(17)}& {\displaystyle }\begin{array}{rl}\underset{\mathit{x}\left(t\right),\mathit{u}\left(t\right),\mathit{\lambda }\left(t\right)}{min}& \underset{{\widehat{t}}_{\mathrm{0}}}{\overset{{\widehat{t}}_{\mathrm{0}}+T_{\text{h}}}{\int }}\left(\right|\left|\mathit{x}\right(t)-{\mathit{x}}^{\text{r}}(t\left)\right|{|}_{Q_{\text{c}}}^{\mathrm{2}}+\left|\right|\mathit{u}\left(t\right)-{\mathit{u}}^{\text{r}}\left(t\right)\left|{|}_{R_{\text{c}}}^{\mathrm{2}}\right)\mathrm{d}t\\ & +\left|\right|\mathit{x}({\widehat{t}}_{\mathrm{0}}+T_{\text{h}})-{\mathit{x}}^{\text{r}}({\widehat{t}}_{\mathrm{0}}+T_{\text{h}})|{|}_{P_{\text{c}}}^{\mathrm{2}}\end{array}\text{(18)}& {\displaystyle }\mathrm{\text{s.t.}}\mathbf{F}\left(\dot{\mathit{x}}\right(t),\mathit{x}(t),\mathit{u}(t),\mathit{\lambda }(t),\mathbf{p})=\mathrm{0},\forall t\in [{\widehat{t}}_{\mathrm{0}},{\widehat{t}}_{\mathrm{0}}+T_{\text{h}}]\text{(19)}& {\displaystyle }{\mathit{h}}_{\text{t}}\left(\dot{\mathit{x}}\right(t),\mathit{x}(t),\mathit{u}(t),\mathit{\lambda }(t),\mathbf{p})\le \mathrm{0},\forall t\in [{\widehat{t}}_{\mathrm{0}},{\widehat{t}}_{\mathrm{0}}+T_{\text{h}}]\text{(20)}& {\displaystyle }\mathit{x}\left({\widehat{t}}_{\mathrm{0}}\right)={\widehat{\mathbf{x}}}_{\mathrm{0}}\end{array}$$

To solve both the flight path generation and the tracking optimal control problems, the problem formulations are transcribed to a nonlinear program (NLP) using direct collocation and solved using an interior-point homotopy (IPH) method. For a detailed description of the solution methods for the optimal control formulations, the reader is referred to De Schutter et al. (2023).

The MPC's sample time is 5 ms, matching the time step in the VWE, and the prediction horizon is T_h=0.1 s, corresponding to 20 sample times. In this case, the MPC is re-evaluated at each time step. The current settings were established empirically to ensure simulation stability.

5 Aero-servo coupling

The previous sections introduced the VWE based on geometry-resolved CFD and the AWE system dynamics and control capabilities from AWEbox. This section details the coupling of these frameworks, referred to as the aero-servo coupling. An explicit coupling strategy is employed, meaning each solver is evaluated only once per time step, with no coupling iterations performed within a single time step. This approach introduces a small offset in the system's states for each solver and could theoretically lead to coupling instabilities, but this has not been observed in the simulations.

The coupling process consists of the following steps:

1. Transmit the rigid-body motion of the aircraft, including the deflection of its control surfaces as calculated by AWEbox, to the VWE.

2. The aircraft's mesh is moved according to this rigid-body motion and the control surface deflection rates in the VWE. The flow solver then computes the resulting airflow around the aircraft, which allows for determining the aerodynamic forces and moments acting on the aircraft.

3. These calculated forces and moments are sent back to AWEbox.

4. The aircraft's movement for the next time step is updated, and subsequently, the controller determines the next control action. Figure 8 illustrates the different coupling steps for time step n.

Steps 1 and 3 are managed by the coupling tool CoCoNuT (Delaissé et al., 2021). Step 2 takes place in the VWE and step 4 in AWEbox. Note that the tether is not part of the VWE, so there is no direct interaction between the tether states and this framework. The tether force is determined by the system dynamics (Eq. 3).

Figure 8 Overview of the explicit aero-servo coupling for time step n to go to time step n+1 (for $n\in [\mathrm{0},n_{\text{end}}-\mathrm{1}]$). The numbers in brackets indicate the steps as explained in the text.

At initialization (step 0), the states and controls at t=0 s from the AWEbox output files, obtained by the flight path generation, are used to initialize the solvers. To set up the VWE, the control surface component grids are deflected with δ₀, and the aircraft is rotated using R₀ to establish the initial attitude. The aircraft is then positioned at q₀, and the flow is initialized with the wind field. At this stage, the aircraft is stationary, so the flow must first develop before it becomes meaningful. The dynamics and control in AWEbox are initialized with the states and controls at time t=0 s, including the motion. The controller is then activated to determine the control action for the first time step.

5.1 Transmit rigid-body motion (step 1)

The mesh motion of the aircraft components – the wing (w), the ailerons (a), the rudders (r), and the elevator (e) – is achieved using the zone motion function in ANSYS Fluent. This function requires the motion to be defined in terms of a translational velocity v_i, with $i\in \mathit{\{}w,e,r,a\mathit{\}}$, a rotational velocity ω_i, a rotation-axis origin o_i, and a rotation-axis direction a_i expressed in the inertial frame (see Fig. 9, left). The rotation-axis origin o_i for each component i is the position of the aircraft’s CG q_n.

Figure 9 (left) Visualization of the mesh motion parameters of the wing and (right) the aircraft's control surface hinges.

Defining the movement of the control surfaces requires special attention, as their deflection motion must be combined with the aircraft's overall motion. The deflection introduces an additional velocity component, ${\mathbf{r}}_i\times {\dot{\mathit{\delta }}}_i{\mathit{h}}_i$, when expressed relative to the aircraft's origin, and the total velocity is given by Eq. (21). The wing deflection rate ${\dot{\mathit{\delta }}}_w$ equals 0, so the second term drops for the wing. For the control surfaces, the vectors r_i and h_i represent the hinge positions and axis of the left aileron, right aileron, elevator, left rudder, and right rudder, respectively, relative to o_i, as shown in Fig. 9 (right). The angular motion due to the deflection rate is also added to the aircraft's angular motion and given by Eq. (22):

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{v}}_i={\dot{\mathbf{q}}}_n+{\mathbf{R}}_n({\mathbf{r}}_i\times {\dot{\mathit{\delta }}}_{i,n}{\mathit{h}}_i)\text{, for }i=\mathit{\{}w,e,r,a\mathit{\}},\text{(22)}& {\displaystyle }{\mathit{\omega }}_{i,\text{I}}={\mathbf{R}}_n({\mathit{\omega }}_n+{\dot{\mathit{\delta }}}_{i,n}{\mathit{h}}_i)\text{, for }i=\mathit{\{}w,e,r,a\mathit{\}}.\end{array}$$

The rotational velocity ω_i and the rotation-axis direction a_i can then be calculated as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{\omega }}_i=\Vert {\mathit{\omega }}_{i,\text{I}}\Vert \text{, for }i=\mathit{\{}w,e,r,a\mathit{\}},\text{(24)}& {\displaystyle }{\mathbf{a}}_i={\displaystyle \frac{{\mathit{\omega }}_{i,\text{I}}}{{\mathit{\omega }}_i}}\text{, for }i=\mathit{\{}w,e,r,a\mathit{\}}.\end{array}$$

Note that each component is moved independently, and they are not linked to one another. Therefore, using the same time integration scheme in both AWEbox and ANSYS Fluent is necessary to prevent drift between the components.

5.2 Solving the flow with grid movement (step 2)

The grid of the components is moved according to the zone motion function defined in the previous step. This grid movement is taken into account in the convective fluxes of the flow’s transport equations. For example, in the continuity equation, this can be expressed as

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\underset{V}{\int }\mathit{\rho }\mathrm{d}V+\underset{S}{\int }\mathit{\rho }(\mathit{v}-{\mathit{v}}_{\text{b}})\cdot \mathit{n}\mathrm{d}S=\mathrm{0}.\end{array}$$

In this equation, v is the flow velocity; v_b is the grid velocity; V and S are the volume and surface of each control volume, respectively; and n is the normal vector to the surface. After solving the flow equations, the pressure distribution is integrated, and the resulting forces ${\mathit{f}}_{\text{I},n+\mathrm{1}}^{\text{f}}$ and moments ${\mathit{m}}_{\text{I},n+\mathrm{1}}^{\text{f}}$ (around the origin of the inertial frame) are exported.

5.3 Feed back the forces and moments (step 3)

The VWE provides the forces and moments in the inertial frame I. While the toolbox AWEbox requires the forces in the inertial frame, it expects the moments in the body-fixed frame B, so the moments are transferred using Eq. (26). For this procedure, the location around which the moment is taken must correspond to the CG. Since the aircraft is first moved in the VWE, its position q_n+1 and attitude R_n+1 for the next time step are computed using the explicit Euler method (Eqs. 27 and 28). These updated values are then used to determine the moments.

$$\begin{array}{}\text{(26)}& {\displaystyle }{\mathit{m}}_{\text{B},n+\mathrm{1}}^{\text{f}}={\mathbf{\text{R}}}_{n+\mathrm{1}}^{\mathrm{T}}({\mathit{m}}_{I,n+\mathrm{1}}^{\text{f}}-{\mathit{q}}_{n+\mathrm{1}}\times {\mathit{f}}_{\text{I},n+\mathrm{1}}^{\text{f}})\text{(27)}& {\displaystyle }{\mathit{q}}_{n+\mathrm{1}}={\mathit{q}}_n+{\dot{\mathbf{q}}}_n\mathrm{\Delta }t\text{(28)}& {\displaystyle }{\mathbf{R}}_{n+\mathrm{1}}={\mathbf{R}}_n+{\dot{\mathbf{R}}}_n\mathrm{\Delta }t={\mathbf{R}}_n+{\mathbf{R}}_n{\mathit{\omega }}_{n,\times }\mathrm{\Delta }t\end{array}$$

5.4 Solving the system dynamics and control (step 4)

During the startup of the simulation (when $n<n_{\mathrm{1}}=\mathrm{220}$), the forces ${\mathit{f}}_{\text{I},n}^{\text{a}}$ and moments ${\mathit{m}}_{\text{B},n}^{\text{a}}$ from the AAM (Eqs. 9 and 11) are used, as the flow is still building up in the CFD solver. When $n_{\mathrm{1}}=\mathrm{220}<n<n_{\mathrm{2}}=\mathrm{440}$, there is a transition period during which a weighted average is taken between the forces ${\mathit{f}}_{\text{I},n+\mathrm{1}}^{\text{f}}$ and moments ${\mathit{m}}_{\text{B},n+\mathrm{1}}^{\text{f}}$ from the VWE and the AAM:

$$\begin{array}{}\text{(29)}& {\displaystyle }{\mathit{f}}_{\text{I}}=(w-\mathrm{1}){\mathit{f}}_{\text{I},n}^{\text{a}}+w{\mathit{f}}_{\text{I},n+\mathrm{1}}^{\text{f}}\text{ and}\text{(30)}& {\displaystyle }{\mathit{m}}_{\text{B}}=(w-\mathrm{1}){\mathit{m}}_{\text{B},n}^{\text{a}}+w{\mathit{m}}_{\text{B},n+\mathrm{1}}^{\text{f}}\text{, for }{n}_{\mathrm{1}}<n<n_{\mathrm{2}}.\end{array}$$

The weight w of the VWE forces and moments varies linearly from 0 at n=n₁ to 1 at n=n₂. After n>n₂, only the forces and moments from the VWE are considered. This time step corresponds to 2.2 s, after which the starting vortex is sufficiently distant and startup effects are considered negligible. Using these forces and moments, the system of differential-algebraic equations (Sect. 4.1) is solved for ${\ddot{\mathit{q}}}_n$ and ${\dot{\mathit{\omega }}}_n$, and these values are filled in the state derivative vector ${\dot{\mathbf{x}}}_n$ = ${\left\{\begin{array}{c}\dot{\mathit{q}},\ddot{\mathit{q}},\dot{\mathbf{R}},\dot{\mathit{\omega }},\dot{\mathit{\delta }},\dot{l},\ddot{l},\textcolor{red}{}\end{array}\right\}}_n$ together with the control inputs ${\dot{\mathit{\delta }}}_n$ and , previously determined from the MPC. The explicit Euler scheme is then used to update the states for the next time step. This can be expressed as

$$\begin{array}{}\text{(31)}& {\mathbf{x}}_{n+\mathrm{1}}={\mathit{x}}_n+{\dot{\mathbf{x}}}_n\mathrm{\Delta }t.\end{array}$$

Finally, the MPC is used to calculate the control inputs ${\dot{\mathit{\delta }}}_{n+\mathrm{1}}$ and for the new states, using the method explained in Sect. 4.4.

This loop (steps 1–4) continues for each new time step, starting again with step 1.

6 Results

We simulate a one-loop crosswind flight with the MegAWES aircraft in a wind field representative of offshore conditions as a demonstration for the aero-servo coupling. First, we present the optimized reference trajectory for this wind condition and aircraft parameters. Then, we present results from a simulation in which the VWE and AWE system dynamics are fully coupled, tracking the reference flight path within the VWE. Finally, we present a qualitative analysis of the flow field for this simulation and summarize the computational time required for the simulation.

6.1 Optimized reference flight path

The optimized reference trajectory, generated by the POCP detailed in Sect. 4.3 and based on the AAM, is illustrated in Figs. 10 and 11. The corresponding optimization parameters are summarized in Table B1. The aircraft starts at the top of the flight path, entering the reel-out phase (red). During this phase, the aircraft flies at an angle of attack of 4°, and the aircraft descends, causing both speed and power output to increase. The power reaches a plateau at the maximum value of 2.5 MW (put as a constraint), where it remains for approximately 5 s until the aircraft reaches the bottom of the flight path. As the aircraft ascends, the angle of attack, speed, and power output decrease, and it transitions into the reel-in phase (blue), where power is consumed. At its peak, the power required to reel in the aircraft is approximately 1.8 MW. This one-loop power cycle takes 20.0 s to complete and produces an average power output of 436 kW.

Figure 10 Visualization of the optimized reference trajectory. The grey aircraft is shown to scale at its initial position, connected by the black tether to the ground station, which is located at the origin of the inertial frame. The inlet wind field is depicted in blue, and the dashed box indicates the simulation domain of the VWE.

Figure 11 Power and airspeed plotted over time.

6.2 Tracking the reference trajectory

In this section, we demonstrate the aero-servo coupling, as outlined in Sect. 5, by tracking the reference trajectory in the VWE using MPC. From Figs. 12 and 13, we observe that the reference trajectory is tracked with high accuracy in terms of position, with a maximum deviation of 4.0 m occurring at the bottom of the loop and a root mean square deviation of 1.6 m over the cycle. The power curve is tracked with moderate accuracy, with the largest deviation of 442 kW occurring during the transition from reel out to reel in. The root mean square deviation over the cycle amounts to 93 kW, and the reduction in average power is 16 kW, which represents 4 % of the reference average power.

Figure 12 The reference trajectory (REF) and the trajectory in the coupled simulation (COSIM).

Figure 13 The reference power (REF) and the resulting power from the coupled simulation (COSIM) over time.

The aerodynamic properties of this coupled simulation (COSIM), more specifically, the angle of attack α (a), the side slip β (b), and the apparent wind speed V_a (c), are plotted in the left column of Fig. 14. The apparent wind speed remains within 0.5 m s⁻¹ from the reference values, while the angle of attack and, in particular, the side-slip angle deviate by up to 2.5° and 5°, respectively, from the reference trajectory. This increased deviation arises because the side-slip angle β is not directly associated with any aircraft state explicitly tracked by the MPC. The resulting aerodynamic forces from the VWE (blue), shown in the right column of Fig. 14 in the body-fixed frame, exhibit oscillations around the reference forces, with a maximum deviation of 0.16 for C_z, 0.02 for C_y, and 0.04 for C_x, at t=9.0 s. Future work could focus on control strategies to minimize these oscillations, thereby reducing structural stresses and extending remaining useful life (RUL).

A comparison is now made between the resulting forces obtained from the VWE and those predicted by the AAM, the model employed within the MPC framework. The resulting forces from the AAM and VWE exhibit a consistent trend. The maximum force deviations are 0.080 for C_z, 0.007 for C_y, and 0.008 for C_x, and the root mean square deviation over the cycle is 0.025 for C_z, 0.002 for C_y, and 0.004 for C_x. The CFD framework used to derive the stability derivatives for the AAM also forms the foundation of the VWE. Therefore, the remaining discrepancies likely arise from the AAM's limitations in capturing nonlinear aerodynamic effects beyond the angle of attack, as well as unsteady aerodynamic phenomena.

Figure 14 (left) The aerodynamic properties from the reference trajectory (REF) and the coupled simulation (COSIM; properties hold for both VWE and AAM). (right) The aerodynamic forces, expressed in the body-fixed frame, resulting from the VWE and predicted by the AAM.

In the left column of Fig. 15, the control surface deflections of (a) the aileron δ_a, (b) the elevator δ_e, and (c) the rudders δ_r are plotted. In addition, the angular rates (d) p, (e) q, and (f) r are shown. The control surface deflections and angular rates in the simulation generally follow the reference trajectory. The largest deviation is observed in the rudder deflection, reaching approximately 8° at t=5.0 s. The rudder deflection also hits the constraint of −10° during the simulation. Note that this constraint for trajectory generation is ±5°. Meanwhile, the maximum deviation for the angular rates amounts to approximately 9° s⁻¹ for the pitch rate at t=8.6 s. The control inputs are more aggressive than the reference, leading to greater oscillations in the aerodynamic moments around the reference value, with a maximum deviation of 0.08 for the pitch moment C_m coefficient at t=8.6 s, as seen in the right column of Fig. 15.

Figure 15 (left) Control surface deflections and (middle) rotation rates from the reference trajectory (REF) and the coupled simulation (COSIM). (right) The aerodynamic moments, expressed in the body-fixed frame, from the VWE and predicted by the AAM.

The moments predicted by the AAM and VWE generally follow the same global trend. The maximum offset between the two models amounts to 0.01 for the roll C_l, 0.03 for the pitch C_m, and 0.004 for the yaw moment coefficient C_n. The root mean square deviation between the two models over the cycle amounts to 0.006 for the roll C_l, 0.01 for the pitch C_m, and 0.001 for the yaw C_n moment coefficient. Because the roll moment coefficient C_l is normalized with the span b, in contrast with the pitch C_m that is normalized with the chord c, the offset in the roll moment is the largest in magnitude. This discrepancy in the moment coefficients arises due to the complex combination of multiple aerodynamic contributions, influenced by the highly dynamic motion of the aircraft. For example, the yaw rate r of the aircraft induces an asymmetric lift distribution, generating a roll moment. The ailerons are deflected to counteract this moment. In the AAM, these aerodynamic contributions are combined linearly. In contrast, the yaw rate r in the VWE also impacts aileron effectiveness, a phenomenon not captured in the AAM. Specifically, the rotational motion increases aileron effectiveness because the outer aileron, which experiences less flow separation (as further discussed in the next section), encounters a higher apparent wind speed. This interaction accounts for the observed offset in the AAM's predictions. Furthermore, the pitch and yaw moments are also significantly influenced by rotational rates, further underscoring the limitations of the AAM. These observations highlight the necessity of using a full CFD-based approach, as implemented in the VWE, to accurately capture the complex aerodynamic interactions.

6.3 Qualitative analysis of the flow field

The VWE provides the velocity and pressure of the flow field throughout the trajectory. These data are valuable for analyzing the design and operation of future AWE systems. In this section, we demonstrate some key flow features. Figure 16 shows the pressure coefficient distribution on the aircraft at two different times: t=5 s and t=15 s, which correspond to the highest and lowest power points in the cycle. By using the pressure distribution on the wing, the local lift distribution (represented by C_z) can be derived, as shown in Fig. 17 at four different time instances. The drop in lift is attributed to the aileron gap. Furthermore, it is evident that the lift distribution is asymmetric, a result of the rotational motion. The ailerons are deflected to compensate for this asymmetry. However, while the left aileron is deflected downward, it does not increase the local lift coefficient due to flow separation on the aileron, as seen in Fig. 18c. The resulting lift distribution deviates from the ideal elliptical lift distribution, which maximizes performance for conventional aircraft. While the ideal lift distribution for AWE systems is not yet well-established, it is clear that improvements are needed for this aircraft design to optimize performance.

Figure 16 Pressure distribution at t=5 s and t=15 s.

Figure 17 Spanwise c_z distribution at four different time instances.

A contour plot of the apparent wind velocity at t=5 s is shown in Fig. 18, providing a visualization of the flow field around all lifting surfaces. As previously mentioned, the left aileron experiences flow separation, with the right aileron also showing slight separation. Additionally, the rear portion of the main wing, located just ahead of the aileron, exhibits separated flow. This highlights the need for an improved aileron design to enhance the performance and controllability of the system and demonstrates the ability of the VWE to assess these flow phenomena for the whole power cycle. The main wing also displays slight flow separation, impacting the elevator, which operates in its wake. The flow around the elevator and rudders remains attached. Despite the equal deflection of the rudders, the flow field around them is not symmetric due to interactions between the rudders and the circular motion of the aircraft.

Figure 18 Contour plots of the apparent velocity magnitude at t=5 s. The plots show the body-fixed x–z plane at (a) y=0 m (covering the wing and elevator), (b) y=15 m (right aileron), and (c) $y=-\mathrm{15}$ m (left aileron) and the (d) x–y plane at z=2 m (rudders).

The wall shear stress in the x_B direction along the aircraft surfaces is examined, where negative values indicate flow reversal and signal regions of flow separation. Figure 19 presents this parameter at t=5 s and t=15 s. It is observed that the trailing edge of the main wing and the entire left aileron exhibit consistent flow separation throughout the trajectory. These findings provide valuable insights for enhancing the aircraft's design to optimize performance across the complete power cycle.

Figure 19 Visualization of negative wall shear stress in the body-fixed x direction (red) as an indication of separated flow regions at t=5 s and t=15 s.

6.4 Notes on computational time

The simulation was performed on a system with 2×20 core Intel Xeon Gold 6242R processors (3.1 GHz) and 187.4 GB of system memory. The peak memory usage during the simulation reached 92.2 GB. A breakdown of the computational time for one power loop cycle is provided in Table 1. Furthermore, it is calculated that 63 % of the Fluent solver calculation time is attributed to the overset technique. It can also be observed that a substantial portion of the total simulation time was spent on data saving, indicating a potential area for optimization in future runs.

The MPC algorithm required approximately 0.6 s of wall time per time step. Although this is relatively efficient in the context of the full simulation, it exceeds the simulation time step of 5 ms, making this control configuration unsuitable for real-time implementation. However, the calculation time of the MPC depends on the implementation and hardware, and faster implementations are available, such as acados (Verschueren et al., 2020).

Table 1 Computational time of a one-loop power cycle.

7 Conclusion and outlook

This study introduces a comprehensive approach that couples a virtual wind environment (VWE), represented by geometry-resolved computational fluid dynamics (CFD), with the airborne wind energy (AWE) system dynamics and control toolbox, AWEbox (De Schutter et al., 2023), to enable aero-servo simulations for AWE systems. The aero-servo coupling is demonstrated by tracking a pre-optimized one-loop reference trajectory for the MegAWES aircraft (Eijkelhof and Schmehl, 2022) in the VWE using the model predictive controller (MPC) from AWEbox. The simulation achieved 96 % of the reference power with a maximum trajectory deviation of 4 m. We compared the resulting forces and moments against predictions from an analytical aerodynamic model (AAM). The comparison revealed consistent trends, although deviations were observed due to aerodynamic effects not captured by the quasi-steady AAM, such as the nonlinear contribution of the rotational motion on the control surfaces' effectiveness and force/moment coefficients. These findings underscore the importance of employing full CFD simulations.

This analysis has highlighted key flow characteristics, such as flow separation during crosswind flight maneuvers, to inform potential design and operational improvements. Enhancing the aileron design could help prevent flow separation, thereby increasing the control authority of these surfaces and boosting the overall system performance. The simulation revealed a significant interaction between the rotational motion of the aircraft, which is common in crosswind flight, and aileron effectiveness. These insights can be further studied and used to refine the AAM, ultimately reducing the model mismatch in the MPC and improving the controller's effectiveness.

Despite its capabilities, the current framework also exhibits certain limitations. First, it should be noted that full-geometry CFD simulations are significantly more computationally demanding, rendering them impractical for rapid design iterations. To limit the computational expense of the simulations presented in this work, the boundary layer at the aircraft surface is not fully resolved; instead, wall functions are employed. While this approach reduces the accuracy in predicting flow separation, it provides a reasonable overall impression of the flow around the aircraft. Nevertheless, the current framework is compatible with further mesh refinement. Furthermore, the current CFD mesh arrangement supports only first-order discretization, which may further limit simulation accuracy. The use of higher-order schemes is recommended to enhance the precision of the CFD results. A validation study is essential to establish the credibility of the simulation outcomes. Secondly, wing deformation and fluid–structure interaction (FSI) effects are neglected in this study, based on the findings of Pynaert et al. (2023), which showed that structural deformations of the investigated aircraft were minimal. Specifically, the FSI analysis indicated only a 1.4 % reduction in aerodynamic loads. Nevertheless, the framework has been designed to support the inclusion of FSI in future studies involving more flexible aircraft configurations. Finally, it is worth noting that the MPC controller is currently executed at every time step; it is recommended to explore lower-frequency MPC evaluations.

This work represents a step forward in the ongoing development of a comprehensive aero-servo-elastic coupling, building upon prior advancements in aeroelastic modeling (Pynaert et al., 2023), which will be integrated into the current approach. Notably, although this work focuses on ground-gen systems, the framework can be adapted for fly-gen systems with necessary modifications.

Appendix A: Stability derivatives

This section outlines the calculation of the MegAWES stability derivatives using CFD. Three distinct setups, illustrated in Fig. A1, are employed for these calculations. Each setup uses identical numerical settings and the same aircraft grid as detailed in Sect. 3. All simulations for stability derivative calculations are conducted with an airspeed V_a of 80 m s⁻¹. The first and second setups involve steady-state calculations, and the third setup involves transient calculations.

The first setup focuses on calculating the stability derivatives associated with the control surface deflections δ_a,e,r and side-slip angle β. In this simulation, the front and sides of the main wing grid serve as the inlet, while the back functions as the outlet. The airspeed of the aircraft, at the specified angle of attack and side-slip angle, is applied at the inlet.

The second setup is used to calculate the stability derivative related to the roll rate p. In this simulation, a small background grid is used, with overset connectivity to the aircraft grid. The roll motion is introduced by applying the corresponding frame motion to all components of the aircraft. The aircraft’s airspeed, at the specified angle of attack, is applied at the inlet.

The third setup is used to calculate the stability derivative related to the main contribution of the angle of attack α, the yaw rate r, and the pitch rate q. This setup is similar to the VWE described in Sect. 3 but with a smaller background and zero wind velocity at the inlet. The aircraft moves according to the method described in step 1 from Sect. 5. Three flight maneuvers are considered: the first is a descending, ascending, and horizontal straight-flight maneuver to simulate a positive, negative, and zero angle of attack, respectively. The aircraft's straight motion with angle of attack α is described by

$$\begin{array}{}\text{(A1)}& \mathit{\omega }=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right],\dot{\mathit{q}}=\left[\begin{array}{c}-V_{\text{a}}\cos \mathit{\alpha }\\ \mathrm{0}\\ -V_{\text{a}}\sin \mathit{\alpha }\end{array}\right].\end{array}$$

The second and third flight maneuvers are a pure yawing flight with yaw rate r and a pure pitching flight with pitch rate q. These maneuvers are illustrated in Fig. A2. The aircraft's motion is described by the following equations for the yaw motion:

$$\begin{array}{}\text{(A2)}& {\mathit{\omega }}_r=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ r\end{array}\right],{\dot{\mathit{q}}}_r=\left[\begin{array}{c}-V_{\text{a}}\cos \mathit{\alpha }\cos \mathit{\psi }\\ -V_{\text{a}}\cos \mathit{\alpha }\sin \mathit{\psi }\\ -V_{\text{a}}\sin \mathit{\alpha }\end{array}\right],\end{array}$$

and the pitch motion (assuming small angles α and θ):

$$\begin{array}{}\text{(A3)}& {\mathit{\omega }}_q=\left[\begin{array}{c}\mathrm{0}\\ q\\ \mathrm{0}\end{array}\right],{\dot{\mathit{q}}}_q=\left[\begin{array}{c}-V_{\text{a}}\cos \mathit{\alpha }\cos \mathit{\theta }\\ \mathrm{0}\\ -V_{\text{a}}(\sin \mathit{\alpha }-\sin \mathit{\theta })\end{array}\right].\end{array}$$

In total, 53 simulations are performed to calculate the stability derivatives. For the main angle of attack contribution C_i,0, the angle of attack is varied between −10° and 10° in steps of 5°. For the other contributions, the angle of attack is varied between −5° and 5° in steps of 5°. For each angle of attack, the side slip β and control surface deflections δ_a,e,r are set to a value of 5° and 10°. Additionally, the elevator deflection δ_e is set to −5° and −10°. A value of 10° s⁻¹ and 20° s⁻¹ is used for the contribution of rotation rates.

Because each simulation considers the variation of only one contribution j, while all other contributions are set to zero, Eq. (11) can be reformulated as follows to calculate the stability derivatives:

$$\begin{array}{}\text{(A4)}& C_{i,j}={\displaystyle \frac{C_i-C_{i,\mathrm{0}}}{j}}\end{array}$$

for j = β and δ_a,e,r,

$$\begin{array}{}\text{(A5)}& C_{i,j}={\displaystyle \frac{\mathrm{2}{V}_{\text{a}}(C_i-C_{i,\mathrm{0}})}{bj}}\end{array}$$

for j = p and r, and

$$\begin{array}{}\text{(A6)}& C_{i,j}={\displaystyle \frac{\mathrm{2}{V}_{\text{a}}(C_i-C_{i,\mathrm{0}})}{cj}}\end{array}$$

for j = q.

Here, the stability derivatives C_i,j collect the contributions of the quantity j to the forces in the i direction and the moments along the i axis. The stability derivatives are subsequently fitted to a second-order polynomial function of α:

$$\begin{array}{}\text{(A7)}& C_{i,j}=\left[\begin{array}{ccc}{c}_{\mathrm{2}}& c_{\mathrm{1}}& c_{\mathrm{0}}\end{array}\right]\left[\begin{array}{c}{\mathit{\alpha }}^{\mathrm{2}}\\ \mathit{\alpha }\\ \mathrm{1}\end{array}\right].\end{array}$$

The resulting coefficients are given in Table A1 for the forces and Table A2 for the moments. AWEbox only uses the values in bold.

Figure A1 Simulation setups to calculate the stability derivatives.

Figure A2 (a) The pure r motion and (b) q motion (Mulder et al., 2013).

Table A1 Aerodynamic force stability derivatives. AWEbox only uses the values in bold.

Table A2 Aerodynamic moment stability derivatives.

Appendix B: AWE system parameters and constraints

The parameters and constraints of the wind profile and AWE system used in the simulations are summarized in Tables B1 and B2, respectively.

Table B1 Wind profile and AWE system parameters.

Table B2 AWE system constraints during path generation and tracking.