Spatio-Temporal Wind Modeling for UAV Simulations
by Kenan Cole and Adam Wickenheiser

This work originally appeared as a chapter in the dissertation entitled
Reactive Trajectory Generation and Formation Control for Groups of UAVs in
Windy Environments [1].

##

This work is licensed under a Creative Commons “Attribution 4.0” license |
---|---

###### Table of Contents

  1.     1. ‣
  2. 1 Introduction
  3. 2 Types of Wind Models
  4. 3 Over-Land Wind Turbulence Modeling
  5. 4 Over-Water Wind Turbulence Modeling
    1. 4.1 Harris Model
    2. 4.2 Forristall Model
    3. 4.3 Ochi and Shin Model
    4. 4.4 Friction Velocity Calculation
    5. 4.5 Spectral Model Comparison
  6. 5 Spatial Wind Distribution
    1. 5.1 Wind Field Propagation
    2. 5.2 Vector Components
  7. 6 Gust Model
    1. 6.1 Amplitude
    2. 6.2 Shape
    3. 6.3 Duration
    4. 6.4 Propagation
    5. 6.5 Occurrence Frequency
    6. 6.6 Lateral and Longitudinal Dissipation
  8. 7 Conclusions

###### List of Figures

  1. 1.1 Example wind data from a wind farm in the Netherlands [4] (Figure 3.3a).
  2. 1.2 Example gusts from literature
  3. 2.1 Spatial wind distribution at $z$ = 2.5 m showing the variation in the wind
  4. 2.2 Vertical wind distribution showing the variation when wind ($u$) and waves
  5. 2.3 Logarithmic wind profile based for $u_{10}$ = 10 m/s
  6. 3.1 Example time series of the six turbulence components
  7. 4.1 Comparison of the PSD for the Ochi and Shin, Forristall, and Harris models
  8. 4.2 Comparison of time series data for the Ochi and Shin, Forristall, and Harris
  9. 5.1 Spatial wavefield for input parameters $u_{10}$ = 12.8 m/s, $c_{p}$ = 7.5 m/s
  10. 5.2 Comparison of Sullivan and Ochi and Shin spectral distributions
  11. 5.3 Coordinate transformation to simulate wind fields in any direction.
  12. 5.4 Comparison of the ratio of the Von Kármán translational components
  13. 5.5 Comparison of the ratio of the Von Kármán rotational components
  14. 6.1 Relationship between the mean speed and the maximum gust amplitude
  15. 6.2 Mean gust shape characterization by comparing non-dimensional parameters
  16. 6.3 Representative gust shape defined by [4, 70]
  17. 6.4 Representative gust shape allowing independent specification
  18. 6.5 Relationship between relative gust amplitude and gust duration
  19. 6.6 Probability of detecting a gust relative to a wind turbine

###### List of Tables

  1. 3.1 Low-altitude Von Kármán model
  2. 3.2  Variable descriptions for the Von Kármán turbulence model
  3. 4.1  Variable descriptions for the Harris wind-over-water turbulence model
  4. 4.2  Variable descriptions for the Forristall wind-over-water turbulence model
  5. 4.3  Variable descriptions for the Ochi and Shin wind-over-water turbulence model

###

Nomenclature

$A_{f}$ |  experimentally determined parameter used by Forristall model (4.2)
---|---
$b$ |  aircraft wingspan, m (or ft) (3)
$B_{f}$ |  experimentally determined parameter used by Forristall model (4.2)
$c_{p}$ |  phase speed of the wind waves, m/s (4.4)
$C$ |  drag coefficient for Harris model (4.1)
$D(\omega,\theta)$ |  spreading function (5)
$D_{0}$ |  coefficient for spreading function (5)
$f$ |  frequency of wind power spectral density functions, Hz (3)
$f^{*}_{f}$ |  non-dimensional frequency used by Forristall model (4.2)
$f^{*}_{h}$ |  non-dimensional frequency used by Harris model (4.1)
$g$ |  gravity, 9.81 m/s2
$g_{1}$ |  gust amplitude coefficient (6.2)
$g_{2r}$ |  gust time scale coefficient (6.2)
$g_{2f}$ |  gust time scale coefficent (6.2)
$g_{3}$ |  gust sigmoid curve error approximation coefficient (6.2)
$g_{4}$ |  pre-gust dip coefficient (6.2)
$g_{5}$ |  post-gust dip coefficient (6.2)
$k$ |  Von Kármán constant, 0.42 (2)
$k(\omega)$ |  wavenumber (5)
$K_{gust,lat}$ |  scaling coefficient for lateral spread of gust (6.6)
$l_{x}$, $l_{y}$ |  length scale for gust probability of propagation (6.6)
$l_{gust,span}$ |  lateral spread of gust, m (6.6)
$L$ |  scaling length for Harris model, m (4.1)
$L_{u}$ |  Von Kármán model longitudinal length scale, ft (3)
$L_{v}$ |  Von Kármán model lateral length scale, ft (3)
$L_{w}$ |  Von Kármán model vertical length scale, ft (3)
$\mathbf{p}$ |  position vector in inertial frame
$\mathbf{p}_{g}$ |  goal position
$\mathbf{p}_{gust}$ |  gust starting position (6.6)
$\mathbf{p}^{w}$ |  position vector in wind frame (5)
$s$ |  spreading parameter (5)
$S(\omega)$ |  power spectral density function, m2/s (4.5)
$t$ |  simulation time, s
$t_{g1}$ |  time span from start of gust to peak, s (6.2)
$t_{g2}$ |  time span from peak to end of gust, s (6.2)
$t_{h}$ |  time span for gust at peak magnitude, m/s (6.2)
$u(z)$ |  mean wind speed at altitude $z$
$u_{10}$ |  mean wind speed at 10m altitude, m/s (2)
$u_{20}$ |  mean wind speed at 20ft altitude, ft/s (3)
$u_{*}$ |  friction velocity, m/s (2)
$v_{air}$ |  mean wind speed, m/s (5.2)
$v_{air,spatial}$ |  scalar wind speed at some ($x$,$y$) position and time $t$, m/s (5)
$v_{gust}$ |  maximum gust magnitude, m/s (6.1)
$v_{gust,nom}$ |  gust magnitude, m/s (6.2)
---|---
$v_{turb}$ |  scalar wind speed at some ($x$,$y$) position and time for any propagation direction and model, m/s (5)
$\mathbf{v}_{air}$ |  translational wind vector in inertial frame, m/s (5.2)
$x_{wa}$ |  wave age (4.4)
$x$ |  position along $\mathbf{x}_{I}$ axis, m (or ft)
$\mathbf{x}_{B}$ |  body frame $x$ axis
$\mathbf{x}_{I}$ |  inertial frame $x$ axis
$\mathbf{x}_{W}$ |  wind frame $x$ axis (5)
$y$ |  position along $\mathbf{y}_{I}$ axis, m (or ft)
$\mathbf{y}_{B}$ |  body frame $y$ axis
$\mathbf{y}_{I}$ |  inertial frame $y$ axis
$\mathbf{y}_{W}$ |  wind frame $y$ axis (5)
$z$ |  altitude, m (or ft)
$\mathbf{z}_{B}$ |  body frame $z$ axis
$\mathbf{z}_{I}$ |  inertial frame $z$ axis
$\mathbf{z}_{W}$ |  wind frame $z$ axis (5)
$z_{0}$ |  sea roughness, m (2)
$z_{0}^{*}$ |  non-dimensional sea roughness (4.4)
$\delta x$, $\delta y$ |  distance from vehicle to gust center (6.6)
$\theta$ |  spreading angle, $\pm\frac{\pi}{2}$ (5)
$\rho$ |  density of water or air, kg/m3
$\sigma$ |  wind speed variance for Forristall model, m2/s2 (4.2)
$\sigma_{air}$ |  standard deviation of wind speed over 10 minute period, m/s (6.1)
$\sigma_{u}$ |  Von Kármán model longitudinal intensity, ft/s (3)
$\sigma_{v}$ |  Von Kármán model lateral intensity, ft/s (3)
$\sigma_{w}$ |  Von Kármán model vertical intensity, ft/s (3)
$\tau$ |  shear stress, MPa (2)
$\tau_{gust}$ |  gust timespan, s (6.3)
$\Phi_{f}$ |  Forristall model power spectral density, m2/s (4.2)
$\Phi_{h}$ |  Harris model power spectral density, m2/s (4.1)
$\Phi_{os}$ |  Ochi and Shin model power spectral density, m2/s (4.3)
$\Phi_{os}^{*}$ |  Non-dimensional power spectral density for Ochi and Shin model (4.3)
$\Phi_{p}$ |  Von Kármán model pitch power spectral density, 1/s (3)
$\Phi_{q}$ |  Von Kármán model roll power spectral density, 1/s (3)
$\Phi_{r}$ |  Von Kármán model yaw power spectral density, 1/s (3)
$\Phi_{u}$ |  Von Kármán model longitudinal power spectral density, ft2/s (3)
$\Phi_{v}$ |  Von Kármán model lateral power spectral density, ft2/s (3)
$\Phi_{w}$ |  Von Kármán model vertical power spectral density, ft2/s (3)
$\psi$ |  randomly generated phase angle from 0 to $2\pi$, rad (4.5)
$\omega$ |  frequency of wind power spectral density functions, rad (3)
$\boldsymbol{\omega}_{air}$ |  rotational wind vector in inertial frame, rad/s (5.2)

##  Chapter 1 1 Introduction

