AERODYNAMIC STABILITY ANALYSIS OF QUADCOPTERS SUBJECT TO UNSTEADY

WIND CONDITIONS

by

Mason Slingluff

A thesis submitted to the faculty of

The University of North Carolina at Charlotte
in partial fulfillment of the requirements
for the degree of Master of Science in
Mechanical Engineering

Charlotte

2021

Approved by:

______________________________
Dr. Mesbah Uddin

______________________________
Dr. Chen Fu

______________________________
Dr. Artur Wolek
 ii

©2021

Mason Slingluff

ALL RIGHTS RESERVED

 iii

ABSTRACT

MASON SLINGLUFF. Aerodynamic stability analysis of quadcopters subject to unsteady wind

conditions. (Under the direction of DR. MESBAH UDDIN)

The quadcopter has become an increasingly utilized tool in the military and aerospace

fields over the last several years. For quadcopters to be used in surveillance, search and rescue

and even delivery applications, the quadcopter must be able to perform in unpredictable and

sometimes harsh wind conditions. The environments where drones would be needed can include

wide open spaces where large gusts of wind are possible as well as urban environments that include

tall buildings and other structures that can create chaotic wind patterns. Understanding the

potential wind conditions as well as the effects on the drone’s ability to fly are essential for drones

to be effectively utilized. A quadcopter’s propellers, motors and design structure are all dependent

on the expected forces and conditions that the quadcopter will be subjected to. Using

Computational Fluid Dynamics (CFD), this study investigates the effects of unsteady wind

conditions on a quadcopter’s ability to maintain a hovering scenario. Preliminary isolated rotor

simulations are conducted using the Improved Delayed Detached Eddy Simulation turbulence

model to validate study parameters. Spectral analysis is used to evaluate the meshes of the isolated

rotor to determine which mesh will be used in full scale simulations. Unsteady wind conditions

are replicated in the simulations by using the von Karman atmospheric disturbance model in

accordance with the Department of Defense’s common practice. The thrust coefficient needed to

maintain a stable hovering position from each propeller is recorded at three different turbulence

levels of the von Karman model. The flow structures of the drones in turbulent conditions are

compared to non-turbulent conditions and the changes in the wakes are presented.
 iv

ACKNOWLEDGMENTS

I would like to express my gratitude to my advisor, Dr. Mesbah Uddin, for his guidance

and mentorship through the research and development of this thesis and the University of North

Carolina at Charlotte for allowing me to utilize the computational resources that were necessary

for my research. I would also like to thank my parents and friends for their support throughout my

graduate program. I can confidently say that without the people around me I would not have been

successful in the writing of this thesis. Finally, I would like to thank my committee members, Dr.

Chen Fu and Dr. Artur Wolek, for their time and constructive criticisms to better improve my

research.
 v

TABLE OF CONTENTS

LIST OF TABLES vii

LIST OF FIGURES viii

CHAPTER 1. INTRODUCTION 1

1.1 Objectives of the thesis 1

1.2 Organization of thesis 2

CHAPTER 2. THEORETICAL BACKGROUND 4

2.1 Governing equations of fluid flow 4

2.2 Navier-Stokes equations 5

2.3 Computational fluid dynamics 6

2.4 Background of turbulence modeling 8

2.5 The beginning of turbulence modeling 9

2.6 Improved delayed detached eddy simulation 13

2.7 Turbulence near the wall 15

2.8 Spectral analysis 16

CHAPTER 3. FLOW PATTERNS IN QUADCOPTER FLIGHT 19

3.1 Vortex structures of single rotor flow 19

3.2 Vortex structure of multi-rotor drones 28

3.3 Von Karman atmospheric disturbance modelling 32

CHAPTER 4. RESULTS 34
 vi

4.1 Preliminary URANS simulations 34

4.2 IDDES isolated rotor simulations 44

4.3 Full drone hover 56

4.4 Full drone unsteady wind conditions 60

CHAPTER 5. CONCLUSIONS 75

REFERENCES 79
 vii

LIST OF TABLES

Table 1: Mechanical power, thrust coefficient and figure of merit predictions of mesh

sizes compared to experimental predictions from [4]. 51

Table 2: Mechanical power, thrust coefficient and figure of merit predictions of time-step

sizes compared to experimental predictions from [4]. 53

Table 3: Hovering case drone component lift values. 60

Table 4: Von Karman velocity profile characteristics. 62

 viii

LIST OF FIGURES

Figure 1: Three-dimensional streamline patterns for a view orthogonal to the surface

of the wing [9]: (a) at Φ=36º and (b) Φ=270º. 20

Figure 2: Chordwise sectional cuts of (a,b) spanwise vorticity ωzC/Vrg, (c,d)

spanwise velocity w/Vrg, (e,f) spanwise vorticity flux wωzC/(Vrg2) and (g,h) vertical

(downwash) velocity v/Vrg at angles of rotation (a,c,e,g) Φ=36º and (b,d,f,h) Φ=270º

from [9]. 21

Figure 3: “Free-run” PIV measurement results of the propeller at hover motion from

[14]: (a) instantaneous vorticity, (b) averaged velocity and (c) averaged vorticity. 22

Figure 4: Velocity distributions of the “Phase-locked” PIV measurement results at

different phase angles [14]: (a) 0 degrees, (b) 60 degrees, (c) 120 degrees, (d) 180

degrees, (e) 240 degrees and (f) 300 degrees. 24

Figure 5: Vorticity distributions of the “Phase-locked” PIV measurement results at

different phase angles from [14]: (a) 0 degrees, (b) 60 degrees, (c) 120 degrees, (d)

180 degrees, (e) 240 degrees and (f) 300 degrees. 25

Figure 6: Dynamic process of tip vortices interaction from [14]: (a) 0 degrees, (b) 60

degrees, (c) 120 degrees, (d) 180 degrees, (e) 240 degrees and (f) 300 degrees. 27

Figure 7: Overset grids for a quadrotor system and overset grids for a quadcopter

vehicle system from [15]: (left) quadrotor system and (right) quadcopter vehicle

system. 28

Figure 8: Q-criterion, pressure, and velocity vectors for a quadrotor in hover. Rotor

separation distance=0.3D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]: Velocity

 ix

is projected on the Q-criterion surface and the pressure is projected on the center

plane. 29

Figure 9: Q-Criterion, pressure, and velocity vectors for a quadrotor in hover. Rotor

separation distance=0.5D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]: Velocity

is projected on the Q-criterion surface and the pressure is projected on the center

plane. 30

Figure 10: Q-Criterion, pressure, and velocity vectors for a quadrotor in hover.

Rotor separation distance=0.5D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]:

Velocity is projected on the Q-criterion surface and the pressure is projected on the

center plane. 31

Figure 11: Domain of isolated rotor CFD simulation 35

Figure 12: Thrust coefficient predictions of URANS simulations vs experimental

data from [4]. 36

Figure 13: Mechanical power predictions of URANS simulations vs experimental

data from [4]. 37

Figure 14: Figure of merit predictions of URANS simulations vs experimental data

from [4]. 37

Figure 15: Leading edge and trailing edge vortices profile at 3600 RPM. 39

Figure 16: Leading edge and trailing edge vortices profile at 5400 RPM. 39

Figure 17: Leading edge and trailing edge vortices profile at 7200 RPM. 40

Figure 18: Mean vorticity of 3600 RPM case. 41

Figure 19: Mean vorticity of 5400 RPM case. 42

Figure 20: Mean vorticity of 7200 RPM case. 43

 x

Figure 21: Mesh domain of 40 million mesh size. 45

Figure 22: (a) vorticity of wake at 20 million cells, (b) vorticity of wake at 40

million cells, (c) vorticity of wake at 60 million cells. 46

Figure 23: Q-criterion of isolated rotor. 47

Figure 24: Energy spectrum in focus region of 60 million cell mesh size. 48

Figure 25: Energy spectrum in focus region of 40 million cell mesh size. 49

Figure 26: Energy spectrum in focus region of 20 million cell mesh size. 50

Figure 27: (a) vorticity of wake at 2-degree time step, (b) vorticity of wake at 3-

degree time step, (c) vorticity of wake at 4-degree time step. 52

Figure 28: Energy spectrum in focus region of 3-degree time-step size. 54

Figure 29: Energy spectrum in focus region of 4-degree time-step size. 55

Figure 30:Mesh domain of full drone simulations. 57

Figure 31: Aerial view of full drone simulations. 58

Figure 32: Vorticity of full-drone simulation in hovering condition. Where left is

the vorticity along the nearest rotor plane, the center is the vorticity along the center

plane of the drone and the right is the vorticity along the furthest moth rotors. 59

Figure 33: Vorticity directly under propellers in hovering condition. 59

Figure 34: Light turbulence von Karman atmospheric turbulence model velocities. 61

Figure 35: Moderate turbulence von Karman atmospheric turbulence model

velocities. 61

Figure 36: Severe turbulence von Karman atmospheric turbulence model velocities. 62

Figure 37: Light turbulence full drone simulation drone moment data. 63

Figure 38: Light turbulence full drone simulation propeller power data. 64
 xi

Figure 39: Light turbulence full drone simulation propeller thrust coefficient data. 65

Figure 40: Vorticity of full-drone simulation in light turbulence condition. Where

left is the vorticity along the nearest rotor plane, the center is the vorticity along the

center plane of the drone and the right is the vorticity along the furthest moth rotors. 65

Figure 41: Vorticity directly under propellers in light turbulence condition. 66

Figure 42: Moderate turbulence full drone simulation drone moment data. 67

Figure 43: Moderate turbulence full drone simulation propeller power data. 68

Figure 44: Moderate turbulence full drone simulation propeller thrust coefficient

data. 68

Figure 45: Vorticity of full-drone simulation in moderate turbulence condition.

Where left is the vorticity along the nearest rotor plane, the center is the vorticity

along the center plane of the drone and the right is the vorticity along the furthest

moth rotors. 69

Figure 46: Vorticity directly under propellers in moderate turbulence condition. 70

Figure 47: Severe turbulence full drone simulation drone moment data. 71

Figure 48: Severe turbulence full drone simulation propeller power data. 72

Figure 49: Severe turbulence full drone simulation propeller thrust coefficient data. 72

Figure 50: Vorticity of full-drone simulation in Severe turbulence condition. Where

left is the vorticity along the nearest rotor plane, the center is the vorticity along the

center plane of the drone and the right is the vorticity along the furthest moth rotors. 73

Figure 51: Vorticity directly under propellers in moderate turbulence condition. 74

 1

CHAPTER 1. INTRODUCTION

The topic of this thesis is the numerical investigation of flow surrounding a quadcopter that

is experiencing oscillating wind conditions, created using the von Karman wind turbulence model,

while maintaining a hovering scenario. Quadcopter drones are being used in the commercial

industry as a hobby for enthusiasts as well as in military applications. Quadcopters are capable of

being utilized in search and rescue missions as well as surveillance and reconnaissance.

Understanding the flight dynamics in realistic flight conditions is important to properly utilize

drones in real world scenarios. Quadcopters can be subject to changing wind conditions and large

wind gusts that can alter the trajectory of the quadcopter. Understanding these potential wind

conditions as well as the effects on the flight dynamics of the quadcopter is the first and most

important step in developing a quadcopter. The crucial aspects of the flow surrounding a

quadcopter are the thrust generated by the propellers as well as the power requirements to maintain

thrust. The quadcopter that will be evaluated is the DJI Phantom 3 which is a drone that has been

extensively evaluated by NASA [1,2].

Efforts to develop quadcopter technologies are being carried out by both the United States

Army and NASA. “The Quadrotor Guidance, Navigation, and Control project was established to

develop and apply these [advanced flight control and obstacle field navigation] advanced

technologies to a small-scale commercial-off-the-shelf (COTS) multi-rotor UAS to support

emerging DoD and industry needs” [3]. NASA has conducted many simulations over the DJI

