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

Georgios Deskos, Alfonso del Carre, Rafael Palacios

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´ arm´n model, the two-dimensional Kaimal model and full three-dimensional wind fields a 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´ arm´n models). a 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

## 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

Email address: r.palacios@imperial.ac.uk (Rafael Palacios)

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 a a 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 a a 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 a a altitudes (e.g. CAT). The low-altitude operating environment of a generic VFA is shown schematically in figure 1. The near-surface turbulence coherent structures (e.g. high/lowspeed 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

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-48c47dea-cd04-48da-b2bf-e9e7c9943f14.jpg)

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 a a 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 onedimensional von K´rm´n, a two-dimensional stochastic formulation Kaimal model (Kaimal a a 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., e 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 ( ∼ hours) and the time spent by a VFA in the lower ABL ( ∼ 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 secondorder 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., n 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 a a 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, coh uu =1. In this study we use the von K´rm´n model. a a 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 _ { \text{ref} } } \cdot \frac { 1 } { \left [ 1 + \left ( 1. 3 3 9 L _ { u } \frac { 2 \pi f } { U _ { \text{ref} } } \right ) ^ { 2 } \right ] ^ { 5 / 6 } }.$$

for the longitudinal velocity and

$$S _ { v, w } ( f ) = \frac { 2 \sigma _ { v, w } ^ { 2 } L _ { v, w } } { \pi U _ { \text{ref} } } \cdot \frac { 1 + \frac { 8 } { 3 } \left ( 2. 6 7 8 L _ { v, w } \frac { 2 \pi f } { U _ { \text{ref} } } \right ) ^ { 2 } } { \left [ 1 + \left ( 2. 6 7 8 L _ { v, w } \frac { 2 \pi f } { U _ { \text{ref} } } \right ) ^ { 2 } \right ] ^ { 1 1 / 6 } }.$$

for the lateral velocities. As our interest is in the low-altitudes, independent turbulence length scales L , L , L u v w are assumed in three directions. They are chosen as in MIL-HDBK1797 (1997) so that 2 L w = h and L u = 2 L v = h/ (0 177+0 000823 . . h ) 1 2 . , where h is the altitude (in feet). The corresponding turbulence intensities are σ , σ u v , σ w . For the simulations, a timesignal of 250 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),

$$\coth _ { u u } ( f, r, U _ { \text{ref} } ) = \exp \left [ - \, a \sqrt { \left ( \frac { r f } { U _ { \text{ref} } } \right ) ^ { 2 } + \left ( \frac { r b } { L _ { u } } \right ) ^ { 2 } } \right ]$$

where a =10 and b =0 12 are the coherence decrement and offset parameters, . L u the streamwise turbulence lengthscale, U 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 inhouse 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 skewsymmetric form

$$\frac { \partial \widetilde { u } _ { i } } { \partial t } + \frac { 1 } { 2 } \left ( \widetilde { u } _ { j } \frac { \partial \widetilde { u } _ { i } } { \partial x _ { j } } + \frac { \partial \widetilde { u } _ { i } \widetilde { u } _ { j } } { \partial x _ { j } } \right ) = - \frac { 1 } { \rho } \frac { \partial \widetilde { p } ^ { * } } { \partial x _ { i } } - \frac { 1 } { \rho } \frac { \partial p _ { \infty } } { \partial x _ { i } } - \frac { \partial \tau _ { i j } } { \partial x _ { j } }, \quad i, j = \{ 1, 2, 3 \},$$

$$\frac { \partial \widetilde { u } _ { i } } { \partial x _ { i } } = 0,$$

where ˜ p ∗ = ˜ p + 1 3 / ˜ ˜ u u i i , and ˜ u i are the filtered components of the modified pressure and velocity fields, respectively, and ρ is the fluid density. The indices (1 2 3) correspond to a , , ( x, y, z ) coordinate system. The sub-filter stresses -∂ τ j ij (residual terms from the filtering process) are calculated using the standard Smagorinsky model (Smagorinsky, 1963)

$$\tau _ { i j } = - 2 ( C _ { S } \Delta ) ^ { 2 } | \widetilde { S } | \widetilde { S } _ { i j }, \quad \widetilde { S } _ { i j } = \frac { 1 } { 2 } \left ( \frac { \partial \widetilde { u } _ { i } } { \partial x _ { j } } + \frac { \partial \widetilde { u } _ { j } } { \partial x _ { i } } \right ),$$

where C S is the Smagorinsky constant which is corrected near the wall using the Mason and Thomson (1992)

$$C _ { S } = \left ( C _ { 0 } ^ { n } + \left \{ \kappa \left ( \frac { z } { \Delta } + \frac { z _ { 0 } } { \Delta } \right ) \right \} ^ { - n } \right ) ^ { - 1 / n }.$$

