import sys import os def add_to_path(p): if p and p not in sys.path:
sys.path.append(p) . . textwidth=433.62pt, textheight=505.89pt

# Vertical Airborne Wind Energy Farms with High Power Density per Ground Area
based on Multi-Aircraft Systems

Jochem De Schutter1, Jakob Harzer1, Moritz Diehl1,2 1Systems Control and
Optimization Laboratory, Department of Microsystems Engineering (IMTEK) and
Department of Mathematics, University of Freiburg, Georges-Koehler-Allee 102,
79110 Freiburg, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>freiburg.de<EMAIL_ADDRESS>of Mathematics,
University of Freiburg, Ernst-Zermelo-Strasse 1, 79104 Freiburg im Breisgau,
Germany.

###### Abstract

This paper proposes and simulates vertical airborne wind energy (AWE) farms
based on multi-aircraft systems with high power density (PD) per ground area.
These farms consist of many independently ground located systems that are
flying at the same inclination angle, but with different tether lengths, such
that all aircraft fly in a large planar elliptical area that is vertical to
the tethers. The individual systems are assigned non-overlapping flight
cylinders depending on the wind direction. Detailed calculations that take
into account Betz’ limit, assuming a cubically averaged wind power density of
7 m/s, give a potential yearly average PD of 43 MW/km2. A conventional wind
farm with typical packing density would yield a PD of 2.4 MW/km2 in the same
wind field. More refined simulations using optimal control result in a more
modest PD of 6 MW/km2 for practically recommended flight trajectories. This PD
can already be achieved with small-scale aircraft with a wing span of 5.5 m.
The simulations additionally show that the achievable PD is more than an order
of magnitude higher than for a single-aircraft AWE system with the same wing
span.

##  I Introduction

Because of the abundant availability of wind and solar energy resources, in
principle only a tiny fraction of the earth’s surface area would suffice to
generate all of humanity’s energy needs. Nevertheless the power density (PD)
per ground surface area of wind and solar power technologies is still a
relevant quantity, since the infrastructure costs of renewable energy farms,
such as grid connection and installation logistics, scale proportionally with
the farm area. The power density of existing wind power farms is estimated to
be around PD = 2 MW/km2. For solar PV farms we have a PD of around 10 MW/km2
[1].

Airborne wind energy (AWE) is an upcoming renewable energy technology which
aims at harvesting the steady and strong high-altitude winds that cannot be
reached by conventional wind technology, at only a fraction of the resources.
AWE developers mainly consider single-aircraft AWE systems (S-AWES), which are
based on the principle of one tethered aircraft flying fast crosswind
maneuvres. However, S-AWES are subject to several limitations that impede the
technology to increase PD with respect to conventional wind.

First, S-AWES are characterized by high tether drag dissipation losses. These
losses are inversely proportional to the aircraft size, which is why large and
heavy (and thus, costly) aircraft are needed to achieve the efficiency needed
for a high PD. Second, large aircraft come with an inherently large turning
radius, which corresponds to a large trajectory footprint on the ground.
Therefore, a very dense geometric spacing of units - possibly with shared
airspace - needs to be achieved for high PD. Third, since the maximum tether
length is limited due to the drag losses, S-AWES in park configuration would
all operate at similar altitudes, which, in combination with close geometric
spacing, would lead to large losses due to wake interaction [2].

$\displaystyle D$$\displaystyle d$$\displaystyle
D\sin\theta_{\mathrm{e}}$$\displaystyle\theta_{\mathrm{e}}$ Figure 1: Top and
side view of a vertical M-AWES wind farm with $N=30$,
$\theta_{\mathrm{e}}=40^{\circ}$, $d/D=0.16$ and $l_{\min}/D=0.42$.

To overcome these limitations, this work proposes (and simulates) the concept
of vertical AWE farms based on multi-aircraft AWE systems (M-AWES), as
depicted in Fig 1. In M-AWES, two or more aircraft fly very tight loops around
a shared, quasi-stationary main tether [3]. Therefore, M-AWES are very
efficient even for small aircraft size while using airspace more effectively
resulting in a lower trajectory footprint on the ground. The proposed vertical
M-AWES farm layout additionally exploits the fact that M-AWES can fly at
arbitrarily high locations above the ground, so that they can operate at
distinct locations in the sky, thereby avoiding wake interaction. A somewhat
similar idea based on networked rotary AWE systems was proposed (but not
simulated) in [4].