Wind affects the stability and maneuverability of UAVs, which can be
particularly dangerous when operating near obstacles or each other. In order
to test the effectiveness of formation control laws and the impact of windy
environments on the vehicles, spatio-temporal wind fields must be modeled.
Each vehicle within the formation experiences unique wind conditions, but
these conditions are correlated to the conditions experienced by the other
vehicles.

A realistic wind field is not static; rather it is should include turbulence
and gusting. Turbulence is defined by the National Weather Service (NWS) as an
“abrupt or irregular movement of air that creates sharp, quick updrafts or
downdrafts” [2]. These perturbations on the mean wind speed are always present
with varying but bounded magnitudes and act in any translational or rotational
direction. Gusting on the other hand is defined by NWS as “a sudden, brief
increase in speed of the wind” [3]. Gusts are not always present, have a much
larger magnitude over the mean wind speed than turbulence, and may propagate
through the wind field before dissipating. Both of these phenomena impact the
vehicle’s ability to follow a trajectory and/or maintain a formation.

While there are many studies that addresses UAV flight control and formations,
there is a much smaller subset that considers the effects of wind on the
vehicles. Of the works that do, their wind models are simple compared to real
wind data that exhibit turbulence and gusting, such as the wind data from a
wind farm in the Netherlands [4], shown in Fig. 1.1. UAV studies such as
[5]-[9] consider wind as a constant disturbance to their controllers.
Similarly, [10] use the maximum velocity with a bounded acceleration as the
worst case scenario for their control law evaluation, and [11] assume slow
velocity changes and therefore model the wind as a piecewise constant function
in acceleration. These models do not reflect the turbulence or gusting that is
present in a real environment. A slightly more complex wind model of the form
$A\sin(\omega t+\phi)$ is used in other studies [12]-[14], but like the
constant disturbance model it is not stochastic, so it is not a good
reflection of realistic turbulence. Moving away from constant and periodic
disturbances, [15, 16] use white noise as a disturbance for UAV simulations.
The power spectral density (PSD) of white noise is uniform over all
frequencies, which does not reflect measured wind data. The Dryden turbulence
model [17] is utilized in [18],[19],[20],[21] for testing vehicle control
laws; however, this is a stationary process so gusting is not considered. Both
turbulence and gusting should be considered in the wind model since they
affect the vehicle response and control design. The controller must
continuously combat turbulence, whereas it only occasionally experiences
gusts, but the gust may exceed the vehicle’s control authority. A controller
that is robust to one may not be robust to the other, and the vehicle should
be capable of handling both.

Figure 1.1:  Example wind data from a wind farm in the Netherlands [4] (Figure
3.3a). Notice the turbulence (small perturbations about the mean values) and
the gusts (large changes in the wind speed above the mean) in the data. UAVs
operating outdoors may be subject to wind disturbances similar to this.

In addition to wind generation, the spatial distribution to model the distinct
yet correlated conditions experienced by each vehicle is needed. Spatial
distributions are considered for a few different applications depending on the
objective of the control algorithm. For example, [22] generates a spatio-
temporal wind field using basis functions and time-varying coefficients to
test UAVs’ abilities to generate trajectories for sampling, estimating, and
building a map of an unknown wind field. Their application is focused on
building a map over a large area over longer time scales than what is
considered in this dissertation. The longer time scales may miss important
turbulence and/or potentially significant gust accelerations that are
important for controller evaluation. More similar to this research are
[23],[24],[25],[26],[27], where trajectories are generated to efficiently move
individual vehicles or formations of vehicles through wind fields. The flow
fields generated are both uniform and nonuniform; however, in all cases no
turbulence is considered. Additionally, of the spatial wind field works
reviewed none consider gusting.

Gust modeling is prominent in wind turbine literature and thermal soaring
applications but is less prevalent in UAV control literature. Studies such as
[28] indicate the need for gust and turbulence modeling for UAV applications
and provide examples for generating a $(1-\cos)$ curve model that is defined
in [17]. The gust and turbulence are generated for application to a single
vehicle with no spatial considerations. Fernández et al. [29] use random data
for the $x$ and $y$ wind velocities, where any increase above the mean is
considered a “gust” by the authors. They use a fan for experimental results
but do not compare either the simulated or fan gust models to any real-world
data. Alternately, [30] use a gust shape that is similar to the gust shape
presented by [4] to model vertical gusts for thermal soaring controllers. In
both of these applications only a single vehicle is considered, and there are
no spatial variations in wind. Examining gusts from a different perspective,
[4] published a study summarizing the statistics of gusts after recording data
over a three year period at a rate of 5 Hz (to capture the high frequency,
short time scale variations) at a wind farm in the Netherlands. This work
provides useful correlations for gust and wind parameters such as magnitude,
duration, and propagation, but does not endorse a single method for spatio-
temporal gust creation. Figure 1.2 shows some of the example gust shapes from
the literature.

Figure 1.2:  Example gusts from the literature. (A) Discrete $(1-\cos)$ gust
model defined in [17]. (B) Experimental gust data from wind farm [4] (Figure
3.2). (C) Random data used to represent wind [29] (Figure 5a). (D) Normalized
wind profile for testing a soaring controller [30] (Figure 6).

Of the works reviewed, there is no one model that generates a spatio-temporal
field that includes turbulence and gusting to use to evaluate UAV controller
performance. This chapter addresses this deficiency by developing a model for
both over-land and over-water environments. In addition to providing a
representative wind field, the model should also be able to run on a personal
computer for evaluating a variety of different simulation scenarios. In the
following sections, the type of model (i.e. analytical, large-eddy simulation,
stochastic) that best satisfies the requirements for this application is
addressed, including the corresponding over-land and over-water turbulence
models. Next, the spatio-temporal spreading functions to address how each
vehicle experiences unique wind conditions that are correlated with what the
other vehicles are experiencing are developed. Lastly, the gusting models are
developed and how they interact with the spatio-temporal field is presented.

##  Chapter 2 2 Types of Wind Models

Environmental wind, turbulence, and gusting models range from very simple
analytical models to computational fluid dynamics (CFD) simulations. There are
several different types of CFD models, including Reynolds-Averaged Navier-
Stokes (RANS), large-eddy simulation (LES), and detached eddy simulation
(DES). In all of these CFD models, the area or volume of the fluid is
discretized and the equations are solved iteratively.

For the case of wind-over-water modeling, there are several groups that have
used large-eddy simulations [31]-[34] to model and understand the wind and
wave interactions. Two examples in Figs. 2.1 and 2.2 highlight example data
from Sullivan, et al. [31] where the spatial spread at a constant altitude and
the variation of the wind with altitude are shown, respectively. These
examples show the capabilities of CFD models to generate realistic and
accurate spatio-temporal wind data at high spatial resolutions. Accurate CFD
solutions, however, are highly dependent on the quality of the mesh (or grid)
as well as the analysis time step. If the problem and initial conditions are
not defined appropriately, then the results may be inaccurate and/or may even
diverge[35, 36]. Additionally, the computational requirements to generate wind
data for evaluating formation and vehicle controllers exceeds the computing
power of a personal computer, which means CFD is not feasible for use in this
work.

Figure 2.1:  Spatial wind distribution at $z$ = 2.5 m showing the variation in
the wind over water for heavier wind and waves (left) and much calmer wind and
waves (right) [31]. Figure 2.2:  Vertical wind distribution showing the
variation when wind ($u$) and waves ($c$) are in the same direction (upper),
and opposite directions (lower) [32].

On the opposite end of the spectrum, there are analytical models that describe
averaged horizontal wind profiles as a function of altitude [37, 38], but
offer no time or spatial variations. One such model is a logarithmic profile
presented in Fig. 2.3, where the dependence with altitude is governed by

$u(z)=\frac{u_{*}}{k}\ln\frac{z}{z_{0}}$ (2.1)

Figure 2.3:  Logarithmic wind profile for $u_{10}$ = 10 m/s (the wind speed at
an altitude of 10 m), $u_{*}=0.16$ m/s, and $z_{0}=0.0005$ m.

where $u_{*}$ is the friction velocity, $k$ is the Von Kármán constant taken
as 0.42, $z$ is the altitude, and $z_{0}$ is the sea roughness. The friction
velocity is the velocity at the boundary layer between the wind and water and
is sometimes called the shear velocity. If shear stress data is available,
then the friction velocity is calculated by

$u_{*}=\sqrt{\frac{\tau}{\rho}}$ (2.2)

where $\tau$ is the shear stress and $\rho$ is the density. If no shear stress
data is available, the friction velocity is solved iteratively with the sea
roughness, as discussed in Sec. 4.4.

The desired model lies between the extremes of the analytical and CFD models
such that enough of the fidelity is preserved to realistically approximate the
wind environment for a group of vehicles, but the computational expense is
reduced so that wind can be generated for vehicle formation simulations
running on a personal computer. A model that achieves this is a stochastic
turbulence model that is expanded into a spatio-temporal distribution by
spreading functions (as discussed in Sec. 5). The benefits of a stochastic
turbulence model are that the frequency content is preserved and realistic, it
is computationally feasible, and it can be combined with a gust model that
propagates through the spatio-temporal field. Since it is a stochastic
process, the flight controllers cannot be tuned to a specific wind time
history. Additionally, since the frequency content is preserved, it can be
used to determine the expected response of the vehicle. There are several
stochastic turbulence models in the literature for both over-land and over-
water environments. These models are examined more closely in Secs. 3 to 4.3
to choose one that is appropriate for this research.

##  Chapter 3 3 Over-Land Wind Turbulence Modeling

