# Low-resource finetuning of foundation models beats
state-of-the-art in histopathology

###### Abstract

To handle the large scale of whole slide images in computational pathology,
most approaches first tessellate the images into smaller patches, extract
features from these patches, and finally aggregate the feature vectors with
weakly-supervised learning. The performance of this workflow strongly depends
on the quality of the extracted features. Recently, foundation models in
computer vision showed that leveraging huge amounts of data through supervised
or self-supervised learning improves feature quality and generalizability for
a variety of tasks. In this study, we benchmark the most popular vision
foundation models as feature extractors for histopathology data. We evaluate
the models in two settings: slide-level classification and patch-level
classification. We show that foundation models are a strong baseline. Our
experiments demonstrate that by finetuning a foundation model on a single GPU
for only two hours or three days depending on the dataset, we can match or
outperform state-of-the-art feature extractors for computational pathology.
These findings imply that even with little resources one can finetune a
feature extractor tailored towards a specific downstream task and dataset.
This is a considerable shift from the current state, where only few
institutions with large amounts of resources and datasets are able to train a
feature extractor. We publish all code used for training and evaluation as
well as the finetuned models 111https://github.com/beneroth13/dinov2.

Index Terms—  Self-supervised learning, foundation models, medical imaging,
histopathology

##  1 Introduction

Fig. 1: We propose finetuning a DINOv2 ViT-S, which yields at least equal
performance compared to CTransPath and RetCCL but in a fraction of domain
specific training time. Performance is measured on three datasets: TCGA &
CPTAC (WSI-level classification) and NCT-CRC (patch-level classification).

Recent progress in natural language processing, in particular the success of
large language models [1], has shown that huge data sources can be leveraged
to build models that generalize well to a wide variety of tasks. These models,
also called foundation models [2], are trained on large-scale datasets without
task-specific supervision and can be adapted for specialized downstream
applications.

Similar approaches have been developed for computer vision, where the most
influential aspects for a good performance lie in large-scale datasets and
self-supervised learning techniques. The core concept is to efficiently
represent complex visual data in a lower-dimensional space. One common
approach is contrastive learning, which involves pushing similar instances
closer in the learned representation space while separating dissimilar ones
[3]. Training networks with contrastive loss often needs a large batch size
and has the risk to collapse to one representation. To account for this, non-
contrastive methods, such as MoCov3 [4], were proposed. MoCov3 creates a
dynamic dictionary of negative examples using a momentum update mechanism. A
similar approach is student-teacher network training, where the student only
gets a masked version of the input and tries to learn similar representations
as the teacher network. The student is updated during training whereas the
teacher is updated with the exponential moving average of the student weights
[5]. This strategy can be used for knowledge distillation, where the knowledge
from a larger model is transferred to a smaller one.

Translating these efforts to medical tasks would be highly beneficial but is
challenging for several reasons: Large-scale and publicly available data
sources are scarce, imaging modalities and downstream tasks are very diverse,
high performance on very specific tasks is preferred over a good performance
on a broad set of tasks, and, finally, the patients’ privacy needs to be
preserved. Nevertheless, a few models were proposed as potential foundation
models that are trained on medical imaging data, but not extensively on
histopathology data [6, 7].

In computational pathology, self-supervised feature extractors are commonly
used to compute a low-dimensional representation of a whole slide image that
was tessellated into smaller patches [8, 9, 10]. These low-dimensional
features are then usually aggregated for a downstream task using weakly-
supervised learning. The largest dataset used for pre-training these feature
extractors comprises 15 million patches [10]. However, this is still
magnitudes smaller than the large-scale pre-training datasets used in
classical computer vision, e.g., 142 million images curated from a dataset of
1.2 billion images for DINOv2 [11]. There has been a recent effort to scale-up
the pre-training for histopathology-specific feature extractors using private
datasets [12, 13, 14] but none of the models is currently publicly available
and could be included in this study.

Chen et al. [12] report using 32 NVIDIA A100 80GB GPUs and a batch size of
3,072 for 125,000 iterations without specifying the exact training time.
CTransPath was trained for 250 hours on 48 NVIDIA V100 GPUs [10]. RetCCL [9]
used 32 NVIDIA V100 32G GPU for 300 hours. Virchow [13] and REMEDIS [6] do not
specify their setup but the scale of the data and the trained model suggest a
similar compute heavy setup. As the compute resources needed to train such a
model is out of reach for many institutions, only publicly available models
such as CTransPath [10] or RetCCL [9] are used.

In our work, we propose a shift towards low-resource fine-tuning on specific
datasets and show its applicability. We present a benchmark of four of the
largest foundation models in classical computer vision on histopathology
downstream tasks and compare the models to the state-of-the-art domain-
specific feature extractors CTransPath and RetCCL. We finetune ViT models of
the best performing self-supervised learning method DINOv2 on a single GPU and
show that finetuning foundation models on task-specific data can outperform
large-scale trained domain-specific models while needing only a fraction of
resources and training time (Figure 1).

##  2 Experiments & Methods

Our comparative analysis includes ResNet50 [15], pre-trained on ImageNet and
its truncated variant as baseline models. We benchmark four vision foundation
models that have been trained on natural images: ImageBind, a multi-modal
model trained on paired data consisting of images and a corresponding second
modality such as text or depth-maps [16]; Segment Anything (SAM), a
segmentation model, trained with supervision on a large amount of images with
segmentation masks [17]; BEiT, a vision transformer (ViT) trained with self-
supervision based on masked image modeling [18]; DINOv2, a self-supervised
teacher-student model trained on a large curated dataset [11]. As
histopathology-specific feature extractors, we evaluate the state-of-the-art
models CTransPath [10] and RetCCL [9], which are specifically designed for
patches from whole-slide images and are trained on 15 million histopathology
patches.

###  2.1 Experimental setup

We evaluate the methods on three datasets of colorectal cancer tissue, two
with annotations for slide-level classification and one with patch-level
annotations.

Slide-level classification. The public repository The Cancer Genome Atlas
(TCGA) [19] contains whole-slide images of 632 patients with colorectal cancer
(CRC) from the cohorts COAD and READ with annotations of microsatellite
instability (MSI). We perform 5-fold cross validation using three folds for
training, one for validation, and one for testing. As an external test set, we
use the CRC cohort of the public database The Clinical Proteomic Tumor
Analysis Consortium, (CPTAC) [20] with 110 patients.

We tessellate the whole slide images at 10$\times$ magnification in patches of
size $256\times 256$ pixels (px) and extract the features222Code for the pre-
processing pipeline available at https://github.com/peng-
lab/HistoBistro/tree/feature_extraction.. Subsequently, we train a slide-level
classifier using transformer-based aggregation following Wagner et al. [21].

Patch-level classification. The NCT-CRC-100K (NCT) dataset [22] contains
100,000 non-overlapping image patches from H&E-stained WSIs from nine classes
of colorectal cancer and normal colon tissue. The images are of size
$224\times 224$ pixels (px) at 0.5 $\mu$m/px magnification and are color-
normalized using Macenko’s method [23]. We use the full dataset for training.
For testing, we use CRC-VAL-HE-7K [22], an independent dataset of 7,180
patches from an external cohort with the same tissue classes.

We extract the features using the introduced models as feature extractors and
evaluate the classification using the sklearn $K$-nearest neighbors
implementation with $K=20$ as well as the sklearn implementation of logistic
regression classifier with 1,000 iterations and $\ell_{2}$-regularization
coefficient $\lambda$ of $\frac{100}{MC}$, where $M$ is the embedding
dimension and $C$ is the number of classes. The choice of hyperparameters and
evaluation methods follows Chen et al. [24].

|  | AUROC
---|---|---
Model | Histo | TCGA | CPTAC
ResNet50 |  | $0.67{\scriptstyle\pm 0.10}$ | $0.65{\scriptstyle\pm 0.05}$
ResNet50 truncated |  | $0.68{\scriptstyle\pm 0.02}$ | $0.68{\scriptstyle\pm 0.02}$
SAM (ViT-B) |  | $0.56{\scriptstyle\pm 0.07}$ | $0.62{\scriptstyle\pm 0.03}$
SAM (ViT-H) |  | $0.55{\scriptstyle\pm 0.05}$ | $0.64{\scriptstyle\pm 0.02}$
BEiT (ViT-B) |  | $0.52{\scriptstyle\pm 0.07}$ | $0.51{\scriptstyle\pm 0.02}$
ImageBind (ViT-h) |  | $0.64{\scriptstyle\pm 0.06}$ | $0.71{\scriptstyle\pm 0.02}$
DINOv2 (ViT-S) |  | $0.73{\scriptstyle\pm 0.07}$ | $0.72{\scriptstyle\pm 0.06}$
DINOv2 (ViT-g) |  | $0.66{\scriptstyle\pm 0.12}$ | $0.60{\scriptstyle\pm 0.03}$
RetCCL (ResNet50) | $\checkmark$ | $0.84{\scriptstyle\pm 0.08}$ | $0.77{\scriptstyle\pm 0.05}$
CTransPath (Swin-T) | $\checkmark$ | $\underline{0.85}{\scriptstyle\pm 0.05}$ | $\underline{0.82}{\scriptstyle\pm 0.04}$
DINOv2 (ViT-S) | $\checkmark$ | $\textbf{0.89}{\scriptstyle\pm 0.05}$ | $\textbf{0.85}{\scriptstyle\pm 0.02}$
DINOv2 (ViT-g) | $\checkmark$ | $0.84{\scriptstyle\pm 0.05}$ | $0.79{\scriptstyle\pm 0.03}$
Table 1: Finetuned DINOv2 outperforms dedicated feature extractors CTransPath and RetCCL on AUROC scores on WSI classification (MSI detection) in CRC tissue from TCGA and CPTAC (external). Also, DINOv2 ViT-S notably outperforms the ViT-g variant. |  | 20-NN | linear probe
---|---|---|---
Model | Histo | ACC | F1 | ACC | F1
ResNet50 |  | 0.78 | 0.81 | 0.86 | 0.88
ResNet50 truncated |  | 0.87 | 0.91 | 0.88 | 0.89
SAM (ViT-B) |  | 0.75 | 0.76 | 0.83 | 0.81
SAM (ViT-H) |  | 0.69 | 0.71 | 0.79 | 0.84
BEiT (ViT-B) |  | 0.63 | 0.66 | 0.59 | 0.65
ImageBind (ViT-h) |  | 0.87 | 0.89 | 0.90 | 0.91
DINOv2 (ViT-S) |  | 0.88 | 0.90 | 0.90 | 0.92
DINOv2 (ViT-g) |  | 0.91 | 0.93 | 0.92 | 0.94
RetCCL (ResNet50) | $\checkmark$ | 0.91 | 0.93 | 0.92 | 0.94
CTransPath (Swin-T) | $\checkmark$ | 0.95 | 0.96 | 0.94 | 0.95
DINOv2 (ViT-S) | $\checkmark$ | 0.94 | 0.95 | 0.93 | 0.94
DINOv2 (ViT-g) | $\checkmark$ | 0.93 | 0.95 | 0.94 | 0.96
Table 2: Finetuned DINOv2 matches the performance of the best-scoring feature
extractor, CTransPath, on balanced accuracy (ACC) and weighted F1-scores (F1)
for CRC tissue classification on NCT-patches.