Phantom 3 and published works that show the effects of rotor RPM and rotor placement relative

to the body of the quadcopter.

1.1 Objectives of the thesis

The objectives of the thesis are listed below:

 2

• To conduct a validation study of isolated rotor flows and compare the results to existing

experimental data performed by NASA [4]. The validation study will compare the thrust

force and the power generated by the rotor as well as evaluate the effects of mesh size and

time-step lengths.

• To investigate the transient flow surrounding a DJI Phantom 3 quadcopter with unsteady

and turbulent inlet conditions that mimic real world environments. Three inlet conditions

are investigated with varying degrees of turbulence. The turbulence modeling technique

that is used is the IDDES model. The thrust and power ranges needed to maintain a

hovering condition will be presented as well as the observed flow structures in the wake of

the quadcopter.

1.2 Organization of thesis

The current chapter highlights the objectives and the methodology of the work done in the

thesis. The next chapters each have their own focused content that is introduced.

Chapter 2 discusses the theoretical background and introduces some basic equations of

fluid motion. The goal of this chapter is to explain to the reader how CFD uses these equations in

a simulation. After discussing the equations, an overview of turbulence and how it is modeled is

discussed.

Chapter 3 is focused on presenting previous works and providing a summary of the current

state of predicting flows over quadcopters and other similar drones. This chapter will include

verification and validation cases using experimental methods of visualizing flow around isolated

rotors and full-scale drones. Important characterizations of the experiments such as thrust

coefficient, power requirement and figure of merit are presented from previous works. The

importance of the thrust coefficient is that the thrust coefficient is a non-dimensional value of
 3

rotors that allows rotors of differing size and angular velocity to be compared [5]. The power

requirement is a measurement of the rate of energy that is required to maintain the angular velocity

of the rotor. The figure of merit is a ratio of the ideal induced power to the actual required power

of the rotor [6]. The figure of merit, like the thrust coefficient, allows for comparison of rotors of

differing sizes and angular velocities. Wake characteristics are investigated and key factors in

typical rotor flow are presented to validate the simulation results of the current work. Finally, the

von Karman turbulence model used to simulate turbulent wind conditions is discussed.

Chapter 4 contains the work that has been done to model unsteady wind conditions over

the DJI Phantom 3 quadcopter. The first section of this chapter describes the parameters of the

simulations that were performed and why those parameters were chosen. The expected level of

accuracy from an IDDES simulation is explained as well as the results are compared to

experimental data.

The final chapter, chapter 5, draws all major conclusion from the thesis. Objectives of the

thesis are revisited and evaluated.

 4

CHAPTER 2. THEORETICAL BACKGROUND

The focus of this chapter is to review the theories used in CFD. The governing equations

of the Navier-Stokes equations are presented and explained.

2.1 Governing equations of fluid flow

The governing equations of fluid flow are comprised of the Navier-Stokes equations. The

equations are derived from the conservation laws of physics. These laws of physics, the

conservation of mass and momentum, can be described using four equations. The first equation is

known as the continuity equation and is shown in equation (1).

𝜕𝜌 𝜕(𝜌𝑢) 𝜕(𝜌𝑣) 𝜕(𝜌𝑤)

+ + + =0 (1)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

In equation (1), the variables p, t, and 𝜌 denote pressure, time and density, respectively, while u,

v, and w represent velocity components in x, y, z directions, respectively. This equation equates

the rate of increase of mass inside an element with net rate of flow inside an element through the

defined faces of an element therefore satisfying the law of conservation of mass. In equation (1)

the leftmost term is used for compressible flow; but for the work done in this thesis, this term is

assumed to be zero thus, the flow is considered incompressible. The next three equations are

known as the momentum equations and are shown in equations (2), (3) and (4).

𝜕(𝜌𝑢) 𝜕(𝜌𝑢2 ) 𝜕(𝜌𝑢𝑣) 𝜕(𝜌𝑢𝑤) 𝜕𝑝 𝜕𝜏𝑥𝑥 𝜕𝜏𝑦𝑥 𝜕𝜏𝑧𝑥

+ + + =− +[ + + ] + 𝑆𝑥 (2)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑣) 𝜕(𝜌𝑢𝑣) 𝜕(𝜌𝑣 2 ) 𝜕(𝜌𝑣𝑤) 𝜕𝑝 𝜕𝜏𝑥𝑦 𝜕𝜏𝑦𝑦 𝜕𝜏𝑧𝑦

+ + + =− +[ + + ] + 𝑆𝑦 (3)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑤) 𝜕(𝜌𝑤𝑢) 𝜕(𝜌𝑤𝑣) 𝜕(𝜌𝑤 2 ) 𝜕𝑝 𝜕𝜏𝑥𝑧 𝜕𝜏𝑦𝑧 𝜕𝜏𝑧𝑧

+ + + =− +[ + + ] + 𝑆𝑧 (4)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

These equations equate the rate of increase of momentum of fluid particles to the sum of the forces

acting on the particles. This satisfies the law of conservation of momentum. Each fluid element
 5

is affected by pressure p, which is a normal stress component, and the nine viscous terms τij, which

are viscous stress components. An additional term is added known as the source term Si. This

term represents the overall effect of body forces that act on the fluid element.

It is important to note that there is one additional equation derived from the law of

conservation of energy. In this equation, the heat transfer through conduction is made using an

analogy of viscous momentum transfer. Because the work in this thesis does not need the heat

transfer portion of the fluid flow, this equation is not reviewed.

2.2 Navier-Stokes equations

The Navier-Stokes equations are achieved by redefining the nine viscous terms that are

present in the momentum equations. The way that the viscous terms are defined in the Navier-

Stokes equations is by expressing them as functions of the strain rate. The strain rates are

comprised of the linear deformation rate and the volumetric deformation rate. The linear

deformation rate can be related to the dynamic viscosity, μ. The volumetric deformation rate can

be related to the viscosity of the fluid, η. The relationship between the stain rate components and

the viscosity constants, μ and η, are shown in equations (5)-(10).

𝜕𝑢 𝜕𝑢 𝜕𝑣 𝜕𝑤
𝜏𝑥𝑥 = 2𝜇 +𝜂( + + ) (5)
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝑣 𝜕𝑢 𝜕𝑣 𝜕𝑤
𝜏𝑦𝑦 = 2𝜇 +𝜂( + + ) (6)
𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝑤 𝜕𝑢 𝜕𝑣 𝜕𝑤
𝜏𝑧𝑧 = 2𝜇 +𝜂( + + ) (7)
𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕𝑢 𝜕𝑣
𝜏𝑥𝑦 = 𝜏𝑦𝑥 = 𝜇 ( + ) (8)
𝜕𝑦 𝜕𝑥

𝜕𝑢 𝜕𝑤
𝜏𝑥𝑧 = 𝜏𝑧𝑥 = 𝜇 ( + ) (9)
𝜕𝑧 𝜕𝑥
 6

𝜕𝑣 𝜕𝑤
𝜏𝑦𝑧 = 𝜏𝑧𝑦 = 𝜇 ( + ) (10)
𝜕𝑧 𝜕𝑦

Substituting all the function in equations (5)-(10) with the strain rate terms in equations (2)-(4),

the Navier-Stokes equations are achieved. The new equations are shown in equations (11)-(13).

𝜕(𝜌𝑢) 𝜕(𝜌𝑢2 ) 𝜕(𝜌𝑣𝑢) 𝜕(𝜌𝑤𝑢) 𝜕𝑝 𝜕 2𝑢 𝜕 2𝑢 𝜕 2𝑢

+ + + =− + 𝜇 [ 2 + 2 + 2 ] + 𝑆𝑥 (11)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑣) 𝜕(𝜌𝑢𝑣) 𝜕(𝜌𝑣 2 ) 𝜕(𝜌𝑤𝑣) 𝜕𝑝 𝜕 2𝑣 𝜕 2𝑣 𝜕 2𝑣

+ + + =− + 𝜇 [ 2 + 2 + 2 ] + 𝑆𝑥 (12)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑤) 𝜕(𝜌𝑢𝑤) 𝜕(𝜌𝑣𝑤) 𝜕(𝜌𝑤 2 ) 𝜕𝑝 𝜕 2𝑤 𝜕 2𝑤 𝜕 2𝑤

+ + + =− +𝜇[ 2 + + ] + 𝑆𝑥 (13)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑧 𝜕𝑥 𝜕𝑦 2 𝜕𝑧 2

The Navier-Stokes equations in this form represent Newtonian fluids. The equations for non-

Newtonian fluids differ from equations (11)-(13) because non-Newtonian fluids do not follow the

assumption that viscous stresses are proportional to rates of deformation like Newtonian fluids.

The fluid that will be simulated in the present work follows the Newtonian fluid assumption.

2.3 Computational fluid dynamics

Computational fluid dynamics is a numerical approach using computational methods to

analyze a given fluid flow. There are many variations of CFD software on the market currently

but the one utilized in the present thesis is Star CCM+ version 2020.2. The respective software’s

solver technique is used to approximate variables such as the velocity and vorticity at given

locations in the fluid flow. The predictions of the numerical approach are then compared to

experimental results or higher fidelity CFD simulations to determine the authenticity and accuracy

of the approximation.

The workflow in developing a CFD simulation include many steps in setting up the domain

of the simulation, choosing the solving method, and visualizing and organizing the data that is

solved. The first step in developing a CFD simulation is creating a region that the computational
 7

domain can solve. In an external fluid flow, like the one studied in the present thesis, a sufficiently

large domain is chosen and given parameters of where the fluid flow will enter from and where

the fluid flow will leave as well as the direction that the fluid flow will be traveling. The domain

that is created is then spit up into discrete elements and is now called the mesh. The size of the

mesh is determined based on many different factors such as the turbulence model utilized, the local

distance from a wall in the domain and the Reynolds number of the flow, to name a few. Once the

mesh is created the simulation begins and the solver starts to calculate the target variables. This

process begins by transforming the partial differential equations that are the governing equations

of fluid flow, equation (11)-(13), into algebraic equations that can be solved rapidly using linear

algebra. There are four methods of converting the partial differential equations into algebraic

equations known as discretization. The four methods are Spectral method, Finite volume method,

finite element method, and finite difference method. The discretization method utilized in the

current thesis was the finite element method. The computer then solves the algebraic equations

continually until the solution is considered converged. When the simulations are completed, the

task of post-processing begins. The objective of post-processing is to organize the data gathered

from the simulations into plots and visualizations that characterize the observed flow.

The workflow described in the last paragraph is a very generalized overview of the process

of creating and simulating a fluid flow. Where research is focused on is development of and

analyzation of turbulence models used in this process. How mesh density affects given solving

methods and turbulence modeling methods is another large area of research in CFD as well. So,

while the overview of the workflow may seem relatively simple, the abundance of research

opportunities in each small area of the workflow is immense.

 8

2.4 Background of turbulence modeling

To understand turbulence modeling, one must first understand what turbulence is.

Turbulence is characterized by random fluctuating motions in fluid flow. Fluid flows that have a

large amount of turbulence are flows with high Reynolds numbers. The Reynolds number is

defined as the ratio of momentum effects to viscous effects. This means that as the Reynolds

number of the flow increases, the effects of momentum on the fluid flow are more dominate and

as the Reynolds number of the flow decreases, the viscous effects begin to dominate the flow.

When the viscous terms have a relatively larger effect on flow, the flow is considered laminar.

When the momentum effects are very dominate, the flow is considered turbulent. Turbulence is

generally created by a disturbance in the flow such as an object that the flow must interact with.

Turbulent fluid flow around an object will begin to create swirling motions that are known as

eddies. These eddies can be characterized by their length, time, and velocity scales. While the

flow is random and chaotic, a trend can be noticed in the transfer of energy between eddies. The

trend that appears is that energy flows from the largest eddies, the ones that get their energy from

the mean flow to the smaller eddies. This continues until the eddies are so small that the energy in

these very small eddies dissipates away into internal energy. Turbulent flows are much more