The over-land case is considered first, where the Dryden and Von Kármán
turbulence models are the standards for modeling wind turbulence. These two
models are defined in the military specification MIL-F-8785C [17] and military
handbook MIL-HDBK-1797 [39] for aircraft that are flying at both high ($>$
2000 ft or 610 m) and low altitudes ($<$ 2000 ft or 610 m). While the models
are similar in their power spectral densities and the length scale parameters,
the Von Kármán turbulence model better matches data of continuous turbulence
[40] and is the preferred model of MIL-F-8785C.

While the Von Kármán and Dryden models were developed for large aircraft,
Patel et al. [41] and Watkins and Vino [42, 43] verified that these models are
still reasonable for small UAVs at lower altitudes. Patel et al. collected
data using a small UAV on what they characterize as a “moderately gusty” day
to use for comparison with the Dryden model. They use the change in energy,
$dE/dt$, as the measure of a “gust” and compare the PSD of this data to the
Dryden PSD. The slopes and amplitudes match well; the only difference is there
are slightly larger amplitudes in the low-frequency range compared to the
Dryden model. Similarly, Watkins and Vino [42, 43] recorded wind data with a
fixed anemometer on a vehicle in moving and stationary experiments to
determine how the turbulence and wind velocities vary for aerial vehicles that
are in the low altitude range ($<$ 610 m). They also found good agreement with
the low-altitude Von Kármán model and, similar to Patel et al., saw slightly
higher amplitudes at the lower frequencies. The frequency range of interest is
0-8 Hz, where the upper limit is cited by Gage [44] based on the highest
frequency for which the Von Kármán and Dryden model transfer functions are
valid according to [39].

Since there is good agreement for both Von Kármán and Dryden for smaller UAVs
and the military specification states that the Von Kármán model shall be used
when feasible (MIL-F-8785C Section 3.7.1), the Von Kármán model is chosen for
over-land simulations. Table 3.1 summarizes the low-altitude Von Kármán model
definition and the definitions for the length scales for the low altitude
model, since the UAVs are operating below 610 m. The variables are defined in
Tbl. 3.2, which differentiates what is known and what is calculated to compute
the power spectral density.

Table 3.1: Low altitude Von Kármán model scaling lengths, intensities, and spectral density | Scaling length | Intensity | Spectral density
---|---|---|---
Lon | $L_{u}=\frac{z}{\left(0.177+0.000823z\right)^{1.2}}$ | $\sigma_{u}=\frac{\sigma_{w}}{\left(0.177+0.000823z\right)^{0.4}}$ | $\Phi_{u}=\frac{2\sigma_{u}^{2}L_{u}}{\pi u_{20}}\frac{1}{\left(1+\left(1.339\frac{L_{u}\omega}{u_{20}}\right)^{2}\right)^{5/6}}$
Lat | $L_{v}=L_{u}$ | $\sigma_{v}=\sigma_{u}$ | $\Phi_{v}=\frac{\sigma_{v}^{2}L_{v}}{\pi u_{20}}\frac{1+8/3\left(1.339\frac{L_{v}\omega}{u_{20}}\right)^{2}}{\left(1+\left(1.339\frac{L_{v}\omega}{u_{20}}\right)^{2}\right)^{11/6}}$
Vert | $L_{w}=z$ | $\sigma_{w}=0.1u_{20}$ | $\Phi_{w}=\frac{\sigma_{w}^{2}L_{w}}{\pi u_{20}}\frac{1+8/3\left(1.339\frac{L_{w}\omega}{u_{20}}\right)^{2}}{\left(1+\left(1.339\frac{L_{w}\omega}{u_{20}}\right)^{2}\right)^{11/6}}$
Roll |  |  | $\Phi_{p}=\frac{\sigma_{w}^{2}}{L_{w}u_{20}}\frac{0.8\left(\frac{\pi L_{w}}{4b}\right)^{1/3}}{1+\left(\frac{4b\omega}{\pi u_{20}}\right)^{2}}$
Pitch |  |  | $\Phi_{q}=\frac{\pm\left(\frac{\omega}{u_{20}}\right)^{2}}{1+\left(\frac{4b\omega}{\pi u_{20}}\right)^{2}}\Phi_{w}(\omega)$
Yaw |  |  | $\Phi_{r}=\frac{\pm\left(\frac{\omega}{u_{20}}\right)^{2}}{1+\left(\frac{3b\omega}{\pi u_{20}}\right)^{2}}\Phi_{v}(\omega)$
Table 3.2:  Variable descriptions for the Von Kármán turbulence model Variable |  Description |  Known/ |  Value
---|---|---|---
|  |  Calculated |
$\Phi_{u}$, $\Phi_{v}$, $\Phi_{w}$ |  PSD for longitudinal, lateral, and vertical turbulence (ft2/s) |  Calculated |  Table 3.1
$\Phi_{p}$, $\Phi_{q}$, $\Phi_{r}$ |  PSD for pitch, roll, and yaw turbulence (1/s) |  Calculated |  Table 3.1
$\sigma_{u}$, $\sigma_{v}$, $\sigma_{w}$ |  Root mean square intensities for longitudinal, lateral, and vertical turbulence velocities (ft/s) |  Calculated |  Table 3.1
$L_{u}$, $L_{v}$, $L_{w}$ |  Scaling lengths for longitudinal, lateral, and vertical turbulence (ft) |  Calculated |  Table 3.1
$\omega=2\pi f$ |  Frequency (rad) |  Known |  $f=$0-8 Hz
$z$ |  UAV altitude (ft) |  Known |  UAV Dependent
$u_{20}$ |  Mean wind speed at 20 ft (6 m) altitude (ft/s) |  Known |  Simulation Dependent
$b$ |  Aircraft wingspan (ft) |  Known |  $b=1.12$ ft (0.34 m)

The Von Kármán model PSD is used to generate a spatio-temporal wind field as
discussed in Sec. 5 for use in the over-land UAV simulations. Figure 3.1 shows
example time series data for the translational and rotational disturbances.

Figure 3.1:  Example time series generation of the six turbulence components
generated by the Von Kármán model with $u_{20}=6$ m/s, $z=10$ m, and $b=0.34$
m. (A) The three translational turbulence components that are added to the
mean wind speed along each direction. (B) The rotational components about each
axis. (C) Comparison between the three translational components and the
translational velocity induced at the motor/propeller assembly from the
rotational turbulence components. Notice that the induced translational
velocity is comparatively small for a quadrotor since the wingspan is very
short.

##  Chapter 4 4 Over-Water Wind Turbulence Modeling

The over-water environment is considered next, where there are several
stochastic models [45]-[51] that are primarily used in the field of offshore
engineering. These models are used to analyze the loads on offshore structures
and consider a wide range of environmental conditions including high winds and
hurricanes. They are based on and/or verified by data taken from various
bodies of water. To choose the best model, popular models such as the Harris
[52], Forristall [45], and Ochi and Shin [48, 53] models are examined in Secs.
4.1 to 4.3 and then compared in Sec. 4.5.

One key assumption in all of the models is a neutrally stable atmosphere. This
means that the air is well mixed so any temperature differences do not have an
effect. Additionally, neutral stability assumes that the temperature of the
water is less than the temperature of the air [54]. If this case were
reversed, then an unstable atmosphere would result. Walmsley [49] presents a
method for incorporating temperature data through use of the Monin-Obukhov
length [55], which requires absolute and potential temperature data. Hsu and
Blanchard [54] also discuss the impact that an unstable atmosphere has on
calculating the profiles through knowledge of the temperature distribution.
Since temperature data would be generated somewhat arbitrarily to feed these
alternate models, there is minimal (if any) improvement in accuracy to the
model, so the assumption of a neutrally stable atmosphere is used.

###  4.1 Harris Model

Published in 1971, the Harris model was one of the earliest over-water
stochastic models proposed [52]. It is used for high velocity wind, and when
compared to hurricane data the Harris spectrum follows the same shape but
shows some inaccuracies at particular frequencies. The Harris model PSD is
given by

$\Phi_{h}=\frac{4Cu_{10}^{2}f^{*}_{h}}{f\left(2+f^{*^{2}}_{h}\right)^{5/6}}$
(4.1)

where $C$ is the drag coefficient, $f^{*}_{h}$ is the non-dimensional
frequency, and $u_{10}$ is the wind speed at 10 m altitude. The drag
coefficient is a constant equal to 0.002 in what is described qualitatively as
“rough” wave conditions and 0.0015 in “moderate” conditions [45]. The non-
dimensional frequency for the Harris model is defined by

$f^{*}_{h}=\frac{fL}{u_{10}}$ (4.2)

where $f$ is the frequency in Hz and $L$ is a scaling length. The scaling
length was initially defined as 1800 m to fit the available data [45, 52].
Harris [52] later corrected this value by determining that the length scales
increase with altitude and defined a length scale range from approximately 50
m to 126 m for altitudes from 18 m to 182 m. This range of length scales is
closer to what Chakrabarti [56] presents in his summary of Harris’s model in
which he defines scaling lengths ranging from 60 m to 400 m.

Interestingly, the model in Eq. 4.1 was an improvement of another spectral
model, the Davenport model [57], and has the same form as the Von Kármán model
summarized in Tbl. 3.1. The coefficients differ as Harris matched experimental
wind data and Von Kármán matched wind tunnel data.

Table 4.1 summarizes the Harris model parameters differentiating what is known
and what is calculated. The Harris PSD can be used to generate spatio-temporal
wind data as described in Sec. 5.