The remainder of this text is structured as follows. Section II describes the
vertical M-AWES park concept and determines an upper PD limit based on a
simplified analysis. Section III states the system model used in the
simulation study while Section IV proposes an optimal control problem (OCP)
formulation to compute PD-optimal power cycles. Section V describes the
numerical results of a case study, where we compute PD-optimal orbits for both
S- and M-AWES, for a small and moderate aircraft size, and where we
investigate the trade-off between PD and wing area efficiency. Section VI
discusses the main conclusions and gives suggestions for future research.

##  II Vertical M-AWES parks

The vertical M-AWES parks proposed in this work are circular and consist of
many independently ground located M-AWES that are all flying at the same
tether inclination angle $\theta_{\mathrm{e}}$, but with different tether
lengths, depending on the wind direction, such that all systems fly in a large
planar elliptical area that is vertical to the tethers, as shown in Fig 1.

The individual systems are assigned non-overlapping “operation cylinders” that
depend on the wind direction and correspond to a circular ground area.
Otherwise, the individual systems are completely independent and can e.g. be
started and landed independently. Each system might be located on a small
tower of e.g. 5 m height in order to minimize its impact on the usability of
the ground area for agriculture.

The distance from one system to all others on the ground is lower bounded by
the maximum diameter $d$ of each system’s ground area circle. The optimal
circle packing density of $N$ circles in a large circle of diameter $D$
depends on $N$ and is assumed to be $\rho_{\mathrm{circle}}=70\%$ here. Thus,
from now on we assume that

$N\frac{\pi d^{2}}{4}=\rho_{\mathrm{circle}}\frac{\pi D^{2}}{4}$ (1)

which means that we choose to build altogether
$N=\rho_{\mathrm{circle}}D^{2}/d^{2}$ individual units in the park. This
abstracts from the individual system’s size $d$ and number $N$, such that we
only need to remember the packing loss factor $\rho_{\mathrm{circle}}$.

Together, the systems form a large inclined elliptical area in the sky, and
all wings fly in this elliptical area. The area is perpendicular to the
tethers, and thus forms an ellipse with a maximum width that is equal to $D$,
but a minimum width of $D\sin\theta_{\mathrm{e}}$. The M-AWES are assumed to
be able to fly ellipses with approximately this aspect ratio so the circle
packing from the ground can be mapped by an affine transformation to the
ellipse packing in the sky. The elliptical area forms an angle of
$\theta_{\mathrm{e}}$ with the vertical, and thus, the effective area of this
“actuation ellipse” is again reduced by a factor $\cos\theta_{\mathrm{e}}$
resulting in a height of the inclined ellipse of only
$D\cos\theta_{\mathrm{e}}\sin\theta_{\mathrm{e}}$.

The shortest tether length $l_{\min}$ defines the location of the lowest point
of the ellipse, which is located $l_{\min}\sin\theta_{\mathrm{e}}$ above the
ground and $l_{\min}\cos\theta_{\mathrm{e}}$ downwind, extending the ground
boundaries of the wind park. This causes an extended park diameter
$D+2l_{\min}\cos\theta_{\mathrm{e}}$, defining a circle above which the AWE
systems can fly. This virtual area enlargement does not lead to increased
infrastructure costs and is therefore neglected.

We assume that the M-AWES wing area is adapted to reach the Betz limit
$\eta=16/27$ on the available flight area assigned to each individual system.
The overall power of the wind park is then proportional to its ground area,
but affected by a variety of losses:

  * •

the circle packing loss $\rho_{\mathrm{circle}}=70\%$;

  * •

the geometric area reduction efficiency that reaches a maximum value of
$\cos\bar{\theta}_{\mathrm{e}}\sin\bar{\theta}_{\mathrm{e}}=0.5$ for
$\bar{\theta}_{\mathrm{e}}=45^{\circ}$;

  * •

the Betz factor $\eta=16/27$;

leading altogether to an effective loss of

$\eta_{\mathrm{tot}}=\eta\cos\bar{\theta}_{\mathrm{e}}\sin\bar{\theta}_{\mathrm{e}}\rho_{\mathrm{circle}}=0.3\rho_{\mathrm{circle}}=21\%\
.$ (2)

The maximum power density is then given by

$\mathrm{PD}_{\mathrm{max}}=\frac{1}{2}\eta_{\mathrm{tot}}\rho_{\mathrm{air}}v^{3}.$
(3)

For example, for a wind speed of $v=7\ \mathrm{m/s}$ and
$\rho_{\mathrm{air}}=1.2\ \mathrm{kg/m}^{3}$ this results in a
$\mathrm{PD}_{\mathrm{max}}$ of 43 MW/km2.