###  2.2 Finetuning and implementation details

To evaluate performance improvements of the foundation models by finetuning,
we focus on DINOv2 in the smallest and largest variants ViT-S/14 (21M
parameters, Table 3) and ViT-g/14 (1.1B parameters), respectively. We finetune
both models using the original DINOv2
implementation333https://github.com/facebookresearch/dinov2 on the NCT and
TCGA dataset. As a starting point, the publicly available pretrained weights
for the backbones were loaded into both the student and the teacher network.
Since the weights for the DINO-heads are not publicly available, a random
initialization is used. As an additional augmentation we add vertical
flipping.

Similar as proposed in [11], we use input images of size $224\times 224$ px
with a local crop size of $98\times 98$ px. In order to adapt the positional
encoding from its original size of $518\times 518$ px to the desired size of
$224\times 224$ px, we employ bicubic interpolation.

In our finetuning experiments, we utilize the provided configuration file from
the training pipeline with a base learning rate of $2\times 10^{-4}$. To
account for the comparatively small datasets, we adjust the number of
iterations per pseudo-epoch and the number of pseudo-epochs. As batch-size we
use 256 for ViT-S and 32 for ViT-g. Using a smaller batch size of 128 for
ViT-S yielded similar results as with batch size 256. In the KoLeoLoss the
hyperparameter $\epsilon$ was changed from $1\times 10^{-8}$ to $1\times
10^{-4}$ to avoid infinite values in the loss. The remaining hyperparameters
are kept consistent with the values specified in the configuration file
provided in the official implementation. Training was done on a single NVIDIA
A100 GPU (80GB).

Slide-level classification. For finetuning on TCGA, we tessellate the whole
slide image into patches of size $512\times 512$ px at $20\times$
magnification. We thereby obtain a training dataset of 2.5M patches, which are
downsampled using RandomResizeCrop to $224\times 224$ px as input to the
vision transformer. We set the iterations per epoch to $500$ and the number of
pseudo-epochs to $200$.

Patch-level classification. We use the full training dataset NCT of 100,000
images for finetuning and set the iterations per epoch to $100$ and the number
of pseudo-epochs to $100$, respectively.

###  2.3 Results

Fig. 2: Performance over time of finetuning a ViT-s with DINOv2: a) on NCT-CRC
and evaluating on the external NCT-CRC testset on patch-level classification
and b) on TCGA and testing on TCGA (5-fold cross-validation) and CPTAC
(external testset) on WSI-level classification.

Slide-level classification. In our experiments, DINOv2 outperforms all other
computer vision foundation models, i.e., SAM, BEiT, and ImageBind, in both the
small (ViT-S) and giant (ViT-g) variant (Table 1). Notably, the giant variant
performs substantially worse than the small variant (e.g., 0.66 vs. 0.73 on
TCGA, 0.6 vs. 0.72 on CPTAC). The performance of the larger network is
surpassed by the truncated version of ResNet50 trained on ImageNet, which
consistently outperforms its full counterpart. Interestingly, the features
from ImageBind lead to a strong generalization on CPTAC, where the external
validation almost matches the performance of DINOv2 (0.71 vs. 0.72) despite a
low performance on the in-domain test set TCGA (0.64). These models are all
outperformed by the histopathology feature extractors, RetCCL and CTransPath,
showing that feature quality benefits from training on histopathology-specific
data. However, finetuning DINOv2 on our 2.5M patches from TCGA further
improves the performance and outperforms the state-of-the-art feature
extractors, where the smaller variant based on ViT-S performs better than the
large version. Compared to CTransPath we only utilize 0.6% (72 h vs. 12,000 h)
of GPU hours to reach peak performance for training the ViT-S models (Figure
1). Peak performance is reached at approximately 30,000 to 40,000 iterations
for ViT-S compared to around 60,000 for ViT-g (Figure 2).

Patch-level classification. As for slide-level classification, DINOv2
outperforms all other computer vision foundation models, where the margin to
ImageBind and the truncated ResNet50 version is small. Histopathology-specific
pre-training improves also in this case, and CTransPath achieves best results
accross all metrics (Table 2). However, RetCLL obtains an almost identical
performance. Finetuning both DINOv2 variants also almost match the performance
of CTransPath with only two GPU hours (0.08% of CTransPath’s) of training
(Figure 1). The peak performance is approximately reached at 2,000 iterations
for ViT-S compared to around 10,000 for ViT-g (Figure 2).

In both tasks, all foundation models show a very similar behavior: ImageBind
and the truncated ResNet50 follow DINOv2 regarding the performance, where
ImageBind has roughly 20 times more parameters (Table 3). Conversely, both
variants of SAM as well as BEiT underperform compared to other models. The
training duration for NCT is significantly shorter than for TCGA due to the
smaller number of slides in the training set. To ensure no overfitting is
happening, we evaluate each task on an external test set, CPTAC and CRC-VAL-
HE-7K, repsecitvely, as well as the in-domain TCGA test set.

Model | #params | feature size
---|---|---
ResNet50 | 25.6 M | 1024
ResNet50 truncated | 8.5 M | 1024
SAM (ViT-B) | 91 M | 256
SAM (ViT-H) | 636 M | 256
BEiT (ViT-B) | 85.7 M | 768
ImageBind (ViT-h) | 630 M | 1024
RetCCL (ResNet50) | 23.5 M | 2048
CTransPath (Swin-T) | 27.5 M | 768
DINOv2 (ViT-S) | 21 M | 384
DINOv2 (ViT-g) | 1.1 B | 1536
Table 3: Parameter count and feature dimension of benchmarked models.

##  3 Conclusion

In our study, we benchmarked four of the most popular publicly available
state-of-the-art foundations models in the field of computer vision as feature
extractors. By fine-tuning the most promising model with the DINOv2 framework,
we achieved notable improvements that are comparable or better than the
current state-of-the-art feature extractor in histopathology CTransPath, while
only using $0.08\%$ or $0.6\%$ of their compute budget, which amounts to two
hours or three days of training on one single A100 GPU depending on the
dataset. Our experiments were conducted with a limited scope, utilizing two
datasets for training and three for testing purposes. To fully conclude a
superior performance compared to established feature extractors, more
experiments are needed on more diverse benchmarks.

##  4 Ethics statement

This research study was conducted retrospectively using human subject data
made available in open access. Ethical approval was not required as confirmed
by the license attached with the open access data.

##  5 Acknowledgments

V.K. and S.J.W. were supported by the Helmholtz Association under the joint
research school “Munich School for Data Science - MUDS”. S.J.W. was supported
by the Add-on Fellowship of the Joachim Herz Foundation. This work was also
supported by the BMBF-funded de.NBI Cloud within the German Network for
Bioinformatics Infrastructure (de.NBI) (031A532B, 031A533A, 031A533B,
031A534A, 031A535A, 031A537A, 031A537B, 031A537C, 031A537D, 031A538A). C.M.
has received funding from the European Research Council under the European
Union’s Horizon 2020 research and innovation program (grant agreement number
866411) and is supported by the Hightech Agenda Bayern. The results shown here
are in part based upon data generated by the TCGA Research Network
444https://www.cancer.gov/tcga.

## References

  * [1] Wayne Xin Zhao, Kun Zhou et al.,  “A survey of large language models,”  arXiv preprint arXiv:2303.18223, 2023.
  * [2] Rishi Bommasani, Drew A. Hudson et al.,  “On the opportunities and risks of foundation models,” 2022.
  * [3] Ting Chen, Simon Kornblith et al.,  “A simple framework for contrastive learning of visual representations,” 2020.
  * [4] Xinlei Chen, Saining Xie et al.,  “An empirical study of training self-supervised vision transformers,” 2021.
  * [5] Chengming Hu, Xuan Li et al.,  “Teacher-student architecture for knowledge learning: A survey,” 2022.
  * [6] Shekoofeh Azizi, Laura Culp et al.,  “Robust and data-efficient generalization of self-supervised machine learning for diagnostic imaging,”  Nature Biomedical Engineering, pp. 1–24, 2023.
  * [7] Michael Moor, Qian Huang et al.,  “Med-flamingo: a multimodal medical few-shot learner,” 2023.
  * [8] Abtin Riasatian, Morteza Babaie et al.,  “Fine-tuning and training of densenet for histopathology image representation using tcga diagnostic slides,”  Medical Image Analysis, vol. 70, pp. 102032, 2021.
  * [9] Xiyue Wang, Yuexi Du et al.,  “Retccl: clustering-guided contrastive learning for whole-slide image retrieval,”  Medical image analysis, vol. 83, pp. 102645, 2023.
  * [10] Xiyue Wang, Sen Yang et al.,  “Transformer-based unsupervised contrastive learning for histopathological image classification,”  Medical Image Analysis, vol. 81, pp. 102559, 2022.
  * [11] Maxime Oquab, Timothée Darcet et al.,  “Dinov2: Learning robust visual features without supervision,” 2023.
  * [12] Richard J. Chen, Tong Ding et al.,  “A general-purpose self-supervised model for computational pathology,” 2023.
  * [13] Eugene Vorontsov, Alican Bozkurt et al.,  “Virchow: A million-slide digital pathology foundation model,”  arXiv preprint arXiv:2309.07778, 2023.
  * [14] Gabriele Campanella, Ricky Kwan et al.,  “Computational pathology at health system scale–self-supervised foundation models from three billion images,”  arXiv preprint arXiv:2310.07033, 2023.
  * [15] Kaiming He, Xiangyu Zhang et al.,  “Deep residual learning for image recognition,”  CoRR, vol. abs/1512.03385, 2015.
  * [16] Rohit Girdhar, Alaaeldin El-Nouby et al.,  “Imagebind: One embedding space to bind them all,” 2023.
  * [17] Alexander Kirillov, Eric Mintun et al.,  “Segment anything,” 2023.
  * [18] Hangbo Bao, Li Dong et al.,  “Beit: Bert pre-training of image transformers,” 2022.
  * [19] Cancer Genome Atlas Networket al.,  “Comprehensive molecular characterization of human colon and rectal cancer,”  Nature, vol. 487, no. 7407, pp. 330, 2012.
  * [20] Nathan J. Edwards, Mauricio Oberti et al.,  “The cptac data portal: A resource for cancer proteomics research,”  Journal of Proteome Research, vol. 14, no. 6, pp. 2707–2713, 2015,  PMID: 25873244.
  * [21] Sophia J. Wagner, Daniel Reisenbüchler et al.,  “Transformer-based biomarker prediction from colorectal cancer histology: A large-scale multicentric study,”  Cancer Cell, vol. 41, no. 9, pp. 1650–1661.e4, 2023.
  * [22] Jakob Nikolas Kather, Niels Halama et al.,  “100,000 histological images of human colorectal cancer and healthy tissue,” May 2018.
  * [23] Marc Macenko, Marc Niethammer et al.,  “A method for normalizing histology slides for quantitative analysis,”  in 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2009, pp. 1107–1110.
  * [24] Richard J. Chen, Tong Ding et al.,  “A general-purpose self-supervised model for computational pathology,” 2023.