complex because of these eddies and how they transfer energy. This makes obtaining a solution

for turbulent flow significant more difficult to predict than a laminar flow.

Moving forward into how turbulence can be modelled, some of the most widely used

assumption and the most important breakthroughs in understanding turbulence should be

examined. The first assumption, dating back to 1877, that is used very often in turbulence

modelling is the Boussinesq assumption. The Boussinesq assumption states that eddies have a

viscosity value that is linearly proportional to the mean strain rates observed in turbulent fluid
 9

flow. This hypothesis allows turbulence models to calculate the eddy viscosity rather than the

strain rate terms to significantly reduce the computational efforts to reach a converged solution.

How the eddy viscosity value is determined varies between different turbulent modelling

techniques. The next breakthrough that was significant was by Kolmogorov in the 1940s.

Kolmogorov introduces the idea of the 2/3 law which gave way to the -5/3 decay rate that is used

to predict the rate of energy dissipation relative to eddy size. Kolmogorov’s theory suggested that

between the large energy producing eddies and the smallest dissipating eddies, there is a sub-range

where energy is cascaded down between eddies at a rate of -5/3. This trend can be noticed in all

turbulent flows and leads to another major revelation in CFD. The theory behind the Kolmogorov

hypothesis then leads to the formation of the Kolmogorov scale. The Kolmogorov scale is the

smallest scale that eddies can be and still be part of the energy cascade. Any scale smaller than

the Kolmogorov scale does not follow continuum physics and does not need to be resolved. These

scales were developed based off dimensional analysis from the first theory that Kolmogorov

hypothesized. The turbulence model known as the k-ω model is based off the Kolmogorov

hypothesis.

2.5 The beginning of turbulence modeling

The amount of research in turbulence models began to increase dramatically in the 70s and

80s as computational resources and techniques began to advance. The first modern day turbulence

model that was developed was the Large-Eddy Simulation (LES). This model was created by

Deardorff in 1970. The next model developed was the Direct Numerical Simulation (DNS)

developed in 1972 by Orszag and Patterson. Unfortunately, both these models were extremely

computationally expensive in nature and were not able to be utilized on complex flow. This holds

true in even the current time of this thesis. This computational cost caused the development of
 10

many variations of RANS models to be developed in the early 1970s era as well. The first of the

RANS models developed was the k-ε model developed by Launder and Spalding in 1972. This

model was improved upon in 1974 to become the standard k-ε model. The next RANS model was

originally developed in 1970 but was improved upon by Wilcox. This model was known as the k-

ω model. The most recent improvement of the k- ω model is still heavily used in current research.

It is known as the Shear Stress Tensor (SST) k-ω model. This model is known to have high

prediction accuracy in the viscous sublayer and in adverse pressure gradients. The k-ε model and

the k-ω models of the early 1970s mark the beginning of CFD use in fluid dynamics research.

Since then, the development of new and improved turbulence models has been ongoing and is

enjoying the benefits of ever-increasing computational resource availability. In the current times,

modeling turbulence can be done using many different types of turbulence models each with its

own advantages and drawbacks. In general, there are four types of turbulence modeling strategies.

They are known as the RANS models, the hybrid RANS-LES models, the LES models and DNS.

The least computationally expensive turbulence modeling strategy is the RANS modeling strategy.

The RANS model can be broken up into two sub-classes. The Reynolds stress models, and

the eddy viscosity models. What both sub-classes have in common is that they use time-averaging

in their modeling technique. This means that the variables that are being calculated in the grid are

being split into two components: the mean component and the fluctuating component. The

breakdown of the variables is shown in equation (14).

𝑢 = 𝑈 + 𝑢′ (14)

In equation (14), U represents the mean component of the velocity and 𝑢′ represents the fluctuating

component of the velocity. Substituting the mean and fluctuating components into the Navier-

Stokes continuity equation (1) yields the RANS continuity equation shown in equation (15).
 11

𝜕𝑈 𝜕𝑉 𝜕𝑊
+ + =0 (15)
𝜕𝑥 𝜕𝑦 𝜕𝑧

The momentum equations for the RANS modeling technique are shown in equation (16), (17)

and (18).

𝜕(𝜌𝑈 2 ) 𝜕(𝜌𝑈𝑉) 𝜕(𝜌𝑈𝑊)

+ +
𝜕𝑥 𝜕𝑦 𝜕𝑧
(16)
𝜕𝑝 2
𝜕 𝑈 𝜕 𝑈 𝜕 𝑈 2 ̅̅̅̅̅̅
𝜕𝑢 2
′ 𝑢′ ̅̅̅̅̅̅
𝜕𝑢 ′𝑣 ′ ̅̅̅̅̅̅
𝜕𝑢 ′𝑤 ′
=− + 𝜇[ 2 + 2 + 2] − 𝜌[ + + ] + 𝑆𝑥
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑈𝑉) 𝜕(𝜌𝑉 2 ) 𝜕(𝜌𝑉𝑊)

+ +
𝜕𝑥 𝜕𝑦 𝜕𝑧
(17)
𝜕𝑝 𝜕 2𝑉 𝜕 2𝑉 𝜕 2𝑉 ̅̅̅̅̅̅
𝜕𝑣 ′ 𝑢′ ̅̅̅̅̅̅
𝜕𝑣 ′𝑣 ′ ̅̅̅̅̅̅
𝜕𝑣 ′𝑤 ′
=− + 𝜇 [ 2 + 2 + 2] − 𝜌 [ + + ] + 𝑆𝑦
𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑧 𝜕𝑥 𝜕𝑦 𝜕𝑧

𝜕(𝜌𝑈𝑊) 𝜕(𝜌𝑉𝑊) 𝜕(𝜌𝑊 2 )

+ +
𝜕𝑥 𝜕𝑦 𝜕𝑧
(18)
𝜕𝑝 𝜕 2𝑊 𝜕 2𝑊 𝜕 2𝑊 ̅̅̅̅̅̅
𝜕𝑤 ′ 𝑢′ ̅̅̅̅̅̅
𝜕𝑤 ′𝑣′ ̅̅̅̅̅̅̅
𝜕𝑤 ′𝑤 ′
=− +𝜇[ 2 + + ] − 𝜌 [ + + ] + 𝑆𝑧
𝜕𝑧 𝜕𝑥 𝜕𝑦 2 𝜕𝑧 2 𝜕𝑥 𝜕𝑦 𝜕𝑧

This means that variables are being constantly averaged to create the mean value of the computed

variables. This time-averaging also creates new terms known as the Reynolds stresses. The

equations can be altered an additional time to create the URANS equations which are the same as

the RANS equations except with the addition of an unsteady term in the momentum equations

shown in equation (19).

𝜕(𝜌𝑈)
(19)
𝜕𝑡

With both the URANS and RANS methods, how these Reynolds stresses are modeled is what

makes the eddy viscosity models and the Reynolds stress models different from one another. The

eddy viscosity models assume that turbulent stress is proportional to the mean rate of strain. What
 12

separates the eddy viscosity models from each other is the transport equations used to derive the

eddy viscosity term. The Reynolds stress models are much more complicated and require seven

new transport equations to determine all the Reynolds stress created from time-averaging. The

additional transport equations increase the computational expense of the Reynolds stress modeling

technique however because the Reynolds stress model does not assume that turbulent stress is

proportional to the mean rate of strain, the model can better predict complex flows that have

anisotropic turbulence. The next turbulence modeling technique follows the same trend of

increasing computational expense but also increasing expected accuracy of prediction.

The hybrid RANS-LES modeling technique uses both the time averaging technique that is

implement in RANS modeling as well as the space averaging technique that is used in LES. Hybrid

RANS-LES modeling implements the less accurate but more computationally inexpensive time-

averaging technique in flow that is near the wall while using the more accurate space filtering

technique in flow that is away from the wall. The more accurate space filtering technique takes

the time-dependent Navier-Stokes equations and filters them into either Fourier space or

configuration space. This creates new terms known as the sub-grid scale stresses. There are

several different ways to solve the sub-grid scale stresses like the one-equation SGS model or the

dynamic Smagorinsky model to name a couple. The benefit of the hybrid RANS-LES modelling

approach is that it combines the most beneficially aspects of the time averaging technique and the

space filtering technique. The areas where the RANS models are used benefit from the reduced

computational expense of the modeling technique and the areas where LES is used, benefit from

increased accuracy of the predictions. This allows a CFD analyst to custom tailor a CFD

simulation to the available computational resources and desired accuracy of results. As the
 13

availability of computational resources become more abundant, switching from a hybrid RANS-

LES turbulence model to an LES model may become warranted.

The LES uses the space filtering method, like the hybrid RANS-LES model, except it uses

the space filtering method in one hundred percent of the domain rather than in just the far-field

regions. The space filtering method separates the small-scale eddies that are modeled using a sub-

grid scale model and the large-scale eddies that are solved directly.

2.6 Improved delayed detached eddy simulation

The turbulence model that will be utilized in the current study is the Improved Delayed

Detached Eddy Simulation (IDDES) model. This model is considered a hybrid RANS-LES

method that uses the SST k-ω model in the flow near the wall and directly solves the flow in the

outer region of the flow. The SST k-ω model is a RANS turbulence model and an eddy viscosity

model.

The SST k-ω model was created by Mentor in 1994. The purpose of the turbulence model

was to take advantage of the standard k-ε turbulence model’s free stream independence in the far

field region of the flow while also utilizing the standard k-ω turbulence model’s accuracy in the

near to the wall region. How this is achieved is by modifying the standard k-ε model into a k-ω

model formulation by substituting the ε term with ω. Then both models are multiplied by a

blending function. The two models are then added together so that when the flow is near the wall

the blending function returns a value of 1 so that all the k-ω model is utilized and as the flow grows

farther from the wall, the value of the blending function becomes 0 so that all the k-ε model is

utilized. The SST k-ω model also uses a damped cross-diffusion derivative term in the vorticity

equation. The way that the SST k-ω model models the eddy viscosity term is shown in equation

(20).
 14

𝛼∗ 𝑎1
𝜇𝑡 = 𝜌𝑘min , (20)
𝜔 2
2√𝑘 500𝑣
𝑆tanh ((max ( ∗ , )) )
𝛽 𝜔𝑑 𝑑 2 𝜔
( )

In equation (20), k represents the turbulent kinetic energy, 𝛼 ∗ represents a model constant, 𝛽 ∗

represents a model constant, 𝜔 represents the and d represents the distance from the wall. Current

literature suggests that, in comparison to other commonly used turbulence models, the SST k-ω

produces better correlated predictions in adverse pressure gradient flow, especially around airfoils

and transonic shocks. The transport equations for the SST k-ω turbulence model are given by

equations (21) and (22) [7].

𝜕𝑘 𝜕𝑘 2 𝜕𝑢𝑘 2 𝜕𝑢𝑖 𝜕 𝜕𝑘
+ 𝑈𝑗 = (𝜇𝑡 (2𝑆𝑖𝑗 − ) − 𝜌𝑘𝛿𝑖𝑗 ) − 𝛽∗ + [(𝜈 + 𝜎𝑘 𝜐𝑡 ) ] (21)
𝜕𝑡 𝜕𝑥𝑗 3 𝜕𝑥𝑘 3 𝜕𝑥𝑗 𝜕𝑥𝑗 𝜕𝑥𝑗

𝜕𝜔 𝜕𝜔 𝛾 𝜕 𝜕Ω 𝜎𝑤2 𝜕𝑘 𝜕2
+ 𝑈𝑗 = − 𝛽 ∗ 𝜔2 + [(𝜐 + 𝜎𝜔 𝜐𝑡 ) ] + 2(1 − 𝐹1 ) (22)
𝜕𝑡 𝜕𝑥𝑗 𝜐𝑡 𝜕𝑥𝑗 𝜕𝑥𝑗 2 𝜕𝑥𝑗 𝜕𝑥𝑗