Table 4.1:  Variable descriptions for the Harris wind-over-water turbulence model Variable |  Description |  Known/ |  Value
---|---|---|---
|  |  Calculated |
$\Phi_{h}$ |  Harris PSD (m2/s) |  Calculated |  Eq. 4.1
$L$ |  Scaling length (m) |  Known |  50 m - 126 m for altitudes 18 m - 182 m
$f^{*}_{h}$ |  Non-dimensional frequency |  Calculated |  Eq. 4.2
$f$ |  Frequency (Hz) |  Known |  $f=$0-8 Hz
$C$ |  Drag coefficient |  Known |  0.002 or 0.0015
$u_{10}$ |  Mean wind speed at 10 ft (3 m) altitude (m/s) |  Known |  Simulation Dependent

###  4.2 Forristall Model

Forristall [45] presents an alternate spectral wind profile with data compiled
from hurricanes in the Gulf of Mexico and gale force winds in the North Sea.
The NWS defines gale force winds from 17.5 - 24.2 m/s (39.1 - 54.1 mph) and
hurricane force winds at $>$ 33 m/s ($>$ 74 mph). Both weather events are
beyond the proposed operating range for the Ascending Technologies Hummingbird
UAVs; however, gusts may reach into the gale force wind range.

Forristall concluded that the data could be fit to a Gaussian distribution
defined by

$\Phi_{f}=\frac{A_{f}f_{f}^{*}\sigma^{2}}{f\left(1+B_{f}f^{*}_{f}\right)^{5/3}}$
(4.3)

where $A_{f}$ and $B_{f}$ are parameters determined from experimental data,
$f^{*}_{f}$ is the non-dimensional frequency (different from the Harris
model), $\sigma^{2}$ is the variance of wind speed, and $f$ is the frequency
of the wind in Hz. The non-dimensional frequency and variance are defined by

$\displaystyle f^{*}_{f}$ $\displaystyle=\frac{fz}{u(z)}$ (4.4)
$\displaystyle\sigma$ $\displaystyle=1.92u_{*}$ (4.5)

where $z$ is the altitude, $u(z)$ is the wind speed at altitude $z$ as
calculated by Eq. 2.1, and $u_{*}$ is the friction velocity either calculated
by Eq. 2.2 or simultaneously with the sea roughness as described in Sec. 4.4.

Table 4.2 summarizes the Forristall model parameters differentiating what is
known and what is calculated. The Forristall PSD can be used to generate
spatio-temporal wind data as described in Sec. 5.

Table 4.2:  Variable descriptions for the Forristall wind-over-water turbulence model Variable |  Description |  Known/ |  Value
---|---|---|---
|  |  Calculated |
$\Phi_{f}$ |  Forristall PSD (m2/s) |  Calculated |  Eq. 4.3
$A_{f}$ |  Constant |  Known |  $42\pm 9.07$ [45]
$B_{f}$ |  Constant |  Known |  $63\pm 13.6$ [45]
$f^{*}_{f}$ |  Non-dimensional frequency |  Calculated |  Eq. 4.4
$\sigma$ |  Wind speed variance (m2/s2) |  Calculated |  Eq. 4.5
$f$ |  Frequency (Hz) |  Known |  $f=$0-8 Hz
$z$ |  UAV altitude (m) |  Known |  UAV Dependent
$u(z)$ |  Wind speed at altitude $z$ (m/s) |  Calculated |  Eq. 2.1
$u_{*}$ |  Friction velocity (m/s) |  Calculated |  Eq. 2.2 or Sec. 4.4
$u_{10}$ |  Mean wind speed at 10 ft (3 m) altitude (m/s) |  Known |  Simulation Dependent

###  4.3 Ochi and Shin Model

The third model for comparison is the Ochi and Shin [48] model that was
published in the late 1980s and is based on fitting the averaged wind data for
seven different data sets. When compared to previous models, it shows much
greater accuracy to the averaged data, particularly at low frequencies [56].
Further supporting this, Myrhaug and Ong [46] show good agreement for the Ochi
and Shin model compared to much more recent data taken by Andersen and Løvseth
[58] off the western coast of Norway in gale and high wind conditions. The
improved accuracy comes from the piece-wise definition based on the frequency.
As with the Forristall model, this model was fit to high wind data sets for
use in offshore engineering structural analysis, so the accuracy of the
extrapolation to lower wind speed is unknown. The PSD for the Ochi and Shin
model is defined as

$\Phi_{os}=\frac{\Phi_{os}^{*}u_{*}^{2}}{f}$ (4.6)

where $\Phi_{os}^{*}$ is the non-dimensional PSD and $f$ is the frequency. The
non-dimensional PSD is defined by