# Assessment of low-altitude atmospheric turbulence models for aircraft
aeroelasticity

Georgios Deskos  Alfonso del Carre  Rafael Palacios<EMAIL_ADDRESS>Department of Aeronautics, Imperial College London, United Kingdom

###### Abstract

We investigate the dynamic response of flexible aircraft in low-altitude
atmospheric turbulence. To this end, three turbulence models of increasing
fidelity, namely, the one-dimensional von Kármán model, the two-dimensional
Kaimal model and full three-dimensional wind fields extracted from large-eddy
simulations (LES) are used to simulate ambient turbulence near the ground.
Load calculations and flight trajectory predictions are conducted for a
flexible high-aspect ratio aircraft, using a fully coupled nonlinear flight
dynamics/aeroelastic model, when it operates in background atmospheric
turbulence generated by the aforementioned models. Comparison of load
envelopes and spectral content, on vehicles of varying flexibility, shows
strong dependency between the selected turbulence model and aircraft
aeroelastic response (e.g. 58% difference in the predicted magnitude of the
wing root bending moment between LES and von Kármán models). This is mainly
due to the presence of large flow structures at low altitudes that have
comparable dimensions to the vehicle, and which despite the relatively small
wind speeds within the Earth boundary layer, result in overall high load
events. Results show that one-dimensional models that do not capture those
effects provide fairly non-conservative load estimates and are unsuitable for
very flexible airframe design.

###### keywords:

Atmospheric turbulence , gust loads , geometrically-nonlinear structures ,
dynamic aeroelasticity , flight mechanics , large-eddy simulation

††journal: Journal of Fluids and Structures

##  1 Introduction

Atmospheric turbulence and its effects on aircraft flight has been a major
concern of the aerospace sciences for more than a century (Hunsaker and
Wilson, 1917). The efforts to characterise the impact of atmospheric
turbulence on aircraft loads can be traced back to the early development and
consolidation of discrete gust and continuous turbulence models (Jones, 1940;
Miles, 1956; Zbrozek, 1965; Eichenbaum, 1971; Houbolt, 1973). Today,
airworthiness certification procedures (e.g. European Aviation Safety Agency
(2007)) have adapted guidelines and methodologies derived from both that early
research and the accumulated experience of decades of in-service flight data.
However, the very recent development of high-altitude pseudo satellites (HAPS)
for an increasingly wide range of applications in remote sensing and
communications (Tozer and Grace, 2001; Yang and Mohammed, 2010) is challenging
the conventional understanding of the interactions between (nonlinear)
aeroelasticity and atmospheric turbulence. The extreme operational demands and
requirements of HAPS have resulted in more slender and much lighter platforms
which inevitably renders them more flexible. This is true for example for
stratospheric solar high-aspect ratio aircraft (NASA Helios, Boeing Oddyseus,
Facebook Aquila, Airbus Zephyr, etc.) (Gonzalo et al., 2018). Poor
understanding of the low-altitude response of such very flexible aircraft
(VFA) has previously led to poorly controlled flight dynamics and even vehicle
loss (Noll et al., 2004). Gaining physical insight into the aeroelastic
interactions while quantifying the uncertainty of flight dynamics is
considered key in further developing these technologies.

VFA aeroelasticity is characterised by large (geometrically-nonlinear)
structural displacements and by strong interactions between the low-frequency
structural modes and the vehicle flight dynamics. Conventional aircraft, on
the other hand, present a certain degree of decoupling between the structure
and the flight dynamics modes. Consequently, while aeroelastic analysis
methods in the latter case are based on linearised frequency-domain solvers
with potentially some higher-order corrections, the linear assumption brings
limitations that can result in poor estimates in the response of VFA (Cesnik
et al., 2014). Afonso et al. (2017) provide a comprehensive survey on the
state of the art of VFA aeroelasticity with focus on commercial aviation.
Indeed, while the current HAPS configurations are already relatively flexible
and feature low-frequency structural modes, they are likely to be even more
flexible as technology matures.

Continuous turbulence design criteria are often based on the von Kármán
spectrum with few empirical observations. Employing such approach to study
atmospheric turbulence interacting with very flexible aircraft aeroelasticity
poses two major issues. First, the von Kármán spectrum is derived from
isotropic turbulence arguments, a condition which is rarely met in atmospheric
turbulence, and never at low altitudes. In the lower part of the planetary
boundary layer (surface layer), turbulence is strongly modulated by the
presence of an aerodynamically rough surface and heat flux due to solar
radiation. These effects combined together can significantly alter the flow
dynamics and therefore deviate from the isotropic turbulence case implied by
the von Kármán-like turbulence that is usually found in higher altitudes (e.g.
CAT). The low-altitude operating environment of a generic VFA is shown
schematically in figure 1.

Figure 1: Schematic representation of the atmospheric turbulence experienced
by the flexible aircraft.

The near-surface turbulence coherent structures (e.g. high/low-speed streaks)
will dominate the flight dynamics of such vehicle. The second shortcoming
relates to the coupling between these large-scale turbulence structures (e.g.
eddies) with the dynamics of very flexible aircraft. As high-wing deformations
are expected under normal trim conditions, aircraft are prone to dynamic
instabilities of significant magnitude (Afonso et al., 2017; Noll et al.,
2004). In particular, high-aspect ratio wings are expected to present a
coupling between rigid-body and elastic motions since the frequency of the
elastic modes are closer to those of the flight dynamics modes. This strong
coupling has been well-documented in, e.g., Schmidt (2012) and Van Schoor and
Von Flotow (1990) for general aircraft aeroelasticity and, more specifically
for high-aspect ratio aircraft, in Patil et al. (2001), Tang and Dowell
(2004), Raghavan and Patil (2009), and Stanford and Beran (2013), among
others. To this day, a limited number of experimental studies on high-aspect
ratio wings exists. They have shown however, that very flexible wings have a
significantly different gust response characteristics under different
deformations and load conditions (Liu et al., 2016). Inherently, additional
analysis is required in order to gain insight into the complex interactions
between wind gusts and aircraft aeroelasticity. This particular study focuses
on investigating the effect of low-altitude turbulence on very flexible high-
aspect ratio aircraft by first conducting a sensitivity analysis of
atmospheric turbulence generation methods including the industry-standard von
Kármán model as well as demonstrating the differences in the dynamic
aeroelastic response between very flexible and stiff aircraft.

To investigate the sensitivity of the different methods used to produce the
continuous turbulent field on flexible aircraft gust loads, we examine the
response of a representative very flexible aircraft model to three atmospheric
turbulence input methods namely, the one-dimensional von Kármán, a two-
dimensional stochastic formulation Kaimal model (Kaimal et al., 1972) enhanced
with a spatial coherence model as well as full-three dimensional wind fields
produced by means of large-eddy simulations (Porté-Agel et al., 2000; Bou-Zeid
et al., 2005; Deskos, 2019). The paper starts with a brief introduction to
atmospheric turbulence modelling techniques in Section 2. The setups used to
simulate the wind fields are presented in Section 3 together with comparing
their statistics and respective velocity spectra. Subsequently, we make use of
a fully-coupled computational aeroelasticity model to make time-domain
predictions of the flight dynamics and loads. The results are presented in
Section 4. Results from the simulations including a sensitivity analysis for
the three turbulence inputs and a cross-comparison between a flexible and a
stiffer aircraft are presented in Section 5. Finally, a summary together with
some key findings are presented in Section 6.

##  2 Atmospheric turbulence modelling

Within the atmospheric boundary layer (ABL), turbulence evolves in response to
the local topography as well as the heating and cooling of the Earth’s
surface. The latter, categorises the ABL in three distinct states, namely, a
neutral boundary layer (NBL), a convective boundary layer (CBL) and a stable
boundary layer (SBL). In the neutral boundary layer, thermal effects remain
unimportant and air parcels move adiabatically solely due to Coriolis forcing
and surface friction. Thermal effects on the other hand are significant in the
convective and stable boundary layers in which temperature stratification
leads to a capping inversion (inversion layer). In the CBL, capping inversion
acts as a lid and turbulence convection is dominated by large plumes
(thermals) that transport and mix air away from the surface. The CBL state is
often found following the sunrise. On the other hand, the SBL is characterised
by a strong temperature stratification which typically leads to a shallower
boundary layer thickness, strong wind shear and much smaller and weaker eddies
than the neutral and convective cases Wyngaard (2010). Nonetheless, all three
states (except for the convective case) are very rarely in a state of
equilibrium and transition from one to another takes place multiple times
during the day. However, provided the disparity between the ABL state
evolution ($\sim$ hours) and the time spent by a VFA in the lower ABL ($\sim$
minutes), the assumption of ABL equilibrium can be considered a valid
assumption.

In simulating atmospheric turbulence, two approaches can be generally
considered. The first approach uses stochastic formulations to generate wind
fields based on first and second-order statistics (spectral/co-spectra).
Stochastic generation of wind fields has been an active area of research over
the last decades (Veers, 1988; Mann, 1994, 1998; Muñoz-Esparza et al., 2015).
For applications in wind engineering (Jeong et al., 2014; García et al.,
2017), the most common stochastic model is that of Kaimal (Kaimal et al.,
1972; Veers, 1988). The second technique on the other hand, utilises detailed
numerical simulations to resolve the fluid flow equations (or a
parametrisation of them) for a given temporal and spatial scale. Such a
popular method is large-eddy simulation (LES). The idea underpinning LES lies
in resolving the large energetic structures of the flow while modelling the
smaller scales. LES produces three-dimensional wind field for a given
computational domain and its quality depends on the resolution and modelling
parameters. LES can be therefore considered as a high-fidelity model, as it
resolves all relevant turbulent scales. Further details on both approaches are
discussed next.