To obtain an estimate for a conventional wind energy farm operating in the
same wind field, we assume operation at Betz’ limit and a circular packing
with a distance of at least 6 rotor diameters between the systems. There are
no geometric area reduction losses, resulting in a total efficiency
$\eta_{\mathrm{tot}}=1.2\%$ and a $\mathrm{PD}_{\mathrm{max}}$ of 2.4 MW/km2
when wake losses are ignored. This PD is a factor 17 lower than the potential
PD of the vertical M-AWES farm.

##  III System model

To make a more realistic assessment, we will compute detailed, PD-optimal
orbits for individual S- and M-AWES in this work. We consider “lift-mode” AWE
systems, where power is produced in a periodic fashion: first, the tether is
reeled out at high tension, driving a winch at the ground station. Then,
tether is reeled back in again at low tension, resulting in a net positive
energy gain. This section presents the M-AWES dynamics used in the optimal
control problem formulation in Section IV and introduces an averaged induction
model to account for the Betz losses.

###  III-A Multi-aircraft dynamics

In the following simulations, we use the multi-aircraft model structure and
model parameters described in [5]. The system dynamics model all six degrees-
of-freedom of the aircraft in the system. The tethers are assumed to be
straight and inelastic, which is a good assumption when tether tension is
high. The dynamics are expressed in non-minimal coordinates and summarized by
the implicit DAE

$F(\dot{x}(t),x(t),u(t),z(t),\theta,a)=0$ (4)

and consistency conditions $C(x(t))=0$.

The system variables consist first of the state vector
$x\in\mathbb{R}^{n_{\mathrm{x}}}$. The control vector
$u\in\mathbb{R}^{n_{\mathrm{u}}}$ consists of the aircraft aileron, elevator
and rudder deflection rates as well as the tether reeling acceleration. The
algebraic state $z\in\mathbb{R}^{n_{\mathrm{z}}}$ consists of the Lagrange
multipliers related to the constraints that define the interlinkage of
aircraft and tethers. The system parameters $\theta\in\mathbb{R}^{n_{\theta}}$
represent parameters that can be optimized over, such as the main tether
diameter and, in the M-AWES case, the secondary tether length and diameter.
The variable $a\in\mathbb{R}$ is the average induction factor that will be
explained in Section III-B.

For the sake of brevity, we refer the reader to [5] for a complete and formal
description of the system variables, dynamics, aerodynamic forces, consistency
conditions, etc. Here we only explicitly discuss those model components
relevant for the wind power availability. We model the wind shear in a
simplified way with a power law approximation:

$u_{\infty}(z)=u_{\mathrm{ref}}\left(\frac{z}{z_{\mathrm{ref}}}\right)^{c_{\mathrm{f}}}\
,$ (5)

with $u_{\infty}(z)$ the freestream velocity at altitude $z$ and
$u_{\mathrm{ref}}$ the reference wind speed measured at altitude
$z_{\mathrm{ref}}=100$ m, with $c_{\mathrm{f}}=0.15$ a surface friction
coefficient typical for flat, open terrain.

The atmospheric density drop with altitude is modeled using the international
standard atmosphere model [6]:

$\rho(z)\coloneqq\rho_{\mathrm{0}}\left(\frac{T_{0}-T_{\mathrm{L}}z}{T_{0}}\right)^{\frac{g}{T_{\mathrm{L}}R}-1}\
,$ (6)

where $R$ is the universal gas constant. The parameters $T_{0}$ and $\rho_{0}$
are the temperature and air density at sea level, and $T_{\mathrm{L}}$ is the
temperature lapse rate.

The model is based on a validated, small aircraft model with wing span
$b_{\mathrm{ref}}=5.5\ \mathrm{m}$, mass $m_{\mathrm{ref}}=36.5\ \mathrm{kg}$
and inertia tensor $J_{\mathrm{ref}}$ given in [7]. In order to be able to
evaluate the dynamics also for larger wing spans $b$, we utilize the following
mass upscaling formula:

$\displaystyle m$
$\displaystyle=m_{\mathrm{ref}}\left(\frac{b}{b_{\mathrm{ref}}}\right)^{\kappa}\
,\quad\text{and}\quad
J=J_{\mathrm{ref}}\left(\frac{b}{b_{\mathrm{ref}}}\right)^{\kappa+2}\ ,$ (7)

with upscaling exponent $\kappa=2.4$.

###  III-B Induction model

