Wake Prediction and Modeling Improvement of an existing model

(1)

Wake Prediction and Modeling

Improvement of an existing model

Author:

Jacques Wild

Supervisor at AVTECH:

David Rytter

Supervisor at KTH:

Arthur Rizzi

Supervisor at ISAE:

Sébastien Massenot

Master Thesis Dissertation

Double-degree with ISAE Supaero, Toulouse and KTH, Royal Institute of Technology, Stockholm

Toulouse, France and Stockholm, Sweden 2010–2014

(2)

List of Figures

1 Typical two phase decay of a wake vortex . . . 10

2 LIDAR velocity . . . 15

3 Wind Spectrum . . . 18

4 VAD scan used by Frehlich in [9] . . . 18

5 Accepted EDR calculation through velocity function . . . 23

6 Accepted EDR calculation through velocity function . . . 23

7 Rejected velocity structure function . . . 24

8 Rejected velocity structure function . . . 24

9 EDR September 3rd . . . 25

10 Mean EDR in september 2013 400 meters away from the LIDAR . . . 26

11 EDR profile in calm atmosphere (Heat flux is negative) . . . 28

12 EDR profile in turbulent atmosphere (Heat flux is positive) . . . 29

13 EDR profile in turbulent atmosphere (Heat flux is positive) . . . 29

14 Wake created by an A321 bouncing on an inversion layer on November 20th, 2013 10:07:54 . . . 32

15 Vertical Velocity computation . . . 32

16 Distribution of the descent velocities . . . 33

17 EDR profile at 2 UTC, stable boundary layer . . . 38

18 EDR profile at 8 UTC, unstable boundary layer . . . 38

19 EDR profile at 14 UTC, unstable boundary layer . . . 38

20 EDR profile at 20 UTC, stable boundary layer . . . 38

21 Longitudinal Position of the wakes . . . 39

22 Altitude of the wakes . . . 39

23 Altitude . . . 40

24 Track in YZ plane . . . 40

25 Strength of the wakes . . . 41

26 High Crosswind, Longitudinal Position of the wakes . . . 42

27 High Crosswind, Altitude of the wakes . . . 42

32 Strength of the wakes . . . 45

(5)

Abstract

This thesis is dedicated to the study of atmospheric turbulence and its effects on wake-vortex decay. Wake-vortices are a by-product of the lift that makes every aircraft fly. Under certain special conditions, the decay-time of a wake-vortex can be high and thereby impose constrains on separation between aircrafts. A precise knowledge of background turbulence will hence improve the ability of predicting the decay of a wake-vortex.

To compute the turbulence, LIDAR measurements have been used to calculate an EDR value (Eddy Dissipation Rate) of the air mass of interest. This parameter indeed quantifies the background turbulence and is of importance in the decay model. The results of this calcula-tion will then be used to improve an overall wake-vortex simulacalcula-tion tool developed by AVTECH.

(6)

Extrait

Ce rapport vient concrétiser mon projet de fin d’études qui avait pour sujet l’étude des turbulences atmosphériques et leur effet sur la dissipation des tourbillons. Ces tourbillons sont un effet secondaire de la portance d’un avion et sont créés lors de la phase de vol. Sous certaines conditions, le temps de dissipation de ces tourbillons peut s’avérer être élevé, ce qui influera la distance de séparation entre avions. Une connais-sance précise de la turbulence atmosphérique permettra donc d’améliorer la prédiction de la dissipation des tourbillons.

Pour calculer la turbulence, des mesures provenant d’un capteur LI-DAR ont été utiliées pour calculer une valeur d’EDR (Eddy Dissipation Rate ou Taux de Dissipation des Tourbillons) de la masse d’air étudiée. Ce paramètre quantifie en effet la turubulence atmosphérique joue un rôle important dans le modèle de dissipation des tourbillons. Les résul-tats de ces calculs ont par la suite été utilisés pour améliorer un outil de simulation plus complet des tourbillons, développé par AVTECH.

(7)

Referat

Denna rapport är tillägnad studerandet av atmosfärisk turbulens och dess effekter på vortexvirvlars (eng. wake-vortex) avmattning. Vor-texvirvlar är en biprodukt från lyftkraftsproduktionen som gör att fly-gplan flyger. Under vissa speciella förhållanden kan tiden det tar för vortexvirvlarna att avta vara lång och därmed påtvinga större avstånd mellan flygplan. En god kunskap om bakgrundsturbulens kan därför för-bättra förmågan att förutsäga tiden det tar för en vortexvirvel att avta. För att beräkna turbulensen har LIDAR-mätningar använts för att räkna ut den relevanta luftmassans EDR-värde (Eddy Dissipation Rate). Denna parameter kan kvantifiera bakgrundsturbulensen och är viktig i modellen som beräknar avmattningen. Resultatet av denna beräkning kommer sedan användas för att förbättra den generella simuleringen av vortexvirvlar med det verktyg som är utvecklat av AVTECH.

(8)

1 Introduction

1.1 Background

With a fast-pace growing air traffic, some of today’s hubs almost reached their maxi-mum capacity acknowledging the current separation standards. For instance, the current issue about the building of a new airport in London is highly related to this lack of landing and take-off slots. There are mainly two ways of tackling this issue. Firstly, it is building of new airports, which is often problematic in dense urban envi-ronments. Secondly, it is reduction of separation between aircrafts. This separation is the compulsory time delay between two aircrafts landing on the same runway to avoid the following plane being caught up in the vortices created by the leader aircraft.