###  2.1 Stochastic formulations for turbulence generation

In analysing the dynamic response of an aircraft, three-dimensional isotropic
turbulence models are often used to generate the turbulence input. They
produce stochastic time signals for the three velocity components ($u,v,w$) by
using either the Dryden or the von Kármán spectrum. The three-component signal
is produced by passing band-limited white noise through forming filters
obtained from spectral factorisation of the semi-empirical models (Ly and
Chan, 1980; Campbell, 1986). Uniform signals in transversal planes are
subsequently generated in the area of interest creating a velocity field with
perfect streamwise coherence, $\text{coh}_{uu}=\,$$1$. In this study we use
the von Kármán model. For a nominal reference velocity $U_{ref}$, it is
defined by the power spectral density,

$S_{u}(f)=\frac{2\sigma_{u}^{2}L_{u}}{\pi
U_{\textnormal{ref}}}\cdot\frac{1}{\bigg{[}1+\bigg{(}1.339L_{u}\frac{2\pi
f}{U_{\textnormal{ref}}}\bigg{)}^{2}\bigg{]}^{5/6}}.$ (1)

for the longitudinal velocity and

$S_{v,w}(f)=\frac{2\sigma_{v,w}^{2}L_{v,w}}{\pi
U_{\textnormal{ref}}}\cdot\frac{1+\frac{8}{3}\bigg{(}2.678L_{v,w}\frac{2\pi
f}{U_{\textnormal{ref}}}\bigg{)}^{2}}{\bigg{[}1+\bigg{(}2.678L_{v,w}\frac{2\pi
f}{U_{\textnormal{ref}}}\bigg{)}^{2}\bigg{]}^{11/6}}.$ (2)

for the lateral velocities. As our interest is in the low-altitudes,
independent turbulence length scales $L_{u},L_{v},L_{w}$ are assumed in three
directions. They are chosen as in MIL-HDBK-1797 (1997) so that $2L_{w}=h$ and
$L_{u}=2L_{v}=\,h/(0.177+0.000823h)^{1.2}$, where $h$ is the altitude (in
feet). The corresponding turbulence intensities are
$\sigma_{u},\sigma_{v},\sigma_{w}$. For the simulations, a time-signal of
$250\text{\,}\mathrm{s}$ is generated and subsequently convected into a full
wind field at speed $U_{ref}$ by means of the “frozen turbulence” hypothesis.

A more accurate representation of low-altitude atmospheric turbulence can also
be sought via the Kaimal model. The Kaimal model can be used to make
predictions of the wind field during the neutral boundary layer state using
empirically calibrated constants. In the present study, Kaimal is used in
conjunction with a spatial coherence model for the streamwise component of the
velocity field, $u(y,z,t)$, known as the International Electrotechnical
Commission (IEC) exponential coherence model (IEC-61400, 2005),

$\text{coh}_{uu}(f,r,U_{\textnormal{ref}})=\exp\bigg{[}-a\sqrt{\bigg{(}\frac{rf}{U_{\textnormal{ref}}}\bigg{)}^{2}+\bigg{(}\frac{rb}{L_{u}}\bigg{)}^{2}}\bigg{]}$
(3)

where $a=\,$$10$ and $b=\,$$0.12$ are the coherence decrement and offset
parameters, $L_{u}$ the streamwise turbulence lengthscale,
$U_{\textnormal{ref}}$ is the mean wind speed at the reference height, $f$
denotes frequency and $r$ is a radial distance in a plane vertical to the
streamwise direction. To generate the 2D wind field we make use of the open-
source code TurbSim (Jonkman and Kilcher, 2012). As TurbSim outputs planar
fields of the velocity field as a function of time we use again Taylor’s
frozen turbulence hypothesis, $u(x,t)=u(0,x-tU)$, to generate a full three-
dimensional velocity field.

###  2.2 High-fidelity turbulent flows via large-eddy simulation

Large-eddy simulations (LES) of atmospheric turbulence are conducted using our
in-house solver WInc3D (Deskos et al., 2018; Deskos, 2019). It is an explicit
LES framework which solves the unsteady, incompressible, filtered Navier-
Stokes equations using their skew-symmetric form

$\frac{\partial\widetilde{u}_{i}}{\partial
t}+\frac{1}{2}\bigg{(}\widetilde{u}_{j}\frac{\partial\widetilde{u}_{i}}{\partial
x_{j}}+\frac{\partial\widetilde{u}_{i}\widetilde{u}_{j}}{\partial
x_{j}}\bigg{)}=-\frac{1}{\rho}\frac{\partial\widetilde{p^{*}}}{\partial
x_{i}}-\frac{1}{\rho}\frac{\partial p_{\infty}}{\partial
x_{i}}-\frac{\partial\tau_{ij}}{\partial x_{j}},\quad i,j=\\{1,2,3\\},$ (4)
$\frac{\partial\widetilde{u}_{i}}{\partial x_{i}}=0,$ (5)

where $\widetilde{p^{*}}=\widetilde{p}+1/3\widetilde{u_{i}}\widetilde{u_{i}}$,
and $\widetilde{u}_{i}$ are the filtered components of the modified pressure
and velocity fields, respectively, and $\rho$ is the fluid density. The
indices $(1,2,3)$ correspond to a $(x,y,z)$ coordinate system. The sub-filter
stresses $-\partial_{j}\tau_{ij}$ (residual terms from the filtering process)
are calculated using the standard Smagorinsky model (Smagorinsky, 1963)

$\tau_{ij}=-2(C_{S}\Delta)^{2}|\widetilde{S}|\widetilde{S}_{ij},\quad\widetilde{S}_{ij}=\frac{1}{2}\bigg{(}\frac{\partial\widetilde{u}_{i}}{\partial
x_{j}}+\frac{\partial\widetilde{u}_{j}}{\partial x_{i}}\bigg{)},$ (6)

where $C_{S}$ is the Smagorinsky constant which is corrected near the wall
using the Mason and Thomson (1992)

$C_{S}=\bigg{(}C_{0}^{n}+\bigg{\\{}\kappa\bigg{(}\frac{z}{\Delta}+\frac{z_{0}}{\Delta}\bigg{)}\bigg{\\}}^{-n}\bigg{)}^{-1/n}.$
(7)

where $z_{0}$ is the roughness lengthscale, $\Delta=\sqrt[3]{\Delta x\Delta
y\Delta z}$ the grid size and choosing $(C_{0},n)=\,(0.14,3)$. Finally,
$-\partial_{i}p_{\infty}$ is a constant pressure gradient applied to flow to
reproduce the effect of geostrophic forcing. It is worth noting that in
equation (4) the viscous term of the Navier-Stokes equation has been
neglected. This is due to the high-Reynolds number and the fact that the mesh
is coarse enough so that the near-ground viscous layer will not have any
impact on the calculations. Another implication of being under-resolved in the
near-ground surface is that applying a no-slip boundary condition is devoid of
any meaning. Instead, we simulate the bulk effect of velocity shear by
replacing it with slip velocity and imposing a local shear stress through a
wall-stress model. Commonly, the wall-stress model of Moeng (1984) is used by
considering the Monin-Obukhov similarity (MOS) theory and in this study we
employ the recent formulations of Bou-Zeid et al. (2005),

$\tau^{\text{wall}}(x,y)=\tau_{w}(x,y)\frac{\widehat{\widetilde{u}}_{i}(x,y,\Delta
z/2)}{\sqrt{\widehat{\widetilde{u}}_{x}^{\,2}(x,y,\Delta
z/2)+\widehat{\widetilde{u}}_{y}^{\,2}(x,y,\Delta z/2)}},$ (8)

and

$\tau_{w}(x,y)=-\left[\frac{\kappa}{\ln\left(\frac{\Delta
z/2}{z_{0}}\right)}\right]^{2}\left[\widehat{\widetilde{u}}_{x}^{\,2}(x,y,\Delta
z/2)+\widehat{\widetilde{u}}_{y}^{\,2}(x,y,\Delta z/2)\right].$ (9)

in which the horizontal-averaged boundary condition is replaced by a twice-
filtered velocity. In addition, we make use of periodic boundary conditions in
lateral directions and free/partial-slip
($\partial_{3}\widetilde{u}_{i}=\widetilde{u}_{3}=0$ for $i\in{1,2}$) for the
bottom $z=0$, and the top boundary $z=H$, respectively. To solve the governing
equations, sixth-order compact finite-difference schemes are used (Laizet and
Lamballais, 2009) on a Cartesian mesh in a half-staggered arrangement (the
same mesh is used for the three velocity components
$(u_{1},u_{2},u_{3})=\,(u,v,w)$, with a different mesh used for pressure $p$),
pressure correction is achieved by a direct spectral solver, while for the
time integration a low-storage third-order Runge-Kutta scheme is employed.
Additionally, explicit filtering is also applied to all fields and at each
time step using discrete sixth-order filters (Gaitonde and Visbal, 1998).
Finally, massive parallelisation of the numerical solver is achieved using MPI
and an efficient 2D pencil domain decomposition approach (Laizet and Li,
2011). The current implementation has been run on over $\mathcal{O}(10^{5})$
number of processors and in multiple platforms and has consistently shown
excellent scaling properties.

##  3 Comparing turbulence generation methods

For a robust comparison between the different synthetic turbulence models (von
Kármán, Kaimal, and LES), we must first ensure that they all yield
statistically similar results. This includes both the time-averaged velocity
profiles as well as velocity spectra. Here, the neutral atmospheric boundary
layer on a flat Earth is considered with three levels of surface roughness,
$z_{0}=\,$$0.0001\text{\,}\mathrm{m}0.001\text{\,}\mathrm{m}0.01\text{\,}\mathrm{m}$.
This defines three flow conditions that will be referred to as _Mild_ ,
_Regular_ and _Severe_ , respectively, and correspond to surface friction
velocities
$u_{*}=\,$$0.174\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}0.235\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}0.274\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$.
The friction velocities have been selected such that the mean velocity at
altitude $10\text{\,}\mathrm{m}$ attains a value of approximately
$U_{10}\approx\,$$5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ in all
cases. Further details are provided in the table 1 below.