The filtering function that determines whether to use a URANS modeling technique versus

directly solving for the flow is like how a standard LES determines whether to use a sub-grid scale

model. The function for the length scale for IDDES formulations is shown in figure (23).

ΔIDDES = min⁡(max(0.15𝑑, 0.15Δ, Δmin ) , Δ) (23)

The filtering function determines whether the SST k-ω model will be used or if the LES model

will be used. The LES model uses the Navier-Stokes equations to directly solve for the flow

outside of the URANS sub-grid scale model region. A more in-depth review of the IDDES model

and its advantages is done by Gritskevich et al [8].

 15

2.7 Turbulence near the wall

The wall has a large impact on the behavior of turbulent flow and therefore constitutes a

unique way of solving this region of flow. A major assumption that is made when dealing with

flow near the wall is that as the distance from the wall decreases, so does the velocity. This

decrease in velocity continues until the flow at the wall is considered to have zero velocity. This

is known as the no-slip wall condition. This no-slip wall condition creates a region of flow that is

dominated by viscous forces known as the boundary layer. The boundary layer is not dependent

on the free-stream parameters of the flow. The parameters that do affect the boundary layer are

the fluid density, fluid viscosity and the wall shear stress. Using dimensional analysis of these

terms’ equations (24) and (25) are derived.

𝑈
𝑢+ =
𝜏 (24)
√ 𝜌𝜔

𝜏
𝜌√ 𝜌𝜔 𝑦
(25)
𝑦+ =
𝜇

In equation (24), 𝑢+ represents the dimensionless velocity value, U represents the velocity vector

and 𝜏𝜔 represents the wall shear stress. In equation (25), 𝑦 + represents the distance in a non-

dimensional form to the wall and y represents the distance to the wall. The formulation of these

terms is known as the law of the wall. Both terms are dimensionless and help determine in which

part of the boundary layer the observed flow is considered. The boundary layer can be broken up

into four main regions. These regions are the viscous sublayer, the buffer layer, the log-law region

and the outer layer. The region closest to the wall is the viscous sublayer. In this region u+ and y+

are equal and this layer ends at y+ equaling five. The flow that is at the wall in the viscous sub-

layer is stationary. The flow in this region does not have any turbulent eddy motions. This region
 16

is also known as the linear sub-layer due to the linear relationship between the velocity of the flow

and the distance from the wall. The next region is the buffer layer. The buffer layer is affected by

both viscous stresses and turbulent stresses. In this region the behavior changes from a linear

relationship between the velocity and the distance from the wall to a log-law relationship that

defines the log-law region of the flow. This transition occurs between y+ equaling five to y+

equaling approximately 60. As the flow surpasses y+ equaling 60 the flow follows the log-law.

The log-law is shown in equation (26).

1
𝑢+ = ln(𝐸𝑦 + ) (26)
𝑘

In equation (26), E is a constant. The log-law is valid for y+ values between 30 and 500; flow

beyond y+>500 is known to be affected by the boundary layer very little and is dominated by

inertial forces and is free from direct viscous effects. This region follows the velocity-defect law

and is also known as the law of the wake. It is common to differentiate between the outer layer of

the boundary layer and all the other regions by grouping all the other regions of the boundary layer

into just the inner region of the flow. This inner region includes the viscous sublayer, the buffer

region and the log-law region.

2.8 Spectral analysis

Spectral analysis of the turbulent kinetic energy of a fluid flow is a popular way of

measuring the complex phenomena that are present on a wide range of scales in the space and time

domain. There is an abundance of experimental research studies that have gathered time histories

of wide ranges of complex flows and then used a form of spectral analysis to decompose the time

histories into modes. The method of decomposition that is of interest in the current thesis is the

discrete Fourier transform. The discrete Fourier transform uses a time history of the turbulent

kinetic energy at a probe, or a group of probes located at a specific point in the domain of the flow
 17

and then converts the time history into an energy spectrum. The energy spectrum that is generated

by the discrete Fourier transform shows a decomposed version of the time history with the x-axis

showing the range of frequencies that are present in the flow and the y axis showing the energy

spectrum at each frequency. The relevance of the energy spectrum is that the energy spectrum

decompositions that have been performed in experimental research show a very distinct profile.

The profile suggests that there are three main frequency groups of a turbulent flow. The low

frequency area of the energy spectrum decomposition shows the energy producing frequencies.

The low frequency magnitude corresponds with eddies of larger size because the period of these

eddies in the time history is much longer requiring a larger amount of space for the eddy. This is

shown by the Energy spectrum being relatively high at the low frequencies, suggesting that the

eddies that are represented by these frequencies have a large amount of kinetic energy. These are

known as the energy producing eddies. The next group of frequencies are in what is known as the

energy cascade. These frequencies are higher than the energy producing eddies but have a negative

trend of energy. The slope of the negative trend has been documented by many experimental

papers as to be -5/3, in agreeance with Kolmogorov’s theory. In the energy cascade region of the

energy spectrum decomposition, the eddies are transferring their energy from the larger eddies,

lower frequency eddies, to the smaller eddies, higher frequency eddies. The energy is transferred

to smaller and smaller eddies until the eddies reach a frequency that is associated with the next

group of frequencies. This group of frequencies is known as the dissipative region. The eddies

that are in this frequency range dissipate their energy away into heat and the energy is converted

from kinetic energy to thermal energy.

Using the profile of the energy spectrum decomposition, the scales at which a CFD

simulation is capturing can be determined. Specifically, for the IDDES turbulence model, the time
 18

history of the flow that is solved for directly in the LES region can be decomposed and the energy

spectrum can be compared to the experimental profile. If higher frequencies of the energy

spectrum are not captured, this can suggest a need for a denser mesh in these regions. If the higher

frequencies are not captured, an underprediction in turbulent kinetic energy will occur and flow

predictions will become inaccurate.

 19

CHAPTER 3. FLOW PATTERNS IN QUADCOPTER FLIGHT

3.1 Vortex structures of single rotor flow

To understand expected flow around rotating rotors, experimental data must be gathered to

compare the simulated data against. Over the last two decades improvements in experimental

measurement devices have allowed substantial advancements on the insight of the vortical

structures of a rotating rotor. Ozen and Rockwell [9] give a good overview of some of the

experimental investigations that have taken place. The techniques that have been used to study

vortex structures include; qualitative smoke, bubble visualization, dye visualization, sectional

imaging involving either chordwise or spanwise planes, and volume imaging with interpretation

of selected sections. Some of the experiments that have been conducted include stereoscopic

particle image velocimetry performed by Zhou et al [10] and Sutkowy et al [11], static load cell

measurement of rotor thrust performed by Ostar-Exel et al [12] and six degree of freedom motion

rig analysis performed by Baris and Britcher [13]. The main flow structure that can be observed

in the abundant number of experiments ranging across many Reynolds numbers are the presence

of a stable leading-edge vortex, vortex breakdown, spanwise velocity magnitude and distribution,

and the history of the rotor’s lift. In Ozen and Rockwell’s [9] experiment, a stable leading-edge

vortex was observed and represented using volume representations. The wing that was observed

was a plexiglass wing that had a chord length of 38.1 millimeters, a span of 78.1 millimeters and

a thickness of 2.8 millimeters. It has an angular velocity of 5.6 rad/s. By using an angular

displacement stereo particle image velocimetry technique, the quantitative imaging was capture.

Streamlines of the lead-edge vortices were captured at 36 degrees and 270 degrees angle of

rotation. These values were chosen through preliminary experimentation to find the most
 20

representative rotation angle that have a constant lift coefficient. The streamlines are shown in

figure 1.

Figure 1: Three-dimensional streamline patterns for a view orthogonal to the surface of the wing
[9]: (a) at Φ=36º and (b) Φ=270º.
The streamlines show a well-defined leading-edge vortex along the span of the rotating wing. At

an angle of rotation of 36 degrees, the swirling motion that represents the leading-edge vortex is

presented along the entire span of the wing as well as the root and tip of the wing. At 270 degrees,

the leading-edge vortex is much larger and is only well defined up to 60 percent of the entire span

of the wing. The root vortex is still present, but the tip vortex is not well defined and does not

have a consistent axis of rotation. Another view of the vortical structure of these two cases is

shown in figure 2.
 21

Figure 2: Chordwise sectional cuts of (a,b) spanwise vorticity ωzC/Vrg, (c,d) spanwise velocity
w/Vrg, (e,f) spanwise vorticity flux wωzC/(Vrg2) and (g,h) vertical (downwash) velocity v/Vrg at
angles of rotation (a,c,e,g) Φ=36º and (b,d,f,h) Φ=270º from [9].
 22

In figure 2(b), the vortical structure near the root of the wing has a short compact profile. As the

distance from the root of the wing increases, the vortical profile begins to elongate and deflect

away from the surface of the wing. The behavior of the leading-edge vortex is a characteristic of

rotating rotor flow that can be observed in all ranges of Reynolds numbers and rotor dimensions.

After the flow has left the surface of the rotor blade it then moves towards the wake of the

rotor. Ning investigated this region of the flow using a high-resolution particle image velocimetry

(PIV) on a rotor in hover [14]. The results of the experiment are shown in figure 3.

Figure 3: “Free-run” PIV measurement results of the propeller at hover motion from [14]: (a)
instantaneous vorticity, (b) averaged velocity and (c) averaged vorticity.
The measurement was done in free-run meaning that the image taking frequency was not the same

frequency as the harmonic motion of the rotor. Figure 3(a) illustrates the instantaneous vorticity

of the wake region. The dark blue circles represent the tip vortices moving downstream with every

revolution. The large red region directly below the motor represents a positive vorticity region
 23

that is due to flow passing the blunt base of the rotor. In figure 3(b), the average velocity of the

flow is represented. In this figure, three major regions can be identified. The first region being

the inflow region located directly above the rotor. In this region, flow is sucked into the rotating

field and then pushed downward into the wake. The difference in velocity between the flow above

the rotor, the inflow region, and the flow directly below the rotor creates a pressure difference that

is the reason why the rotor is producing thrust. The region directly below the inflow region is the

induced flow region. This region is characterized by an increase in velocity and a decrease in

pressure. The decrease in pressure causes a contraction of the wake region and as the flow moves

farther downstream the wake begins to expand again due to dissipation. In figure 3(b), the flow

directly below the tip of the wing has a velocity that is upward right next to the induced flow region

that has a strong downward velocity and a contraction due to the increase in velocity. This

difference in velocity is what creates the tip vortex. The area within the induced flow region with

the highest velocity starts directly below the wing tip at the beginning of the wake and moves

inward as the flow moves downstream. The region outside of the inflow region and the induced

flow region is known as the quiescent region. This region maintains a velocity of less than two

meter per second. In figure 3(c), the averaged vorticity is shown. This figure shows the trail and

dissipation of the leading-edge vortices and the vortices due to the base of the rotor. Both vortices,

the leading-edge and the base, have a very high rate of dissipation that is represented by the

expansion of the vortical region and the dramatic reduction in magnitude as the flow moves

downward. To further analyze the wake region, the velocity and vorticity distributions were

measured at phase-locked intervals. Six intervals were gathered at 0 degrees, 60 degrees, 120

degrees, 180 degrees, 240 degrees and 300 degrees. The distributions are shown in figure 4 and

5.
 24

Figure 4: Velocity distributions of the “Phase-locked” PIV measurement results at different

phase angles [14]: (a) 0 degrees, (b) 60 degrees, (c) 120 degrees, (d) 180 degrees, (e) 240 degrees
and (f) 300 degrees.
 25

Figure 5: Vorticity distributions of the “Phase-locked” PIV measurement results at different

phase angles from [14]: (a) 0 degrees, (b) 60 degrees, (c) 120 degrees, (d) 180 degrees, (e) 240
degrees and (f) 300 degrees.
The major difference that can be observed between the free-run measurements and the phase-

locked measurements is the presence of low velocity pockets shown by the dotted circle in figure

