Detached Eddy Simulations of a simplified Tractor-Trailer Geometry
Except where reference is made to the work of others, the work described in this thesis
is my own or was done in collaboration with my advisory committee. This thesis does
not include proprietary or classified information.
Harshavardhan Ghuge
Certificate of Approval:
Brian Thurow
Assistant Professor
Aerospace Engineering
Christopher Roy, Chair
Assistant Professor
Aerospace Engineering
Anwar Ahmed
Associate Professor
Aerospace Engineering
Joe F. Pittman
Interim Dean
Graduate School
Detached Eddy Simulations of a simplified Tractor-Trailer Geometry
Harshavardhan Ghuge
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
May 10, 2007
Detached Eddy Simulations of a simplified Tractor-Trailer Geometry
Harshavardhan Ghuge
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at their
expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Harshavardhan Ghuge, son of Arvind and Vaishali Ghuge was born on July 02, 1981,
in Pune, India. He attended Maharashtra Institute of Technology, Pune and graduated
in June 2004 with the degree of Bachelor of Mechanical Engineering. He joined the
Masters program in the department of Aerospace Engineering at Auburn University in
August 2004.
iv
Thesis Abstract
Detached Eddy Simulations of a simplified Tractor-Trailer Geometry
Harshavardhan Ghuge
Master of Science, May 10, 2007
(B.E., Maharashtra Institute of Technology, Pune University, 2004)
116 Typed Pages
Directed by Christopher J. Roy
This thesis presents a comparison of the results from the Detached Eddy Simulation
(DES) turbulence model used for modeling flow over a 1/8th scale simplified model of a
tractor-trailer geometry (or the Ground Transportation System) with the experimental
data. The main aim of this study is to assess the DES model for validation purpose.
From the results it is concluded that the DES model gives very good agreement with
the experimental data for the global quantities like drag but does not predict the details
of the wake structure correctly. The reason for the inaccuracy in the prediction may
be originating from the lack of knowledge about the transfer of information between
the Reynolds-Averaged Navier-Stokes region (where the turbulence is modeled) and the
Large Eddy Simulation region (where the turbulence is simulated).
v
Acknowledgments
I would like to acknowledge the U.S. Department of Energy FreedomCAR and Ve-
hicle Technologies (CAR stands for Cooperative Automotive Research) for funding this
work at Auburn University. I would like to thank my adviser Dr. Christopher J. Roy,
Assistant Professor, Aerospace Engineering Department for his guidance and support
throughout the research and toward the completion of this thesis.
Many thanks to Dr Matthew F. Barone of the Aerosciences and Compressible Fluid
Mechanics Department from Sandia Laboratories for his insightful guidance during this
research. I would also like to thank Dr. Brian Thurow, Assistant Professor, Aerospace
Engineering Department, for his assistance during my research and for serving as a
committee member. I am thankful to Dr. Anwar Ahmed, Associate Professor, Aerospace
Engineering Department for serving as a committee member.
Most importantly, I am grateful to my mother for patiently supporting me through-
out my education.
vi
Style manual or journal used LATEX: A Document Preparation System by Leslie
Lamport (together with the style known as ?aums?)
Computer software used TEX (specifically LATEX), MATLAB 7.0.4, Tecplot and
the departmental style-file aums.sty
vii
Table of Contents
List of Figures x
List of Tables xiv
Nomenclature xv
1 Introduction 1
1.1 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Introduction to Turbulence Models . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Reynolds Averaged Navier-Stokes Equations . . . . . . . . . . . . . 4
1.2.2 Introduction to Large Eddy Simulation (LES) . . . . . . . . . . . . 7
1.2.3 Introduction to Hybrid RANS/LES Models . . . . . . . . . . . . . 8
1.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.1 Experimental Research . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3.2 Computational Research . . . . . . . . . . . . . . . . . . . . . . . . 13
2 Problem Description 24
2.1 Wind Tunnel Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1.1 Surface Pressures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.1.2 Drag, Lift and Side Force Measurements . . . . . . . . . . . . . . . 27
2.1.3 Oil-film Interferometry (OFI), Hot-film Anemometry and Particle
Image Velocimetry (PIV) . . . . . . . . . . . . . . . . . . . . . . . 27
3 Simulation Approach 28
3.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 The DES Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 The SACCARA Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4 Discretization of the Convective Terms . . . . . . . . . . . . . . . . . . . . 33
3.4.1 Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4.2 Upwind Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5 Discretization of the Diffusion Terms . . . . . . . . . . . . . . . . . . . . . 35
3.6 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.7 Grid Generation Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.8 Domain Decomposition and Load Balancing . . . . . . . . . . . . . . . . . 36
3.9 Boundary Conditions used for the Simulations . . . . . . . . . . . . . . . 37
3.10 Characteristic Scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
viii
4 Numerical Accuracy 43
4.1 Statistical Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2 Initial Transients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 Effect of Number of Sub-iterations used for Iterative Convergence . . . . . 57
4.4 Effect of Time Step (??s) . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.5 Effect of the Mesh Size on DES . . . . . . . . . . . . . . . . . . . . . . . . 63
5 Results and Discussion 67
5.1 Contour and Streamline Plots . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2 Velocity Profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.3 Pressure Distribution Profile . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.4 Drag, Lift and Side Force Coefficients . . . . . . . . . . . . . . . . . . . . 79
5.5 Vortex Core Locations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.6 Turbulent Eddies in the Wake . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.7 Power Spectra of the Unsteady Pressure Signal . . . . . . . . . . . . . . . 83
6 Conclusion 88
7 Recommendations 91
Bibliography 92
Appendices 97
A Fortran Code for creating animations using Tecplot 98
ix
List of Figures
2.1 The GTS model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.1 Reference static pressure comparison between the coarse mesh simulations
and the experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2 Coarse mesh structure near the truck base . . . . . . . . . . . . . . . . . . 41
3.3 Medium mesh structure near the truck base . . . . . . . . . . . . . . . . . 42
4.1 Time windows used for calculating statistically converged results . . . . . 45
4.2 Location of pressure and velocity profiles used for statistical convergence . 46
4.3 Pressure distribution along line T and corresponding to different time
windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4 Velocity profiles along line A and corresponding to different time windows 48
4.5 Error in the pressure distribution obtained from smaller time windows
with respect to the pressure distribution obtained from the biggest time
window . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.6 Error in the u/Vref velocity profiles obtained from smaller time windows
with respect to the u/Vref velocity profile obtained from biggest time
window . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.7 Error in the v/Vref velocity profiles obtained from smaller time windows
with respect to the v/Vref velocity profile obtained from biggest time window 49
4.8 Error in the urms/Vref velocity profiles obtained from smaller time win-
dows with respect to the urms/Vref velocity profile obtained from biggest
time window . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.9 Time windows used for estimating period of initial transients . . . . . . . 52
4.10 Pressure distribution along line T and corresponding to different time
windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
x
4.11 Velocity profiles along line A and corresponding to different time windows 54
4.12 Error in the pressure distribution from all time windows with respect to
the pressure distribution obtained from the window farthest along time-axis 55
4.13 Error in the u/Vref profiles from all time windows with respect to the
u/Vref profile obtained from the window farthest along time-axis . . . . . 55
4.14 Error in the v/Vref profiles from all time windows with respect to the
v/Vref profile obtained from the window farthest along time-axis . . . . . 56
4.15 Error in the urms/Vref profiles from all time windows with respect to the
urms/Vref profile obtained from the window farthest along time-axis . . . 56
4.16 Pressure distribution along line T and corresponding to different number
of sub-iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.17 Velocity profiles along line A and corresponding to different number of
sub-iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.18 Pressure distribution along line T and corresponding to different time-steps 61
4.19 Velocity profiles along line A and corresponding to different time-steps . . 62
4.20 Medium mesh vertical streamwise contour plot . . . . . . . . . . . . . . . 64
4.21 Coarse mesh vertical streamwise contour plot . . . . . . . . . . . . . . . . 64
4.22 Pressure distribution along line T obtained from the coarse and medium
mesh simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.23 Velocity profiles along line A obtained from the coarse and medium mesh
simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.1 Planes used for comparison of velocity data . . . . . . . . . . . . . . . . . 67
5.2 Vertical streamwise PIV window in the experiment . . . . . . . . . . . . . 68
5.3 Vertical streamwise contour plot from the medium mesh simulations . . . 68
5.4 Horizontal streamwise PIV window in the experiment (y/W = 0.696) . . 69
5.5 Horizontal streamwise contour plot from the medium mesh simulations
(y/W = 0.709) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
xi
5.6 Horizontal streamwise PIV window in the experiment (y/W = 1.044) . . 70
5.7 Horizontal streamwise contour plot from the medium mesh simulations
(y/W = 1.05) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.8 Comparison of u/Vref profiles from the experiment and simulations at
y/W = 0.348 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.9 Comparison of u/Vref profiles from the experiment and simulations at
y/W = 0.696 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.10 Comparison of u/Vref profiles from the experiment and simulations at
y/W = 1.044 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.11 Comparison of v/Vref profiles from the experiment and simulations at
y/W = 0.348 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.12 Comparison of v/Vref profiles from the experiment and simulations at
y/W = 0.696 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.13 Comparison of v/Vref profiles from the experiment and simulations at
y/W = 1.044 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.14 Comparison of w/Vref profiles from the experiment and simulations at
y/W = 0.348 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.15 Comparison of w/Vref profiles from the experiment and simulations at
y/W = 0.696 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.16 Comparison of w/Vref profiles from the experiment and simulations at
y/W = 1.044 and z/W = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.17 Comparison of pressure distribution along base of the truck from the ex-
periment and the medium mesh simulations . . . . . . . . . . . . . . . . . 77
5.18 Comparison of pressure distribution on the front surface of the GTS model 78
5.19 Comparison of pressure distribution on the top surface of the GTS model 78
5.20 Comparison of pressure distribution along the right side of the GTS model 78
5.21 Comparison of pressure distribution on the bottom surface of the GTS
model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
xii
5.22 Coarse mesh mean and instantaneous drag coefficient . . . . . . . . . . . . 80
5.23 Medium mesh mean and instantaneous drag coefficient . . . . . . . . . . . 80
5.24 Location of vortex cores in vertical streamwise plane (z/W = 0.0) . . . . . 81
5.25 Location of vortex cores in horizontal streamwise plane (y/W = 0.696) . . 81
5.26 Location of vortex cores in horizontal streamwise plane (y/W = 1.044) . . 81
5.27 Vortical structures resolved using LES in the wake of the GTS on the
medium mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.28 Fast Fourier transform of the pressure signal from experimental data (x/W
= 7.647, y/W = 0.63 and z/W = -0.46) . . . . . . . . . . . . . . . . . . . 85
5.29 Fast Fourier transform of the pressure signal from coarse mesh simulations
(x/W = 7.647, y/W = 0.63 and z/W = -0.46) . . . . . . . . . . . . . . . 86
5.30 Fast Fourier transform of the pressure signal from medium mesh simula-
tions (x/W = 7.647, y/W = 0.63 and z/W = -0.46) . . . . . . . . . . . . 86
5.31 Fast Fourier transform of the pressure signal in the horizontal streamwise
plane from the coarse mesh simulations . . . . . . . . . . . . . . . . . . . 86
5.32 Fast Fourier transform of the pressure signal in the horizontal streamwise
plane from the coarse mesh simulations . . . . . . . . . . . . . . . . . . . 86
5.33 Fast Fourier transform of the pressure signal from the vertical streamwise
plane in the coarse mesh and near the bottom surface of the GTS model . 87
5.34 Fast Fourier transform of the pressure signal from the vertical streamwise
plane in the coarse mesh and near the top surface of the GTS model . . . 87
xiii
List of Tables
1.1 Effect of drag-reduction concepts on wind-averaged drag [5]. . . . . . . . 11
2.1 Wind tunnel conditions [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Experimental values of pressure coefficient at different z/W locations
along the y/W axis of the trailer base [4]. . . . . . . . . . . . . . . . . . . 26
3.1 Coarse and medium mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2 Reference pressure comparison . . . . . . . . . . . . . . . . . . . . . . . . 38
3.3 Characteristics scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.1 Absolute error with respect to biggest time window . . . . . . . . . . . . . 50
4.2 Absolute error with respect to time window farthest along the time-axis . 56
5.1 Total uncertainty values of the velocity and pressure from the simulations 74
5.2 Comparison of the drag, side and lift force coefficients . . . . . . . . . . . 79
xiv
Nomenclature
English Symbols
b distance between square cylinder bottom face and the channel floor (m)
c cth processor
cb1, cb2, cv1, cw1
cw2, cw3 closure coefficients in Spalart-Allmaras model
CD coefficient of drag
CDES detached eddy simulation constant
CDmean mean value of coefficient of drag
Ck coefficient of thermal conductivity (W/m/K)
CL coefficient of lift
CLmean mean value of coefficient of lift
Cp time-averaged coefficient of pressure
cp specific heat at constant pressure (J/kg/K)
(Cp)error error in the Cp profile
CS coefficient of side force
xv
Chapter 1
Introduction
Class 8 tractor-trailers carry the major load of transportation in the trucking indus-
try and are used for transportation of most of the commodities across the United States
of America. The weight of a modern Class 8 tractor-trailer can be up to 80,000 pounds.
At a common highway speed of 70 mph 65% of the total engine output is consumed in
overcoming the aerodynamic drag. Also the energy used in overcoming the aerodynamic
drag increases as the vehicle speed increases [1]. During late 1970?s and early 1980?s
improvement in truck aerodynamics to reduce fuel consumption led to the development
of the following technologies [2]:
1. Cab shaping
2. Cab-mounted deflectors
3. Trailer front-end fairings
4. Cab side extenders
5. Body front-edge rounding
6. Tractor-trailer gap seals
7. Trailer side skirts
8. Rear-boat-tailing
1
The first five items are called the first generation drag reduction devices and are already
in use. The last three items are called second generation drag reducing devices and are
not used widely. A typical Class 8 tractor-trailer has a wind-averaged drag coefficient of
CD = 0.6. A wind-averaged drag coefficient is a drag coefficient that takes into account
the effects of wind blowing from all the directions as is the case when a vehicle is moving
on the road. The present day truck drag coefficients can be reduced by as much as 50%
[1].
According to the statistics presented in [3] there were 2.2 million trucks on the
US roads in the year 2003. Each truck traveled approximately 62,900 miles with an
average fuel efficiency of 5.3 miles/gallon. A 50% reduction in truck drag coefficient
would give nearly 25% reduction in fuel consumption at highway speeds. This could
save fuel worth $1.5 billion per year [1]. Lower fuel consumption will also result in a
reduction in pollution emissions and a reduced dependence on foreign oil.
A significant reduction in cost can be achieved by replacing some of the road and
wind-tunnel testing required for the aerodynamic design and analysis of the Class 8
tractor-trailer with computational simulations. At the same time, one should view the
computational results with caution because not all the flow physics may be completely
captured by the computational fluid dynamics (CFD) simulations and the drag prediction
may differ from the true values. The US Department of Energy has developed a multi-
year project in which it has funded a number of organizations and universities with the
aim to establish a clear understanding of the drag producing flow phenomena over a Class
2
8 tractor-trailer through joint experiments and computations. Eventually this research
will help in the design of ?smart? drag reducing devices for this class of tractor-trailers [1].
One aim of this project is to provide high-quality aerodynamic data for direct validation
of computational tools used for simulation of the flow over a Class 8 tractor-trailer. Use
of the validated computational tools will help reduce the experimental cost in the future.
As a first step of the validation study a very simplified model of a Class 8 tractor-trailer
called the Ground Transportation System (Figure 2.1) was studied experimentally at
NASA Ames [4]. A detail explanation of this experimental study isprovided inChapter 2.
To perform validation studies, computational simulations of the Ground Transportation
System (GTS) were performed using the Detached Eddy Simulation turbulence model [8].
The simulation approach is described in Chapter 3. Comparison with the experimental
data and the results will be discussed in Chapters 4 and 5.
1.1 Research Objective
This work will compare the results from the computational simulations of the sim-
plified tractor-trailer model (or the GTS model) with the experimental data to assess
the Detached Eddy Simulation model as a CFD tool to simulate the flow over a Class
8 tractor-trailer. The comparison will include the drag coefficient of the flow over the
truck surface, details of the flow structure in the wake of the truck, pressure coefficient
distribution at various locations on the truck surface and the frequency of the vortices
shed from the truck base.
3
1.2 Introduction to Turbulence Models
1.2.1 Reynolds Averaged Navier-Stokes Equations
From the engineering point of view we are usually interested in average or mean
quantities while analyzing turbulent flows. Since turbulence can be characterized by
random chaotic motion of the fluid, we rely on statistical approaches to calculate the
averaged quantities. The Navier-Stokes equations are mass-averaged using a statistical
approach called Favre (mass) averaging [6]. A Favre-averaged velocity vector is given by
?vi = 1?? lim?t??
integraldisplay ?t
0
?vidt (1.1)
where ?vi is the Favre-averaged velocity vector, ?? is the mean density and vi is the instan-
taneous velocity vector. So the Favre-averaged continuity equation (tensor notation) in
conservation form is given by
???
?t +
?(???vi)
?si = 0, (1.2)
the Favre-averaged momentum equation in conservation form is given by
?(???vi)
?t +
?(???vj ?vi)
?sj =?
??p
?si +
?
?sj
bracketleftBig
?tji ??v??j v??i
bracketrightBig
(1.3)
4
and the Favre-averaged energy equation in conservation form is given by
?
?t
bracketleftBigg
??
parenleftbigg
?e+ ?vi?vi2
parenrightbigg
+ ?v
??
i v
??
i
2
bracketrightBigg
+ ??s
j
bracketleftBigg
???vj
parenleftbigg
?h+ ?vi?vi
2
parenrightbigg
+ ?vj?v
??
i v
??
i
2
bracketrightBigg
= ??v
j
bracketleftbigg
?qLj ??v??j h?? +tijv??i ? 12?v??j v??i v??i
bracketrightbigg
+ ??s
j
bracketleftBig
?vi
parenleftBig
?tij ??v??i v??j
parenrightBigbracketrightBig
(1.4)
Equations 1.2 to 1.4 are called the Reynolds-Averaged Navier-Stokes (RANS) equations
[6]. It should be noted that they are obtained by performing Favre-averaging of the
Navier-Stokes equations, but it is conventional to call them the RANS equations. The
mean pressure is given by using the mean density and temperature in the perfect gas
law
?p = ??R?T (1.5)
There are some unknown terms in the RANS equations. Any term in the RANS
equations which involves the product of two or more Favre fluctuating variables (v??i and
h??) is an unknown term because the exact correlation between two Favre fluctuating
variables is not known. These terms are approximated by solving additional partial
differential and algebraic equations. Depending on the number of partial differential
equations used to approximate these unknown terms, the RANS models are divided into
four categories [6].
5
Algebraic Models
The algebraic models use algebraic equations to solve for the unknown terms in
the RANS. Some examples of these kind of models are the Cebeci-Smith model and the
Baldwin-Lomax model [6].
One-Equation Models
One-equation models use one partial differential equation and some algebraic equa-
tions to approximate the unknown terms. Some examples of these kinds of models are
the Spalart-Allmaras model and the Baldwin-Barth model [6].
Two-Equation Models
Two-equation models use two partial differential equations and some algebraic equa-
tions to approximate the unknown terms. Of all the RANS models two-equation models
are the most popular for the industrial applications. The k?? and k?epsilon1 models are the
most popular two-equation models [6].
Second-Order Closure Models
The effects of stream line curvature, sudden changes in strain rate, secondary mo-
tion, etc. cannot be accommodated by two equation models. These features are handled
by second order closure models. For three-dimensional flows, a second order closure
model introduces seven new partial differential equations. Some examples of these kind
6
of models are the Launder-Reece-Rodi model and the Wilcox multiscale model. How-
ever, because of the large number of equations and complexity involved in second-order
closure models they are not as popular as the two-equation model and the algebraic
models [6].
1.2.2 Introduction to Large Eddy Simulation (LES)
The RANS equations model the effect of the turbulent eddies on the flow, but they
do not resolve the turbulent eddies. Large Eddy Simulation (LES) is other turbulence
modeling technique in which the larger eddies are resolved and the smaller ones are mod-
eled. Resolving means the Navier-Stokes equations are numerically solved without any
kind of approximation. LES is computationally expensive but is usually more accurate
than RANS. Two important steps in LES are as follows [6]:
1. Filtering
2. Sub-grid scale modeling
Filtering separates the large eddies (which are resolved) from the smaller eddies (which
are modeled). It introduces a scale ? which represents the smallest turbulent scale or
eddy that can be resolved by the LES model. Any scale or eddy smaller than ? is modeled
using a sub-grid scale model. Commonly used filters in LES are the box function, the
Fourier cutoff function and the Gaussian function [6]. Some examples of sub-grid scale
models are the Smagorinsky eddy viscosity model [41] and the Dynamic one-equation
model [40].
7
1.2.3 Introduction to Hybrid RANS/LES Models
The biggest drawback of LES simulations for wall-bounded flows is that they are
prohibitively expensive because of the very fine mesh requirements of the LES model
in the boundary layer. To overcome this drawback the LES model is combined with a
RANS model to give the hybrid RANS/LES modeling approach. In this approach RANS
modeling is used in the attached boundary layer region and LES is used away from the
boundary layer and in the highly separated flow regions [20]. Examples of the hybrid
RANS/LES models are the Detached Eddy Simulation model (DES) [8] and the Linear
Numerical Scale model (LNS) [7]. A detailed explanation of the DES model is given in
Section 3.2.
1.3 Literature Review
In this section a review of some experimental and computational research performed
on bluff body flows such as the GTS model, a square cylinder and a surface mounted
cube is presented. The cases of the flow over a square cylinder and a surface mounted
cube are lower Reynolds number cases and the case of the flow over the GTS model is
a high Reynolds number case. This section will also compare the results from various
CFD models used for modeling the flow over bluff bodies. Some new hybrid RANS/LES
models are also reviewed.
8
1.3.1 Experimental Research
Cooper [2] has reviewed the past work done in the field of truck aerodynamics for
drag reduction purpose. He has described the design aspects and the potential of drag
reduction in the second generation drag reducing devices. Based on the data available
he has concluded in his work that the combination of first and second generation devices
can reduce fuel consumption by over 4,000 US gallons for a tractor-trailer operating
125,000 miles/year. A net reduction in wind-averaged drag coefficient based on a highway
speed of 65 mph ? ?CD(65) = 0.183 is possible with a combination of first and second
generation devices. The second generation devices are mechanically more complex than
the first-generation aerodynamic package and their implementation will require efficient
mechanical design and government action to encourage their use.
High Reynolds Number Flow Cases
The Reynolds number for a Class 8 tractor-trailer traveling on a highway at an
average speed of 70 mph is approximately 6 million. So the flow over tractor-trailer
is considered a high Reynolds number case. Effect of Reynolds number on the drag
coefficient of a tractor-trailer is negligible for Reynolds number greater than 2 million.
McCallen et al. [1] have stressed the importance of the understanding of the flow
structure in the wake and the flow in the gap between the tractor and the trailer for opti-
mized design of the drag reducing devices. A brief summary of the suitability of various
turbulence models for the computational simulations of the flow over the tractor-trailer
9
is discussed showing the limitations of the RANS models. Results from the experimen-
tal investigations of the drag reducing devices are discussed. It is concluded that these
devices can provide up to 11% to 12% of fuel savings. The authors have mentioned the
need to make the trailer manufacturers realize the potential of the drag reducing devices
fitted on the trailers and a need to design these devices so that their maintenance cost
is minimum.
In an experimental paper by Storms et al. [5] two studies were performed using the
Generic Conventional Model (GCM). A GCM is a model of a Class 8 tractor-trailer with
all the details of the shape of the tractor-trailer body.
1. The effect of variation of Reynolds number on the aerodynamic drag acting on the
tractor-trailer.
2. The variation in the aerodynamic drag with use of various drag reducing devices.
In the first study it was found that the effect of the Reynolds number was most evident
at high yaw angles (? > 8? and ? < -8 ?) where there was a significant reduction in
drag at lower Reynolds numbers (Reynolds number was varied from 1 to 6.5 million).
In the second study six drag reduction devices were studied experimentally. They are
listed below.
1. Side and top extenders
2. Boattail plates
3. Base flaps
10
4. Trailer belly box
5. Trailer side skirts
6. Trailer full skirt
Table 1.1 shows the effectiveness of drag-reduction devices measured relative to the
extender configuration. The positive sign indicates an increase in the wind averaged
drag and the negative sign indicates a decrease in the drag value. The wind averaged
drag reported in the paper was computed for a highway speed of 55 mph.
Table 1.1: Effect of drag-reduction concepts on wind-averaged drag [5].
Flow ? ?CD, 55 mph % difference
No Extenders +0.156 +37.0
Extenders 0 0
Boattail plates -0.058 -13.7
Base flaps -0.082 -19.4
Trailer belly box -0.050 -11.8
Trailer side skirts -0.026 -6.2
Trailer full skirt +0.016 +3.8
Low Reynolds Number Flow Cases
The lowReynoldsnumber flow caseshavethe Reynoldsnumber lessthan onemillion.
The effect of wall proximity on flow in the wake of a square cylinder was analyzed by
Martinuzzi et al. [10] at a Reynolds number of 18,900 based on the diameter. They
observed four different flow regimes which are related to changes in the gap flow between
the cylinder bottom face and the neighboring wall. For b/D > 0.9 the back pressure,
the drag coefficient and the strength of the shed vortices in the wake of the cylinder are
11
insensitive to the gap height. For 0.6 < b/D < 0.9, the strength of the shed vortices
and fluctuating lift decrease while the base pressure increases as b/D is reduced. For 0.3
< b/D < 0.6 the flow reattaches intermittently on the bottom faces, which causes the
shedding of vortices to become increasingly irregular. For b/D < 0.3 the lower shear layer
reattaches permanently to the bottom face and the periodic fluctuations are completely
suppressed.
Martinuzzi and Tropea [11] studied the flow structures around obstacles with dif-
ferent aspect ratios placed in a flow channel. The Reynolds number based on the height
of the channel was varied from 80,000 to 120,000. For aspect ratios (w/H) > 6 the
cross-stream component of velocity and pressure gradient are negligible in the middle of
the wake showing that this region can be treated as two-dimensional. It is shown that
the total pressure loss due to the presence of the object is proportional to the blockage
ratio. It is concluded that the recovery length behind a three dimensional obstacle is
shorter than in the case of a two-dimensional flow.
Lyn et al. [12] studied experimentally the flow over a square cylinder at Reynolds
number of Re = 21,400 based on the cylinder diameter. Ensemble-averaged statistics
of the flow around the square cylinder are obtained using two component laser-Doppler
measurement. A distinction is made between the flow in the near wake region where the
vortices are ?mature? and have a distinct identity and the base region where the vortices
grow to maturity and are then shed.
12
1.3.2 Computational Research
High Reynolds Number Flow Cases
Roy et al. [21] have performed RANS simulations on the GTS model at 0? yaw angle
using the baseline Menter k?? model. The quantities compared with the experimental
data [4] include vehicle drag, surface pressure and the time-averaged velocities in the
trailer near wake. The prediction of pressure and the time-averaged velocity in the trailer
wake is in poor agreement with the experimental data due to the inability of the RANS
models to accurately capture the near-wake vortical structure. This is one motivation
for using the DES model (which actually resolves the larger eddies) to simulate the flow
over the GTS model.
Maddox et al. [18] performed DES of the GTS model at 0? and 10? yaw angles
and compared his results with experiment performed by Storms et al. [4]. The compu-
tational domain extended from approximately 15W upstream of the GTS up to 13.5W
downstream of the GTS. The dimensionless time step (non-dimensionalized using the
width W of the GTS model and upstream velocity 93.9 m/s) used was 0.02. For the
10? yaw angle case the reattachment of the separation bubble on the leeward side of
the model occurred substantially further along the GTS than indicated by the pressure
measurements made by Storms [4]. The statistics were acquired for a non-dimensional
time of 81 units. The present work does not agree with this time-window used by Mad-
dox because from the research presented in this thesis the smallest time-window required
for getting statistically converged results is ??w ? 180 units. For the 0? yaw angle the
13
drag coefficient was calculated to be 0.279 as against the experimental value of 0.249.
For the 10? yaw angle case the drag coefficient from the simulations was 1.379 and the
experimental value was 1.253. No details of the wake flow have been shown or described.
The authors have attributed the discrepancy in the force predictions at 10 degree yaw
to the transition from laminar to turbulent flow in the boundary layers.
Unaune et al. [29] have performed DES on the simplified tractor-trailer at 0? yaw
angle using commercial modeling software Fluent 6.1 and compared their results with
the experiment performed by Storms et al. [4]. Unstructured mesh was used in the
simulations. Wall functions were used to model the boundary layer and the first cell
center was 0.75 mm from the truck surface. The grid consisted of 13.8 million cells.
The time step used in the simulations were 0.001 s for first 0.08 s and 0.0005 s for the
remaining simulation period. The estimated period of initial transients was 0.2 s and the
statistics were collected for a time window of 0.2 s to 0.55 s. The research presented in
this thesis does not agree with this and shows that the initial transience period is nearly
equal to 1 s and the time-window required for obtaining statistically converged results is
nearly equal to 0.6 s. The details of the wake given by Unaune show that the results are
not statistically converged. The drag coefficient computed from the simulations was CD
= 0.253 while the experimental value of the drag coefficient is CD = 0.249. In the power
spectra of the unsteady pressure sensor located in the base of the trailer there were no
obvious peaks showing any particular shedding frequency.
14
Baysal and Bayraktar [30] have used the renormalization group (RNG) k?epsilon1 model
for the simulation of flow over a Generic Conventional Model (GCM). The computations
were performed on the left half of the computational domain with proper boundary
conditions at the symmetry plane. The fine mesh consisted of about 5 million cells.
On the ground boundary the velocity of the flow was set equal to the flow at the inlet
boundary to create a moving ground boundary condition. The flow over the tractor-
trailer without tires and with stationary ground plane was also computed. The flow
in the gap between the tractor and trailer was altered with removal of the tires. The
computed drag value is about 6% less as a result of removing tires. For the stationary
ground case the boundary layer on the ground thickened to alter the entire undercarriage
flow. Also the trailer wake become skewed and was driven towards the ground. The
computed drag value was 9% more as a result of making the ground stationary.
Low Reynolds Number Cases
Ramesh et al. [13] have performed 3D unsteady RANS simulations of the flow over
a square cylinder at Reynolds number of Re = 22,000 based on the cylinder diameter.
The RANS model used was a non-linear k ?epsilon1 model [48]. Unlike the standard k ?epsilon1
model the non-linear model uses non-linear terms in the calculation of Reynolds stresses.
Although the results of the bulk parameters like the reattachment length, mean drag and
the Strouhal number are inferior to results from LES simulations by Rodi [14] and by
Nakayama and Vengadesan [15] and they do not produce same level of agreement with
experimental data the prediction of the mean quantities and the unsteady phenomenon
15
is better than the linear eddy viscosity models. On the whole when comparing the
computational time required by the present model and by the LES simulations, the
accuracy achieved by the unsteady non-linear k?epsilon1 model is significant and hence this
model can be used to simulate 3D unsteady complex engineering flows and reasonable
success can be achieved in its application.
Roy et al. [20] have examined the flow over a square cylinder using steady state
RANS methods and DES. The RANS modeling was performed on the vertical symme-
try plane of the square cylinder. For the DES the numerical issues examined include
statistical convergence, iterative convergence, temporal discretization error and spatial
discretization error. All the RANS models greatly over predicted the length of the recir-
culation zone. The models gave poor predictions of the mean and fluctuating velocities
in the wake. The results from the DES gave good agreement with the experimental data
for the global quantities and the wake mean and rms velocities.
Rodi et al. [14] had organized a workshop on LES simulations of a square cylinder at
a Reynolds number ofRe = 22,000 and a surface mounted cube at Reynolds number ofRe
= 3000 and Re = 22,000. In the conclusion the authors have mentioned that the results
from the simulations of the flow over a cube were in better agreement with each other
and with the experimental data because the flow over a cube is fully turbulent, unlike the
square cylinder case which involved transitional flow. The predictions produced by the
use of different sub-grid scale models are different and the differences are larger at high
Reynolds number. The dynamic one equation model gave better spatial distribution of
16
the sub-grid eddy viscosity, but there is no guarantee of getting overall better results
relative to Smagorinsky model.
L?ubcke et al. [22] have devised a RANS model called the explicit algebraic stress
model (EASM) and applied it to the square cylinder at Reynolds number of Re = 22,000
and the circular cylinder case at a Reynolds number of Re = 3900 and Re = 140,000. The
EASM model shows improvements in prediction of unsteady flows at low computational
cost as compared to other RANS models. The results predicted by the EASM matched
fairly well when compared with the LES simulations and the experimental data.
Krajnovi?cand Davidson[23] performedLESsimulations ofthe surface mounted cube
at a Reynolds number of Re = 40,000. The sub-grid scale models used were the dynamic
one equation model (OEM) proposed by Davidson [24] and the localized dynamic sub-
grid scale kinetic energy equation model proposed by Menon and Kim [25]. A series
of time-averaged velocities and turbulent stresses are computed and compared with the
experimental data. The results of velocities obtained from the simulation are generally in
better agreement with the experiment than the results obtained for the Reynolds stresses.
Predictions made without any model give poor agreement with the experimental data
but the two sub-grid models give good agreement with the experiment. It was observed
that the grid refinement has greater effect on the separation length than the reattachment
length. The Strouhal number is obtained from the Fourier transformation of the side
force signal. The Strouhal number for the fine mesh using the OEM model was 0.146
and best agreed with the experimental result of 0.145. The authors concluded that LES
17
proposed in the paper using simple inlet boundary conditions, on a relatively coarse grid
with the two sub-grid scale models gave accurate results at an acceptable computational
cost.
Sohankar et al. [26] examined results from different sub-grid scale models used in
LES simulations of flow around a square cylinder at Reynolds number of Re = 22,000.
The sub-grid scale models used were the Smagorinsky model (SSM) [41], the standard
dynamic model (DSM) [42] and the dynamic one-equation model [40]. When blockage
was taken into account the DSM model gave the worst results when compared with the
experiment. The SSM model produced similar results to those produced by the dynamic
one-equation model except for the Strouhal number and the root mean square (rms)
drag. Effects of spatial and temporal resolution and the computational span wise length
were analyzed using the dynamic one-equation model. Finer spatial resolution improved
agreement with the experiment. When the spanwise dimension was increased from four
to seven diameters there was a 6% reduction in sectional rms drag. When the time
resolution was increased by a factor of 2 the rms lift reduced by 5%. The one-equation
sub-grid model was successful in accounting for the backscattering phenomenon. Among
the three models the lowest level of all components of Reynolds stress is predicted by
the DSM. It was concluded that lower Reynolds stresses correspond to higher pressure
region in the wake flow which leads to lower drag force. The results produced by the
dynamic one-equation model give better agreement with experiment than the other two
subgrid models.
18
Barone and Roy [28] used the DES model for simulation of the bluff body wake
flow. The authors have used a combination of Yee?s [32] second order symmetric total
variation diminishing (STVD) scheme along with characteristic-based filtering for solving
the Navier-Stokes equations. The characteristic-based filtering was obtained using the
artificial compression method (ACM) switch of Harten [33]. The first case analyzed
was flow in the wake of a square cylinder at a Reynolds number of Re = 21,400. Out
of the coarse, medium and fine grids, results from the fine grid are mostly superior
than other grids, except the rms drag fluctuation and prediction of the stream-wise
velocity fluctuation in the near wake region. The DES model with the present numerical
scheme is not able to give accurate predictions of Reynolds shear stress in near wake.
In the summary for this case the fine grid DES results are found competitive with LES
calculations in the prediction of global quantities. The results of the coarse mesh using
the low dissipation scheme (combination of STVD and ACM switch) gave nearly same
results as the fine mesh using the STVD scheme. It is concluded that the DES model
succeeds where RANS models often fail in predicting the mean flow and global flow
quantities.
Rodi [31] has used RANS and LES to calculate flow over bluff bodies and the results
are compared with the experimental data. The bluff bodies used are a long square
cylinder (Re = 22,000) and a surface mounted cube (Re = 40,000). Calculations were
carried out using k?epsilon1 model. k?epsilon1 model with the Kato Launder [34] correction was
also used. The Kato Lauder correction eliminates spurious turbulence production from
19
the k?epsilon1 model in the stagnation region. A two-layer approach was used in connection
with k?epsilon1 model for both cases. In the two-layer approach [35] the viscous sub-layer is
resolved with a simpler one-equation model. Smagorinsky?s SGS model and Germano?s
dynamic model were used as sub-grid scale models.
For the square cylinder, results provided by the k?epsilon1 model are poor. The Kato-
Launder modification yields improved results. Further improvement in the results is
obtained with the combination of the two-layer approach resolving the near wall region
and the Kato-Launder correction. In all RANS calculations the turbulent fluctuations
are severely under-predicted. It was found that none of the LES results are uniformly
good and entirely satisfactory, and there were large differences between the individual
calculations which are difficult to explain. For the surface mounted cube the k ? epsilon1
model leads to poor prediction over the top. The result is improved with Kato-Launder
correction. A combination of the two layer approach and the Kato-Launder correction
gives better results for the global quantities except the separation length which worsens
when we shift from standard model to the combination. LES gives better results overall
because it simulates all the complex features of the 3-dimensional flow. The difference
between various LES calculations is much smaller for the surface mounted cube when
compared with those from the square cylinder.
20
Hybrid RANS/LES Models
Spalart and Squires [16] have addressed the precautions which should be taken to
obtain better results from simulation of flow over bluff bodies using the DES model.
They are listed below.
1. Very long time samples are essential to converge the statistics, as vortex shedding
has strong modulations in DES. Simulating only a few cycles of such shedding is
unsafe.
2. The authors have mentioned that the error resulting from dissipation from numer-
ical discretization of the Navier-Stokes equations combined with the DES model
can be large and requires great care in controlling it. One of the solution to this
problem is increasing the grid density by refining the grid and varying the time
step involved in the simulations.
3. AprecautionarycheckshouldbeperformedbytheusertoconfirmthattheRANS/LES
interface is not deep inside the boundary layer because it leads to inaccuracies in
the solution.
The authors have also highlighted some areas of research regarding the use of DES
model that need to be investigated.
1. Research on the flexible design of DES simulation which will involve LES treatment
within the boundary layer particularly upstream of the separation line because this
has given better results in some cases [17].
21
2. Another issue is the area in which a shear layer, after separation, needs to generate
?LES content? which it did not possess as a boundary layer upstream. If the grid
is not constructed properly the DES model may not create any LES content. This
issue may be important for the research presented in this thesis because shear layers
separated from the boundary layers formed on the GTS surfaces are released in
the LES region and are not correctly simulated.
3. A source of ambiguity exists in use of DES for vortex-dominated flows because of
the effect of the DES limiter. The risk involved is that the solution may not have
LES behavior in the intended LES region.
Nichols [27] used three hybrid RANS/LES models to simulate the flow over a circular
cylinder (Re = 8 ? 106). The models were:
1. Spalart-Allmaras DES model (SA-DES model)
2. Shear Stress Transport DES model (SST DES model)
3. SST-Multiscale DES model
The SST DES model and SST-Multiscale model transition from RANS to LES as a
function of local turbulent length scale and the local grid spacing unlike the Spalart-
Allmaras DES model. The Spalart-Allmaras DES model transitions from RANS to
LES as a function of local grid spacing only. The SA-DES model did not reach grid
convergence for the cylinder average drag in the refinement study made by the author.
The SST-based models are less sensitive to grid density variations than the SA-DES
22
model. The grid sensitivity of the SA-DES model indicates that a filter function that
includes some measure of the local turbulence length scales along with the local grid scale
is superior to a filter function that includes only the local grid scale. A computational
mesh with a grid resolution resulting in a ratio of the turbulent length scale to grid
length scale greater than two produced a reasonable simulation with all three hybrid
models for the square cylinder case.
Camarri et al. [7] have devised a hybrid RANS/LES model called as the Limited
Numerical Scale (LNS) model. They have applied it to the simulation of the flow around
a square cylinder at a Reynolds number of Re = 22,000. The model uses local blending of
two eddy viscosities. One eddy viscosity comes from Smagorinsky?s sub-grid scale model
in the LES and other comes from the k?epsilon1 RANS model. In LNS, unlike DES, there is
no a-priori separation between regions to be treated with LES and RANS models. Also,
DES is based on the Spalart-Allmaras RANS model while the LNS strategy is directly
applicable to most available RANS models. Simulations carried out on a grid sufficiently
refined for LES show that the LNS results are almost identical to those given by LES.
For a coarser resolution the results from LNS are better than LES and RANS because
of its blending criteria.
23
Chapter 2
Problem Description
2.1 Wind Tunnel Experiment
A 1/8-scale simplified model of the tractor-trailer geometry (or the GTS model) was
studied experimentally in the NASA Ames 7 ft?10 ft wind tunnel. A complete set of the
experimental data and the surface definition of the model are included on a CD-ROM
for further analysis and comparison. The wind tunnel has 4.57 m long test section and a
constant height of 2.13 m. The width of the test section is 3.05 m with 1% wall divergence
[4]. Figure 2.1 shows the GTS geometry used during the experiment. The length of the
model is 2.48 m (x/W = 7.647), width W = 0.324 m and the height is 0.45 m (y/W
= 1.392). This model is without any wheels and there is no tractor-trailer gap, which
simplifies the grid generation. The dimensions of the wind tunnel and the GTS model
were non-dimensionalized by the trailer width W = 32.38 cm. Because of the interest in
wake details, the GTS model was located 13.33 cm downstream of the beginning of the
test section of the wind tunnel and the bottom of the model was located 7.6 cm above
the wind-tunnel floor. A body-axis coordinate system was used in the experiment. Force
and pressure measurements were made at yaw angles (?) ranging from -14? to 14?. Yaw
angle is the angle made by the longitudinal axis of a vehicle with the direction of the air
flow. Detailed data sets were obtained at Reynolds numbers of Re = 2?106 and Re =
7.4?105. These Reynolds numbers are based on the width of the trailer. These numbers
24
Figure 2.1: The GTS model
correspond to wind-tunnel velocities of 205 mph and 75 mph respectively [4]. It should
be noted that the computational simulations are carried out in this research are only for
the case of Reynolds number of Re = 2?106 and a yaw angle of ? = 0?. Table 2.1 shows
the wind tunnel flow conditions for the Re = 2 million experiment.
Table 2.1: Wind tunnel conditions [4].
?, ? 0
Re, million 2.021
M 0.279
V?, m/s 94.48
q, N/m2 5266.82
pt, N/m2 102492.48
ps or pref, N/m2 97216.07
Tt, K 285.76
Ts, or Tref, K 276.59
??, kg/m3 1.180
25
2.1.1 Surface Pressures
79 surface pressure taps were located on the model and 43 additional pressure taps
were located on the right wall of the test section (looking upstream). These pressure
taps were used to study the pressure distribution around the GTS model. In addition,
one unsteady pressure transducer was located on the rear of the trailer (x/W = 7.647,
y/W = 0.63 and z/W = -0.46 ). The unsteady pressure signal was used to obtain
power spectra of the unsteady pressure. The reference static pressure tap was located
at x/W = 4.47, y/W = 2.588 and z/W = -4.7. This location is on the wall of the
wind tunnel. See figures in [4] for the location of the pressure taps. The reference static
pressure is used in calculation of the pressure coefficient and the drag coefficient. The
uncertainty in the pressure coefficient calculated from the experiment was ?0.002. This
uncertainty estimate includes only the measurement precision adjusted for the run-to-
run repeatability [4]. Table 2.2 shows the experimental pressure coefficient distribution
along the y/W axis of the trailer base at three z/W locations. These values are available
in the CD-ROM. The values of these experimental pressure coefficients will be used for
comparison with computational results.
Table 2.2: Experimental values of pressure coefficient at different z/W locations along
the y/W axis of the trailer base [4].
y/W location ?Cp at z/W = 0 ?Cp at z/W = 0.2206 ?Cp at z/W = 0.44120
1.3333 -0.080177 -0.087359 -0.082658
1.0147 -0.101592 -0.107468 -0.104073
0.6961 -0.162834 -0.162051 -0.153302
0.3775 -0.202139 -0.198483 -0.187775
0.0588 -0.155130 -0.152519 -0.143639
26
2.1.2 Drag, Lift and Side Force Measurements
The drag measurements are accurate to within ?1 lbf (4.45 N). The repeatability
of the drag coefficient measurements was ?0.001 for ? ? -5? and ?0.01 for ? < -5?.
These error bands include measurement resolution and point-to-point repeatability [4].
The accuracy of the side and lift forces is not mentioned. The experimental values of
drag, lift and side force coefficient will be compared with simulation results in Section 5.
2.1.3 Oil-film Interferometry (OFI), Hot-film Anemometry and Particle Im-
age Velocimetry (PIV)
Oil-film interferometry was done on the top and right side of the GTS model to
measure the skin friction on those surfaces. The uncertainty in this technique is ?5%.
Hot-film measurements were made on the right side of the cab to determine the existence
and extent of flow seperation in the vicinity of the front corner radius. Three-component
PIV was conducted in several planes in the wake of the trailer to measure the velocities
in those planes [4].
27
Chapter 3
Simulation Approach
3.1 Governing Equations
The Favre-averaged Navier-Stokes equations (equations 1.2 to 1.4) are solved in the
simulations. The Mach number at which the experiment was performed was M = 0.279
and the solver used for the simulations was a compressible fluid flow solver. Equations
1.2 to 1.4 are in differential form. It should be noted that the differential and the
integral forms of the Navier-Stokes equations are equivalent, and the CFD solver used
in this research uses integral form of the Navier-Stokes equations (i.e., the finite-volume
method).
3.2 The DES Model
The DES model used in this work is a hybrid RANS/LES model proposed by Spalart
and co-workers [8]. It uses the Spalart-Allmaras one-equation model for modeling the
attached boundary layer region and LES for the simulations of highly separated flow
regions. The Spalart-Allmaras model [6] uses one partial differential equation for calcu-
lation of the working variable ??
???
?t +vi
???
?si = cb1
?S???cw1fw
parenleftbigg??
d
parenrightbigg2
+ 1? ??s
i
parenleftbigg
(? + ??) ????s
i
parenrightbigg
+ cb2? ????s
i
???
?si (3.1)
28
The working variable is used for the calculation of the eddy viscosity
?t = ??f?1 (3.2)
Other than the main partial differential equation used for calculating the working variable
the model uses following auxiliary equations
fv1 = ?
3
?3 +cv13 (3.3)
fv2 = 1? ?1+?f
v1
(3.4)
fw = g
parenleftbigg 1+c
w36
g6 +cw36
parenrightbigg1/6
(3.5)
? = ??? (3.6)
g = r +cw2parenleftbigr6 ?rparenrightbig (3.7)
r = ???S?2d2 (3.8)
29
?S = S +fv2 ??
?2d2 (3.9)
S = radicalbig2?ij?ij (3.10)
The tensor ?ij is the rotation tensor and is given by
?ij = 12
parenleftbigg?v
i
?xj ?
?vj
?xi
parenrightbigg
(3.11)
The closure coefficients [6] are given by
cb1 = 0.1355
cb2 = 0.622
cv1 = 7.1
? = 23
cw1 = cb1?2 + (1+cb2)?
cw2 = 0.3
cw3 = 2
? = 0.41 (3.12)
The transition from RANS modeling in the boundary layer to LES outside the boundary
layer is performed by changing the definition of the distance d of a point to the nearest
30
wall in the equation 3.9 [8]. The distance d is replaced with ?d.
?d = min(d,CDES?) (3.13)
The constant CDES is given by Spalart [8] and is equal to 0.65. Far from the wall, the
value of ?d becomes
?d = CDES? (3.14)
where ? is the local grid spacing and is equal to maximum mesh spacing in the three
coordinate directions
? = max(?x,?y,?z) (3.15)
When the production term is balanced with the dissipation term we get the dynamic
eddy viscosity
?t = ?(CSm?)2radicalbigSijSij (3.16)
and the value of the constant CSm
CSm = CDES
radicalBigg
cb1fv1
cw1fw (3.17)
31
where Sij is the strain rate tensor
Sij = 12
parenleftbigg??v
i
?xj +
? ?vj
?xi
parenrightbigg
(3.18)
In the outer part of the boundary layer CSm asymptotes to CSm = 0.29CDES = 0.19.
The DES model thus asymptotes to a Smagorinsky-type LES model in the bluff-body
wake assuming sufficient mesh refinement [8].
3.3 The SACCARA Code
The SACCARA (Sandia Advanced Code for Compressible Aerothermodynamics
Research and Analysis) code is a compressible fluid flow solver at the Sandia Laboratories
and was used for the DES computations performed during this research. It was developed
from the INCA code [37] written by Amtec Engineering. It is used to solve the Navier-
Stokes equations and the turbulence transport equations in conservation form. Code
verification has been performed on the SACCARA code by code-to-code comparison
with other Navier-Stokes codes [38] and with Direct Simulation Monte Carlo method [39].
The code uses a Lower-Upper Symmetric Gauss-Seidel scheme to solve the equations.
This scheme is based on the works of Yoon et al. [44] and Peery and Imlay [46] and has
excellent scalability up to thousands of processors.
32
3.4 Discretization of the Convective Terms
3.4.1 Finite Volume Method
The discretization approach used for convective terms is the finite volume method
(FVM). The FVM relies on obtaining the governing equations in integral form. It can
also be used for unstructured grids unlike the finite difference methods. The SACCARA
code is based on cell-centered finite-volume approach. The cell-centered finite-volume
approach involves evaluation of exact integrals of properties over a cell by approximating
them by the values at the center of the cell.
3.4.2 Upwind Schemes
An upwind scheme was used for the discretization of the convective terms. In the
upwind scheme the convective term can be expressed as a central difference term plus a
numerical diffusion term. These methods provide a physically-based method for intro-
ducing numerical diffusion. They can capture sharp gradients/discontinuities and they
minimize dispersion. Three types of upwind schemes are the flux vector splitting (FVS)
approach, the flux difference splitting (FDS) approach and hybrid upwind schemes. The
SACCARA code has available both FVS and FDS schemes, and a FDS scheme was used
for the DES.
The current simulations employ Yee?s symmetric total variation diminishing (TVD)
scheme. The scheme is modified to incorporate Harten?s artificial compression switch
(ACM) [33] to further reduce the dissipation, along with a characteristic based filtering
33
method [32]. Yee?s scheme uses the Roe interface flux, which may be written as
Zm+1/2 = Zm +Zm+12 + 12R1m+1/2?m+1/2 (3.19)
Here 12 [Zm +Zm+1] is the central differencing portion of the numerical flux Zm+1/2 and
R1m+1/2?m+1/2 is the nonlinear dissipation term. R1m+1/2 is the right eigenvector of
the matrix ?Zm+1/2/?Uj+1/2 where U is the conservation flow vector. The non-linear
dissipation term is combined with Harten?s switch and this switch is used to decide the
amount of non-linear dissipation in a non-dissipative scheme. The modified dissipation
term looks like
Z?m+1/2 = 12R1m+1/2??m+1/2 (3.20)
where the elements of ??m+1/2 are given by
?l?m+1/2 = K?lm+1/2?lm+1/2 (3.21)
The?lm+1/2 are the elements of the vector ?m+1/2 in the symmetric TVD scheme and
are given by
?lm+1/2 =?
vextendsinglevextendsingle
vextendsingle?lm+1/2
vextendsinglevextendsingle
vextendsingle
parenleftBig
?lm+1/2 ?Qlm+1/2
parenrightBig
(3.22)
34
where ?lm+1/2 is the eigenvalue of ?Zm+1/2/?Uj+1/2 and Qlm+1/2 is the minmod limiter.
The dissipation can be scaled down globally by reducing the value of the numerical
dissipation constant K. The value of K varies between 0 to 1. ?lm+1/2 serves as Harten?s
switch and it can be used to further scale down the numerical dissipation in regions of
smooth flow
?lm+1/2 =
vextendsinglevextendsingle
vextendsingle?lm+1/2 ??lm?1/2
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle?l
m+1/2 +?
l
m?1/2
vextendsinglevextendsingle
vextendsingle
(3.23)
where ?lm+1/2 are the elements of R1?1j+1/2 (Um+1 ?Um).
3.5 Discretization of the Diffusion Terms
The viscous terms were discretized using central differencing.
3.6 Temporal Discretization
Sub-iterative procedure is used to obtain second-order accuracy during temporal
discretization. In this the summation of the discretized temporal derivative and the
steady-state residual at the n + 1 time level is iterated until it reaches the specified
tolerance value.
3.7 Grid Generation Approach
The grid generation approach used during this research was the structure grid ap-
proach using the commercial grid generation software Gridgen. In this approach the grid
35
is laid in a regular repeating pattern called block. For a 3-dimensional mesh the struc-
tured grid uses hexahedral elements in a computationally rectangular array. The grid
can be stretched and twisted to fit it according to the body shape. The multi-block grid
generation approach which is used in this research allows several blocks to be connected
together to construct the whole domain. The requirements of block to block connectiv-
ity can vary according to the software used. The SACCARA code uses point to point
connectivity where the blocks must match topologically and physically at the boundary.
Structured grids provide a very high degree of control over the grid but the time
and expertise required to create a structured grid for the entire model is a challenging
task [49].
3.8 Domain Decomposition and Load Balancing
In dealing with geometrically large complicated systems the mesh domain can be
decomposed into subdomains and each subdomain solved on a separate processor of
a cluster or a supercomputer. This method is called domain decomposition method
(DDM). The main advantages of DDM include efficiency of solvers, savings in compu-
tational storage conducive to parallel processing and saving in computational time [45].
In a DDM method the domain ?(t) is expressed as a union of subdomains
?(t) =
np(t)uniondisplay
c=1
?c(t) (3.24)
36
An important consideration while solving large grids on distributed memory archi-
tectures is the distribution of the mesh across the memory of the machine at runtime
so that the calculated load is evenly balanced and the amount of interprocessor com-
munication is minimized [45]. For the present case the coarse and medium meshes were
sub-divided into 98 sub-domains, and these sub-domains were distributed over 86 proces-
sors for both grids. Table 3.1 shows the number of cells in coarse and medium mesh the
number processors used for the simulations and the average number of cells/processor.
Table 3.1: Coarse and medium mesh
Mesh Number of cells np average number of cells/processor
Coarse 3.9 ? 106 86 45,349
Medium 1.32 ? 107 86 153,488
3.9 Boundary Conditions used for the Simulations
The inflow boundary employs stagnation values for pressure pt = 102,653 N/m2
and temperature Tt = 282.1 K and enforces inflow normal to the boundary. The outflow
boundary used a fixed static pressure of ps = 97,700 N/m2. The back pressure was chosen
so that the tunnel wall reference pressure located at x/W = 4.47, y/W = 2.588 and z/W
= -4.7 matched with the experiment. The reference pressure from the experiment and the
average reference pressure from the coarse mesh simulations are compared in figure 3.1.
The experimental reference pressure is 97,336 N/m2 and the average reference pressure
from the coarse mesh simulations is 97,295 N/m2. It can be seen in the figure that the
reference pressure from the simulations matches pretty well with reference pressure from
37
the experiment. Table 3.2 compares average reference pressure values from the coarse
and medium mesh with the experimental value. Slip conditions on velocity are employed
on the top and side walls of the wind tunnel to save the computational cost, while the
floor of the wind tunnel, the GTS surface and the support posts employ no-slip velocity
conditions and assume an adiabatic wall. The freestream dynamic eddy viscosity is set
at ?t = 1?10?5 Ns/m2. Solid wall boundary condition for the turbulence model can be
found in [47].
Table 3.2: Reference pressure comparison
Experiment Coarse Mesh Medium Mesh
97,336 N/m2 97,295 N/m2 97,249 N/m2
3.10 Characteristic Scales
Time has been non-dimensionalized by a reference time scale defined as
? = W/V? = 0.0034s (3.25)
Where W is the trailer width and V is the freestream velocity. The value of V? was
calculated from the reference static pressure and stagnation pressure and it came out
to be 93.8 m/s1. Characteristic length scales for the medium and coarse meshes are
given in Table 3.3 (Note that the length scales are non-dimensionalized by the trailer
width). Extremely fine grid spacing is required near the wall because the RANS model
is integrated to the wall. The maximum y+ value near the wall from the coarse mesh
38
is 1.4. The maximum grid spacing in the trailer wake leads to approximately 26 points
across the trailer width for the coarse grid case and 38 points for the medium grid case.
Figures 3.2 and 3.3 show the coarse and the medium mesh structure near the truck base.
The medium mesh has higher grid density than the coarse mesh.
Table 3.3: Characteristics scales
Mesh ?wall ?wake
Coarse 6.17?10?5 0.12
Medium 3.08?10?5 0.081
39
?
p(
N/
m2
)
0 200 400 60095000
96000
97000
98000
99000
100000
Coarse Mesh - Instantaneous
Coarse Mesh - Mean
Experiment
Reference Static Pressure from the Coarse Mesh
and the Experiment
Figure 3.1: Reference static pressure comparison between the coarse mesh simulations
and the experiment
40
x/W
y/W
7 8 9 10 11
0
1
2
3
4
Coarse Mesh near the Truck Base,
z/W = 0.0
Truck
base
Figure 3.2: Coarse mesh structure near the truck base
41
x/W
y/W
7 8 9 10 11
0
1
2
3
4
Medium Mesh near the Truck Base,
z/W = 0.0 plane
Truck
base
Figure 3.3: Medium mesh structure near the truck base
42
Chapter 4
Numerical Accuracy
Numerical accuracy is concerned with difference between numerical solution and
exact solution to a differential equation. The Navier-Stokes equations for general com-
pressible flows do not have an exact solution. In this Chapter we will assess the numerical
accuracy of the solution by comparing it with the best available solution. In Section 4.1
an effort is made to predict the length of time window required for obtaining statistically
converged results. In Section 4.2, the time period of existence of initial transients in the
solution is estimated. In Section 4.3, the results are analyzed to find out the number of
sub-iterations required for obtaining numerically accurate results. In Section 4.4 effect
of the temporal discretization on the results is studied. In Section 4.5, the effect of the
mesh size on the solution is studied. It should be noted that discretization error cannot
be fully assessed in this case because PDEs solved in the DES model change with a
change in the mesh size. Hence the whole model is changed when it is applied from the
coarse to the medium mesh.
4.1 Statistical Convergence
DES is time dependent and to assess mean properties from these simulations we
need to average the data from the simulations over some time window. The required
size of the time window is determined by performing a statistical convergence analysis.
43
In this research the pressure and velocities in the base of the truck were used to perform
this analysis. The reason for selecting these variables near the base region was that this
is the unsteady flow region where the eddies are formed and shed. The flow over the
front, top, side and bottom of the GTS model is mostly steady flow without any massive
separation. Hence this flow is modeled using RANS. The pressure and velocity data
were time averaged over a successively increasing time window. When the profiles from
two successive time windows matched each other within a specific tolerance then the
smaller time window out of the two windows is said to be the time window over which
the data should be averaged to get statistically converged results. Figure 4.1 shows the
unsteady pressure signal in the base of the truck (x/W = 7.647, y/W = 0.63 and z/W
= -0.46) plotted against non-dimensional time ?. The width of the rectangular windows
correspond to the successively increasing time windows used for averaging the pressure
and velocities. Note that the averaging starts in the reverse direction along the time-
axis starting from ? = 571. This is done in order to stay away from the initial transient
period. In figure 4.2 the line T which is parallel to the y/W axis shows the location
along which the pressure profiles are time-averaged. Line A, which is parallel to x/W
axis, shows the location along which the velocity profiles are time averaged to obtain
statistical convergence.
Figure 4.3 shows the the pressure coefficient (Cp) distribution plotted along the line
T and corresponding to the various time windows shown in figure 4.1. In the figure
it can be seen that only the Cp profile corresponding to ?start = 527 or the profile
44
?
p(
N/
m2
)
0 200 400 60095000
95500
96000
96500
97000
97500
98000
98500
99000
99500
100000
?start = 301, ?end = 571
?start = 347, ?end = 571
?start = 394, ?end = 571
?start = 439, ?end = 571
?start = 483, ?end = 571
?start = 527, ?end = 571
Calculation of time window for
statistical convergence
Figure 4.1: Time windows used for calculating statistically converged results
45
X
Y
Z
A
T
Figure 4.2: Location of pressure and velocity profiles used for statistical convergence
corresponding to smallest time window is not statistically converged. All other profiles
seem to be statistically converged. Figure 4.4 shows the u/Vref, v/Vref and urms/Vref
velocity profiles plotted along the line A and corresponding to the various time windows
shown in figure 4.1. In this figure a look at the v/Vref and u/Vref profiles suggests
that all the profiles except the one corresponding to ?start = 527 seem to to statistically
converged. But the u/Vref profiles corresponding to the three biggest time windows are
different from three profiles corresponding to the three smaller windows at some places.
Figure 4.5 shows the absolute errors in all the Cp profiles corresponding to the various
time windows shown in figure 4.1 with respect to the Cp profile corresponding to the
biggest time window in figure 4.1. The biggest time window gives the best statistically
converged results because it is averaged over longest period of time. Figures 4.6, 4.7
and 4.8 show the absolute errors in the u/Vref, v/Vref and urms/Vref velocity profiles,
46
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.3: Pressure distribution along line T and corresponding to different time win-
dows
47
x/W
u/V
re
f,
v/V
re
f,
u rm
s/
V r
ef
7.5 8 8.5 9 9.5 10-0.4
-0.2
0
0.2
0.4
0.6
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3 , nIt = 4
urms/Vref
v/Vref
u/Vref
Figure 4.4: Velocity profiles along line A and corresponding to different time windows
48
(Cp)error
y/W
-0.01 -0.005 0 0.005
0
0.2
0.4
0.6
0.8
1
1.2
1.4
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.5: Error in the pressure distri-
bution obtained from smaller time win-
dows with respect to the pressure distri-
bution obtained from the biggest time
window
x/W
(u
/V
re
f) e
rro
r
8 8.5 9 9.5-0.01
0
0.01
0.02
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.6: Error in the u/Vref veloc-
ity profiles obtained from smaller time
windows with respect to the u/Vref ve-
locity profile obtained from biggest time
window
x/W
(v/
V r
ef)
er
ro
r
8 8.5 9 9.5-0.04
-0.03
-0.02
-0.01
0
0.01
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.7: Error in the v/Vref veloc-
ity profiles obtained from smaller time
windows with respect to the v/Vref ve-
locity profile obtained from biggest time
window
x/W
(u
rm
s/
V r
ef)
er
ro
r
8 8.5 9 9.5
-0.02
-0.01
0
0.01
0.02
0.03
0.04
?start = 527, ?end = 571
?start = 483, ?end = 571
?start = 439, ?end = 571
?start = 394, ?end = 571
?start = 347, ?end = 571
?start = 301, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.8: Error in the urms/Vref
velocity profiles obtained from smaller
time windows with respect to the
urms/Vref velocity profile obtained from
biggest time window
49
respectively, corresponding to the various time windows shown in figure 4.1 with respect
to the biggest time window. Table 4.1 shows the largest absolute errors with respect
to the biggest time window for various profiles seen in figures 4.5, 4.6, 4.7 and 4.8 (the
relative percentage error is not calculated because the denominator in some cases may be
zero). From the plots it is concluded that the smallest time window showing statistically
converged results is the window with ??w = 177 because it has reasonably small errors
with respect to the biggest time window. Other observation made from the various
plots and the table is that the Cp profiles converges the fastest and the v/Vref profile
tends to converge the slowest. Although this analysis is based on coarse mesh results, it
is assumed to hold true for the medium mesh results as well. Henceforth, any analysis
done with the coarse and medium mesh is by averaging the data over time ??w = 177 for
achieving statistically converged results. This non-dimensional time window corresponds
to a physical time of ?tw = 0.611 s.
Table 4.1: Absolute error with respect to biggest time window
(Cp)error (u/Vref)error (v/Vref)error (urms/Vref)error ??w
0.01041 0.02694 0.03916 0.04119 44
0.00448 0.0222 0.02063 0.019286 88
0.00488 0.02457 0.01425 0.01342 132
0.00470 0.00665 0.01103 0.01025 177
0.00359 0.00303 0.01015 0.00557 224
0 0 0 0 270
50
4.2 Initial Transients
Since the simulations are initialized with freestream values for velocity pressure and
temperature the flow over the GTS during some initial time period oscillates because of
the waves generated due to sudden exposure of the truck to the freestream flow. Figure
4.9 shows the plot of the unsteady pressure signal from the truck base region plotted
against time. The large fluctuations in the pressure signal near ? = 0 denotes the period
of initial transients. As the simulations are reduced in time the initial transients die
out. The data collected during the initial transient period should be neglected, thus
the time period of existence of initial transients needs to be calculated. Figure 4.9
shows the time window of width ??w ? 177 (this width was obtained from statistical
convergence analysis in section 4.1) moved along the time-axis. Figures 4.10 and 4.11
show the time-averaged pressure and velocity profiles corresponding to the time windows
shown in figure 4.9. The Cp profiles seen in 4.10 show that the profiles corresponding
to first four time windows (starting from ? = 0) are in the initial transience period.
The last three pressure profiles seem to overlap each other and it can be said that initial
transients have minimized in this period. Almost the same conclusion can be drawn after
observing the time averaged u/Vref, v/Vref and urms/Vref velocity profiles in figure 4.11.
Figure 4.12 shows the absolute errors in the Cp profiles corresponding to various time
windows shown in figure 4.9 with respect to the Cp profile corresponding to the time
window starting at ?start = 394 (or the final time window). The final time window
is furthest along the time-axis from the initial transience period hence it is considered
51
?
p(
N/
m2
)
0 200 400 60095000
95500
96000
96500
97000
97500
98000
98500
99000
99500
100000
?start = 394, ?end = 571
?start = 348, ?end = 527
?start = 301, ?end = 483
?start = 232, ?end = 417
?start = 162, ?end = 348
?start = 95, ?end = 278
?start = 23, ?end = 208
Estimation of initial transience period
Figure 4.9: Time windows used for estimating period of initial transients
52
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
?start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.10: Pressure distribution along line T and corresponding to different time win-
dows
53
x/W
u/V
re
f,
v/V
re
f,
u rm
s/
V r
ef
7.5 8 8.5 9 9.5 10-0.4
-0.2
0
0.2
0.4
0.6
?start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
urms/Vref
v/Vref
u/Vref
Figure 4.11: Velocity profiles along line A and corresponding to different time windows
54
(Cp)error
y/W
-0.005 0 0.005 0.01 0.0150
0.2
0.4
0.6
0.8
1
1.2
?start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.12: Error in the pressure dis-
tribution from all time windows with
respect to the pressure distribution ob-
tained from the window farthest along
time-axis
x/W
(u
/V
re
f) e
rro
r
8 8.5 9 9.5-0.1
-0.05
0
0.05
?start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.13: Error in the u/Vref pro-
files from all time windows with respect
to the u/Vref profile obtained from the
window farthest along time-axis
the ideal window. Figures 4.13 to 4.15 show the absolute errors in the u/Vref, v/Vref
and urms/Vref profiles corresponding to all the time windows shown in figure 4.9 with
respect to the final time window. It is very obvious that the pressure and velocity profiles
corresponding to window starting at ?start = 23 will show maximum errors because of
the initial transients present in the results. Table 4.2 shows the biggest absolute errors
in the pressure and velocity profiles from all the windows with respect to pressure and
velocity profiles from the window corresponding to ?start = 394 (or the final window). It
is observed from the various plots and the table that the Cp profile converges the fastest
and the urms/Vref profile converges the slowest. This analysis helps in estimating the
period of initial transients. Examination of the pressure and velocity errors from table
4.2 suggests that the errors corresponding to windows with ?start = 301 and ?start = 348
55
x/W
(v/
V r
ef)
er
ro
r
8 8.5 9 9.5-0.04
-0.02
0
0.02
0.04 ?
start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, Kappa = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.14: Error in the v/Vref pro-
files from all time windows with respect
to the v/Vref profile obtained from the
window farthest along time-axis
x/W
(u
rm
s/
V r
ef)
er
ro
r
8 8.5 9 9.5
-0.02
-0.01
0
0.01
0.02
?start = 23, ?end = 208
?start = 95, ?end = 278
?start = 162, ?end = 348
?start = 232, ?end = 417
?start = 301, ?end = 483
?start = 348, ?end = 527
?start = 394, ?end = 571
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3, nIt = 4
Figure 4.15: Error in the urms/Vref pro-
files from all time windows with respect
to the urms/Vref profile obtained from
the window farthest along time-axis
are acceptable and match closely with the profiles corresponding to the window with
?start = 394 (or the final window). Hence we conclude that the initial transients period
exists from ? = 0 to ? = 301. Therefore any data in the region from ? = 0 to ? = 301
will be neglected. This corresponds to a physical time period of ?tw ? 1 s. All the
comparisons with the experimental data will be made by using the data obtained after
? = 301.
Table 4.2: Absolute error with respect to time window farthest along the time-axis
Time window (Cp)error (u/Vref)error (v/Vref)error (urms/Vref)error
?start = 23, ?end = 208 0.01951 0.09295 0.04087 0.01897
?start = 95, ?end = 278 0.01743 0.07368 0.03312 0.02639
?start = 162, ?end = 348 0.00822 0.05606 0.03014 0.02241
?start = 232, ?end = 417 0.00902 0.04203 0.02519 0.02383
?start = 301, ?end = 483 0.00246 0.01557 0.01524 0.01934
?start = 348, ?end = 527 0.00664 0.00850 0.00992 0.01533
?start = 394, ?end = 571 0 0 0 0
56
4.3 Effect of Number of Sub-iterations used for Iterative Convergence
The number of sub-iterations to be performed for achieving iterative convergence
during each time step is not exactly known for this complicated turbulent flow. The time
steps used for the simulations are small (??s1 = 1.162?10?3 and ??s2 = 2.325?10?3)
and for such a small time steps, four to six sub-iterations performed for each time-step
is believed to give an iteratively converged solution. The effect of the number of sub-
iterations performed during each time-step was studied by performing one, four and eight
sub-iterations per time step and the pressure and velocity profiles from these cases are
compared in figures 4.16 and 4.17. The comparison of the pressure profile shows that
the profiles corresponding to nIt = 4 and nIt = 1 are equidistant from the nIt2 profile,
which should be the most accurate profile because of the larger number of sub-iterations.
Hence no conclusion can be made from the pressure plot. The effect of number of sub-
iterations is clearly visible while comparing the u/Vref velocity profiles. Comparison of
the velocity profiles infigure 4.17 shows that the u/Vref profile corresponding to nIt =
1 shows the biggest absolute error of approximately 0.03 with respect to the nIt = 8
case. The u/Vref velocity profile corresponding to nIt = 4 shows much smaller error
with respect to the nIt = 8 case. It is concluded from the velocity profile plot that four
sub-iterations are sufficient for achieving iterative convergence when the time step of
??s1 = 1.162 ? 10?3 is used in the simulations.
57
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
?start = 406, ?end = 583, nIt = 1
?start = 394, ?end = 571, nIt = 4
?start = 406, ?end = 583, nIt = 8
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3
Figure 4.16: Pressure distribution along line T and corresponding to different number
of sub-iterations
58
x/W
u/V
re
f,
v/V
re
f,
u rm
s/
V r
ef
7.5 8 8.5 9 9.5 10-0.4
-0.2
0
0.2
0.4
0.6
?start = 406, ?end = 583, nIt = 1
?start = 394, ?end = 571, nIt = 4
?start = 406, ?end = 583, nIt = 8
DES, Coarse Mesh, K = 0.2, ??s = 1.16e-3
urms/Vref
v/Vref
u/Vref
Figure 4.17: Velocity profiles along line A and corresponding to different number of
sub-iterations
59
4.4 Effect of Time Step (??s)
The effect of time step (??s) used for discretizing the Navier-Stokes equations was
analyzed by using two time steps for the coarse mesh. As mentioned previously the
time-steps used were ??s1 = 1.162 ? 10?3 (?ts1 = 4 ? 10?6 s) and ??s2 = 2.325 ?
10?2 (?ts2 = 8 ? 10?6 s). Figure 4.18 shows the pressure distribution at the z/W = 0
along the vertical axis of the truck base (Line T in Figure 4.2) obtained from the two
time-steps used. Figure 4.19 shows the velocity profiles (along line A in Figure 4.2) of the
u, v and urms velocity profile obtained from the two time steps used for discretization.
There is a significant difference in the pressure and velocity profiles obtained from the
two time steps. One has to be cautious while comparing the profiles from these two time
steps. The number of sub-iterations used for iterative convergence at the time step with
??s1 case is 4 and for the ??s2 case is 6. As mentioned in section 4.3 the number of sub-
iterations that need to be performed at each time step to guarantee iterative convergence
of the solution is not known. So even after performing 6 sub-iterations for the ??s2 case
it is not know that the solution iteratively converged to the same level as in the ??s1
case with 4 sub-iterations. So it cannot be concluded whether the difference between
the pressure and velocity profiles corresponding to the two time steps used is a result of
temporal discretization error or an error due to iterative convergence. We can conclude
that the difference between the pressure and velocity profiles from the two time-steps
used for the simulations is a combination of the error due to iterative convergence, error
due temporal discretization and presence of some transients in the pressure and velocity
60
profiles for the ??s2 case. The presence of transients in the ??s2 case comes into the
picture because this case was obtained by changing the time step from ??s1 to ??s2 at
time ? = 394 during the simulations.
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
?start = 394, ?end = 571
?start = 422, ?end = 599
DES, Coarse Mesh, K = 0.2, nIt = 4
??s1 = 1.16e-3, ??s2 = 2.32e-3
??s1 = 2.32e-3
??s2 = 1.16e-3
Figure 4.18: Pressure distribution along line T and corresponding to different time-steps
61
x/W
u/V
re
f,
v/V
re
f,
u rm
s/
V r
ef
7.5 8 8.5 9 9.5 10-0.4
-0.2
0
0.2
0.4
0.6
?start = 394, ?end = 571 (??s1 = 1.16e-3)
?start = 422, ?end = 599 (??s2 = 2.32e-3)
DES, Coarse Mesh, K = 0.2, nIt = 4
??s1 = 1.16e-3, ??s2 = 2.32e-3
urms/Vref
v/Vref
u/Vref
Figure 4.19: Velocity profiles along line A and corresponding to different time-steps
62
4.5 Effect of the Mesh Size on DES
The results from the medium and coarse mesh are compared in this section to ex-
amine the effect of the mesh size on the results from the simulations. Figure 4.20 and
4.21 show the u/Vref velocity contours and streamlines in the vertical streamwise plane
(red colored plane in figure 5.1) from the medium and the coarse mesh simulations. The
upper vortex from the coarse mesh simulations is much more dominant than the lower
vortex, while the vortices are almost equal in size for the medium mesh results. Figure
4.22 compares the pressure distribution obtained from the two grids. The pressure distri-
bution obtained from the medium mesh is much more flat than the pressure distribution
from the coarse mesh. The reason for this is the almost equal dominance of the upper
and lower vortices in the medium mesh results in figure 4.20 and the larger size of the
upper vortex in the coarse mesh simulations in figure 4.21. Figure 4.23 compares the
velocity profiles from the two grids. The u/Vref and urms/Vref velocity profiles from
the grids match to each other but the v/Vref profiles from the two grids are different.
The reason for this can again be attributed to the difference in size of the upper vortex
obtained from the coarse and medium grid results.
63
x/W
y/W
8 9 10
0
0.5
1
1.5
2
2.5
u/Vref
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
DES, Medium Mesh, nIt = 5
K = 0.2, ??s = 2.32e-3, ?start = 304, ?end = 481
Figure 4.20: Medium mesh vertical stream-
wise contour plot
x/W
y/W
7.5 8 8.5 9 9.5 10
0
0.5
1
1.5
2
2.5
u/Vref
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
DES, Coarse Mesh, nIt = 4
K = 0.2, ??s = 1.16e-3, ?start = 301, ?end = 479
Figure 4.21: Coarse mesh vertical stream-
wise contour plot
64
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Coarse Mesh, ?start = 301, ?end = 479
Medium Mesh, ?start = 304, ?end = 481
DES, Medium and Coarse Mesh, K = 0.2, nIt medium = 5,
??s medium = 2.32e-3, nIt coarse = 4, ??s coarse = 1.16e-3
Figure 4.22: Pressure distribution along line T obtained from the coarse and medium
mesh simulations
65
x/W
u/V
re
f,
v/V
re
f,
u rm
s/
V r
ef
7.5 8 8.5 9 9.5 10-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Coarse Mesh, ?start = 301, ?end = 479
Medium Mesh, ?start = 304, ?end = 481
DES, Medium and Coarse Mesh, K = 0.2, nIt medium = 5,
??s medium = 2.32e-3, nIt coarse = 4, ??s coarse = 1.16e-3
urms/Vref
v/Vref
u/Vref
Figure 4.23: Velocity profiles along line A obtained from the coarse and medium mesh
simulations
66
Chapter 5
Results and Discussion
The results from the simulations will be compared with the experimental data in
this chapter. The reason for this being that the medium mesh results are closer to the
experiments than the coarse mesh results.
5.1 Contour and Streamline Plots
Figure 5.1 shows the various planes in the wake of the GTS model from which the
velocity data from the experiments and the simulations are compared. Figures 5.2 to 5.7
X
Y
Z
Figure 5.1: Planes used for comparison of velocity data
show a comparison of contour and streamline plots from experiment and medium mesh
simulations. The rectangular window shown in the contour plots from the simulations is
67
x/W
y/W
7.5 8 8.5 9 9.5 10
0
0.5
1
1.5
2
2.5
u/Vref
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
Experimental PIV window, u/Vref velocity contours
and streamlines, vertical streamwise cut: z/W = 0.0
Figure 5.2: Vertical streamwise PIV
window in the experiment
x/W
y/W
7.5 8 8.5 9 9.5 10
0
0.5
1
1.5
2
2.5
u/Vref
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
DES, Medium Mesh, u/Vref velocity contours
and streamlines, vertical streamwise cut: z/W = 0.0
Figure 5.3: Vertical streamwise contour
plot from the medium mesh simulations
the window corresponding to the PIV data plane in the experiment. All the data planes
shown in the figures are near the base of the GTS model. The gray colored rectangular
block seen in all the contour plots is the base of the GTS model. Figure 5.2 shows the
vertical streamwise PIV data plane (the red colored plane in the figure 5.1) from the
experiment. Figure 5.3 shows the corresponding contour plot from the simulations of
the medium mesh. Theses two figures show the u/Vref velocity contours and streamlines
in the vertical streamwise plane at z/W = 0.0. The lower vortex is located at x/W
? 8.0 and y/W ? 0.38 in the experiment and at x/W ? 8.6 and y/W ? 0.2 in the
simulations. Thus the experimental data shows the lower vortex located much closer to
the base of the GTS model and at a higher location than the simulations. The upper
vortex is not seen in the experimental PIV window but the nature of the streamlines in
the window suggests that it may be located near the upper right corner of the window.
68
The distribution of the vortices in the medium mesh simulations is much more symmetric
as seen in figure 5.3. This difference in location of the lower and the upper vortices from
the experiment and the simulations gives a major difference in the pressure distribution
along the base of the GTS model.
x/W
z/W
8 8.5 9 9.5 10-1
-0.5
0
0.5
1
v/Vref
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.05
-0.06
-0.08
-0.1
-0.14
Experimental PIV window, v/Vref velocity contours and
streamlines, horizontal streamwise cut: y/W = 0.696
Figure 5.4: Horizontal streamwise PIV
window in the experiment (y/W =
0.696)
x/W
z/W
8 8.5 9 9.5 10-1
-0.5
0
0.5
1
v/Vref
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
-0.1
-0.12
-0.14
DES, Medium Mesh, v/Vref velocity contours and
streamlines, horizontal streamwise cut: y/W= 0.71
Figure 5.5: Horizontal streamwise con-
tour plot from the medium mesh simu-
lations (y/W = 0.709)
Figure 5.4 shows a horizontal streamwise data plane (horizontal violet colored plane
in figure 5.1) from the experiment at y/W = 0.696 showing the v/Vref velocity contours
and streamlines. The corresponding contour plot from the simulations is shown in the
figure 5.5. The v/Vref velocity near the base of the GTS model (area near x/W = 7.7 and
y/W = 0) is going down in the experimental PIV window whereas the simulations show
an opposite trend. As we move away from the GTS model (area near x/W = 8.5 and z/W
= 0) the PIV window shows the v/Vref velocity coming up but the simulations show an
opposite trend. Thus the vertical flow direction from the experiment and the simulations
69
does not match in this horizontal streamwise plane. This is due to the difference in
locations of the upper and lower vortices in the experiment and the simulations as seen
in figures 5.2 and 5.3.
x/W
z/W
8 8.5 9 9.5 10-1
-0.5
0
0.5
1
v/Vref
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
-0.1
-0.12
-0.14
Experimental PIV window, v/Vref velocity contours
and streamlines, vertical streamwise cut: z/W = 1.044
Figure 5.6: Horizontal streamwise PIV
window in the experiment (y/W =
1.044)
x/W
z/W
8 8.5 9 9.5 10-1
-0.5
0
0.5
1
v/Vref
0.1
0.08
0.06
0.04
0.02
0
-0.02
-0.04
-0.06
-0.08
-0.1
-0.12
-0.14
DES, Medium Mesh, v/Vref velocity contours and
streamlines, horizontal streamwise cut: y/W = 1.05
Figure 5.7: Horizontal streamwise con-
tour plot from the medium mesh simu-
lations (y/W = 1.05)
Figure 5.6 shows the v/Vref velocity contours and streamlines in a horizontal stream-
wise plane at y/W = 1.04 (horizontal blue colored plane in figure 5.1) from the experi-
ment and figure 5.7 shows the corresponding plot from the simulations. Even in this case
there is a mismatch of the v/Vref velocity flow directions near and away from the GTS
model when a comparison is made between the experimental plot and the simulations.
70
5.2 Velocity Profiles
Figures 5.8 to 5.16 show a comparison of the various velocity profiles at different
locations in the wake between the experiment and the medium mesh simulations. Un-
certainty bars have been placed on the profiles from the simulations. These uncertainty
bars take into account the numerical uncertainty due to the statistical convergence, ini-
tial transients, effect of sub-iterations (nIt) and effect of the time step used for the
simulations (??s). Table 5.1 shows the total uncertainty associated with the various
flow variables. It should be noted that the banded profiles seen for all the plots from
5.8 to 5.16 show the upper and lower limit due to the numerical uncertainty. Experi-
mental uncertainty is not available for any data. The profiles from the simulations are
compared with experimental data available from two PIV windows. One of the data sets
is obtained from the vertical PIV window, in the z/W = 0 plane (red colored plane in
figure 5.1) and the other dataset is obtained from one of the horizontal strreamwise PIV
windows corresponding to either y/W = 0.348, y/W = 0.696 or y/W = 1.044 (Green,
Violet or Blue colored plane in the figure 5.1, respectively). Ideally the data from the
two PIV windows should overlap each other, but it can be seen in all the plots from
5.8 to 5.16 that there is a discrepancy in the experimental data used from the two PIV
windows. The reason for this discrepancy is not known. Figure 5.8 to 5.10 compare the
u/Vref velocity profiles at different locations in the wake. At no location is the u/Vref
profile from the simulations matching with either of the experimental datasets. Also
71
generally the recovery of the u/Vref velocity from the experiments starts earlier than for
the simulations.
x/W
u/V
re
f
7.5 8 8.5 9
-0.3
-0.2
-0.1
0
0.1
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of u/Vref velocity profile: y/W = 0.348
and z/W = 0
Figure 5.8: Comparison of u/Vref pro-
files from the experiment and simula-
tions at y/W = 0.348 and z/W = 0
x/W
u/V
re
f
7.5 8 8.5 9
-0.3
-0.2
-0.1
0
0.1
0.2 Medium MeshMedium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horzontal PIV Window
Comparison of u/Vref velocity profile: y/W = 0.696
and z/W = 0
Figure 5.9: Comparison of u/Vref pro-
files from the experiment and simula-
tions at y/W = 0.696 and z/W = 0
Figures 5.11 to 5.13 compare the v/Vref velocity profiles from the simulations with
the experimental data from the PIV windows. The profile from the simulations shows a
completely different trend than found in the data. The main reason of the mismatch of
the experimental data with the velocity profiles is the difference in locations of the lower
and upper vortices seen in figures 5.2 and 5.3 from the experiment and the simulations
respectively.
Figures 5.14 to 5.16 compare the w/Vref velocity profiles from the simulations with
the data from the PIV windows. Since these plots are in the z/W = 0 plane which is
the symmetry plane of the GTS model the cross flow velocity (w/Vref) should be ideally
72
x/W
u/V
re
f
7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2
-0.3
-0.2
-0.1
0
0.1
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of u/Vref velocity profile: y/W = 1.044
and z/W = 0
Figure 5.10: Comparison of u/Vref profiles from the experiment and simulations at y/W
= 1.044 and z/W = 0
x/W
v/V
re
f
7.5 8 8.5 9
-0.3
-0.2
-0.1
0
0.1
0.2
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experimant - Horizontal PIV Window
Comparison of v/Vref velocity profile: y/W = 0.348
and z/W = 0
Figure 5.11: Comparison of v/Vref pro-
files from the experiment and simula-
tions at y/W = 0.348 and z/W = 0
x/W
v/V
re
f
7.5 8 8.5 9
-0.3
-0.2
-0.1
0
0.1
0.2
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of v/Vref velocity profile: y/W = 0.696
and z/W = 0
Figure 5.12: Comparison of v/Vref pro-
files from the experiment and simula-
tions at y/W = 0.696 and z/W = 0
73
x/W
v/V
re
f
7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2
-0.3
-0.2
-0.1
0
0.1
0.2
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of v/Vref velocity profile: y/W = 1.044
and z/W = 0
Figure 5.13: Comparison of v/Vref profiles from the experiment and simulations at y/W
= 1.044 and z/W = 0
zero. The experimental data from the PIV windows is generally in the uncertainty limit
of the profile from the simulations.
Table 5.1: Total uncertainty values of the velocity and pressure from the simulations
Source of uncertainty Cp u/Vref v/Vref w/Vref
Initial transients 0.00246 0.01557 0.01524 0.0212
Statistical convergence 0.0047 0.00665 0.01103 0.01396
Number of sub-iterations (nIt) 0.007 0.01915 0.0150 0.02296
Time step (??s) 0.0175 0.05116 0.0430 0.00564
Total uncertainty 0.02 0.0572 0.0493 0.0346
Table 5.1 shows the various uncertainties in the velocity and pressure profiles from
the simulations. These uncertainties are the biggest absolute errors calculated by finding
the difference between the current solution and the best available solution. The total
uncertainty is obtained by summing the square the various uncertainties and finding the
74
x/W
w/
V r
ef
7.5 8 8.5 9
-0.1
0
0.1
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of w/Vref velocity profile: y/W = 0.348
and z/W = 0
Figure 5.14: Comparison of w/Vref pro-
files from the experiment and simula-
tions at y/W = 0.348 and z/W = 0
x/W
w/
V r
ef
7.5 8 8.5 9
-0.1
0
0.1
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Vertical PIV Window
Experiment - Horizontal PIV Window
Comparison of w/Vref velocity profile: y/W = 0.696
and z/W = 0
Figure 5.15: Comparison of w/Vref pro-
files from the experiment and simula-
tions at y/W = 0.696 and z/W = 0
x/W
w/
V r
ef
7.6 7.8 8 8.2 8.4 8.6 8.8 9 9.2
-0.1
-0.05
0
0.05
0.1
Medium Mesh
Medium Mesh Upper & Lower Limit
Experiment - Verical PIV Window
Experiment - Horizontal PIV Window
Comparison of w/Vref velocity profile: y/W = 1.044
and z/W = 0
Figure 5.16: Comparison of w/Vref profiles from the experiment and simulations at y/W
= 1.044 and z/W = 0
75
root of the sum.
Utotal =
radicalBig
U2stat +U2tran +U2It +U2??s (5.1)
5.3 Pressure Distribution Profile
Figure 5.17 compares the pressure distribution along the base of the GTS model at
z/W = 0.0. The reason for the discrepancy between the pressure distribution from the
experiment and the simulations can again be traced back to the difference in location of
the lower and the upper vortices near the base of the GTS model (see figures 5.2 and
5.3) from the experiment and simulations.
In the experiment the lower vortex is closer to the base of the GTS model and
hence the experiment shows a large negative pressure in the lower base region. The
upper vortex is located much farther from the base of the GTS model hence the pressure
in the upper part of the base is higher. The medium mesh simulations on the other
hand have the lower and upper vortices almost of equal size and equidistant from the
base of the model. Therefore the pressure distribution from the simulations is almost
a flat curve having a minor bulge in the negative direction at the locations where the
vortices are located. We conclude that the overall magnitude of the Cp predicted by
the simulations matches that of the Cp predicted by the experiment but the profiles do
not match. It should be noted here that the scale of the Cp profile in the figure 5.17
is enlarged to highlight the differences. Unaune et al. [29] has compared the Cp profile
76
along the base of the GTS from the simulations and the experiment on a very large scale
and hence the plot doesn?t show any details of comparison.
Cp
y/W
-0.2 -0.15 -0.1 -0.05 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Medium Mesh
Medium Mesh
Upper & Lower Limit
Experiment
Comparison of the pressure distribution along truck
base: z/W = 0
Figure 5.17: Comparison of pressure distribution along base of the truck from the ex-
periment and the medium mesh simulations
Figures 5.18 to 5.21 compare the pressure distribution over the front, top, side and
bottom surfaces of the GTS model from the simulations with the experimental data. It
77
Cp
y/W
-2 -1 0 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Medium Mesh
Experiment
Pressure distribution along the front surface at
z/W=0.0 from experiments and medium mesh
Figure 5.18: Comparison of pressure
distribution on the front surface of the
GTS model
x/W
C p
0 2 4 6 8
-0.4
-0.2
0
0.2
0.4
Medium Mesh
Experiment
Pressure distribution along the top surface of the GTS
at z/W = 0.0 from the experiments and medium mesh
Figure 5.19: Comparison of pressure
distribution on the top surface of the
GTS model
x/W
C p
-2 0 2 4 6 8-0.6
-0.4
-0.2
0
0.2
Medium Mesh
Experiment
Pressure distribution along the side of the GTS (z/W=0.5)
at y/W = 0.696 from the experiments and medium mesh
Figure 5.20: Comparison of pressure
distribution along the right side of the
GTS model
x/W
C p
0 2 4 6 8
-1.5
-1
-0.5
0
Medium Mesh
Experiment
Pressure distribution along the bottom surface of GTS
at z/W =0.0 from the experiment and the medium mesh
Figure 5.21: Comparison of pressure
distribution on the bottom surface of
the GTS model
78
should be noted that the flow over these surfaces is attached, unlike the massively sepa-
rated near wake flow. The pressure distribution obtained from the simulations matches
very well with the experiment in these attached boundary layer regions.
5.4 Drag, Lift and Side Force Coefficients
Table 5.2 compares the drag, side and lift force coefficients predicted by the coarse
and medium mesh DES and other drag coefficient values available in the literature with
the experimental data. The wind tunnel and the GTS model constructed in the simula-
tions is symmetric about the z/W = 0.0 plane so the side force predicted by the coarse
and the medium mesh simulations should ideally be zero. The small non-zero values
represent statistical error. Figures 5.22 and 5.23 from the coarse and medium meshes
respectively, show the variation of the instantaneous and the mean drag with time. Note
the higher amplitude oscillations of the medium mesh compared to the coarse mesh.
This may be due to the better spatial resolution provided by the medium mesh.
Table 5.2: Comparison of the drag, side and lift force coefficients
(CDmean) (CSmean) (CLmean)
Experiment 0.249?0.01 0.007 -0.158
Coarse Mesh DES 0.2565 -1e-4 -0.1187
Medium Mesh DES 0.2439 -3e-5 -0.1224
Unaune et al. [29] DES 0.253 NA NA
Maddox et al. [18] DES 0.279 NA NA
Roy et al. [21] (RANS Menter k??) 0.298 NA NA
Roy et al. [21] (RANS Spalart-Allmaras) 0.413 NA NA
79
?
C D
350 400 450 500 5500.1
0.15
0.2
0.25
0.3
0.35
0.4
CDmean
CD
Coarse Mesh Drag Coefficient
Figure 5.22: Coarse mesh mean and in-
stantaneous drag coefficient
?
C D
300 350 400 450 5000.1
0.15
0.2
0.25
0.3
0.35
0.4
CDmean
CD
Medium Mesh Drag Coefficient
Figure 5.23: Medium mesh mean and
instantaneous drag coefficient
5.5 Vortex Core Locations
Figures 5.24, 5.24 and 5.26 compare the location of the vortex cores in vertical and
horizontal streamwise plane between the experiment and the simulations. The location
of the upper vortex cannot be seen in the PIV window in the vertical streamwise plane
(red colored plane in figure 2.5) from the experiment but it may be located in the upper
right corner outside the PIV window shown in figure 5.2. While the prediction of the
location of the vortex cores improves from the coarse mesh to the medium mesh in the
vertical streamwise plane (z/W = 0), the same cannot be said about the horizontal
streamwise planes. However it should be noted that the vortex cores in the horizontal
planes from the medium and coarse mesh are already very close to the experiment as
compared to what we see in the vertical streamwise plane.
80
x/W
y/W
7.5 8 8.5 9 9.5
0
0.5
1
1.5
Coarse Mesh
Medium Mesh
Experiment
Comparison of the location of vortex cores from
Experiment, Coarse and Medium mesh: z/W = 0
Figure 5.24: Location of vortex cores in
vertical streamwise plane (z/W = 0.0)
x/W
z/W
7.5 8 8.5 9 9.5-1
-0.5
0
0.5
1
Coarse Mesh
Medium Mesh
Experiment
Comparison of the location of vortex cores from
Experiments, Coarse and Medium mesh : y/W = 0.69
Figure 5.25: Location of vortex cores
in horizontal streamwise plane (y/W =
0.696)
x/W
z/W
7.5 8 8.5 9 9.5-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Coarse Mesh
Medium Mesh
Experiment
Comparison of the location of vortex cores from
Experiment, Coarse and Medium mesh: y/W=1.044
Figure 5.26: Location of vortex cores in horizontal streamwise plane (y/W = 1.044)
81
5.6 Turbulent Eddies in the Wake
Figure 5.27 shows the instantaneous z-vorticity in the vertical streamwise plane
(top) and the isosurfaces of z-vorticity in the wake of the GTS model (bottom) for the
medium grid. The z-vorticity is a measure of the rotation of the flow about the z-axis,
Figure 5.27: Vortical structures resolved using LES in the wake of the GTS on the
medium mesh
red counter-clockwise rotation and blue denotes the clockwise rotation. Figure 5.27 gives
82
an indication of the size and nature of the turbulent structures in the wake resolved by
the DES model on the medium mesh (13.2 million grid cells).
5.7 Power Spectra of the Unsteady Pressure Signal
In the experiment there was one high-frequency pressure transducer located on the
base of the GTS model (x/W = 7.647, y/W = 0.63 and z/W = -0.46). The power spectra
for the unsteady pressure signal is obtained by performing a fast Fourier transformation
(FFT) of the normalized pressure signal from the experiment and the simulations. FFT
is an efficient way to perform a discrete Fourier transform of a signal which converts
the signal from the time domain into the frequency domain. The power spectra serves
the purpose of giving an idea about the frequency of the turbulent structures in the
region near the base of the GTS model. The frequency corresponding to the peak in
the power spectra is the frequency containing the most turbulent energy. FFT of the
pressure signals at other locations in horizontal streamwise plane (y/W = 0.696) and
vertical streamwise plane (z/W = 0) from the coarse mesh simulations is also performed.
A comparison of the power spectra from the experiment and from the simulations will
give an idea about the frequency of the turbulent eddies from the experiment and the
simulations. Figures 5.28, 5.29 and 5.30 show the power spectra obtained from the exper-
iment and the corresponding spectra from coarse mesh and medium mesh simulations,
respectively. The y axis in the plots denotes the energy associated with the unsteady
83
turbulent wake in units of decibel (dB). The power spectra obtained from the experi-
ment shows a peak frequency at nearly 24Hz. The red line shown in all three plots is
the best fit curve of the FFT of the experimental data. It shows the general trend of the
FFT of the experimental data. The FFT plots from the simulations (Figures 5.29 and
5.30) show three peaks at approximately 80 Hz, 110 Hz and 130 Hz. Also some peaks
are seen in the range of 10 to 30 Hz which are highlighted by the circle and appear to
roughly correspond to the dominant frequencies found in the experiment. Figures 5.31
to 5.34 show power spectra of the pressure signals at four locations in the wake from
the coarse mesh simulations. These locations at x/W ? 8.15 (approximately W/2 away
from the trailer base) in the regions where the attached boundary layers on the side, top
and bottom surfaces of the GTS model separate from the trailer base. Figures 5.31 and
5.32 show the spectra obtained from the pressure signals in the separated side-wall shear
layers; however, there are no distinct peaks seen in these two plots. Figures 5.33 and 5.34
give the power spectra in the shear layers which separate from the bottom and top of the
truck, respectively and show distinct peaks which roughly correspond to the peaks seen
in the figures 5.29 and 5.30. Animations of the turbulent flow in the wake were performed
and they showed the mixing of the shear layers from the separated boundary layers on
the side, top and bottom surfaces of the GTS with the LES region in the wake. The
peaks seen in figure 5.33 and 5.34 appear to correspond to the large-scale unsteadiness
of the shear layers seen in the animations. The peaks seen in the FFT plots from the
simulations thus seem to correspond to the dominant shedding frequency of the shear
84
layers generated by the boundary layers on the top, side and bottom surfaces of the GTS
model. This problem does not appear in the experimental data and highlights one of
basic drawbacks of the DES model, namely that there is no mechanism for converting the
modeled turbulence from the attached RANS boundary layers into resolved fluctuations
in the LES wake region. Strouhal number based on the width of the trailer was also
calculated. The FFT performed on the normalized drag plot gave a peak frequency of ?
132 Hz. Since the FFT plots of the pressure signal also showed peak frequency of ? 132
Hz this frequency is used in the calculation of Strouhal number. The Strouhal number
from the computations is ? 0.45. The Strouhal number from the experimental data is
0.08.
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Experiment FFT
Experiment FFT - best fit curve
Power spectra of unsteady pressure from experiments
Figure 5.28: Fast Fourier transform of the pressure signal from experimental data (x/W
= 7.647, y/W = 0.63 and z/W = -0.46)
85
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Coarse mesh
Experiment FFT - best fit curve
Power spectra of unsteady pressure from Coarse
mesh
Figure 5.29: Fast Fourier transform of
the pressure signal from coarse mesh
simulations (x/W = 7.647, y/W = 0.63
and z/W = -0.46)
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Medium mesh FFT
Experiment FFT - best fit curve
Power spectra of unsteady pressure from Medium
mesh
Figure 5.30: Fast Fourier transform of
the pressure signal from medium mesh
simulations (x/W = 7.647, y/W = 0.63
and z/W = -0.46)
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Power spectra of pressure signal in horizontal stream-
wise plane at x/W = 8.158, y/W = 0.696 and z/W=0.5034
Figure 5.31: Fast Fourier transform of
the pressure signal in the horizontal
streamwise plane from the coarse mesh
simulations
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Power spectra of pressure signal in horizontal stream-
wise plane at x/W = 8.158, y/W = 0.696 and z/W = -0.5034
Figure 5.32: Fast Fourier transform of
the pressure signal in the horizontal
streamwise plane from the coarse mesh
simulations
86
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Power spectra of pressure signal in vertical streamwise
plane at x/W = 8.158, y/W = -0.036 and z/W = 0
Figure 5.33: Fast Fourier transform of
the pressure signal from the vertical
streamwise plane in the coarse mesh and
near the bottom surface of the GTS
model
Frequency (Hz)
En
er
gy
(d
B)
0 50 100 150 200 250 300
-40
-20
0
20
Power spectra of pressure signal in vertical streamwise
plane at x/W = 8.158, y/W = 1.392 and z/W = 0
Figure 5.34: Fast Fourier transform of
the pressure signal from the vertical
streamwise plane in the coarse mesh and
near the top surface of the GTS model
87
Chapter 6
Conclusion
The initial transient period of ? ? 300 (t ? 1 s) was quite long, thus making the
simulations computationally expensive than expected. The time window required for
achieving statistically converged results was found to be ??w ? 180 (?tw ? 0.6 s).
These estimates for the initial transient period and time window required for statistical
convergence will help in future DES simulations. Based on this finding we concluded
that similar studies made by Unaune et al. [29] and Maddox et al. [18] have not
marched long enough in the physical time during their computations to get out of this
initial transient period and neither have they collected the statistics for a long enough
period to get statistically converged results. Analysis of the number of sub-iterations
(nIt) required to achieve an iteratively converged solution suggested that the results
from the case with 4 sub-iterations are reasonably close to the solution obtained from 8
sub-iterations; however computations performed with 1 sub-iteration did not give good
results. It is concluded that 4 sub-iterations are sufficient for achieving accurate results
when a non-dimensional time step of ??s = 1.16 ? 10?3 (?ts = 4 ? 10?6 s) is used.
It should be noted that this analysis was performed using the coarse mesh simulations.
It is computationally expensive to use the medium mesh for this analysis, so we have
assumed that the conclusions drawn for the required number of sub-iterations hold true
for the medium mesh results as well. Analysis of the the time step (??s) that should be
88
used for the simulations suggested that the results from the ??s = 2.32 ? 10?3 (?ts =
8 ? 10?6 s) case were significantly different than the results from the ??s = 1.16 ? 10?3
(?ts = 4 ? 10?6 s) case. Thus a maximum time-step of ??s = 1.16 ? 10?3 should be
used instead of ??s = 2.32 ? 10?3 to get accurate results. The grid size has a significant
effect on the results from the simulations. The medium mesh results are much better
than the coarse mesh simulations.
The drag value predicted by the simulations matches very well with the experimental
value. The pressure values predicted by the DES in the attached boundary layer regions
also match very well with the experiment. The range of the pressure values predicted by
the DES approach in the base of the GTS model matches with the range of the values
of pressure from the experiments. The actual pressure profile, however, does not match
with the experimental data. The prediction of the location of the vortices near the base
of the model does not match with the experimental data. Here the DES model fails
to accurately capture the details of the turbulent flow in the near base region. The
inability of these DES simulations to accurately predict the details of the turbulent wake
structure suggests that this model may not be able to accurately predict the drag when
drag reduction devices are added to the base. The failure of the current DES simulations
to predict wake details appear to be related to one of the main drawbacks of the hybrid
RANS/LES methods, namely that there is no mechanism to transfer information from
the RANS region (i.e., the boundary layer) to the LES region (i.e., the wake). Hence
when the shear layers in the trailer base enter the wake region governed by the LES
89
model, there is no way to convert the large, modeled eddy viscosity into large-scale
turbulent fluctuations required for the LES model. Animations of these RANS-based
shear layers confirmed that the large turbulent energy peaks seen in the 80-150 Hz range
correlate with the frequency of the shear layer flapping instability. The experimental
data did not show these features. The DES model is cheaper than the LES approach
and the Direct Numerical Simulations, but is still computationally expensive.
90
Chapter 7
Recommendations
It is expected that the fine mesh simulations will give better comparison with ex-
periments, hence fine mesh simulations should be performed. Research on the transfer of
information from the modeled turbulence obtained from RANS equations in the bound-
ary layer to the LES model in the separated region should improve the predictions of
the wake flow details using DES on coarser grids. The front portion of the coarse mesh
should be truncated and the results from the truncated coarse mesh should be compared
with the results from the full geometry coarse mesh. If these simulations show good
agreement then it can be concluded that the front portion of the GTS does not play a
significant role in the turbulent flow in the near-base region. Then the front portion of
the GTS can be truncated to make the fine grid simulations computationally less expen-
sive. The effect of drag reduction devices such as the boattail plates which are attached
to the base of the model can then be studied using DES on the truncated mesh. Finally,
truncating the front portion of the GTS mesh will allow the attached RANS boundary
layers to be perturbed so they contain the proper large-scale turbulent fluctuations when
they reach the LES region.
91
Bibliography
[1] MacCallen, R., Salari, K., Ortega, J., Browand, F., Hammache, M., Hsu, T.,
Leonard, A., Chatelain, P., Rubel, M., Roy, C. J., DeChant, L., Hassan, B., Hei-
neck, J. T., Yaste, D., Storms, B., Englar, R., and Pointer, D., 2004, ?DOE?s Effort
to Reduce Truck Aerodynamic Drag - Joint Experiments and Computations Lead
to Smart Design,? 34th AIAA Fluid Dynamics Conference, June 28 - July 1.
Portland, OR.
[2] Cooper, K. R., 2003, ?Truck Aerodynamics Reborn - Lessons from the Past,? SAE
Technical Paper Series 2003-01-3376.
[3] US DOE Transportation Energy Data Book: Edition 23, 2003.
[4] Storms, B. L., James, J. C., Heineck, J. T., Walker, S. M., Driver, D. M., and Zilliac,
G. G., 2001, ?An Experimental Study of the Ground Transportation System (GTS)
Model in the Nasa Ames 7 - by 10 - ft Wind Tunnel,? NASA/TM-2001-209621.
[5] Storms, B. L., Satran, D. R., Heineck, J. T., and Walker, S. M., 2004, ?A Study
of Reynolds Number Effects and Drag-Reduction Concepts on a Generic Tractor-
Trailer,? AIAA 2004-2251.
[6] Wilcox, D. C., ?Turbulence Modeling for CFD? DCW Industries Inc, 5354 Palm
Drive, La Ca?nada, CA, November 1994, pp. 5-19, 322-327.
[7] Camarri, S., Salvetti, M. V., Koobus, B., and Dervieux, A., 2004 ?A hybrid
RANS/LES approach applied to Bluff-body flow simulation,? European Congress
on Computational Methods in Applied Sciences and Engineering, Jyvaskyla, Fin-
land, July 24-28.
[8] Spalart, P. R., Jou, W. H., Strelets, M., and Allmaras, S., 1997, ?Comments on the
Feasibility of LES for Wings and on a Hybrid RANS/LES Approach,? Advances in
DNS/LES, Columbus,OH, Greyden Press.
[9] Vino, G., Watkins, S., Mousley, P., Watmuff, P., and Prasad, S., 2005 ?Flow Struc-
tures in the near-wake of the Ahmed model,? Journal of Fluids and Structures,
Vol. 20, pp. 673-695.
[10] Martinuzzi, R. J., Bailey, S. C., and Knopp, G. A., 2003 ?Influence of Wall Proximity
on Vortex Shedding from a Square Cylinder,? Experiments in Fluids, Vol. 34, pp.
585-596.
92
[11] Martinuzzi, R. J., and Tropea, C., 1993 ?The Flow Around Surface-Mounted, Pris-
matic Obstacles Placed in a Fully Developed Channel Flow,? Journal of Fluids
Engineering, Vol. 115, pp. 85-92.
[12] Lyn, D. A., Einav, S., Rodi, W., and Park, J., 1995 ?A Laser-Doppler Velocimetry
Study of Ensemble-averaged Characteristics of the Turbulent near Wake of a Square
Cylinder,? Journal of Fluids Mechanics, Vol. 304, pp. 285-319.
[13] Ramesh, V., Vengadesan, S., and Narasimhan J. L., 2005 ?3D Unsteady RANS Sim-
ulations of Turbulent Flow Over Bluff Body by Non-linear Model,? International
Journal of Numerical Methods for Heat & FluidFlow, Vol. 16, No. 6, pp. 660-
673.
[14] Rodi, W., 1997 ?Status of Large Eddy Simulation: Results of a Workshop,? Journal
of Fluids Engineering, Vol. 119, pp. 248-262.
[15] Nakayama, A., and Vengadesan, S. N., 2002 ?On the Influence of Numerical Schemes
and Sub-grid Stress Models on Large-eddy Simulation of Turbulent Flow Past
Square Cylinder,? International Journal of Numerical Methods for Heat &
Fluid Flow, Vol. 38, No. 3, pp. 227-253.
[16] Spalart, P. R., and Squires, K. D., 2004 ?The Status of Detached-Eddy Simulation
for Bluff Bodies,? The Aerodynamics of Heavy Vehicles: Trucks, Buses, and Trains,
R. McCallen, F. Browand and J. Ross, eds., Springer, 19.
[17] Nikitin, N. V., Nicoud, F., Wasistho, B., Squires K. D., and Squires, K. D., 2000
?An Approach To Wall Modeling in Large Eddy Simulations,? Physics of Fluids,
Vol. 12, pp. 7-10.
[18] Maddox, S., Squires, K. D., Wurtzler, K. E., and Forsythe, J. R., 2004 ?Detached-
Edy Simulation of the Ground Transportation System,? The Aerodynamics of Heavy
Vehicles: Trucks, Buses, and Trains, R. McCallen, F. Browand and J. Ross, eds.,
Springer, 19.
[19] Squires, K. D., Forsythe, J. R., and Spalart, P. R., 2005 ?Detached-Eddy Simula-
tion of the Separated Flow Over a Rounded-Corner Square,? Journal of Fluids
Engineering, Vol. 127, pp. 959-966.
[20] Roy, C. J., DeChant, L. J., Payne, J. L., and Blottner, F. G., 2003 ?Bluff-Body
Flow simulations Using Hybrid RANS/LES,? AIAA 2003-3889.
[21] Roy, C. J., Payne, J. L., and Payne, M. M, 2006 ?RANS Simulations of a Simplified
Tractor/Trailer Geometry,? Journal of Fluids Engineering, Vol. 128, pp. 1083-
1089.
93
[22] L?ubcke, H., Schmidt, St., Rung, T., and Thiele, F., and 2001 ?Comparison of LES
and RANS in Bluff-body Flows,? Journal of Wind Engineering and Industrial
Aerodynamics, Vol. 89, pp. 1471-1485.
[23] Krajnovi?c, S., and Davidson, L., 2002 ,?Large-Eddy Simulation of the Flow Around
a Bluff Body,? AIAA Journal, Vol. 40, No. 5 pp. 927-935.
[24] Davidson, L., 1997 ,?Large-Eddy Simulation: A Dynamic One-Equation Sub-grid
Model for Three-Dimensional Recirculating Flow,? 11th International Symposium
on Turbulent Flows in Complex Geometry, Vol. 3, pp. 26.1-26.6.
[25] Menon, S., Kim, W.-W., 1996, ?High Reynolds Number Flow Simulations Using
the Localized Dynamic One-Equation Subgrid Model,? AIAA 96-0425.
[26] Sohankar, A., Davidson, L., and Norberg, C., 2000,?Large Eddy Simulation of Flow
Past a Square Cylinder: Comparison of Different Subgrid Scale Models,? Journal
of Fluids Engineering, Vol. 122, pp. 39-47.
[27] Nichols, R. H., 2006 ,?Comparison of Hybrid Turbulence Models for a Circular
Cylinder and a Cavity,? AIAA Journal Vol. 44, No. 6. pp. 1207-1219.
[28] Barone, M. F., Roy, C. J., 2005 ,?Evaluation of Detached Eddy Simulation for
Turbulent Wake Applications,? AIAA Paper 2005-0504.
[29] Unaune, S. V., Sovani, S. D., and Kim, S., 2005 ,?Aerodynamics of a Generic Ground
Transportation System: Detached Eddy Simulation,?The Aerodynamics of Heavy
Vehicles: Trucks, Buses, and Trains, R. McCallen, F. Browand and J. Ross, eds.,
Springer, 19.
[30] Baysal, O., and Bayraktar, I., 2000 ,?Computational Simulations for the External
Aerodynamics of Heavy Trucks,? SAE Technical Paper Series 2000-01-3501.
[31] Rodi, W., 2000 ?Comparison of LES and RANS Calculations of the Flow Around
Bluff Bodies,? Journal of Wind Engineering and Industrial Aerodynamics, Vol.
69-71, pp. 55-75.
[32] Yee, H. C., Sandham, N. D., and Djomehri, M. J., 1998 ?Low-dissipative High-
Order Shock-Capturing Methods Using Characteristic-Based Filters,? Journal of
Computational Physics, Vol. 150, pp. 199-238.
[33] Harten, A., 1984 ?On Class of High Resolution Total-Variation-Stable Schemes,?
SIAM Journal of Numerical Analysis, Vol. 21, No. 1, pp. 1-23.
[34] Kato, M., and Launder, B., 1993 ?The modeling of Turbulent Flow around Station-
ary and Vibrating Square Cylinders,? 9th International Symposium on Turbulent
Shear Flows, Kyoto.
94
[35] Norris, H., and Reynolds, W., 1975 ?Turbulent Channel Flow with a Moving Wavy
Boundary,? Department of Mechanical Engineering. Report FM/10, Stanford Uni-
versity.
[36] Wong, C. C., Blottner, F. G., Payne, J, L., and Soetrisno, M., 1995 ?Implementation
of a Parallel Algorithm for Thermo-Chemical Non equilibrium Flow Solutions,?
AIAA Paper 95-0152.
[37] INCA User?s Manual, Version 2.0, 1995, Amtec Engineering, Inc., Bellevue, WA.
[38] Payne, J, L., and Walker, M. A., 1995 ?Verification of Computational Aerodynamic
Predictions for Complex Hypersonic Vehicles using the INCA Code,? AIAA Paper
95-0762.
[39] Roy, C. J., Gallis, M, A., Bartel, T. J., and Payne, J. L., 2003 ?Navier-Stokes
and DSMC predictions for Laminar Hypersonic Shock-Induces Separations,? AIAA
Journal Vol. 41, No. 6. pp. 1055-1063.
[40] Germano, M., Piomelli, U., Moin, P., and Cabot, W., 1991 ?A Dynamic Sub-grid
scale Eddy Viscosity Model,? Physics of Fluids Vol. 3, No. 7. pp. 1760-1765.
[41] Smagorinsky, J., 1963 ?General Circulation Experiments with the primitive Equa-
tions, Part I: The Basic Experiment,? Monthly Weather Review Vol. 91, pp. 99-165
[42] Ghosal, S., Lund, T., Moin, P., and Akselvoll, K., 1995 ?A Dynamic Localiza-
tion Model for Large-Eddy Simulation of Turbulent Flows,? Journal of Fluid
Mechanics., Vol. 286, pp. 229-225
[43] Spalart, P. R., and Allmaras, P. R., 1994 ?A One-Equation Turbulence Model for
Aerodynamic Flows,? La Recherche Aerospatiale Vol. 1, pp. 5-21.
[44] Yoon, S., and Jameson, A., 1987 ?An LU-SSOR Scheme for the Euler and Navier-
Stokes Equations,? AIAA Paper 87-0600.
[45] Chung, T. J., ?Computational Fluid Dynamics? 2002 Cambridge University Press.
[46] Peery, K. M., and Imlay, S. T., 1986 ?An Efficient Implicit Method for Solving
Viscous Multi-Stream Nozzle/Afterbody Flow Fields,? AIAA Paper 86-1380.
[47] Roy, C. J., and Blottner, F, G., 2003 ?Methodology for Turbulence Model Valida-
tion: Application to Hypersonic Flows,? Journal of Spacecrafts and Rockets Vol.
40, No. 3. pp. 313-325.
[48] Shih, T. H., and Lumley, J, L., 1997 ?A New Reynolds Stress Algebraic Equation
Model,? Computational Methods in Applied Mechanical Engineering Vol. 125,
No. 3. pp. 287-302.
[49] http://www.cfdreview.com/article.pl?sid=01/04/28/2131215
95
[50] http://en.wikipedia.org/wiki/Fast Fourier transform
96
Appendix
97
Appendix A
Fortran Code for creating animations using Tecplot
Following is a Fortran code used during creating animations using Tecplot
Program DES-Coarse-Mesh-OPLTF-Animat
implicit double precision(a-h,o-z)
OPEN(UNIT=13,FILE=?OPLTF11?,STATUS=?unknown?)
OPEN(UNIT=14,FILE=?Reduced-OPLTF11?,STATUS=?unknown?)
write(14,*)?TITLE =?INCA Version 2 Input File Created by Gridge ??
write(14,*)?VARIABLES = ?X(M)? ?Y(M)? ?Z(M)? ?U(M/S)? ?V(M/S)?
& ?W(M/S)? ?P(N/M2)? ?T(K)? ?R(KG/M3)? ?Mach??TVisc??
i1=1
do itstep=1,5
read (13,*)
write(*,*)itstep
read(13,*)
read(13,*)
i1=i1+3
do nzone=1,46
read(13,11)zone,ii,jj
i1=i1+1
if ((mod(ii,5)).eq.0)then
ians1=0
ians2=ii/5
else
ians1=mod(ii,5)
ians2=(ii-ians1)/5
end if
read(13,*)
read(13,*)
i1=i1+2
do nvar=1,11
do j=1,jj
do i=1,ians2
read(13,*)
i1=i1+1
end do
if(ians1.ne.0)then
98
read(13,*)
i1=i1+1
end if
end do
end do
end do
end do
iin1=0
write(*,*)i1
do itstep=6,100000000
write(*,*)itstep
read(13,*,end=26)
read(13,*)
read(13,*)
i1=i1+3
do nzone=1,46
read(13,11,end=26)zone,ii,jj
read(13,*)
read(13,*)
i1=i1+3
if(itstep.eq.(6+5*iin1))then
write(14,20)zone,ii,jj
20 format(?ZONE T=?K=?,I7,?, ?, I= ?,I7,?, J= ?,I7,?,
F=BLOCK ?)
write(14,*)?DT=(SINGLE SINGLE SINGLE SINGLE SINGLE SINGLE SING
LE SINGLE SINGLE SINGLE SINGLE)?
end if
if ((mod(ii,5)).eq.0)then
ians1=0
ians2=ii/5
else
ians1=mod(ii,5)
ians2=(ii-ians1)/5
end if
do nvar=1,11
do j=1,jj
do i=1,ians2
! write(*,*)itstep
! write(*,*)?reached here?
! write(*,*)i1
read(13,21)v1,v2,v3,v4,v5
if(itstep.eq.(6+5*in1))then
99
write(14,21)v1,v2,v3,v4,v5
end if
end do
if(ians1.ne.0)then
if(ians1.eq.1)then
read(13,22)v6
if(itstep.eq.(6+5*in1))then
write(14,22)v6
end if
end if
if(ians1.eq.2)then
read(13,23)v6,v7
if(itstep.eq.(6+5*in1))then
write(14,23)v6,v7
end if
end if
if(ians1.eq.3)then
read(13,24)v6,v7,v8
if(itstep.eq.(6+5*in1))then
write(14,24)v6,v7,v8
end if
end if
if(ians.eq.4)then
read(13,25)v6,v7,v8,v9
if(itstep.eq.(6+5*in1))then
write(14,25)v6,v7,v8,v9
end if
end if
end if
end do
end do
end do
iin1=iin1+1
end do
11 format(13X,I2,9X,I2,7X,I2)
21 format(5(1x,E15.8))
22 format(1x,E15.8)
23 format(2(1x,E15.8))
24 format(3(1x,E15.8))
25 format(4(1x,E15.8))
26 continue
100
end program
101