Table 1:  Atmospheric turbulence flow parameters Name |  $z_{0}$ [m] |  $u_{*}$ [m/s] |  $U_{10}$ [m/s] |  $\mathcal{I}$ %
---|---|---|---|---
Mild | 0.0001 | 0.174 | 5.05 | 8.87
Regular | 0.001 | 0.235 | 5.41 | 10.74
Severe | 0.01 | 0.274 | 4.74 | 12.73

Before proceeding to presenting the results it is also instructive to report
the simulation configurations and some key parameters used for generating the
data. Starting with large-eddy simulations, we consider a region of
$L_{x}\>\times\>L_{y}\>\times\>L_{z}=$1000\text{\,}\mathrm{m}$\>\times\>$500\text{\,}\mathrm{m}$\>\times\>$500\text{\,}\mathrm{m}$$
discretised by $N_{x}\times N_{y}\times N_{z}=\,1024\times 512\times 513$ mesh
nodes. The time step was selected to be equal to $\Delta
t=\,$$0.05\text{\,}\mathrm{s}$ so that a Courant-Friedrichs-Lewy (CFL) number
of approximately $0.25$ is maintained throughout. Each simulation was run for
50 large-eddy turnover times ($\delta/u_{*}$, where
$\delta=\,$$1000\text{\,}\mathrm{m}$ is the boundary layer height) a spin-up
time that has been previously used to obtain converged statistics (Bou-Zeid et
al., 2005; Xie et al., 2015). Restricted by the CFL number, the total number
of time steps required during the precursor simulation amounts to $2$ million.
All three cases were run using $4096$ processors and required approximately
$1.5$ million CPU hours each. For the simulations involving the aeroelastic
coupling, a small number of snapshots were extracted (every 3 seconds) and
interpolated to obtain wind fields for a total of 10 minutes. For the Kaimal
model a smaller mesh of $N_{y}\>\times\>N_{z}=\,$$32$ $\times$ $32$ nodes was
used for a domain of
$L_{y}\>\times\>L_{z}=$50\text{\,}\mathrm{m}$\>\times\>$50\text{\,}\mathrm{m}$$,
obtaining a resolution of approximately $1.5\text{\,}\mathrm{m}$. The
calculations were made for a total analysis time of $250\text{\,}\mathrm{s}$
using a time step $\Delta t=\,$$0.25\text{\,}\mathrm{s}$. To reproduce the
desired results TurbSim (Jonkman and Kilcher, 2012) was used in a customised
mode by applying a logarithmic profile and setting the turbulence level at the
characteristic altitude, surface frictions and coherence parameters so that we
closely match the LES results. Lastly, for the von Kármán model, the standard
deviations of the three velocity components $\sigma_{u}$, $\sigma_{v}$ and
$\sigma_{w}$ together with the length scale $L_{u}$ extracted from the LES
simulations were used to generate zero-mean time series before added to the
background characteristic velocity $U_{10}$. A comparison of the vertical
profiles and power spectral densities (PSD) are shown in figure 2 and 3,
respectively. The vertical mean velocity profiles agree well for the Kaimal
and LES data whereas the von Kármán can only match the velocity at the
reference altitude.

Figure 2: Mean velocity and turbulence intensity profiles for the three cases

Larger discrepancies can be found between the the models in the streamwise and
transverse velocity variances. In particular, LES predicts a linear decrease
of the streamwise variance, a result that agrees with the observed data as
well as theoretical scaling laws (Wyngaard, 2010), whereas the Kaimal and von
Kármán models yield uniform distributions for all components.

Figure 3: Power spectral density of velocity fluctuations for each velocity
component ($u^{\prime}$, $v^{\prime}$ and $w^{\prime}$) and all three models
(1D von Kármán, 2D Kaimal and 3D LES).

Similarly, PSDs from the three models, shown in figure 3, appear to collapse
together for the streamwise component but exhibit larger discrepancies in the
lateral ones. Specifically, in the vertical velocity component PSDs, LES data
plots significantly deviate from those of Kaimal and von Kármán. Those two
models assume flow isotropy which results in nearly identical PSDs over the
$v^{\prime}$ and $w^{\prime}$ components.

Figure 4: Instantaneous snapshots of a _mild_ streamwise velocity, $u(t)$, for
the Kaimal and LES models.

Nevertheless, while our comparison of the flow characteristics is ended here
by comparing at most second-order statistics (velocity variance, spectra) more
differences between the LES, which produces a non-gaussian distribution of
velocity, and the stochastically-generated turbulence can be further sought.
This is for example demonstrated qualitatively in figure 4. The two snapshots
of the LES and Kaimal, respectively, are both taken from the _Mild_ case and
exhibit large differences in the formation of the near-wall large coherent
structures. LES flow structures consist of large contra-rotating streamtubes
generating large structures which extend as high as $25\text{\,}\mathrm{m}$
while Kaimal consists of smaller structures extending only a few metres along
the y-direction. The von Kármán model on the other hand is not plotted as only
a uniform distribution will be assumed in a 2D planar snapshot. These
qualitatively observed differences between the three models which are not
captured by either the flow statistics will be re-assessed later on during the
aeroelastic simulations in Section 5.

##  4 Nonlinear time-domain aeroelastic solver

The aeroelastic model consists of a structural solver based on a non-linear
geometrically exact beam theory (GEBT) and an aerodynamic model based on the
unsteady vortex lattice method (UVLM), as originally proposed by Murua et al.
(2012). The 1D structure is closely coupled to convergence at each FSI time
step to the discrete-time UVLM using a coinciding spanwise panel
discretisation. In particular, quadratic finite elements are used for the
GEBT, with a discretisation based on nodal displacements and rotations, which
are given in a body-attached, moving frame of reference $A$ (this frame moves
and rotates with respect to an inertial frame, $G$) as shown in figure 5.

Figure 5: Frames of reference used in the structural solver and
parametrisation of the instantaneous geometry.

The beam dynamics description is derived from the equations of a curved,
geometrically nonlinear beam using the local element frame of reference $B$.
Subsequently, the discrete system of equation (after applying a quadratic
finite elements approximation) is expressed in terms of $\bm{\eta}$, the state
variable containing the nodal displacements and rotations, and
$[\bm{v}_{A},\bm{\omega}_{A}]$, the state variable of velocities and rotations
of frame $A$, as (Hesse and Palacios, 2016)

$\bm{\mathcal{M}}(\bm{\eta})\left\\{\begin{array}[]{c}\dot{\bm{v}}_{A}\\\
\dot{\bm{\omega}}_{A}\end{array}\right\\}+\bm{\mathcal{Q}}_{\text{gyr}}(\bm{\eta},\dot{\bm{\eta}},\bm{v}_{A},\bm{\omega}_{A})+\bm{\mathcal{Q}}_{\text{stiff}}(\bm{\eta})=\bm{\mathcal{Q}}_{\text{ext}}$
(10)

where $\bm{\mathcal{M}}$ is the mass matrix and
$\bm{\mathcal{Q}}_{\text{gyr}}$, $\bm{\mathcal{Q}}_{\text{stiff}}$ and
$\bm{\mathcal{Q}}_{\text{ext}}$ are the discrete gyroscopic, stiffness, and
external generalised forces respectively. This equation is then solved
iteratively using a fixed-point iteration scheme while time integration is
carried out with an explicit Newmark-$\beta$ scheme formulated in incremental
form. The applied forces are given in the material frame of reference ($B$)
allowing for follower force input.

The aerodynamic solver is based on the unsteady vortex lattice method (UVLM),
a three-dimensional nonlinear potential computational technique for lifting
surfaces. To this end, vortex rings are placed on the wings of the aircraft
and by applying the no-penetration condition, Biot-Savart’s law and the Kutta-
Joukowski condition we are able to calculate the bound circulation
$\bf{\Gamma}$ vector used to calculate the force increment which is split into
steady and unsteady (apparent mass) contributions. These two quantities can be
integrated to yield the aerodynamic forces at the centre point of the segments

$\partial\bm{F}_{\text{steady}}=\rho_{\infty}||\bf{\Gamma}||(\bm{v}\times\partial\bm{l})$
(11)

where $k\in[0,4)$ is the segment counter and $\bm{v}$, $\partial\bm{l}$ the
structure segment’s velocity and length. The unsteady force at the centre
point of the panel is calculated via,

$\partial\bm{F}_{\text{unsteady}}=\rho_{\infty}\dot{||\bf{\Gamma}||}(t)(\hat{\bm{v}}\times\partial\bm{l})$
(12)

where $\hat{\bm{v}}$ is a panel-weighted velocity and
$\dot{||\bf{\Gamma}||}(t)$ is the rate of change of the magnitude of bound
vorticity with respect to an inertial observer. To model the wake of the
aircraft, a background-flow-convected wake model is used. In this way, we do
not account for the mutual interactions between the free (shed) vortex rings,
which significantly reduces the computational cost during each iteration (Del
Carre and Palacios, 2019). Instead, the wake is allowed to convect by the
ambient flow and interact with the structure without exhibiting the near tip
roll-up. Studies conducted by the authors have determined that an ambient
flow-convected wake yields near identical results when compared to the free-
force wake for the problems of interest, while the savings in computational
cost are remarkable. The present mid-fidelity aeroelastic solver deals with
the coupled behaviour of the structural/aerodynamic effects of flexible
structures and it has been implemented within the framework SHARPy, written in
python, with low level libraries in C++ and Fortran, and parallelisation using
OpenMP.

###  4.1 3D velocity fields as aeroelastic inputs

A final aspect of the present model is the methodology for combining the
turbulence input data with the dynamic aeroelastic solver. Atmospheric
turbulence fields either in a two-dimensional form (vertical planes) or three-
dimensional sub-domains of a flow realisation where generated for a time
period of $120\text{\,}\mathrm{s}$. Subsequently, the data was stored and read
by SHARPy at each time step. The LES input data was found to impose practical
constraints in the reading and interpolation process as each snapshot amounts
to more than 3 Gb of binary memory (for $1024\times 512\times
513\sim$270\text{\,}\,$$million grid points). This poses a bottleneck which
may be approached in two ways. A direct solution is to read and store in RAM
only the two necessary snapshots. This is suitable for distributed computers,
where permanent storage latency is high, and nodes can be devoted to a small
number of simulations each. However, shared memory computers usually have a
lower permanent storage latency and the RAM available per processor node is
lower. The solution is then to map the binary file stored in the relatively
fast hard drive into an array-like structure, so data is directly read from
the permanent storage. This solution was seen to drastically reduce the memory
overhead with a 20% increase of the interpolation time cost. The 3D velocity
field reader implemented here supports both methods, and it implements a cache
where data fields are read only once during each simulation. In the end, the
decision of which method is used for the 3D field read is a combination of the
architecture the software is being run in, the number of cases running
simultaneously and the time between LES snapshots. Once the 3D field is mapped
or copied into an array, a multivariable interpolation is conducted. First,
the volumetric data is given as a structured regular grid. This makes
trilinear interpolation a suitable solution, while the quicker closest
neighbour interpolation is also available for finer 3D field discretisations.
Results presented here are obtained using a trilinear interpolation in space.
The time-dependency of the 3D velocity fields is considered in this analysis.
When multiple snapshots of the field are given, a linear interpolation is
performed between the two closest snapshots to the simulation time. A
visualisation of the combined 3D velocity fields and aircraft mode is shown in
figure 6 below.