State-of-the-art induction models for AWE are typically a variation of Betz’
analysis for conventional wind turbines, thus based on a steady-state analysis
[8, 9]. However, detailed wind field simulations [2] show that for lift-mode
AWE systems, induction is inherently time-dependent: at the beginning of the
reel-out phase, induction starts to build up, reaching its peak when
transitioning into the reel-in phase, after which it starts to decline. When
the power cycle re-starts, the wind is almost “fresh” again.

While the induction model proposed in [10] accounts for the dynamic
variability of the aircraft trajectories over a power cycle, it still assumes
an instantaneous build-up of induction. Therefore, in this work, we propose
the following model, based on a momentum balance applied to averaged flight
quantities.

We first compute the annular swept area for each aircraft $k\in\mathcal{K}$,
with $\mathcal{K}$ the index set of all aircraft in the system. We integrate
over the reel-out phase the norm of the aircraft’s flight speed $\dot{q}$
multiplied with the wing span, weighted with the local dynamic pressure:

$A_{\mathrm{s},k}\coloneqq\int\limits_{0}^{T_{\mathrm{ro}}}\frac{1}{2}\rho(q_{\mathrm{z},k}(t))u^{2}_{\mathrm{\infty}}(q_{\mathrm{z},k}(t))b\lVert\dot{q}_{k}(t)\rVert\mathrm{d}t\
,$ (8)

where, by including the dynamic pressure inside the integral, we account for
variability of wind speed and air density along the trajectory. The parameter
$T_{\mathrm{ro}}$ is the reel-out phase duration, and $q_{\mathrm{z},k}$ is
the vertical position component of aircraft $k$.

We assume that the force acting on this annulus is the main tether force, and
that during the reel-in phase, this force is zero. The average tether force
over one power cycle of period $T$ is then given by:

$\bar{F}_{\mathrm{t}}\coloneqq\frac{1}{T}\int\limits_{0}^{T_{\mathrm{ro}}}F_{\mathrm{t}}(t)\mathrm{d}t\
,$ (9)

with the expression of the main tether force $F_{\mathrm{t}}(t)$ in [5].

Momentum conservation applied to these average quantities then gives an
algebraic equation for the average induction factor $a$, i.e.:

$\bar{F}_{\mathrm{t}}=4a(1-a)\sum\limits_{k\in\mathcal{K}}A_{\mathrm{s},k}\ .$

The apparent wind speed that each aircraft $k$ experiences is then given by

$u_{\mathrm{a},k}\coloneqq(1-a)u_{\infty}(q_{\mathrm{z},k})e_{\mathrm{x}}-\dot{q}$
(10)

with $e_{\mathrm{x}}\coloneqq\begin{bmatrix}1&0&0\end{bmatrix}^{\top}$ the
unit vector in the $x$-direction.

Similar to the steady-state models proposed in [8] and [9], this induction
model assumes that both flight annulus and tether force are perpendicular to
the wind speed vector. This is a crude assumption, given that the main tether
elevation angle is intrinsically non-zero. Nevertheless, this averaged
approach gives a first-order account of the time-dependency of the induction
for lift-mode systems, which is necessary in particular for lift-mode orbits
with short reel-out phases.

By adding the values $A_{\mathrm{s},k}$ in the M-AWES case, we assume that the
swept areas of the different aircraft do not overlap. To avoid double-
counting, a no-overlap condition will be enforced as a constraint in the
optimal control problem. For S-AWES, we will check a posteriori that the swept
area does not self-overlap.

##  IV Problem formulation

This section discusses the path constraints used for the simulations in this
work, and introduces a constraint which expresses the available flight
cylinder as a function of the radius of the corresponding circular ground
area. Then, this section presents a periodic OCP formulation to compute PD-
optimal power cycles based on the dynamics presented in the previous section.

###  IV-A Flight envelope

Path constraints need to be enforced along the trajectory, to avoid flight
envelope violations and to preserve structural integrity of the airframes and
the tethers. More concretely, we impose the following constraints:

  * •

Tether stress should not exceed the material yield stress with a safety factor
3.

  * •

The tether force should be strictly positive to avoid tether sag and to
preserve model validity.

  * •

Aircraft roll and pitch angles should be smaller than 90∘ to avoid collision
with the tether. Note that for real-world trajectories, a larger safety factor
would be appropriate.

  * •

The angle-of-attack and side-slip angle of all aircraft are bounded to avoid
stall and preserve model validity.

  * •

The aileron, elevator and rudder control surfaces and their rates are bounded.

  * •

The aircraft should remain above the ground with a safety distance of 100 m.

  * •

The tether length is bounded from above by a value $l_{\mathrm{t,max}}$ = 700
m. For the M-AWES, this constraint is active at the end of the reel-out phase.

  * •