The minimum separation aircrafts have to respect is a consequence of the wake vortices produced by every airplane that flies. A wake vortex is a wind eddy that can be hazardous if a plane flies into it, creating an unexpected roll moment that can lead to a loss of control of the aircraft in the most extreme cases.

The current standards have so far proven themselves to be safe. Nevertheless, they were designed to hold with every meteorological condition. This implies that in case of favourable weather conditions, the separation between aircrafts could potentially be reduced without loss of safety.

The work undertaken in this master thesis is addressing this problem by modelling air turbulence and its effects on wake vortex decay. The long term aim is to help AVTECH develop a reliable wake-vortex decay tool that can help airports to deal with wake vortex issues.

1.2 The project

AVTECH Sweden AB is a company that has been involved in modelling wake vortex decay since it became an issue for currently saturated airports. The company got involved with wake-vortex issues in 2010 and has completed projects involving Monte-Carlo simulations as a tool to analyse wake-vortex decay.

Several master theses were written on this project : David Rytter in [14] started the project by implementing the P2P-model described by Holzäpfel in [11]. This was the basis for all future work done on wake vortex at AVTECH. Erik Peldan continued this work in [12] recently by improving the solving process of the differential equation system.

AVTECH recently got access to wind data produced by two Lockheed-Martin LIDAR device that can observe wake decay directly at one airport in the middle East. The first of the two LIDARs can observe wakes out of ground effect (see Peldan in [12] for more explanation about wakes in ground effect and out of ground effect) whereas the second one focuses on in ground effect vortices. These LIDARs measure the radial wind velocity, from whom the EDR can be calculated to directly test the model that has been developed. Furthermore, a LIDAR can be equipped with wake-tracking software and observe wakes directly. The analysis of the LIDAR data and its exploitation is what this thesis will be focusing upon.

1.3 Acknowledgment

(9)

to thank my academic tutors, Arthur Rizzi from KTH and Sébastien Massenot from ISAE who accepted to supervise my work.

I would like to thank AVTECH Sweden for giving me the opportunity to partici-pate in a project where I learned a great deal about aviation safety and how a project is managed and run.

2 The Model

This section describes the physical equations that are used in the Monte-Carlo simu-lation. For more information on this topic, Peldan gave many details in [12].

2.1 The decay model

AVTECH used a wake-vortex decay model developed by Holzäpfel at the German DLR as a milestone for the whole wake-vortex project. This article (see [11] ) is the basis for predicting the decay of wake vortices. This model establishes two different phases in a wake vortex decay. The first phase, the diffusion phase has a slow decay whereas during the second phase, the rapid decay phase, the wakes decay much faster. The decay for these two phases can be found in the following equations:

Γ(r, t) = A − exp − r 2 4ν1(t − T1) (1)

for the first phase, whereas as the faster decay from the second phase is described by the following equation:

Γ(r, t) = A − exp −r2 ν1(t − T1) − exp −r2 ν2(t − T2) (2)

An example of the two-phase decay can be seen on figure 1. LIDAR observations of such decay can be found in the annex.

Γ∗ is an average circulation 5 to 15 meters away from the core radius. This is the quantity that is usually used in wake vortex studies such as [11] or [13]. These equations are the basis for the wake behaviour in the Monte-Carlo simulations carried out by AVTECH and will be used in this project every time a wake is simulated.

The parameters T1 is the time where the diffusion phase starts (shortly after

creation of the wake vortex) and the time T2is the time where the decay phase starts.

These parameters can be adjusted depending on the EDR value, and this was one of the motivations for this thesis. An excellent track where T1 and T2 are visible

can be observed in section 9.3 on figure 25. The parameters ν1 and ν2 are viscosity

(10)

Figure 1: Typical two phase decay of a wake vortex

3 LIDAR

The data that is used in this master thesis comes from a LIDAR observing wakes created oncoming aircrafts out of ground effect and another LIDAR observing wakes in ground effect.

Both LIDARs are equipped with a wake-vortex tracking software and their primary purpose is to be able to observe wake-vortices. Nevertheless, it will foremost be used to estimate the EDR in this master thesis, which is a measurement of the background turbulence. It has to be controlled that no plane-induced turbulence is hanging in the air when a background turbulence is to be made.

3.1 Physical quantities used in LIDAR physics

A LIDAR (Light Detection and Ranging) is a measurement device which has a wide range of application and is famous for having measured the distance Earth-Moon with a good precision. A LIDAR device analyses the backscattered signal it receives. It can for instance be used to measure distances or the properties of the air the LIDAR beam goes through. This last case is the case that will actually be used in this thesis. A LIDAR measures the radial velocity of the wind velocity field. Nevertheless, the velocity cannot be measured at every point in space since a LIDAR has a definite spatial resolution. It can therefore deliver only one measurement of the wind every ∆x meters. The letter R will denote the radial coordinate for the measurement, and t the time of the measurement in the following. The radial velocity measured by the LIDAR at distance R and time t is denoted v(R, t). What velocity is measured is explained in more details in section 3.4.

(11)

• PRF is the Pulse Repetition Frequency (in Hz). It is the number of pulses per second of the laser. For the LIDAR data I have access to: P RF = 750 Hz. • The sampling rate SR

SR = 2.5 ∗ 108Hz (3)

• The range per sample defined by the following equation rangeP erSample = c

2× SR; (4)

• τ is the FWHM of the pulse envelope, or the laser pulse direction. For this LIDAR, τ = 3 × 10−7 s.