Figure 6: Visualisation of the 3D LES-generated velocity field (severe case)
and its combination with the aircraft model are shown.

##  5 Numerical results

As discussed in Section 1, in the design process of aircraft the 1D von Kármán
spectrum is typically used to simulate continuous gusts. Having already
presented the differences between the three turbulence-generation models (1D
von Kármán, 2D Kaimal and 3D LES), in this section we shall attempt to
quantify their effects on aircraft aeroelasticity from short simulated flights
in open loop (that is, without feedback control on the vehicle dynamics) at
the reference height $H_{\textnormal{ref}}=\,$$10\text{\,}\mathrm{m}$.

###  5.1 Test case definition

Simulations will be carried out on the High-Altitude Long-Endurance (HALE)
T-tail configuration of Murua et al. (2012), with minor modifications. This
test case provides a simple, easily reproducible and representative HALE
aircraft example, which has been previously used for aeroelastic simulations
by Murua et al. (2012) and Hesse and Palacios (2016), among others. The
vehicle features a high-aspect ratio wing with relatively low stiffness in
torsion and out-of-plane bending and a conventional T-tail configuration, with
elevator and rudder as control surfaces. All aircraft structures are flexible,
although fuselage and tail are significantly stiffer than the wing. Finally,
all aerodynamic surfaces of the aircraft are represented by a symmetric
aerofoils. Two aircraft configurations are considered, which will be referred
to as the _stiff_ and _flexible_ vehicles. They exhibit exactly the same
geometry and inertial characteristics with the former stiffness properties
scaled by a factor of 100. Details are presented in table 2 for both
configurations. The only active control surface in the simulations below is
the elevator, and once the aircraft is trimmed for cruise flight, it is left
fixed at that deflection.

Figure 7: Schematic representation of the T-tail aircraft geometry. Table 2:
Stiffness and mass parameters for the aircraft model.

Component |  $EA$ $\left[\text{N}\cdot\text{m}\right]$ |  $GJ$ $[\text{N}\cdot\text{m}^{2}]$ |  $EI_{y}$ $GJ$ $[\text{N}\cdot\text{m}^{2}]$ |  $EI_{z}$ $[\text{N}\cdot\text{m}^{2}]$ |  $\bar{m}$ [kg$\cdot$m-1] |  $\bar{J}$ [kg$\cdot$m]
---|---|---|---|---|---|---
|  |  | Flexible Aircraft |  |  |
Wing | 1.5e4 | 1.5e4 | 6.0e4 | 4.2e6 | 0.75 | 0.1
Fuselage | 1.5e5 | 1.5e5 | 6.0e5 | 4.2e7 | 0.1 | 0.01
Tail | 1.5e5 | 1.5e5 | 6.0e5 | 4.2e7 | 0.1 | 0.01
|  |  | Stiff Aircraft |  |  |
Wing | 1.0e6 | 1.0e6 | 2.0e6 | 1.4e8 | 0.75 | 0.1
Fuselage | 1.0e7 | 1.0e7 | 2.0e7 | 1.4e9 | 0.1 | 0.01
Tail | 1.0e7 | 1.0e7 | 2.0e7 | 1.4e9 | 0.1 | 0.01

###  5.2 Flexible aircraft dynamics for each atmospheric turbulence model

To examine the relative importance of input turbulence data in predicting both
loads and flight trajectories for a flexible aircraft we undertake a number of
simulations for all three turbulence intensity (TI) cases. Three-dimensional
atmospheric turbulence fields are extracted from the three models such that
they exhibit similar flow statistics, as discussed in section 3, with
simulations corresponding to $120\text{\,}\mathrm{s}$ of real time and $\Delta
t=\,$$0.025\text{\,}\mathrm{s}$. The produced time-series for the aircraft
loads and trajectories for each TI and turbulence-generation model are
subsequently used in our assessments. For the LES-generated background flows,
multiple aeroelastic simulations were run for each turbulent intensity case in
order to explore the spanwise flow variations of the simulated atmospheric
turbulence. A summary of all aeroelastic simulations is shown in table 3
below. Each case is defined by the average reference velocity and turbulence
intensity as well as the initial spanwise position in the turbulent flow
field, $Y$, of the aircraft mid-wing point. Values of
$Y=\,$$0\text{\,}\mathrm{m}10\text{\,}\mathrm{m}50\text{\,}\mathrm{m}200\text{\,}\mathrm{m}-200\text{\,}\mathrm{m}$
are selected based on the lateral turbulence integral length scale, which is
calculated to be between $150\text{\,}\mathrm{m}200\text{\,}\mathrm{m}$ for
the three turbulence intensity cases.

Table 3:  Summary of the various simulations along with some key parameters Case |  $U_{10}$ (m/s) |  $\mathcal{I}_{10}$ (%) | Y$(t=0)$ (m) | Turbulence model
---|---|---|---|---
L11 | 5.05 | 8.87 | 0 | LES
L12 | 5.05 | 8.87 | 10 | LES
L13 | 5.05 | 8.87 | 50 | LES
L14 | 5.05 | 8.87 | 200 | LES
L15 | 5.05 | 8.87 | -200 | LES
L21 | 5.41 | 10.74 | 0 | LES
L22 | 5.41 | 10.74 | 10 | LES
L23 | 5.41 | 10.74 | 50 | LES
L24 | 5.41 | 10.74 | 200 | LES
L25 | 5.41 | 10.74 | -200 | LES
L31 | 4.74 | 12.73 | 0 | LES
L32 | 4.74 | 12.73 | 10 | LES
L33 | 4.74 | 12.73 | 50 | LES
L34 | 4.74 | 12.73 | 200 | LES
L35 | 4.74 | 12.73 | -200 | LES
K1 | 5.05 | 8.87 | - | Kaimal
K2 | 5.41 | 10.74 | - | Kaimal
K3 | 4.74 | 12.73 | - | Kaimal
VK1 | 5.05 | 8.87 | - | von Kármán
VK2 | 5.41 | 10.74 | - | von Kármán
VK3 | 4.74 | 12.73 | - | von Kármán

Aircraft trajectories are tracked at the origin of the vehicle axes shown in
figure 7 and the trajectories for all simulations in table 3 are shown in
figure 8. We may observe that the aircraft exhibits a similar vertical
trajectory for all turbulence inputs by maintaining an altitude around
$10\text{\,}\mathrm{m}$ with “dives” and “jumps” of around
$1\text{\,}\mathrm{m}3\text{\,}\mathrm{m}$ occurring throughout. A significant
deviation from this pattern is observed in the response to the von Kármán
model for both the regular and severe cases. Perfect spanwise coherence in the
background flow implies generation of net lift resultant that results in
substantial inertia relief and an overall gain of altitude. The transverse
component of the trajectory exhibits larger deviations between the three
turbulence-generation models as well as within LES itself when the aircraft is
initially placed at different spanwise locations with respect to the same
background flow
($Y=\,$$0\text{\,}\mathrm{m}10\text{\,}\mathrm{m}50\text{\,}\mathrm{m}200\text{\,}\mathrm{m}-200\text{\,}\mathrm{m}$).
Interactions with the large scale flow features in the LES-generated
turbulence results in significant variation of the aircraft trajectory, which
can deviate both right and left while in some cases it may experience a full
U-turn and continue flying in the opposite direction. The other two models
(Kaimal and von Kármán), which are characterised by less or no variability in
the transverse direction, result in nearly straight flight trajectories.

Figure 8: Flight trajectory of the flexible configuration and all three
turbulence generation methods

To demonstrate the effect of spanwise velocity coherence experienced by the
aircraft wing as well as its relationship with the aerodynamic behaviour (e.g.
spanwise lift distribution) and aircraft trim shape, figure 9 presents
ensembled-averaged spanwsie velocity and lift distributions for the aircraft
wing from an LES-generated turbulence case. This corresponds to the regular TI
and the aircraft initiated at $Y=\,$$50\text{\,}\mathrm{m}$ spanwise distance
from the reference location (centre of the LES computational domain). For this
case, a U-turn of the aircraft occurs during the simulation.

Figure 9: Flight trajectory and ensemble-averaged lift/velocity distribution
of the flexible configuration for the regular LES case

Clearly, as the aircraft comes across large coherent structures it experiences
an incident velocity with large variations across its span (see figure 9). In
turn, the non-uniform velocity distribution gives rise to asymmetric loads
which through the aircrafts flexibility lead to a gradual veer to the right,
which leads to a vehicle U-turn. It should be noted that the aircraft altitude
does not change significantly during the U-turn trajectory and that ensemble
averaged lift distributions reduced further this variability. This is a result
of a passive control mechanism inherent to the aircraft by its design which
makes the aircraft navigate through turbulence without having full control
over its trajectory. This of course can be addressed by a flight control
system which can alter the aircraft trim and thus have an impact on loads as
well. As far as the root loads are concerned, similar variability is observed
for all three atmospheric turbulence inputs. In figure 10 we observe that
torsion fluctuates around a mean value of
$3.5\text{\,}\mathrm{N}\text{\,}\mathrm{m}$ for all three models with larger
amplitudes found for the regular and severe Kaimal cases. Out-of-plane (OOP)
and in-plane (IP) wing root moments, on the other hand, attain a mean value of
around $-1700\text{\,}\mathrm{N}\mathrm{m}$ and
$370\text{\,}\mathrm{N}\mathrm{m}$, respectively, with the fluctuations under
the Kaimal and LES inputs being more pronounced in the time series. It is
worth mentioning here that the aeroelastic simulations which make use of the
Kaimal model exhibit relatively high-frequency fluctuations. These
fluctuations which appear to dominate the loads after about
$12.5\text{\,}\mathrm{s}$ coincide with the instance that the aircraft’s
trajectory deviates from a straight path (the one aligned with the mean flow).
This effect can be attributed to the lack of streamwise coherence of the
Kaimal model as well as the assumption of “frozen turbulence” which generates
non-physical longitudinal structures once the 2D planes are projected to a
three-dimensional wind field. Once the aircraft is misaligned from the mean
flow, it is incident to a non-uniform velocity field with large velocity
differences across its wing span, increasing thus the torsion experienced by
the root. Similar fluctuations can also be observed for the “severe” case and
the LES inputs. This however is not an artefact of the LES method but rather
the effect of the near-surface large coherent flow structures.