The induction factor $a$ should be positive and smaller than 0.5 to avoid flow
acceleration or reversal.

We refer the reader to [5] for those numerical bound values not mentioned
explicitly in this text.

###  IV-B Cylindrical flight constraint

A central feature of vertical M-AWES parks is that each system is assigned an
individual, tilted flight cylinder. The intersection of this flight cylinder
with the ground gives a circular area with radius $R$ (= $d/2$). The flight
cylinder has an ellipsoidal cross-section with a major axis length $R$ and a
minor axis length $R\sin\theta_{\mathrm{e}}$, with $\theta_{\mathrm{e}}$ the
elevation angle of the cylinder. This elevation angle is included as an
optimization variable.

We express the flight cylinder constraint for all aircraft in the system in
the following way. Note that for M-AWES, the constraint is also imposed on the
juncture node between main tether and secondary tethers. First, we rotate the
aircraft position into the ellipsoidal cylinder frame:

$\displaystyle\hat{q}_{\mathrm{y},k}$ $\displaystyle\coloneqq
q_{\mathrm{y},k}$ (11) $\displaystyle\hat{q}_{\mathrm{z},k}$
$\displaystyle\coloneqq
q_{\mathrm{z},k}\cos(\theta_{\mathrm{e}})-q_{k,\mathrm{x}}\sin(\theta_{\mathrm{e}})\
.$ (12)

Then, we define the constraining ellipse axes as

$\displaystyle\hat{R}_{\mathrm{y}}$ $\displaystyle\coloneqq R-b/2$ (13)
$\displaystyle\hat{R}_{\mathrm{z}}$ $\displaystyle\coloneqq
R\sin\theta_{\mathrm{e}}-b/2\ ,$ (14)

which ensures that the entire wing with span $b$ remains in the ellipse. The
constraint then reads as

$\frac{\hat{q}_{\mathrm{z},k}^{2}}{\hat{R}_{\mathrm{z}}^{2}}+\frac{\hat{q}_{\mathrm{y},k}^{2}}{\hat{R}_{\mathrm{y}}^{2}}\leq
1\ .$ (15)

These flight constraints, together with the constraints mentioned Section IV-A
are summarized by the expression $h(x(t),u(t),z(t),\theta,a)\geq 0$.

In the M-AWES case, we need to ensure that the swept areas of the two
connected aircraft do not overlap during the reel-out phase, to preserve model
validity. For the purposes of this work, we therefore propose to pre-structure
the M-AWES OCP so that the solution consists of one single loop: half a loop
for the reel-out phase, and half a loop for the reel-in phase. During the
reel-out phase, each aircraft is assigned one half of the flight cylinder.
During the reel-in phase, the two aircraft switch flight regions.

Formally, we express the no-overlap constraint for aircraft $k$ as

$\displaystyle
h_{\mathrm{no},k}(x,\theta_{\mathrm{e}},\phi_{0})\coloneqq\hat{q}_{\mathrm{z},k}\cos(\phi_{0})-\hat{q}_{\mathrm{y},k}\sin(\phi_{0})\
.$ (16)

The angle $\phi_{0}$ rotates the intersecting half-plane that divides the
flight cylinder in two and can be chosen freely by the optimizer.

For a dual-aircraft system, with aircraft nodes $k\in\mathcal{K}=\\{2,3\\}$,
the no-overlap condition for the two aircraft is combined with a phase-fixing
constraint on the tether reel-out speed $\dot{l}_{\mathrm{t}}$:

$h_{\mathrm{no}}(x,\theta_{\mathrm{e}},\phi_{0})\coloneqq\begin{bmatrix}h_{\mathrm{no},2}(x,\theta_{\mathrm{e}},\phi_{0})\\\
-h_{\mathrm{no},3}(x,\theta_{\mathrm{e}},\phi_{0})\\\
\dot{l}_{\mathrm{t}}\end{bmatrix}\ .$ (17)

This constraint is greater or smaller than zero depending on the phase.The
constraints mentioned in this section are summarized by the expression
$h(x(t),u(t),z(t),\theta,a)\geq 0$.

###  IV-C Optimal control problem

We can now directly compute periodic flight trajectories that optimize the
power density, by solving the following periodic OCP:

$\displaystyle\underset{\displaystyle\begin{subarray}{c}x(\cdot),u(\cdot),z(\cdot)\\\
\theta,a,T_{\mathrm{ro}},T_{\mathrm{ri}}\\\
R,\theta_{\mathrm{e}},\phi_{0}\end{subarray}}{\mathrm{min}}\quad-\frac{1}{T}\int\limits_{0}^{T}\rho_{\mathrm{circle}}\frac{P(t)}{\pi
R^{2}}\mathrm{d}t$ $\displaystyle F(\dot{x}(t),x(t),u(t),z(t),\theta)$
$\displaystyle=0,$ $\displaystyle\ \forall t\in[0,T],$ $\displaystyle
h(x(t),u(t),z(t),\theta,a)$ $\displaystyle\geq 0,$ $\displaystyle\ \forall
t\in[0,T],$ $\displaystyle h_{\mathrm{no}}(x(t),\theta_{\mathrm{e}},\phi_{0})$
$\displaystyle\geq 0,$ $\displaystyle\ \forall t\in[0,T_{\mathrm{ro}}],$
$\displaystyle-h_{\mathrm{no}}(x(t),\theta_{\mathrm{e}},\phi_{0})$
$\displaystyle\geq 0,$ $\displaystyle\ \forall t\in(T_{\mathrm{ro}},T],$
$\displaystyle x(0)-x(T)$ $\displaystyle=0,$
$\displaystyle\bar{F}_{\mathrm{t}}-4a(1-a)\sum\limits_{k\in\mathcal{K}}A_{\mathrm{s},k}$
$\displaystyle=0,$

where the overall time $T$ is defined as the sum of the reel-out time and the
reel-in time, which are free optimization variables: $T\coloneqq
T_{\mathrm{ro}}+T_{\mathrm{ri}}$. The initial and final state of the
trajectory are free, but must be equal. The cost function is chosen so as to
maximize the average power output divided by the circular ground area occupied
by the system. This PD is multiplied with $\rho_{\mathrm{circle}}$ to account
for packing losses.

The M-AWES OCP is discretized using direct collocation with 40 intervals, and
Radau polynomials of degree 4. For S-AWES, the no-overlap conditions in (17)
are omitted, but the phase-fixing constraint is retained. Since the optimal
time period is larger for this problem, we increase the number of collocation
intervals to 100.

The NLP is formulated in Python using the open-source AWE optimal control
framework AWEbox [11], which builds on the symbolic framework for algorithmic
differentation and nonlinear optimization CasADi [12]. AWEbox solves the NLP
with IPOPT [13] and the linear solver MA57 [14].

##  V Numerical results

This section presents and discusses PD-optimal periodic orbits for M-AWES and
S-AWES, both for a small and moderate aircraft size. In a second step,
periodic orbits for each variant are computed for a range of fixed values for
the ground circle radius $R$, in order to investigate the trade-off between
power density and wing area efficiency. We conclude with a critical discussion
of the obtained results in the light of the modeling assumptions.

###  V-A Optimal power density solutions

The periodic OCP is solved for both a S-AWES and a M-AWES with two aircraft.
First we use the small wing span ($b=5.5$ m) of the original aircraft model.
Then we do the same for an upscaled version of the same model ($b=26$ m).
Table I summarizes the optimal results for the different variants. All of the
results have a similar optimal elevation angle
$\theta_{\mathrm{e}}^{\ast}\approx 40^{\circ}$. This is close to the
theoretically optimal value of $45^{\circ}$ obtained in Section II.

Fig. 2 shows the optimal trajectories for the small-size S-AWES and M-AWES.
The S-AWES has a rather large optimal ground circle radius $R^{\ast}=46.5$ m,
and very low power output because of the large tether drag losses. Hence the
power density is also impractically low (0.2 MW/km2). The optimal radius for
the M-AWES is a factor 2.5 smaller, at a value only slightly larger than three
wing spans. The dual-aircraft configuration thus allows the system to fly
extremely tight circles. Combined with the efficiency gain due to reduced
tether drag and the increased flying altitude, this results in a power density
that is 35 times higher (7.3 MW/km2) than for S-AWES. Optimizing for power
density drives this system to make very efficient use of the available
airspace: the optimal induction factor $a^{\ast}=0.21$ is close to the
theoretically optimal value of 1/3. For the S-AWES, induction is almost
negligible.