where z 0 is the roughness lengthscale, ∆ = 3 √ ∆ ∆ ∆ x y z the grid size and choosing ( C , n 0 ) = (0 14 . , 3). Finally, -∂ p i ∞ 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 ) },$$

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 ].$$

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 ( ∂ u 3 ˜ i = ˜ u 3 = 0 for i ∈ 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 , u 1 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 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, a a 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 0001m, 0 001m and 0 01m. . . . This defines three flow conditions that will be referred to as Mild , Regular and Severe , respectively, and correspond to surface friction velocities u ∗ =0 174ms . -1 , 0 235 m s . -1 and 0 274 m s . -1 . The friction velocities have been selected such that the mean velocity at altitude 10 m attains a value of approximately U 10 ≈ 5 m s -1 in all cases. Further details are provided in the table 1 below. Before pro-

| Table 1:   | Atmospheric turbulence flow parameters   | Atmospheric turbulence flow parameters   | Atmospheric turbulence flow parameters   | Atmospheric turbulence flow parameters   |
|------------|------------------------------------------|------------------------------------------|------------------------------------------|------------------------------------------|
| Name       | z 0 [m]                                  | u ∗ [m/s]                                | U 10 [m/s]                               | 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                                    |

ceeding 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 × L y × L z = 1000m × 500 m × 500 m discretised by N x × N y × N z = 1024 × 512 × 513 mesh nodes. The time step was selected to be equal to ∆ t =0 05s 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 ( δ/u ∗ , where δ =1000m 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 × N z =32 × 32

nodes was used for a domain of L y × L z = 50m × 50 m, obtaining a resolution of approximately 1 5 m. The calculations were made for a total analysis time of 250 s using a time step . ∆ =0 25s. t . 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 a a three velocity components σ u , σ v and σ 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 a a at the reference altitude. Larger discrepancies can be found between the the models in the

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-c426e8db-c472-467a-8507-f3cb9e06db28.jpg)

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 a a distributions for all components. 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 a a isotropy which results in nearly identical PSDs over the v ′ and w ′ components. 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 m while Kaimal consists of smaller structures extending only a few metres along the y-direction. The von K´ arm´n model on the other hand is not plotted as only a a uniform distribution will be assumed in a 2D planar snapshot. These qualitatively observed

Figure 3: Power spectral density of velocity fluctuations for each velocity component ( u ′ , v ′ and w ′ ) and all three models (1D von K´rm´n, 2D Kaimal and 3D LES). a a

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-dee7a1d5-3543-4e95-9afb-4d7a0802f064.jpg)

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-99e62743-ad97-4ecd-aa70-e32896550b22.jpg)

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,

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-b18d0fd8-e8ec-428f-b761-2bd9d02a3017.jpg)

with a discretisation based on nodal displacements and rotations, which are given in a bodyattached, moving frame of reference A (this frame moves and rotates with respect to an inertial frame, G ) as shown in figure 5.

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 η , the state variable containing the nodal displacements and rotations, and [ v A , ω A ], the state variable of velocities and rotations of frame A , as (Hesse and Palacios, 2016)

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

where M is the mass matrix and Q gyr , Q stiff and Q 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β 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 Γ 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 F _ { \text{steady} } & = \rho _ { \infty } | | \Gamma | | ( \mathbf v \times \partial l ) \\ & \quad \quad \quad$$

where k ∈ [0 , 4) is the segment counter and v , ∂ l the structure segment's velocity and length. The unsteady force at the centre point of the panel is calculated via,

$$\partial F _ { \text{ unsteady} } = \underset { \cdot \dots } { \rho } _ { \infty } | | \dot { \mathbf r } | | ( t )$$

where ˆ is a panel-weighted velocity and v ˙ || Γ || ( ) is t 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 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 × 512 × 513 ∼ 270 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-815e4c3f-e69f-4c1c-82e3-3e6aa8e49be5.jpg)

## 5. Numerical results

As discussed in Section 1, in the design process of aircraft the 1D von K´rm´n spectrum a a 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), a a 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 ref =10m.

## 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.

30%

Figure 7: Schematic representation of the T-tail aircraft geometry.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-9c739620-2a25-4951-b203-0f36619c2af2.jpg)

Table 2: Stiffness and mass parameters for the aircraft model.