4(a). The cause of these pockets is the rotor blade cutting through the rotation plane and creating
 26

a periodic change in velocity. The area between the pockets has an increase in velocity and then

another decrease when the next pocket is encountered. This pattern of low velocity pockets

continuously moves downstream until the pattern dissipates away. Analyzing the pattern in a

three-dimensional sense, the low velocity pockets act as a sheet that follows the rotor as it rotates.

These low velocity pockets have a lower vorticity than the tip vortices and the vortices due to the

base of the rotor. To better understand the behavior of the wake and the vortex sheet, figure 5 was

indexed with indexes “A” representing the vortices due to the first blade of the rotor and “B”

representing the vortices due to the second blade. The indexes are also numbered numerically

with “1” representing the vortices created from the current rotation, “2” representing vortices from

the previous rotation and “3” representing vortices from the rotation before “2”. Looking at the

tip vortex and the sheet vortex generated by A.1, as the flow moves farther downstream, the vortex

sheet travels faster than the tip vortex. This is due to the vortex sheet being in the induced flow

region and being accelerated downstream to generate thrust. The tip vortex moves downstream

much slower and eventually will combine with the tip vorticity from B.1 like how the tip vortices

from B.3 and A.3 have combined from previous rotations. The velocity at which the vortex sheet

travels downstream is approximately twice that of the velocity of the tip vorticity. Following the

tip vorticity, the slip stream of the wake can be outlined. At 0.7 times the radius of the rotor in the

downstream direction, the flow changes from being dominated by the acceleration of the flow to

being dominated by diffusion. The location where the change from acceleration dominance to

dissipative dominance is known as the threshold plane. The region before the threshold plane is

known as the near-region and the region after the threshold is known as the far-field region. The

far-field region is where the tip vortices mix and dissipate away. The way that the vortices

dissipate is by a roll-up motion. The roll-up motion is shown in figure 6.

 27

Figure 6: Dynamic process of tip vortices interaction from [14]: (a) 0 degrees, (b) 60 degrees, (c)
120 degrees, (d) 180 degrees, (e) 240 degrees and (f) 300 degrees.
The roll up motion is caused by the interaction between the tip vortex and the vortex sheet. The

vortex sheet pulls the tip vortex inward and past the previous tip vortex. This begins the roll-up

motion and contributes to the dissipation in the far-field region. This roll-up behavior is

observed only in hovering scenarios. In an environment where multiple rotors are in close

proximity this type of wake behavior becomes more complex as the wakes of adjacent rotors

interact with one another.

 28

3.2 Vortex structure of multi-rotor drones

When investigating quadcopter wake structures, the location of each rotor, the fuselage,

and the wings play a very crucial role in efficiency and performance. Yoon, Lee and Pulliam

investigated the role that each of the structures play on the performance of a quadcopter by

performing a Detached Eddy Simulation on a XV-15 quadcopter [15]. The XV-15 quadcopter has

rotor diameters of 25 feet. The distance between rotors is 0.2867 diameters from one another but

this aspect of the quadcopter was manipulated to observe the effect of rotor placement. The overset

grid system of the quadcopter is shown in figure 8.

Figure 7: Overset grids for a quadrotor system and overset grids for a quadcopter vehicle system
from [15]: (left) quadrotor system and (right) quadcopter vehicle system.
The simulation was run to thirty rotations to ensure that the flow had reached a quasi-steady state.

The first configuration of the quadcopter was with rotors placed 0.3 diameters away from one

another. The first phenomenon that was different from the isolated rotor flows was the interaction

of vortices above the rotors. This is illustrated in figure 9.

 29

Figure 8: Q-criterion, pressure, and velocity vectors for a quadrotor in hover. Rotor separation
distance=0.3D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]: Velocity is projected on the Q-
criterion surface and the pressure is projected on the center plane.
Vortices from the tip of the rotors are being pulled up and towards the center of the quadcopter.

The velocity vectors of these vortical structure indicate that flow is counter rotational. A center

plane of the pressure distribution is also shown in figure 9. The pressure distribution shows a

propagation of pressure waves. The propagation waves were not present in isolated rotor flow.

This same phenomenon was noticed even when the rotors were separated by 0.5 diameters, shown

in figure 10.
 30

Figure 9: Q-Criterion, pressure, and velocity vectors for a quadrotor in hover. Rotor separation
distance=0.5D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]: Velocity is projected on the Q-
criterion surface and the pressure is projected on the center plane.
The sizes of the vortices were reduced, however. There was still a slight interaction between the

wakes as well as the pressure distribution still showing wave propagations. When the rotors were

separated by one full diameter from one another, the vortical structures vanished, shown in figure

11.
 31

Figure 10: Q-Criterion, pressure, and velocity vectors for a quadrotor in hover. Rotor separation
distance=0.5D, Mtip=0.69, ReCtip=4.9 million, θ=10º from [15]: Velocity is projected on the Q-
criterion surface and the pressure is projected on the center plane.
The only interaction between the rotors was a slight pressure reduction between the rotors. When

the rotor distance was increased to two diameters, the pressure wave propagations disappeared

altogether as well as the wake effects. This illustrates the negative effects on the efficiency of the

rotors relative to the distance from another on the overall behavior of the wake of the quadcopter.

As the distance from the rotors decreased, so did the efficiency of the quadcopter system as a

whole.

The breakdown of the vertical forces of a quadcopter shows an interesting phenomenon as

well. The thrust created by the four rotors while the quadcopter body is present in the domain was

3% higher than that of the thrust generated by four isolated with no quadcopter body in the domain.

The net lift on the quadcopter was lower than the net lift of the isolated four rotors only due to the

negative lift created by the body of the quadcopter. This suggests that the presence of the body of
 32

the quadcopter creates a more desirable wake for the four rotors and increases the thrust that the

rotors can generate. The body generates extremely unsteady fluctuations of pressure which may

be a potential cause of the loss in lift associated with the body of the quadcopter being present in

the domain.

3.3 Von Karman atmospheric disturbance modelling

The von Karman atmospheric disturbance model is the preferred method of replicating real

world atmospheric conditions using an energy spectrum as a basis to calculate the turbulence

velocities. The method of applying the von Karman atmospheric disturbance model is described

in the Department of Defenses’ interface standard flying qualities of piloted aircraft [16]. While

this is the preferred method of the Department of Defense there are several other ways to replicate

turbulence. Davoudi et al reviews several different ways to replicate turbulence in the atmospheric

boundary layer that are more stringent that the von Karman model while Li et al and Cybyk et al

have performed numerical investigation on drone subject to simple low Reynolds number flows

[17,18,19]. The von Karman turbulence spectrum is created by using forming filters and passing

band-limited white noise through the filters. It is used to create continuous turbulent flows that

has qualities that are consistent with the comparable structural analyses of flight. The distribution

of the velocities in the x, y and z directions are dependent upon the altitude and the direction and

speed at which the object of significance is moving through the atmosphere. The energy spectrums

for the velocities are given by equations (27), (28) and (29) [16].

2𝐿𝑢 1
𝜑𝑢𝑔 (Ω) = 𝜎𝑢2 5 (27)
𝜋
[1 + (1.339𝐿𝑢 Ω)2 ]6

8 (28)
2𝐿 𝑣
1 + 3 (2.678𝐿𝑣 Ω)2
𝜑𝑣𝑔 (Ω) = 𝜎𝑣2 11
𝜋
[1 + (2.678𝐿𝑣 Ω)2 ] 6
 33

8 (29)
2𝐿𝑤
1 + 3 (2.678𝐿𝑤 Ω)2
𝜑𝑤𝑔 (Ω) = 𝜎𝑢2 11
𝜋
[1 + (2.678𝐿𝑤 Ω)2 ] 6

In Equation (27), 𝜑𝑖𝑔 represents the spectra function, Ω represents the spatial frequency, σi

represents the turbulent intensity and Li represents the turbulent length scale. How the signal is

generated in the correct characteristics is by passing white noise through the forming filters that

are given in equations (30), (31) and (32) [16].

2 𝐿 𝐿
𝜎𝑢 √𝜋 ∙ 𝑉𝑢 (1 + 0.25 𝑉𝑢 𝑠)
(30)
𝐿 𝐿 2
1 + 1.357 𝑉𝑢 𝑠 + 0.1987 ( 𝑉𝑢 ) 𝑠 2

1 𝐿 𝐿 𝐿 2
𝜎𝑣 √𝜋 ∙ 𝑉𝑣 (1 + 2.7478 𝑉𝑣 𝑠 + 0.3398 ( 𝑉𝑣 ) 𝑠 2 )
(31)
𝐿 𝐿 2 𝐿 3
1 + 2.9958 𝑉𝑢 𝑠 + 1.9754 ( 𝑉𝑣 ) 𝑠 2 + 0.1539 ( 𝑉𝑣 ) 𝑠 3

1 𝐿 𝐿 𝐿 2
𝜎𝑤 √ ∙ 𝑤 (1 + 2.7478 𝑤 𝑠 + 0.3398 ( 𝑤 ) 𝑠 2 )
𝜋 𝑉 𝑉 𝑉 (32)
2
𝐿 𝐿 𝐿 3
1 + 2.9958 𝑉𝑤 𝑠 + 1.9754 ( 𝑉𝑤 ) 𝑠 2 + 0.1539 ( 𝑉𝑤 ) 𝑠 3

In equations (30), (31) and (32), V represents the speed at which the aircraft moves through the

turbulence field and s represents the white noise signal. One of the notable benefits of using the

von Karman method over the Dryden method, which is a common substitute of the von Karman

method, is that the spectral form of the auto correlation function in the von Karman method allows

for a finite microscale. The finite microscale creates a higher proportion of spectral energy at

higher frequencies when compared to other atmospheric disturbance models.

 34

CHAPTER 4. RESULTS

4.1 Preliminary URANS simulations

To better understand the scale of the flow around the DJI Phantom 3 drone, three

simulations were performed on an isolated rotor. The simulations were ran at three different RPMs

on Siemens Star CCM+ platform to determine whether rotor velocity had an impact on the

prediction accuracy of the simulations. The first simulation ran at 3600 RPM, the second was ran

at 5400 RPM, and the third was ran at 7200 RPM. These rates were chosen because of the

abundance of experimental and high-fidelity literature that is available for the DJI phantom rotors

at these rates. The turbulence model selected was the SST k-ω turbulence model. The mesh size

of the simulation was 5.4 million cells for all three of the RPMs. The URANS mesh size is

substantially smaller than the mesh that was generated for the IDDES simulations to utilize the

efficiency of the SST k-ω model. The mesh domain was a cylinder surrounding the isolated rotor,

shown in figure 12.

 35

Figure 11: Domain of isolated rotor CFD simulation

The diameter of the domain was 3 rotor diameters wide and 6 rotor diameters long. The boundary

conditions of the simulation were a stagnation inlet at the top of the domain, the sides of the

cylinder were walls with a slip condition so that a boundary layer does not form along the wall,

and the bottom of the domain was considered a pressure outlet.

The thrust coefficient, power and figure of merit were compared to experimental data

gather by Russel et al [4]. The power requirement of the rotor is obtained by measuring the torque

of the rotor and multiplying the torque by the angular velocity to give the power in Watts. The

power requirement is measured with each time step and then average over the last 20 rotations to

give the final measurement. The thrust is measure by the upward force that is generated on the
 36

rotor. The thrust is measured at each time step and then averaged over the last 20 rotations to give

the final thrust value. The thrust coefficient can then be calculated using equation (33) [20]

𝑇
𝐶𝑇 = (33)
1 2
2 𝜌𝐴(Ω𝑅)

In equation (33), T represents the thrust of the rotor, Ω represents the angular velocity of the rotor,

R represents the radius of the rotor and A represents the reference area of the rotor. The thrust

coefficient is averaged over the last 20 rotations to get the final nominal value. The figure of merit

is found using equation (34) [4].