$\displaystyle{-20}$$\displaystyle{0}$$\displaystyle{20}$$\displaystyle\hat{z}$
[m]S-AWES$\displaystyle{-50}$$\displaystyle{-25}$$\displaystyle{0}$$\displaystyle{25}$$\displaystyle{50}$$\displaystyle\hat{y}$
[m]$\displaystyle{-10}$$\displaystyle{0}$$\displaystyle{10}$$\displaystyle\hat{z}$
[m]M-AWESreel-out Figure 2: PD-optimal flight trajectories in ellipsoidal
coordinates for $b=5.5$ m, with $\theta^{\ast}_{\mathrm{e}}\approx
40^{\circ}$.
$\displaystyle{-100}$$\displaystyle{0}$$\displaystyle{100}$$\displaystyle\hat{y}$
[m]$\displaystyle{-50}$$\displaystyle{0}$$\displaystyle{50}$$\displaystyle\hat{z}$
[m]S-AWES$\displaystyle{-100}$$\displaystyle{0}$$\displaystyle{100}$$\displaystyle\hat{y}$
[m]M-AWES - reel-out Figure 3: PD-optimal flight trajectories in ellipsoidal
coordinates for $b=26.0$ m, with $\theta_{\mathrm{e}}^{\ast}\approx
40^{\circ}$.

Fig. 3 shows the optimal trajectories for the moderate-size S-AWES and M-AWES.
The S-AWES solution improves on two fronts compared to the small-size results.
First, since the relative tether drag contribution to the total system drag
decreases with increasing aircraft size, the moderate-size system is more
efficient and flies at a higher altitude, resulting in a significantly higher
power output. Second, for a more than four times larger aircraft, the ground
area radius only increases by a factor of 2.5, the induction factor increases
by a factor of 5, and the PD (1.7 MW/km2) increases by a factor larger than 6.

For the M-AWES, the ground radius increases with a factor 5, but the large
aircraft size also reduces the impact of the secondary tether drag, which is
why the power density increases slightly to 8.4 MW/km2, which is still almost
a factor 5 larger than for the moderate-size S-AWES.

Note that while in the envisioned vertical wind farms, M-AWES need an
individual flight cylinder to avoid collisions, this need not be the case for
S-AWES. In fact, S-AWES can be packed closer together than done in this work,
even with overlapping flight cones when synchronized properly, as proposed in
[15]. From this perspective, the obtained PD results are too pessimistic.

However, the model used in this work neglects wake interaction effects within
the farm, which is a good assumption for a vertical M-AWES farm, but not for
densely and horizontally packed S-AWES, in particular as they grow larger. The
induction factors obtained in this study are small ($<$ 0.05) but they might
still have a non-negligible effect on the power output of downstream systems.
In this sense, the proposed packing densities might not be overly
conservative.

TABLE I: PD-optimal (A-D) and practically recommended (E) solution parameters and outputs. Label | System |  $b$ [m] |  $a$ [-] |  $R$ [m] | PD [$\frac{\text{MW}}{\text{km}^{2}}$] |  $\bar{P}$ [kW]
---|---|---|---|---|---|---
A | S-AWES | 5.5 | 0.01 | 46.5 | 0.2 | 2.0
B | M-AWES | 5.5 | 0.21 | 18.4 | 7.3 | 11.0
C | S-AWES | 26.0 | 0.05 | 116.8 | 1.7 | 107.0
D | M-AWES | 26.0 | 0.09 | 96.0 | 8.4 | 347.6
E | M-AWES | 5.5 | 0.13 | 30.4 | 5.9 | 24.5

###  V-B Trade-off between ground and wing area

Optimizing for power density results in a suboptimal solution in terms of
average power output for a given wing area
$\bar{P}_{\mathrm{S}}\coloneqq\bar{P}/(\lvert\mathcal{K}\rvert S)$, with $S$
the aerodynamic surface of a single aircraft in the system. In practice, a
trade-off between these two objectives needs to be found: for a given power
output, we want to both minimize the trajectory footprint and the required
wing area.

Starting from the PD-optimal solution, the Pareto front between PD and
$\bar{P}_{\mathrm{S}}$ is constructed by re-solving the OCP for fixed and
increasing values of $R$. Fig 4 shows the result of this parametric sweep for
all variants. After a certain value of $R$, the $\bar{P}_{\mathrm{S}}$-optimal
solution is reached and the cylindrical flight constraint becomes inactive for
all larger $R$.

For increasing $R$, all system variants are able to increase power output in
two ways. First, the systems fly at lower, more power-optimal elevation
angles, down to a value $\theta_{\mathrm{e}}\approx 25^{\circ}$ for all
variants. The larger value of $R$ compensates the minor ellipse axis reduction
with $\sin\theta_{\mathrm{e}}$. Second, the systems fly trajectories with a
larger harvesting area, which results in a lower induction factor but overall
in a net increase in power.