• ∆r is the pulse range, i.e. the spatial extension of the LIDAR pulse, the region illuminated by the LIDAR pulse. By definition, this quantity is equal to (c is the speed of light):

∆r = c∆T

2 = 45 m (5)

• ∆p is referred to as the length to a range gate. It is the distance the pulse is moving during the measurement (estimation of the velocity).

∆p =M Tsc

2 = 76.8 m (6)

M is the amount of measured samples per gate. This quantity for the undertaken measurements is 128.

• ∆R is the range resolution of the LIDAR. It is defined as the sum of the distance the pulse is moving between two measurements and its spatial width: ∆R = ∆r + ∆p = 121.8 m.

• TS is the sampling interval of the data. In this case, we have access to the

sampling rate (parameter set for the LIDAR measurements) which is 2.5 × 108

Hz and yields a sampling time of about 4 × 10−9 seconds.

The gates the LIDAR uses are highly overlapping since the distance between two consecutive gates is 9.6 meters according to the LIDAR manual and the preceding LIDAR parameters. In the remaining part of this thesis, capital R will denote the radial distance between the LIDAR and the centre of one range gate. The lower case letter r will denote the radial distance between the LIDAR and the current point.

3.2 What a LIDAR is measuring

(12)

In the rest of the thesis, the first EDR calculation will be held with the metho-dology described in [8]. This paper indeed uses a vertically staring LIDAR to get its data. This situation is the closest one to the scan pattern that is used by the LIDAR. It can be considered that a LIDAR scanning the same plane and passing every 8 seconds through the same angle is the same as a LIDAR pointing all the time to the same direction but with a 8 seconds delay between two measurements.

3.3 LIDAR measurements

This section focuses on the different calculations that need to be done before any turbulence estimation. The following quantities are requested to calculate the EDR value.

The fluctuating or turbulent component of the radial velocity is defined as:

v0(r, t) = v(r, t) − hv(r)i (7)

The brackets h.i indicate a time average in the following.

The EDR calculation is then based on the following quantities, calculated with the estimation of the turbulence component of the radial velocity:

• The covariance: Bv(s, r) = D v0r −s 2 v0r +s 2 E (8) • The structure function:

Dv(s, r) = D v0r −s 2 − v0r + s 2 E (9)

Different models for the structure function have been used in previous studies regarding air turbulence. Frehlich describes in [8] that in case of isotropic turbulence and high Reynolds number, the structure function can be modelled by:

Dv(s, r) = Cv2/3(r)s2/3 (10)

An empirical model for the structure function used by Frehlich in [9] and widely used by the scientific community working with EDR is:

Dv(s, r) = 2σ2v(r)Λ[s/L0(r)] (11)

The parameters s, which is the separation (in meters) of the two observation points and the parameter r which is the centre of the observation points play an important role in the EDR calculation.

(13)

3.4 Velocity estimated by a LIDAR

According to [6], several different velocities can be computed from a LIDAR signal: the mean linear velocity, the pulse weighted velocity and the conditional ensemble average velocity. The desired measured velocity is named m and can be chosen to be either one of those.

Furthermore, the key point is, as explained in [7]: «The accuracy of a velocity estimates ˆv requires a definition of the desired measurement which is typically chosen as the conditional ensemble average mave(R, t) of the velocity estimate for a given

realisation of the atmosphere v(R, t). »

The estimate of the velocity given by the LIDAR is by definition an averaged quantity, because the only quantity the LIDAR has the capacity to measure is the average of the velocity around the range gate denoted by R. According to [6], three different velocities can be chosen for the random variable that has to be estimated.

The three quantities are:

mlin= 1 ∆p Z R+∆P2 R−∆P 2 v(r, t)dr mave= hˆv[R, t | v(r)]i mwgt= Z ∞ −∞ v(r, t)W (R − r)dr

According to [6] and [7], the ensemble average mave velocity is approximated by

the pulse-weighted velocity mwgtwell. However, these references do not explain under

which hypothesis this approximation is valid, so for our purpose, it is assumed that this can be used.

The velocity mwgtis the one that is used estimated by the LIDAR used for this

thesis.

Once the desired measurement has been chosen an effective and reliable velocity estimator has to be used to get an estimation of the velocity. The estimator then estimates the velocity from the LIDAR signal. The quality of the estimation depends on the chosen estimator. The performances of different velocity estimators have been compared in [5].

The function W is related to the laser emission profile of the laser. It is defined by: W (s) = 1 ∆p Z ∆p/2 −∆p/2 In(r + s)ds (12)

This last function (lidar-pulse-range-weighting function for a range gate centred at distance R) is defined by In(r) = 1 √ πrP exp −(R − r) 2 r2 P (13) The velocity estimate delivered by the LIDAR is denoted ˆv. It is defined by the following equation, according to [6].

ˆ

(14)

Nevertheless, ˆv is an unbiased estimator and it can therefore be represented as : ˆ

v(R, t) = mwgt(R, t) + e(R, t) (15)

Hence the importance of the random error e and its statistical description. This choice of the velocity estimator is of importance, since it will determine the statistics of the error and the variance which the EDR will be computed from. A good understanding of what an estimator is required here. The estimation in itself depends as much on the choice of the estimator as on the sample that is used.

3.5 Maximum-Likelihood Estimator

The estimator that is used to estimate the velocity is the maximum likelihood esti-mator. It will be briefly described here. An estimator is used to infer the value of a parameter. It maps the sample space with a possible value that the parameter could take. It is as well a random variable and its statistics can be looked at as well. The statistics of estimators is therefore quite important to describe a good estimator, or at least an estimator with the desired properties.