| Component         | EA [N · m]        | GJ [N · m 2 ]     | EI y GJ [N · m    | EI z [N · m 2     | ¯ m [kg · m - 1   | ¯ J [kg · m]      |
|-------------------|-------------------|-------------------|-------------------|-------------------|-------------------|-------------------|
| Flexible Aircraft | Flexible Aircraft | Flexible Aircraft | Flexible Aircraft | Flexible Aircraft | Flexible Aircraft | 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    | Stiff Aircraft    | Stiff Aircraft    | Stiff Aircraft    | Stiff Aircraft    | Stiff Aircraft    | 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 s of real time and ∆ t =0 025s. The pro-. duced 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 =0m, 10m, 50m, 200m and -200 m are selected based on the lateral turbulence integral length scale, which is calculated to be between 150 m to 200 m for the three turbulence intensity cases.

Table 3: Summary of the various simulations along with some key parameters

| Case   |   U 10 (m/s) |   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 a a     |
| VK2    |         5.41 |      10.74 | -               | von K´rm´n a a     |
| VK3    |         4.74 |      12.73 | -               | von K´rm´n a a     |

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 m with 'dives' and 'jumps' of around 1 m to 3 m occurring throughout. A significant deviation from this pattern is observed in the response to the von K´rm´n a a 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 =0m, 10m, 50m, 200m and -200 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´ arm´n), a which are characterised by less or no variability in the transverse direction, result in nearly straight flight trajectories. To demonstrate the effect of spanwise

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-8ecc6ab3-7544-4f5f-a6ff-41fb1822e1ff.jpg)

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 =50m 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. 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

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-910a5232-c36d-414b-a953-6503e7282b2b.jpg)

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 N 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 Nm and 370 Nm, 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 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. Nonetheless, in order

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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-aa6fd7dc-4e79-40dd-b2c9-991540dc7ffa.jpg)

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. 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

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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-de33abb3-965e-4d3e-b1b3-2661040a062a.jpg)

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. a a 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 a a 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. Another important difference that is captured by the a a 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 (∆ u = u t ( +∆ ) t -u t ( ), ∆ M = M t ( +∆ ) 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 illustrates the existence of both intermittent velocity difference and loads, particularly for the 'severe' case in which a clear heavy tail

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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-2970ef09-a301-429e-9fbe-8ef510781335.jpg)

distributions, | ∆ U ∞ | &gt; . 0 1 m s -1 and | ∆ M OOP | &gt; . 7 5 N 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 | ∆ 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 a a is characterised by significantly less variability and thus, its PDF appears much narrower. 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 a a 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 ∼ 0 001 . in the high frequencies. In addition, all three cases exhibit a peak around f =1 3Hz which . corresponds to the mode II (natural torsion frequency). The peak around mode II is more

Figure 13: PDF of the incident velocity magnitude∆ | U ∞ | and out-of-plane moment M OOP for all models and turbulence intensity cases. Normal distribution is included for comparison.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-f52e06d1-ffdb-4574-a505-7ec32029a4ca.jpg)

pronounced for the severe turbulence case while modes I ( f =0 4Hz) and V ( . f =3 92Hz) . 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 Hz to 1 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 a a about an order of magnitude smaller than their respective counterparts generated by LES and Kaimal turbulence inputs.

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-6f719d6e-adf4-4969-9e07-d69d99dce046.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-2ed519ef-3541-4109-8f53-723ebdc6d421.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-579c6ae3-b39a-4339-8170-60c59a6e2f72.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-0c8fb252-9cbb-4195-b424-a1d6d1c39ff8.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-c914e7bc-a796-4400-bcf5-b120ea3af8dc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-226ef062-e455-42db-9d69-8f4918527f31.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-4e4a90dc-4f9e-4c73-b567-2ff0373c6e9d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-3fc03cf9-8803-401c-9630-b9e6cf782b0d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-c18f00b6-39bf-4749-80db-63aaa3846e4c.jpg)

## 5.3. Aeroelastic effects assessment

Having demonstrated the differences between three turbulence input models (von K´rm´n, a a 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. Aircraft flexibility is found to fundamentally change

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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-3fb2b5de-9cb0-4209-9138-3af5d1de56ba.jpg)

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 N m to -2500 N m while experiencing only positive torsion values throughout. This shift in the loads response is also reflected in the

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00372v1-figure-67eb1b85-2366-4514-8a88-caf171fb2adc.jpg)

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 threedimensional 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´ arm´n model was the a worst performing method of those considered, particularly when used in a moving threedimensional 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

- J. C. Hunsaker, E. B. Wilson, Report on Behavior of Aeroplanes in Gusts, Technical Report, NACA Technical Report 1, NACA-TR-1, 1917.
- R. T. Jones, The unsteady lift of a wing of finite aspect ratio, Technical Report, National Advisory Committee for Aeronautics, 1940. NACA-TR-681.