Figure 10: Root loads from the three turbulence generation models for the
regular case. In the figure the torsion, out-of-plane (OOP) and in-plane (IP)
moments are shown.

Nonetheless, in order to fully appreciate the magnitude of the loads it is
worth plotting the load envelops for all three TI cases and root loads of
interest, namely torsion, OOP and IP.

Figure 11: Root loads envelopes for all three models and cases. In all figures
the out-of-plane (OOP) root moment is plotted against torsion.

The load envelops which are constructed by the convex hull obtained from the
maximum instantaneous loads exerted on the aircraft, provide information
regarding the loading regime expected for each TI case. To this end, from both
figures 11 and 12 we may observe the gradual expansion of the load envelopes
as the input turbulence intensity increases. The shape of the load envelopes
is also consistent in all three turbulence intensity cases. In the case of the
OOP plotted against torsion, the load envelopes are stretched along the OOP
direction signifying the significant increase of OOP with increasing levels of
turbulence. Similarly, the shape of the IP versus torsion envelopes appears to
be elongated with the larger in-plane moments occurring at negative root
torsion. Nonetheless, in both cases (OOP and IP) the choice of the turbulence
input model appears to be more significant. The LES data provide by far the
largest load envelope with its boundaries twice as big as those of Kaimal. On
the other hand, the Kaimal model yields smaller load envelopes with the
difference between the two models being increased with the magnitude of
turbulence intensity. Surprisingly, the smallest load envelopes are formed by
the von Kármán model input. This is true for both the IP and OOP moments. We
argue that this difference is due to the fact that in the von Kármán model,
the spanwise coherence is by definition equal to unity, leading to a uniform
spanwise load distribution. As a result, the dynamic loads give rise to a
lifting of the overall structure through “inertial relief” (applied work by
gust forces is transformed into rigid body kinetic energy). This is confirmed
by figure 8 in all three turbulence intensity cases by the vertical trajectory
of the von Kármán model.

Figure 12: Wing root loads envelopes for all three models obtained from
combining all available time histories. In all figures the in-plane (IP) root
moment is plotted against torsion.

Another important difference that is captured by the load envelopes is the
loads’ frequency of occurrence. We may observe that the boundaries of all
Kaimal and LES load envelopes are due to intermittent wind gusts which
increase with increasing levels of turbulence. To determine the potential
effects of intermittency on the aircraft loads, we have computed the
probability distribution function (PDF) of the incident (background) velocity
and load increments ($\Delta u=u(t+\Delta t)-u(t)$, $\Delta M=M(t+\Delta
t)-M(t)$) and examined the characteristics of its “tail”. Intermittency in
both the incident velocity and the loads can be captured by heavy tail
distributions that significantly deviate form the standard Gaussian
distribution.

Figure 13: PDF of the incident velocity magnitude$\Delta|U_{\infty}|$ and out-
of-plane moment $M_{OOP}$ for all models and turbulence intensity cases.
Normal distribution is included for comparison.

Figure 13 illustrates the existence of both intermittent velocity difference
and loads, particularly for the “severe” case in which a clear heavy tail
distributions, $|\Delta
U_{\infty}|>\,$$0.1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$ and $|\Delta
M_{OOP}|>\,$$7.5\text{\,}\mathrm{N}\text{\,}\mathrm{m}$ in both at positive
and negative. It is worth noting that in the incident velocity difference PDFs
the Kaimal model appears to be more intermittent. This is again due to the
frozen turbulence assumption in combination with the deviation of the aircraft
from a straight path. The wind field intermittency also translates into the
observed loads. More specifically, the PDFs for $|\Delta M_{OOP}|$ show that
the aircraft experiences large variation of out-of-plane loads sporadically
with probability density of less that $0.01$, particularly in the “severe”
case. In the other two “milder” cases the two models (Kaimal and LES) appear
to be in a better agreement. Von Kármán turbulence is characterised by
significantly less variability and thus, its PDF appears much narrower.

Figure 14: PSDs of the root loads (Torsion, OOP and IP) for all three
turbulence input cases (LES, Kaimal and von Kármán).

Comparing the OOP moment distribution with the other two models as well as
with that of the incident velocity for the von Kármán model we may not be able
to infer the origin of such a large difference. However, following our earlier
argument we may again justify such a large difference in the PDF occurs by the
stronger inertia relief in this case. Finally, it is worth presenting the
power spectral density (PSD) functions of the root loads from all the three
models in figure 14. Again, comparing the PSDs for the three models a similar
behaviour can be obtained. For instance, the PSD of torsion is characterised
by a plateau region over of low-frequency region, followed by a steep slope
and low-energy content $\sim$$0.001$ in the high frequencies. In addition, all
three cases exhibit a peak around $f=\,$$1.3\text{\,}\mathrm{Hz}$ which
corresponds to the mode II (natural torsion frequency). The peak around mode
II is more pronounced for the severe turbulence case while modes I
($f=\,$$0.4\text{\,}\mathrm{Hz}$) and V ($f=\,$$3.92\text{\,}\mathrm{Hz}$) are
less recognisable in the PSD plots. Nonetheless, the effect of increasing the
level of turbulence intensity on each PSD can be seen in all three cases by
“pushing” the distributions up. This effect is clearer over the lower
frequencies ($0.1\text{\,}\mathrm{Hz}1\text{\,}\mathrm{Hz}$) and in particular
for the simulations using the LES-generated turbulence input. A similar trend
is followed by the IP moment. On the other hand, the OOP follows a smoother
slope and lacks a “plateau” region, confirming that its response varies across
the range of scales. Nonetheless, it is again worth noting that again the von
Kármán simulations yield PSD distributions with an amplitude of about an order
of magnitude smaller than their respective counterparts generated by LES and
Kaimal turbulence inputs.

###  5.3 Aeroelastic effects assessment

Having demonstrated the differences between three turbulence input models (von
Kármán, Kaimal and LES), in this section we present a cross-comparison between
the _Flexible_ and the _Stiff_ aircraft configurations. For this analysis we
only use the LES turbulence inputs. The two aircraft configurations are
expected to exhibit a different response to dynamic gust loading. The flight
dynamics of the former (_Flexible_) aircraft configuration are characterised
by its coupling with the low-frequency structural modes Afonso et al. (2017).
On the other hand, the _Stiff_ aircraft configuration is expected to exhibit a
substantial degree of decoupling between the two. To address the question of
how much does stiffness affect the flight dynamics and structural loads of an
aircraft subject to atmospheric turbulence, we present and compare the load
envelopes and PSD distributions of the two configurations. To this end, the
two aircraft configurations have been simulated for the same number of cases
by maintaining all other parameters the same. We begin with the presentation
of the OOP vs torsion load envelope in figure 15.

Figure 15: Root loads envelopes for the “Stiff” and “Flexible” cases for all
levels of TI. In all figures the out-of-plane (OOP) root moment is plotted
against torsion.

Aircraft flexibility is found to fundamentally change the dynamic gust
response of the aircraft. The stiff aircraft is found to be prone to larger
OOP root moment loads and suffer less from large variations in torsion. More
specifically, OOP moment attains values ranging from
$-500\text{\,}\mathrm{N}\text{\,}\mathrm{m}$ to
$-2500\text{\,}\mathrm{N}\text{\,}\mathrm{m}$ while experiencing only positive
torsion values throughout.

Figure 16: Comparison of PSD distributions for the root loads (Torsion, OOP
and IP) between the stiff and flexible aircraft.

This shift in the loads response is also reflected in the respective PSDs of
the loads. Figure 15 shows a comparison between the two cases for all computed
root loads. As in the load envelopes, torsion fluctuations of the stiff case
are an order-of-magnitude) smaller than those of the flexible one with the
differences being more pronounced in the low-frequency region. This is a well-
known steady aeroelastic effect that cascades in – lift twists the wing thus
reducing the bending moments and increases torsion. Conversely, the OOP PSD
shows a significant increase in the energy over the mid-range frequencies
(inertial subrange). Finally smaller difference can be observed for the IP
moment with the stiff case exhibiting smaller fluctuations. This transfer of
energy from torsion and IP moment fluctuations to OOP signifies the
fundamental differences in the turbulence/loads transfer function between the
flexible and the stiff aircraft. The former, transforms the OOP loads into
torsion which in turn changes its aerodynamic behaviour and therefore the
resultant root loads. This can be seen as a passive load alleviation mechanism
which targets a reduction of the OOP moment. The stiff aircraft on the other
hand, by experiencing a certain degree of decoupling between structural and
flight dynamics, remains more responsive to incident turbulence (the $-5/3$
scaling is maintained in the inertial sub-range) implying it may experience
larger fatigue loads during its operation.

##  6 Conclusions

In the present study we have conducted aircraft aeroelasticity simulations
using turbulence input fields which are assumed to be representative of their
operating environment. The aeroelastic coupling and dynamic loading on a very
flexible aircraft was found to be highly modulated by the coherent structures
and flow features of low-altitude atmospheric turbulence. This may be of
critical importance for some very light air vehicles currently under
consideration.

The choice of the turbulence generation method was found to have a significant
impact on both loads and flight dynamics of the aircraft. LES and Kaimal
models provide similar results with LES being able to predict larger load
envelopes. Kaimal was found to be problematic when it is used within our
moving frame of reference as the application of the frozen turbulence
hypothesis used to transform temporal 2D planes into spatial full three-
dimensional wind fields does not generate three-dimensional coherence in the
wind field. As a result, small deviations of the aircraft from the dominant
wind direction generated sudden “load jumps” which affected the load envelopes
in our simulations. LES, on the other hand, simulates a more realistic
turbulent environment with large scales extending in the streamwise direction
and generating flow coherence. The von Kármán model was the worst performing
method of those considered, particularly when used in a moving three-
dimensional frame of reference. The uniform planar velocity fields from this
approach result in spanwise loading which substnunderpredicts the magnitude of
the loads and overpredicts the energy transferred into rigid body motion.
Therefore its use for loads and aeroelastics at low altitudes is not advised.