𝑇 𝑇
𝐹𝑀 = √ (34)
𝑃 2𝜌𝜋𝑅 2

In equation (34), P represents the power of the rotor. The figure of merit is averaged over the last

20 rotations to get the final nominal value. The results of the simulations compared to the

experimental results are shown in figures 12, 13 and 14.

0.025

0.02
Thrust Coefficient

0.015

0.01

0.005

0
3600 4100 4600 5100 5600 6100 6600 7100
RPM

URANS Experimental

Figure 12: Thrust coefficient predictions of URANS simulations vs experimental data from [4].
 37

90
80
70
60
PM (W)

50
40
30
20
10
0
3600 4100 4600 5100 5600 6100 6600 7100
RPM

URANS Experimental

Figure 13: Mechanical power predictions of URANS simulations vs experimental data from [4].

1.2
1
0.8
FM

0.6
0.4
0.2
0
3600 4100 4600 5100 5600 6100 6600 7100
RPM

URANS Experimental

Figure 14: Figure of merit predictions of URANS simulations vs experimental data from [4].
The prediction of the URANS simulations on the thrust coefficient show a dependance on rotor

angular velocity. The lower angular velocity simulation showed an underprediction of thrust while

the middle angular velocity simulation showed an overprediction. The highest angular velocity

simulation predicted the thrust accurately. The power predictions of the URANS simulations show

another dependency on rotor angular velocity. As the simulations increased in angular velocity,
 38

the predicted power requirement began to become underpredicted. This caused an underprediction

in the highest angular velocity case of 24 percent. Due to the errors in prediction of the thrust and

power requirements the figure of merit also suffered. The highest angular velocity case had a 47

percent error in overprediction. This is due to the thrust coefficient matching the experimental

case but the power requirement being underpredicted by a significant amount.

The inaccuracies of the thrust and power prediction are assumed to be mainly due to two

main reasons. The first being the expected inaccuracies of using a lower fidelity turbulence model

such as URANS model. This model is solving for the mean flow and cannot predict the smaller

fluctuations that may affect the power measurements. The second reason is due to the assumption

that the rotor is rigid. This assumption is not physically correct as the rotor will bend due to the

torque of the motor and the resistance of the fluid. It has been observed in previous experimental

studies that the bending of the rotor wings can cause unsteady moments on the rotor that cannot

be predicted by the simulations in the present work. Even with the URANs turbulence model

inaccuracies and the incorrect assumption of rotor rigidity, the simulations were able to maintain

a prediction accuracy of below 24 percent in thrust coefficient predictions and below 47 percent

in figure of merit predictions. The effects of the rotor bending will be present in all simulations

and cannot be eliminated without using a structural model of the rotor and simulating bending.

Adding this aspect to the simulations would greatly increase computational expense. Because of

the increase in computational expense the assumption of the rotor being rigid is maintained

throughout all simulations presented in the current work.

To further evaluate the performance of the URANS predictions, the spanwise leading edge

vortex is compared to figure 1 and figure 2 of the experimental papers. The development of the
 39

leading-edge vortex and the trailing edge vortex in the 3600, 5400 and 7200 RPM ranges are shown

in figures 15, 16 and 17.

Figure 15: Leading edge and trailing edge vortices profile at 3600 RPM.

Figure 16: Leading edge and trailing edge vortices profile at 5400 RPM.
 40

Figure 17: Leading edge and trailing edge vortices profile at 7200 RPM.
The profile of the leading-edge vortex near the root of the rotor is close to the surface of

the rotor and as the location of the leading-edge vortex moves towards the tip of the rotor, the

leading-edge vortex begins to get larger and moves away from the surface. This profile growth is

also seen in figures 1 and 2 of the experimental data of the rotating wing. The trailing edge of the

rotor also begins to develop near the tip of the rotor similarly to the experimental data. This

suggests that the URANS model is capturing the leading-edge vortex and the trailing edge vortex

profiles correctly.

The vorticity structure of the wake of the URANS simulations were captured and the mean

results of the vorticity in the “in-the-page” direction are shown for the three RPMs that were

simulated in figures 18, 19 and 20.

 41

Figure 18: Mean vorticity of 3600 RPM case.

 42

Figure 19: Mean vorticity of 5400 RPM case.

 43

Figure 20: Mean vorticity of 7200 RPM case.

The results of the vorticity in the wake region show a strong vorticity magnitude underneath the

wing tip where the wing tip vorticity is located, like the experimental results shown in figure 3(c).

The contraction of the wake is shown accurately as well as the presence of the vortex sheets that

are being forced downstream in the induced flow region. The low-velocity pockets between the

vortex sheet layers are captured throughout the RPM range as well. The vorticity underneath the
 44

base of the rotor in the “out-of-the-page” direction is capture accurately with a strong magnitude

directly underneath the base of the rotor. With the URANS results resembling the experimental

wake structures that were investigated by Russel et al’s experimental studies, it is assumed that

the wake structure of the rotor is being captured correctly [4]. The next step is to move into higher

fidelity simulations of the isolated rotor case to study the simulation parameters that will be used

in the full-scale drone simulations.

4.2 IDDES isolated rotor simulations

After the URANS simulations were performed the mesh of the IDDES isolated rotor

simulations could be created. The mesh was constructed based on the guidance of Spalart’s paper

[21]. The mesh was broken up into three main regions of flow. The Euler region, the URANS

region and the LES region. The Euler region consisted of the region of flow that was outside of

the rotor’s wake and the rotor itself. This region of the flow had a coarse mesh because there was

not much turbulence in this region and a finer mesh would be computationally inefficient. The

URANS region of the flow consisted of the boundary layer of the flow and some of the viscous

region of the flow. This region of the flow is close enough to the wall of the rotor for the IDDES

turbulence model to utilize URANS modeling rather than LES modeling. The URANS region was

mostly made up of the prism layer which was constructed to meet a wall y+ of approximately 0.1

and a growth rate of 1.2. Once the flow was far enough away from the wall of the rotor the flow

entered the LES region. The LES region had the highest density mesh because in this region there

was a large amount of turbulence. The LES region included the entirety of the wake and eventually

transferred into the Euler region where there was significantly reduced turbulence far downstream

of the wake. This template was used to construct three IDDES meshes with sizes of
 45

approximately 20 million, 40 million and 60 million cells. The mesh with 40 million cells is shown

in figure 21.

Figure 21: Mesh domain of 40 million mesh size.

The three meshes were simulated at 4860 RPM and compared to one another to determine which

mesh size was most computationally effective. The time histories were gathered at each distinct

mesh region to determine how well the scales of flow were being captured at each region’s

location.

The vorticity of the wake regions with each mesh density are shown in figure 22.
 46

(a) (b) (c)

Figure 22: (a) vorticity of wake at 20 million cells, (b) vorticity of wake at 40 million cells, (c)
vorticity of wake at 60 million cells.
The vorticity scale in figure 22 is -500 to 500 s-1. The vorticity predictions show good agreement

with the experimental results that were described in Ning’s experimental paper, figure 5 [14]. The

tip vorticity is well defined in all the meshes that were tested. The wake becomes contracted

directly after the rotor, suggesting that the low-pressure contraction is being accurately captured

by all the meshes as well. The behavior of the tip vortices colliding and dissipating away further

in the wake is also captured. The main difference between the three mesh sizes in the prediction

of the vorticity is that the denser meshes, 40 million and 60 million, capture smaller scale

dissipation of the tip vortices. As the tip vortices collide and the vortices dissipate away, the denser

mesh shows the vorticity stretch and separate into smaller eddies. This is characteristic of eddy

dissipation and is expected in dissipation of energy. The behavior of the colliding of the tip

vortices is shown three-dimensionally in figure 23 of the Q-criterion iso-surface.

 47

Figure 23: Q-criterion of isolated rotor.

The scale for the q-criterion for figure 23 is 100,000 s-2. The vorticity scale in figure 23 is -500 to

500 s-1. As the tip vortex created by the left and right wings move downward, they begin to collide

and form one large vortex that then begins to dissipate away. The magnitude of vorticity decreases,

shown by the dark blue color, as the large vortex continues to move downstream from the rotor.

The energy spectra for the three mesh cases in the focus region of the mesh are shown in

figures 24, 25 and 26.

 48

Figure 24: Energy spectrum in focus region of 60 million cell mesh size.
 49

Figure 25: Energy spectrum in focus region of 40 million cell mesh size.
 50

Figure 26: Energy spectrum in focus region of 20 million cell mesh size.
The 60 million mesh has the highest resolution and should therefore capture the largest number of

scales of motion. As the resolution of the solution becomes coarser with the 40 million and 20

million cell meshes, the solutions capture less scales of motion. The number of scales of motion

that are captured is illustrated using the energy spectrum graphs. The densest mesh at 60 million

cells shows higher turbulent kinetic energy in the higher frequency range and therefore indicates

capturing of smaller scales of motion. The next densest mesh, the 40 million cell size mesh, shows

a decrease in the energy spectrum magnitude around the 800 hertz frequency point, indicating that

scales of motion that are at the 800 hertz frequency and higher are not being captured to the same
 51

extent as the 60 million cell size mesh. At the coarsest mesh, the 20 million cell size mesh, the

magnitude begins to drop even earlier than the 40 million cell size mesh. Again, indicating that

less scales of motion are being captured. Two important frequency peaks are captured at 81 hertz

and 41 hertz. Both frequencies are associated with the rotation rate of the rotor. At 4860 RPM,

the rotation rate is 81 rotations per second. This signifies that a wing tip vorticity is being created

above the energy spectrum probe 81 times per second per wing. The 41 hertz frequency signifies

the wing tip vortices of the left and right wings passing the probes. The 81 hertz frequency

signifies the energy that is measured when the wing tip vortices of each wing collide and begin to

dissipate away. The presence of these peaks at all mesh sizes shows that all the mesh sizes are

capturing the wing tip vortices and the dissipation of the vortices at the same rate.

The final criteria in choosing which mesh size to implement in the full drone simulation

were the prediction of the power, thrust coefficient and figure of merit of the three mesh sizes.

The predictions of each of these characteristics are shown in table 1.

Table 1: Mechanical power, thrust coefficient and figure of merit predictions of mesh sizes
compared to experimental predictions from [4].

Experimental 20 Million 40 Million 60 Million

Mechanical Power (W) 22.5 25.1 24.6 25.6

Thrust Coefficient 0.015 0.015 0.015 0.016

Figure of Merit 0.65 0.66 0.66 0.67

The predictions of the simulations were the within 2 percent of the experimental data for the thrust

coefficient and the Figure of Merit. For the Mechanical power, the predictions ranged from 9.3-

13.3 percent of the experimental data. The larger error in the mechanical power predictions is

assumed to be due to some of the simulation assumptions such as the propeller not having mass as
 52

well as the propeller being considered rigid. In the experiment done by Russel et al, the propeller

has inertia due to its mass and the propeller wings bend slightly during rotation [4]. The bending

and inertia are not captured in the simulations and are a source of error in predictions that may lead

to overpredictions in mechanical power. The mesh size that was chosen was the 20 million cell

size due to the accurate thrust and figure of merit predictions. The negative aspect of the 20 million

cell size mesh is that some of the smaller scale motions will not be captured in the full drone

simulations but the drone characteristics and the major flow patterns around the propeller will still

be captured properly.

With the mesh size captured, the next aspect of the simulation to be tested was the time

step. The 20 million cell size mesh was simulated two more times with a 3-degree time-step and

a 4-degree time-step. This was done to compare the effects of the time-step duration on the

accuracy of the simulations. The vorticity scenes of the 2-degree, 3-degree and 4-degree time

steps are shown in figure 27.

(a) (b) (c)

Figure 27: (a) vorticity of wake at 2-degree time step, (b) vorticity of wake at 3-degree time step,
(c) vorticity of wake at 4-degree time step.
 53

The vorticity scale in figure 27 is -500 to 500 s-1. As the time-step grew larger, the tip vortices