The estimator is that is implemented in the LIDAR is the maximum likelihood estimator. For a set of n observations noted x1, x2, . . . , xn, the density function is:

f (θ|x1, x2, . . . , xn) = f (x1|θ) × f (x2|θ) × . . . × f (xn|θ) (16)

θ can be a scalar or the vector here. The idea behind the ML estimator is that it reverses this dependence and look at this function with fixed x1, x2, . . . xn. It then

computes the θ which minimizes this likelihood function.

L(θ|x1, . . . , xn) = n

Y

j=1

f (xj|θ) (17)

The post processing software uses this estimator to calculate the velocity based on the signal that is read by the LIDAR. Every velocity that is produced by the LIDAR is derived from the Maximum-Likelihood estimator. Frehlich compares the performances of this estimator to others in [5].

An example of a wind plot from the LIDAR can be seen on figure 2.

(15)

Figure 2: LIDAR velocity

4 Coding of software in C++ to extract relevant

data

The output file of a LIDAR data contains a large amount of measurements. To compute the EDR, only the radial velocity is relevant. Since the only output file available from the LIDAR was a raw data file, encrypted in a binary file, software aimed at extracting the relevant data had to be coded. Since the documentation about the organisation of the data was sometimes vague, this turned out to be a time-consuming process.

4.1 Binary/Hexadecimal convertor

The manufacturer of the LIDAR delivered a manual explaining how the data is orga-nised. This manual explains how to read the data in hexadecimal format. This is the reason why a binary to Hexadecimal function has to be created.

First, software that could read the binary file and create an hexadecimal output out of it. A UNIX function called "hexdump" already exists but the function had to be integrated in the produced software to calculate the EDR so that it needed to be coded as well.

Every line is 64 bytes long. The C++ function "read" was used instead of the function "getline" to store every line in memory since the delimiters "

(16)

The ".prd" file that is the extension for the binary file has a size of about 400 MB for ten minutes of LIDAR data recording. After extraction to hexadecimal format, it is about 1.3 GB megabytes.

4.2 Format of the output file

The organisation of the data is explained in a document. The raw LIDAR file contains velocity data, filtered velocity data, SNR data, filtered SNR data, spectral data, time and position for every data, configuration data, status data. Every different data type is identified by a header in hexadecimal format. This header is then followed by the size of the measurement. This can be summarised in the following table.

Content of the block Header Version number Size of the block Size of the block 2 bytes 2 bytes 4 bytes

4.2.1 Velocity measurement header

If the 2 bytes corresponding to a precise header are found in the data, it has to be checked that it is indeed a header and that the following bytes indeed contain a version number and a size. Several functions were coded to achieve the preceding and to assure that if two headers identification came up on the same line, both are read and only the actual header is selected.

Thereafter, the size of the block is stored, and the software reads the correct amount of bytes after that to get all the velocity measurements.

4.2.2 Record Header

For every measurement header comes a mesaurement info header that indicates the time and the position of the velocity measurement that will follow. The relevant information can be found in the following table:

Content of the block Header Version number Length of the Header Length of the whole record

Size of the block 2 bytes 2 bytes 4 bytes 4 bytes

Content of the block Year Month Day Hour Minute Second Nanoseconds Size of the block 2 bytes 1 byte 1 byte 2 bytes 1 byte 1 byte 4 bytes

4.3 Extracting relevant information

(17)

5 EDR Calculation

This section describes in the subsections 5.1.1 and 5.1.2 two different ways of calcula-ting a raw estimate from LIDAR data for the velocity structure function. Once this raw estimate is calculated, an estimate for the velocity structure function is given by the equation (32). The method that is used to calculate the EDR value is the one described in subsection 5.1.1. This choice will be justified at this point.

The fit between this velocity structure function and the theoretical velocity struc-ture function defined in section 5.2 gives the EDR.

5.1 The raw longitudinal structure velocity function

5.1.1 Calculation with a fixed beam Paper from 1998 [8]

This paper explains how to compute EDR from a LIDAR pointing in one direction only. The LIDAR data we have access to is a scan, but the LIDAR is pointing at the same direction every eight seconds, so that the difference with this article lays in the time interval between two measurements (one second in [8]). Nevertheless, the methodology described here is adapted to the set of data from the LIDAR.

For all the data presented here, ∆t is equal to 8 seconds . The velocity structure function given by [8] is:

ˆ D(R1− R2, R0) = 1 N N −1 X l=1 [ˆv0(R1, l∆t) − ˆv 0 (R2, l∆t)]2 ! (18)

It must be kept in mind, even if it is not explicitly written, that this structure function is used to calculate the EDR in R0. This point is in the middle of the two

points described by R1 and R2

The paper states the following: "This estimate is a true spatial estimate and does not require Taylor’s frozen hypothesis to map a temporal statistic to a spatial statistic. In practice, these estimates are averaged over a few range-gates."

An example of a wind spectrum can be seen on figure 3.

The wind spectrum often correlates with the Kolmogorov scaling. It has to be re-membered that the hypothesis under which the turbulent velocity temporal spectrum satisfies the Kolmogorov scaling is that the atmospheric turbulence is locally isotropic (see [8]).

5.1.2 Calculation with a varying beam Paper from 2005 [9]

(18)

Figure 3: Wind Spectrum

Figure 4: VAD scan used by Frehlich in [9]

(19)

