Iowa State University
Digital Repository @ Iowa State University
Graduate heses and Dissertations
Graduate College
2010
Conigurational efect on dust cloud formation and
brownout
Sayan Ghosh
Iowa State University
Follow this and additional works at: htp://lib.dr.iastate.edu/etd
Part of the Aerospace Engineering Commons
Recommended Citation
Ghosh, Sayan, "Conigurational efect on dust cloud formation and brownout" (2010). Graduate heses and Dissertations. Paper 11489.
his hesis is brought to you for free and open access by the Graduate College at Digital Repository @ Iowa State University. It has been accepted for
inclusion in Graduate heses and Dissertations by an authorized administrator of Digital Repository @ Iowa State University. For more information,
please contact hinefuku@iastate.edu.
Configurational effect on dust cloud formation and brownout
by
Sayan Ghosh
A thesis submitted to the graduate faculty
in partial fulfillment of the requirements for the degree of
MASTER OF SCIENCE
Major: Aerospace Engineering
Program of Study Committee:
R. Ganesh Rajagopalan, Major Professor
Thomas J. Rudolphi
Bong Wie
Iowa State University
Ames, Iowa
2010
Copyright c Sayan Ghosh, 2010. All rights reserved.
ii
TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iv
NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v
ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vii
CHAPTER 1. OVERVIEW . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.1
Current Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.1.1
Turbulence models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.1.2
Rotor model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.1.3
Dust transport model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
CHAPTER 2. TURBULENT ROTOR FLOW MODELING . . . . . . . . .
7
2.1
Reynolds Average Navier-Stokes (RANS) Equation
. . . . . . . . . . . . . . .
7
2.2
Eddy Viscosity-based Turbulence Model . . . . . . . . . . . . . . . . . . . . . .
9
2.3
2.2.1
Baldwin-Lomax Model . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2.2
Spalart-Allmaras Model . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2.3
k − ǫ Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
2.2.3.1
Standard k − ǫ model . . . . . . . . . . . . . . . . . . . . . . .
14
2.2.3.2
Realizable k − ǫ model . . . . . . . . . . . . . . . . . . . . . .
16
Rotor Source Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
18
2.3.1
18
Rotor source calculation . . . . . . . . . . . . . . . . . . . . . . . . . . .
CHAPTER 3. DUST TRANSPORT MODELING . . . . . . . . . . . . . . .
3.1
Particle Characteristics
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
23
iii
3.2
3.3
3.1.1
Equivalent particle diameter . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.1.2
Forces on particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
3.1.3
Particle terminal velocity . . . . . . . . . . . . . . . . . . . . . . . . . .
29
3.1.4
Particle threshold friction velocity . . . . . . . . . . . . . . . . . . . . .
31
Modes of Particle Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
3.2.1
Dust emission mechanism and modeling . . . . . . . . . . . . . . . . . .
34
3.2.2
Vertical particle flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
Dust Transport Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
3.3.1
Particle Stokes number
. . . . . . . . . . . . . . . . . . . . . . . . . . .
38
3.3.2
Particle volume fraction . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
3.3.3
Particle bulk density (Dust cloud density) . . . . . . . . . . . . . . . . .
39
3.3.4
Particle loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
3.3.5
Dilute versus dense flows and phase coupling . . . . . . . . . . . . . . .
39
3.3.6
Rotorcraft Brownout Model . . . . . . . . . . . . . . . . . . . . . . . . .
41
3.3.6.1
Dust Transport Equation . . . . . . . . . . . . . . . . . . . . .
43
3.3.6.2
Dust Particle Velocity . . . . . . . . . . . . . . . . . . . . . . .
45
CHAPTER 4. NUMERICAL METHOD AND ALGORITHM . . . . . . . .
4.1
4.2
4.3
46
General Convection-Diffusion Equation . . . . . . . . . . . . . . . . . . . . . . .
46
4.1.1
Scalar variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
Momentum Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.2.1
Pressure correction equation
. . . . . . . . . . . . . . . . . . . . . . . .
55
4.2.2
The pressure equation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
4.2.3
The SIMPLER algorithm . . . . . . . . . . . . . . . . . . . . . . . . . .
59
Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
59
4.3.1
Inlet boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . .
60
4.3.2
Outlet boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . .
60
4.3.3
Wall boundary condition
. . . . . . . . . . . . . . . . . . . . . . . . . .
61
4.3.4
Computational process . . . . . . . . . . . . . . . . . . . . . . . . . . . .
61
iv
CHAPTER 5. RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1
63
Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
5.1.1
3D incompressible solver . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
5.1.2
Turbulence models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
5.1.3
Rotor model validation
. . . . . . . . . . . . . . . . . . . . . . . . . . .
70
5.1.4
Brownout model validation . . . . . . . . . . . . . . . . . . . . . . . . .
71
5.2
Turbulence Model Comparison for Rotor in Ground Effect . . . . . . . . . . . .
73
5.3
Fuselage Effect on Wake and Ground Signature . . . . . . . . . . . . . . . . . .
75
5.4
Fuselage Effect on Brownout . . . . . . . . . . . . . . . . . . . . . . . . . . . .
89
5.5
Fuselage Shape Effect on Wake and Ground Signature . . . . . . . . . . . . . .
95
5.6
Multi-Rotor Configuration Effect on Brownout . . . . . . . . . . . . . . . . . . 102
CHAPTER 6. CONCLUSION, REMARKS, AND FUTURE WORK
. . .
131
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
133
APPENDIX A. WALL FUNCTIONS . . . . . . . . . . . . . . . . . . . . . . .
139
A.1 Standard Wall Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
A.2 Launder Spalding Wall Function . . . . . . . . . . . . . . . . . . . . . . . . . . 140
A.3 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
A.4 Wall Function for Rough Wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
v
LIST OF TABLES
Table 2.1
Closure coefficient for Baldwin-Lomax Model . . . . . . . . . . . . . .
11
Table 2.2
Closure coefficient for Spalart-Allmaras Model . . . . . . . . . . . . . .
14
Table 2.3
Closure coefficient for Standard k − ǫ Model . . . . . . . . . . . . . . .
16
Table 2.4
Closure coefficient for Realizable k − ǫ Model . . . . . . . . . . . . . .
18
Table 4.1
Coefficient and source term for scalar transport equations . . . . . . .
51
Table 5.1
Rotor properties for helicopter in hover . . . . . . . . . . . . . . . . . .
75
Table 5.2
Rotor properties for helicopter with different fuselage shape . . . . . .
95
Table 5.3
Properties of rotor for multi-rotor configurations . . . . . . . . . . . . . 103
vi
LIST OF FIGURES
Figure 2.1
Aerodynamic forces on blade section . . . . . . . . . . . . . . . . . . .
19
Figure 3.1
Particle size characterization . . . . . . . . . . . . . . . . . . . . . . . .
24
Figure 3.2
Particle size distribution in different regions . . . . . . . . . . . . . . .
25
Figure 3.3
Aerodynamic drag and lift on spherical particle . . . . . . . . . . . . .
26
Figure 3.4
Aerodynamic drag coefficient, CD , as a function of Reynolds number Rep 27
Figure 3.5
Particle terminal velocity, wt , Shao (2000) [24] . . . . . . . . . . . . . .
30
Figure 3.6
Particle threshold friction velocity u∗t Iversen et al. (1976) [26] . . . .
31
Figure 3.7
Modes of particle transport . . . . . . . . . . . . . . . . . . . . . . . .
33
Figure 3.8
Dust emission mechanism from Shao (2000) [24] . . . . . . . . . . . . .
35
Figure 3.9
Dilute vs dense flows Crowe et al. (1998) [25] . . . . . . . . . . . . . .
40
Figure 3.10
Mass concentration of various particle size ranges at rotor tip location
during Sandblaster-2 field test on rotorcraft brownout (from Midwest
Research Institute (MRI) technical report (Chatten (2007) [2])) . . . .
42
Figure 4.1
The control volume for scalar variable φ. . . . . . . . . . . . . . . . . .
47
Figure 4.2
Control volumes of u and v momentum . . . . . . . . . . . . . . . . . .
53
Figure 4.3
Flowchart of Rotorcraft Brownout Model
. . . . . . . . . . . . . . . .
62
Figure 5.1
u and v velocity at center plane in driven cavity at Re = 100 . . . . .
64
Figure 5.2
u and v velocity at center plane in driven cavity at Re = 400 . . . . .
64
Figure 5.3
Velocity magnitude contour plot at center plane of driven cavity at Re
Figure 5.4
= 100 and 400 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
Skin friction coefficient with Reynolds number on flat plate
66
. . . . . .
vii
Figure 5.5
Impinging jet simulation set up . . . . . . . . . . . . . . . . . . . . . .
67
Figure 5.6
Impinging jet velocity profile . . . . . . . . . . . . . . . . . . . . . . . .
67
Figure 5.7
Velocity magnitude contour plot (impinging jet) . . . . . . . . . . . . .
68
Figure 5.8
Reynolds Shear Stress u′ w′ in impinging jet . . . . . . . . . . . . . . .
69
Figure 5.9
Aerodynamic blade load comparisons with Rabbott’s experiment at various Θ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 5.10
70
Particle transport by rotor downwash in wind tunnel experiment (Nathan
et al. (2008) [44]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
Figure 5.11
Dust particle distribution predicted by brownout model
. . . . . . . .
72
Figure 5.12
Ground velocity profile for hovering CH-53E rotor . . . . . . . . . . . .
74
Figure 5.13
Ground friction velocity at rotor height H = 0.5R . . . . . . . . . . . .
76
Figure 5.14
Ground friction velocity at rotor height H = 0.75R . . . . . . . . . . .
77
Figure 5.15
Ground friction velocity at rotor height H = 1.0R . . . . . . . . . . . .
77
Figure 5.16
Ground friction velocity at rotor height H = 1.5R . . . . . . . . . . . .
78
Figure 5.17
Ground friction velocity at rotor height H = 2.0R . . . . . . . . . . . .
78
Figure 5.18
Ground friction velocity at rotor height H = 2.5R . . . . . . . . . . . .
79
Figure 5.19
Ground friction velocity at rotor height H = 3.0R . . . . . . . . . . . .
79
Figure 5.20
Ground friction velocity at rotor height H = 4.0R . . . . . . . . . . . .
80
Figure 5.21
Ground friction velocity at rotor height H = 8.0R . . . . . . . . . . . .
80
Figure 5.22
Velocity magnitude contours at H = 0.5R . . . . . . . . . . . . . . . .
82
Figure 5.23
Velocity magnitude contours at H = 1.0R . . . . . . . . . . . . . . . .
83
Figure 5.24
Velocity magnitude contours at H = 1.5R . . . . . . . . . . . . . . . .
84
Figure 5.25
Velocity magnitude contours at H = 2.0R . . . . . . . . . . . . . . . .
85
Figure 5.26
Ground velocity profile at different radial location for different rotor
heights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
86
Figure 5.27
Wall jet velocity at different radial locations . . . . . . . . . . . . . . .
88
Figure 5.28
Dust cloud formation after 5s for rotor-only and rotor-fuselage configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
viii
Figure 5.29
Visualizing dust cloud characteristics . . . . . . . . . . . . . . . . . . .
Figure 5.30
Dust cloud characteristics seen through rotor center cut plane for rotoronly and rotor-fuselage configurations . . . . . . . . . . . . . . . . . . .
Figure 5.31
92
93
Particle path lines at rotor mid plane for rotor-only and rotor-fuselage
configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
94
Figure 5.32
Four fuselage shapes used . . . . . . . . . . . . . . . . . . . . . . . . .
95
Figure 5.33
Ground friction velocity comparison at H = 0.6R . . . . . . . . . . . .
96
Figure 5.34
Ground friction velocity comparison at H = 1.0R . . . . . . . . . . . .
97
Figure 5.35
Ground friction velocity of different fuselage cross-section at H = 1.5R
97
Figure 5.36
Ground friction velocity comparison at H = 2.0R . . . . . . . . . . . .
98
Figure 5.37
Ground friction velocity comparison at H = 3.0R . . . . . . . . . . . .
98
Figure 5.38
Velocity magnitude and streamlines for different fuselage shapes at H =
0.6R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure 5.39
99
Velocity magnitude and streamlines for different fuselage shapes at H =
1.0R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
Figure 5.40
Velocity magnitude and streamlines for different fuselage shapes at H =
1.5R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
Figure 5.41
Rotor configurations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Figure 5.42
Streamlines in lateral plane for tandem rotor hovering at different heights104
Figure 5.43
Streamlines in longitudinal plane for tandem rotor hovering at different
heights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Figure 5.44
Streamlines in lateral plane for tilt rotor hovering at different heights . 107
Figure 5.45
Streamlines in longitudinal plane for tilt rotor hovering at different heights108
Figure 5.46
Streamlines in lateral planes for quad rotor hovering at different heights 109
Figure 5.47
Streamlines on longitudinal plane for quad rotor hovering at different
heights . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Figure 5.48
Ground friction velocity comparison at H = 1.0Rref
. . . . . . . . . . 113
Figure 5.49
Ground friction velocity comparison at H = 1.5Rref
. . . . . . . . . . 113
ix
Figure 5.50
Ground friction velocity comparison at H = 2.0Rref
. . . . . . . . . . 113
Figure 5.51
Ground friction velocity comparison at H = 2.5Rref
. . . . . . . . . . 113
Figure 5.52
Dust cloud comparison for multiple rotor configurations at H = 1.0Rref 114
Figure 5.53
Dust cloud comparison for multiple rotor configurations at H = 1.5Rref 115
Figure 5.54
Dust cloud comparison for multiple rotor configurations at H = 2.0Rref 116
Figure 5.55
Dust cloud comparison for multiple rotor configurations at H = 2.5Rref 117
Figure 5.56
Dust cloud characteristics through fuselage lateral mid plane at H =
1.0Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
Figure 5.57
Dust cloud characteristics through fuselage lateral mid plane at H =
1.5Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
Figure 5.58
Dust cloud characteristics through fuselage lateral mid plane at H =
2.0Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
Figure 5.59
Dust cloud characteristics through fuselage lateral mid plane at H =
2.5Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
Figure 5.60
Dust cloud characteristics through fuselage longitudinal mid plane at
H = 1.0Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
Figure 5.61
Dust cloud characteristics through fuselage longitudinal mid plane at
H = 1.5Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
Figure 5.62
Dust cloud characteristics through fuselage longitudinal mid plane at
H = 2.0Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
Figure 5.63
Dust cloud characteristics through fuselage longitudinal mid plane at
H = 2.5Rref . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
Figure 5.64
Dust cloud evolution after 2s at H = 1.0Rref
. . . . . . . . . . . . . . 123
Figure 5.65
Dust cloud evolution after 3s at H = 1.0Rref
. . . . . . . . . . . . . . 124
Figure 5.66
Dust cloud evolution after 5s at H = 1.0Rref
. . . . . . . . . . . . . . 125
Figure 5.67
Dust cloud evolution after 7s at H = 1.0Rref
. . . . . . . . . . . . . . 126
Figure 5.68
Dust cloud evolution after 10s at H = 1.0Rref . . . . . . . . . . . . . . 127
Figure 5.69
Comparison of tandem rotors with and without fuselage at H = 1.0Rref 128
x
Figure 5.70
Comparison of tilt rotors with and without fuselage at H = 1.0Rref . . 129
Figure 5.71
Comparison of quad rotors with and without fuselage at H = 1.0Rref . 130
xi
NOMENCLATURE
b
rotor blade span
Cd
drag coefficient
Cl
lift coefficient
c
airfoil chord
d m , dp
particle equivalent diameter
g
acceleration due to gravity
k
turbulent kinetic energy
M
local Mach number
p
pressure
p′
pressure correction
Re
Reynolds number
Rep
particle Reynolds number
i
Srotor
rotor momentum source term
Stv
particle Stokes number
ui
fluid velocity components
upi
particle velocity components
Ui
mean fluid velocity components
u′i
turbulent fluid fluctuating velocity components
u′ , v ′
velocities correction
û, v̂
pseudo velocity components
ur
particle relative velocity
u∗
fluid friction velocity
xii
u∗t
particle threshold terminal velocity
V∞
freestream velocity
~
V
velocity vector
~f lap
V
rotor flapping velocity
wt
particle terminal velocity
Greek Symbols
α
angle of attack
αp
particle volume fraction
β
rotor flapping angle
ǫ
turbulent dissipation
Γd
dust diffusion coefficient
µ
fluid viscosity
µt
turbulent eddy viscosity
µpt
turbulent particle eddy viscosity
νt
turbulent kinematic eddy viscosity
φ
general, scalar transport variable (in rotor modeling used as twist angle)
ρ
fluid density
ρp
particle density
ρd
dust cloud density
σp
particle density ratio
θr
rotor blade pitch angle at radius r
τij
Reynolds stress tensor −ρu′i u′j
τv
particle velocity response time
τF
flow field time characteristics
τf
fluid relaxation time
τw
ground shear stress
xiii
ACKNOWLEDGEMENTS
I would like to acknowledge my advisor, Dr. Ganesh Rajagopalan, for his support and
encouragement throughout the entirety of my graduate studies. I would like to thank Air
Force Office of Scientific Research (AFOSR) for providing the financial support to carry out
my research and graduate studies. In addition, I appreciate all the help and discussion with
Dr. William Warmbrodt and Dr. Marvin Moulton and other MURI teams, during the course
of research. I express my deepest gratitude to Dr. James D. Iversen for his suggestion and
guidance. I would also like to thank my other committee members: Dr. Thomas J. Rudolphi
and Dr. Bong Wie. Many thanks are due to the Department of Aerospace Engineering at Iowa
State University and its staff for administrative help and computing resources.
I would also like to thank my friends and colleagues, Angela Lestari, Mark Lohry, Janani
Murallidharan, and Kanchan Guntupalli for their valuable suggestions and support. Most
importantly, I would like to express my deepest appreciation for the support, confidence and
love provided to me by my parents. Without them I would not have made it to where I am
today.
1
CHAPTER 1.
OVERVIEW
Rotorcraft brownout is in-flight visibility restriction caused by clouds of sand and dust
particles during landing, take off, and near ground flight operations in arid desert terrain.
This complex phenomenon is caused by entrainment of dust, sand, and ground particles by
rotor downwash and is compounded by fuselage geometry and its orientation with respect
to the ground. Nonlinear forces and moments on the fuselage, and highly unsteady wind
velocities are common in near ground operations and play a significant role in the behavior of
the particulate clouds that create the brownout condition.
To date, systematic study of brownout is sparse in literature. Some flight test data on
dust distribution for different aircraft are available in Rogers [1] and Chatten [2]. However, in
order to computationally simulate this phenomenon, it is first imperative to correctly predict
the rotor downwash in different flight conditions. Also important is understanding the physics
of particle entrainment and distribution in the brownout dust cloud. Some knowledge can
be obtained from the field of riverine and aeoline sedimentology. It has been found that the
empirical models used in describing the behavior of suspended particles in water or air can be
applied to the helicopter brownout problem also.
Two approaches are popular in studying particle entrainment processes. One is the Lagrangian approach by Keller [3] in which the trajectories of individual particles are tracked and
these particles are taken to represent a dust cloud. The dynamics of each individual particle
is modeled. Thus, for the overall flow to be reliably represented, a large number of particles
must be modeled, thereby making this approach computationally very expensive. In the second approach, the so-called Eulerian approach, overall dynamics of the particle distribution in
the air is modeled using suitable transport equations. This approach has been applied to the
2
brownout problem by Ryerson [4], and Haehnel [5]. All of the above mentioned research use a
two-phase model to represent the particle dynamics and assume a one-way coupling between
the fluid and the particles, i.e., the fluid affects the particle but not the other way around.
Another documented study of the brownout problem was conducted by Phillips and Brown
[6]. They used the Vorticity Transport Model (VTM) to solve the rotor flowfield and coupled
it with an Eulerian-based dust Transport equation to represent the dust entrainment process.
Brown’s dust transport model uses empirical correlations borrowed from sedimentology and
studies the evolution of particle density distribution and its relation with flowfield velocity and
vorticity. This study found the interaction of the tip vortices with the ground is by far the most
important contributor to the formation of the initial dust cloud. Therefore, the vorticity and
velocity profile in the rotor wake are crucial. Phillips et al [7] determined them to be affected
by parameters like solidity and blade twist. From these studies it can be concluded that the
geometric properties of the rotor have a great influence on the type and extent of the dust
cloud formed. Although VTM is found to be a good model in predicting the rotor downwash,
no fuselage effect is considered in this study. Since it is based on empirical correlations, the
dust transport model has room for further investigation.
Another approach, ([8], [9]) employs a real-time free wake model, coupled with a dust and
debris model, that uses analytical methods in a Lagrangian framework. This model results in
a real-time brownout simulator that gives pilots a feel for the brownout cloud formed under
different flight conditions, ground cover and type of aircraft. A visual obscuration model, based
on a light scattering method, is also used to visualize brownout and train pilots accordingly.
However, the Lagrangian approach makes this approach computationally expensive.
One of the important factors in developing a dust transport model that can be used in
conjuction with a Computational Fluid Dynamics (CFD) code is the development of a suitable
entrainment function to describe the entrainment flux. Haehnel [10] studied the physics of
particle entrainment under the influence of an impinging jet and determined the driving force
for particle entrainment is the Reynolds stress associated with the turbulent fluctuations of
the flow. Since an impinging jet is similar to the wake of the rotor in ground effect, it can
3
be deduced that solving correctly for the turbulent fluctuations in the rotor’s flowfield is an
important step in solving the brownout problem. Thus, a complete brownout model that
incorporates all the components of brownout as observed in reality is yet to be developed.
1.1
Current Research
The goal of this current research is to develop a computational test bed to study fuselage
and rotor configuration effects on brownout. The emphasis of the simulation approach has
been on
• global flow field rather than on the flow near the rotor blades.
• body influence on the brownout and on the inception of particle entrainment.
Brownout is a complex phenomenon with multiple scales and physics dictating the physics
of the simulation. The complex phenomenon involves aerodynamic processes like ground effect
on rotor wake, turbulent burst, surface wall jet, tip vortices, and their interactions. These
aerodynamic interactions, coupled with the interaction of the particles, become a complex
multi-phase phenomenon. In addition to the rotor-induced vortical wake, there are unsteady
effects, due to the bluff body underneath the rotor as well as the dynamics of the particles that
complicate the simulation. Hence, to study this behavior accurately as well as economically,
judicious compromises and an engineering approach to the complete physics are required.
Accuracy of the simulation, particularly close to the ground where particle entrainment
takes place, is highly dependent upon the turbulence characteristics and ground shear predicted
by the mathematical models used. Since the phenomenon is inherently three-dimensional, a
three-dimensional flow solver has been developed for studying various turbulence models. Many
variants of turbulence models are considered and implemented in the solver before settling on
the final choice for studying the fuselage configuration effects on rotorwash and brownout. In
addition, a solution procedure for dust cloud evolution and dynamics has been developed.
To accurately and economically predict the turbulent flow field, incompressible Reynolds
Averaged Navier Stokes (RANS)-based solver has been developed in the Cartesian coordinate
4
system. In the solver, the governing flow equations are solved using a finite-volume based
method known as SIMPLER ([11]). In this algorithm, the flow field is determined by solving
for primitive variables, namely the static pressure and the velocity vector, directly from the
mass and momentum conservation equations.
To solve for the Reynolds turbulent stresses, additional relations or partial differential
equations (PDE) are solved, known as turbulence models. In general, these turbulence model
can be broadly classified [12] into Zero Equation models, One Equation models and Two
Equation models. These models are also known as linear eddy viscosity models, as they assumes
a linear relationship between turbulent eddy viscosity and the mean strain rate. Other complex
RANS-based models, such as Reynolds Stress models [13], Non-linear Models [14] etc. can also
be found in the literature. Since there has not been much research on the applicability of these
models on rotorcraft brownout flows, it has been decided to use linear eddy viscosity models
for the current research because of their relatively low complexity and low computational cost.
Although these models are most widely used for many industrial flow problems, one of the
drawbacks of these models is that it assumes isotropic character of turbulence, which may be
erroneous for the anisotropic turbulence condition found in brownout.
1.1.1
Turbulence models
Zero Equation models, as the name denotes, do not solve any partial differential equation.
Instead, turbulent viscosity is solved using an algebraic equation. On the other hand, One
Equation and Two Equation models solve one and two partial differential equations, respectively. In general, complexity and accuracy of linear eddy viscosity models increases with the
number of equations solved. A candidate model from each of these categories has been selected,
implemented, and tested.
Zero Equation Model (Baldwin Lomax Model)
This model implementation is
straightforward for simple single body flows. It gives acceptable results for simple attached
flow with mild separation and re-attachment. The eddy viscosity is calculated as a function
of turbulent length scale (lt ) and turbulent velocity scale (vt ), both are semi-empirically mod-
5
eled. Length scale is of the order of the boundary layer thickness calculated using vorticity
distribution; thus, eliminating the need to calculate the boundary layer thickness explicitly.
However the model does not perform well where turbulence properties are not proportional
to the mean flow length scale, such as in flows with strong separation and recirculation. The
model is computationally very efficient.
One Equation Model (Spalart Allmaras Model)
One equation model overcomes
some of the problems of the Zero Equation Model by solving the transport equation for turbulent viscosity-like variables. Since a transport equation is solved, it captures the effects of
convection and diffusion of turbulence. Thus, the flow downstream is affected by the upstream
turbulence. However, this model still needs to semi-empirically model the length scale to
complete the equations. This is achieved by making the source term (destruction term) as a
function of distance from the nearest wall. This model gives better results for wall-bounded
flows, separated flows, and boundary layer flows with adverse pressure gradients. However,
calculation of an additional equation increases the computational requirement and storage of
additional variables.
Two Equation Model (Standard and Realizable k − ǫ Model)
Solving two equa-
tions avoids the dependence on empirical and semi-empirical algebraic formulations for both
turbulent velocity and length scales. The two coupled partial differential equation are solved
for turbulent kinetic energy, k, and turbulent dissipation, ǫ. This model is widely used for a
variety of turbulent flows involving separation and re-attachment, and flows with complex turbulent characteristics. However, solving two additional equations requires more computational
time and memory.
1.1.2
Rotor model
To model the rotor, the method developed by Rajagopalan and Fanucci, Rajagopalan and
Lim, and Rajagopalan and Mathur [15, 16, 17] has been used, where the rotor’s influence is
considered as a momentum source which imparts momentum to the fluid flow around it. The
6
momentum source is time averaged and its magnitude depends upon the rotor’s geometry and
flow characteristics around the rotor. The advantage of this model is that there is no need
of a body fitted grid around the rotor; thus, reducing the memory and computational time
requirements. In the time averaged model of the rotor used in this study the effect of discrete
blades are not captured.
1.1.3
Dust transport model
The transport of dust has been modeled by the multiphase species transport equation.
Since the number of dust particles during brownout is immense, dust particle transport is
assumed to be in a continuum phase and is solved in the Eulerian frame of reference. The
advantage of modeling dust in such a manner reduces the computational requirements and
memory. Also, the coupling between particle and flow phase has been considered as a one-way
coupled system, where the particles are affected by the flow and the particles do not have
any effect on the flow. The modeling of dust lift-off criteria, entrainment, and saltation are
accomplished through empirical formulation used in geophysical research. The particle velocity
field is calculated through an Algebraic slip model, which calculates the relative particle velocity
using an algebraic expression.
In summary, all of the aforementioned models have been integrated into a single program,
Rotorcraft Brownout Model (RBM). In the present research, RBM has been used to study the
following problems:
1. Study of various turbulence models on rotorcraft problems in ground effect.
2. Effects of fuselage on rotor wake, ground signature and brownout.
3. Different rotor configurations and ground height effect on rotor wake, ground signature,
and brownout.
7
CHAPTER 2.
2.1
TURBULENT ROTOR FLOW MODELING
Reynolds Average Navier-Stokes (RANS) Equation
The equation for conservation of mass and momentum of incompressible rotor flow field is
given by
∂ui
=0
∂xi
(2.1)
∂tij
∂ρui ∂ρui uj
∂p
i
+
=−
+
+ Srotor
∂t
∂xj
∂xi
∂xj
(2.2)
where ρ is density, xi is position, ui is velocity, p is pressure, and t is time. The term tij is
viscous stress tensor given by
tij = 2µSij
(2.3)
where µ is molecular viscosity and sij is the strain rate tensor given by
Sij =
1
2
∂uj
∂ui
+
∂xj
∂xi
(2.4)
i
Srotor
is the rotor momentum source term used to model the rotor (discussed in section 2.3).
The above set of equations form the incompressible Navier-Stokes equation and can be
used to solve for laminar as well as turbulent flows. However, most of the practical turbulent
flows contain cascades of eddies with a wide range of time and length scales. Direct Numerical
Simulation (DNS) is a numerical method, which resolves all the length and time scales of
turbulent motion using very fine grids. Large Eddy Simulation (LES) is another numerical
technique, in which the larger eddies are resolved; whereas, the smaller turbulent structures
are modeled. Nevertheless, in case of complex flows, such as rotor flows, both of the above
methods need an impractical degree of grid refinement and simulation time for design use.
8
An alternate approach is to solve for the time averaged mean flow as introduced by
Reynolds. In this concept, the instantaneous velocity, ui , in an inhomogeneous turbulence
can be expressed in terms of mean velocity, Ui , and fluctuating velocity, u′i .
ui = Ui + u′i
(2.5)
The mean velocity Ui is the time averaged velocity defined as
Ui = lim
T →∞
1
2T
Z
T
ui dt
−T
(2.6)
where T is a time interval large enough compared to the turbulent length scale, but smaller
than the time in which the time-averaged mean flow varies. Also, the time average of the
fluctuating velocity is zero.
lim
T →∞
1
2T
Z
T
u′i
dt = 0
−T
(2.7)
Replacing the instantaneous velocity, ui , in the Navier-Stokes equation with the sum of
time averaged velocity, Ui , and fluctuating velocity, u′i , and then taking the time average of
the equation yields Reynolds Average Navier-Stokes (RANS) equations.
∂Ui
=0
∂xi
(2.8)
∂p
∂
∂ρUi ∂ρUi Uj
i
+
=−
+
2µSij − ρu′i u′j + Srotor
∂t
∂xj
∂xi ∂xj
(2.9)
where −ρu′i u′j is the time average of product of fluctuating velocities and are called Reynolds
Stresses, denoted by τij .
τij = −ρu′i u′j
(2.10)
wherein τij is a symmetric tensor, i.e., τij = τji .
The RANS equations are not closed form equations. The Reynolds stresses, τij , are six
additional unknowns, that need to be calculated to complete the RANS equations. Additional
relationships or equations are required to relate or calculate τij in terms of known or calculated
9
variables. These additional relations or equations are the main essence of turbulence modeling.
2.2
Eddy Viscosity-based Turbulence Model
According to Boussinesq’s hypothesis, the turbulent Reynolds stress is related linearly to
the mean strain-rate tensor of the fluid flow. Analogous to viscous stress, the proportionality factor is called turbulent viscosity or eddy viscosity, µt ,. Boussinesq’s hypothesis for
incompressible flow is
2
τij = 2 µt Sij − ρkδij
3
(2.11)
where the term 23 ρkδij is a linear constitutive relation used while solving two-equation models
or models which solve for turbulent kinetic energy, k. Turbulent kinetic energy is defined as
k=
1 ′2
′
′
u1 + u22 + u32
2
(2.12)
Using Boussinesq’s hypothesis, the RANS (Eq. 2.9) equation becomes
2 ∂k
∂Ui ∂Uj
∂ρUi ∂ρUi Uj
∂p
∂
i
(µ + µt )
−
+
=−
+
+
+ Srotor
∂t
∂xj
∂xi ∂xj
∂xj
∂xi
3 ∂xi
(2.13)
Thus, using the above assumption, the problem of finding Reynolds stress reduces to determination of eddy viscosity. Unlike molecular viscosity, µ, eddy viscosity, µt , depends on local
flow characteristics, mainly the local turbulent length scale and velocity scale. The eddy viscosity is calculated using various eddy viscosity-based turbulence models which can be broadly
categorized, depending upon the number of transport equations being solved.
• Zero Equation Models (Cebeci-Smith Model, Baldwin-Lomax Model)
• One Equation Models (Baldwin-Barth Model, Spalart-Allmaras Model)
• Two Equation Models (k − ǫ Model, k − ω Model)
In the current research, one representative model from each category, i.e., Baldwin-Lomax
Model, Spalart-Allmaras Model, and k − ǫ has been implemented and tested.
10
2.2.1
Baldwin-Lomax Model
The Baldwin-Lomax model [12, 18] is a two-layer algebraic zero-equation eddy viscosity
model, which gives the eddy viscosity, µt , as a function of the local boundary layer velocity
profile. The earlier zero-equation models require determination of the outer edge of boundary
layer to calculate the length scale for boundary layer and wakes. This problem is eliminated
in this model as it uses the vorticity distribution to compute the required length scales. One
of the major disadvantages of this model is it is unable to capture the time history effect of
turbulence, such as convection and diffusion, because no transport equation is being solved.
This model is quite robust and useful for simpler flow geometries.
Equations: The turbulent viscosity is given as
µt =
µ
t inner
if y ≤ ym
µt outer
if y > ym
(2.14)
where ym is the smallest distance from the surface, where µt inner is equal to µt outer . The inner
region is given by the Prandtl - Van Driest formula:
2
|ω|
µt inner = ρlmix
(2.15)
where mixing length lmix is given by
lmix = ky 1 − e
−y +
A+
where y + is non-dimensional height calculated as y + =
|ω| is magnitude of vorticity given by
|ω| =
(2.16)
p
τw /ρ y/ν. A+ is a constant.
p
2ωij ωij
(2.17)
where
1
ωij =
2
∂uj
∂ui
−
∂xj
∂xi
(2.18)
11
The outer region is given by:
µt outer = ρ K Ccp Fwake FKleb (y)
(2.19)
where K and Ccp are constants. Fwake is given as
Fwake = M IN
ymax Fmax ; Cwk ymax
u2dif
Famx
!
(2.20)
where Cwk is a constant.
Fmax =
1
[maxy (lmix |ω|)]
κ
(2.21)
where ymax is y at which lmix |ω| is a maximum.
FKleb is the intermittency factor given by:
"
FKleb (y) = 1 + 5.5
y CKleb
ymax
6 #−1
(2.22)
where CKleb is a constant.
udif is the difference between maximum and minimum speed in the profile. For boundary
layers, the minimum is always set to zero.
√
√
udif = max( ui ui ) − min( ui ui )
The values of closure coefficient used in this model is given in Table 2.1.
Table 2.1
Closure coefficient for Baldwin-Lomax Model
A+
Ccp
CKleb
Cwk
κ
K
26
1.6
0.3
0.25
0.4
0.0168
(2.23)
12
2.2.2
Spalart-Allmaras Model
Spalart-Allmaras Model [19] is a one-equation eddy viscosity model which solves a transport
equation for a viscosity-like variable, ν̃, used to calculate kinematic turbulent viscosity, νt .
Equations: Kinematic turbulent eddy viscosity, νt , is related to ν̃ as:
νt = ν̃fv1
(2.24)
where
fv1 =
χ3
3 ,
χ3 + Cv1
χ=
ν̃
ν
(2.25)
The intermediate variable ν̃ is calculated using the transport equation:
∂ρuj ν̃
∂ρν̃
+
= Sdif f usion + Sproduction − Sdestruction + Strip
∂t
∂xj
(2.26)
where Sdif f usion is a diffusion term, Sproduction is the source term due to turbulence production,
Sdestruction is the sink term due to turbulence destruction, and Strip is a trip term used for
turbulence transition.
The diffusion term is defined as
Sdif f usion =
∂ ν̃
∂ ν̃ ∂ ν̃
ρ ∂
(ν + ν̃)
+ Cb2
σ ∂xj
∂xj
∂xi ∂xi
(2.27)
The production source term depends upon vorticity and allows good modeling near the
wall.
Sproduction = ρCb1 [1 − ft2 ]S̃ ν̃
(2.28)
where Cb1 is a constant. S̃ is given by
S̃ = S +
ν̃
fv2 ,
κ2 d 2
fv2 = 1 −
χ
1 + χfv1
(2.29)
13
d is distance from nearest wall and S is magnitude of vorticity given by
S=
p
2Ωij Ωij ,
∂uj
1 ∂ui
−
)
Ωij = (
2 ∂xj
∂xi
(2.30)
ft2 is given by
ft2 = Ct3 exp(−Ct4 χ2 )
(2.31)
where Ct3 and Ct4 are constants.
The destruction term depends upon the nearest wall distance, d.
Sdestruction
Cb1
= ρ Cw1 fw (r) − 2 ft2
κ
2
ν̃
d
(2.32)
The term fw (r) calibrates the outer region of boundary layer and depends upon the characteristic length r, and is given by:
6
1 + Cw3
fw = g 6
6
g + Cw3
1/6
,
g = r + Cw2 (r6 − r),
r=
ν̃
S̃κ2 d2
(2.33)
where Cw2 and Cw3 are constants.
The trip term is modeled as
Strip = ρft1 ∆U 2
(2.34)
where ∆U 2 is a norm between the velocity at transition and field point considered. ft1 is given
by
ω2
ft1 = Ct1 gt exp −Ct2 t 2 [d2 + gt2 d2t ]
∆U
(2.35)
In the current research, the tripping terms ft1 and ft2 have been neglected (set to zero).
The values of closure coefficient is given in Table 2.2
2.2.3
k − ǫ Models
k − ǫ model is a two-equation turbulence model which solves two transport equations for
turbulent kinetic energy, k, and turbulent dissipation, ǫ. Unlike zero-equation and one-equation
14
Table 2.2
Closure coefficient for Spalart-Allmaras Model
σ
Cb1
Cb2
κ
Cw1
Cw2
Cw3
Cv1
Ct1
Ct2
Ct3
Ct4
2/3
0.1355
0.622
0.41
Cb1 /κ + ( 1 + Cb2 )/κ
0.3
2.0
7.1
1.0
2.0
1.1
2.0
models, the two-equation model is also called a “complete” turbulence model since it solves for
both turbulence length scale and velocity scale, instead of relating it to some flow properties.
The turbulence length scale is taken as k 3/2 /ǫ and velocity scale is approximated as k 1/2 . Using
these, the turbulent viscosity can be calculated as
µt = C µ ρ
k2
ǫ
(2.36)
where Cµ is a model constant in Standard k − ǫ model with value of 0.09.
Many variant of k−ǫ model are discussed in the literature. These formulations mainly differ
in the method of calculating turbulent velocity, calculation of source term (both generation
and destruction term) of ǫ-equation and model constants. In this research two variants of the
k − ǫ model have been used. The first is the original and the most popular model, developed
by Launder and Spalding [20] well known as the Standard k − ǫ model. The other model is
relatively new, proposed by Shih [21], known as Realizable k − ǫ model.
2.2.3.1
Standard k − ǫ model
The Standard k − ǫ model is a semi-empirical model in which the equation for turbulent
kinetic energy, k, is derived from the exact equation; whereas, the transport equation for ǫ is
more of a mathematical expression obtained from physical reasoning. This model is only valid
15
for fully turbulent flows.
Equations: Transport equation for turbulent kinetic energy is given by
∂
∂
∂
(ρk) +
(ρkui ) =
∂t
∂xi
∂xj
µt
µ+
σk
∂k
+ Pk − ρǫ
∂xj
(2.37)
where Pk is the production term for turbulent kinetic energy and ρǫ is the destruction term
due to turbulent dissipation.
The transport equation for turbulent dissipation is given by
∂
∂
∂
(ρǫ) +
(ρǫui ) =
∂t
∂xi
∂xj
µt
µ+
σǫ
ǫ
∂ǫ
ǫ2
+ C1ǫ Pk − C2ǫ ρ
∂xj
k
k
(2.38)
2
where C1ǫ kǫ Pk is the production and C2ǫ ρ ǫk is the destruction of turbulent dissipation.
The production term Pk is modeled as
Pk = −ρu′i u′j
∂uj
∂ui
= τij
∂xi
∂xj
Using Boussinesq’s hypothesis, Pk can be related to mean rate of strain tensor S =
as
Pk = µ t S 2
(2.39)
p
2Sij Sij
(2.40)
Calculating turbulent production term using Eq.2.40 may cause overproduction in turbulent
kinetic energy in highly accelerated or decelerated regions and the stagnation regions. To overcome this problem, Kato and Launder [22] proposed an ad-hoc modification to the production
term of turbulent kinetic energy equation. Kato and Launder proposed to replace one of the
strain rate, S, with vorticity, Ω, to modify the production term as
Pk = µt SΩ
(2.41)
s
∂uj 2
1 ∂ui
Ω=
−
2 ∂xj
∂xi
(2.42)
where Vorticity Ω is given as
16
The values of coefficient used are listed in Table 2.3.
Table 2.3
Closure coefficient for Standard k − ǫ Model
C1ǫ
C2ǫ
Cµ
σk
σǫ
2.2.3.2
1.44
1.92
0.09
1.0
1.3
Realizable k − ǫ model
The Standard k−ǫ model may sometimes violate two mathematical constraints on Reynolds
Stresses based on the physics of turbulent flows. Using Boussinesq hypothesis, normal Reynolds
stress for an incompressible flow can be written as
ρu′2 = −2µt
∂U
2
+ k
∂x
3
(2.43)
The normal Reynolds Stress ρu′2 , a positive quantity, may become negative or non-realizable
if the strain is high as
k ∂U
1
=
≈ 3.7
ǫ ∂x
3Cµ
(2.44)
Another defect of the Standard k − ǫ model (Shih (1995) [21]) is it may violate the Schwarz
2
′2
inequality for shear stresses u′α u′β ≤ u′2
α uβ at a high strain rate.
To overcome these deficiencies in the Standard k − ǫ model, two modifications have been
adopted in the Realizable k − ǫ model as follows:
1. A new transport equation for turbulent dissipation (new production and destruction
term) based on the dynamic equation of the mean-square vorticity fluctuation.
2. A new eddy viscosity formulation in which Cµ is a variable.
Equations: The transport equation for turbulent kinetic energy is given by
∂
∂
∂
(ρk) +
(ρkui ) =
∂t
∂xi
∂xj
µ+
µt
σk
∂k
+ Pk − ρǫ
∂xj
(2.45)
17
The equation for turbulent dissipation is given by
∂
∂
∂
(ρǫ) +
(ρǫuj ) =
∂t
∂xj
∂xj
µt
µ+
σǫ
ǫ2
∂ǫ
√
+ ρ C1 Sǫ − ρ C2
∂xj
k + νǫ
(2.46)
2
where ρ C1 Sǫ is the new production term and ρ C1 Sǫ − ρ C2 k+ǫ√νǫ is the new destruction term
for turbulent dissipation. The term C1 is given as
η
,
C1 = max 0.43,
η + 5.0
k
η=S ,
ǫ
S=
p
2Sij Sij
(2.47)
The transport equation for turbulent kinetic energy is the same as the Standard k−ǫ model.
The only difference is in the calculation of the source term of dissipation equation. One of the
remarkable features noted in the ǫ equation is it does not contain the Pk in production source
term. It is believed this helps to better represent the spectral energy transfer [23]. Also, the
singularity problem of destruction term in the ǫ equation of the Standard k − ǫ model (when
k = 0) does not exist because of the modified destruction term.
The turbulent viscosity is modeled the same as Eq. 2.36 but with the variable, Cµ . Cµ is
given by
Cµ =
q
U ≡ Sij Sij + Ω̃ij Ω̃ij ,
∗
1
∗
A0 + As kUǫ
Ω̃ij = Ωij − 2ǫijk ωk ,
(2.48)
Ωij = Ωij − ǫijk ωk
(2.49)
Ωij is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular
velocity ωk .
The constants A0 and As are given by:
A0 = 4.04, As =
√
6 cos φ
(2.50)
where
√
Sij Sjk Ski
1
φ = cos−1 ( 6W ), W =
,
3
S̃ 3
S̃ =
p
Sij Sij ,
1
Sij =
2
∂uj
∂ui
+
∂xi
∂xj
(2.51)
18
The values of the coefficient used in the Realizable k − ǫ Model are given in Table 2.4.
Table 2.4
Closure coefficient for Realizable k − ǫ Model
C1ǫ
C2ǫ
σk
σǫ
2.3
1.44
1.9
1.0
1.2
Rotor Source Modeling
In the present research, the method developed by Rajagopalan and Fanucci, Rajagopalan
and Lim, and Rajagopalan and Mathur [15, 16, 17], where the rotor has been treated as a
i
momentum source is adopted. The rotor is modeled as the momentum source term, Srotor
, and
is added to the momentum equation at the point where the blades pass. The magnitude of the
source term is a function of rotor geometry, blade cross-sectional aerodynamic characteristics,
and the local flow field.
i
i
~ , Ω, x, y, z, ρ, µ, Re, M, b, c)
Srotor
= Srotor
(Cl , Cd , α, α̇, V
(2.52)
The source term is time averaged and individual blades are not recognized. Explicit details of
the chordwise flow over rotor blade are not resolved.
To calculate rotor source term, the rotor is discretized into a spanwise element. The blade
properties, like chord length, airfoil thickness, plane deflection, and cross-sectional area are
assumed to be constant over the element. The source term is added to the grid cells, where the
blade pass while the rotor spins. A brief description of rotor source calculation is shown here.
More details can be found in Rajagopalan and Lim [16] and Rajagopalan and Mathur,[17].
2.3.1
Rotor source calculation
The rotor source term is calculated by computing the rotor force using the local angle of
attack and Mach number at each spanwise element. This is achieved by finding the Cl and
CD values from a lookup table for the given airfoil at that particular angle of attack and Mach
19
number.
~ = Vx êx + Vy êy + Vz êz is the velocity vector at position (x, y, z), where the
If V
Figure 2.1
Aerodynamic forces on blade section
rotor passes, then it can be transferred to the blade coordinate system (n, θ, s), such that the
velocity vector in rotor plane becomes
V~ ′ = vn ên + vθ êθ + Vs ês
(2.53)
If the rotor blade is flapping about the flapping hinge, then the tangential velocity of the
center of the blade element, with respect to the flapping hinge, is given by
~f lap = (β̇Rf )ên
V
(2.54)
where Rf is the radial distance from the hinge offset and β̇ is the flapping rate. In general, the
flapping motion of the blade is given by the flapping angle, β, represented as a Fourier series
of azimuthal angle ψ as
β = a 0 + ΣN
n=1 (an cos nψ + bn sin nψ)
(2.55)
From the above equation the flapping rate can be calculated as
β̇ =
∂β
ω
∂ψ
(2.56)
where ω is the angular velocity of rotor. Using the Eq. 2.53 and 2.54 the flow velocity relative
20
to the blade is given as
~rel = V
~ − (V
~bl + V
~f lap )
V
= vn′ ên + vθ′ êθ + Vs′ ês
(2.57)
~bl is the velocity of the blade in (n, θ, s) coordinate.
where V
The local angle of attack, α, at any spanwise position, s, is given by
α=φ−ǫ
(2.58)
ǫ = arctan(−vn′ /vθ′ )
(2.59)
where ǫ is inflow angle given by
The twist angle, φ, at that section can be determined from the rotor blade geometry twist
and from cyclic pitch variation of the blade root section θr as
θ r = A0 + Σ N
n=1 (An cos nψ + Bn sin nψ)
(2.60)
The local Mach number is calculated from
M ′ = v ′ /a
(2.61)
where a is the speed of sound and v ′ is the relative velocity seen by the airfoil section which
q
is v ′ = v ′ 2n + v ′ 2θ . After calculating the local Mach number and local angle of attack, the
aerodynamic coefficients, Cl and Cd , are obtained from the airfoil lookup table. Determining
the aerodynamic coefficients in this manner implicitly allows the compressibility effect on the
aerodynamic characteristics of the blade section.
21
Hence, the lift and drag in the directions, ên and êθ , can be calculated as
fn = L cos ǫ − D sin ǫ
fθ = L sin ǫ + D sin ǫ
(2.62)
where L and D are given by
L =
D =
1 ′2
ρv Cl c ds
2
1 ′2
ρv Cd c ds
2
(2.63)
The resultant force can be transformed to the computational domain to obtain the resultant
force acting on the blade F~ . The force acting on the fluid is −F~ which is added to the
momentum source term. However, in the steady (time-averaged) rotor source model, this force
is distributed to the computational cells passed by the rotor blade. The fraction of the force
added to the cells is calculated as follows.
If trev is the time period of rotation of rotor, i.e., trev = 2π/ω, then the time the blade
element spends in a control volume with width ∆θ is given by
t∆θ =
∆θ
ω
(2.64)
Then, the fraction of the time blade element spends in a cell is
tf rac =
∆θ
2π
(2.65)
If there are Nb blades, then the time averaged source term is given by
~ = Nb ∆θ (−F )
S
2π
(2.66)
22
CHAPTER 3.
DUST TRANSPORT MODELING
Dust particle entrainment during rotorcraft brownout is a process of soil particle movement
by the aerodynamic force of impinging rotor downwash with distinct phases of particle lift-off,
transport, and deposition. There are many factors which affect the physics of the process,
which includes strength of rotor downwash, ground signature, aerodynamic roughness length,
land-surface topography, soil texture, soil composition, etc. Although the phenomenon of
particle entrainment and transport takes place in most of the surface types, the intensity and
magnitude of entrainment is intense in desert arid region, where the percentage of dust in the
soil is high and surface vegetation does not restrict the motion of particles.
One of the closely related physical phenomenon to rotorcraft brownout can be found in
the geophysical process of aeolian sedimentology and particle transport. There have been a
number of studies completed in the geophysical research group, which deals with the behavior
of particle motion and transport by air as well as water. Although the basic physics of particle
motion, both in geophysical process and rotorcraft brownout, remains the same, the condition
in which this process takes place is different. In the geophysical level, flow is mainly parallel
to the ground; whereas, in the brownout scenario, the flow impinges on the ground with
many other complex aerodynamic processes. Development of mathematical formulations of
wind erosion problem under rotor flow is scarce in literature. The present research uses a
mathematical formulation based on empirical and semi-empirical relations from the geophysical
research group. Most of the mathematical relations discussed in this chapter are related to
wind erosion, Shao (2000) [24]. To model the transport of particles in the rotor flow field,
the particle motion is modeled as a multi-phase flow. The governing equation for particle
transport has been taken from the multi-phase research group, which formulates air as the
23
carrier medium and particle as the discrete medium ([25]).
3.1
Particle Characteristics
The particle transport process involves various complex interactions, like particle-flow interaction, particle-particle interaction, and particle-surface interaction. The behavior of these
interactions are influenced by various physical properties of particles, like shape, size, density,
terminal velocity, threshold velocity, etc.
3.1.1
Equivalent particle diameter
The particle size is one of the major factors determining particle transport behavior and
characterization. In general, particle size is closely related to particle shape. The shape
of a particle can be very irregular and varies from sphere to plates, very angular to wellrounded, and from rough to smooth in their texture. These properties can be defined by
various morphological parameters, like sphericity, roundness, and surface roughness. However,
it can be extremely difficult to precisely determine the size of the particle unless the shape is
well defined. A simple and efficient way to characterize the particle is through equivalent
particle size or diameter. Equivalent particle size of a particle is the diameter of the
equivalent sphere which retains certain aerodynamic or optical properties as identical to the
particle. Equivalent particle diameter can be based on mass, volume terminal velocity, drag,
etc. For example, mass equivalent and volume equivalent diameter are:
dm =
6m
πρp
1/3
dv =
6V
π
1/3
(3.1)
where m is particle mass, V is the particle volume, and ρp is the particle density. The choice
of the type of particle equivalent diameter depends upon the problem studied. For the current
research, equivalent diameter, based on drag and terminal velocity, is the most appropriate.
Based on the size, a particle can be characterized among roughly four categories as shown
in Figure 3.1, i.e., gravel (2000µm < d ≤ 2m), sand (60 < d ≤ 2000µm), slit (2 < d ≤ 60µm),
24
and clay (d ≤ 2µm) (based on MIT and British Standards Institute). For simplicity, in the
current work the particles have been divided in two categories - sand and dust. Sand particles
(d > 100µm) are those particles, which are relatively heavier and cannot be suspended in the
air for long time. They are deposited back in the ground, due to the gravitational effect, but
they play a major role in dust entrainment by the process known as saltation. Dust particles
(d < 100µ) are those, which are relatively smaller and play a major role in forming dust
clouds. Because of lower terminal velocity of dust particles, they are carried by the flow for
longer distance and time, particularly by the effect of turbulent eddies.
Figure 3.1
Particle size characterization
Generally soil contains distribution of particle sizes which is represented by particle size
distribution density function p(d), Figure 3.2. If P (d) is probability (or mass fraction) of
particle with diameter less than d, then
P (d) =
Z
d
p(d)δd
(3.2)
0
Given the particle size distribution function of a given soil, the mass fraction of any range of
size can be calculated. For example, mass fraction of a particle size between da and db is given
25
Figure 3.2
Particle size distribution in different regions
by
P (db ) − P (da ) =
3.1.2
Forces on particle
Various forces acting on a particle include
• Gravitational Force Fg
• Aerodynamics Drag Fd
• Aerodynamics Lift Fl
• Buoyancy Force Fb
• Magnus Force Fm
• Electrostatic Force Fe
Z
db
p(d)δd
da
(3.3)
26
• Particle-Particle Interaction Fp
Gravitational force is due to gravity, which is simply the weight of the particle acting vertically.
Fg = −mg
(3.4)
where g is the acceleration due to gravity.
Aerodynamic drag is produced when a particle moves relative to the surrounding fluid,
creating a pressure difference ∆p between frontal and wake regions, Figure 3.3. Additional
force is also exerted by shear σij of the fluid on the surface.
Figure 3.3
Aerodynamic drag and lift on spherical particle
Fdi = −
Z
pni dS +
S
Z
σij nj dS
(3.5)
S
The above equation for aerodynamic drag can be quantified by using drag coefficient CD
by
1
Fd = ρCD A |ur | Ur
2
(3.6)
where ur = u − up is the relative velocity of particle, A is the representative area, and up is
particle velocity.
27
Aerodynamic drag coefficient CD is the function of the particle Reynolds number, Rep =
ρur d/µ, Figure 3.4. There has been numerous studies completed on the relationship between
CD and Rep . In general, if Rep << 1, then it is assumed to be Stokes flow. For this region,
CD varies inversely with Rep . With increasing Rep , CD reaches a constant value around 0.45.
Known as the inertial region, this region is between 750 < Rep < 3.5 × 105 . At a much higher
Reynolds number, Rep > 3 × 105 , there is a sudden decrease in the CD value to around 0.1,
due to turbulence.
Figure 3.4
Aerodynamic drag coefficient, CD , as a function of Reynolds
number Rep
At Stokes flow region, the flow is a creeping flow, where the inertial effect can be neglected.
The governing equation in this region can be given by
∂ 2 ui
∂p
=µ
∂xi
∂xi ∂xj
(3.7)
By solving the above equation for flow over a spherical particle, force due to pressure can be
found as
Fp = πµdp U
(3.8)
28
where U is free stream velocity. Similarly, the shear force can be calculated as
Fτ = 2πµdp U
(3.9)
Thus, total drag force in Stokes region can be given as
Fd = 3πµdp U
(3.10)
Fd = 3πµdp ur
(3.11)
and in terms of relative velocity,
Using the Fd , the drag coefficient for Stokes flow can be found as
CD =
24
Rep
(3.12)
With an increase in Reynolds number, the inertial forces become more significant and
the drag coefficient is higher than Stokes drag. In the current work, a simpler expression of
CD (Rep ) is used
CD =
24
fdrag
Rep
(3.13)
where fdrag is known as the drag function and can calculated as [25]
fdrag =
1 + 0.15Re0.687
p
0.0183Re
if Rep ≤ 1000
(3.14)
if Rep > 1000
Similar to aerodynamic drag, aerodynamic lift is produced due to a difference in pressure
above and below the particle, due to difference in velocity (Bernoulli’s effect) as shown in
Figure 3.3. Aerodynamic lift can also be approximately calculated [24] through
1
Fl = Cl ρA(∇U 2 )d
2
(3.15)
29
where ∇U 2 is a gradient of U and Cl is the aerodynamic lift coefficient roughly proportional
to CD , as Cl = 0.85CD .
The buoyancy force, Fb , on the particle is proportional to the density ratio of the particle
to fluid. In the present case, the density of particle is of the order of 103 ; whereas, the density
of fluid is of the order of 1. Since the density ratio is of the order of 103 , Fb has been neglected.
Magnus force, Fm , is force caused by rotation of the particle and acts perpendicular to both
the direction of motion and rotation. Electrostatic force, Fe , is cased by electrostatic charge
on particle and electric field near the surface. In the present work, Fl , Fm and Fe have been
neglected.
Particle-particle Fp interaction can occur for two reasons. The first is due to particleparticle collision, which happens when the particle density is high. This also is neglected in
the current research. The other interaction is caused by a cohesive force, due to inter-molecular
attraction. This force is dominant when particle size is small.
3.1.3
Particle terminal velocity
Particle terminal velocity is the particle to fluid relative velocity, when the aerodynamic
forces balance the gravitational force on the particle [24]. At this instance, the particle experiences zero acceleration dupi /dt = 0. Neglecting all external forces, except the aerodynamic
drag force, Fd , and gravitation force, Fg , the simplified equation for particle motion can be
written as
m
dupi
= 3πµdp fdrag uri + δi3 mg
dt
(3.16)
Using the expression for fdrag from Eq. (3.13), the equation of particle motion can be written
as
fdrag
dupi
=
uir + δi3 g
dt
τv
(3.17)
where τv is the particle velocity response time,
τv =
ρp d2p
18µ
(3.18)
30
Applying dupi /dt = 0 in Eq. (3.17), horizontal components of terminal velocity are found
to be zero. The vertical component, wt , is
wt =
gτv
fdrag
(3.19)
Using Eq. (3.13) and Eq. (3.18), the vertical terminal velocity can be written as
wt =
4gdp σp
3CD (Rept )wt
(3.20)
where σp = ρp /ρ is the particle density ratio. Rept = wt dp /ν is particle Reynolds number
at terminal velocity. In general, the above equation is an implicit equation in wt , since CD
is the function of Rept , which, in turn, is a function of wt . This equation is solved using
numerical iteration. In the special case for flow around Stokes region, CD = 24/Rep and,
thus, wt = σp gd2p /18ν. At a very high Reynolds number, CD = 0.48 and wt = 1.66(σp gdp )(1/2) .
Plot showing the variation of terminal velocity with particle diameter at three different density
ratios are shown in Figure 3.5.
Figure 3.5
Particle terminal velocity, wt , Shao (2000) [24]
31
3.1.4
Particle threshold friction velocity
Particle threshold friction velocity, u∗t , is the minimum friction velocity, u∗ =
p
τw /ρ,
required to give the initial motion to the particle on the ground against the retarding gravitational force and cohesive forces.
Figure 3.6
The forces which drive the initial motion are mainly
Particle threshold friction velocity u∗t Iversen et al. (1976) [26]
the aerodynamic drag Fd and lift Fl , which are directly proportional to wall shear near the
ground. Particle threshold velocity depends upon capacity of the surface to resist wind erosion,
essentially a function of surface properties. In general, threshold velocity depends on various
factors like
• Particle diameter (dp )
• Particle size distribution (p(d))
• Particle density (ρp )
• Fluid density (ρ)
• Impact of impinging particle (air,water)
• Particle friction Reynolds number (Re∗ )
32
There have been many theories developed to calculate the threshold friction velocity. One of
the earliest theories, known as Bagnold Scheme [27] or Bagnold Formulae, was proposed
by Ralph Alger Bagnold (founder and first commander of the British Army’s Long Range
Desert Group during World War II).
u∗t = AB
p
σp gdp
(3.21)
where the coefficient AB = AB (Re∗t ) depends on particle threshold Reynolds number (Re∗t =
u∗i dp /ν). It was found empirically to be in the range of 0.1 and 0.2. A problem with this
scheme was it predicts a lower threshold friction velocity for particles below the viscous sub
layer. Later Greeley and Iversen [28] hypothesized that AB is not constant but is a function
of aerodynamic drag and the inter particle cohesive force.
AB = A1 F (Re∗t )G(d)
(3.22)
where G(d) is the function which takes cohesive force into account. The threshold friction
velocity for different particle types are plotted in Figure 3.6.
Very recently Lu and Shao [29] proposed that under ideal conditions, u∗t can be described
as a function of particle property alone
1
u∗t =
κ
s
a1
ρp
a2
gdp +
ρ
ρdp
(3.23)
where κ represents a surface roughness factor taken to be 1.0 in this current research. From
experiments the value of the coefficient a1 and a2 are found to be approximately 0.0123 and
3 × 10−4 kgs−2 , respectively. In the current research the above semi-empirical expression for
u∗t has been used.
33
3.2
Modes of Particle Transport
Based on different categories of particle and flow condition, the particle adopts a different
type of transport mechanism [24] (Figure 3.2).
Figure 3.7
Modes of particle transport
• Suspension
This type of particle transport is common among relatively smaller particles, referred
to as dust particles here. Once the particle has lift-off from the ground, it stays in the
air for a longer time and travels longer distance. Due to its low terminal velocity, it
is easily influenced by turbulent fluctuating velocities. Thus, for this type of particle,
turbulent diffusion plays a very important role in the transport mechanism. For most
geophysical processes, particles below . 70µm are considered to have a suspension type of
motion. Also, if the terminal velocity of the particle is very small compared to the vertical
turbulent velocity, then the particle follows the suspension mode. If the vertical terminal
velocity is approximated as κu∗ , then the particle goes to suspension, if wt /κu∗ << 1.
• Saltation
When fluid friction velocity is higher than the threshold friction velocity, the relatively
larger sand particle is entrained into the flow. Since the terminal velocity is large for
these particles, they fall to the ground, due to gravitational force. While in flight these
particles extract momentum from the flow and hit the ground with higher energy. This
34
impact further intensifies the emission of more particles and dust into the flow. Thus,
this mechanism plays a very important role in dust particle entrainment.
The layer in which saltating particles moves is called saltation layer. The thickness of
this layer, know as saltation layer thickness (zm ), is the maximum height the saltating
particle attains. The height of the saltation layer depends upon the friction velocity. As
predicted by Owen [30], saltation layer thickness can be given by
zm = c
u2∗
2g
(3.24)
where c is a constant assumed to be 1.6. Owen also hypothesized the saltation layer
behaves as a rough surface with roughness length, zos , proportioal to the saltation layer
thickness, zm , for the flow above the saltation layer.
• Creep
Larger particles > 1, 000µm are too heavy to be lifted by aerodynamic force. These
particles roll (or creep) along the ground. One of the major causes of surface creep is the
momentum imparted to the heavy particles from the impact of saltating particles.
3.2.1
Dust emission mechanism and modeling
So far in this chapter, various characteristics and properties of all sizes of particle have been
described. Since much anecdotal evidence suggests that in the case of brownout, relatively
small, fine, powder-like particles are the main ingredient of the brownout cloud [6], focus of
this research will be on dust particles. In general, there are three mechanism of dust emission
as depicted in Figure 3.8.
1. Aerodynamic Lift This occurs when aerodynamic lift is large enough to pick the dust
particles from the ground. Generally, these particles are very small, where the gravitational and aerodynamic lift are much smaller; whereas, the cohesive force can be
dominant. Dust emission, due to this mechanism, is small.
35
2. Saltation Bombardment As discussed before, when the saltating particle falls on the
ground, due to high momentum content, it transfers energy to the dust particles, which
helps overcome the binding forces. The dust emission, due to this mechanism, is an order
in magnitude higher than aerodynamic lift.
3. Disaggregation This happens when sand particles coated with dust disintegrate under
strong wind erosion.
Figure 3.8
3.2.2
Dust emission mechanism from Shao (2000) [24]
Vertical particle flux
Experiments conducted by Nalpanis [31] showed vertical particle flux, F , decreases exponentially with height.
F = Fo e
−λgz
u2
∗
(3.25)
where Fo is vertical particle flux at ground, z is the height from ground, and λ is a parameter
which around 50%. Through various field measurement with different kinds of soils accom-
36
plished by Gillette [32], it was found that Fo is proportional to un∗ , where n varies between
2.9 and 4.4 depending on soil type. Based on this relation, there has been a number of dust
emission schemes proposed by various researchers. Although a generalized scheme has not been
developed due to the large numbers of parameters involved, it has been found that dust emission is directly linked with intensity of saltation and indirectly related to the binding strength
of dust particles. For the current research, a scheme proposed by Marticorena and Bergametti
[33] has been used. The vertical particle flux depends upon saltation or horizontal particle flux
Q as
Fo = Qe13.4f −6.0
(3.26)
where f is percentage of clay in soil. For this research f = 0.1 was used.
Similar to vertical particle flux, saltation flux also depends on friction velocity u∗ . The
present research uses a model developed by Bagnold [27] to calculate the saltation flux Q.
Q=
where co = 0.25 +
0.33wt
u∗ .
0
if u∗ < u∗t
co ρ u3∗ 1 −
g
u2∗t
u2∗
(3.27)
if u∗ ≥ u∗t
As discussed earlier, particle entrainment happens when ground friction velocity, u∗ , is
greater than threshold friction velocity, u∗t . The friction velocity is calculated from ground
p
stress as u∗ =
(τw /ρ). In most geophysical research, like uniform steady parallel flow,
the ground stress is calculated from the gradient of mean flow velocity (the first moment of
velocity field). But, in the experiment conducted by Haehnel [10] on particle entrainment under
impinging jet, it was found that there was significant particle entrainment in the jet centerline
stagnation zone, where the wall stress based on gradient of mean flow velocity is zero. Instead
of wall stress based on first moment of velocity field, Haehnel found the driving force of the
particle entrainment is the Reynolds stress, v ′ w′ (second moment of velocity field), due to
turbulent velocity fluctuation. Although, Reynolds stress is same as average shear stress based
37
on the mean velocity field for steady parallel flows, it does not go to zero at the stagnation
region, as in the case of impinging jet. Haehnel also suggested the relationship of the Reynolds
stress with the turbulent kinetic energy as v ′ w′ = 0.2k. The current research uses Reynolds
stress, as suggested by Haehnel, to calculate the dust flux.
3.3
Dust Transport Modeling
A system or flow consisting of a mixture solid, liquid and gas phases is called a multiphase
system or flows. Multiphase flows can be classified into different categories, depending upon
the nature of the system or the flow regime. Based on the phases involved, multiphase flows
can be classified [34, 25] into various categories.
• Gas liquid flows
• Gas solid flows
• Liquid solid flows
• Three phase flows
In the case of rotorcraft brownout, flow can be considered as a solid-gas multiphase turbulent dispersed flow, where air is the carrier phase and dust is the discrete phase. Both phases
have separately defined volume fraction and velocity field. Conservation equations for the flow
of each phase are solved to get the density and velocity field. Transfer of mass, momentum,
and energy between the phases are taken care through coupling between the phases.
The two different theoretical models used to solve multiphase flows are Eulerian-based
multi-fluid model and Lagrangian-based particle tracking model. In both models the fluid
is solved as a continuum phase. In the Eulerian-based model, a cloud of dust is also solved
as continuum, where as in the Lagrangian-based model each dust particle is treated as an
independent entity and the trajectories are calculated individually.
The advantages and disadvantages of both of these models have been described by many
researchers (Durst et al (1984) [35], Shirolkara et al. (1996) [36], Crowe et al. (1996) [37]]. In
38
the Lagrangian-based model, since each particle is tracked separately, this model gives good
insight into particle dynamics. A disadvantage with this model is it becomes computationally
expensive when the number of particles is large. The Eulerian-based approach for particle
transport is easy to interpret and is computationally more efficient, since it solves an additional
transport equation for the discrete phase.
Mathematical modeling of multiphase flows depends on various properties of dispersed
phase flows. Some of the important properties which apply specifically to particle-gas flows
are described next.
3.3.1
Particle Stokes number
Stokes number is a dimensionless number which gives an idea of particle behavior in fluid
particle flows. The Stokes number of a particle can be defined as
Stv =
τv
τF
(3.28)
where τv is particle velocity response time given in Eq. 3.18 and τF is time characteristic of
flow field which can be given by
τF =
Lref
Uref
(3.29)
where Uref is the reference velocity of the flow field and Lref is characteristic reference length.
If Stv << 1, then the response time for particles to react to a change in flow behavior is
small. So particles will have enough time to respond to the change in flow and the particle
will have almost the same velocity as that of fluid (velocity equilibrium). If Stv >> 1, then
particle response time is high and relative velocity of the particle with respect to the flow will
be significant. An approximate relation between particle velocity, fluid velocity, and Stokes
number is given by
up
1
∼
u
1 + Stv
(3.30)
It can be seen as Stv → 0, particle velocity attains fluid velocity; whereas, when Stv → ∞,
the particle velocity approaches zero.
39
3.3.2
Particle volume fraction
Particle volume fraction ,αp , is defined as volume of the particle per unit volume.
αp =
lim
δV →δV o
δVp
δV
(3.31)
where δVp is the volume of the particle in the volume.
3.3.3
Particle bulk density (Dust cloud density)
Particle bulk density is the mass of particle per unit volume of mixture. In this current
work, it has also been reported as Dust Cloud Density, ρd .
ρd =
lim
δV →δV o
δMp
δV
(3.32)
where δMp is the mass of particle in volume. Dust cloud density is related to particle
volume fraction, αp , and particle density, ρp , as
ρd = αp ρp
3.3.4
(3.33)
Particle loading
The ratio of mass flux of particle to that of mass flux of fluid is called loading Z.
Z=
Ṁp
Ṁf
(3.34)
where Ṁp is the mass flow rate of particle and Ṁf is the mass flow rate of fluid. For small
Stokes number it can be approximated as Z = ρd /ρ.
3.3.5
Dilute versus dense flows and phase coupling
Dilute flow is a flow where the motion of a particle is governed by the forces acting upon
them by the fluid, whereas in a dense flow particle motion is governed by the collisions. In
40
general, if particle collision time, τc , is very large compared to particle response time, τv , then
the flow is called dilute flow, otherwise, it is a dense flow.
Figure 3.9
Dilute vs dense flows Crowe et al. (1998) [25]
Another way to differentiate dilute flow from dense flow [25] can be given in terms of
particle diameter as
dp <
3µ
Zρur
(3.35)
This relation is depicted in the Figure 3.9. For turbulent flows the relation can be given in
terms of standard deviation, σ, of particle fluctuation velocity as
dp <
1.33µ
Zρσ
(3.36)
The exchange of mass, momentum, and energy between phases is governed by phase coupling parameter. Since in the rotorcraft brownout case there is no phase change of particle,
mass coupling is not considered. Similarly, energy coupling is also neglected in this present
research. Three types of momentum coupling can be considered.
1. One-way coupling: Particle is influenced by fluid through aerodynamic drag; whereas
41
fluid is not being affected by the presence of particles.
2. Two-way coupling: Both particle and fluid affect each other through aerodynamic drag
and turbulence transfer.
3. Four-way coupling: In addition to two-way coupling, particle-particle collision is also
taken into consideration.
The type of coupling required for modeling muliphase flows varies from problem-to-problem.
In general, one-way coupling is straightforward and efficient. Requirement of two- or four-way
coupling is estimated through parameter magnitude estimation, which generally depends upon
factor-like particle loading Z, particle Stokes number Stv , particle bulk density ratio C = ρd /ρ,
etc.
One of the commonly used parameters to decide upon the coupling is known as momentum
coupling parameter, Πmom . An approximate definition of momentum coupling parameter is
given by Crowe et al. (1998) [25]
Πmom =
C
1 + Stv
(3.37)
If Πmom << 1, then coupling of momentum between phases is insignificant, and one way
coupling is used.
3.3.6
Rotorcraft Brownout Model
From the field studies achieved by Chatten [2] (Midwest Research Institute), on brownout
cloud characteristics on six rotorcraft configurations at the Yuma Proving Ground (YPG), it
was found that brownout dust cloud contains particles ranging from 2µm to larger than 350µm
(Figure 3.10). The mass concentration (mg/m3 ) at rotor tip location for various particle sizes
ranges from 10mg/m3 to 700mg/m3 .
Based on these test data and other anecdotal reports, it can be inferred that the major
contributor to brownout clouds is a relatively smaller particle (0 − 62µm), since number of
particle per unit volume (number density) is highest. High particle number density causes
42
Figure 3.10 Mass concentration of various particle size ranges at rotor tip
location during Sandblaster-2 field test on rotorcraft brownout
(from Midwest Research Institute (MRI) technical report
(Chatten (2007) [2]))
more scattering of light, thus causing visibility problems. Therefore, the focus of this current
research is on relatively small dust particles.
Particle Stokes Number in brownout: Considering particle density of ρp = 2800kg/m3
and air viscosity of µ = 1.7 × 105 , the particle response time varies from 9 × 10−4 s . τv .
3.5 × 10−2 s for particle range of 10µm to 62µm. Time characteristic of fluid τF can be found
in the range of 0.15s . τF . 0.6s, using Lref = 6 − 12m (rotor radius) and Uref = 20 − 40m/s
(wall jet velocity or rotor-induced velocity). Based on the time characteristics of particle and
fluid, the particle Stokes number varies between 1.5 × 10−3 . Stv . 0.23. This shows the
Stokes number is Stv < 1 for all dust particle sizes and Stv << 1 for smaller dust particles.
Particle loading in brownout: At a lower Stokes number, particle loading is approximated as Z = ρd /ρ. As seen from the test data, particle loading is of the order of Z ∼ 10−3 ,
the brouwnout flow condition can be considered as a dilute flow.
Phase coupling in brownout: With particle bulk density ratio C = ρd /ρ < 7−4 , momentum coupling parameter Πmom << 1. Accordingly, flow can be considered as one-way
coupled. However, the particle bulk density ratio has been found to be comparatively large
43
near the ground, which is better suited for two-way coupling. Irrespective of this, since this
zone (saltation and entrainment zone) of high bulk density ratio is generally much smaller
compared to the overall flow field, two-way coupling has been ignored for this region.
Considering the above properties of particle discrete phase during brownout and the fact
the number of particles involved is large, the Eulerian-based dust transport model has been
used to model and calculate the dust transport for the present case. With very low particle
volume fraction, volume fraction of fluid is almost equal to unity. With the one-way coupling
assumption, the fluid flow is not affected by the particle and the conservation equations (given
in CHAPTER 2) for turbulent rotor flow are solved independent of particle flow.
3.3.6.1
Dust Transport Equation
The mass transport equation of a particle can be calculated by summing the mass change
of the particle in the control volume to net efflux of mass from the control surface and the rate
of dust generation in this volume.
∂
∂
∂
(αp ρp ) +
(αp ρp upi ) −
(αp ρp wt ) = Sd
∂t
∂xi
∂x3
where up is particle velocity and Sd is dust source term. The third term,
(3.38)
∂
∂x3 (αp ρp wt ),
takes care of the settling effect on a particle, due to gravity. The dust source term, Sd , is equal
to dust generation near the ground due to dust emission. The magnitude of Sd is equal to
vertical dust flux given by Eq. 3.25.
The above equation can be written in terms of particle dust density, using the relation
ρd = αp ρp as
∂
∂
∂
(ρd ) +
(ρd upi ) −
(ρd wt ) = Sd
∂t
∂xi
∂x3
(3.39)
In turbulent flows, the particle velocity upi and dust density ρd can be divided into two
components, i.e., average component and fluctuating component (similar to derivation of RANS
equation) as
44
upi = ūpi + u′pi
(3.40)
ρd = ρ̄d + ρ′d
(3.41)
Replacing this in the transport equation and time averaging the equation [25], the equation
can be rewritten as
∂
∂
∂
∂
(ρ̄d ) +
(ρ̄d ūpi ) −
(ρ̄d wt ) = Sd −
(ρ′ u′ )
∂t
∂xi
∂x3
∂xi d pi
(3.42)
With reference to Fick’s Law, the dispersion of the dust in turbulent flows can be modeled
as
− ρ′d u′pi = Γd
∂ ρ̄d
∂xi
(3.43)
where Γd is dispersion or diffusion coefficient of dust in turbulent flows. The magnitude of Γd
depends upon the flow characteristics and can be determined through experiments. One of the
semi-empirical formulations used in geophysical research ([24]) is
Γd = ΓF
β 2 w2
1+4 2t
σ
(3.44)
where β is a dimensionless coefficient, ΓF is the fluid diffusion coefficient, and σ is the standard
deviation of turbulent velocity. In the current research, a relatively simple expression for Γd ,
mostly used in multiphase research [36], has been used.
Γd =
µpt
σtp
(3.45)
where µpt is turbulent particle viscosity and σtp is turbulent particle Schmidt number, which
varies from 0.34−0.7. Turbulent particle viscosity µpt is related to turbulent fluid eddy viscosity
µt by
µpt =
µt
1 + ττfv
where τf is fluid relaxation time given by τf = 0.2k/ǫ [36].
(3.46)
45
3.3.6.2
Dust Particle Velocity
The relative velocity between the dust particle and fluid is defined as
~ur = ~up − ~u
(3.47)
For dilute flows with low particle response time (0.001 < τv < 0.01s), an algebraic slip formulation [38, 23] is used to calculate particle relative velocity. Assuming that local equilibrium
has been reached between the particle and fluid, the relative velocity can be given by
~ur =
τv (ρp − ρ)
~a
fdrag
ρp
(3.48)
where ~a is acceleration given by the mixture velocity (fluid velocity for low particle volume
fraction) as
~a = f~ − (~u.∇)~u −
∂~u
∂t
(3.49)
where f~ is the acceleration due to body force. Since no electrostatic, magnetic effect, or other
body forces have been considered and the effect of gravitational force has already been taken
into account in the particle transport equation using the terminal velocity, wt , f = 0 has been
used to calculate particle relative velocity.
In turbulent flow, a diffusion term is added to the relative velocity to accommodate the
effect of turbulent fluctuations.
~ur =
τv (ρp − ρ)
~a + ~utr
fdrag
ρp
(3.50)
where
~utr = µpt
1
∇ρd
ρd
(3.51)
46
CHAPTER 4.
NUMERICAL METHOD AND ALGORITHM
The transport equation for all the flow variables are solved using Finite Volume Method
(FVM). Structured Cartesian grid has been used to discretize the domain.
4.1
General Convection-Diffusion Equation
The general Convection-Diffusion equation in three-dimensional Cartesian system can be written as:
∂
∂
(ρφ) +
∂t
∂x
∂φ
ρuφ − Γ
∂x
∂
+
∂y
∂φ
ρvφ − Γ
∂y
∂
+
∂z
∂φ
ρwφ − Γ
∂z
=S
(4.1)
where S is the external source term, and u, v and w denote the velocity components in the x,
y and z directions, respectively.
The equation can be simplified as:
∂(ρφ) ∂Jx ∂Jy
∂Jz
+
+
+
=S
∂t
∂x
∂y
∂z
(4.2)
where Jx , Jy and Jz represent the total (convection plus diffusion) fluxes defined by:
∂φ
∂x
∂φ
Jy ≡ ρvφ − Γ
∂y
∂φ
Jz ≡ ρwφ − Γ
∂z
Jx ≡ ρuφ − Γ
(4.3)
(4.4)
(4.5)
For a control volume given in Figure 4.1, the scalar variable φ is located at the control volume
47
center, P . The discretized form of this equation is obtained upon integration over the control
volume and is given by the following equation.
[(ρφ) − (ρφ)0 ]
∆x∆y∆z
+ α(Jφ−e − Jφ−w + Jφ−n − Jφ−s + Jφ−t − Jφ−b ) +
∆t
0
0
0
0
0
0
(1 − α)((Jφ−e
− Jφ−w
+ Jφ−n
− Jφ−s
+ Jφ−t
− Jφ−b
)=
α(SC + SP φ) + (1 − α)(SC0 + SP0 φ0 )∆x∆y∆z
(4.6)
The subscripts e, w, etc. refer to the values computed at the control volume faces, as seen
in Figure 4.1. The source term S has been linearized as S = SC + SP φ. For the unsteady
term, ρ and φ are assumed to prevail over the entire control volume. The ’old’ values from the
previous time step are denoted by subscript o.
Figure 4.1
The control volume for scalar variable φ.
For incompressible flows, the continuity equation reduces to:
(Fe − Fw + Fn − Fs + Ft − Fb ) = 0
(4.7)
(Fe0 − Fw0 + Fn0 − Fs0 + Ft0 − Fb0 ) = 0
(4.8)
Multiplying the continuity equation (Eq. 4.8) by αφ for the current values and (1 − α)φ0 for
the ‘old’ values and subtracting them from the discretized Convection-Diffusion Equation, we
48
get:
(φ − φ0 )ρ0
∆x∆y∆z
+ α[(Jφ−e − Fe φ) − (Jφ−w − Fw φ) +
∆t
(Jφ−n − Fn φ) − (Jφ−s − Fs φ) +
(Jφ−t − Ft φ) − (Jφ−b − Fb φ)] +
0
0
(1 − α)[(Jφ−e
− Fe0 φ0 ) − (Jφ−w
− Fw0 φ0 ) +
0
0
(Jφ−n
− Fn0 φ0 ) − (Jφ−s
− Fs0 φ0 ) +
0
0
(Jφ−t
− Ft0 φ0 ) − (Jφ−b
− Fb0 φ0 )]
= α(SC + SP φ) + (1 − α)(SC0 + SP0 φ0 )∆x∆y∆z
(4.9)
The terms, such as (Jφ−e − Fe φ), can be rewritten as follows:
Jφ−e − Fe φ = aE (φ − φE )
(4.10)
Jφ−w − Fw φ = aW (φW − φ)
(4.11)
Jφ−n − Fn φ = aN (φ − φN )
(4.12)
Jφ−s − Fs φ = aS (φS − φ)
(4.13)
Jφ−t − Ft φ = aT (φ − φT )
(4.14)
Jφ−b − Fb φ = aB (φB − φ)
(4.15)
where the coefficients are obtained by the following equations:
aE
= De A(|Pe |) + [[−Fe , 0]]
(4.16)
aW
= Dw A(|Pw |) + [[Fw , 0]]
(4.17)
aN
= Dn A(|Pn |) + [[−Fn , 0]]
(4.18)
aS = Ds A(|Ps |) + [[Fs , 0]]
(4.19)
aT
(4.20)
= Dt A(|Pt |) + [[−Ft , 0]]
aB = Db A(|Pb |) + [[Fb , 0]]
(4.21)
49
The subscripts E, W , etc. denote the neighboring control volume center values.
The diffusion terms are computed from:
De =
Dw =
Dn =
Ds =
Dt =
Db =
Γe ∆y∆z
(δx)e
Γw ∆y∆z
(δx)w
Γn ∆x∆z
(δy)n
Γs ∆x∆z
(δy)s
Γt ∆y∆x
(δz)t
Γb ∆y∆x
(δz)b
(4.22)
(4.23)
(4.24)
(4.25)
(4.26)
(4.27)
and the Peclet numbers from:
Pe =
Pw =
Pn =
Ps =
Pt =
Pb =
Fe
De
Fw
Dw
Fn
Dn
Fs
Ds
Ft
Dt
Fb
Db
(4.28)
(4.29)
(4.30)
(4.31)
(4.32)
(4.33)
For the Power Law scheme, the function A(|P |) is computed from:
A(|P |) = [[0, (1 − 0.1|P |)5 ]]
(4.34)
50
and the equation becomes
φ
ρ∆x∆y∆z
+ α(aE + aW + aN + aS + aT + aB − SP ∆x∆y∆z)φ
∆t
= α(aE φE + aW φW + aN φN + aS φS + aT φT + aB φB ) +
(1 − α)[a0E φ0E + a0W φ0W + a0N φ0N + a0S φ0S + a0T φ0T + a0B φ0B ] −
(1 − α)[a0E + a0W + a0N + a0S + a0T + a0B − SP0 ∆x∆y∆z]φ0 +
φ0
ρ∆x∆y∆z
+ ∆x∆y∆z(αSC + (1 − α)SC0 )
∆t
(4.35)
This equation can be rewritten as:
ρ
∆x∆y∆z
φ + α aP φ = α(aE φE + aW φW + aN φN + aS φS + aT φT + aB φB ) + btotal (4.36)
∆t
where:
aP
= aE + aW + aN + aS + aT + aB − SP ∆x∆y∆z
btotal = α b + bo
(4.37)
(4.38)
b = SC ∆x∆y∆z
(4.39)
bo = (1 − α)[a0E φ0E + a0W φ0W + a0N φ0N + a0S φ0S + a0T φ0T + a0B φ0B ]
−(1 − α)a0P φ0 + ρ
4.1.1
∆x∆y∆z 0
φ + (1 − α)b0
∆t
(4.40)
Scalar variables
The scalar variables φ, diffusion coefficient Γ, and source term (S) for each scalar transport
equation is given in Table 4.1.
It should be noted that for dust transport equation particle density ρp and particle velocities
upi is used in place of fluid density ρ and fluid velocities ui respectively.
Source term linearization is completed in the following manner.
51
Table 4.1
Coefficient and source term for scalar transport equations
φ
Spalart-Allmaras
ν̃
Γ
ρ
σ (ν
S
ρCb1 [1 − ft2 ]S̃ ν̃ + σρ Cb2 |∇ν|2
+ ν̃)
h
+ρ Cw1 fw −
Standard k − ǫ
Realizable k − ǫ
Dust Transport
i
ν̃ 2
d
+ ρft1 ∆U 2
µt
σk
Pk − ρǫ
µ+
µt
σǫ
C1ǫ kǫ Pk − C2ǫ ρ ǫk
µ+
µt
σk
Pk − ρǫ
µt
σǫ
ρ C1 Sǫ − ρ C2 k+ǫ√νǫ
k
ǫ
k
ǫ
αp
Cb1
f
κ2 t2
µ+
µ+
2
2
µtp
σpt
F
Spalart-Allmaras
For ν̃ equation
SC
=
ρ
Cb2 |∇ν|2 + ρft1 ∆U 2
σ
Cb1
+ν̃ max 0, ρCb1 [1 − ft2 ]S̃ + ρ Cw1 fw − 2 ft2
κ
SP
Cb1
= min 0, ρCb1 [1 − ft2 ]S̃ + ρ Cw1 fw − 2 ft2
κ
ν̃
d2
ν̃
d2
(4.41)
52
Standard k − ǫ
For k equation
SC
= Pk
SP
= −ρ
cµ ρk
µt
(4.42)
For ǫ equation
SC
SP
ǫ
= C1ǫ Pk
k
ǫ
= −C2ǫ ρ
k
(4.43)
Realizable k − ǫ
For k equation
SC
= Pk
SP
= −ρ
cµ ρk
µt
(4.44)
For ǫ equation
SC
= ρ C1 Sǫ
SP
= −ρ C2
k+
ǫ
√
νǫ
(4.45)
Dust Transport
For αP equation
SC
= F
SP
= 0
(4.46)
53
4.2
Momentum Equations
The discretized momentum equations for the velocities can be derived similarly from the general
Convection-Diffusion Equation with the φ term replaced by the velocity components u, v and
w. In this research, a staggered mesh is used. This means the velocity components u, v and w
are positioned at different control volume faces as seen in Figure 4.2. Therefore, the equations
must be adjusted accordingly.
(a) u-momentum control volume
Figure 4.2
(b) v-momentum control volume
Control volumes of u and v momentum
If the source term (btotal ) of u, v, and w momentum equations in the discretized momentum
are denoted as btotal
, btotal
, and btotal
u
v
w , then the discretized form of the momentum equation for
u becomes:
ρ
∆x∆y∆z
uP + αauP uP = α(aE uE + aW uW + aN uN + aS uS + aT uT + aB uB ) + btotal
(4.47)
u
∆t
54
where
auP
= aE + aW + aN + aS + aT + aB − SPu ∆x∆y∆z
= αbu + buo + (pW − pE )∆y∆z
btotal
u
bu = SCu ∆x∆y∆z
(4.48)
(4.49)
(4.50)
buo = (1 − α)[a0E u0E + a0W u0W + a0N u0N + a0S u0S + a0T u0T + a0B u0B ]
−(1 − α)(auP uP )0 + ρ
∆x∆y∆z 0
u + (1 − α)b0u
∆t
(4.51)
Similar expressions can be written for v- and w- momentum equations as:
ρ
∆x∆y∆z
(4.52)
vP + αavP vP = α(aE vE + aW vW + aN vN + aS vS + aT vT + aB vB ) + btotal
v
∆t
where
avP
= aE + aW + aN + aS + aT + aB − SPv ∆x∆y∆z
btotal
= αbv + bvo + (pS − pN )∆x∆z
v
bv = SCv ∆x∆y∆z
(4.53)
(4.54)
(4.55)
0
0
0
0
bvo = (1 − α)[a0E vE
+ a0W vW
+ a0N vN
+ a0S vS0 + a0T vT0 + a0B vB
]
−(1 − α)(avP vP )0 + ρ
∆x∆y∆z 0
v + (1 − α)b0v
∆t
(4.56)
and
ρ
∆x∆y∆z
total
wP + αaw
(4.57)
P wP = α(aE wE + aW wW + aN wN + aS wS + aT wT + aB wB ) + bw
∆t
55
where
aw
P
= aE + aW + aN + aS + aT + aB − SPw ∆x∆y∆z
= αbw + bw
btotal
o + (pB − pT )∆x∆y
w
bw = SCw ∆x∆y∆z
(4.58)
(4.59)
(4.60)
0
0
0
0
]
+ a0N wN
+ a0S wS0 + a0T wT0 + a0B wB
+ a0W wW
bw
= (1 − α)[a0E wE
o
0
−(1 − α)(aw
P wP ) + ρ
4.2.1
∆x∆y∆z 0
w + (1 − α)b0w
∆t
(4.61)
Pressure correction equation
The discretized momentum can be solved similar to discretized scalar equations, if the
pressure field is given. Since the correct pressure field is not known initially, a guessed pressure
field is used to solve for the velocity field. Since the pressure filed is not correct, the calculated
velocity field may not satisfy the continuity equation. Thus, a correction (pressure) is applied
to correct the velocity filed such that the mass conservation is satisfied in each computational
cell. The pressure correction is calculated by solving the pressure correction equation.
For the pressure control volume Pp , the discretized momentum equation of east face u
velocity can be written similar to Eq. 4.51 as
ae ue = Σanb unb + b + (pP − pE )Ae
(4.62)
where subscript nb denotes neighboring point and b is the source term without the pressure
term. If the guessed pressure filed is denoted by p∗ , then the imperfect velocity is denoted as
u∗ , v ∗ , and w∗ . These imperfect velocity fields are solution of the following discretized equation
ae u∗e = Σanb u∗nb + b + (p∗P − p∗E )Ae
(4.63)
∗
an vn∗ = Σanb vnb
+ b + (p∗P − p∗N )An
(4.64)
∗
at we∗ = Σanb wnb
+ b + (p∗P − p∗T )At
(4.65)
56
Solving the velocities using the discretized momentum equation with guessed pressure field
p∗ in general will not satisfy the continuity equation. If p′ is the pressure correction required
to get the correct pressure field then
p = p∗ + p′
(4.66)
The velocity is corrected as
u = u∗ + u′
v = v∗ + v′
w = w∗ + w′
(4.67)
Subtracting Eq. 4.65 from Eq. 4.62, and using Eq. 4.67, we have
ae u′e = Σanb u′nb + (p′P − p′E )Ae
(4.68)
ae u′e = (p′P − p′E )Ae
(4.69)
u′e = de (p′p − p′E )
(4.70)
by dropping the Σanb u′nb term
where
de ≡
Ae
.
ae
(4.71)
The Eq. 4.70 is called the velocity-correction formula. The corrected velocities are given as
ue = u∗e + de (p′p − p′E )
vn = vn∗ + dn (p′p − p′N )
wt = wt∗ + dt (p′p − p′T )
(4.72)
The discretized equation of p′ is obtained from continuity equation. The continuity equation
57
is
∂ρ ∂(ρu) ∂(ρv) ∂(ρw)
+
+
+
=0
∂t
∂x
∂y
∂z
(4.73)
Integrating the equation over the pressure control volume yields
(ρP − ρ0P )∆x∆y∆z
+ [(ρu)e − (ρu)w ] ∆y∆z + [(ρv)n − (ρv)s ] ∆z∆x
∆t
+ [(ρw)t − (ρw)b ] ∆x∆y = 0
(4.74)
Substituting the velocity component by the velocity correction expression (Eq. 4.72) gives
aP p′P = aE p′E + aW p′W + aN p′N + aS p′S + aT p′T + aB p′B + b
(4.75)
where
aE
= ρe de ∆y∆z
aW
= ρw dw ∆y∆z
aN
= ρn dn ∆z∆x
aS = ρs ds ∆z∆x
aT
= ρt dt ∆x∆y
aB = ρb db ∆x∆y
aP
= aE + aW + aN + aS + aT + aB
b =
(ρP − ρ0P )∆x∆y∆z
+ [(ρu∗ )e − (ρu∗ )w ] ∆y∆z
∆t
+ [(ρv ∗ )n − (ρv ∗ )s ] ∆z∆x + [(ρw∗ )t − (ρw∗ )b ] ∆x∆y
4.2.2
(4.76)
The pressure equation
Since, in the velocity correction equations, the neighbor-point velocity correction is not
considered, the pressure correction has the entire responsibility to correct the velocity field.
This leads to severe pressure-correction fields, which is efficient in correcting the velocity field,
58
but pressure field may not be calculated properly. To overcome this problem, the pressure is
calculated through pressure equation and pressure-correction is only used for correcting the
velocities.
The momentum equation can be re-written as
ue =
Σanb unb + b
+ de (pP − pE )
ae
(4.77)
Σanb unb + b
ae
(4.78)
Defining pseudo-velocity ûe by
ûe =
The momentum equation can be written as
ue = ûe + de (pP − pE )
(4.79)
Similarly other velocities can be written as
vn = v̂n + dn (pP − pN )
(4.80)
wt = ŵt + dt (pP − pT )
(4.81)
Replacing the above definition of velocity in discretized continuity equation (Eq. 4.74) we have
a P p P = a E p E + a W p W + a N pN + a S p S + a T p T + a B p B + b
(4.82)
where aE , aN , aN , aS , aT , and aB are similar to pressure correction coeffiecient as given in Eqn.
4.76 and b is given by
b =
(ρP − ρ0P )∆x∆y∆z
+ [(ρû)e − (ρû)w ] ∆y∆z
∆t
+ [(ρv̂)n − (ρv̂)s ] ∆z∆x + [(ρŵ)t − (ρŵ)b ] ∆x∆y
(4.83)
59
4.2.3
The SIMPLER algorithm
The SIMPLER algorithm [11] consist of solving the pressure equation to obtain the pressure
field and uses the pressure-correction filed only to correct the velocity field. The algorithm
follows :
1. Start with a given initial flow field.
2. Calculate the unsteady portion of the center coefficients.
3. Calculate the unsteady portion of the source terms.
4. Calculate the coefficients of the momentum equations and the pseudo-velocities.
5. Using the pseudo-velocities, calculate the source term for the pressure equation.
6. Calculate the coefficients for the pressure equation and solve the pressure equation to
obtain the pressure field.
7. Using the calculated pressure field, solve the momentum equations to obtain the velocity
field (u, v, and w).
8. Calculate the source terms of the pressure correction equation and solve for the pressure
corrections.
9. Correct the velocities using the velocity correction equations.
10. Return to Step 4 and repeat until convergence.
11. Go to Step 2 and start with a new time level.
4.3
Boundary conditions
To solve the discretized equations, boundary conditions are needed for each transport
equation. Mainly three boundary conditions, inlet, outlet, and wall boundary conditions, have
been used in almost all the simulations carried out in this current research.
60
4.3.1
Inlet boundary conditions
Momentum Equation: Inlet velocity is provided.
Spalart-Allmaras Model: Turbulent viscosity ratio, νt /ν, is provide at the inlet. In
general, the turbulent viscosity ratio varies from 3.0 to 10.0.
k − ǫ Model: Turbulent intensity (I) and turbulent length scale (l) are provided, based
on the turbulent kinetic energy and turbulent dissipation calculated.
3
kinlet = (U I)2
2
ǫinlet
3
k2
= Cµ
l
3
4
(4.84)
(4.85)
Turbulent intensity of 1 to 10% has been used. Turbulent length scale of the order rotor
radius(l ≈ R), has been used.
Dust Transport Model: Dust concentration at inlet is zero αp = 0.
4.3.2
Outlet boundary conditions
Momentum Equation: Outlet velocity is calculated, based on mass conservation, such
that
ṁout = ṁin
(4.86)
where ṁout is total mass outflow from all outlet boundaries and ṁin is the inflow from inlet
boundaries.
Spalart-Allmaras Model: Neumann boundary condition is used such that
∂ ν̃
=0
∂n
(4.87)
61
k − ǫ Model: Neumann boundary condition is used such that
∂k
∂n
∂ǫ
∂n
= 0
= 0
(4.88)
Dust Transport Model:Neumann boundary condition is used such that
∂αp
=0
∂n
4.3.3
(4.89)
Wall boundary condition
Momentum Equation: No-slip boundary condition is used such that velocity on the
face of the wall is zero. In the case of k − ǫ model, µe is calculated on the wall faces using
the function approach (discussed in APPENDIX-1). For the cells which lie inside the body
geometry, the diffusion is set to infinite.
Spalart-Allmaras Model:
ν̃ = 0
(4.90)
k − ǫ Model: k − ǫ model cannot be integrated to the wall. To bridge this, wall function
(discussed in APPENDIX-1) is used to provide the boundary condition.
Dust Transport Model: Dust concentration at the wall has been taken to be zero.
αp = 0
4.3.4
(4.91)
Computational process
During the computational process, RANS equation is solved first at each time step. While
solving the RANS, the rotor source model is invoked at each iteration to add the rotor source
terms. From the converged flowfield, the turbulent properties are solved by solving turbulence
models. After the calculation of flow field variables, the dust tranport model is started, which
is solved for dust concentration. The entire algorithm is given in the flow chart shown in Figure
62
4.3.
Figure 4.3
Flowchart of Rotorcraft Brownout Model
63
CHAPTER 5.
RESULTS
A three-dimensional Cartesian grid-based incompressible laminar solver has been developed
from an existing two-dimensional solver. This three-dimensional solver has been used as a
testbed for implementing, validating, and testing various turbulence models and the brownout
dust transport model. For solving the rotor flow, the existing rotor source model has been
integrated with the turbulent solver.
5.1
5.1.1
Validation
3D incompressible solver
To validate the three-dimensional laminar incompressible solver, a standard test case of lid
driven cavity has been used. In this test case, the cavity is cube shaped with 1m length(l) and
the lid moves at constant speed (Ulid ) of 1m/s. The test has been completed at two Reynolds
numbers of 100 and 400, where the Reynolds number is defined as
Re =
ρUlid l
µ
(5.1)
The wall of the cavity has been modeled as viscous no slip boundary condition. The domain
is discretized into a 42 × 42 × 42 non-uniform grid.
Results of center plane velocity distributions are compared with Flux Corrector Method
(FCM) used by Wirogo [39]. At Re = 100, horizontal velocity, u, and vertical velocity, v, at
center plane are plotted in Figure 5.1.
For Re = 400, horizontal velocity (u) and vertical
velocity (v) at center plane are plotted in Figure 5.2. The contour plot of velocity magnitude
at center plane for both Reynolds numbers is plotted in Figure 5.3
64
Figure 5.1
u and v velocity at center plane in driven cavity at Re = 100
Figure 5.2
u and v velocity at center plane in driven cavity at Re = 400
At both the Reynolds numbers, u velocity profiles match well with FCM. In the v velocity
profile, the present scheme (Power-Law) predicts a higher maximum velocity, due to the high
numerical diffusivity caused by the first order scheme.
5.1.2
Turbulence models
Four turbulence models have been integrated with the 2D and 3D incompressible solvers.
The first test is achieved by simulating turbulent flow over the 2D flat plate. Skin friction
2 ) is plotted against Reynolds number (Re = ρV x/µ) for all
coefficient (Cf = 2τw /ρVinf
inf
turbulence models in Figure 5.4. Results are compared with empirical formulation for skin
65
Figure 5.3
Velocity magnitude contour plot at center plane of driven cavity
at Re = 100 and 400
friction over the square plate for turbulent flows by Schlichting [40]
Cf =
0.059
(1/5)
Rex
(5.2)
The Baldwin-Lomax and Spalart-Allmaras models give good comparisons for most of the
Reynolds number range. Both k − ǫ models over-predict skin friction at a relatively low
Reynolds number (< 3.0 × 105 ). This is because both k − ǫ models used are for high Reynolds
number fully turbulent flows. Among both k − ǫ models, Realizable k − ǫ gives a better results.
Another case of 3D circular impinging jet is used for comparing the turbulence models.
Circular impinging jet closely resembles rotorcraft in ground effect, which have various features
like recirculation, stagnation zone, and wall jet.
The simulation set-up is shown in Figure 5.5. Set-up and flow condition are similar to
experiments carried out by Cooper et al. [41]. In this present case, the H/D ratio is taken
to be 2 and Reynolds number is 70, 000, where Reynolds number is based on impinging bulk
velocity and diameter of nozzle.
Velocity profiles at three different radial locations (r/R = 1.0, 2.0 and 3.0) are plotted in
Figure 5.6 for all turbulent models and have been compared with the experimental results.
Velocity magnitude contour plots for all the cases are shown in Figure 5.7.
66
Figure 5.4
Skin friction coefficient with Reynolds number on flat plate
It can be seen that at r/R = 1.0, Laminar and the Baldwin-Lomax model give the best
result. Further downstream in the wall jet region, both laminar and Baldwin-Lomax model
over-predict the maximum wall jet velocity and at higher heights under-predict the profile
velocity. In the wall jet region, one-equation model and two-equation model give better comparisons.
Turbulent shear stress u′ w′ is also plotted at various radial locations in Figure 5.8. It can
be seen that none of the linear eddy viscosity-based model predicts the turbulent shear stress
accurately. The reason, as given by Craft et al. [42], is due to anisotropic behavior of turbulent
fluctuating velocity in the impinging jet case, where turbulent kinetic energy is dominated by
normal straining and fluctuating velocity normal to the wall is larger than that parallel to the
wall at the impingement zone. Nevertheless, after comparing all the results, the two-equation
model seems to give the best qualitative results.
67
Figure 5.5
Impinging jet simulation set up
(a) r/R = 1.0
(b) r/R = 2.0
(c) r/R = 3.0
Figure 5.6
Impinging jet velocity profile
68
(a) Baldwin-Lomax
(b) Spalart-Allmaras
(c) Standard k − ǫ
(d) Realizable k − ǫ
(e) Laminar
Figure 5.7
Velocity magnitude contour plot (impinging jet)
69
(a) r/R = 1.0
(b) r/R = 2.0
(c) r/R = 3.0
Figure 5.8
Reynolds Shear Stress u′ w′ in impinging jet
70
5.1.3
Rotor model validation
To validate rotor source modeling, simulation has been accomplished on a rotor similar
to experiments completed by Rabbott [43]. A two-bladed rotor has been used with 15f t
diameter. Simulations were carried out at various collective pitches (Θ). For the current
simulation, NACA0012 airfoil was used. Chord length was 0.153 times the radius and there
was zero twist along the blade. The variation of aerodynamic blade load with radial locations
at various collective pitches is plotted in Figure 5.9.
Figure 5.9
(a) Θ = 3.0o
(b) Θ = 4.5o
(c) Θ = 8.5o
(d) Θ = 9.2o
Aerodynamic blade load comparisons with Rabbott’s experiment at various Θ
The rotor source model with tip correction, is able to perfectly capture the trend of the
71
aerodynamic load along the blade radius.
5.1.4
Brownout model validation
It is very difficult to validate the brownout model. There have not been many quantitative
experimental data on brownout. The present model has been qualitatively validated with wind
tunnel experiments conducted by Nathan and Green [44] at University of Glasgow. In their
experiment, a small model rotor was used and was placed at one radius above the ground.
The wind tunnel speed was set so the thrust-normalized advanced ratio µ∗ is equal to 0.65.
Thrust-normalized advanced ratio was given by
µ
µ∗ = p
(CT /2)
(5.3)
where µ is rotor advance ratio and CT is coefficient of thrust.
µ = V cos α/(ΩR)
(5.4)
where α is tilt angle, Ω is angular velocity, and R is the radius of rotor. CT is given by
CT =
T
ρ(ΩR)2 A
(5.5)
where T is rotor thrust and A is rotor disc area. In their experiment, very fine, powder-like
particles were used as dust. The motion of the particle transported by rotor flow was captured
by high speed camera (Figure 5.10).
Simulation was achieved on a similar rotor and flow
condition using the current brownout model. Dust particle density distribution and velocity
are shown in Figure 5.11. It can be seen that the numerical model was able to capture the
same extent of the brownout cloud. Another remarkable feature captured by the present model
is the stagnation region of the dust particle on the ground at about one radius upstream.
72
Figure 5.10 Particle transport by rotor downwash in wind tunnel experiment (Nathan et al. (2008) [44])
(a) Dust particle Density
(b) Particle velocity vector
Figure 5.11 Dust particle distribution predicted by brownout model
73
5.2
Turbulence Model Comparison for Rotor in Ground Effect
For the selection of an appropriate turbulence model for further studies, simulation was
carried out with all turbulence models to solve for turbulent rotor flow. The focus of this study
was to verify the ability of various turbulence model to predict rotor outwash and wall jet of a
hovering rotor in ground effect. This study used the CH-53E rotor, hovering at a height of two
rotor radius at thrust coefficient CT = 0.007. The field test data are documented by Preston
[45].
Comparisons of various turbulence models (Laminar, Baldwin-Lomax, Spalart Allmaras,
Standard k−ǫ, and Realizable k−ǫ) with flight test data and Vorticity Transport Model (VTM)
used by Phillips and Brown [46] for ground velocity profiles at different radial locations are
shown in Figure 5.12.
As seen in these plots, all models under-predict the maximum wall jet velocity particularly
close to the rotor, at radial position r/R = 1.0, 1.25 and 1.5. The reason might be due to
the presence of the fuselage in flight test data not modeled in the current simulations. In the
wall jet region at r/R = 1.75 and above, all models closely predict maximum wall jet velocity.
Although at farther radial locations, models predict lower velocity at higher height. This might
be due to numerical diffusion caused by first order schemes and the turbulence models used.
Near the wall jet region, Standard k − ǫ and Realizable k − ǫ have the best comparisons.
Based on the these outcomes and results from the circular impinging jet case, the Realizable
k − ǫ model was found to be most suitable and was used for all subsequent simulations.
74
(a) r/R = 1.0
(b) r/R = 1.25
(c) r/R = 1.50
(d) r/R = 1.75
(e) r/R = 2.0
(f) r/R = 3.0
Figure 5.12 Ground velocity profile for hovering CH-53E rotor
75
5.3
Fuselage Effect on Wake and Ground Signature
To study the effects of body (fuselage) on wake and ground signature of the rotor, single
rotor helicopter was used for simulation at different heights in ground effect (IGE). Two sets of
simulations were completed, first without the fuselage (named rotor configuration) and second
with fuselage (named rotor-fuselage configuration). Height of the rotor center was varied from
H = 0.5R to 8R. The properties of the rotor are given in Table 5.1.
Table 5.1
Rotor properties for helicopter in hover
Rotor radius(R)
Number of blades
Twist
Root cutout
Hinge offset
Chord
12m
7
−17o (linear)
0.325
0.0633
0.0617R
The blade collective pitch was maintained so that the coefficient of thrust was 0.007 when
hovering at H = 1.0R.
Simulation was in full turbulence mode, with Realizable k − ǫ model. The mesh size of
103 and 84 in x and y direction, respectively, was used, with a refined grid around the body.
Approximately 50 cells were used along the rotor diameter were used to properly capture the
rotor flow. The mesh size in the z direction varied from 75 to 94 for rotors at different heights
with a refined grid around the body and rotor. The refined grid with stretching was used near
the ground to capture the rotor outwash properly. Simulation time for each configuration was
between 15 to 25s, depending upon the time required by the wall jet to reach the steady state.
p
The ground signature was represented as the distribution of friction velocity (u∗ = τw /ρ)
along the ground. This gave an idea of erosion capability of the rotor downwash and outwash.
The erosion and dust entrainment took place in the region where the ground friction velocity
was higher than the particle threshold friction velocity. Ground signature of rotor wake and
oil flow pattern, with and without fuselage, at different heights are shown in Figures 5.13 to
5.21. The black color circle shows the position of the rotor.
As seen, effects of fuselage are significant at most heights. At lower height of H = 0.5 and
76
0.75R, the strength and extent of ground signature outside the rotor perimeter is not much
different for rotor and rotor-fuselage configuration, except near the tail of fuselage. There
is a fuselage footprint of low friction velocity zone just below the rotor. The effect is more
pronounced between height of H = 1.0 and 3R. A zone of high ground friction velocity around
the perimeter of rotor was observed for rotor-fuselage configuration between these heights.
This is due to an increase in downwash velocity in the presence of a body in the rotor wake.
The extent of the ground signature is almost the same between both configurations except it
is a little distorted due to the tail of the fuselage in the rotor-fuselage configuration.
An interesting observation was found in the pattern of oil flow distribution for rotor-fuselage
configuration. With an increase in rotor height, the shadow of the fuselage (tadpole-like shape
made by the oil flow pattern) can be seen rotating counter-clock wise. This is due to the
swirling flow nature of the rotor wake.
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.13 Ground friction velocity at rotor height H = 0.5R
77
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.14 Ground friction velocity at rotor height H = 0.75R
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.15 Ground friction velocity at rotor height H = 1.0R
78
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.16 Ground friction velocity at rotor height H = 1.5R
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.17 Ground friction velocity at rotor height H = 2.0R
79
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.18 Ground friction velocity at rotor height H = 2.5R
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.19 Ground friction velocity at rotor height H = 3.0R
80
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.20 Ground friction velocity at rotor height H = 4.0R
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.21 Ground friction velocity at rotor height H = 8.0R
81
The contour plot of velocity magnitude for rotor and rotor-fuselage configuration at rotor
height of H = 0.5 to 2R is shown in Figures 5.22 to 5.25. This plot is shown in the longitudinal
plane (along the fuselage) and four cross planes (slices) as show in the figures. A front view of
all planes (slice A to D) is shown in four thumbnails below each figure.
It is observed that at height 0.5R, the rotor wake and outwash are not much affected by
the presence of the fuselage. The same observation can be found in the ground signature
discussed earlier. At heights H = 1.0 to 2.0R, the effect of fuselage can be clearly seen in all
slices. In addition to the unsteady behavior below the fuselage, two major differences are the
increase in the bulk velocity of rotor wake and increase in the wall jet velocity. Presence of
fuselage decreases the flow area, thus, the fuselage increases the downwash velocity. Increase
in downwash velocity increases the wall jet velocity.
The ground velocity profile at three different lateral positions from rotor center (y = 1,
2, and 3R) for different rotor heights are shown in Figure 5.26. The comparison has been
made between rotor configurations and rotor-fuselage configurations. As seen in the figure, the
velocity profile for both configurations are not much different at lower heights (H = 0.5 and
0.75R). As observed earlier, fuselage effect is most significant between H = 1.0 to 3R, after
which the effect starts to decrease. The maximum wall jet velocity at three radial locations for
different heights is plotted in Figure 5.27. At the lateral position y = 1.0R, the wall jet velocity
first decreases with height for heights below 3R and then increases for higher heights for both
configurations. At all the heights, the wall jet velocity for the rotor-fuselage configuration is
always higher at this position. At the lateral position y = 2.0R, the wall jet velocity decreases
and becomes asymptotic to a constant value. At lower heights, H < 2.5R, wall jet velocity
for the rotor-fuselage configuration is higher than the rotor configuration. There is a crossover
at H = 2.5R. The wall jet velocity for the rotor is higher. This same trend is observed at
y = 3.0R, but the crossover takes place at H = 2R. The difference in wall jet velocity between
both configurations follows the same trend at all the locations, which reaches a maximum
around H = 1.5R.
82
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.22 Velocity magnitude contours at H = 0.5R
83
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.23 Velocity magnitude contours at H = 1.0R
84
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.24 Velocity magnitude contours at H = 1.5R
85
(a) Rotor only
(b) Rotor and Fuselage
Figure 5.25 Velocity magnitude contours at H = 2.0R
86
(a) H = 0.5R
(b) H = 0.75R
(c) H = 1.0R
(d) H = 1.5R
(e) H = 2.0R
(f) H = 2.5R
Figure 5.26 Ground velocity profile at different radial location for different
rotor heights
87
(a) H = 3.0R
(b) H = 4.0R
(c) H = 8.0R
Figure 5.26 (continued) Ground velocity profile at different radial location for different rotor heights
88
(d) y = 1.0R
(e) y = 2.0R
(f) y = 8.0R
Figure 5.27 Wall jet velocity at different radial locations
89
5.4
Fuselage Effect on Brownout
To study the fueslage effect of the brownout, the particle transport model has been used
for two configurations, rotor-only configuration and rotor-fuselage configuration, at different
heights. Rotor properties are the same as in the last section with N ACA0012 airfoil. For
all subsequent simulations, particle size of dp = 10µm is used with a particle density of ρd =
2800kg/m3 . The roughness length on the ground is taken as zo = 250µm.
The brownout cloud after t = 5s (≈ 15 rotor revolution) is shown in Figure 5.28 for both
configurations. Dust clouds shown in the images are the iso-surface of dust cloud density (ρd )
of 50mg/m3 . At a height of 0.5R, substantial differences in dust cloud patterns are observed.
A fountain of dust can be seen coming out around the region of rotor center in rotor-only
configuration. In rotor-fuselage configurations, the rise in dust is blocked by the presence of
the body. To visualize the dust characteristics near the rotorcraft, the dust cloud in the rear
of the rotor center lateral plane (y-z plane) is removed, as shown in Figure 5.29. The various
layers (iso surface) of dust density can be seen through the removed side as seen in Figure 5.30,
where the x-axis is normal to the images. These images also give an idea of the extent of dust
cloud height forming ahead of the rotor. As the height of the rotor increases, the height of the
dust fountain at the rotor’s center decreases for rotor configurations, due to the effect of rotor
downwash. In contrast, for the rotor-fuselage configurations, the height of the dust fountain
just below the fuselage increases and attempts to catch up with the fuselage. This is caused
by the recirculating flows in the wake of the fuselage, which lift the dust higher, as seen in the
particle path lines plots in Figure 5.31. This effect of the recirculating flows and the height of
the fountain below the fuselage decreases after a height of 2.0R. At all these heights, the dust
density in the dust fountain is much higher in the rotor-fuselage simulations. Even in the rotor
outwash region, higher density clouds (larger iso-surface value) are found in the rotor-fuselage
simulations.
Although the core of dust clouds in the rotor-fuselage configuration are determined denser,
the average height achieved by the dust cloud is almost the same for both configurations at
each height.
90
Rotor Only
Fuselage and Rotor
(a) H = 0.5R
(b) H = 0.5R
(c) H = 1.0R
(d) H = 1.0R
(e) H = 1.5R
(f) H = 1.5R
Figure 5.28 Dust cloud formation after 5s for rotor-only and rotor-fuselage
configuration
91
Rotor Only
Fuselage and Rotor
(a) H = 2.0R
(b) H = 2.0R
(c) H = 2.5R
(d) H = 2.5R
Figure 5.28 (continued) Dust cloud formation after 5s for rotor-only and rotor-fuselage configuration
92
Figure 5.29 Visualizing dust cloud characteristics
93
Rotor Only
Fuselage and Rotor
(a) H = 0.5R
(b) H = 0.5R
(c) H = 1.0R
(d) H = 1.0R
(e) H = 1.5R
(f) H = 1.5R
(g) H = 2.0R
(h) H = 2.0R
(i) H = 2.5R
(j) H = 2.5R
Figure 5.30 Dust cloud characteristics seen through rotor center cut plane
for rotor-only and rotor-fuselage configurations
94
Rotor Only
Fuselage and Rotor
(a) H = 0.5R
(b) H = 0.5R
(c) H = 1.0R
(d) H = 1.0R
(e) H = 1.5R
(f) H = 1.5R
(g) H = 2.0R
(h) H = 2.0R
(i) H = 2.5R
(j) H = 2.5R
Figure 5.31 Particle path lines at rotor mid plane for rotor-only and rotor-fuselage configurations
95
5.5
Fuselage Shape Effect on Wake and Ground Signature
To study the effect of fuselage shape on wake and ground signature, four simplified geometries with circular, elliptical, square, and rectangular cross-sections have been use for the
fuselage. The width of each fuselage is 9m. The height of elliptical and rectangular fuselage is
7m. The rotor properties are given in Table 5.2. Collective pitch has been set so the coefficient
of thrust at H = 1.0R is 0.007 for all the fuselage configurations. This is done by running
the simulation for each rotor configuration at different collective pitches to obtain a relation of
thrust coefficient and collective pitch, at a particular height. Using this relation, the collective
pitch is calculated to get the desired coefficient of thrust.
Figure 5.32 Four fuselage shapes used
Table 5.2
Rotor properties for helicopter with different fuselage shape
Rotor radius(R)
Number of blades
Twist
Root cutout
Hinge offset
Chord
Airfoil N ACA0012
8.02m
4
−12o (linear)
0.2
0.03
0.0617R
The ground signature (ground friction velocity) of each fuselage shape at various heights is
96
plotted in Figures 5.33 to 5.37. As seen, no significant difference has been found in the ground
signature for different fuselage shapes, if the projection of the fuselage on the ground is same.
Contour plot of velocity magnitude and streamline on lateral and longitudinal rotor mid
plane are shown in Figures 5.38 to 5.40 for all fuselage shapes at heights of 0.6, 1.0, and
1.5R. At height 0.6R, flow is mostly attached below the fuselage for all configurations, except
there is a recirculation just below the fuselage nose for the square, rectangular, and elliptical
configurations. At heights of 1.0R, the region below the fuselage becomes highly unsteady for
all configurations. Recirculating flows can be seen in the front view of all configurations. Flow
is highly un-symmetric in all cases, due to unsteadiness and the swirling nature of the rotor’s
downwash. The recirculating region below the fuselage nose is also found in rectangular and
square fuselages. Similar trends are observed for all higher heights (H = 1.5R and above).
(a) Circular Fuselage
(b) Square Fuselage
(c) Rectangular Fuselage
(d) Elliptical Fuselage
Figure 5.33 Ground friction velocity comparison at H = 0.6R
97
(a) Circular Fuselage
(b) Square Fuselage
(c) Rectangular Fuselage
(d) Elliptical Fuselage
Figure 5.34 Ground friction velocity comparison at H = 1.0R
(a) Circular Fuselage
(b) Square Fuselage
(c) Rectangular Fuselage
(d) Elliptical Fuselage
Figure 5.35 Ground friction velocity of different fuselage cross-section at
H = 1.5R
98
(a) Circular Fuselage
(b) Square Fuselage
(c) Rectangular Fuselage
(d) Elliptical Fuselage
Figure 5.36 Ground friction velocity comparison at H = 2.0R
(a) Circular Fuselage
(b) Square Fuselage
(c) Rectangular Fuselage
(d) Elliptical Fuselage
Figure 5.37 Ground friction velocity comparison at H = 3.0R
99
(a) Circular Fuselage (Front View)
(b) Circular Fuselage (Side View)
(c) Square crosss-section(Front View)
(d) Square crosss-section (Side View)
(e) Rectangular crosss-section(Front View)
(f) Rectangular crosss-section (Side View)
(g) Elliptical crosss-section(Front View)
(h) Elliptical crosss-section (Side View)
Figure 5.38 Velocity magnitude and streamlines for different fuselage
shapes at H = 0.6R
100
(a) Circular Fuselage (Front View)
(b) Circular Fuselage (Side View)
(c) Square crosss-section (Front View)
(d) Square crosss-section (Side View)
(e) Rectangular crosss-section (Front View)
(f) Rectangular crosss-section (Side View)
(g) Elliptical crosss-section (Front View)
(h) Elliptical crosss-section (Side View)
Figure 5.39 Velocity magnitude and streamlines for different fuselage
shapes at H = 1.0R
101
(a) Circular Fuselage (Front View)
(b) Circular Fuselage (Side View)
(c) Square crosss-section (Front View)
(d) Square crosss-section (Side View)
(e) Rectangular crosss-section (Front View)
(f) Rectangular crosss-section (Side View)
(g) Elliptical crosss-section (Front View)
(h) Elliptical crosss-section (Side View)
Figure 5.40 Velocity magnitude and streamlines for different fuselage
shapes at H = 1.5R
102
5.6
Multi-Rotor Configuration Effect on Brownout
Full rotorcraft configurations with single rotor, tandem rotor, tilt rotor, and quad rotor
have been used for this study (Figure 5.41). The rotor properties used are given in Table
5.3. N ACA0012 airfoil has been used for the rotor blades. The collective pitch of each rotor
is adjusted so the CT /σ for all the rotorcrafts are the same, while hovering at OGE (out of
ground effect). The mechanism of calculating the collective pitch is the same as discussed in
last section. σ is the solidity of the rotor given by σ = Nb c/πR, where Nb is the number of
rotor blades and c is the rotor chord.
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.41 Rotor configurations
The simulation was carried out for four rotor heights 1.0, 1.5, 2.0, and 2.5Rref , where Rref
is the rotor radius of a single rotor helicopter. Corresponding non-dimensional heights (H/R)
for the tandem rotor are 1.33, 2, 2.66, and 3.66. For tilt and quad rotor, the non-dimensional
103
Table 5.3
Single rotor
Tandem rotor
Tilt rotor
Quad rotor
Properties of rotor for multi-rotor configurations
Rotor Radius
Rotor Blades (per rotor)
Twist
σ
CT
12m
9m
6m
6m
5
3
3
3
12o
12o
16o
16o
0.095
0.083
0.105
0.105
0.0071
0.0063
0.0076
0.0077
heights (H/R) are 2, 3, 4, and 5.
Velocity streamlines at three different lateral planes – front rotor mid plane, fuselage mid
plane, and rear rotor mid plane – are plotted for tandem and quad rotors at various heights
in Figures 5.42 and 5.46. The tilt rotor streamlines are plotted on lateral rotor mid plane in
Figure 5.44. Streamlines along longitudinal plane are also plotted in Figures 5.43, 5.45, and
5.47.
In tandem rotor, large recirculating regions are observed underneath the fuselage at rotor
heights of H = 1.0, 1.5, and 2.0Rref . At fuselage mid plane, rotor wake are found to be shifted
on the starboard side. This is due to the effect of mixing two counter rotating and swirling
rotor downwashes. The rotor outwash and the wall-jet thickness in the longitudinal plane are
determined to be similar to the single rotor configuration. The recirculating region near the
ground at height 2.5Rref has diminished.
Similar to tandem rotor, tilt and quad rotor configurations show a heavy recirculating
zone below the fuselage at a height of 1.0, 1.5, and 2.0Rref . In particular, the region below
the fuselage wings is found highly unsteady. One of the major differences between tilt and
quad rotor with respect to tandem rotor is the rotor outwash and wall jet thickness along the
longitudinal plane. At lower altitudes, the outwash and the wall jet of tilt and quad rotor are
found thicker compared to tandem and single rotor. This plays a very important role in dust
particle transport in front of the fuselage.
104
(a) H = 1.0Rref
(b) H = 1.5Rref
Figure 5.42 Streamlines in lateral plane for tandem rotor hovering at different heights
105
(a) H = 2.0Rref
(b) H = 2.5Rref
Figure 5.42 (continued) Streamlines in lateral plane for tandem rotor hovering at different heights
106
(c) H = 1.0Rref
(d) H = 1.5Rref
(e) H = 2.0Rref
(f) H = 2.5Rref
Figure 5.43 Streamlines in longitudinal plane for tandem rotor hovering at
different heights
107
(a) H = 1.0Rref
(b) H = 1.5Rref
(c) H = 2.0Rref
(d) H = 2.5Rref
Figure 5.44 Streamlines in lateral plane for tilt rotor hovering at different
heights
108
(a) H = 1.0Rref
(b) H = 1.5Rref
(c) H = 2.0Rref
(d) H = 2.5Rref
Figure 5.45 Streamlines in longitudinal plane for tilt rotor hovering at different heights
109
(a) H = 1.0Rref
(b) H = 1.5Rref
Figure 5.46 Streamlines in lateral planes for quad rotor hovering at different heights
110
(a) H = 2.0Rref
(b) H = 2.5Rref
Figure 5.46 (Continued) Streamlines in lateral planes for quad rotor hovering at different heights
111
(c) H = 1.0Rref
(d) H = 1.5Rref
(e) H = 2.0Rref
(f) H = 2.5Rref
Figure 5.47 Streamlines on longitudinal plane for quad rotor hovering at
different heights
112
To compare the dust erosion capability of different rotor configurations, the ground signature of each configuration at different heights is plotted in Figures 5.48 to 5.51. A significant
difference of ground friction velocity has been found between single and tandem rotor with
respect to tilt and quad rotors at the same dimensional height. In particular, tandem rotor is
determined to give strong ground signature for all heights. Even at the same non-dimensional
height (H/R) of 2.0, (comparing Figures 5.50(a), 5.49(b), 5.48(c), and 5.48(d)), the tandem
rotor has the most intense ground friction.
Dust clouds (iso-surface of 50mg/m3 ) after simulation time of 5s for different rotor configurations at different heights are shown in Figures 5.52 to 5.55. At all heights the spread of
the dust cloud is more for the single and tandem rotors. Dust characteristics along the rotorcraft’s lateral plane is visualized by removing the dust cloud on the rear side of the fuselage
mid lateral plane. Dust characteristic seen through fuselage’s center lateral plane (front view)
are shown in Figures 5.56 to 5.59, where the x-axis is normal to the images. Similarly, for
dust characteristics along the longitudinal plane (side view), the dust cloud is removed from
the port side of the fuselage’s longitudinal plane (x-z plane) and seen through this plane in
Figures 5.60 to 5.63, where the y-axis is normal to the plane. Although in all scenarios single
and tandem rotor form more denser clouds, in terms of front visual reference, these two are
determined as better. As seen in the front view at heights 1.0 and 1.5Rref , huge dust clouds
form in front of the tilt and quad rotor, thus restricting the pilot’s front visibility. This is
mainly caused by transport of dust particles by stronger and thicker rotor outwash in front of
the fuselage in tilt and quad rotor configuration. In terms of side view perspective, tilt and
quad rotor are determined as better. There is an unusually huge dust cloud on the starboard
side of the tandem rotor at greater heights. This is caused by dust carried by the shifted rotor
downwash on the starboard side as discussed earlier.
Another interesting phenomenon observed is the dust fountain just below the fuselage tilt
and quad rotor is not as intense as compared to single and tandem rotors. This might be due
to a smaller rotor in tilt and quad rotors, thus, smaller downwash area hitting the ground as
shown in the ground signature.
113
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.48 Ground friction velocity comparison at H = 1.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.49 Ground friction velocity comparison at H = 1.5Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.50 Ground friction velocity comparison at H = 2.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.51 Ground friction velocity comparison at H = 2.5Rref
114
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.52 Dust cloud comparison for multiple rotor configurations at
H = 1.0Rref
115
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.53 Dust cloud comparison for multiple rotor configurations at
H = 1.5Rref
116
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.54 Dust cloud comparison for multiple rotor configurations at
H = 2.0Rref
117
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.55 Dust cloud comparison for multiple rotor configurations at
H = 2.5Rref
118
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.56 Dust cloud characteristics through fuselage lateral mid plane
at H = 1.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.57 Dust cloud characteristics through fuselage lateral mid plane
at H = 1.5Rref
119
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.58 Dust cloud characteristics through fuselage lateral mid plane
at H = 2.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.59 Dust cloud characteristics through fuselage lateral mid plane
at H = 2.5Rref
120
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.60 Dust cloud characteristics through fuselage longitudinal mid
plane at H = 1.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.61 Dust cloud characteristics through fuselage longitudinal mid
plane at H = 1.5Rref
121
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.62 Dust cloud characteristics through fuselage longitudinal mid
plane at H = 2.0Rref
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.63 Dust cloud characteristics through fuselage longitudinal mid
plane at H = 2.5Rref
122
Time evolution of dust for all configurations at height 1.0Rref is shown in Figures 5.64
- 5.68. After 2 s, the dust cloud height formed by tandem rotor is significantly higher than
for other configurations. It can be clearly seen that the major dust entrainment is from the
starboard side of the tandem rotors. After 3 s, all configurations show considerable size of dust
cloud formation. Single rotor shows dust formation all along the perimeter of the rotor. Some
effects of the fuselage tail are also found in the rear zone. The tandem rotor also shows dust
formation from all sides, but dust cloud heights are larger on the lateral sides. Dust cloud size
in the tilt rotor are comparatively small with a major evolution of dust coming from front and
rear. The quad rotor also shows similar trends with significant dust evolution from the sides
as well. After 5 s, all the configurations are almost engulfed in the dust cloud and there is no
visual reference for pilots. The single rotor still has some near rotor visual reference. After
10 s, all rotorcraft are engulfed in the cloud, without any visual reference. The dust cloud is
almost axisymmetric for a single rotor. In the tandem rotor, a huge dust cloud is seen on the
starboard side; whereas, for tilt and quad rotors, it is seen in the front and rear.
Comparison of fuselage effect on multi rotor configurations is shown in Figures 5.69 - 5.71.
Plots are shown for ground friction and dust clouds for multi-rotor configurations, with and
without fuselage, hovering at 1.0Rref . Similar to the helicopter, there is a significant effect of
fuselage on ground signature. The strength of ground friction velocity is found higher in the
case with fuselage for all configurations. Average height of the dust for rotor and rotor-fuselage
configurations are found to be almost the same, but the dust cloud density for the fuselage
rotor configurations are higher.
123
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.64 Dust cloud evolution after 2s at H = 1.0Rref
124
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.65 Dust cloud evolution after 3s at H = 1.0Rref
125
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.66 Dust cloud evolution after 5s at H = 1.0Rref
126
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.67 Dust cloud evolution after 7s at H = 1.0Rref
127
(a) Single rotor
(b) Tandem rotor
(c) Tilt rotor
(d) Quad rotor
Figure 5.68 Dust cloud evolution after 10s at H = 1.0Rref
128
(a) Ground friction velocity for tandem (b) Ground friction velocity for tandem
rotors only
rotors with fuselage
(c) Tandem rotors only (dust cloud)
(d) Tandem rotors with fuselage (dust cloud)
(e) Tandem rotors only (Pilot’s front view)
(f) Tandem rotors with fuselage (Pilot’s front view)
(g) Tandem rotors only (Pilot’s side view)
(h) Tandem rotors with fuselage (Pilot’s side view)
Figure 5.69 Comparison of tandem rotors with and without fuselage at
H = 1.0Rref
129
(a) Ground friction velocity for tandem (b) Ground friction velocity for tandem
rotors only
rotors with fuselage
(c) Tilt rotors only (dust cloud)
(d) Tilt rotors with fuselage (dust cloud)
(e) Tilt rotors only (Pilot’s front view)
(f) Tilt rotors with fuselage (Pilot’s front view)
(g) Tilt rotors only (Pilot’s side view)
(h) Tilt rotors with fuselage (Pilot’s side view)
Figure 5.70 Comparison of tilt rotors with and without fuselage at
H = 1.0Rref
130
(a) Ground friction velocity for tandem (b) Ground friction velocity for tandem
rotors only
rotors with fuselage
(c) Quad rotors only (dust cloud)
(d) Quad rotors with fuselage (dust cloud)
(e) Quad rotors only (Pilot’s front view)
(f) Quad rotors with fuselage (Pilot’s front view)
(g) Quad rotors only (Pilot’s side view)
(h) Quad rotors with fuselage (Pilot’s side view)
Figure 5.71 Comparison of quad rotors with and without fuselage at
H = 1.0Rref
131
CHAPTER 6.
CONCLUSION, REMARKS, AND FUTURE WORK
In the current research, a computational test bed has been developed for studying rotorcraft
brownout. Various turbulence models have been implemented and tested for impinging type
flows, from which the Realizable k − ǫ model has been selected as more suitable for rotor flows.
Momentum source, based time averaged rotor model, has been integrated to solve for rotor
flow. The dust transport model based on Eulerian frame of reference has been used for solving
the brownout dust transport. The dust entrainment on the ground has been modeled through
empirical relations used in geophysical research.
Based on the simulation completed on rotor and fuselage-rotor configurations it has been
determined that presence of body has a significant effect in brownout conditions. This effect is
particularly pronounced between heights of 1.0 to 3.0 times the rotor radius. The same effect
is also observed in multiple rotor configurations.
In the simulations with different fuselage shapes, it has been found that the cross-sectional
shape of the fuselage does not have much effect on ground signature if the shadow (the projection on the ground plane) of the body is the same. Although there is not much difference in the
wake, separation zones have been observed near the fuselage nose for square and rectangular
shaped fuselage because of sharp corners.
In the study of multiple rotor configuration, it has been found the tandem rotor has the
most significant ground signature. Although the amount of dust entrainment in tandem is
high, from the pilot’s front view perspective, has the least visual restriction. Both tilt and
quad rotor show lower visual performance because of large dust cloud evolution in front of the
fuselage.
132
Future work
Although the ability of the present computational model to predict and
calculate the brownout scenario looks promising, further work is needed to model an unsteady
rotor, which can capture the effect of an individual rotor.
All the test cases studied in the current research are completed for hovering rotors at
different heights; whereas in an actual scenario, the height of the rotorcraft decreases with
time. Also, current tactics used by the pilot while landing in an arid desert region is by flying
at very low altitudes parallel to the ground and suddenly land like a fixed wing aircraft. By
achieving this, the pilot prevents the rotorcraft from engulfing itself in the dust cloud and,
thus, extend the visual reference for a longer time. To model this actual scenario, a moving
body method must be developed.
133
BIBLIOGRAPHY
[1] Rogers. S. J. (1968). Evaluation of The Dust Cloud Generated by Helicopter Rotor Downwash. Technical Report 67-81, U.S. Army Aviation Material Laboratories.
[2] Chatten, C. (2007). Sandblaster 2 Support of See-Through Technologies for Particulate Brownout. Task 5 Technical Report,U.S. Army Aviation and Missile Command
Under Contract No. W31P4Q-07-C-0215 Midwest Research Institute (MRI),Project No.
110565.1.001
[3] Keller, J. D., Whitehouse, G. R., Wachspress, D. A. ,Teske, M. E., and Quackenbush,
T. R. (2006). A Physics-Based Model of Rotorcraft Brownout for Flight Simulation. 62nd
Annual American Helicopter Society Forum, Phoenix, AZ.
[4] Ryerson, C. C., Hachnel, R. B., Koenig, G. G., and Moulton, M. A. (2005).Visibility
Enhancement in Rotorwash Clouds. AIAA Aerospace Sciences Meeting and Exhibit, Reno,
Nevada.
[5] Haehnel, R. B., Moulton, M. A., Wenren, Y., and Steinhoff, J. (2008).A Model to Simulate
Rotorcraft-Induced Brownout. 64th Annual Forum of the American Helicopter Society,
Montreal.
[6] Phillips, C. and Brown, R. E. (2008). Eulerian simulation of the fluid dynamics of helicopter brownout. American Helicopter Society 64th Annual Forum, Montreal, Canada,
April 29 - May 1, 2008.
[7] Phillips, C., Kim, H. W. and Brown, R. E. (2009). The Effect of Rotor Design on the
Fluid Dynamics of Helicopter Brownout. 35th European Rotorcraft Forum, Germany.
134
[8] Wachspress, D. A., Whitehouse, G. R., Keller, J. D., Yu, K., Gilmore, P., Dorsett, M.,
and McClure, K. (2009). A high fidelity brownout model for real time flight simulation
and trainer. American Helicopter Society 65th Annual Forum, Grapevine.
[9] Keller, J. D., Whitehouse, G. R., Wachspress, D. A., Milton, E. T., Quackenbush, T. R.
(2006). A Physics-Based Model of Rotorcraft Brownout for Flight Simulation Applications.
American Helicopter Society 62nd Annual Forum, Phoenix, AZ.
[10] Haehnel, R. and Dade, W. B.. (2008). Physics of Particle Entrainment Under the Influence
of an Impinging Jet. Proceedings Army Science Conference, Orlando, FL.
[11] Patankar, S. V. (1980). Numerical heat transfer and fluid flow. Washington, DC, Hemisphere Publishing Corp..
[12] Wilcox, D.C., (1998) Turbulence Modeling for CFD. 2nd Ed., DCW Industries, Inc..
[13] Launder, B. E., Reece, G. J. and Rodi, W. (1975). Progress in the Development of a
Reynolds-Stress Turbulent Closure. Journal of Fluid Mechanics, 68 (3), 537–566.
[14] Craft, T. J., Launder, B. E., and Suga, K. (1997). Prediction of turbulent transitional
phenomena with a nonlinear eddy-viscosity model. International Journal of Heat and
Fluid Flow, 18 (1), 15–28.
[15] Rajagopalan, R. G. and Fanucci, J. B. (1985). Finite difference model for vertical axis
wind turbines. Journal of Propulsion and Power, 1 432–436.
[16] Rajagopalan, R. G. and Lim, C. K. (1991). Laminar Flow Analysis of a Rotor in Hover.
Journal of American Helicopter Society, 36 (1) 12–23.
[17] Rajagopalan, R. G. and Mathur, S, R. (1993). Three Dimensional Analysis of a Rotor in
Forward Flight. Journal of American Helicopter Society, 38 14–25.
135
[18] Baldwin, B. S. and Lomax, H. (1978). Thin Layer Approximation and Algebraic Model
for Separated Turbulent Flows. American Institute of Aeronautics and Astronautics,
Aerospace Sciences Meeting, 78–257.
[19] Spalart, P. R. and Allmaras, S. R. (1992). A One-Equation Turbulence Model for Aerodynamic Flows. American Institute of Aeronautics and Astronautics, 92–0439.
[20] Launder, B. E. and Spalding, D. B. (1974). The numerical computation of turbulent flows.
Computer Methods in Applied Mechanics and Engineering, 3 (2), 269–289.
[21] Shih, T. H., Liou, W. W., Shabbir, A., Yang, Z. and Zhu, J. (1995). A new eddy-viscosity
model for high Reynolds number turbulent flows–Model development and validation. Computers Fluids, 24 (3), 227–238
[22] Kato, M. and Launder, B. E. (1993). The Modeling of Turbulent Flow Around Stationary
and Vibrating Square Cylinders. Proc. 9th Symposium on Turbulent Shear Flows, Kyoto,
10.4.1–10.4.6.
[23] Fluent 6.3 (2006). Fluent User’s Guide. ANSYS, Inc.
[24] Shao, Y. (2000). Physics and modelling of wind erosion. Boston : Kluwer Academic,
Dordrecht.
[25] Crowe, C., Sommerfeld, M. and Tsuji, Y. (1998). Multiphase flows with droplets and
particles. CRC Press (Boca Raton, Fla) .
[26] Iversen, J. D., Greeley, R. and Pollack, J. B., (1976). Wind blown dust on Earth, Mars
and Venus. Journal of Atmospheric Science, 33, 2425-2429.
[27] Bagnold, R. A. (1941). The physics of blown sand and desert dunes. Journal of Atmospheric Science, London: Methuen, 265 pp.
[28] Greeley, R. and Iversen, J. D. (1985). Wind as a Geological Process. Cambridge University
Press, 99-100.
136
[29] Lu, H. and Shao, Y. (2001). Towards quantitative prediction of dust storms: An integrated
wind erosion modelling system and its application. Environmental Modelling and Software,
16 (3), 233–249.
[30] Owen, P. R. (1964). Saltation of uniform grains in air. Journal of Fluid Mechanics, 20 (2),
225-242.
[31] Nalpanis, P., Hunt, J. C. R.and Barrett, C. F. (1993). Saltating particles over flat beds.
Journal of Fluid Mechanics, 251, 661–685.
[32] Gillette, D. A. (1977). Fine particulate emissions due to wind erosion. Transactions of the
American Society of Agricultural Engineers, 20, 890–897.
[33] Marticorena, B. and Bergametti, G. (1995).Modeling the atmospheric dust cycle: 1. Design of a soil-derived dust emission scheme. Journal of Geophysical Research, 100 (D8),
415–430.
[34] Ishii, M. (1975). Thermo-fluid dynamic theory of two-phase flow. Paris: Evrolles
[35] Durst, F., Miloevic, D. and Schnung, B. (1984). Eulerian and Lagrangian predictions
of particulate two-phase flows: a numerical study. Journal of Instructional Computing
Research, 8 (2), 101–115.
[36] Shirolkara, J. S., Coimbrab, C. F. M. and Queiroz McQuay, M. (1984).Fundamental aspects of modeling turbulent particle dispersion in dilute flows. Progress in Energy and
Combustion Science, 22 (4), 363–399.
[37] Crowe, C. T., Troutt, T. R. and Chung, J. N. (1984). Numerical models for two-phase
turbulent flows. Annual Review of Fluid Mechanics, 28, 11–43.
[38] Manninen, M., Taivassalo, V., and Kallio, S. (1996). On the mixture model for multiphase
flow. VTT Publications 288, Technical Research Centre of Finland, 1996.
137
[39] Wirogo, S. (1997). Flux corrected method: An accurate approach to fluid flow modeling.
Ph.D Thesis, Iowa State University. AAT 9737770
[40] Schlichting, H., Gersten, K. and Gersten, K. (1979). Boundary-Layer Theory. McGraw
Hill, New York. ISBN: 978-3-540-66270-9.
[41] Cooper, D., Launder, D. C. and Liao, G. X. (1993)]. Impinging jet studies for turbulence
model assessment. Part I: Flow-field experiments. International Journal of Heat Mass
Transfer, 36, 2675–2684.
[42] Craft, T. J., Graham, L. J. W., and Launder, B. E., (1993). Impinging jet studies for turbulence model assessment. Part 2: An examination of the performance of four turbulence
models. International Journal of Heat and Mass Transfer, 36, 2685–2697.
[43] Rabbott, J. P. (1956). Static-Thrust Measurements of the Aerodynamic Loading on a Helicopter Rotor Blade. Technical Report, National Aeronautics and Space Administration,
Washington DC.
[44] Nathan, N. D. and Green, R. B. (2008). Measurements of a rotor flow in ground effect and
visualization of the brown-out phenomenon. Helicopter Society Annual Forum, Montreal,
Canada.
[45] Preston, J. R. (1994). VTOL Downwash/Outwash Operational Effects Model. Annual
Forum Proceedings, American Helicopter Society, 50 (V2),1101.
[46] Phillips, C. and Brown, R. E. (2008). The effect of helicopter configuration on the fluid
dynamics of brownout. 34th European Rotorcraft Forum, 16th - 19th Sept 2008, Liverpool,
England.
[47] Tennekes, H. and Lumley, J. L. (1972). A First Course in Turbulence. MIT Press, Cambridge, MA.
138
[48] Jones, W. P. and Launder, B. E (1972). The Prediction of Laminarization with a TwoEquation Model of Turbulence. International Journal of Heat and Mass Transfer, 15,
301–314
139
APPENDIX A.
WALL FUNCTIONS
At high Reynolds number, a very steep gradient of velocity exists for the flow near the
wall. Based on non-dimensional distance y + =
u∗ y
ν ,
the region near the wall can be divided
into three different layers [47].
1. viscous sub-layer (0 < y + < 5)
2. buffer layer (5 < y + < 30)
3. inertial sub-layer (30 < y + < 200)
To accurately capture the velocity gradient, a very refined mesh is required near the wall, which
can capture the viscous sub-layer. This increases computational requirements tremendously for
the high Reynolds number flow, where the viscous sub-layer is very thin. Also the k − ǫ Model
(and other high Reynolds number turbulence models) do not consider the effect of near-wall in
their modelling. So, these models cannot be integrated into the wall since they give inaccurate
predictions. To overcome this, the first node of the mesh is kept in an inertial sub-layer and
the other required variables are calculated using law of the wall or wall function.
A.1
Standard Wall Function
In the attached turbulent boundary layer, velocity profile in inertial sub-layer is given by
logarithmic profile [12] known as wall of law.
U
1
= ln(Ey + )
u∗
κ
(A.1)
140
where κ is van Karman constant equal to 0.41 and E is empirical constant equal to 9.793. If
Up is the velocity at the first node point in the inertial sub-layer, then wall friction can be
calculated as
τw =
ρu∗ Up
ln(Ey + )
(A.2)
The wall friction calculated through the wall function approach is implemented in the
existing discretized equation by modifying the wall viscosity, µe , which can be calculated as
τw = µ
ρu∗ Up
ρu∗ yp κ
Up
Up
∂U
= µe
⇒ µe
=
⇒ µe =
+
∂y
yp
yp
ln(Ey )
ln(Eyp+ )
A.2
(A.3)
Launder Spalding Wall Function
One of problems with the standard wall function is the u∗ goes to zero at separating and
re-attachment point, thus causing singularity. Launder and Spalding [20] proposed to use
1/4 √
identity u∗ = Cµ
k in calculating effective wall viscosity.
1/4 p
ρCµ
k p yp κ
µe =
ln(Eyp∗ )
where
yp∗
=
1/4 p
yp C µ
kp
ν
(A.4)
(A.5)
In the current research, the wall function is applied for the wall adjacent node with yp∗ >
11.225 as used in commercial software, Fluent [23].
A.3
Turbulence
Based on the local equilibrium hypothesis in the wall adjacent cell, production of turbulent
kinetic energy and dissipation are assumed equal. Accordingly, turbulent production in the
wall adjacent cell is calculated as [23]
Pk ≈ τ w
τw
∂U
= τw
p
1/4
∂y
κρCµ
k p yp
(A.6)
141
and the turbulent dissipation in wall adjacent node is calculated as
3/4 3/2
ǫp =
A.4
Cµ kp
κyp
(A.7)
Wall Function for Rough Wall
Effective viscosity for the rough wall is given by
1/4 p
k p yp κ
ρCµ
µe =
∗
ln(Eyp ) − ∆B
(A.8)
where ∆B is the roughness function which depends upon non-dimensional roughness height
Ks+ . Ks+ is given by
Ks+ =
ρKs u∗
µ
(A.9)
where Ks is the physical roughness height.
The wall is considered smooth when Ks+ < 2.25. In this case, ∆B = 0. For a fully rough
regime, when Ks+ > 90,
∆B =
1
ln(1 + Cks Ks+ )
κ
where Cks is roughness constant which varies from 0.5 to 1.0.
(A.10)