predictions became more smeared, and the dissipation of the vortices began to occur earlier. The

wake contraction predictions also suffer as the time-step becomes larger. The wake prediction

results indicate that time-step size has a large impact on predicting the structure of the wake. The

results of the mechanical power, the thrust coefficient and the figure of merit for each time-step

size is shown in table 2.

Table 2: Mechanical power, thrust coefficient and figure of merit predictions of time-step sizes
compared to experimental predictions from [4].

Experimental 2-Degree 3-Degree 4-Degree

Mechanical Power (W) 22.5 25.1 25.1 25.1

Thrust Coefficient 0.015 0.015 0.015 0.015

Figure of Merit 0.65 0.66 0.66 0.65

There was no difference in the mechanical power and thrust coefficient predictions of the three

step sizes with a 1.5 percent difference in Figure of Merit prediction. Thus, capturing of the tip

vortices did not significantly impact the characteristic prediction of the propeller performance.

The energy spectrum of the 3-degree time-step and the 4-degree time-step are shown in figures 28

and 29.
 54

Figure 28: Energy spectrum in focus region of 3-degree time-step size.

 55

Figure 29: Energy spectrum in focus region of 4-degree time-step size.

Both the 3-degree time-step and the 4-degree time-step maintain the same profile as the 2-degree

time-step indicating that the time-step range did not have a dramatic impact on the scales that were

captured. The overall effects of the time-step changes were that the prediction of the wake

structure suffered greater with little to no impact on the energy spectrum and the propeller

performance characteristics. Because of the importance of the wake structure, the smallest time-

step size of 2 degrees per second was selected. Following the mesh size studies and the time-step

studies, the full-scale drone mesh was ready to be created.

 56

4.3 Full drone hover

Before simulating the full drone in turbulent conditions, the drone was simulated in a

hovering scenario with no turbulence present. The hovering scenario is imposed by having the

body of the drone in a fixed location in the domain with no movement in the x-, y- or z-directions.

The propellers are also in a fixed location except for the rotation axis on the upward direction.

This locked location assumption is slightly different from real world scenarios as the drone would

tilt slightly forward to counteract the wind movement in the domain. Using the coarsest mesh size

template of 20 million cells and the smallest time-step size of 2 degrees per second, the mesh for

the full-scale drone was created and consisted of approximately 50 million cells. The mesh is

shown in figure 28.

 57

Figure 30:Mesh domain of full drone simulations.

The mesh uses the same template as the isolated rotor cases used with focus regions and departure

regions identified for region of refinement. The refinement regions were extended downstream of

the drone to better capture some of the wake of the drone that will be cause by the von Karman

wind turbulence entering from the left side of the domain. The overset regions of the propellers

were reduced compared to the isolated rotor cases so that they do not overlap and cause errors in

the simulations as well as allow for more space in the background mesh to predict the high vorticity

regions that are between the rotor blades. The aerial view of the mesh is shown in figure 29.
 58

Figure 31: Aerial view of full drone simulations.

The drone body is encircled by the focus region mesh and the wake of the drone and the propeller

has an extended departure region. The hovering simulation was ran for 5400 time-steps. Each

time step is 2 degrees per second giving the entire hovering simulation a total simulation time of

30 propeller rotations. The first 10 rotations were simulated using only 1 inner iteration time-step

to quickly initialize the solution while the last 20 rotations used an inner iteration time-step of 7

iteration so that the solution would accurately converge. The 7-iteration time-step was selected

based on preliminary testing that showed that the 7-iteration time-step was the most efficient size.

The results of the vorticity in the wake of the drone are shown in figure 32.
 59

Figure 32: Vorticity of full-drone simulation in hovering condition. Where left is the vorticity
along the nearest rotor plane, the center is the vorticity along the center plane of the drone and
the right is the vorticity along the furthest moth rotors.
The vorticity scale in figure 32 is -500 to 500 s-1. The tip vortices are well defined under the

outer wing rotation location. The vortices move downward and combine with the other vortices

as they dissipate away like the isolated rotor cases. The body of the drone creates small,

detached vortices of their own in the section of the body directly under the propeller. Once the

wing tip vortices had combined, the wake broke down and became chaotic. The vorticity

directly under the propellers of the drone is shown in figure 33.

Figure 33: Vorticity directly under propellers in hovering condition.

 60

The vorticity scale in figure 33 is -500 to 500 s-1. The area between the four rotor shows a large

amount of chaotic vorticity as the leading and trailing tip vortices collide and mix with one

another. This area is directly above the drone body and all the turbulent flow created in this

region is confined between the propellers and the drone body. The drone body also causes drag

which affects the lift of the drone. The propeller lifts and the drone body lift are shown in table

2.

Table 3: Hovering case drone component lift values.

Drone Body (N) Propeller 1(N) Propeller 2 (N) Propeller 3 (N) Propeller 4 (N)

Hover -1.00 3.82 3.82 3.81 3.82

The propeller thrusts are lower than the thrust generated in the isolated rotor case. This is

assumed to be because of the body of the drone interfering in the wake of the rotors. The drone

body also creates negative lift which also lowers the amount of total lift of the drone body and

propellers. This trend was also present in Russel et al experimental paper [4].

4.4 Full drone unsteady wind conditions

The full drone simulations used three degrees of turbulence; a 15-knot light turbulence, a

30-knot moderate turbulence and a-45 knot severe turbulence condition. The graphs of the von

Karman atmosphere turbulence models that were used at shown in figures 34, 35 and 36.
 61

Figure 34: Light turbulence von Karman atmospheric turbulence model velocities.

Figure 35: Moderate turbulence von Karman atmospheric turbulence model velocities.
 62

Figure 36: Severe turbulence von Karman atmospheric turbulence model velocities.
The three turbulence models have the same profile with differing degrees of velocity fluctuations.

The mean velocities, RMS velocities and the turbulent intensity of the three velocity profiles are

shown in table 4.

Table 4: Von Karman velocity profile characteristics.

Mean Mean Mean

RMS X RMS Y RMS Z Turbulent
X-Velocity Y-Velocity Z-Velocity
(m/s) (m/s) (m/s) Intensity
(m/s) (m/s) (m/s)

Light
0.205 0.216 0.308 0.305 0.378 0.489 0.837
Turbulence

Moderate
0.410 0.431 0.615 0.611 0.754 0.977 0.930
Turbulence

Severe
0.614 0.645 0.920 0.915 1.130 1.465 0.730
Turbulence
 63

Turbulent intensity is given in equation (35).