ˆ Draw(k∆s) = 1 NT(NR− k) × NT X l=1 NR−k X j=1 ˆ v0[R1+ (j − 1)∆s, φ1+ (l − 1)∆φ, θ] − ˆv 0 [R1+ (j + k − 1)∆s, φ1+ (l − 1)∆φ, θ] (19) It takes into account both a time average and a spatial average. Problem might be that NR is the "number of range gates for each azimuth angle" and NT the "number

of velocity measurements for the given VAD scan". NT therefore is the amount of

different angles. This approach might be used as well. To be able to use it with the data set, one must assume that the result stands if the elevation angle and the azimuth are swapped. The data set indeed has a fixed azimuth but a varying elevation angle.

The approach that is appropriate is therefore the approach described in the ar-ticle [8].

These two different way of calculating the raw function are suited for different purposes. The approach described in [8] is used for one beam calculation whereas the approach taken into [9] is more adapted to VAD scans or the glide-path like calculation described in [1].

5.2 The theoretical structure function

The function Λ is defined in the von Karman model as the following function:

Λ(x) = 1 − 2 2 3 Γ(1/3)x 1/3 K1/3(x) (20)

According to Frehlich in [9], "If the fluctuations of the radial velocity are homoge-neous and isotropic over the plane defined by the fixed elevation angle θ", then "the second-order statistics can be described by the longitudinal structure function".

Still according to [9], "for a VAD lidar scan, the extent of the effective sensing volume transverse to the lidar beam is given by:"

∆h = R sin ∆φ (21)

If ∆h is small compared to ∆p, Dwgt is related to Dv(s). Still according to [9],

"for a Gaussian transmitted LIDAR pulse with a von Karman model for the turbulent wind field". This information has to be kept in mind if this approach is to be followed in the future. The azimuth angle ∆φ will probably be replaced by ∆θ.

The theoretical structure function is then defined by:

Dwgt(s, σ, L0) = 2σ2G(s/∆p, µ, ξ) (22)

The function G is defined by: G(m, µ, ξ) =

Z ∞

−∞

(20)

The function F is defined by:

F (x, µ) = 1