The M-AWES power output increases by a factor up to 2.9 for the small-size
system and by up to a factor of 1.7 for the moderate-size system, at the cost
reducing the PD with a factor of 2.6 in both cases. A good compromise for the
small-size system might be PD $\approx$ 6 MW/km2 and
$\bar{P}_{\mathrm{S}}\approx$ 4 kW/m2, marked by the label “E” in Fig. 4 and
summarized in Table I. Interestingly, the small-size M-AWES completely
dominates the moderate-size S-AWES by a large margin. M-AWES based on small
aircraft can thus be efficiently deployed both as a single unit for small-
scale applications, as well as in AWE farms for utility-scale electricity
generation.

$\displaystyle{0}$$\displaystyle{1}$$\displaystyle{2}$$\displaystyle{3}$$\displaystyle{4}$$\displaystyle{5}$$\displaystyle\bar{P}_{\mathrm{S}}$
[kW/m2]$\displaystyle{0}$$\displaystyle{2}$$\displaystyle{4}$$\displaystyle{6}$$\displaystyle{8}$PD
[W/m2]ABCDEM-AWES, b = 5.5 mM-AWES, b = 26 mS-AWES, b = 5.5 mS-AWES, b = 26 m
Figure 4: Pareto efficiency front between power density and power per wing
area for small- and moderate-size S-AWES and M-AWES.

##  VI Conclusion

In this paper we proposed vertical M-AWES farms with high PD per ground area.
We determined the theoretical potential of these farms and computed and
compared detailed PD-optimal flight trajectories for both M-AWES and S-AWES of
different sizes. The achieved PD of the recommended small-size M-AWES design
“E” is significantly lower than the theoretical estimate, by a factor of 7. A
big loss factor is the fact that the optimal flight annulus only covers part
of the elliptical cross-section, and thus does not exploit the total available
harvesting area. Future work should explore M-AWES trajectories that use more
area of the elliptical cylinders in order to achieve power densities that are
closer to what is theoretically possible.

## Acknowledgements

This research was supported by DFG via Research Unit FOR 2401 and project
424107692 and by the EU via ELO-X 953348.

## References

  * [1] J. Van Zalk, P. Behrens, The spatial extent of renewable and non-renewable power generation: A review and meta-analysis of power densities and their application in the u.s, Energy Policy 123 (2018) 83–91.
  * [2] T. Haas, J. De Schutter, M. Diehl, J. Meyers, Large-eddy simulation of airborne wind energy farms, Wind Energy Science 7 (3) (2022) 1093–1135.
  * [3] M. Zanon, S. Gros, J. Andersson, M. Diehl, Airborne wind energy based on dual airfoils, IEEE Transactions on Control Systems Technology 21 (2013) 1215–1222.
  * [4] R. Read, Kite networks for harvesting wind energy, in: Airborne Wind Energy: Advances in Technology Development and Research, Springer Singapore, 2018.
  * [5] J. De Schutter, R. Leuthold, T. Bronnenmeyer, R. Paelinck, M. Diehl, Optimal control of stacked multi-kite systems for utility-scale airborne wind energy, in: Proceedings of the IEEE Conference on Decision and Control (CDC), 2019.
  * [6] C. Archer, An introduction to meteorology for airborne wind energy, in: Airborne Wind Energy, Springer Berlin / Heidelberg, 2013.
  * [7] E. C. Malz, J. Koenemann, S. Sieberling, S. Gros, A reference model for airborne wind energy systems for optimization and control, Renewable Energy 140 (2019) 1004–1011.
  * [8] M. De Lellis, R. Reginatto, R. Saraiva, A. Trofino, The betz limit applied to airborne wind energy, Renewable Energy 127 (2018) 32–40.
  * [9] M. Kheiri, V. S. Nasrabad, F. Bourgault, A new perspective on the aerodynamic performance and power limit of crosswind kite systems, Journal of Wind Engineering and Industrial Aerodynamics.
  * [10] R. Leuthold, J. De Schutter, E. C. Malz, G. Licitra, S. Gros, M. Diehl, Operational regions of a multi-kite awe system, in: European Control Conference (ECC), 2018.
  * [11] https://github.com/awebox (2022).
  * [12] J. Andersson, J. Gillis, G. Horn, J. Rawlings, M. Diehl, CasADi – a software framework for nonlinear optimization and optimal control, Mathematical Programming Computation 11 (1) (2019) 1–36.
  * [13] A. Wächter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical Programming 106 (1) (2006) 25–57.
  * [14] HSL, A collection of Fortran codes for large scale scientific computation. (2011).
  * [15] P. Faggiani, R. Schmehl, Airborne Wind Energy: Advances in Technology Development and Research, Springer Singapore, 2018, Ch. Design and Economic of a Pumping Kite Wind Park, pp. 391–411.