√1 (𝑢𝑥′2 + 𝑢𝑦′2 + 𝑢𝑧′2

3 (35)
𝜎=
√𝑈𝑥2 + 𝑈𝑦2 + 𝑈𝑧2

In equation (35), 𝑢𝑖′ represents the fluctuating component of the velocity and Ui represents the

mean velocity component. The turbulent flow is imposed on the left side of the domain after 10

initial revolutions of the drone are performed to properly initialize the flow. After the initial 10

revolutions, the von Karman velocity profiles are introduced into the domain. The propellers

maintain a rotor velocity of 5400 RPM. Each simulation ran for a total of 1.75 seconds. The

power, moment and thrust coefficients are measured throughout the simulation.

The drone moments from the point when the light turbulent flow is introduced to the

domain is shown in figures 37.

Figure 37: Light turbulence full drone simulation drone moment data.
 64

The moments that are created on the drone stay below 0.03 Nm. Both the clockwise and

counterclockwise axis maintain a moment that is fluctuating around 0±0.04 Nm. The reason for

the fluctuations is that as the rotors rotate, the rotor blades move in and out of the turbulent air that

is trapped between all four the rotors, at the drone’s center, and the air that is free to move

downstream, away from the drone body. The power and thrust coefficients of the propellers as the

drone is subjected to the light turbulence is shown in figure 38 and 39.

Figure 38: Light turbulence full drone simulation propeller power data.
 65

Figure 39: Light turbulence full drone simulation propeller thrust coefficient data.
The power requirements and the thrust coefficient stays relatively stable with the turbulent air

created by the light turbulence case. The propeller power maintained a power requirement of 36±2

W. The thrust coefficient maintained a value of 0.016±0.0015. The effects on the wake flow

structure at the end of the light turbulence case is shown in figure 40.

Figure 40: Vorticity of full-drone simulation in light turbulence condition. Where left is the
vorticity along the nearest rotor plane, the center is the vorticity along the center plane of the
drone and the right is the vorticity along the furthest moth rotors.
 66

The vorticity scale in figure 40 is -500 to 500 s-1. The light turbulence has little to no effect on the

wake structure of the drone when compared to hovering scenario in figure 32. The formation of

the wing tip vortices is not impact by the light turbulence either. The wake is also not shifted

downstream and maintains its shape directly underneath the drone. The vorticity directly under

the propellers is shown in figure 41.

Figure 41: Vorticity directly under propellers in light turbulence condition.

The vorticity scale in figure 41 is -500 to 500 s-1. Again, the propellers are not affected by the

light turbulence condition imposed and resembles the hovering scenario shown in figure 33.

Because the drone rotors experienced little impact from the turbulent air, there was no significant

moments on the drone and had no changes in wake structure. The drone is considered stable for

the span of the simulation in the light turbulence case.

The next simulation that was ran was the moderate turbulence case at a wind speed

magnitude of 30 knots. The angular velocity and drone moments from the point when the moderate

turbulent flow was introduced to the domain is shown in figures 42.

 67

Figure 42: Moderate turbulence full drone simulation drone moment data.
The moments that are created on the drone stay below 0.03 Nm like the light turbulence case. The

only major difference is towards the end of the simulation, the clockwise moment dips at an all-

time low to around -0.04 Nm. The power and thrust coefficients of the propellers as the drone is

subjected to the moderate turbulence is shown in figures 43 and 44.

 68

Figure 43: Moderate turbulence full drone simulation propeller power data.

Figure 44: Moderate turbulence full drone simulation propeller thrust coefficient data.
The power requirements and the thrust coefficient stays relatively stable. The effects on the wake

flow structure at the end of the moderate turbulence case is shown in figure 45.
 69

Figure 45: Vorticity of full-drone simulation in moderate turbulence condition. Where left is the
vorticity along the nearest rotor plane, the center is the vorticity along the center plane of the
drone and the right is the vorticity along the furthest moth rotors.
The vorticity scale in figure 45 is -500 to 500 s-1. The moderate turbulence has little to no effect

on the wake structure of the drone like the light turbulence case and has the same structure as the

hovering scenario. The formation of the wing tip vortices is not impact by the moderate turbulence.

The wake is also not shifted downstream and maintains its shape directly underneath the drone.

The vorticity directly under the propellers is shown in figure 46.

 70

Figure 46: Vorticity directly under propellers in moderate turbulence condition.

The vorticity scale in figure 46 is -500 to 500 s-1. The propellers are not affected by the moderate

turbulence condition imposed and resembles the hovering scenario shown in figure 33. The results

of the moderate case were very similar to the light turbulence case suggesting that the turbulence

increase from the light turbulence case to the moderate turbulence case was not enough to affect

the drone’s performance. Because of these results, the drone is considered stable in the moderate

turbulence case as well.

The final simulation that was run was the severe turbulence case with a wind speed

of 45 knots. The drone moments from the point when the severe turbulent flow was introduced to

the domain is shown in figure 47.

 71

Figure 47: Severe turbulence full drone simulation drone moment data.
The large magnitude of the clockwise and counterclockwise moments indicate that the drone was

much more effected by the severe turbulence compared to the light and moderate turbulence levels.

The moment that was created by the turbulence reached a peak of 0.15 Nm along the

counterclockwise axis and a peak of -0.1 Nm along the clockwise axis. The overall moments had

an average of 0±0.15 Nm. The power and thrust coefficients of the propellers as the drone is

subjected to the severe turbulence is shown in figures 48 and 49.

 72

Figure 48: Severe turbulence full drone simulation propeller power data.

Figure 49: Severe turbulence full drone simulation propeller thrust coefficient data.
Both the power and thrust coefficients are affected by the turbulence as well. The power

requirement increased in range of value to 36±5 W. The thrust coefficient range increased with
 73

an average of 0.016±0.003. The increases in the range shows that the propeller’s performance was

affected by the presence of the turbulent wind conditions. The effects on the wake flow structure

at the end of the moderate turbulence case is shown in figure 50.

Figure 50: Vorticity of full-drone simulation in Severe turbulence condition. Where left is the
vorticity along the nearest rotor plane, the center is the vorticity along the center plane of the
drone and the right is the vorticity along the furthest moth rotors.
The vorticity scale in figure 50 is -500 to 500 s-1. The severe turbulence case shows changes in

the downstream portion of the wake with turbulent air moving through the propeller wake and

shifting the lower portion of the wake downstream. Some turbulence moves outside of the induced

wake zone and is present up to 0.5 diameters downstream of the drone. The defined edges of the

wing tip vortices are also dissipated away earlier in the wake when compared to the hovering, light

turbulence and moderate turbulence cases. The vorticity directly under the propellers is shown in

figure 51.
 74

Figure 51: Vorticity directly under propellers in moderate turbulence condition.

The vorticity scale in figure 51 is -500 to 500 s-1. The presence of vorticity along the outer region

of the rotor’s diameter suggests that the wing tip vortices are impacting the formation of the wing

tip vortices and are most likely the reason to the changes in the rotor’s performance. This correlates

with the thrust coefficients being affected and the power requirements of the rotors beginning to

vary with a higher degree. Because there is a larger moment that was created on the drone, the

drone is slightly unstable with effects of the turbulent air still present past the wake zone of the

drone. When turbulence levels reach the wind speed of 45-knot range, the phantom 3 drone is

more unpredictable and is affected by the severe turbulence. The zone surrounding the drone

where turbulent air is present also increases in size because of the changes in drone performance.
 75

CHAPTER 5. CONCLUSIONS

The flight dynamics of a quadcopter experiencing unsteady wind conditions has been

evaluated and presented. Preliminary low fidelity simulations have been performed to validate

computational fluid dynamic simulation techniques follow by small scale high fidelity

computational fluid dynamic simulations and finally full-scale simulations. Background on fluid

dynamics and turbulence modeling techniques that are used in the present paper have been

explained. Previous literature surrounding experimental studies on quadcopter in hovering

scenarios has also been presented and used as a means of validation for the current work’s results.

URANS simulations using the SST k-ω turbulence model were performed on the small-

scale single propeller. The results showed accurate wake structures when compared to

experimental data. The effects of propeller angular velocity were evaluated using three RPMs of

3600, 5400 and 7200. The results showed that higher angular velocity had a negative effect on

prediction accuracy when compared to experimental studies. The maximum error of the prediction

on thrust coefficient, rotor power and figure of merit was 47 percent. There was a clear presence

of wing tip vortices. The wake of the rotor also demonstrated a contraction in the induced flow

region and showed signs of vortex sheets that matched the wake structures of experimental studies.

The URANS simulations validated the domain conditions and the methodology of the simulations

by providing accurate propeller characteristics and wake behavior. Because of the success of the

URANS simulations, the higher fidelity simulations were able to be performed next.

The higher fidelity single rotor simulations utilized the IDDES turbulence model.

Five simulations were performed to evaluate the effects of mesh density and time-step size. The

mesh sizes tested were 20 million, 40 million and 60 million. The densest mesh was able to predict

the dissipation of the wing tip vortices with the most accuracy by capturing the smaller eddies that
 76

were created as the large eddy is transferring energy down the energy cascade. All the mesh sizes

were able to capture the wake structures accurately when compared to the experimental data. The

energy spectrum was captured in all the simulations and showed the critical frequencies of 81 hertz

and 41 hertz that correspond with the rotor’s angular velocity of 4860 RPM. The figure of merit,

thrust coefficient and rotor powers all increased in accuracy when compared to the lower fidelity

simulations. The figure of merit and the thrust coefficients matched almost exactly to the

experimental data and the rotor power measurements maintained an error range of 9.3-13.3 percent

between the three simulations. The larger error in rotor power measurements is due to the

assumption of rotor rigidity and the assumption of 0 mass in the rotor. Because of these

assumptions the inertia of the rotor and the rotor wing bending phenomenon cannot be accurately

captured, affecting the power measurements. The three-mesh study simulation were all performed

at a time step length of 2 degrees per step. Two additional simulations were performed at 3 and 4

degrees per step to evaluate the effects of time step length on the simulation predictions. The

larger time step did not affect the predictions of the figure of merit, rotor power and thrust

coefficient. However, the dissipation of the wing tip vortices was greatly affected. The wing tip

vortices experienced early and dramatic dissipation and smearing, losing their shape and travel

path. Because of the large effect of time step size, the smallest time-step of 2 degrees per second

is recommended and used in the full-scale simulations. The coarsiest mesh was selected as the

template for the full-scale mesh for simulation efficiency and because at this mesh size, the

important frequencies of the energy spectrum are still able to be captured. If computational

resources were more abundant it would be reasonable to use the denser meshes of 40 or 60 million

cells but the simulation times will increase exponentially with the increase in domain cell size so

simulations time duration should always be considered when evaluating mesh sizes. With the
 77

mesh template chosen and the time-step size chosen the full-scale simulations were able to be

performed.

The full-scale simulations included four simulations. The first being a baseline hovering

simulation with no turbulence being subjected to the drone. The follow three simulations subjected

the hovering drone to three degrees of turbulence using the von Karman turbulence model. The

three turbulence levels were a light turbulence level with a wind speed of 15 knots, a moderate

turbulence level with a wind speed of 30 knots and a severe turbulence level with a wind speed of

45 knots. Each of the turbulent cases was initialized with 10 rotor revolutions with no turbulence

in the domain before the turbulent wind was initiated on the left inlet wall of the domain. The light

and moderate turbulent cases showed little effect on the drone’s ability to maintain hover. The

largest moment that was generated during the two simulations was a moment of 0.03 Nm. The

required power levels did not increase or decrease and matched the power levels of the hovering

case. The thrust coefficient did not change in any of the turbulence cases. This shows that the

rotor efficiency was not negatively affected by the turbulent gusts. The severe turbulence case

showed a different result. The largest moment that was created on the drone was 0.15 Nm. The

range of moments that was observed throughout the simulation was increased by 275 percent. The

thrust coefficient was also affected by showing an increase in the range of observed values by 100

percent. The power requirements were affected by the severe turbulence case by an increase in

range of 150 percent. The wake structure was affected by the severe wind gusts by moving slightly

downstream of the wind gusts with small levels of turbulence still present up to 0.5 diameters

downstream of the drone. This suggests that turbulent air can go through the wake of the drone

and affect other drone or aircraft that are behind it. The aerial view of the wake shows vorticity

affecting the development of the wing tip vortices. The vorticity is this area is most likely the
 78

reason for changes in the rotor performance characteristics. Because of the larger moment, the

drone is expected to be unstable in the 45 knots severe turbulence case. With these results the DJI

Phantom 3 drone is predicted to be able to handle light and moderate turbulence level but will

suffer in performance in an environment where severe degrees of turbulence are expected in a real-

world environment.
 79

REFERENCES

[1] Yoon, S., Lee, H.C. and Pulliam, T.H., 2016. Computational analysis of multi-rotor flows.

In 54th AIAA aerospace sciences meeting (p. 0812).

[2] Theodore, C.R., 2018. A Summary of the NASA Design Environment for Novel Vertical

Lift Vehicles (DELIVER) Project.

[3] Cheung, K.K., Wagster IV, J.A., Tischler, M.B., Ivler, C.M., Berrios, M.G., Berger, T. and

Lehmann, R.M., 2017, May. An overview of the US Army Aviation Development

Directorate quadrotor guidance, navigation, and control project. In Vertical Flight Society

Forum (Vol. 73, p. 2).

[4] Russell, C., Jung, J., Willink, G. and Glasner, B., 2016, May. Wind tunnel and hover

performance test results for multicopter UAS vehicles. In AHS 72nd annual forum (pp. 16-

19).

[5] Eraslan, Y., Özen, E. and OKTAY, T., 2020. A Literature Review on Determination of

Quadrotor Unmanned Aerial Vehicles Propeller Thrust and Power Coefficients. In EJONS

X–INTERNATIONAL CONFERENCE

[6] Jimenez Garcia, A. and Barakos, G.N., 2016. CFD analysis of hover performance of rotors

at full and model-scale conditions. Aeronautical Journal, 120(1231), pp.1386-1424.

[7] Patel, Y., 2010. Numerical investigation of flow past a circular cylinder and in a staggered

tube bundle using various turbulence models.

[8] Gritskevich, M.S., Garbaruk, A.V., Schütze, J. and Menter, F.R., 2012. Development of

DDES and IDDES formulations for the k-ω shear stress transport model. Flow, turbulence

and combustion, 88(3), pp.431-449.

 80

[9] Ozen, C.A. and Rockwell, D., 2012. Three-dimensional vortex structure on a rotating

wing. Journal of Fluid Mechanics, 707, p.541.

[10] Zhou, W., Ning, Z., Li, H. and Hu, H., 2017. An experimental investigation on

rotor-to-rotor interactions of small UAV propellers. In 35th AIAA applied aerodynamics

conference (p. 3744).

[11] Sutkowy Jr, M.L., 2018. Relationship between Rotor Wake Structures and

Performance Characteristics over a Range of Low-Reynolds Number Conditions (Doctoral

dissertation, The Ohio State University).

[12] Ostar-Exel, L., 2019. The effects of varying diameter on coaxial propellers for the

propulsion of multirotor systems (Doctoral dissertation, Rutgers University-School of

Graduate Studies).

[13] Baris, E., Britcher, C. P., & Altamirano, G. 2019. Wind Tunnel Testing of Static

and Dynamic Aerodynamic Characteristics of a Quadcopter. AIAA Aviation 2019 Forum.

http://doi.org/10.2514/6.2019-2973

[14] Ning, Z., 2018. Experimental investigations on the aerodynamic and aeroacoustic

characteristics of small UAS propellers.

[15] Yoon, S., Lee, H.C. and Pulliam, T.H., 2016. Computational analysis of multi-rotor

flows. In 54th AIAA aerospace sciences meeting (p. 0812).

[16] Chalk, C.R., Neal, T.P., Harris, T.M., Pritchard, F.E. and Woodcock, R.J.,

1969. Background Information and User Guide for MIL-F-8785B (ASG),'Military

Specification-Flying Qualities of Piloted Airplanes'. CORNELL AERONAUTICAL LAB

INC BUFFALO NY.

 81

[17] Davoudi, B., Taheri, E., Duraisamy, K., Jayaraman, B. and Kolmanovsky, I., 2020.

Quad-rotor flight simulation in realistic atmospheric conditions. AIAA Journal, 58(5),

pp.1992-2004.

[18] Li, C., Dong, H. and Cheng, B., 2020. Tip vortices formation and evolution of

rotating wings at low Reynolds numbers. Physics of Fluids, 32(2), p.021905.

[19] Cybyk, B.Z., McGrath, B.E., Frey, T.M., Drewry, D.G., Keane, J.F. and Patnaik,

G., 2014. Unsteady airflows and their impact on small unmanned air systems in urban

environments. Journal of Aerospace Information Systems, 11(4), pp.178-194.

[20] Ventura Diaz, P. and Yoon, S., 2018. High-fidelity computational aerodynamics of

multi-rotor unmanned aerial vehicles. In 2018 AIAA Aerospace Sciences Meeting (p.

1266).

[21] Spalart, P.R. and Streett, C., 2001. Young-person's guide to detached-eddy

simulation grids (pp. 1-23). National Aeronautics and Space Administration, Langley

Research Center.