2√πµexp[−µ

2_{(x + 1)}2_{] + exp[−µ}2_{(x − 1)}2_{] − 2 exp(−µ}2_x2₎₊

1

2(x + 1)Erf[µ(x + 1)] + (x − 1)Erf[µ(x − 1)] − 2Erf(µx) (24)

Erf is the error function defined by: Erf(z) =√2

π Z z

0

exp(−t2dt) (25)

The parameters µ and χ are defined by

µ = (2 ln 2)1/2∆p

∆r (26)

χ = ∆p L0

(27)

The functions detailed here are used to calculate the theoretical structure function in (22).

5.3 The azimuth structure function

It is defined the following way in [9]. It is stated that "If the transverse dimension ∆ h of the LIDAR sensing volume for each radial velocity estimate is much less than the effective range resolution ∆R, the mean velocity vwgt can be approximated as

men-tioned beforehand. The Doppler LIDAR azimuth structure function then becomes: Dwgt(s, σ, L0) = 2σ2Gφ _s ∆p, µ, χ (28) In this case s = R(φ1− φ2) Gφ(m, µ, χ) = 2 Z ∞ 0 F (x, µ) Λχpm2_{+ x}2_{− Λ(χx) + Λ} D χpx2_{+ m}2 m 2 m2_{+ x}2 dx (29) The function Λ is defined in (20). Since different studies have shown that there is a strong correlation between the results obtained with these two different functions, only the azimuthal structure function will be considered in the following.

The function ΛD is defined by:

ΛD=

x4/3

21/3_Γ(1 3)

K2/3(x) (30)

(21)

ˆ Draw(R, kR∆φ, θ) = 1 NT − k Nt−k X j=1 ˆ v0[R, φ1+ (j − 1)∆φ, θ] − ˆv 0 [R, φ1+ (j + k − 1)∆φ, θ] (31)

5.4 Calculation of the structure function

The structure function estimate is simply derived from the raw estimate of the velocity structure function by the following equation:

ˆ

Dwgt(R1− R2) = ˆDraw(R1− R2) − E(R1− R2) (32)

The subsection 5.4.2 explains how to calculate the error function denoted by E(R1− R2), and ˆDraw(R1− R2) has been derived in the section 5.

5.4.1 Hypothesis

The turbulence of the air is well described by the azimuthal or longitudinal structure function when the van Karman model is valid. This model states that the turbulence field is isotropic and homogeneous as mentioned in [9].

5.4.2 Calculation of the error

In order to calculate an estimate of the structure function, according to [9], an un-biased estimate is given by the following formula:

ˆ

Dwgt(R1− R2) = ˆDraw(R1− R2) − E(R1− R2) (33)

E(R1− R2) is an unbiased correction for the contribution from the estimation

error e(R, φ, θ) of the velocity estimates. In [3], it is stated that a spectral approach is preferred since they have no access to the raw data. This error is computed trough the following formula:

ˆ σ2_∆e(R1, R2) = 1 (N/2 − jT + 1)∆t N/2 X j=jT ˆ Φ∆v(j∆f, R1, R2) (34)

The letter ˆΦ∆v represents the spectrum of the difference in velocities between

two gates. The error e(R, φ, θ) is equal to the quantity ˆσ2_∆e(R1, R2) which can be

calculated thanks to the spectrum of the difference in velocities. The spectral method is described more precisely in [8]

(22)

5.5 Matching of the structure functions and EDR calculation

5.5.1 Minimisation of cost function

Once the structure function has been calculated, it is compared to the theoretical one. This comparison is done by minimizing the following function.

χ2= 1 NT NT X k=1 h ˆ_D_wgt_{(k∆s) − D}_wgt_{(k∆s, σ, L}₀₎i2 D2 wgt(k∆s, σ, L0) (35)

Once L0 and σ have been calculated through the minimization of this function,

the EDR is finally computed with these values through the following formula:

= ₂4/3_π √ 3Γ(1/3)Γ(4/3)Cv 3/2 _σ3 L0 = 0.933668σ 3 L0 (36)

(23)

5.5.2 Results

The following two figures show in red the calculated Structure function with mea-sured LIDAR winds. In blue, the best match for the model prediction Dwgt and

the structure function is displayed. The software minimizes the function defined in equation (35) and gives the corresponding value of L0 and σ.

Figure 5: Accepted EDR calculation through velocity function

(24)

The following two figures show calculated structure functions (in red) that were not taken into account since the turbulence scale L0 is smaller than 50 meters which

is not enough for a good precision according to [3].

Figure 7: Rejected velocity structure function

Figure 8: Rejected velocity structure function

All the EDR measurements that have a small turbulence scale L0 usually don’t

(25)

5.5.3 EDR curves

The EDR is the measure of the atmospheric background turbulence. The following plots show the EDR calculated by LIDAR velocity measurements. The LIDAR that is used to track down wake-vortices is also use to calculate the background turbulence. A particular attention is therefore given not to measure wake turbulence that has been created by an aircraft and that would add turbulence to the background turbulence. The atmospheric turbulence described by the EDR varies during the day. Usually, atmospheric turbulence is very low at night, and increases during the day when the sun is warming up the lower layer of air.

The following plot shows the LIDAR measured variation over one day: September, 3rd 2013.

Figure 9: EDR September 3rd

The following curve shows the mean variation of the EDR over a day over the period September 2nd-September 15th 2013. One measurement per hour has been derived from LIDAR data. One measurement per hour was the best compromise between the quality of the measurement and the need for precision.

(26)

Figure 10: Mean EDR in september 2013 400 meters away from the LIDAR

6 EDR profile

The LIDAR measurement enables to compute the EDR at certain points in space. An EDR profile can be created following the methodology developed by Han in [10]. This article exposes a very simple way to establish a full EDR profile based on two EDR measurements. In [13], this model is used because of the need to create a low computational time demanding EDR profile.

In this section, an EDR profile will be build based on two LIDAR measurements. This section uses Han’s article [10] with the results from the LIDAR. It has to be mentioned that the results derived in [10] are based on similarity theory: the Monin-Obukhov similarity and the mixed-layer similarity theories.

6.1 Description of the algorithm

In our case, this model will be used to improve the quality of the simulations. The EDR value indeed has, as mentioned before, a crucial role in the wake-vortex decay and this algorithm gives an idea of what the EDR can look like. This computation is based on similarity theory.

This model uses quantities that are easy to compute and that give a plausible EDR profile if the atmosphere stays "standard". It won’t hold if an inversion layer hovers above the altitude where the EDR measurements are taken. But this model is simple to use and it is the first time that it will be implemented using LIDAR based EDR.

(27)

if the atmospheric boundary layer is stable or unstable. A stable boundary layer corresponds to the state of the atmosphere at night : calm air and few turbulence. The unstable boundary layer describes the turbulent air rising under the action of the heat flux coming from the ground. The Obukhov length is defined by :

L = − u 3 ∗ k ∗ (g/T0)( ¯w 0 θ_v0)s (37)

k is the von Karman constant, g is the gravitational acceleration and T0 is the

reference temperature. The Obukhov length, if a comparison with fluid mechanics is dared, describes the state of the atmospheric boundary layer. The Earth indeed has several layers composing the atmosphere and the atmospheric boundary layer is the one that "sticks" to it. The Obukhov length describes the state and the length of this layer.

u∗ is the friction velocity and ( ¯w0θv0)s is the surface heat flux. According to [10],

the friction velocity can be adjusted to get a continuous EDR profile at z = 40. The value of u∗ might therefore differ a little from the value indicated in (38). These

parameters can be calculated the following way:

u∗= k∆¯u φm(ζm) ln(z2/z1) (38) ( ¯w0 θ0 v)s= _k2_∆u∆θ v φm(ζm)φh(ζm) ln(z2/z1)2 (39)

These parameters are the parameters that have effect on the atmospheric tur-bulence. The surface heat flux changes sign during the day. At night, the ground radiates heat towards the cold air. During daytime, the ground absorbs the heat that is produced by the sun.

6.2 EDR formulas

The different EDR formula for the different types of atmospheric boundary layers are described here. The type of boundary layer is determined by the value of the preceding parameters.

6.2.1 Neutral or stable boundary layer

If the boundary layer is neutral or stable (_Lz ≥ 0) and if the current altitude is in the surface layer (z ≤ 0.1h), then the EDR has the following expression:

(z) = u 3 ∗ kz 1.24 + 4.3z L (40) h is the boundary layer height. It is given by the following equation for a neutral or stable boundary layer:

(28)

If the current altitude is in the outer layer (z ≥ 0.1h), then, according to this model, the EDR takes the following expression:

(z) = u 3 ∗ kz 1.24 + 4.3z L 1 − 0.85z h 1.5 (42)

6.2.2 Unstable boundary layer

If the current altitude is in the surface layer (z ≤ 0.1h), then the EDR is given by:

(z) = u 3 ∗ kz 1 + z L 2/32/3 (43)

On the contrary, if the current altitude verifies z ≤ 0.1h, then: (z) = w 3 ∗ h 0.8 − 0.3z h (44)

6.3 Results

The output of the software, with EDR calculated at 70 meters by the LIDAR is given on the following plot for a stable atmospheric boundary layer.

Figure 11: EDR profile in calm atmosphere (Heat flux is negative)

(29)

Figure 12: EDR profile in turbulent atmosphere (Heat flux is positive)

This plot summarizes the two preceding plots and shows the different EDR profiles over the day. The black line at 70 meters shows the altitude where the LIDAR input is given. According to this model, the EDR decreases with the altitude.

(30)

(31)

7 Initial vortex separation distance

7.1 Methodology

To improve the wake-vortex model further, the modeling of the initial distance bet-ween wake-vortices can be improved. b0 is the initial vortex separation in meters. It

is different from the wingspan since the vortices need some "roll-up" time. So far, it was assumed that this distance was computed by the following equation:

b0= B × s (45)

where

s = π 4

B = wingspan of the considered aircraft

This equation is valid for an elliptical distribution, which was assumed for every plane. This value is going to be tested for the Boeing 777-300ER aircraft, which is a common one at the considered airport.

There are two ways to verify if this hypothesis stands. Both were explored by Delisi in [4]. Since I have access to LIDAR data for this thesis, the initial distance between the two wake vortices can be calculated directly. Some bias on the measurement is to be expected because of the bias that stems from the LIDAR measurement and because of the fact that the LIDAR only records a position every 8 seconds. The wake therefore had time to move away between its creation and the first LIDAR recording. The second way of calculating b0is through the calculation of V0, the initial vertical

descent velocity of the wake vortices. V0=

gM 2πρU b2

0

(46)

U is the approach speed. Even if it depends on the landing weight of the aircraft, without flight recorded data it will be assumed to have 85% of the Maximum Landing Mass as in [4]. If flight recorded data is available, a very precise calculation can be carried out. For a 777-300ER aircraft this velocity is around 149 kts or 276 km/h depending on the weight of the aircraft.

ρ is the standard atmosphere density and is taken to be 1.2754kg.m−3. The following V0 calculation only stands for a Boeing 77W. For this aircraft type, the

wingspan is 64.8 meters and the maximum landing weight is 251 290 kg.

7.2 Weather assumptions

(32)

Figure 14: Wake created by an A321 bouncing on an inversion layer on November 20th, 2013 10:07:54

Following Delisi, the descent rate is regarded as the slope of the following linear function that represents the position of the vortices against time. This track was produced by a B77W landing in the considered airport on November, 1st 2013. The plot displayed in figure 15 shows how the initial descent velocity is calculated.

(33)

7.3 Results

Each 777 track detected by the LIDAR has its descent rate calculated. This enables to take the median descent velocity over a month. The value b0 is then deduced by

rewriting equation (46). The chosen descent velocity is chosen to be the median value and for november 2013 V0= −1.75m/s. The distance to the median velocity can be

seen on the figure 16.

Figure 16: Distribution of the descent velocities

b0= s gM 2πρU V0 (47) b0= _{9.81 × 213600} 2 × 3.14 × 1.2754 × 76.6 × 1.75 12 = 44.2 m (48)

The value that can be derived for s becomes :

s =44.2

64.8 = 0.68 (49)

(34)

• The precise value of the mass M of the aircraft.

• The value of the approach speed U that is closely linked to the mass of the aircraft.

• The value of the atmospheric pressure ρ.

(35)

8 Conclusion

8.1 Final thoughts

The education I received both at KTH and ISAE-Supaéro was fully relevant to prepare me to carry out engineering projects such as the one presented here. The lectures I attended in aerodynamics, informatics, laser physics and modeling (from a physical and mathematical point of view) were necessary to be able to carry out this master thesis. I want to thank both institutions as well as the European Union education programs that enabled me to get this education.

I particularly enjoyed working on this project because it involved the use of dif-ferent skills. I acquired programming knowledge as well as modeling skills. By parti-cipating in an engineering project I learned how to work under time pressure as well as predict what can go wrong and what can be time-consuming in such a project. In the scope of this project, I learned to use software developed by engineers and researchers working on wake-vortex issues. I met experts working on aviation safety and regulation.

An important achievement of this master thesis was the coding of software that extracts relevant data. One chooses a direction of the beam, a distance away from the LIDAR, and the software extracts the data one chooses (radial velocity to compute the EDR or even other data one might be interested in). While programming this software I acquired knowledge in C++ that I did not possess before.

I want to thank AVTECH Sweden again for all the opportunities I had working on this project. I now feel confident towards working as an engineer in the aviation field or in the aerospace engineering field.

8.2 Further Work

The work undertaken during this master thesis opened several new questions regar-ding this project. To be tuned further, the model will need to be compared against flight recorded data, so that the theoretical strength of the wake can be known more precisely. The evaluation of the strength is a critical quantity if the separation between airplanes is to be lowered.

All the required elements toward a safety analysis are now available. The LIDAR tracks enable to analyse how hazardous the wake can potentially be for the following aircraft and find the most dangerous weather conditions as well as the weather condi-tions where wakes do not pose any danger. Defining a safety case, establishing one or several severity metrics is part of the work that can now be led to use the full potential of this project.

The influence on the value of the parameters T1, T2, ν1and ν2 defined in section

from the EDR can be studied further and the modeling of their value improved based on the LIDAR measurements.

(36)

9 Annex

9.1 Code of the hexdump function

This function is used to extract binary data and convert them to hexadecimal data. The input argument is a pointer to the data to be extracted, the output argument is a string containing 32 hexadecimal characters. This amount can easily be changed.

string hexdump(const char *pAddressIn) {

unsigned char *pTmp,ucTmp;

unsigned char *pAddress = (unsigned char *)pAddressIn; pTmp = pAddress;

int size = 16; char temp[100];

/* Change every two byte to hexadecimal format */ try

{

for(int lIndex = 0; size > 0; size--, lIndex += 2) {

ucTmp = *pTmp++;

sprintf(temp + lIndex, "%02X",(unsigned short)ucTmp); }

}

catch(exception e) {

cout << e.what() << endl; }

int nmber = numberOfCharsInArray(temp); if(nmber == 32)

{

char temp2[33];

/* Change endianness */

(37)

}

/* Turn char to String*/

string temp2bis = turnChar32ToString(temp2); return temp2bis;

} else {

cout << "EOF" << endl; return 0;

(38)

9.2 EDR profile at different times of the day

Figure 17: EDR profile at 2 UTC, stable boundary layer Figure 18: EDR profile at 8 UTC, unstable boundary layer

(39)

9.3 Different LIDAR Tracks

The following tracks are the output of the wake-tracking software. For each plane passing in the LIDAR, wakes are created and the

9.3.1 Example Track

These wake vortices were created by a A340 plane on November 6th. It can be observed that the initial wake strength is about 500 m2_{/s, which is the mean value for a wake created by a heavy (wake category H) aircraft.}

On this track, the two phase decay can be seen perfectly (on figure 25).

Figure 21: Longitudinal Position of the wakes

(40)

Figure 23: Altitude

(41)

(42)

9.3.2 Crosswind

In this case, the wind blows the wakes away the wakes and they are not hazardous for the following plane.

Figure 26: High Crosswind, Longitudinal Position of the wakes

(43)

9.3.3 Tail Wind

This track was produced by a A380-800 on November 4th. It can be seen on figure 32 that the initial strengh of this A380-800 is about 700 m2_{/s but can actually rise until 800 m}2_/s.

(44)

Figure 30: Altitude

(45)

Figure 32: Strength of the wakes

(46)

9.3.4 Bouncing

This track was produced by a A332 on November 2013, 16th. Figure 37 shows that the initial strength of this aircraft is about 400m2_{/s, which is regular for a heavy aircraft.}

(47)

Figure 35: Altitude

(48)

(49)

Bibliography

[1] Pak Wai Chan. Lidar-based turbulence intensity calculation using glide-path scans of the doppler light detection and ranging (lidar) systems at the hong kong international airport and comparison with flight data and a turbulence alerting system. Meteorologische Zeitschrift, 19(6):549–563, 2010.

[2] PW Chan. Generation of eddy dissipation rate map at the hong kong internatio-nal airport based on doppler lidar data. In 12th Conference on Aviation, Range, and Aerospace Meteorology, 2006.

[3] PW Chan, KK Lai, and CM Shun. Measurement of turbulence in terrain-disrupted airflow at the hong kong international airport using a doppler lidar. Hrvatski meteorološki časopis, 40(40):503–506, 2005.

[4] Donald P Delisi, Matthew J Pruis, Frank Y Wang, and David Y Lai. Estimates of the initial vortex separation distance, bo, of commercial aircraft from pulsed lidar data. In 51st AIAA Aerospace Science Meeting, pages 2013–0365, 2013. [5] RG Frehlich and MJ Yadlowsky. Performance of mean-frequency estimators

for doppler radar and lidar. Journal of atmospheric and oceanic technology, 11(5):1217–1230, 1994.

[6] Rod Frehlich. Estimation of velocity error for doppler lidar measurements. Jour-nal of Atmospheric and Oceanic Technology, 18(10):1628–1639, 2001.

[7] Rod Frehlich and Larry Cornman. Estimating spatial velocity statistics with co-herent doppler lidar. Journal of Atmospheric and Oceanic Technology, 19(3):355– 366, 2002.

[8] Rod Frehlich, Stephen M Hannon, and Sammy W Henderson. Coherent dop-pler lidar measurements of wind field statistics. Boundary-Layer Meteorology, 86(2):233–256, 1998.

[9] Rod Frehlich, Yannick Meillier, Michael L Jensen, Ben Balsley, and Robert Shar-man. Measurements of boundary layer profiles in an urban environment. Journal of applied meteorology and climatology, 45(6):821–837, 2006.

[10] Jongil Han, S Pal Arya, Shaohua Shen, and Yuh-Lang Lin. An estimation of tur-bulent kinetic energy and energy dissipation rate based on atmospheric boundary layer similarity theory. National Aeronautics and Space Administration, Langley Research Center, 2000.

[11] Frank Holzäpfel and Robert E Robins. Probabilistic two-phase aircraft wake vortex model: application and assessment. Journal of Aircraft, 41(5):1117–1126, 2004.

[12] Erik Peldan. Computation of air-vortices based on gpu technology: Optimizing and parallelizing a model for wake-vortex prediction using opencl. 2013. [13] Fred H Proctor and David W Hamilton. Evaluation of fast-time wake vortex

prediction models. In 47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition, pages 2009–344, 2009.

Wake Prediction and Modeling Improvement of an existing model