$\Phi_{os}^{*}=\left\\{\begin{array}[]{l l}583f^{*}_{f}&f_{f}^{*}\leq 0.003\\\
\frac{420f_{f}^{*^{0.7}}}{\left(1+f_{f}^{*^{0.35}}\right)^{11.5}}&0.003<f_{f}^{*}\leq
0.1\\\
\frac{838f_{f}^{*}}{\left(1+f_{f}^{*^{0.3}}\right)^{11.5}}&f_{f}^{*}>0.1\end{array}\right.$
(4.7)

where the non-dimensional frequency, $f^{*}_{f}$, is defined by Eq. 4.4. Table
4.3 summarizes the Ochi and Shin model parameters differentiating what is
known and what is calculated. The Ochi and Shin PSD can be used to generate
spatio-temporal wind data as described in Sec. 5.

Table 4.3:  Variable descriptions for the Ochi and Shin wind-over-water turbulence model Variable |  Description |  Known/ |  Value
---|---|---|---
|  |  Calculated |
$\Phi_{os}$ |  Ochi and Shin PSD (m2/s) |  Calculated |  Eq. 4.6
$\Phi_{os}^{*}$ |  Non-dimensional PSD |  Calculated |  Eq. 4.7
$f^{*}_{f}$ |  Non-dimensional frequency |  Calculated |  Eq. 4.4
$\sigma$ |  Wind speed variance (m2/s2) |  Calculated |  Eq. 4.5
$f$ |  Frequency (Hz) |  Known |  $f=$0-8 Hz
$z$ |  UAV altitude (m) |  Known |  UAV Dependent
$u(z)$ |  Wind speed at altitude $z$ (m/s) |  Calculated |  Eq. 2.1
$u_{*}$ |  Friction velocity (m/s) |  Calculated |  Eq. 2.2 or Sec. 4.4

###  4.4 Friction Velocity Calculation

The friction velocity, $u_{*}$, is used in all three models described in Secs.
4.1 to 4.3. It is required for the drag coefficient (Eq. 4.1) for the Harris
model, the variance (Eq. 4.5) for the Forristall model, and the non-
dimensional PSD (Eq. 4.7) for the Ochi and Shin model. The calculation of the
friction velocity is done in two primary ways: the Charnock model or the
Volkov model [46]. Both calculations use the rearrangement of the logarithmic
wind profile in Eq. 2.1 to solve for the friction velocity at an altitude of
10 m as given by

$u_{*}=\frac{0.4u(z)}{\ln\frac{z}{z_{0}}}=\frac{0.4u_{10}}{\ln\frac{10}{z_{0}}}$
(4.8)

This equation is dependent on the sea roughness length, $z_{0}$, which is
computed differently for the Charnock and Volkov models. Regardless of the
model used, the sea roughness is dependent on the friction velocity, so the
two variables are solved iteratively.

The Charnock model was proposed in the mid 1950s and is only dependent on the
friction velocity as follows:

$z_{0}=0.0144u_{*}^{2}g$ (4.9)

where $g=9.81$ m/s2.

The Charnock model is considered less accurate than the Volkov model because
it does not take into account the wave age, $x_{wa}=c_{p}/u_{*}$, where
$c_{p}$ is the phase speed of the wind waves. The wave age is “a measure of
the time the wind has been acting on a wave group” [59]. Since the Charnock
model does not take wave age into account, Brown and Swail [47] consider the
Charnock model to be “inadequate to accurately model the wind drag over
water”, indicating the Volkov model should be used.

The Volkov model was proposed in 2001 and is a piecewise function based on the
wave age. For realistic conditions, the wave age reaches a maximum when the
wind and waves are moving at the same rate. The model defines a non-
dimensional sea roughness as follows:

$z_{0}^{*}=\left\\{\begin{array}[]{l l}0.0185&x_{wa}\leq 0.35\\\
0.03x_{wa}e^{-0.14x_{wa}}&0.35<x_{wa}<35\\\ 0.008&x_{wa}\geq
35\end{array}\right.$ (4.10)

where the sea roughness length is a function of the non-dimensional sea
roughness length:

$z_{0}=\frac{z_{0}^{*}u_{*}^{2}}{g}$ (4.11)

To iteratively solve for the sea roughness and the friction velocity, an
initial estimate of one of the parameters is needed. Brown and Swail [47]
present an approximation for the drag coefficient that is only dependent on a
known parameter, $u_{10}$, given by

$C=1000\left(0.65+0.067u_{10}\right)$ (4.12)

The estimated drag coefficient is used to calculate an initial estimate of the
friction velocity as follows:

$u_{*}=\sqrt{C}u_{10}$ (4.13)

To solve for the sea roughness, an estimate for the wave phase speed, $c_{p}$,
is also needed. The realistic limits presented in [46] of $0.03\leq
c_{p}/u_{10}\leq 1.0$ are used to choose a value of $c_{p}$. The simulations
use $c_{p}=0.5u_{10}$.

Using the initial estimate of $C$ to solve for $u_{*}$ and the estimate of
$c_{p}$, the sea roughness is solved from Eq. 4.11. This sea roughness value
updates the friction velocity in Eq. 4.8. This process is repeated until the
solutions converge. Algorithm 1 gives the pseudocode to solve for the friction
velocity, where $\varepsilon=0.001$ for the simulations.

Algorithm 1 Determine the friction velocity

1:procedure $[u_{*},z_{0},C]$ = GetFrictionVelocity($u_{10}$)

2: $C\leftarrow 1000\left(0.65+0.067u_{10}\right)$

3: $u_{*}\leftarrow\sqrt{C}u_{10}$

4: $c_{p}\leftarrow 0.5u_{10}$

5: $e_{u_{*}}\leftarrow 1$

6: while $e_{u_{*}}>\varepsilon$ do

7: $z_{0}[i]\leftarrow\frac{z_{0}^{*}u_{*}[i]^{2}}{g}$

8: $u_{*}[i+1]\leftarrow\frac{0.4u_{10}}{\ln\frac{10}{z_{0}[i]}}$

9: $e_{u_{*}}\leftarrow|u_{*}[i+1]-u_{*}[i]|$

10: $i\leftarrow i+1$

11: end while

12: $C_{d}\leftarrow u_{*}^{2}/u_{10}^{2}$

13:end procedure

The Volkov model is chosen based on the results of Brown and Swail [47] and
Myrhaug [46], who both identify the importance of including the wave age in
the friction velocity calculation. In general, the Volkov model increases the
PSD of the turbulence over the Charnock model [46], which is more constraining
for the control design. The Volkov model is used when computing the spectral
model comparisons in Sec. 4.5.

###  4.5 Spectral Model Comparison

The three spectral models (Harris, Forristall, and Ochi and Shin111The
measured wind data from [48] is given as non-dimensional data that is
converted using $\Phi_{data}=\Phi_{data}^{*}u_{*}^{2}/f$ and $f=f^{*}u(z)/z$
for $z=10$ m, $u(z)=10$ m/s, $c_{p}=5$, and $f=0$ to $8$ Hz.) defined in Secs.
4.1 to 4.3 are compared to one another and the average of measured spectrum
from the seven data sets in [48]. Figure 4.1 shows the comparison. As
expected, the Ochi and Shin model matches the average data very well. The
Forristall model is also accurate but underestimates the low frequency PSD.
The Harris model is very different from the data and other models and
underestimates the PSD.

Figure 4.1:  Comparison of the PSD for the Ochi and Shin, Forristall, and
Harris models compared to average measured wind data from [48]. As expected,
the Ochi and Shin model matches the data the closest. The Forristall model is
also close over a wide frequency range, but underestimates the low frequency
components.

A comparison of the time series data for the three models is shown in Fig.
4.2. The Ochi and Shin model shows the largest velocity magnitude, and both
the Ochi and Shin and Forristall models show high frequency content for the
wind model. The Harris model shows the smoothest time series because its PSD
is lower at higher frequencies. The time series are generated by

$y(t)=\sum_{i=1}^{n}\sqrt{S(\omega_{i})\delta\omega}\cos\left(\omega_{i}t+\psi_{i}\right)$
(4.14)

where $S(\omega_{i})$ is the PSD, $\omega$ is the frequency, and $\psi$ is a
randomly generated phase angle between 0 and $2\pi$ radians [40, 60, 61].

Figure 4.2:  Comparison of time series data for the Ochi and Shin, Forristall,
and Harris models. The Ochi and Shin model shows the largest velocity
magnitudes. The time series data was generated for $z=10$ m, $u(z)=10$ m/s,
and $f=0$ to $8$ Hz. The solution to the friction velocity uses $c_{p}=5$,
producing $u_{*}=0.45$ m/s.

Based on the comparison of the PSD and time series data for the three models,
the Ochi and Shin model is chosen for its accuracy to the averaged
experimental data. It also produces the largest magnitudes for the wind
speeds, which are more constraining for the control design.

##  Chapter 5 5 Spatial Wind Distribution

The generation of a 3D spatial wavefield distribution from directional spectra
is addressed by several groups [31, 53, 62, 63, 64] due to the availability of
data for directional spectra. Of the works reviewed, there were none that
specifically addressed the generation of a 3D spatial windfield from spectral
models; however, the principles for generating 3D wavefields can be applied
with minor modifications.

In order to account for the propagation of the wind in two dimensions, a
spreading function [56, 65] is used. There are several spreading functions;
one of the most common is the $\cos^{2s}$ function [56, 64] given by

$D(\omega,\theta)=D_{0}\cos^{2s}\theta$ (5.1)

where $s$ is a spreading parameter that is a function of frequency. The
spreading function must satisfy

$\int_{-\pi}^{\pi}D(\omega,\theta)d(\theta)=1$ (5.2)

so that there is no change in the energy of the spectrum as a result of the
spreading function. Typically, the spreading function is assumed to be
symmetric about the mean propagation direction, where $\theta=\pm\pi/2$, and
outside of this range the function is 0. Therefore, for the case $s=1$, the
coefficient $D_{0}$ is solved using

$D_{0}=\frac{1}{\int_{-\pi/2}^{\pi/2}\cos^{2s}\theta d\theta}=\frac{2}{\pi}$
(5.3)

The spreading parameter, $s$, can be specified as a constant or fit to
experimental data [66]. While $s$ can be solved as a function of a Fourier sum
of the angular spread away from the mean propagation direction from the
experimental data [67, 68], it is typically taken to be a constant [66]. In
general, increased values of $s$ tend to narrow the distribution (i.e. there
is more spreading at $\theta=0$ than $\theta=\pi/2$), so modeling of swells
tends to use larger $s$. Typical values of $s$ are in the range of 1 to 10
[56] or 1 to 16 [67]. For the applications considered in this dissertation,
$s=1$ is chosen as this is a common value [64] and does not arbitrarily narrow
the propagation.

The spreading function is multiplied by the wind spectra to form the two-
dimensional spectrum as a function of frequency and spreading angle as given
by

$S(\omega,\theta)=S(\omega)D(\omega,\theta)=S(\omega)D(\theta)=S(\omega)D_{0}\cos^{2}\theta$
(5.4)

where $S(\omega)$ is the wind PSD. The superposition used to create the
unidirectional time series in Eq. 4.14 is expanded to Eq. 5 where there is a
double summation over the frequency range and the spreading angle. The
resulting wind magnitude is

$\displaystyle v_{air,spatial}(x,y,t)=$
$\displaystyle~{}~{}\sum_{i=1}^{n}\sum_{j=i}^{m}\sqrt{2D_{0}S(\omega_{i})(\cos^{2}\theta_{j})\Delta\omega\Delta\theta}\cos\left(k_{i}x\cos\theta_{j}+k_{i}y\sin\theta_{j}-\omega_{i}t+\psi_{ij}\right)$
(5.5)

where $\omega_{i}$ varies from the lowest, $\omega_{0}$, to highest,
$\omega_{\infty}$, frequencies that are applicable for the distribution,
$\theta_{j}$ varies from $\pi/2$ to $\pi/2$, $k_{i}$ is the wavenumber, and
$\psi_{ij}$ is the phase angle. The wavenumber is the spatial frequency of the
wave propagation. The wavenumber for the wavefields is obtained from the
linear dispersion relation for ocean waves, which is $\omega^{2}=gk$ [31]. The
wave age calculated in Sec. 4.4 shows that the wind and waves do not
necessarily propagate at the same rate. The wavenumber is estimated from the
frequency range of the wind (based on the literature) and the wind speed,
$u_{10}$. This gives the spatial frequency as $k$ = $\omega/u_{10}$ [69]. The
phase angle, just as in Eq. 4.14, is a random value between 0 and $2\pi$.

An example of the spatial distribution generated by Eq. 5 is shown in Fig. 5.1
using the Ochi and Shin model. A qualitative comparison of the spreading
function wind field generation and the large eddy simulation high fidelity
model is shown in Fig. 5.2 for the same set of starting parameters. The
relative periodicity and spatial size is fairly well re-created, confirming
Eq. 5 is a reasonable method for spatio-temporal wind field generation.

Figure 5.1:  Spatial wavefield for input parameters $u_{10}$ = 12.8 m/s,
$c_{p}$ = 7.5 m/s, at an altitude at $z$ = 2.5 m. Figure 5.2:  Comparison of
(A) Sullivan (Figure 5 right panel from [31]) (B) Ochi and Shin spectral
spatial distributions over a 100 m x 100 m area for starting parameters
$u_{10}=12.8$ m/s, $z=2.5$ m, and $c_{p}/u_{*}\approx 1.4$ m/s.

###  5.1 Wind Field Propagation

The definition in Eq. 5 generates a spatial wind field that propagates with
time in the $+\mathbf{x}_{I}$ direction, as shown in Fig. 5.3. If the
prevailing wind is in a different direction, the field must be rotated using a
coordinate transformation. The new points in the frame, $\mathbf{x}_{w}$ and
$\mathbf{y}_{w}$, are defined as
$\mathbf{p}_{1}^{w}=[x^{w},~{}y^{w},~{}z^{w}]$. This is illustrated in Fig.
5.3 and used to determine the turbulence magnitude by:

$\displaystyle v_{turb}(x,y,t)$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{m}\sqrt{2S(\omega_{i},\theta_{j})\Delta\omega\Delta\theta}\cos\left(k_{i}x^{w}\cos\theta_{j}+k_{i}y^{w}\sin\theta_{j}-\omega_{i}t+\psi_{ij}\right)$
(5.6)

Figure 5.3:  Coordinate transformation to simulate wind fields in any
direction. The prevailing wind is in the $\mathbf{x}_{w}$ direction, and the
wind magnitude is computed for any point $\mathbf{p}_{1}$.

This rotation is important since MIL-F-8785C states that for the low altitude
model, the longitudinal component shall be aligned with the prevailing wind
direction.

###  5.2 Vector Components

In addition to the rotation of the wind field to align with the prevailing
wind direction, the scalar magnitude produced by Eq. 5 must be translated into
a wind vector. The Von Kármán model has four PSD components that are all
independent; however, the Ochi and Shin model has only a single PSD function.

In the Von Kármán model the PSD functions from Tbl. 3.1 are taken to be the
independent vector components. Since the longitudinal component is aligned
with the prevailing wind direction, this defines the relationship between the
longitudinal and lateral components in the inertial frame, as shown in Fig.
5.3. Equations 5.10 and 5.14 summarize the vector implementation of the full
wind velocity experienced by the vehicle.

$\displaystyle\mathbf{v}_{air}$
$\displaystyle=\left[\begin{array}[]{c}\left(v_{turb,lon}(x,y,t)+v_{air}+v_{gust}\right)\cos(\theta_{w})+v_{turb,lat}(x,y,t)\sin(\theta_{w})\\\
\left(v_{turb,lon}(x,y,t)+v_{air}+v_{gust}\right)\sin(\theta_{w})+v_{turb,lat}(x,y,t)\cos(\theta_{w})\\\
v_{turb,vert}(x,y,t)\end{array}\right]$ (5.10)
$\displaystyle\boldsymbol{\omega}_{air}$
$\displaystyle=\left[\begin{array}[]{c}\omega_{air,roll}\\\
\omega_{air,pitch}\\\ \omega_{air,yaw}\end{array}\right]$ (5.14)

where $v_{turb,lon}$, $v_{turb,lat}$, $v_{turb,vert}$, $\omega_{air,roll}$,
$\omega_{air,pitch}$, $\omega_{air,yaw}$ are computed by Eq. 5.6 using the
respective Von Kármán PSD functions in Tbl. 3.1, $v_{air,mean}$ is the mean
wind speed, and $v_{gust}$ is the magnitude of the gust (as defined in Sec.
6).

The Ochi and Shin model is not quite as straightforward since there is only a
single PSD. Hsu and Blanchard [54] provide an initial starting point for
scaling the turbulence intensities for the lateral and vertical components,
but the Von Kármán components provide insight into the relative magnitude
ratios as well. Assuming that the Ochi and Shin model is the analog of the
longitudinal component of the Von Kármán model, then the ratio of the other 5
components in the Von Kármán model to the longitudinal component are used to
scale the Ochi and Shin model PSD. In addition to magnitude, the correlation
between the translational and rotational spectra is also considered. Again
using the Von Kármán model, there are four independent components
(longitudinal, lateral, vertical, roll) and two dependent components
(pitch/vertical and yaw/lateral). This correlation is implemented by using the
same random sequence of values for $\psi$ in Eq. 5.6 for the components that
are correlated.

Figures 5.4 to 5.5 show the relative magnitudes for the translational and
rotational components, respectively, for various wind speeds. Notice that the
rotational components vary minimally with the altitude. These relationships
are used to define the components of the Ochi and Shin model for use in Eqs.
5.10 and 5.14.

Figure 5.4:  Comparison of the ratio of the Von Kármán (A) lateral PSD
($\Phi_{v}$) and (B) vertical PSD ($\Phi_{w}$) to the longitudinal PSD
($\Phi_{u}$) for varying mean wind speeds, $v_{air}$, and altitudes, $z$.
Figure 5.5:  Comparison of the ratio of the Von Kármán (A) roll PSD
($\Phi_{p}$), (B) pitch PSD ($\Phi_{q}$), and (C) yaw PSD ($\Phi_{r}$) to the
longitudinal PSD ($\Phi_{u}$) for varying mean wind speeds, $v_{air}$, and
altitudes, $z$. Note that dependence on altitude is very minimal.

##  Chapter 6 6 Gust Model

The final component to the spatio-temporal wind field is the gust model. Gusts
are only considered for the translational components of wind as, to the best
of the author’s knowledge, no rotational gust models have been proposed to
date. The gusting is based on wind farm statistical gust data from [4, 70],
which is limited to measuring only 2D gusts and does not account for boundary
layer effects near the ground. These studies suggest the following five
parameters when defining a gust:

  1. 1.

Amplitude - maximum amplitude of the gust

  2. 2.

Shape - amplitude vs time profile of the gust

  3. 3.

Duration - time span from gust start to finish

  4. 4.

Propagation speed - speed at which the gust moves in the prevailing wind
direction

  5. 5.

Occurrence frequency - time between discrete gust events

  6. 6.

Longitudinal and Lateral spread/dissipation - spatio-temporal modeling of the
gust

Each of these parameters is considered for discrete gust events and is
discussed briefly in the following subsections. The combination of these
parameters generates the final gust amplitude to use in Eq. 5.10.

###  6.1 Amplitude

The gust amplitude is strongly correlated to the mean wind speed [4]. The plot
in Fig. 6.1 is taken from [4] to show the strong correlation between mean wind
speed and gust amplitude. A linear fit is used to define the maximum gust
amplitude as a function of mean wind speed for the “POT_A=2” line from Fig.
6.1. This line uses the peak over threshold (POT) method, which identifies a
gust if the wind magnitude increases by the threshold amount (in this case 2
m/s for POT_A=2) over the mean wind speed. For the vehicles considered in this
research, 2 m/s changes could be significant. Additionally, the “POT_A=2” line
fits closely with the “IECCorr” lines, which use the correlation method that
identifies gusts by matching gust shape, and the “VelocIncr” lines, which
measure the acceleration from the first to last point of a moving window to
determine if a gust is present. The resulting linear fit is given by

$v_{gust}=\frac{52}{43}v_{air}-\frac{1}{8}$ (6.1)

Figure 6.1:  Relationship between the mean speed and the maximum gust
amplitude. The relationship is nearly linear for all the different gust
detection methods [4] (Figure 4.3).

###  6.2 Shape

It is difficult to define a generic gust shape based on the data that was
collected at the wind farm [4]. No clear spatio-temporal gust shape emerges
from the data. Instead, the authors characterize the trends of the gust shapes
relative to the gust peak by examining the data before and after the gust
peak. The resulting plot is in terms of the relative duration (in this case
data $\pm 10$ seconds around the peak with $dt=1$), compared to the non-
dimensional parameter, $(v_{gust}-v_{air})/\sigma_{air}$, where $\sigma_{air}$
is the standard deviation over a 10 minute period that is also used to
calculate the mean wind speed, $v_{air}$. Figure 6.2 shows what the authors
and other works [71] call the “mean gust shape”. While the experimental data
and theory agree, there is not an obvious way to generate a gust shape in the
time domain from this data.

Figure 6.2:  Mean gust shape characterization by comparing non-dimensional
parameters of relative duration (data before and after the gust peak, in this
case $\pm 10$ seconds with $dt=1$) and change in wind speed over standard
deviation [4] (Figure 4.12). While there is good agreement of experiment and
theory, this does not provide insight for generating time series gust data.

For this research, the gust shape considered is a discrete event that produces
a temporary increase in wind speed over the mean that then returns to its
original value. Consequently, gusts that change the overall mean wind speed
are considered separately and use the $1-\cos$ shape. A generic shape for the
discrete gust event, such as that shown in Fig. 6.3, was suggested by [70] and
is used as the basic gust shape for this dissertation.

Figure 6.3:  Representative gust shape defined by [4, 70] provides the basis
for the generic gust definition used for this dissertation.

Branlard [4] provides the equation for a symmetric curve of similar shape to
that in Fig. 6.3. To allow greater control, the definition from [4] is
generalized as shown in Fig. 6.4 and defined as follows:

$v_{gust,nom}(t)=\left\\{\begin{array}[]{ll}g_{1}\left(1-\left(g_{2r}t-g_{3}\right)^{2}\right)e^{-\left(g_{2r}t-g_{3}\right)^{2}/g_{4}},&0<t<\frac{1}{2}t_{g1}\\\
g_{1},&\frac{1}{2}t_{g1}\leq t\leq\frac{1}{2}t_{g1}+t_{h}\\\
g_{1}\left(1-\left(g_{2f}t^{*}-g_{3}\right)^{2}\right)e^{-\left(g_{2f}t^{*}-g_{3}\right)^{2}/g_{5}},&t>\frac{1}{2}t_{g1}+t_{h}\end{array}\right.$
(6.2)

where $g_{1}$ is the gust amplitude obtained from Eq. 6.1,
$g_{2r}=2g_{3}/t_{g1}$ and $g_{2f}=2g_{3}/t_{g2}$, which scale the time for
the rising and falling sides of the curve respectively, $g_{3}=6$ is the time,
$t$, where $(1-t^{2})e^{-t^{2}/2}$ is within 1% (or less) of zero, $g_{4}>0$
impacts the magnitude of the “dip” before the gust, $g_{5}>0$ impacts the
magnitude of the “dip” after the gust, $t_{h}$ is the time duration at the
peak gust magnitude, $t^{*}=t+\frac{1}{2}(t_{g2}-t_{g1})-t_{h}$, $t_{g1}$ is
the time from gust start to peak for the rising side, and $t_{g2}$ is the time
from peak to gust end for the falling side.

Figure 6.4:  Representative gust shape allowing independent specification for
the rising, falling, and peak hold portions of the gust.

###  6.3 Duration

There is a fair amount of variation in the gust duration as presented in [4]
and shown in Fig. 6.5. The linear region of the plot is used since precise
accuracy of the gust duration is not required, and the relative gust magnitude
in the linear region is consistent with the wind speeds examined for this
research. To be consistent with the gust amplitude linear fit, the “POT_A=2”
line is used again. The linear fit is given by

$\tau_{gust}=0.71\left(v_{gust}-v_{air}\right)+3.51$ (6.3)

Figure 6.5:  Relationship between relative gust amplitude and gust duration
where the relative gust amplitude is defined as $v_{gust}-v_{air}$. Reasonable
wind speeds for this dissertation generally correspond to the linear region of
the plot [4] (Figure 4.10).

###  6.4 Propagation

The data from [4] indicates a reasonable amount of variation in the data for
correlating gust propagation. There are a two parameters that stand out for
consideration: (1) the presence of the expected gust (i.e. did the gust
actually propagate as expected, and, if so, how did the shape of the gust
compare to what was experienced at the first measurement), and (2) the
propagation speed.

There is not a good model that is gleaned from this data to realistically
modify the gust propagation so that either the shape changes or the gust
dissipates out. As such, the only modifications to the gust in the prevailing
wind direction come from the changing turbulence as calculated by either the
Von Kármán or Ochi and Shin turbulence model. Given that the formation
distributions are not spanning several hundred meters or more, this assumption
is reasonable (discussed further in Sec. 6.6).

The propagation speed is the mean wind speed, $v_{air}$, as supported by [4].
It is assumed that the propagation speed is constant over the simulation,
which is reasonable for simulations on the order of tens of minutes as the
mean wind speed is typically computed over 10 minute intervals for wind data
analysis, and the 10 minute averages change slowly [4].

###  6.5 Occurrence Frequency

Branlard [4] presents a summary of the number of gusts per month as detected
by the various algorithms to determine how frequently gusts occur. Again
considering the peak over threshold model for a rise of 2 m/s, there are 9,777
gusts detected per month. This extrapolates to 13.5 gusts per hour or a gust
roughly every 4 minutes. This estimate is used as an average value for
generating gusts in the simulations.

###  6.6 Lateral and Longitudinal Dissipation

To determine the lateral and longitudinal spread of the gust, the authors of
[4] go through the data manually to determine how individual gusts propagate.
In general, there is very strong correlation ($\approx$85%) for the
longitudinal presence of a gust at a point in space 200 m away from but in
line with the initial detection location, as shown in Fig. 6.6. This is well
within the assumed operating area so it is reasonable to seek a model that
shows strong correlation longitudinally.

The lateral spread of the gust drops off significantly around points 250 m
away from and perpendicular to the initial detection location based on the
wind farm data as shown in Fig. 6.6. Since the assumed extents of the
formations are less than 250 m in any direction, the model should show some
lateral spread of gusts to other vehicles.

Figure 6.6:  Probability of detecting a gust relative to a wind turbine that
just detected a gust [4] (Figure 6.9).

The assessment of the propagation is somewhat dependent on the locations of
the wind turbines. The probabilities and correlation functions are given at
specific distances and then curve fit for all distances. Additionally, since
the data is measured on wind turbine towers, there is a minimum distance
required to avoid the wake stream of one turbine affecting the measurement of
another. The proposed curve fit [4] for probability of detecting a gust
laterally and longitudinally is

$P=e^{-\frac{\delta x}{l_{x}}}e^{-\frac{|\delta y|}{l_{y}}}$ (6.4)

where $l_{x}$ and $l_{y}$ are characteristic length parameters associated with
the gust and $\delta x$ and $\delta y$ are the relative distance between the
“gust center” and the detection point. The characteristic length is calculated
by fitting Eq. 6.4 to the data.

Equation 6.4 also provides the basis for scaling the spreading and dissipation
of the gust, where the probability is used as the gust magnitude scaling
parameter, $K_{gust,lat}=P$. The characteristic length parameters are based on
the expected wind environment and $\delta x$ and $\delta y$ are the distance
between the gust center and the vehicle location along the $\mathbf{x}_{w}$
and $\mathbf{y}_{w}$ axes, respectively. The gust center propagates at a speed
of $v_{air}$ in the $\theta_{w}$ direction, and the resulting final gust
magnitude is

$v_{gust}=K_{gust,lat}v_{gust,nom}(t^{*})$ (6.5)

where

$t^{*}=t-\frac{1}{v_{air}}\left((\mathbf{p}(t)-\mathbf{p}_{gust})\cdot[\cos\theta_{w},~{}\sin\theta_{w},~{}0]^{T}\right)$
(6.6)

$\mathbf{p}_{gust}$ is the “starting” position of the gust and $t$ is the
simulation time. This gust magnitude in Eq. 6.5 is used in Eq. 5.10 to fully
define the wind model at any point $\mathbf{p}$.

##  Chapter 7 7 Conclusions

The spatio-temporal wind model defined in this chapter addresses a deficiency
in the literature for realistic spatio-temporal wind field modeling for UAV
formations that includes turbulence and gusting. The model meets the
requirements of preserving realistic features of the environmental turbulence
and gusting as well as computational feasibility. The Von Kármán model is
chosen for operation in over-land environments where it has been shown that
this model is realistic even for smaller UAVs. The Ochi and Shin model is one
of the most accurate wind-over-water turbulence models and is chosen for the
over-water scenarios. To create the spatio-temporal distribution a spreading
function uses the chosen PSD functions to generate a correlated wind field
where each vehicle experiences unique conditions. Lastly, the spreading
function definition is also suitable for superposing with a 2D gust model.

This environmental wind model is intended for use in simulations of groups of
vehicles to simulate the unique but correlated wind conditions experienced by
the vehicles. The stochastic basis of the model allows multiple wind fields to
be used for the verifying control laws.

## Bibliography

  * [1] K. Cole.  Reactive Trajectory Generation and Formation Control for Groups of UAVs in Windy Environments.  PhD thesis, The George Washington University, 2018.
  * [2] National weather service.  http://www.crh.noaa.gov/Image/lmk/Brian_S/
NWS_2006_Turbulence_Briefing.pdf, 2006\.

  * [3] National weather service.  https://graphical.weather.gov/definitions/defineWindGust.html, 2018\.
  * [4] E. Branlard.  Wind energy: On the statistics of gusts and their propagation through a wind farm.  ECN-Wind-Memo-09-005, 2009.
  * [5] S. Jackson, J. Tisdale, M. Kamgarpour, B. Basso, and J.K. Hedrick.  Tracking controllers for small uavs with wind disturbances: theory and flight results.  47th IEEE Conference on Decision and Control, pages 564–569, 2008\.
  * [6] D.R. Nelson, D.B. Barber, T.W. McLain, and R.W. Beard.  Vector field path folloiwng for small unmanned air vehicles.  American Control Conference, Minneapolis, MN, 2006.
  * [7] A. Ratnoo, P.B. Sujit, and M. Kothari.  Adaptive optimal path following for high wind flights.  18th IFAC World Congress, Milano, Italy, 2011.
  * [8] A. Brezoescu, P. Castillo, and R. Lozano.  Straight line path following in windy conditions.  Conference on Unmanned Aerial Vehicles, Geomatics, Zurich, Switzerland., 2011.
  * [9] D. Cabecinhas, R. Cunha, and C. Silvestre.  Experimental validation of a globally stabilizing feedback controller for a quadrotor aircraft with wind disturbance rejection.  In American Control Conference, Washington, DC, pages 1024–1029, 2013.
  * [10] N. Bezzo, K. Mohta, C. Nowzari, I. Lee, V. Kumar, and G. Pappas.  Online planning for energy-efficient and disturbance-aware uav operations.  IEEE/RSJ International Conference on Intelligent Robots and Systems, 2016.
  * [11] J. Zhang, F. Meng, Y. Zhou, G. Lu, and Y. Zhong.  Decentralized formation control of multi-uav systems under wind disturbances.  Proceedings of the 34th Chinese Control Conference, 2015.
  * [12] T.G. McGee and J.K. Hedrick.  Path planning and control for multiple potin surveillance by an unmanned aircraft in wind.  American Control Conference, Minneapolis, MN, 2006.
  * [13] J. Escareno, S. Salazar, H. Romero, and R. Lozano.  Trajectory control of a quadrotor subject to 2d wind disturbances.  Journal of Intelligent Robot Systems, 70:51–63, 2013.
  * [14] I. Tanyer, E. Tatlicioglu, E. Zergeroglu, M. Deniz, A. Bayrak, and B. Ozdemirel.  Robust output tracking control of an unammed aerial vehicle subject to additive state dependent disturbance.  IET Control Theory and Applications, 2016.
  * [15] Derek Miller.  Open loop system indentificationof a micro quadrotor helicopter from closed loop data.  Master’s thesis, University of Maryland, College Park, 2011.
  * [16] F. Chen, R. Jiang, K. Zhang, B. Jiang, and G. Tao.  Robust backstepping sliding-mode control and observer-based fault estimation for a quadrotor uav.  IEEE Transactions on Industrial Electronics, 2016.
  * [17] U.S. Military Specification.  Mil-f-8785c.  1980\.
  * [18] S. Waslander and C. Wang.  Wind disturbance estimation and rejection for quadrotor position control.  In AIAA Unmanned Unlimited, Seattle, WA, 2009.
  * [19] C. Hancer, K.T. Oner, E. Sirimoglu, E. Cetinsoy, and M. Unel.  Robust hovering control of a quad tilt-wing uav.  IEEE, 978:1615–120, 2010.
  * [20] N. Sydney, B. Smyth, and D. Paley.  Dynamic control of autonomous quadrotor flight in an esimated wind field.  IEEE Conference on Decision and Control, pages 3609–3616, 2013\.
  * [21] D. Kun and I. Hwang.  Linear matrix inequality-based nonlinear adaptive robust contorl of quadrotor.  Journal of Guidance, Control, and Dynamics, 2016.
  * [22] X. Lan and M. Schwager.  Rapidly exploring random cycles: Persistent estimation of spatiotemporal fields with multiple sensing robots.  IEEE Transactions on Robotics, 2016.
  * [23] C. Peterson and D. Paley.  Multi-vehicle coordination an estimated time-varying flowfield.  Journal of Guidance, Control, and Dynamics, 2010.
  * [24] C. Peterson and D. Paley.  Distributed estimation for motion coordination in an unknown spatially varying flowfield.  Journal of Guidance, Control, and Dynamics, 2013.
  * [25] C. Peterson and Barton J.  Virtual structure formation of cooperating uavs using wind-compensation command generation and generalized velocity obstacles.  IEEE Aerospace Conference, 2015.
  * [26] R. Mellish and D. Paley.  Backstepping control design for motion coordination of self-propelled vehicles.  IEEE Conference on Decision and Control, 2010.
  * [27] N. Sydney and D. Paley.  Multivehicle coverage control for a nonstationary spatiotemporal field.  Automatica, 2014.
  * [28] J. Etele.  Overview of wind gust modemodel with application to autonomous low-level uav control.  Defence R&D Canada, 2006.
  * [29] R. Fernández, S. Dominguez, and P. Campoy.  $\mathcal{L_{1}}$ adaptive control for wind gust rejection in quad-rotor uav wind turbine inspection.  International Conference on Unmanned Aircraft Systems, 2017.
  * [30] J. Langelaan and G. Bramesfield.  Gust energy extraction for mini- and micro-uninhabited aerial vehicles.  46th Aerosciences Conference, Reno, NV, 2008.
  * [31] P.P. Sullivan, J.C. McWilliams, and T. Hristov.  A large eddy simulation model of high wind marine boundary layers above a spectrum of resolved moving waves.  National Center for Atmospheric Research.
  * [32] P. Sullivan, J. B. Edson, T. Hristov, and J. McWilliams.  Large eddy simulations and observations of atmospheric marine boundary layers above nonequilibrium surface waves.  Journal of Atmospheric Sciences, 65:1225–1245, 2007.
  * [33] N. Suzuki, T. Hara, and P. Sullivan.  Turbulent airflow at young sea states with frequent wave breaking events large eddy simulation.  Journal of the Atmospheric Sciences, 68:1290–1305, 2010.
  * [34] R.R. Harcourt and E.A. D’Asaro.  Large-eddy simulation of langmuir turbulence in pure wind seas.  Journal of Physical Oceanography, 38, 2007.
  * [35] H.C. Yee, J.R. Torczynski, S.A. Morton, M.R. Visbal, and P.K. Sweby.  On spurious behavior of cfd simulations.  International Journal for Numerical Methods in Fluids, 1999.
  * [36] J. Blazek.  Computational Fluid Dynamics: Principles and Applications.  Elsevier Science and Technology, 2005.
  * [37] K. Ruggles.  The vertical mean wind profile over the ocean for light to moderate winds.  Journal of Applied Meteorology, 9:389–395, 1970.
  * [38] J. Tambke and J.A.T. Bye.  Modelling offshore wind profiles using inertially coupled wave boundary layers.  ForWind - Center for Wind Energy Research, Oldenburg, Germany, 2001\.
  * [39] U.S. Military Handbook.  `mil-hdbk-1797`.  1997\.
  * [40] Frederic M. Hoblit.  Gust Loads on Aircraft: Concepts and Applications.  AIAA, 1988, p. 29-50.
  * [41] C. Patel, H.T. Lee, and I. Kroo.  Extracting energy from atmospheric turbulence.  XXIX OSTIV Congress, Lüsse-Berlin, Germany, 2008.
  * [42] S. Watkins and G. Vino.  The turbulent wind environment of birds, insects, and mavs.  15th Australasian Fluid Mechanics Conference, Sydney, Australia, 2004.
  * [43] S. Watkins, J. Milbank, B. Loxton, and W. Melbourne.  Atmospheric wind and their implications for microair vehicles.  AIAA Journal, Vol. 44, No. 11, 2006.
  * [44] S. Gage.  Creating a unified graphical wind turbulence model from multiple specifications.  AIAA Modeling and Simulation Technologies Conference and Exhibit, Austin, TX, 2003.
  * [45] G.Z. Forristall.  Wind spectra and gust factors over water.  Offshore Technology Conference, 1988.
  * [46] D. Myrhaug and M.C. Ong.  Effect of wave age on wind gust spectra over wind waves.  Journal of Offshore Mechanics and Arctic Engineering, 131:034501–1–6, 2009.
  * [47] R.D. Brown and V.R. Swail.  Over-water gust factors.  Ocean Engineering, 18:363–394, 1991.
  * [48] M.K. Ochi and Y.S. Shin.  Wind turbulent spectra for design consideration of offshore structure.  Offshore Technology Conference, 1988.
  * [49] J. Walmsley.  On theoretical wind speed and temperature profiles over the sea with applications to data from sable island, nova scotia.  Atmosphere - Ocean, 26:203–233, 1988.
  * [50] N. Islam, M. Zaheer, and S. Ahmed.  Double hinged articulated tower interaction with wind and waves.  Journal of Wind Engineering and Industrial Aerodynamics, 97:287–197, 2009.
  * [51] M.M. Zaheer and N. Islam.  Dynamic response of articulated tower platforms to random sea environment.  Proceedings of the ASME Conference on Ocean, Offshore, and Arctic Engineering, pages 327–332, 2009.
  * [52] R. I. Harris.  The nature of the wind.  In The modern design of wind-sensitive structures. Proceedings of the seminar held on 18 June 1970 at the Institution of Civil Engineers, Great George Street, London, SW1., 1970.
  * [53] T. Fossen.  Guidance and Control of Ocean Vehicles.  John Wiley and Sons, 1994.
  * [54] S.A. Hsu and B.W. Blanchard.  Estimating overwater turbulence intensity from routine gust factor measurements.  Journal of Applied Metereology, 43:1911–1916, 2004.
  * [55] A.M. Obukhov.  Turbulence in an atmosphere with a non-uniform temperature.  The Institute of Theoretical Geophysics of the Academy of Sciences of the U.S.S.R., 1971.
  * [56] S. Chakrabarti.  Nonlinear Methods in Offshore Engineering.  Elsevier, 1990.
  * [57] A.G. Davenport.  The spectrum of horizontal gustigust near the ground in high winds.  Quarterly Journal of the Royal Meteorological Society, 87(372):194–211, April 1961.
  * [58] O.J. Andersen and J. Løvseth.  The frøya database and maritime boundary layer wind description.  Marine Structures, 19:173–192, 2006.
  * [59] Ams glossary - wave age.  glossary.ametsoc.org/wiki/Wave_age, 2018.
  * [60] A. Abrous, R. Wamkeue, and E. Berkouk.  Modeling and simulation of a wind model using a spectral representation method.  Renewable and Sustained Energy Conference (IRSEC), 2015 3rd International, 2015.
  * [61] J.D. McMinn.  Extension of a kolmogorov atmospheric turbulence model for time-based simulation implementation.  AIAA Paper 97-3532, 1997.
  * [62] M. Kelly and J.C. Wyngaard.  Two dimensional spectra in the atmospheric boundary layer.  Journal of the Atmospheric Sciences, 63:3066–3070, 2006.
  * [63] B. Venas and L.R. Saetran.  Modeling of spatial wind-wave correlations.  AIAA Journal, 36:888–894, 1998.
  * [64] Q. Guo and Z. Xu.  Simulation of deep-water waves based on jonswap spectrum and realization by matlab.  IEEE, 2011.
  * [65] G. Forristall and K. Ewans.  Worldwide measurements of directional wave spreading.  Journal of Atmospheric and Oceanic Technology, 15:440–469, 1997\.
  * [66] Directional wave spectra using cosine-squared and cosine 2s spreading functions.  Coastal Engineering Technical Note, 1985.
  * [67] H. Mitsuyasu, F. Tasai, T. Suhara, S. Mizuno, M. Ohkusu, T. Honda, and K. Rikiishi.  Observations of the directional spectrum of ocean waves using a cloverleaf buoy.  Journal of Physical Oceanography, 5, 1975.
  * [68] D.E. Hasselmann, M. Dunckel, and J.A. Ewing.  Directional wave spectra observed during jonswap 1973.  Journal of Physical Oceanography, 10, 1980.
  * [69] Randall D. Knight.  Physics for Scientists and Engineers.  Pearson, 2004.
  * [70] J.F. Manwell, J.G. McGowan, and A.L. Rogers.  Wind Energy Explained.  J. Wiley and Sons, LTD, 2003.
  * [71] G. C. Larsen, W. Bierbooms, and Hansen. K. S.  Mean gust shapes.  DTU Library, Riso National Laboratory, Roskilde, Denmark, 2003.