- J. W. Miles, The aerodynamic force on an airfoil in a moving gust, Journal of the Aeronautical Sciences 23 (1956) 1044-1050.
- J. K. Zbrozek, Atmospheric gusts: Present state of the art and further research, The Journal of the Royal Aeronautical Society 69 (1965) 27-45.
- F. D. Eichenbaum, A general theory of aircraft response to three-dimensional turbulence, Journal of Aircraft 8 (1971) 353-360.
- J. C. Houbolt, Atmospheric Turbulence, AIAA Journal 11 (1973) 421-437.

European Aviation Safety Agency, Certification Specifications for Certification Specifications for Normal-Category Aeroplanes CS-23, Technical Report, European Aviation Safety Agency, 2007.

- T. C. Tozer, D. Grace, High-altitude platforms for wireless communications, Electronics Communication Engineering Journal 13 (2001) 127-137.
- 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.
- J. Gonzalo, D. L´pez, D. Dom´ ınguez, A. Garc´ ıa, A. Escapa, o On the capabilities and limitations of high altitude pseudo-satellites, Progress in Aerospace Sciences 98 (2018) 37 56.
- 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.
- C. E. S. Cesnik, R. Palacios, E. Reichenbach, Re-examined Structural Design Procedures for Very Flexible Aircraft, Journal of Aircraft 51 (2014) 1580-1591.
- F. Afonso, J. Vale, E. Oliveira, F. Lau, A. Suleman, ´ A review on non-linear aeroelasticity of high aspect-ratio wings, Progress in Aerospace Sciences 89 (2017) 40-57.
- D. K. Schmidt, Modern Flight Dynamics, Mc-Graw Hill, New York, NY, USA, 2012.
- M. C. Van Schoor, A. H. Von Flotow, Aeroelastic characteristics of a highly flexible aircraft, Journal of Aircraft 27 (1990) 901-908.
- 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.
- 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.

| 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.                                                                                                       |
|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| 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.                                                                    |
| 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.                                                                                               |
| 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.                                                                             |
| F. Port´-Agel, e 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.                                                 |
| 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.                                                                                   |
| G. Deskos, Numerical simulations of wind turbine wakes, Ph.D. thesis, Imperial College London, 2019.                                                                                                                                                          |
| J. C. Wyngaard, Turbulence in the Atmosphere, Cambridge University Press, 2010.                                                                                                                                                                               |
| P. S. Veers, Three-Dimensional Wind Simulation, Technical Report, Sandia National Labo- ratories, 1988. SAND88-0152 UC-261.                                                                                                                                   |
| J. Mann, The spatial structure of neutral atmospheric surface-layer turbulence, Journal of Fluid Mechanics 273 (1994) 141-168.                                                                                                                                |
| J. Mann, Wind field simulation, Probabilistic Engineering Mechanics 13 (1998) 269 - 282.                                                                                                                                                                      |
| D. Mu˜oz-Esparza, n 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. |
| 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.                         |
| J. Garc´a, ı J. Mu˜oz-Paniagua, n 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.                                 |
| 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.                                                                    |

- C. W. Campbell, Monte Carlo Turbulence Simulation Using Rational Approximations to von K´ arm´n Spectra, AIAA Journal 24 (1986) 62-66. a
- MIL-HDBK-1797, U.S. Military Handbook MIL-HDBK-1797, Technical Report, Department of Defense, USA, 1997.
- 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.
- B. J. Jonkman, L. Kilcher, TurbSim User's Guide: Version 1.06.00, Technical Report, National Renewable Energy Laboratory, 2012.
- G. Deskos, S. Laizet, M. D. Piggott, Development and validation of the higher-order finitedifference wind farm simulator, winc3d, in: 3rd International Conference on Renewable Energies Offshore (RENEW2018), Lisbon, Portugal, 8-10 October 2018.
- J. Smagorinsky, General circulation experiments with the primitive equations, Monthly Weather Review 91 (1963) 99-164.
- P. J. Mason, D. J. Thomson, Stochastic backscatter in large-eddy simulations of boundary layers, Journal of Fluid Mechanics 242 (1992) 51-78.
- 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.
- 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.
- 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.
- S. Laizet, N. Li, Incompact3d: A powerful tool to tackle turbulence problems with up to o(10 ) computational cores, 5 International Journal for Numerical Methods in Fluids 67 (2011) 1735-1757.
- 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.
- 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.
- H. Hesse, R. Palacios, Dynamic Load Alleviation in Wake Vortex Encounters, Journal of Guidance, Control, and Dynamics 39 (2016) 801-813.

- 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 20192038.
Georgios Deskos Alfonso del Carre Rafael Palacios
Assessment of low-altitude atmospheric turbulence models for aircraft aeroelasticity
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{\'a}rm{\'a}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{\'a}rm{\'a}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.