The comparison between the flexible and the stiff cases has demonstrated the
extent of aeroelastic effects in the previous simulations. Analysis of the two
aircraft flight dynamics shows that stiffer aircraft are prone to
instabilities from high-frequency perturbations (gusts) leading to large OOP
moment loads. More flexible vehicles, on the other hand, have natural load
alleviation mechanisms through structural (elastic) and rigid body motion
coupling. Results have shown smaller OOP loads as the energy contained in gust
loading is distributed between rigid body motion, torsion, OOP and IP moments.
Overall, as the flight dynamics of high-aspect ratio flexible aircraft were
found to be fundamentally different from that of the stiff one, more emphasis
should be given to non-linear aeroelastic phenomena which might occur under
turbulent conditions and are currently not well-captured by existing
computational models.

## Acknowledgements

The authors would like to acknowledge funding from EPSRC (grant EP/R007470/1)
and Airbus Defence and Space.

## References

## References

  * Hunsaker and Wilson (1917) J. C. Hunsaker, E. B. Wilson, Report on Behavior of Aeroplanes in Gusts, Technical Report, NACA Technical Report 1, NACA-TR-1, 1917\.
  * Jones (1940) R. T. Jones, The unsteady lift of a wing of finite aspect ratio, Technical Report, National Advisory Committee for Aeronautics, 1940. NACA-TR-681.
  * Miles (1956) J. W. Miles,  The aerodynamic force on an airfoil in a moving gust,  Journal of the Aeronautical Sciences 23 (1956) 1044–1050.
  * Zbrozek (1965) J. K. Zbrozek,  Atmospheric gusts: Present state of the art and further research,  The Journal of the Royal Aeronautical Society 69 (1965) 27–45.
  * Eichenbaum (1971) F. D. Eichenbaum,  A general theory of aircraft response to three-dimensional turbulence,  Journal of Aircraft 8 (1971) 353–360.
  * Houbolt (1973) J. C. Houbolt,  Atmospheric Turbulence,  AIAA Journal 11 (1973) 421–437.
  * European Aviation Safety Agency (2007) European Aviation Safety Agency, Certification Specifications for Certification Specifications for Normal-Category Aeroplanes CS-23, Technical Report, European Aviation Safety Agency, 2007.
  * Tozer and Grace (2001) T. C. Tozer, D. Grace,  High-altitude platforms for wireless communications,  Electronics Communication Engineering Journal 13 (2001) 127–137.
  * Yang and Mohammed (2010) Z. Yang, A. Mohammed,  Wireless communications from high altitude platforms: Applications, deployment and development,  in: 2010 IEEE 12th International Conference on Communication Technology, pp. 1476–1479.
  * Gonzalo et al. (2018) J. Gonzalo, D. López, D. Domínguez, A. García, A. Escapa,  On the capabilities and limitations of high altitude pseudo-satellites,  Progress in Aerospace Sciences 98 (2018) 37 – 56.
  * Noll et al. (2004) T. E. Noll, J. M. Brown, M. E. Perez-Davis, S. D. Ishmael, G. C. Tiffany, M. Gaier, Investigation of the Helios Prototype Aircraft Mishap, Technical Report, NASA, 2004. Mishap Report.
  * Cesnik et al. (2014) C. E. S. Cesnik, R. Palacios, E. Reichenbach,  Re-examined Structural Design Procedures for Very Flexible Aircraft,  Journal of Aircraft 51 (2014) 1580–1591.
  * Afonso et al. (2017) F. Afonso, J. Vale, É. Oliveira, F. Lau, A. Suleman,  A review on non-linear aeroelasticity of high aspect-ratio wings,  Progress in Aerospace Sciences 89 (2017) 40–57.
  * Schmidt (2012) D. K. Schmidt, Modern Flight Dynamics, Mc-Graw Hill, New York, NY, USA, 2012\.
  * Van Schoor and Von Flotow (1990) M. C. Van Schoor, A. H. Von Flotow,  Aeroelastic characteristics of a highly flexible aircraft,  Journal of Aircraft 27 (1990) 901–908.
  * Patil et al. (2001) M. J. Patil, D. H. Hodges, C. E. S. Cesnik,  Limit-Cycle Oscillations of High-Aspect-Ratio Wings,  Journal of Fluids and Structures 15 (2001) 107–132.
  * Tang and Dowell (2004) D. Tang, E. Dowell,  Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings,  Journal of Fluids and Structures 19 (2004) 291 – 306.
  * Raghavan and Patil (2009) B. Raghavan, M. J. Patil,  Flight dynamics of high aspect-ratio flying wings: Effect of large trim deformation,  Journal of Aircraft 46 (2009) 1808–1812.
  * Stanford and Beran (2013) B. Stanford, P. Beran,  Direct flutter and limit cycle computations of highly flexible wings for efficient analysis and optimization,  Journal of Fluids and Structures 36 (2013) 111 – 123.
  * Liu et al. (2016) Y. Liu, C. Xie, C. Yang, J. Cheng,  Gust response analysis and wind tunnel test for a high-aspect ratio wing,  Chinese Journal of Aeronautics 29 (2016) 91 – 103.
  * Kaimal et al. (1972) J. C. Kaimal, J. C. Wyngaard, Y. Izumi, O. R. Coté,  Spectral characteristics of surface-layer turbulence,  Quarterly Journal of the Royal Meteorological Society 98 (1972) 563–589.
  * Porté-Agel et al. (2000) F. Porté-Agel, C. Meneveau, M. B. Parlange,  A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer,  Journal of Fluid Mechanics 415 (2000) 261–284.
  * Bou-Zeid et al. (2005) E. Bou-Zeid, C. Meneveau, M. Parlange,  A scale-dependent lagrangian dynamic model for large eddy simulation of complex turbulent flows,  Physics of Fluids 17 (2005) 025105.
  * Deskos (2019) G. Deskos, Numerical simulations of wind turbine wakes, Ph.D. thesis, Imperial College London, 2019.
  * Wyngaard (2010) J. C. Wyngaard, Turbulence in the Atmosphere, Cambridge University Press, 2010\.
  * Veers (1988) P. S. Veers, Three-Dimensional Wind Simulation, Technical Report, Sandia National Laboratories, 1988\. SAND88-0152 UC-261.
  * Mann (1994) J. Mann,  The spatial structure of neutral atmospheric surface-layer turbulence,  Journal of Fluid Mechanics 273 (1994) 141–168.
  * Mann (1998) J. Mann,  Wind field simulation,  Probabilistic Engineering Mechanics 13 (1998) 269 – 282.
  * Muñoz-Esparza et al. (2015) D. Muñoz-Esparza, B. Kosović, J. van Beeck, J. Mirocha,  A stochastic perturbation method to generate inflow turbulence in large-eddy simulation models: Application to neutrally stratified atmospheric boundary layers,  Physics of Fluids 27 (2015) 035102.
  * Jeong et al. (2014) M.-S. Jeong, S.-W. Kim, I. Lee, S.-J. Yoo,  Wake impacts on aerodynamic and aeroelastic behaviors of a horizontal axis wind turbine blade for sheared and turbulent flow conditions,  Journal of Fluids and Structures 50 (2014) 66–78.
  * García et al. (2017) J. García, J. Muñoz-Paniagua, A. Crespo,  Numerical study of the aerodynamics of a full scale train under turbulent wind conditions, including surface roughness effects,  Journal of Fluids and Structures 74 (2017) 1–18.
  * Ly and Chan (1980) U. Ly, Y. Chan,  Time-domain computation of aircraft gust covariance matrices,  in: 6th AIAA Atmospheric Flight Mechanics Conference, Danvers, Massachussetts, USA. AIAA Paper No 1980-1615.
  * Campbell (1986) C. W. Campbell,  Monte Carlo Turbulence Simulation Using Rational Approximations to von Kármán Spectra,  AIAA Journal 24 (1986) 62–66.
  * MIL-HDBK-1797 (1997) MIL-HDBK-1797, U.S. Military Handbook MIL-HDBK-1797, Technical Report, Department of Defense, USA, 1997\.
  * IEC-61400 (2005) IEC-61400, IEC-61400-1 Wind turbines-Part 1: Design requirements, 2005. The European Standard EN 61400-1:2005 has the status of a British Standard.
  * Jonkman and Kilcher (2012) B. J. Jonkman, L. Kilcher, TurbSim User’s Guide: Version 1.06.00, Technical Report, National Renewable Energy Laboratory, 2012\.
  * Deskos et al. (2018) G. Deskos, S. Laizet, M. D. Piggott,  Development and validation of the higher-order finite-difference wind farm simulator, winc3d,  in: 3rd International Conference on Renewable Energies Offshore (RENEW2018), Lisbon, Portugal, 8–10 October 2018.
  * Smagorinsky (1963) J. Smagorinsky,  General circulation experiments with the primitive equations,  Monthly Weather Review 91 (1963) 99–164.
  * Mason and Thomson (1992) P. J. Mason, D. J. Thomson,  Stochastic backscatter in large-eddy simulations of boundary layers,  Journal of Fluid Mechanics 242 (1992) 51–78.
  * Moeng (1984) C.-H. Moeng,  A large-eddy-simulation model for the study of planetary boundary-layer turbulence,  Journal of the Atmospheric Sciences 41 (1984) 2052–2062.
  * Laizet and Lamballais (2009) S. Laizet, E. Lamballais,  High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy,  Journal of Computational Physics 228 (2009) 5989 – 6015.
  * Gaitonde and Visbal (1998) D. V. Gaitonde, M. R. Visbal, High-order schemes for Navier-Stokes equations: Algorithm and Implementation into FDL3DI, Technical Report, Air Force Research Laboratory, 1998.
  * Laizet and Li (2011) S. Laizet, N. Li,  Incompact3d: A powerful tool to tackle turbulence problems with up to o(105) computational cores,  International Journal for Numerical Methods in Fluids 67 (2011) 1735–1757.
  * Xie et al. (2015) S. Xie, N. Ghaisas, C. L. Archer,  Sensitivity issues in finite-difference large-eddy simulations of the atmospheric boundary layer with dynamic subgrid-scale models,  Boundary-Layer Meteorology 157 (2015) 421–445.
  * Murua et al. (2012) J. Murua, R. Palacios, J. M. R. Graham,  Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics,  Progress in Aerospace Sciences 55 (2012) 46 – 72.
  * Hesse and Palacios (2016) H. Hesse, R. Palacios,  Dynamic Load Alleviation in Wake Vortex Encounters,  Journal of Guidance, Control, and Dynamics 39 (2016) 801–813.
  * Del Carre and Palacios (2019) A. Del Carre, R. Palacios,  Low-altitude dynamics of very flexible aircraft,  in: Proceedings of the AIAA Scitech 2019 Conference, San Diego, California, USA. AIAA Paper No 2019-2038.