University of Southern California Dissertations and Theses

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

(USC Thesis Other)

Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NUMERICAL MODELING OF SEPARATED FLOWS AT MODERATE
REYNOLDS NUMBERS APPROPRIATE FOR TURBINE BLADES AND
UNMANNED AERO VEHICLES
by
Giacomo Castiglioni
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(AEROSPACE ENGINEERING)
August 2015
Copyright 2015 Giacomo Castiglioni
Dedication
to my grandparents
ii
Acknowledgments
I want to thank my PhD advisor Professor Domaradzki for his support during this
long journey. His work ethic, his knowledge of turbulence models and healthy skepticism
towards `dissipation-less' schemes transformed my approach to CFD. His dry wit and
sense of humor kept things light at all times. Finally I would like to thank him for giving
me the chance to participate in the summer programs in Munich, which I very much
enjoyed. On that note, I would like to thank Stefan for letting me use INCA and his
continued help with the code, and all the other nice people I met at TUM.
Bogdan, a friend and a mentor. The HVAC team at Tesla Motors, in Particular Jo,
Chris, and Andrew for showing me the `real world CFD analyst' work and the trade
secrets of Star-CCM+.
Many former and current `inmates' made the journey more bearable. Informal tech-
nical and non-technical discussions with the `Indians from India' and the Greek need a
particular mention. Prabu, Stelios, Jagan, Vyaas, Dean, DJ, and Orlando; I will always
treasure our countless Friday beer meetings. Trevor, for introducing me to the intricacies
of the USC cluster and for teaching me an invaluable lesson `Never trust what someone
says, including myself'. Ramtin, for our daily conversations from the rst semester till
the end. Friends from the Old World that are always eager to meet me (possibly not
iii
more than twice a year!) whenever I am back home: Marta, Ezra, Samuele, Paolino, Ale,
Brio, Davide, Crema, and Frizzo. New friends from the Old World that I feel like I have
known forever: Andrea and Matteo.
My family for being an innite source of mental, emotional (and economic!) stability;
il papa' e la mamma sono insostituibili. Michele, Alessandra, Mattia, Simone, Emanuele,
Daniela e Luigi; fratelli e sorelle, cugini, zii che mi hanno sempre fatto sentire amato.
Martina, my girlfriend, for patiently listening to my daily rants, always believing in me,
and sharing her life with mine.
iv
Table of Contents
Dedication ii
Acknowledgments iii
List Of Symbols vii
List Of Figures xii
List Of Tables xvi
Abstract xvii
Chapter 1: Introduction 1
Chapter 2: Physical problem 5
2.1 Laminar Separation Bubble . . . . . . . . . . . . . . . . . . . . . . . . . . 5
Chapter 3: Numerical problem 11
3.1 Direct Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Large Eddy Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Reynolds Averaged Navier-Stokes . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Wall modeling in Large Eddy Simulations . . . . . . . . . . . . . . . . . . 24
3.5 Types of numerical error . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.6 Interplay between numerical dissipation and SGS model . . . . . . . . . . 30
Chapter 4: Research goals 34
Chapter 5: Numerical methods 36
5.1 INCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.1.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5.1.2 Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.1.3 Finite Volume Discretization . . . . . . . . . . . . . . . . . . . . . 39
5.1.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.2 Star-CCM+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2.1 explicit LES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.2 Dynamic Smagorinsky model . . . . . . . . . . . . . . . . . . . . . 45
v
5.2.3 WALE model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Chapter 6: Simulation results 48
6.1 Benchmark
ow conguration . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.2 Computational setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.2.1 INCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.2.2 Star-CCM+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.3 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.4 Evaluation metrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6.5 Two-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.5.1 INCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.5.2 Star-CCM+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.6 Three-dimensional results . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.6.1 INCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.6.2 Star-CCM+ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.7 Two-dimensional sweep of angles of attack . . . . . . . . . . . . . . . . . . 69
Chapter 7: Numerical dissipation 74
7.1 Concept . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
7.2 Equations and method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.2.1 Analytical form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.2.2 Discretized form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.3 Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.4 Taylor Green vortex
ow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.4.1 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.4.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
7.4.3 Full periodic domain . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.4.4 Sub-domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.4.5 Surface versus Volume Integrals . . . . . . . . . . . . . . . . . . . . 89
7.4.6 Eects of Compressibility . . . . . . . . . . . . . . . . . . . . . . . 93
7.5 Discussion of the validation case . . . . . . . . . . . . . . . . . . . . . . . 93
7.6 NACA-0012 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
7.6.1 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Chapter 8: Conclusions 106
8.1 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
8.2 Numerical dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
8.3 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
References 111
vi
List Of Symbols
Abbreviations
1-D Three-dimensional
2-D Three-dimensional
3-D Three-dimensional
APG Adverse Pressure Gradient
ALDM Adaptive Local Deconvolution Method
CFL Courant-Friedrichs-Lewy number
DNS Direct Numerical Simulations
FV Finite Volume
HCF High Cycle Fatigue
HLLC Harten, Lax, van Leer, and contact
IB Immersed Boundary
LES Large Eddy Simulations
LODI Local Associated One-dimensional Case
LSB Laminar Separation Bubble
MAV Micro Air Vehicle
NSCBC Navier-Stokes Characteristics Boundary Conditions
vii
RANS Reynolds Averaged Navier-Stokes Equations
R-K Runge-Kutta
SGS Sub Filter Scale
TVD Total-variation Diminishing
UAV Unmanned Aerial Vehicle
uDNS under-resolved Direct Numerical Simulations
WALE Wall-adapting local eddy viscosity
WENO Weighted essentially non-oscillatory
Greek Letters
Angle of attack
Volume fraction of
uid
Smagorinsky model kernel at test lter level
Smagorinsky model kernel at grid lter level
Identity tensor

Specic heats ratio
Dynamic viscosity

T
Turbulent eddy dynamic viscosity
Kinematic viscosity
Generic eld variable
Mean eld variable
SGS stress tensor

dil
SGS pressure dilatation
Density
viii
Re
ection parameter
Viscous stress tensor
Roman Letters
a Speed of sound
B
i
SGS terms, with i = 1, . . . , 6
C
f
Skin friction coecient
c
p
Specic heat capacity at constant pressure
C
p
Pressure coecient
C
s
Smagorinsky constant
c
v
Specic heat capacity at constant volume
C
w
WALE constant
e
kin
Specic kinetic energy for unit mass
e Specic internal energy for unit mass
E Total energy per unit mass
F Inviscid
ux tensor
G Viscous
ux tensor
G Convolution kernel
G Filter function
h Grid spacing
h Specic enthalpy for unit mass
H Total enthalpy per unit mass
k Thermal conductivity
K Rudy and Strikwerda constant
ix
L Germano identity tensor
L Reference length
L Ingoing wave amplitude
M Jacobian of the least square method
M
max
Expected maximum Mach number in the domain
P Lumped resolved pressure
p Pressure
q Heat
ux vector
Q SGS heat
ux vector
S Strain rate tensor
S Square of the velocity tensor
T Temperature
^
T Test ltered SGS stress tensor
u Velocity vector
U Conserved variables vector
V Volume
x Cartesian coordinate vector
Subscripts
[ ]
c
Quantity based on cord length
[ ]
i
Component of a vector
[ ]
ij
Component of a tensor
[ ]
1
Free stream reference quantity
[ ]
sgs
SGS quantity
x
Superscripts
[ ]
d
Deviatoric part of a tensor
g
[ ] Favre ltered eld variable
g
[ ] Favre averaged eld variable
[ ] Filtered eld variable
d
[ ] Test ltered eld variable
[ ]
+
Non-dimensional quantity
[ ]
0
Fluctuating component
[ ]
0
Non-resolved component
[ ]
00
Favre averaging remainder component

[ ] Resolved eld variable
[ ]
T
Transpose
Non-dimensional Numbers
Ma Mach number
Pr Prandtl number
Re Reynolds number
St Strouhal number
Symbols
hi Averaged eld variable
4 Filter cut-o
4n Cell size in the direction normal to the boundary
4x Cell size in the x-direction
4tg Cell size in the direction tangent to the boundary
xi
List Of Figures
2.1 Features of a laminar separation bubble
ow [14]. . . . . . . . . . . . . . . 8
2.2 Experimental polar plots of a E387 airfoil from Spedding and McArthur
[61]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Sketch of an averaged in time pressure coecient for a generic airfoil with a
laminar separation bubble (black dashed line) and without (red solid line).
Disregarding the contribution from the skin friction, the area enclosed by
the curve is proportional to the lift coecient and it is evident this area is
bigger in the presence of the laminar separation bubble. . . . . . . . . . . 10
3.1 Solution of the linear advection equation for dierent numerical schemes.
Computed solutions (blue line), exact solutions (red dots). . . . . . . . . . 28
3.2 Pressure plot of a simulation exhibiting aliasing error, 2h spurious waves
can be recognized. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5.1 Volume fraction of
uid for the celli;j. The light shaded area represents
the
uid domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.1 Grid around the airfoil, only every 5th line is shown (INCA, coarse grid). 50
6.2 Immersed boundary grid close-up (coarse grid). . . . . . . . . . . . . . . . 51
6.3 Grid around the airfoil (Star-CCM+ ). . . . . . . . . . . . . . . . . . . . . 52
6.4 Grid transition between hexahedral cells in the freestream and body tted
cells around the airfoil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.5 The spanwise vorticity eld (10 levels over the range75 !
+
). [14] . . . . 55
6.6 Iso-surface of Q-criterion (Q
+
= 50) colored with x-velocity (Star-CCM+,
UDNS). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
xii
6.7 Two-dimensional Pressure and friction coecients at the wall. DNS of
Jones et al. [40] (black dots), 5th order WENO scheme with ne grid
(solid red line), ALDM with ne grid (dashed dotted blue line), ALDM
with coarse grid (dashed green line), and ALDM with very ne grid (long
dashed black line). [14] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.8 Results of the dynamic mode decomposition for WENO5 and ALDM sim-
ulations. [14] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.9 Two-dimensional Pressure and friction coecients at the wall. DNS of
Jones et al. [40] (black dots), UDNS (solid red line), explicit LES with
Dynamic Smagorinsky model (dashed dotted blue line), and explicit LES
with WALE model (dashed green line). [14] . . . . . . . . . . . . . . . . . 63
6.10 Three-dimensional Pressure and friction coecients at the wall. DNS of
Jones et al. [40] (black dots), UDNS with 5th order WENO scheme with
ne grid (solid red line), DNS with 5th order WENO scheme (short dashed
black line), ALDM with ne grid (dashed dotted blue line), and ALDM
with coarse grid (long dashed green line). [14] . . . . . . . . . . . . . . . . 66
6.11 Three-dimensional Pressure and friction coecients at the wall. DNS of
Jones et al. [40] (black dots), UDNS (solid red line), explicit LES with
Dynamic Smagorinsky model (dashed dotted blue line), and explicit LES
with WALE model (dashed green line). [14] . . . . . . . . . . . . . . . . . 67
6.12 Two-dimensional Pressure and friction coecients at the wall for dierent
angles of attack (Star-CCM+, uDNS). . . . . . . . . . . . . . . . . . . . . 71
6.13 Two-dimensional Pressure and friction coecients at the wall for dierent
angles of attack (Star-CCM+, uDNS). . . . . . . . . . . . . . . . . . . . . 72
6.14 Two-dimensional Lift and Drag polar . . . . . . . . . . . . . . . . . . . . . 73
7.1 Isosurface of Q-criterion (Q = 0) colored with velocity magnitude at initial
and nal time, inviscid case. . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.2 Time-evolution of numerical dissipation rate "
n
tot
, Re =1. Star-CCM+
with 64
3
cells (blue line), Star-CCM+ with 128
3
cells (red dashed line) and
INCA from [58] (black circles). . . . . . . . . . . . . . . . . . . . . . . . . 84
7.3 Time-evolution of normalized enstrophy
(t)=
(0),Re =1. Star-CCM+
with 64
3
cells (blue line), with 128
3
cells (red dashed line), WENO (black
squares) and spectral (black circles) from [59]. . . . . . . . . . . . . . . . . 85
xiii
7.4 Time-evolution of normalized kinetic energy E
kin
(t)=E
kin
(0), Re =1.
Star-CCM+ with 64
3
cells (blue line), with 128
3
cells (red dashed line)
and t
1:2
law (black dots). . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.5 Time-evolution of numerical dissipation rate "
n
tot
, nite Re. Star-CCM+
(blue line) and INCA from [58] (black circles). . . . . . . . . . . . . . . . . 87
7.6 Time-evolution of numerical dissipation rate "
n
tot
. Re =1 (blue line),
Re = 3000 (red dashed line),Re = 1600 (black dash-dotted line),Re = 800
(magenta dotted line), Re = 400 (solid blue line with markers), Re = 200
(red dashed line with markers), andRe = 100 (black dash-dotted line with
markers) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
7.7 Time-evolution of numerical dissipation rate "
n
sub
, Re =1, sub-volume
CSD1. Numerical dissipation rate "
n
sub
(blue line), rate of change of ki-
netic energyE
kin
sub
=t (black dot-dashed line), kinetic energy
uxF
kin
sub
(blue dot-dashed line with markers), acoustic
uxF
ac
sub
(magenta dotted
line) and INCA from [58] (black circles). . . . . . . . . . . . . . . . . . . . 89
7.8 Time-evolution of numerical dissipation rate "
n
sub
, nite Re, sub-volume
CSD1. Numerical dissipation rate "
n
sub
(blue line), rate of change of ki-
netic energyE
kin
sub
=t (black dot-dashed line), kinetic energy
uxF
kin
sub
(blue dot-dashed line with markers), acoustic
uxF
ac
sub
(magenta solid
line with markers), viscous
ux F
vis
sub
(black dashed line with markers) vis-
cous dissipation rate"
vis
sub
(red dashed line) and INCA from [58] (black
circles). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
7.9 Time-evolution of physical and numerical viscosity, nite Re, sub-volume
CSD1. Total viscosity +
n
(blue line), physical viscosity (red dashed
line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
7.10 Time-evolution of numerical dissipation rate "
n
sub
, Re =1, sub-volume
NCSD2, including work due to pressure. Numerical dissipation rate "
n
sub
(blue line), rate of change of kinetic energyE
kin
sub
=t (black dot-dashed
line), kinetic energy
uxF
kin
sub
(blue dot-dashed line with markers), acous-
tic
uxF
ac
sub
(magenta solid line with markers),W
p
sub
(red dashed line with
markers) and INCA from [58] (black circles). . . . . . . . . . . . . . . . . 92
7.11 Time-evolution of numerical dissipation rate "
n
sub
, Re =1, sub-volume
CSD1 with pressure gradient transport term. Numerical dissipation rate
"
n
sub
(blue line), rate of change of kinetic energyE
kin
sub
=t (black dot-
dashed line), kinetic energy
uxF
kin
sub
(red dashed line), pressure trans-
portpt
sub
(magenta solid line with markers) and INCA from [58] (black
circles). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
xiv
7.12 Time-evolution of numerical dissipation rate"
n
sub
for a compressible solver,
Ma = 0:0845, nite Re, sub-volume CSD1. Numerical dissipation rate
"
n
sub
(blue line), viscous dissipation rate"
vis
sub
(red dashed line) and INCA
from [58] (black circles). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
7.13 Sub-domains location for NACA-0012. . . . . . . . . . . . . . . . . . . . . 95
7.14 Time-evolution of numerical and viscous dissipation rates. "
n
sub
(blue solid
line) and "
vis
sub
(red dashed line). . . . . . . . . . . . . . . . . . . . . . . . . 99
7.15 Time-evolution of physical and numerical viscosities. Total viscosity +
n
(blue solid line), physical viscosity (red dashed line) and average total
viscosity +
n
(dot-dashed black line). . . . . . . . . . . . . . . . . . . . 100
7.16 Time-evolution of single terms of the kinetic energy balance, Eq. (7.14).
Plotted are the values for sub-domain SD1 (solid blue line) against SD1a
(red dashed line). (a) change in time of kinetic energy, (b) kinetic energy

ux, (c) acoustic
ux, (d) viscous
ux. . . . . . . . . . . . . . . . . . . . . 101
7.17 Time-evolution of single terms of the kinetic energy balance, Eq. (7.14).
Plotted are the values for sub-domain SD1 (solid blue line) against SD1a
(red dashed line). (e) pressure work, (f) viscous dissipation rate, (g) dissi-
pation function and (h) numerical dissipation rate. . . . . . . . . . . . . . 102
7.18 Time-evolution of physical, numerical and sub-grid-scale viscosities with
dierent values of C
s
. Total viscosity +
n
(blue line), physical viscosity
(red dashed line), average total viscosity +
n
(dot-dashed black line)
and physical viscosity plus sub-grid scale viscosity
SGS
(blue line with
markers) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
xv
List Of Tables
6.1 Summary of two-dimensional time averaged laminar separation bubble
characteristics. N is the number of grid points, x
s
the separation point,
x
r
the reattachment point and l
b
the length of the bubble. [14] . . . . . . 68
6.2 Summary of three-dimensional time averaged laminar separation bubble
characteristics. N is the number of grid points, x
s
the separation point,
x
r
the reattachment point and l
b
the length of the bubble. [14] . . . . . . 68
6.3 Summary of two-dimensional time averaged laminar separation bubble
characteristics for the sweep of angles of attack. N is the number of grid
points,x
s
the separation point,x
r
the reattachment point andl
b
the length
of the bubble. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.1 Summary of averaged in time dissipation rates and viscosities, where
tot
=
+
n
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
xvi
Abstract
Flows over airfoils and blades in rotating machinery, for unmanned and micro-aerial
vehicles, wind turbines, and propellers consist of a laminar boundary layer near the
leading edge that is often followed by a laminar separation bubble and transition to
turbulence further downstream. Typical Reynolds averaged Navier-Stokes turbulence
models are inadequate for such
ows. Direct numerical simulation is the most reliable, but
is also the most computationally expensive alternative. This work assesses the capability
of immersed boundary methods and large eddy simulations to reduce the computational
requirements for such
ows and still provide high quality results. Two-dimensional and
three-dimensional simulations of a laminar separation bubble on a NACA-0012 airfoil at
Re
c
= 510
4
and at 5

of incidence have been performed with an immersed boundary code
and a commercial code using body tted grids. Several sub-grid scale models have been
implemented in both codes and their performance evaluated. For the two-dimensional
simulations with the immersed boundary method the results show good agreement with
the direct numerical simulation benchmark data for the pressure coecient C
p
and the
friction coecient C
f
, but only when using dissipative numerical schemes. There is
evidence that this behavior can be attributed to the ability of dissipative schemes to
damp numerical noise coming from the immersed boundary. For the three-dimensional
xvii
simulations the results show a good prediction of the separation point, but an inaccurate
prediction of the reattachment point unless full direct numerical simulation resolution is
used. The commercial code shows good agreement with the direct numerical simulation
benchmark data in both two and three-dimensional simulations, but the presence of
signicant, unquantied numerical dissipation prevents a conclusive assessment of the
actual prediction capabilities of very coarse large eddy simulations with low order schemes
in general cases. Additionally, a two-dimensional sweep of angles of attack from 0

to 5

is performed showing a qualitative prediction of the jump in lift and drag coecients due
to the appearance of the laminar separation bubble. The numerical dissipation inhibits
the predictive capabilities of large eddy simulations whenever it is of the same order
of magnitude or larger than the sub-grid scale dissipation. The need to estimate the
numerical dissipation is most pressing for low-order methods employed by commercial
computational
uid dynamics codes. Following the recent work of Schranner et al., the
equations and procedure for estimating the numerical dissipation rate and the numerical
viscosity in a commercial code are presented. The method allows for the computation of
the numerical dissipation rate and numerical viscosity in the physical space for arbitrary
sub-domains in a self-consistent way, using only information provided by the code in
question. The method is rst tested for a three-dimensional Taylor-Green vortex
ow in
a simple cubic domain and compared with benchmark results obtained using an accurate,
incompressible spectral solver. Afterwards the same procedure is applied for the rst time
to a realistic
ow conguration, specically to the above discussed laminar separation
bubble
ow over a NACA 0012 airfoil. The method appears to be quite robust and its
application reveals that for the code and the
ow in question the numerical dissipation
xviii
can be signicantly larger than the viscous dissipation or the dissipation of the classical
Smagorinsky sub-grid scale model, conrming the previously qualitative nding.
xix
Chapter 1
Introduction
Reynolds numbers for
ows in rotating machinery, for unmanned aerial vehicles (UAV),
micro air vehicles (MAV), wind turbines, and propellers are usually low or moderate.
Based on wing/blade chord they are typically less than 2 10
6
and often only on the
order of few 10
4
10
5
. By comparison, civilian airplanes are characterized by the
Reynolds numbers ranging from a few million to 80 10
6
for the Boeing 747 at cruising
velocity. Recently the use of UAVs has been proposed as possible Mars explorers, an
application for which high-lift, low Reynolds number airfoils are necessary due to the low
cruising speed required to take pictures [42]. There is a long history of research on high
Reynolds number aerodynamics, driven by its importance in the aerospace industry, and
in time robust numerical tools have been developed for the prediction of lift and drag
forces for streamlined bodies. The recent experimental investigations of low and mod-
erate Reynolds number aerodynamics [38, 35, 61] reveal new features of such
ows that
complicate their prediction compared with high Reynolds number
ows.
1
Aerodynamics at moderate and low Reynolds numbers is less developed, partly be-
cause the eld was less important before, but mostly because of a unique physical di-
culty: such
ows are often dominated by the eects of
ow separation. At low Reynolds
numbers a signicant portion of the boundary layer remains laminar before transitioning
to turbulence. The laminar portion of the boundary layer in certain cases separates due
to an adverse pressure gradient (APG). The
ow separation and reattachment locations
can greatly in
uence lift and drag, and thus aects
ight stability of UAVs, eciency of
wind turbines, and unsteadiness in turbine
ows which is instrumental in predicting high
cycle fatigue (HCF) for turbomachinery components. The physical origin of the observed
behavior is qualitatively well understood: it is caused by a sudden appearance of the
laminar separation bubble (LSB) that changes the
ow pattern around a wing and thus
its properties. In order to produce ecient designs or control schemes to reduce separa-
tion eects we need prediction tools for such
ows. However, while there is experimental
data on both the time-averaged lift and drag forces and the instantaneous
ow elds over
wings, the agreement between dierent experiments is often poor, indicating that the
separation and reattachment processes can be sensitive to details of experimental setups.
The agreement between experiments and computations is even harder to nd, indicat-
ing potential diculties in developing predictive tools. This is because the separation pro-
cess, especially driven by APG as opposed to geometry, is essentially in non-equilibrium,
involving subtle interactions between viscous, advective, and pressure eects that can
be captured only by full Navier-Stokes equations. Indeed, at present the only reliable
numerical results for such
ows have to be obtained using direct numerical simulations
(DNS). However, such DNS require substantial computational resources, long wall-clock
2
runs and long analysis times, e.g., Jones et al. [40] used over 170 million grid points
for a relatively simple 3-D airfoil conguration. In this case the airfoil studied was a
NACA-0012 airfoil at 5

angle of attach and Re
c
= 5 10
4
.
Such numerical simulations of full Navier-Stokes equations are not feasible if a number
of congurations and angles of attack are to be investigated. Therefore, other simulation
options must be considered. One option is to employ Reynolds Averaged Navier-Stokes
equations (RANS) models, modied to account for the reduction of the eddy viscosity
around the separation region. This is an approach commonly used and optimized for high
Reynolds number, turbulent
ows, but is not as useful for the separated
ows of interest.
Another option is to employ Large Eddy Simulations (LES) techniques. For instance,
Yang and Voke [70] reported LES results obtained with the dynamic Smagorinsky model
in a good agreement with experiments for boundary-layer separation and transition caused
by surface curvature atRe = 3; 450. Yet even for this relatively low Reynolds number, the
critical issue in getting the agreement was good numerical resolution (472x72x64 mesh
points) and a high order numerical method, requirements dicult to satisfy in simula-
tions of practical
ows often performed with low order nite dierence or nite volume
methods (e.g., commercial codes). Similarly, Eisenbach and Friedrich [25] performed LES
of
ow separation on an airfoil at high angle of attack atRe = 10
5
using Cartesian grids.
This case also required very high resolutions between 50 and 100 million mesh points.
Furthermore, this specic simulations were done at high angles of attack (18 deg) where
the separation location is less dicult to predict. A rare example of a low resolution LES
is given by Almutairi et al. [2]. In that study results for a laminar separation bubble
3
over an airfoil at Re = 5 10
5
are in a good agreement with Jones et al. [40] DNS data
at 4:5% of the DNS resolution.
In this thesis we report on a study of a separation problem for a
ow around a wing
at moderate Reynolds numbers using LES. The specic geometrical setting is that of a
3-D airfoil at incidence for which detailed DNS results were obtained in [40, 41]. The goal
of this research is to reproduce the laminar separation bubble on a NACA-0012 airfoil
at Re
c
= 5 10
4
at 5

of incidence with resolution reduced drastically from that used
in DNS and typical LES performed for this problem. The reduction target is a factor of
100. To simulate this case two codes have been employed. The rst code and principal
tool has been developed at the technical university of Munich (TUM) [36, 32]. It uses
a Cartesian nite volume discretization (FV) and an immersed boundary method (IB);
from now on this code will be referred to as INCA. The second code employed is the
commercial code Star-CCM+ from CD-Adapco that uses an unstructured FV method.
The rationale of using a second, commercial code to study the same problem is twofold;
rstly Star-CCM+ is used as a second opinion to evaluate dierences that might arise
between the DNS reference case and INCA, secondly comparing Star-CCM+ and the
reference DNS will provide a good sample of the current status of LES models capabilities
in commercial codes.
4
Chapter 2
Physical problem
2.1 Laminar Separation Bubble
Modern understanding of aerodynamics stems from the revolutionary work of Ludwig
Prandtl in 1904 (published in a paper in 1905) [6]. Prandtl was the rst to suggest the
presence of a thin viscous boundary layer around objects in a
ow. His intuition directly
derived from the hypothesis of a no slip condition at the wall and the fact that, in a low
viscosity
uid, the
ow should rapidly recover free stream velocity away from the wall.
Inside this boundary layer, the rapid change in velocity would lead to a signicant velocity
gradient which would in turn induce high shear stresses. Therefore, contrary to what was
believed until the early 19th century, even in a low viscosity
uid the friction forces at
the wall are not negligible. Prandtl thus suggested that the
ow around an object can be
divided into two regions: a thin boundary layer close to the object where friction forces
are dominant, and an inviscid region outside the boundary layer where friction forces are
negligible. In his short paper Prandtl went as far as giving a rst description of boundary
layer separation due to an APG, explaining the mechanism behind separation. The
uid
5
in the boundary layer, close to the wall, facing an APG and having its energy depleted
by frictional losses at a certain point (separation point) has not enough momentum to
stay attached to the wall and separates. At certain low Reynolds numbers, typically
in the range Re
c
= 10
4
10
6
, the laminar separated shear layer undergoes instability
and transition to turbulence. In this range the energized boundary layer has enough
momentum to reattach to the body. Thanks to the concept of boundary layers introduced
by Prandtl, for the rst time it was possible to quantitatively calculate lift and drag of
airfoils. Separation of the boundary layer can have dramatic eects. Once the shear layer
separates it can create large scale disturbances of the mean
ow creating a dramatic loss
in lift.
Separated
ows can be divided into two categories:
ows with geometry induced sepa-
ration and APG induced separation [3]. Geometry induced separation is somewhat easier
to deal with because the
ow separates due to a sharp geometry gradient and therefore
the separation point is xed with respect to the Reynolds number. APG separation is
more complex because not only do both separation and reattachment points depend on
several parameters (pressure gradient, curvature, Reynolds number, free stream turbu-
lence levels, etc.) but the
ow under certain values may not even separate. In both
categories the
ow might reattach to the wall. In this project APG separation with
turbulent reattachment is considered.
Nowadays, as mentioned in the introduction, the physical behavior behind the laminar
separation and transition to turbulence due to APG is qualitatively well understood.
For instance, the attached laminar boundary layer developing on a wing at incidence is
subjected to an APG due to the wing curvature that causes the shear layer to separate.
6
Immediately behind the separation point there is an eectively stagnant
ow region,
the so-called \dead air" region or \bubble", followed by a reverse
ow vortex. The main
features of the laminar separation bubble are illustrated in Figure 2.1 taken from [46]. The
interface between the separated
ow moving away from the wing and the recirculating
ow
in the vicinity of the wing results in a shear layer with an in
ectional mean velocity prole.
This shear layer experiences Kelvin-Helmoholtz instabilities, developing into turbulence,
after rst generating characteristic spanwise vortices. Further downstream, the separated
turbulent
ow reattaches and gradually evolves into the classical turbulent boundary
layer. The above description emerges from numerous experimental investigations, e.g.,
[38, 35, 61, 34] as well as from DNS results [46, 69, 60, 1, 40, 41]. Note that despite the
fact that in classical literature the term \laminar" is used, the actual
ow is a mixture
of regions where the
ow is a laminar, transitional, non-equilibrium turbulent boundary
layer, and an equilibrium turbulent boundary layer.
The appearance of a laminar separation bubble can have a dramatic eect on lift and
drag coecients as can be seen from Figure 2.2 where the polar plot of an E387 airfoil from
the experimental results of Spedding and McArthur [61] is shown. Each point represents
an angle of attack, the dierent lines are dierent Reynolds numbers. From a certain
Reynolds number (Re = 30; 000) and at a certain angle of attack a jump in lift and drag
coecients can be observed. This jump is due to the appearance of a laminar separation
bubble. To presence of the bubble eectively changes the shape of the airfoil that the
free stream experiences and the jump in the lift coecient can be more easily understood
observing the sketch in Figure 2.3. In Figure 2.3 the pressure coecient of a generic
airfoil is sketched in the presence of a laminar separation bubble (red solid line) and in
7
Figure 2.1: Features of a laminar separation bubble
ow [14].
8
Figure 2.2: Experimental polar plots of a E387 airfoil from Spedding and McArthur [61].
the absence of it (black dashed line). If we ignore the contribution from the skin friction,
the area enclosed by the curve is proportional to the lift coecient and it is evident that
the area under the curve is bigger in the presence of the laminar separation bubble.
9
Figure 2.3: Sketch of an averaged in time pressure coecient for a generic airfoil with a laminar
separation bubble (black dashed line) and without (red solid line). Disregarding the contribution
from the skin friction, the area enclosed by the curve is proportional to the lift coecient and it
is evident this area is bigger in the presence of the laminar separation bubble.
10
Chapter 3
Numerical problem
Historically low Reynolds number
ows have been less studied in CFD, on one hand
because research eort was concentrated on higher Reynolds number
ows due to the in-
terests of the aerospace industry and on the other hand because such
ows are dominated
by the eects of
ow separation that are challenging to simulate. In order to produce
ecient designs or control schemes to reduce separations eects, numerical prediction
tools for such
ows are needed. Reliable numerical simulations are dicult to obtain be-
cause the complexity of such
ows is a challenge for turbulence models: in the same
ow
regions are present where the
ow is a laminar, transitional, non-equilibrium turbulent
boundary layer and an equilibrium turbulent boundary layer. In addition, the separation
problem, especially driven by APG is essentially non-equilibrium involving subtle inter-
action between viscous, advective and pressure eects that can be fully captured only
using full Navier-Stokes equations. Ideally DNS would be the perfect tool for parametric
and optimization studies, however it is well known that nowadays its computational cost
is prohibitive. DNS studies are limited to very simple 3-D geometries, usually geomet-
rically homogenous in the spanwise direction and having small physical domains. As
11
an example, in the reference studies of Jones et al. [40, 41] the spanwise extent of the
domain is only 20% of the cord and the simulation required almost 170 million points.
Such computational costs do not allow for parametric studies and there is a need to look
at other simulation options.
One option is to employ RANS models. This technique is signicantly less expensive
than DNS. RANS modeling is commonly used and optimized for high Reynolds number,
turbulent
ows, but not as useful for the separated
ows of interest. RANS simula-
tions cannot accurately predict a priori the separation point, and both separation and
reattachment points strongly depend on free stream levels of turbulence [20]. Several
modications have been proposed to overcome the limits of RANS models and allow for
correct a priori prediction of the separation point, but with very limited success [16].
Another option that promises to be the right balance between accuracy of results
and computational time, is to use LES techniques. Few good quality LES studies of
LSB caused by curvature induced APG are present in literature and usually a DNS-like
resolution is used. For instance Yang and Voke simulated a
at plate with a blunt leading
edge and used almost 2 million points for a Reynolds number of only 3450. For a case
similar to the present study Eisenbach and Friedrich [25] used around 50 and 100 million
points to perform LES simulations of
ow separation over an airfoil at a high angle of
attack. In this case a Cartesian immersed boundary method was used at Re
c
= 10
5
.
Another important point should be made. LES studies are not always compared to
experimental results or DNS and, even when they are, there is very rarely any validation
of skin friction and pressure coecient. LSB are extremely thin and make it very dicult
12
to measure skin friction coecients, thus it is very dicult to have a detailed quantitative
validation between experiments and LES simulations.
In the following sections a more detailed description of the three above mentioned
modeling approaches (DNS, LES and RANS) is given, described in increasing level of
approximation.
3.1 Direct Numerical Simulations
Conceptually the DNS modeling approach is the simplest, it consists of solving the full
Navier-Stokes equations and resolving all scales of motion. In DNS a
ow is distinguished
from another only by the equation of state, material properties, initial and boundary
conditions. This method is very expensive because when dealing with turbulence there
is the need to resolve a wildly separated range of length and time scales. The order of
magnitude of largest scales is the same as for the characteristic length of the domain
L, while the smallest scales decrease as Re
1=2
L
, for the Kolmogorov time-scales, and
Re
3=4
L
for the Kolmogorov length-scales [53]. Re
L
is the Reynolds number based on the
characteristic length of the domain and it is the ratio between the forces of inertia and
viscous forces. For a time-dependent 3-D turbulent simulation the computational work
required is therefore proportional to Re
1=2
L
Re
9=4
L
=Re
11=4
L
. For a turbulent
ow DNS the
solved 3-D and time dependent
ow eld is only one realization of the possible outcomes,
thus since in engineering applications we are interested in the mean quantities there is a
need for temporal averaging. In the aerospace industry the Reynolds numbers of interest
are in the rangeRe
L
= 10
5
10
8
. From the above mentioned arguments it is easy to see
13
how DNS simulations for moderate and high Reynolds numbers quickly become extremely
expensive to compute. This approach is not used in the current study, but it is used in
the reference study of Jones et al. [41].
As reference, the compressible Navier-Stokes in conservative form are
@
@t
+
@ (u
i
)
@x
i
= 0 (3.1)
@ (u
i
)
@t
+
@ (u
i
u
j
)
@x
j
=
@p
@x
i
+
@
ij
@x
j
(3.2)
@E
@t
+
@
@x
j
[(E +p)u
j
] =
@ (u
i

ij
)
@x
j

@q
j
@x
j
(3.3)
where u
i
are the components of the velocity vector, p the pressure, the density and E
the total energy per unit mass. The constitutive relation between stress and strain
rate for a Newtonian
uid is

ij
=

@u
i
@x
j
+
@u
j
@x
i

2
3
@u
k
@x
k

ij

(3.4)
the heat
ux q
i
is dened as
q
i
=k
@T
@x
i
(3.5)
where is the dynamic viscosity,k is the thermal conductivity, andT is the temperature.
The denition of total energy is
E =e +
1
2
u
i
u
i
(3.6)
14
where e is the specic energy per unit mass. Alternatively, instead of the total energy
other primary variables can be used, for instance temperature or enthalpy. In case the
temperature is used as primary variable, the thermal energy equation (with e = c
v
T ,
where c
v
is the specic heat capacity at constant volume) replaces equation (3.3)
@ (c
v
T )
@t
+
@ (u
j
c
v
T )
@x
j
=p
@u
j
@x
j
+
ij
@u
i
@x
j

@q
j
@x
j
(3.7)
In case the enthalpy is used as primary variable (with h = c
p
T , where c
p
is the specic
heat capacity at constant pressure), the total energy equation (3.3) is replaced by
@ (h)
@t
+
@ (u
j
h)
@x
j
=
@p
@t
+u
j
@p
@x
j
+
ij
@u
i
@x
j

@q
j
@x
j
(3.8)
3.2 Large Eddy Simulations
The LES modeling technique is based on the idea of solving only for the large 3-D
turbulent scales, while discarding and modeling the small scales in an attempt to reduce
the vast separation of scales typical of turbulent
ows thus reducing the computational
cost. This approach is justied by the fact that the large scales are the energy containing
scales and are the most aected by boundary conditions while the small scales are sup-
posed to be weakly dependent on the boundary conditions and therefore universal. Also
this technique, as DNS, requires averaging in time to get mean quantities of interest.
There are four conceptual steps in LES modeling: [53]
i Filtering Navier-Stokes equations
ii Derivation of LES equations with sub lter scales (SGS)
1
terms
15
iii Modeling of residual SGS stress tensor
iv Numerical solution of ltered (LES) equations
The ltering concept was introduced by Leonard [44]

(x
i
;t) =
Z
G(x
i
;x
0
i
;t)(x
0
i
;t)dx
0
i
(3.9)
where G is the lter function and for consistency it is required that
Z
G(x
i
;x
0
i
;t)dx
0
i
= 1 (3.10)
Often a convolution form is used for the lter [56]

(x
i
;t) =
Z
G(x
i
x
0
i
;t)(x
0
i
;t)dx
0
i
(3.11)
where G is the convolution kernel. The non-resolved eld is then dened as

0
(x
i
;t) =(x
i
;t)

(x
i
;t) (3.12)
There are a couple of points associated with the ltering idea that are noteworthy
a) even though formally the lter is only applied to spatial scales, it is indirectly also a
temporal lter since the dynamics of Navier-Stokes allows association of temporal scales
with spatial scales; b) there are actually three ranges of scales: represented resolved scales,
represented non-resolved (corrupted) scales and non represented scales; c) the eective
1
Even though the term sub lter scales is more appropriate because it is more general and does not
need to introduce the concept of a grid, in literature this term is usually referred to as sub grid scales.
16
lter that acts on a simulation is a combination of three lters: an analytical lter to
express the ltered Navier-Stokes equations, a lter associated with the computational
grid and a lter induced by the numerical scheme itself.
When dealing with compressible
ows, most authors opt for a change of variables
that simplies the derived LES equations, where the ltered variables are weighted by
the density [28]
~
=

(3.13)
this operation is called Favre ltering and does not apply to pressure and density elds.
Also it should be noted that Favre ltering does not commute with partial dierential
operators. This change of variables is a purely mathematical operation and it allows
for the unaltered form of the ltered continuity equation. The Favre ltered continuity
equation becomes
@
@t
+
@ ( ~ u
i
)
@x
i
= 0 (3.14)
To avoid the appearance of SGS terms in the equation of state another change of variables
is needed. Dierent approaches are suggested in literature; following is presented one
proposed by Vreman (System I) [67]. The idea is to solve for the resolved total energy

E which is dened as

E = ~ e +
1
2
~ u
i
~ u
i
(3.15)
17
One of the advantages of this systems is that it does not require any change in the
thermodynamic variables. Following this choice, the momentum equation becomes
@ ( ~ u
i
)
@t
+
@ ( ~ u
i
~ u
j
)
@x
j
=
@ p
@x
i
+
@~
ij
@x
j
+
@
@x
j
(
ij
+
ij
~
ij
) (3.16)
where
ij
is the SGS stress tensor

ij
= (g u
i
u
j
~ u
i
~ u
j
) (3.17)
and
ij
~
ij
comes from the change of variables due to Favre ltering and is generally
neglected, i.e. it is assumed
ij
= ~
ij
. The energy equation becomes
@

E
@t
+
@
@x
j
h

E + p

~ u
j
i
=
@~ q
j
@x
j
+
@ (~ u
i
~
ij
)
@x
j
B
1
B
2
B
3
+B
4
+B
5
B
6
(3.18)
where B
i
are the SGS terms and are dened as follows
B
1
=
1

1
@
@x
j
(pu
j
p~ u
j
) =
@
@x
j
(C
v
Q
j
) (3.19)
B
2
=p
@u
j
@x
j
p
@~ u
j
@x
j
=
dil
(3.20)
B
3
= ~ u
i
@
ij
@x
j
=
@~ u
i

ij
@x
j

ij
@~ u
i
@x
j
=B
3a
B
3b
(3.21)
B
4
=
ij
@u
i
@x
j

ij
@~ u
i
@x
j
(3.22)
B
5
=
@
@x
j
(
ij
~ u
i
~
ij
~ u
i
) (3.23)
18
B
6
=
@
@x
j
( q
j
~ q
j
) (3.24)
where Q
j
=

g
Tu
j

~
T ~ u
j

is the SGS heat
ux. Formally all the SGS terms B
i
should
be modeled. Vreman [66] did a thorough a priori investigation of all the SGS terms, i.e.
computed the SGS terms on a DNS eld and ltered them on a LES grid. The conclusion
of the study was that the terms B
3b
,B
4
,B
5
,B
6
are small or negligible and therefore are
not modeled. In practice the only modeled terms are the SGS stress tensor
ij
and the
SGS heat
ux Q
j
and sometimes the pressure dilatation term
dil
.
As shown in equations (3.14) -(3.18) the ltering operation produces extra SGS terms
in the Navier-Stokes equations. There are two fundamentally dierent approaches to
model these extra terms:
i the explicit LES approach consists of modeling these extra forcing terms through
a subgrid model. Furthermore we can distinguish between functional modeling and
structural modeling [56]. Functional modeling consists of modeling eects of the SGS
terms based on the ltered quantities, without trying to model the terms themselves.
It assumes a universal character of the small scales. Structural modeling aims to
reproduce the eigenvectors of the SGS terms; it requires knowledge of the structure
of small scales
ii the implicit LES add no extra terms to the governing equations, but the implicit
lters coming from the discretization scheme and the grid are tailored to act as SGS
model, i.e. the scheme is built in such a way that the numerical error acts as SGS
model.
19
Given the above denition of explicit and implicit LES modeling approaches and the
observation that the eective acting lter is always a combination of an analytical, dis-
cretization scheme and grid lters, it is interesting to note that in reality all schemes that
use explicit LES models are in fact hybrid LES. This observation is supported by the work
of Vreman [67] who observed that if the ratio of the lter width to grid spacing is around
unity, even for a 4th order scheme, the discretization error is of the same order as the
SGS terms. Further proof is given by the recent work of Cadiuex et al. [8] where a high
order scheme simulating a LSB gives closer agreement to DNS result in an under-resolved
DNS (uDNS) test case (which can be seen as untailored/non-calibrated implicit LES)
than in the test case where explicit (hybrid) LES was performed. This can be justied
by the fact that the numerical error is acting as additional, undesired SGS dissipation
corrupting the simulation. It can be further suggested that in some sense implicit LES is
better justied (especially in compressible
ows) since no terms are neglected; the eects
of all the SGS terms are globally modeled through the discretization scheme.
Most of the simulations done with INCA in the present study were carried with an
implicit LES scheme and some in uDNS conguration. For details about the method refer
to section 5.1 in the numerical method chapter, while for a discussion about the interac-
tion between the truncation error of the scheme and the SGS model refer to section 3.6.
Similarly control simulations done with Star-CCM+ were carried out testing all the ex-
plicit LES models available in the software and also in uDNS mode. In Star-CCM+ there
is the option to solve for the temperature equation (3.7) or enthalpy equation (3.8). For
the present simulations the temperature formulation has been used. The LES equations
solved by Star-CCM+ follow the derivation of Martin et al. [48, 54]. As most codes, in
20
the ltered temperature and enthalpy equations the only SGS term modeled is the SGS
heat
ux term.
3.3 Reynolds Averaged Navier-Stokes
The RANS modeling technique is by far the least computationally expensive, but
can be fairly involved, especially for compressible
ows. RANS modeling is based on a
statistical approach justied by the idea that turbulence consists of random
uctuations.
Reynolds (1895) proposed to express all quantities in the
ow eld as the sum of two
components, a mean and a
uctuating part
= +
0
(3.25)
where is the mean variable
=

(3.26)
where the bar symbol : represents a Reynolds average operator

=

and
0
repre-
sents the
uctuating component. Since the equations are solved only for the mean eld,
resolution requirements are not very demanding; in addition since only the mean
ow
motions are resolved there is no need to average multiple realizations in time. As we will
see later in this section, averaging of Navier-Stokes equations leads to RANS equations
and again to the so-called closure problem, i.e., the averaging procedure gives rise to an
additional term in the equations (Reynolds stress tensor) that needs to be modeled. First
we need to introduce some hypotheses that are needed to derive the compressible RANS
21
equations. In 1962 Morkovin observed that for wall bounded
ows the compressibility
eects on the turbulent eddies are small (except for hypersonic
ows and
ows with large
temperature gradients) therefore he suggested to assume that density
uctuations have
small eects on turbulence provided thatj
0
j<< [68]. In addition, to simplify the com-
pressible RANS equations and make them similar to the incompressible RANS equation
(similarly to the idea of Favre ltering expressed in equation (3.13) ), it is necessary to
introduce the Favre averaging procedure; this is a mass-averaged operations dened as
~ u
i
=
1

lim
4t!1
1
4t
Z
t
2
t
1
(x
k
;)u
i
(x
k
;)d (3.27)
thus
~ u
i
=u
i
= U
i
+
0
u
0
i
(3.28)
Now the instantaneous velocities u
i
, temperature T , specic internal energy e and
specic enthalpy h are decomposed with Favre Averaging, for instance
u
i
= ~ u
i
+u
00
i
(3.29)
while density, pressurep and heat
uxq
i
are decomposed with regular ensemble average.
Decomposition of the conservative form of Navier-Stokes leads to the compressible RANS
equations or Favre (mass) averaged mean conservation equations
@
@t
+
@ ( ~ u
i
)
@x
i
= 0 (3.30)
22
@ ( ~ u
i
)
@t
+
@ ( ~ u
i
~ u
j
)
@x
j
=
@ p
@x
i
+
@
@x
j

ij
u
00
i
u
00
j

(3.31)
@ ( E)
@t
+
@ ( ~ u
j
H)
@x
j
=
@
@x
j

q
j
u
00
j
h
00
+
ij
u
00
i

1
2
u
00
j
u
00
i
u
00
i

+
@
@x
j
h
~ u
i

ij
u
00
i
u
00
j
i
(3.32)
where E is the total energy per unit mass
E = ~ e +
1
2
~ u
i
~ u
i
+
1
2
u
00
i
u
00
i
(3.33)
and H the total enthalpy per unit mass
H =
~
h +
1
2
~ u
i
~ u
i
+
1
2
u
00
i
u
00
i
(3.34)
where h = e +p= is the specic enthalpy per unit mass. The Favre averaged conti-
nuity (3.14) and momentum equations (3.31) are identical to the laminar counterpart
except for an additional term, the so called Favre averaged Reynolds stressu
00
i
u
00
j
. This
additional term (and other terms in the energy equation (3.32) that will not be analyzed
for sake of brevity) needs to be modeled to close the system of equations (i.e. the closure
problem). For virtually all RANS models the Reynolds stress tensor is modeled based on
the Boussinesq hypothesis (for analogy with the stress tensor) and is dened as
u
00
i
u
00
j
=
T

@~ u
i
@x
j
+
@~ u
j
@x
i

2
3
@~ u
k
@x
k

ij

2
3
k
ij
(3.35)
23
where
T
is the turbulent (eddy) viscosity and k =
1
2
u
00
l
u
00
l

ij
is the turbulent kinetic
energy per unit volume. Several more or less sophisticated RANS models have been
developed (zero, one or two equation models) to specify the local value of
T
, but they
all share the Boussinesq approximation deciencies, i.e. the models will perform poorly
for
ows in which the assumption that the exchange of turbulent momentum behaves
similarly to the exchange of molecular momentum is incorrect. In particular RANS
models are expected to fail for
ows with sudden change in mean strain rate,
ows over
curved surfaces, rotating
ows, 3-D
ows and separated
ows. Because of the proven
inability of RANS to correctly capture LSB this modeling technique is not considered in
this study.
3.4 Wall modeling in Large Eddy Simulations
Whenever realistic geometries, at Reynolds numbers relevant to the aerospace in-
dustry, are modeled with LES the near wall resolution is the limiting factor; in this
context there is the need to distinguish between wall-modeled LES and wall-resolved
LES. Wall-bounded
ows, due to the presence of boundary layers, have strong resolution
requirements. There are two possible approaches, the boundary layer can be fully re-
solved or most of the boundary layer is modeled to reduce the numerical stiness of the
problem. In the rst case is referred to as wall-resolved LES while the latter is referred
to as wall-modeled LES. For the above reason there is a signicant dierence in mesh
size scaling with Reynolds number between the two approaches: wall-resolved LES grid
resolutions scale as Re
13=7
, while wall-modeled LES scales as Re [19]. Given the
24
above scaling laws, it is easy to show that for realistic geometries and moderate to high
Reynolds numbers (for instance an airliner at cruising speed) for now and in the near fu-
ture it is not feasible to use wall-resolved LES. On the other hand in the present study we
focus to moderate to low Reynolds numbers therefore we can aord to use wall-resolved
LES. Furthermore, if wall-modeled LES were to be used in this study it would become
extremely challenging to distinguish between the modeling errors coming from the wall
model and the ones coming from the LES model.
3.5 Types of numerical error
So far we have been talking about numerical dissipation in a rather abstract manner.
The numerical dissipation is not the only type of numerical error; together with the
dispersive error makes up the truncation error. The term truncation refers to the fact
that we use a truncated Taylor expansion to represent a derivative. For linear equations
we need only to worry about the truncation error, while in non-linear equations another
source of error may appear, namely the aliasing error. The aliasing error comes from
evaluation of a function at discrete points. When we calculate a non-linear term from a
discrete solution new high wave-numbers are generated. These high wave-numbers, if the
resolution is not sucient, have no modal support on the computational grid and appear
as aliased low wave-numbers polluting the numerical solution. The character of the error
that will be dominant strongly depends on the numerical scheme used. We are going to
render these concepts less abstract with few examples.
25
Consider the linear scalar advection equation
@u
@t
+
@u
@x
= 0; (3.36)
with periodic boundary conditions and appropriate initial condition, 0 x L and
t 0. The numerical derivative in space for a equally spaced grid of N points can be
calculated with dierent schemes, for instance central dierence
@u
i
@x

u
i+1
u
i1
24x
; (3.37)
or upwind
@u
i
@x

u
i
u
i1
4x
; (3.38)
or a spectral collocation method
@u
i
@x

N1
X
j=0
D
ij
u
j
; (3.39)
where D
ij
is the spatial derivative matrix dened as
D
ij
=
8
>
>
<
>
>
:

L
(1)
i+j
cot

(ij)
N

if i6=j;
0 if i =j:
(3.40)
These three spatial schemes have been chosen as they are representative of schemes with
truncation error predominantly dispersive (central dierence), diusive (upwind) and
schemes in which the truncation error is small (collocated spectral) and therefore the
26
dominant character of the error might be aliasing (when applied to non-linear equations).
This latter behavior is typical of high-order schemes. Since we want to limit our discussion
to the spatial error, in the following examples a high order explicit time integrator (the
classic 4th order Runge-Kutta) and a very small time step have been used.
The dierence in behavior among the schemes is striking (see Figure 3.1). In Fig-
ure 3.1(b) we see the `unpacking' of the initial condition due to dispersive error, i.e.
modes with higher wave-numbers have a larger dispersive error and travel slower than
the modes with lower wave-numbers. In Figure 3.1(d) we see the eect of numerical diu-
sion, in this case all wave-numbers travel at the right speed, but the amplitude is severely
diminished. In Figure 3.1(f) we see that neither dispersive error nor diusive error are
signicant and the exact solution is very well approximated. It seems that high-order
methods are the clear winners, on the other hand they are very dicult to implement
for complex geometries and often require ltering to be stable. Filtering is equivalent to
adding articial viscosity. In Figure 3.2 it is shown an example of a
ow eld corrupted
by aliasing error. Aliasing typically appears as an oscillation which wave-length is 2h,
whereh is the grid spacing, and slowly grows in time. Details about the implementation
and the analysis of the numerical error of the above scheme can be found in [26] for
the upwind and central nite dierence schemes and in [11] for the spectral collocation
method.
27
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: cd Time= 0.000

computed
exact
(a) central dierence, t = 0
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: cd Time= 6.000

computed
exact
(b) central dierence, t = 6
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: up Time= 0.000

computed
exact
(c) upwind dierence, t = 0
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: up Time= 6.000

computed
exact
(d) upwind dierence, t = 6
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: sp Time= 0.000

computed
exact
(e) collocated spectral, t = 0
0 2 4 6 8 10
−0.2
−0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
x
u
Scheme: sp Time= 6.000

computed
exact
(f) collocated spectral, t = 6
Figure3.1: Solution of the linear advection equation for dierent numerical schemes. Computed
solutions (blue line), exact solutions (red dots).
28
Figure 3.2: Pressure plot of a simulation exhibiting aliasing error, 2h spurious waves can be
recognized.
29
3.6 Interplay between numerical dissipation and SGS model
As introduced in section 3.2, in LES the number of degrees of freedom is reduced by
means of a spatial lter that suppresses the eects of small scales at the cost of introducing
sub-grid scale (SGS) unknowns (i.e. for the incompressible Navier-Stokes the SGS stress
tensor) which must be explicitly modeled [53, 56, 28].
An alternative approach is to use the numerical dissipation originating from the dis-
cretization of the Navier-Stokes equations as an implicit LES model. The strategy of using
the numerical dissipation as an implicit model relies on the hypothesis that the eect of
the SGS terms on the resolved scales is primarily dissipative [28]. This approach was
proposed by Boris et al. [7] who utilized a Flux-Corrected Transport (FCT) scheme. The
FCT method assures the monotonicity of the solution, which is why the original approach
was dubbed monotonically integrated LES (MILES). This approach has been expanded
to other monotonic and non-monotonic schemes and has been more generally renamed
as ILES. An example of a non-monotonic implicit LES is the compact nite dierence
scheme stabilized by lters [27]. The MILES approach has been controversial and as such
it has been the object of rigorous investigations [29, 22]. These studies have not been
particularly encouraging. Even when MILES appears to reproduce qualitatively the dy-
namics of Navier-Stokes equations, a more in-depth, quantitative investigation has shown
that this is not the case. Broadly speaking, these studies present two scenarios. Either
the numerical dissipation is excessive with respect to the correct SGS dissipation leading
to poor results both in implicit LES and explicit LES congurations [29], or the scheme is
under-dissipative (with respect to the correct SGS dissipation) leading to good results for
30
short time integrations and poor results for long time integrations due to accumulation of
energy in the high wave numbers [22]. The latter case can potentially be adjusted either
by ltering or with the addition of an explicit SGS model [63]. Both situations show
a fundamental deciency of the MILES approach: there is no mechanism embedded in
the numerics to ensure the correct amount of SGS dissipation. Nevertheless, the implicit
LES approach has been gaining more popularity due its simplicity and good qualitative
results for specic
ows, even if lacking a general justication or detailed comparisons
with DNS results. In implicit LES the numerics are rarely designed in such a way that the
numerical dissipation matches, a priori, the physical SGS dissipation. For a properly de-
signed implicit LES such a matching should be attempted at least for canonical turbulent

ows, e.g., isotropic turbulence or the turbulent channel
ow. Whenever the numerics
are not constrained to reproduce the correct amount of numerical dissipation, implicit
LES may not provide even basic quantities correctly such as the log-law of the wall or the
skin friction. An example of a proper implicit LES implementation is the adaptive local
deconvolution method (ALDM) of Hickel et al. [36], where the discretization is based
on a solution-adaptive deconvolution operator which allows for control of the truncation
error so that the numerical viscosity matches the values predicted for isotropic turbulence
by turbulence theories. Early explicit LES studies pointed out that low-order numerical
methods are not suitable in the explicit LES framework as the interaction between nu-
merical dissipation and the SGS dissipation [43] would negatively aect the results, while
for high-order/spectral methods the leading source of error is aliasing of the non-linear
term. Recently it has been shown that at very coarse resolutions even formally high-order
31
methods can suer from the interaction between numerical dissipation and SGS dissipa-
tion making the addition of an explicit LES model detrimental to the performance of the
code [9].
Despite the issues mentioned above,
ow simulations without explicit LES models are
becoming more and more popular. The reason for this trend and the attractiveness of
the approach are due to several reasons: its simplicity; a resulting qualitative behavior
that mimics the dynamics of Navier-Stokes equations; and the lack of a universal explicit
SGS model that could guarantee clearly superior results in all situations. Furthermore
implicit LES results are often validated with experimental results which may suer from
a high degree of uncertainty for many complex
ows of interest in engineering. Due to
the increasing popularity of the implicit LES approach it is prudent to make a formal
distinction between implicit LES schemes that are designed to provide the correct amount
of SGS dissipation and those that are simply run without an explicit model, with no
constraints imposed on the numerical dissipation. We will refer to implicit LES schemes
in the former case and to under-resolved DNS in the latter.
Recently, the need to assess the numerical dissipation in such simulations has been
addressed by [58] who developed the methodology that allows for the quantication of the
numerical dissipation for an arbitrary CFD code in a self consistent way, i.e., using only
information about the
ow eld provided by the code being analyzed, viz., the solver can
be treated as a black box without the need to know the details of the implementation. This
numerical dissipation quantication tool provides a rigorous method to judge a posteriori
the quality of a given simulation allowing for an impartial assessment of the impact of the
32
numerical dissipation and it will be used in Chapter 7 to analyze the results from Star-
CCM+ simulations. This work was presented in a paper [12] were this methodology was
applied for the rst time to a non-trivial
ow conguration solved using a commercial,
low-order compressible CFD code.
33
Chapter 4
Research goals
The main goal of this project is to assess the capabilities of the LES modeling tech-
nique and IB methods to accurately reproduce a LSB due to APG induced by geometry
curvature at a moderate Reynolds number with a grid resolution of the order of 1% with
respect to the required DNS resolution. The 1% resolution target is commonly achieved
for wall-resolved LES of fully turbulent
ows, for instance in channel
ows [63]. A LSB

ow is a particularly dicult test for LES models because it contains a laminar boundary
layer, separation, transition to turbulence of the separated shear layer, a turbulent non
equilibrium boundary layer and a turbulent equilibrium boundary layer. It has to be
emphasized that the challenge comes not only from the nature of the
ow but also from
the very low resolution goal. In literature when LES is used to simulate LSB usually
DNS-like resolution is used, but as explained in the subsection 3.6 at low resolution the
error of the discretization scheme can be of the same order as the SGS terms, making the
1% resolution goal a severe test for LES models. A secondary objective is to assess the
current capabilities of commercially available low-order code (Star-CCM+) to perform
34
coarse LES studies, this will be done by comparing the results with INCA and with the
DNS reference study of Jones et al. [40].
35
Chapter 5
Numerical methods
5.1 INCA
5.1.1 Governing equations
The numerical code solves the compressible three-dimensional Navier-Stokes equations
in conservative form
Z
V

@U
@t
+r [F (U) +G (U)]

dV = 0; (5.1)
where the vector U is dened as
U = [;u
1
;u
2
;u
3
;E]
T
; (5.2)
and the conserved variables are density , momentum u
i
and total energy E; u
i
are
the velocity components. The inviscid
ux tensor F is dened as
36
F
i
(U) =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
u
i
u
i
u
1
+
i1
p
u
i
u
2
+
i2
p
u
i
u
3
+
i3
p
u
i
E +u
i
p
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
; (5.3)
and G is the viscous
ux tensor
G
i
(U) =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
0

i1

i2

i3
u
k

ik
q
i
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
; (5.4)
where
ij
is the viscous stress tensor dened in equation (3.4), q
i
is the heat
ux as given
in equation (3.5) and the denition of total energy is
E =
1

1
p +
1
2
u
i
u
i
: (5.5)
The above equations are solved in non-dimensional form. The
uid is assumed to be
a perfect gas with Prandtl number of Pr = c
p
=k = 0:72 and a specic heats ratio of

= c
p
=c
v
= 1:4. The above equations are non-dimensionalized using the cord length
c and the free stream values: the velocity U
1
, the density
1
, and the temperature
T
1
. The pressure is non-dimensionalized using
1
U
2
1
. The governing non-dimensional
37

ow parameters are the Reynolds number Re =
1
U
1
c=
1
and the Mach number
Ma =U
1
=a
1
, where a
1
is the speed of sound based on T
1
and
1
is the viscosity in
the free stream. Pressure and temperature are then determined by the non-dimensional
ideal gas equation of state
p =
1

Ma
2
T: (5.6)
A non-dimensional power law is assumed for the temperature dependence of viscosity
=
1
Re
T
0:75
; (5.7)
and the non-dimensional thermal conductivity is expressed by
k =
1
(
1)Ma
2
Pr
: (5.8)
5.1.2 Grid
The generation of suitable grids for LES of complex
ows can be time-consuming and
dicult. Contradictory requirements, such as adequate local resolution and minimum
number of grid points, can deteriorate the grid quality and adversely aect accuracy and
numerical convergence properties. In this case Cartesian grids are used which facilitate
grid generation. Cartesian grids also imply fewer computational operations per grid point
than body tted or unstructured grids. On the other hand, geometric boundaries do not
necessarily coincide with grid lines, so that boundary conditions have to be applied at
the subcell level. INCA implements a conservative immersed boundary (IB) method
38
for representing sharp interfaces between a
uid and a rigid body on Cartesian grids,
see Refs. [49, 32]. A level-set technique is used for describing the interface geometry.
The level-set function is the signed distance between each point in the computational
domain and the immersed interface, which is positive within and negative outside of the

uid domain. The zero contour is the interface between the
uid and the obstacle. The
intersection of the obstacle with the Cartesian grid produces a set of cells that are cut
by the interface. The underlying FV discretization is modied locally in these cut cells
in such a way that it accounts only for the
uidic part of a cell. Required interface
normals, face apertures and FV fractions of cut cells can be computed eciently from
the level-set function. The viscous stresses at the
uid-solid interface are approximated
by linear dierencing schemes. The interface pressure is obtained by solving the one-
sided (symmetric) Riemann problem in the interface-normal direction. As the scheme
operates on
uxes only, this cut-cell FV method maintains accuracy and ensures mass
and momentum conservation. Discrete conservation and a sharp representation of the

uid-solid interface render this method particularly suitable for LES of turbulent
ows.
5.1.3 Finite Volume Discretization
The Navier-Stokes equations are numerically solved on a grid that is too coarse to
resolve all scales representing the turbulent
uid motion. However, the unavailable small-
scale information is crucial for the proper evolution of large-scale structures in the
ow.
Therefore, the eect on the resolved scales of their nonlinear interactions with the unre-
solved SGS has to be represented by a turbulence model. SGS models generally operate
on
ow scales that are marginally resolved by the underlying discretization scheme, i.e.,
39
the truncation error of common approximations for the convective terms can outweigh
the eect of even physically sound models. Representing a full merger of discretization
and SGS model, the Adaptive Local Deconvolution Method (ALDM) has shown to be
a reliable, accurate, and ecient method for implicit LES of Navier-Stokes turbulence
[36, 37]. The implicit SGS model provided by this discretization can be interpreted as
a combination of eddy-viscosity and scale similarity modeling. Model parameters were
determined by a spectral-space analysis of the eective eddy viscosity in isotropic tur-
bulence in the Reynolds numbers innite limit by constraining the numerical dissipation
to the physical SGS dissipation obtained from the analysis of nonlinear interactions in
turbulence. ALDM is used for discretizing the hyperbolic
ux (equation (5.3)), whose
discretization causes the SGS eects of interest, whereas standard centered dierences are
used for discretizing the viscous
ux (equation (5.4)). For comparative purposes other
classic discretization schemes for the inviscid
uxes have been used in under-resolved DNS
mode, i.e., without any explicit or implicit LES model. The other schemes used are 3rd
and 5th order WENO with HLLC
uxes and 4th and 6th order central dierence methods.
These schemes were used because they are known to be either more dissipative (WENO
schemes) or less dissipative (central dierence schemes) than the ALDM discretization.
For time integration, the conditionally stable 3rd order Runge-Kutta (R-K) method of
Gottlieb and Shu [31] is used. This time-discretization scheme is total-variation diminish-
ing (TVD) for CFL 1, provided the underlying spatial discretization is TVD, whereas
the linear stability bound is larger.
As cut cell methods can generate cells with a very small
uid volume fraction, a
special treatment of such cells is necessary when excessively small time steps according to
40
=α
i,j
Figure 5.1: Volume fraction of
uid for the celli;j. The light shaded area represents the
uid
domain.
the CFL stability criterion are to be avoided. To increase the computational eciency,
the code uses a conservative mixing procedure, where the conservative variables of small
cells are mixed with larger neighbors based on the volume fraction, i.e. if the volume
of small cells is smaller than a mixing threshold (see Figure 5.1) [49, 32]. During the
simulations, the time step size is adjusted dynamically according to CFL = 0:5 based on
the full cell size.
5.1.4 Boundary Conditions
In general high order methods simulating compressible
ows require a way to control
of wave re
ection from the boundaries. In particular, in order to limit the size of the
computational domain there is the need for non re
ective boundary conditions to avoid
41
unphysical spurious wave re
ections that could strongly in
uence the entire
ow-eld.
When a steady state solution can be obtained, the behavior of the boundaries could
be ignored. In the present study this is not the case: there will be turbulent unsteady
structures crossing the out
ow boundary. The method used is the Navier-Stokes charac-
teristic boundary conditions (NSCBC) of Poinsot and Lele [52, 47]. The NSCBC method
is based on characteristic analysis [65]. The boundary conditions in a computation can
be divided in physical boundary conditions and numerical boundary conditions. Know-
ing the physical boundary conditions is not sucient to solve the numerical problem.
Whenever the number of physical boundary conditions at the boundary is less than the
number of primitive variables, extra numerical boundary conditions are required. The
idea behind NSCBS is to x only the physical boundary conditions while calculating the
characteristic wave propagation of the other quantities. In the method implemented, the
value of the wave amplitudes in the multidimensional viscous case is inferred by exam-
ining the local associated one-dimensional (LODI) case. Afterwards these approximated
wave amplitudes are used in the full system of conservation equations to calculate the
primitive variables at the boundary for the numerical boundary conditions. Outgoing
waves are calculated with one-sided stencils, while ingoing waves are approximated using
LODI relations to expresses the unknown amplitudes (ingoing) as functions of the known
amplitudes (outgoing). One-dimensional NSCBC formulation is exact only for the 1-D
case, for 2-D and 3-D cases only waves normal to the boundary are correctly transmitted.
For non-re
ecting subsonic out
ow a modication of the general NSCBC formulation is
needed to avoid a drift of the mean pressure in the computational domain. To solve this
42
problem a numerical condition imposing the mean pressure at innity is added. The
ingoing wave amplitude is prescribed
L =K (pp
1
); (5.9)
where p is the local pressure, p
1
is the target mean pressure at innity and K is a
constant calculated following the relation given by Rudy and Strikwerda [55]
K =
a
L

1M
2
max

; (5.10)
where is a constant, that will be called the re
ection parameter, a is the local speed of
sound,L is a domain reference length andM
max
is the expected maximum Mach number
in the domain. NSCBC are applied at the far eld domain boundaries where the only
disturbances likely to reach the boundaries will be in the form of linear waves. At the
downstream out
ow boundary, which will be subject to the passage of nonlinear
uid
structures, a low value for the re
ection parameter has been set.
5.2 Star-CCM+
Star-CCM+ solves the compressible three-dimensional integral Navier-Stokes equa-
tions in conservative form (equation (5.1)). The equations are solved in a preconditioned
dimensional form on an unstructured grid [17]. The time is advanced through a dual
time-stepping implicit scheme. The inviscid
uxes (equation (5.3)) are evaluated by
43
using the Weiss-Smith preconditioned Roes
ux-dierence splitting scheme which is for-
mally at best second order accurate. The viscous
uxes (equation (5.4)) are evaluated by
a standard central dierence scheme. Two dierent explicit models are available in Star:
dynamic Smagorisky [30, 45] and wall-adapting local eddy viscosity model (WALE)[50].
Both have been used in the present work.
5.2.1 explicit LES
Both WALE and Smagorinsky are functional models based on the eddy-viscosity or
Boussinesq hypothesis, i.e. the energy transfer between the small resolved eddies and the
SGS eddies is thought to be analogous to the transfer of energy in molecular viscosity.
With this hypothesis the SGS terms found in section 3.2 are modeled as follows.
The SGS stress tensor is decomposed in deviatoric and isotropic parts

ij
=
d
ij
+
1
3

kk
: (5.11)
For incompressible and low Mach number compressible
ows the isotropic part of the
SGS stress tensor is negligible and it is lumped together with the resolved pressure term
through a change of variable
P = ~ p +
1
3

kk
(5.12)
and we model the deviatoric part of the SGS stress tensor as

d
ij
=
sgs

@~ u
i
@x
j
+
@~ u
j
@x
i

2
3
@~ u
k
@x
k

ij

= 2
sgs

~
S
ij

1
3
~
S
kk

ij

(5.13)
44
where S
ij
=
1
2

@u
i
@x
j
+
@u
j
@x
i

is the strain rate tensor. The SGS heat
ux is
Q
j
=

sgs
Pr
sgs
@
~
T
@x
j
(5.14)
where
sgs
is a subgrid scale viscosity that needs to be calculated. How
sgs
is calculated
is what dierentiates explicit functional models from one another. The SGS Prandtl
number Pr
sgs
is usually kept as a constant Pr
sgs
= 0:6, or can be calculated through a
dynamic procedure. [28]
5.2.2 Dynamic Smagorinsky model
The simplest version of the Smagorinsky model prescribes the SGS viscosity as

sgs
=C
2
s
4
2
j
~
S
ij
j (5.15)
wherej
~
S
ij
j =

2
~
S
ij
~
S
ij

1=2
, C
s
is a constant and4 is a length scale associated with the
lter cut-o and is usually dened as4 =V
1=3
whereV is the volume of the considered
cell. In the dynamic model following the work of Germano et al. [30] and the subsequent
modication suggested by Lilly [45] the constant C
s
is adjusted dynamically for every
time step for each point in space. The idea is to use a second test lter and that the SGS
tensors at both levels can be modeled by the same constant. The test ltering procedure
for an unstructured grid is dened as
^

=
1
P
N
n=0
V
n
N
X
n=0

n
V
n
(5.16)
45
where 0 is the considered cell andN the total number of neighbors. The test ltered SGS
stress tensor is
^
T
ij
=
\
u
i
u
j

1
^

c
u
i
c
u
j
(5.17)
the dierence between the two SGS stress tensors is
L
ij
=
^
T
ij
^
ij
=
\
~ u
i
~ u
j

1
^

c
~ u
i
c
~ u
j
: (5.18)
Assuming that both SGS stress tensors can be modeled by the same constant (scale-
invariance assumption invoked by Germano) then we can write

ij

1
3

kk

ij
=C
2
s

ij
(5.19)
^
T
ij

1
3
^
T
kk

ij
=C
2
s

ij
(5.20)
where and are model-dependent kernels for the dynamic procedure

ij
=2
^
4
2
^
j
^
~
S
ij
j

^
~
S
ij

1
3
^
~
S
ij

ij

(5.21)

ij
=24
2
j
~
S
ij
j

~
S
ij

1
3
~
S
ij

ij

(5.22)
where
^
4 is the test lter width.
Following least square method to determine the constant proposed by Lilly
C
2
s
=
hM
ij
L
ij
i
hM
ij
M
ij
i
(5.23)
46
where M
ij
=
ij

^

ij
andhi usually represents an averaging procedure in the direction
of statistical homogeneity. Since for an unstructured grid this averaging direction cannot
be dened, the averaging is done among all the neighboring cells.
5.2.3 WALE model
The WALE model proposed by Nicoud and Ducros [50] is one of the few SGS models
to exhibit the correct asymptotic behavior of
SGS
close to the wall without the need for
wall functions [28]. The SGS viscosity is dened as

sgs
=C
2
w
4
2

S
d
ij
S
d
ij

3=2

~
S
ij
~
S
ij

5=2
+

S
d
ij
S
d
ij

5=4
(5.24)
with
S
d
ij
=
1
2

@~ u
i
@x
l
@~ u
l
@x
j
+
@~ u
j
@x
l
@~ u
l
@x
i

1
3
@~ u
m
@x
l
@~ u
l
@x
m

ij
(5.25)
and in Star-CCM+ the WALE constant is set at C
w
= 0:544.
47
Chapter 6
Simulation results
6.1 Benchmark
ow conguration
The specic geometrical setting used as a benchmark
ow is that of a NACA-0012
three-dimensional airfoil atRe
c
= 510
4
at 5

of incidence and Mach numberMa = 0:4.
For the same conguration detailed DNS results were obtained by [40, 41] using a fourth
order central dierence scheme. The scheme used entropy splitting of the inviscid
ux
terms to increase stability but neither articial viscosity nor ltering were necessary.
In our simulations the computational domain is a rectangular prism with the following
dimensions: x = [6; 6], y = [7; 7] and z = [0; 0:2], where x is the streamwise direction
and theyz plane is normal to the incoming
ow direction, with z along the wingspan.
The airfoil trailing edge is located at (x;y) = (1; 0).
48
6.2 Computational setup
6.2.1 INCA
The geometry of the NACA-0012 airfoil at 5

of incidence and the computational
grid are shown in Figure 6.1. The computational domain is a rectangular prism with the
following dimensions: x = [6; 6], y = [7; 7] and z = [0; 0:2], where x is the streamwise
direction and the yz plane is normal to the incoming
ow direction, with z along the
wingspan. The airfoil trailing edge is located at (x;y) = (1; 0). For the simulations with
INCA several two-dimensional Cartesian grids have been generated. Here we present only
three dierent resolutions which will be referred to as \coarse",\ne" and \very ne" grid.
The coarse grid has around 240; 000 cells, the ne grid has around 1:6 million cells and the
very ne grid has around 2:3 million cells. For comparison, the baseline two-dimensional
case of Jones et al. [40] used 1:6 million grid points.
Due to the diculties of the IB method to handle very sharp edges a blunt trailing
edge has been adopted. The thickness of the blunt trailing edge is less 0.1% of the chord
while the thickness of the turbulent boundary layer at the trailing edge is around 10%
of the cord, it is thus expected that the presence of the blunt trailing edge will have
negligible eects on the
ow-eld. This local modication is common practice and in
similar studies found in literature has been shown to have negligible in
uence [21].
The meshes are built to satisfy the following requirements at the airfoil boundary:
4n
+
1 and4tg
+
< 25, where4n
+
is the cell size in the wall normal direction and
4tg
+
is the cell size in the tangential direction, both evaluated at the boundary in non-
dimensional wall units (4x
+
= u

4x=, where u

is the wall-shear velocity). Due to
49
x
y
0 0.2 0.4 0.6 0.8 1
0.1
0
0.1
0.2
Figure 6.1: Grid around the airfoil, only every 5th line is shown (INCA, coarse grid).
limitations in the possible stretching of the cells, in practice4tg
+
25. It should be
noted that for the ne grid, even though the total number of cells is the same as in
Jones, the actual distribution of cells/points in the domain is quite dierent due to the
distinct numerical approaches, i.e., Cartesian grid with IB versus generalized curvilinear
coordinates.
To create three-dimensional LES grids the two-dimensional grids are extended in the
spanwise direction using 32 cells within the boundary layer, close to the airfoil and in
the wake, reducing up to 4 cells away from the airfoil where the
ow is laminar to reduce
the computational eort. Same technique has been used to generate a three-dimensional
DNS grid, in this case the number of spanwise cells used is 96. The three-dimensional
coarse grid has around 7:2 million cells, the ne grid has around 46:2 million cells and the
DNS grid has 138:7 million cells. For comparison, the three-dimensional grid of [40] used
170:7 million grid points. For each run, starting from a converged solution, the
ow has
been simulated for around 12 time units,c=U
1
, wherec is the airfoil cord length andU
1
the free stream velocity. The equations of motion have been non-dimensionalized using
50
x
y
0.115 0.12 0.125 0.13 0.135
0.12
0.125
0.13
0.135
Figure 6.2: Immersed boundary grid close-up (coarse grid).
c and U
1
, therefore in the results shown c = 1 and U
1
= 1. Because of the low Mach
number, the time step in the coarse grid had to be set to a low value of t = 0:2 10
4
,
which is a factor of 5 smaller than used by [40]. An even smaller time step is used for the
ner grids.
6.2.2 Star-CCM+
The same geometry as for INCA has been used in simulations performed with the
commercial code Star-CCM+ , with the exception of the sharp trailing edge, which could
be retained. The same domain size as for INCA has been used. In order to quickly
generate the grid, the unstructured mesh generator within Star-CCM+ has been used.
To have a grid similar to the grid used in INCA a hexahedral grid for the outer
ow is used
51
Figure 6.3: Grid around the airfoil (Star-CCM+ ).
with a boundary tted grid around the airfoil (Figure 6.3). The same constraints of4n
+
and4tg
+
as for the grid for INCA have been imposed. In gures 6.3 and 6.4 can be seen
the freestream hexahedrals and the transition to the body tted grid. The total number
of cells used in two-dimensional cases is around 54; 000. For the three-dimensional cases
the mesh is simply extended in a homogeneous fashion in the spanwise direction with 30
cells. Since Star-CMM+ uses an implicit time stepping, a time step needs to be chosen.
According to the work of [18] for a turbulent channel
ow the maximum time step in wall
units (4t
+
=4tu
2

=) that allows for accurate results is4t
+
max
= 0:4; a value below
4t
+
max
has been used for all the simulations.
6.3 Numerical simulations
These simulations had been initiated during the rst SFB-TRR40 summer program
[13] by performing two-dimensional low resolution simulations using INCA. During the
second SFB-TRR40 summer program [15] several runs with increasing resolution were
performed to obtain the best match with the benchmark data. The challenges of obtaining
52
X
Y
0.05 0.1 0.15 0.2 0.25
0.12
0.14
0.16
0.18
0.2
Figure 6.4: Grid transition between hexahedral cells in the freestream and body tted cells
around the airfoil.
good results at very low resolutions come partially from the diculties in using Cartesian
grids stretched in the streamwise direction at curved boundaries. Additionally, two-
dimensional simulations with the ne grid and four dierent discretization schemes have
been performed; in particular we have chosen two schemes that are more dissipative than
ALDM (3rd and 5th order WENO schemes) and two that are less dissipative (4th and 6th
order central dierence schemes). Several two-dimensional runs with the ALDM scheme
and with increasing resolution have been made. It should be noted that, by design,
ALDM provides a consistent SGS turbulence model for three-dimensional turbulence. In
two-dimensional or laminar
ows ALDM does not provide a turbulence model but merely
acts as a slightly dissipative, 2nd order accurate centered discretization. Therefore, the
two-dimensional results presented should be considered as essentially DNS results, while
the three-dimensional results should be considered as an implicit LES. Three-dimensional
simulations for the coarse and ne grid have been carried out only with the ALDM and
the 5th order WENO schemes.
53
With Star-CCM+ the original goal of a reduced grid resolution has been met, so there
was no need for studying ner grid resolutions. Six simulations have been performed, three
two-dimensional cases and the respective three-dimensional cases. An under-resolved
DNS (without an LES model) simulation and two simulations with explicit LES models
available in the code: WALE and Dynamic Smagorinsky model. All the cases presented
here are summarized in Table 6.1 for the two-dimensional simulations and in Table 6.2
for the three-dimensional simulations.
6.4 Evaluation metrics
As a qualitative comparison the instantaneous spanwise vorticity eld of the simulated

ow is shown in Figure 6.5 which can be compared with the results of Jones et al. [40].
The present two-dimensional and three-dimensional simulations capture well qualitative
features observed in high resolution DNS, including the presence of a separated
ow. In
particular the characteristics of the laminar separation bubble described in Ch. 2 can be
observed: laminar boundary layer, separated shear layer, break up of the shear layer,
transition to turbulence, and turbulent reattachment. Only the vorticity elds for the
coarse grid cases with the ALDM scheme are shown because all the other simulations
look qualitatively similar.
Similarly we can observe the same features in Figure 6.6 where the three-dimensionality
of the problem can be appreciated, here it is shown an iso-surface of the Q-criterion col-
ored with x-velocity. The Q-criterion is the second invariant of the velocity gradient
tensor and it is used to isolate vortical structures in the
ow [24].
54
(a) Isocontours of vorticity for two-dimensional coarse grid simulation (INCA, ALDM).
(b) Isocontours of vorticity for three-dimensional coarse grid simulation (INCA, ALDM).
Figure 6.5: The spanwise vorticity eld (10 levels over the range75 !
+
). [14]
55
Figure 6.6: Iso-surface of Q-criterion (Q
+
= 50) colored with x-velocity (Star-CCM+, UDNS).
As a principal quantitative metrics to evaluate our simulations we use quantities that
are of paramount importance for the design of wings or blades: pressure coecient
C
p
=
pp
1
1
2
u
2
1
(6.1)
and skin friction coecient
C
f
=

w
1
2
u
2
1
(6.2)
where P
1
is the free stream pressure, U
1
is the free stream velocity, p is the pressure
at the surface and
w
is the wall shear stress. The pressure and friction coecients
from our simulations are compared with the corresponding data of [40] for their high
resolution DNS case. All quantities are averaged in time for at least 10 time units. Also,
for a comparison among dierent simulations the location of the separation point, x
s
, the
56
reattachment point, x
r
, and the bubble length, l
b
, are listed in Table 6.1 and 6.2 for the
two-dimensional and three-dimensional cases, respectively. This quantitative information
about the laminar separation bubble is extracted from the analysis of the skin friction
coecient plots (see Figs. 6.7(b), 6.9(b), 6.10(b), and 6.11(b)). In particular, at the
suction side the point where the skin friction becomes negative denes x
s
and where the
skin friction becomes positive again denes x
r
.
6.5 Two-dimensional results
6.5.1 INCA
For simulations performed with INCA the comparison for the pressure coecient
(Figure 6.7(a)) is quite good for all cases, especially on the pressure side of the airfoil;
this is expected since on the pressure side the
ow is fully laminar. The prediction of
C
p
on the suction side is not fully satisfactory for the cases with the ALDM scheme:
the plateau formed due to the presence of the laminar separation bubble is captured,
but the secondary peak at the end of the plateau is not. Substantially increasing the
resolution, up to the benchmark DNS resolution, improves the results. However, no
perfect agreement with the DNS reference is attained, and the dierences between the
ne grid and the very ne gird results are negligible. The same conclusions can be drawn
from the skin friction results (see Figure 6.7(b)). There is overall agreement with the
DNS reference, especially at the pressure side, and at the suction side all simulations
accurately predict both the separation and the reattachment point. However, for the
ALDM cases the magnitude and features of the secondary peak in C
f
are not correctly
57
captured even with the very ne grid. After a preliminary three-dimensional run (coarse
grid with ALDM, see Figure 6.10) it was observed that the results were closer to the
forced case of Jones et al. than to the unforced one. The forcing was added by Jones
et al. as a means to promote faster transition to turbulence in the three-dimensional
simulations. This observation suggested that perhaps the IB procedure was generating
some numerical noise at the airfoil boundary and that this noise was acting as a forcing
promoting early transition to turbulence. In order to explore this idea additional runs
with the ne grid were performed with schemes that are known to be more dissipative
than ALDM, hoping that the additional numerical dissipation would damp out the noise
originating from the IB method. The rst run was made with a 3rd order WENO scheme.
This brought an improvement in the results, notably the secondary peak after the bubble
was now resolved, but the added numerical dissipation was excessive resulting in a shorter
bubble length. The second run with a 5th order WENO scheme (denoted as WENO5 in
the following) provided the right amount of numerical dissipation and all the features of
the laminar separation bubble were fully captured as can be seen in Figure 6.7. In an
attempt to further clarify the role of numerical dissipation in the dierence between the
results obtained with WENO5 and ALDM, two additional runs with the ne grid were
made. This time 4th and 6th order central dierence schemes were used. Both runs went
unstable after a little more than one time unit, suggesting that these are under-resolved
DNS or that the noise coming from the IB fully corrupted the simulation after growing
unbounded. Finally one additional run at higher resolution with the ALDM scheme (very
ne grid) was made, but the dierences with the ne grid results are negligible.
58
Flow visualizations show an apparently identical
ow evolution for the WENO5 and
ALDM simulations. In order to identify and isolate the dierences that lead to the
discrepancy in the time averaged C
f
and C
p
distributions, we performed a Dynamic
Mode Decomposition (DMD) for a sub-domain 0:3 x 0:85 that fully encloses the
region where
ow-separation and vortex shedding takes place. DMD is a technique that
allows for a modal analysis of any data sequence [57]. It is used here to identify modes
and frequencies that express the dominant dynamic behavior captured in the simulated
time series. Figure 6.8 shows the resulting frequencies and corresponding amplitudes.
The dierence between the two simulations is striking. For the WENO5 simulation,
the dynamic behavior is fully governed by a single low-frequency dynamic mode with
Strouhal number St = 3:36 and its higher harmonics. The Strouhal number is dened
as St = f=U
1
where f is the frequency. For the ALDM time series, DMD identies
two dominant low-frequency modes with St = 3:0 and St = 2:78, and a broad spectrum
of less energetic modes that also contribute to the
ow evolution. All relevant dynamic
modes are very well resolved on the computational grid and have reached a saturated
state. We conclude that the pronounced secondary peaks of C
f
and C
p
in the WENO5
results are due to vortex shedding taking place at a constant streamwise coordinate and
with constant frequency, which results in a single dynamic mode and its harmonics. The
ALDM solution on the other hand is essentially dominated by two dynamic modes with
slightly dierent frequencies. Though the shed vortices have the same intensity as in the
WENO5 run, this leads to a slight upstream-downstream oscillation of the vortex shedding
location and a more smeared C
f
and C
p
distribution in the temporal average. We were
unable to clearly identify the origin of the observed discrepancies between these two cases.
59
Possible explanations are that the additional dynamic modes result from
ow instabilities
with small growth rates or from disturbances due to the geometry representation, which
are suppressed if a more dissipative numerical method is used.
6.5.2 Star-CCM+
To speed up the initialization of the
ow a multigrid approach has been used, where
results of a coarse simulation were interpolated on a ner grid. Two steps with a very
coarse and a coarse grid, run with a RANS model have been performed to obtain an
initial solution for the nal LES grid. Quantitatively, all Star-CCM+ simulations behave
similarly showing a respectable performance on both the pressure and the suction side.
In particular, the pressure coecient plot (Figure 6.9(a)) shows that the peak at the end
of the plateau on the suction side is well captured in all cases. Similarly, the skin friction
plot (Figure 6.9(b)) shows that the secondary peak at mid-cord on the suction side is
well captured, even by the under-resolved DNS case. The prediction of separation and
reattachment points is excellent in all cases. The only signicant dierence with the DNS
reference case is the wavy pattern on the suction side after reattachment. This dierence
could be explained by high frequency ltering due to the coarse grid and/or the much
bigger time step allowed by the implicit time stepping.
60
x
Cp
0 0.2 0.4 0.6 0.8
1.5
1
0.5
0
0.5
1
Jones et al.
Fine WENO 5
Fine ALDM
Coarse ALDM
Very Fine ALDM
(a) Pressure coecient comparison (INCA)
x
Cf
0 0.2 0.4 0.6 0.8
0.01
0
0.01
0.02
0.03
Jones et al.
Fine WENO 5
Fine ALDM
Coarse ALDM
Very Fine ALDM
(b) Friction coecient comparison (INCA)
Figure 6.7: Two-dimensional Pressure and friction coecients at the wall. DNS of Jones et al.
[40] (black dots), 5th order WENO scheme with ne grid (solid red line), ALDM with ne grid
(dashed dotted blue line), ALDM with coarse grid (dashed green line), and ALDM with very ne
grid (long dashed black line). [14]
61
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0 5 10 15 20 25
ALDM
WENO5
Figure 6.8: Results of the dynamic mode decomposition for WENO5 and ALDM simulations.
[14]
62
x
Cp
0 0.2 0.4 0.6 0.8
1.5
1
0.5
0
0.5
1
Jones et al.
UDNS
Dyn. Smag.
WALE
(a) Pressure coecient comparison (Star-CCM+ )
x
Cf
0 0.2 0.4 0.6 0.8
0.01
0
0.01
0.02
0.03
Jones et al.
UDNS
Dyn. Smag.
WALE
(b) Friction coecient comparison (Star-CCM+ )
Figure 6.9: Two-dimensional Pressure and friction coecients at the wall. DNS of Jones et al.
[40] (black dots), UDNS (solid red line), explicit LES with Dynamic Smagorinsky model (dashed
dotted blue line), and explicit LES with WALE model (dashed green line). [14]
63
6.6 Three-dimensional results
6.6.1 INCA
Four three-dimensional simulations have been performed using INCA. Two were per-
formed with the ALDM model (coarse and ne grid) and two with the 5th order WENO
scheme (ne grid with 32 and 96 cells in the spanwise direction). It should be noted that
the simulations with the ne grid have to be considered as LES because even though the
resolution in thexy plane is DNS-like, the resolution in the spanwise direction (32 cells)
is typical of a coarse LES. The 5th order WENO scheme is known to have high levels
of the numerical dissipation [62]. The high numerical dissipation acts as an uncontrolled
implicit SGS model and if the simulation were to be performed with an explicit SGS
model the unknown numerical dissipation would contaminate the eects of the model
[10]. However the best agreement for the two-dimensional simulations was obtained with
this scheme on the ne grid and for this reason it has been used in a three-dimensional
case as well. The comparison for the pressure coecient (Figure 6.10(a)) is fairly good
overall, although the bubble for both ALDM cases is visibly shorter compared to the DNS
reference. Dierently, for the WENO5 case the bubble is longer. More sensitive is the
skin friction comparison shown in Figure 6.10(b). Similarly to the two-dimensional case
the pressure side is well reproduced. On the suction side, the separation point is fairly
well captured. The reattachment point varies substantially between the coarse and ne
ALDM cases and the bubble in both cases is signicantly shorter than in the benchmark
DNS. We may speculate that the noise from the immersed boundary procedure promotes
early transition to turbulence and consequently early reattachment. With the ner grid
64
the forcing is weaker and the IB solution approximates better the reference DNS data.
For the WENO5 case the separation point is better predicted, but the reattachment oc-
curs too late suggesting that the uncontrolled numerical dissipation provided by WENO
scheme is slightly too low to behave as a correct implicit SGS model. To support this
claim an additional simulation with the full DNS resolution also in the spanwise direction
has been run. As predicted this allow for a good match of the reattachment point with
the DNS of Jones et al.
6.6.2 Star-CCM+
For Star-CCM+ all two-dimensional cases were ran also in three dimensions. As
in the two-dimensional case the C
p
and C
f
plots agree well with the DNS of Jones et
al. The agreement of the pressure coecient (Figure 6.11(a)) is excellent and there
is no signicant dierence between the under-resolved DNS case and the two explicit
LES cases. Looking at the skin friction plot (Figure 6.11(b)) it can be seen that all
cases slightly overestimate the mid-cord peak on the suction side, but separation and
reattachment points are correctly predicted. The good performance of the under-resolved
DNS case and the fact that there is no substantial dierence with the explicit LES cases
are interesting. This behavior suggests that the numerical dissipation provided by the
scheme is dominant therefore the additional SGS dissipation added by the explicit SGS
models must be negligible. This issue will be addressed in dept in a chapter 7.
65
x
Cp
0 0.2 0.4 0.6 0.8 1
1.5
1
0.5
0
0.5
1
Jones et al.
Fine WENO 5
Fine WENO 5 96
Fine ALDM
Coarse ALDM
(a) Pressure coecient comparison (INCA)
x
Cf
0 0.2 0.4 0.6 0.8 1
0.01
0
0.01
0.02
0.03
Jones et al.
Fine WENO 5
Fine WENO 5 96
Fine ALDM
Coarse ALDM
(b) Friction coecient comparison (INCA)
Figure 6.10: Three-dimensional Pressure and friction coecients at the wall. DNS of Jones et
al. [40] (black dots), UDNS with 5th order WENO scheme with ne grid (solid red line), DNS
with 5th order WENO scheme (short dashed black line), ALDM with ne grid (dashed dotted
blue line), and ALDM with coarse grid (long dashed green line). [14]
66
x
Cp
0 0.2 0.4 0.6 0.8 1
1.5
1
0.5
0
0.5
1
Jones et al.
UDNS
Dyn. Smag.
WALE
(a) Pressure coecient comparison (Star-CCM+ )
x
Cf
0 0.2 0.4 0.6 0.8 1
0.01
0
0.01
0.02
0.03
Jones et al.
UDNS
Dyn. Smag.
WALE
(b) Friction coecient comparison (Star-CCM+ )
Figure 6.11: Three-dimensional Pressure and friction coecients at the wall. DNS of Jones
et al. [40] (black dots), UDNS (solid red line), explicit LES with Dynamic Smagorinsky model
(dashed dotted blue line), and explicit LES with WALE model (dashed green line). [14]
67
Simulation Model N 10
5
x
s
x
r
l
b
Jones et al. [40] DNS 16 0.151 0.582 0.431
INCA (coarse) ALDM 2.4 0.177 0.583 0.406
INCA (ne) ALDM 16 0.164 0.585 0.421
INCA (very ne) ALDM 23 0.171 0.597 0.426
INCA (ne) UDNS (WENO 3) 16 0.129 0.562 0.433
INCA (ne) UDNS (WENO 5) 16 0.155 0.586 0.431
INCA (ne) UDNS (CD 4) 16 - - -
INCA (ne) UDNS (CD 6) 16 - - -
Star-CCM+ UDNS 0.54 0.142 0.590 0.448
Star-CCM+ WALE 0.54 0.151 0.581 0.430
Star-CCM+ Dyn. Smag. 0.54 0.153 0.581 0.428
Table6.1: Summary of two-dimensional time averaged laminar separation bubble characteristics.
N is the number of grid points, x
s
the separation point, x
r
the reattachment point and l
b
the
length of the bubble. [14]
Simulation Model N 10
6
x
s
x
r
l
b
Jones et al.[40] DNS 170.7 0.099 0.607 0.508
Jones et al. forced [40] DNS 170.7 0.128 0.500 0.372
Almutairi et al.[2] MTS 7.6 0.114 0.604 0.490
INCA (coarse, 32 pts) ALDM 7.2 0.125 0.462 0.337
INCA (ne, 32 pts) ALDM 46.2 0.138 0.549 0.411
INCA (ne, 32 pts) UDNS (WENO 5) 46.2 0.107 0.641 0.534
INCA (ne, 96 pts) DNS (WENO 5) 138.7 0.112 0.623 0.511
Star-CCM+ UDNS 1.6 0.090 0.601 0.511
Star-CCM+ Dyn. Smag. 1.6 0.090 0.598 0.508
Star-CCM+ WALE 1.6 0.090 0.596 0.506
Table 6.2: Summary of three-dimensional time averaged laminar separation bubble characteris-
tics. N is the number of grid points, x
s
the separation point, x
r
the reattachment point and l
b
the length of the bubble. [14]
68
6.7 Two-dimensional sweep of angles of attack
Given the good performance of Star-CCM+ for the 5

angle of attack case a sweep
of angles of attacks for the two-dimensional case has been performed. The same meshing
criteria and the same domain has been used to create meshes for the following angles of
attack: 0

, 1

, 1:5

, 2

, 3

and 4

. The rationale behind this sweep was to investigate the
capabilities of LES to capture the appearance of the LSB. Only a two-dimensional sweep
has been performed to reduce the computational cost. To avoid excessive line cluttering in
the plots the results have been plotted in two dierent sets, the rst set contain the angles
from 0

to 2

and the second from 3

to 5

. The skin friction and pressure coecients plots
can be seen in Figure 6.12 and Figure 6.13. Additionally separation point, reattachment
point and bubble length have been calculated for each case and collected in Table 6.3.
In Table 6.3 also two additional two-dimensional cases from the PhD thesis of Jones are
reported [39]. The low angles of attack cases, namely from 0

to 3

have been run for
40 non-dimensional time units because in 10 time units the averaged skin friction and
pressure coecients were not fully converged. This is due to the low frequency shedding
happening at these low angles of attack.
There is good quantitative agreement with the results of Jones and the typical trend
for LSB is captured: the reattachment point moves upstream with increasing angle of
attack; similarly the reattachment point moves upstream with increasing angle of attack.
However, the length of the bubble has a maximum, around 3

, after which the reattach-
ment points moves upstream faster than the separation points leading to a shorter and
shorter bubble.
69
Simulation Angle of attack N 10
5
x
s
x
r
l
b
Star-CCM+ (uDNS) 0

0.54 0.656
Jones et al.(DNS) [39] 0

16 0.658
Star-CCM+ (uDNS) 1

0.54 0.482
Star-CCM+ (uDNS) 1:5

0.54 0.470 0.977 0.507
Star-CCM+ (uDNS) 2

0.54 0.441 0.960 0.519
Star-CCM+ (uDNS) 3

0.54 0.311 0.831 0.520
Jones et al.(DNS) [39] 3

16 0.315 0.850 0.535
Star-CCM+ (uDNS) 4

0.54 0.215 0.708 0.493
Star-CCM+ (uDNS) 5

0.54 0.142 0.590 0.448
Jones et al.(DNS) [39] 5

16 0.151 0.582 0.431
Table6.3: Summary of two-dimensional time averaged laminar separation bubble characteristics
for the sweep of angles of attack. N is the number of grid points, x
s
the separation point, x
r
the
reattachment point and l
b
the length of the bubble.
Integral lift and drag coecients have been constructed via integration the skin friction
and pressure coecient curves at each angle of attack. Comparing the resulting polar
plot (Figure 6.14) and the experimental results from Spedding and McArthur (Figure2.2)
we can observed the same qualitative trend. Note that due to the symmetric nature of
our airfoil only positive angles have been considered. A jump in lift and drag coecients
can be observed in the numerical polar plot as well as in the experimental one. For the
NACA-0012 atRe
c
= 510
5
the jump happens between 1:5

and 3

. A close inspection
of the skin friction coecient plot for the 1:5

case (Figure 6.12(b)) reveals that at this
angle of attack a stable LSB is already formed, while at 1:0

there is only a partial,
unstable reattachment.
70
x
Cp
0 0.2 0.4 0.6 0.8
1
0.5
0
0.5
1
0 deg
1 deg
1.5 deg
2 deg
(a) Pressure coecient comparison (Star-CCM+ )
x
Cf
0 0.2 0.4 0.6 0.8
0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
0 deg
1 deg
1.5 deg
2 deg
(b) Friction coecient comparison (Star-CCM+ )
Figure 6.12: Two-dimensional Pressure and friction coecients at the wall for dierent angles
of attack (Star-CCM+, uDNS).
71
x
Cp
0 0.2 0.4 0.6 0.8
1.5
1
0.5
0
0.5
1
3 deg
4 deg
5 deg
(a) Pressure coecient comparison (Star-CCM+ )
x
Cf
0 0.2 0.4 0.6 0.8
0.01
0.005
0
0.005
0.01
0.015
0.02
0.025
0.03
3 deg
4 deg
5 deg
(b) Friction coecient comparison (Star-CCM+ )
Figure 6.13: Two-dimensional Pressure and friction coecients at the wall for dierent angles
of attack (Star-CCM+, uDNS).
72
0 0.005 0.01 0.015 0.02 0.025 0.03
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
C
d
C
l
Figure 6.14: Two-dimensional Lift and Drag polar
73
Chapter 7
Numerical dissipation
The results from chapter 6 presented in subsections 6.5.2 and 6.5.2 show that adding an
explicit model has negligible in
uence on our gures of merit (skin friction and pressure
coecients). This is indicative of the presence of substantial numerical dissipation in
the commercial code. The recent work of Schranner et al. [58] is adapted to quantify
the numerical dissipation in Star-CCM+. This chapter is organized as follow: rst the
conceptual reasoning behind the method is outlined, then a procedure to applied the
method to Star-CCM+ is presented. Subsequently the procedure is tested for a known

ow (Taylor Green Vortex
ow) and nally it is applied to the simulations presented in
chapter 6.
7.1 Concept
The method allows for computation of the numerical dissipation rate and numerical
viscosity in the physical space for arbitrary sub-domains in a self-consistent way using
only information provided by the code in question. The method is inspired by the work
74
of Domaradzki and Radhakrishnan [23] who developed a way to quantify numerical dissi-
pation in a given code with the external help of a spectral solver. The limit of the original
procedure is that an external spectral code is needed. Schranner et al. [58] generalized
the original concept to obtain a self contained method that does not rely on an external
spectral code.
Consider the kinetic energy evolution equation applied to an exact solution of Navier-
Stokes equations.
@e
kin
@t
+N (;u
i
;p) = 0; (7.1)
where e
kin
=
1
2
u
i
u
i
;
@
@t
is the analytic time derivative andN is the analytic operator of
the kinetic energy evolution equation. Now consider the kinetic energy evolution equation
using the analytical operators and applied a posteriori to a general numerical time-space
solution generated with a given numerical scheme
@e
kin
(n)
@t
+N (
(n)
;u
i(n)
;p
(n)
) ="
n
(n)
; (7.2)
where the subscript [ ]
(n)
refers to the numerical solution at a given time level; "
n
is
the residual. It is important to stress that the analytical solution and its numerical
counterpart are not equal, i.e. (;u
i
;p)6= (
(n)
;u
i(n)
;p
(n)
). On the right-hand side
of Eq. (7.2) there is a residual because the constraint represented by Eq. (7.1) is not
enforced by the numerical solver directly. In other words, the kinetic energy evolution
equation is not being solved for numerically, therefore in general the sum of the terms
75
on the left-hand side will not be equal to zero. An estimate of the residual of the kinetic
energy evolution equation is obtained by approximating the exact operators as follows
"
n
(n)

@e
kin
(n)
@t

CD
+N (
(n)
;u
i(n)
;p
(n)
)

HO
; (7.3)
where [ ]

CD
represents a central dierence operator in time and [ ]

HO
represents the
highest-order available spatial operator. In principle the estimation of the second term
in Eq. (7.3) depends on the discretization order of the various terms and
uxes in the
balance, in practice Schranner et al. [58], using as a reference a spectral-Fourier code and
the original procedure of Domaradzki and Radhakrishnan [23], show that the estimate is
quite robust and fairly independent of the discretization applied to compute the balance.
7.2 Equations and method
7.2.1 Analytical form
The transport energy equation for the compressible Navier-Stokes equations is
@E
@t
+
@
@x
j
[(E +p)u
j
] =
@u
i

ij
@x
j

@q
j
@x
j
; (7.4)
where u
i
are the components of the velocity vector, p the pressure, the density and E
the total energy per unit mass. The constitutive relation between stress and strain rate
for a Newtonian
uid was given in Eq. (3.4). The denition of the total energy per unit
mass is
E =e +e
kin
; (7.5)
76
where e is the specic internal energy per unit mass and e
kin
=
1
2
u
i
u
i
is the kinetic
energy per unit mass. Following the procedure of Schranner et al. [58] the transport
equation (3.3) can be separated into the internal energy equation
@e
@t
+
@eu
j
@x
j
=p
@u
j
@x
j
+
ij
@u
i
@x
j

@q
j
@x
j
; (7.6)
and the kinetic energy equation
@e
kin
@t
+
@e
kin
u
j
@x
j
=u
j
@p
@x
j
+u
i
@
ij
@x
j
: (7.7)
The viscous term in equation (7.7) can be rewritten as
u
i
@
ij
@x
j
=
@u
i

ij
@x
j

ij
@u
i
@x
j
; (7.8)
where the rst term on the right hand side represents the viscous work. Thus equa-
tion (7.7) becomes
@e
kin
@t
+
@e
kin
u
j
@x
j
=u
j
@p
@x
j
+
@u
i

ij
@x
j

ij
@u
i
@x
j
: (7.9)
For nite volume methods (FV) where the integral form of the Navier-Stokes equations
is solved, it is convenient to rewrite the pressure term as
u
j
@p
@x
j
=
@u
j
p
@x
j
p
@u
j
@x
j
: (7.10)
77
Therefore the integral form of equation (7.9) can be written as
Z
V

@
@t
e
kin
+
@e
kin
u
j
@x
j
+
@u
j
p
@x
j
p
@u
j
@x
j

@u
i

ij
@x
j
+
ij
@u
i
@x
j

dV =
@
@t
Z
V
e
kin
dV +
Z
A

e
kin
u
j
+u
j
pu
i

ij

n
j
dA +
Z
V

p
@u
j
@x
j
+
ij
@u
i
@x
j

dV =
@
@t
E
kin
+F
kin
+F
ac
F
vis
W
p
+"
vis
= 0; (7.11)
whereF
kin
,F
ac
,F
vis
are the kinetic energy, acoustic and viscous
uxes, n
j
the outward
unit vector normal to the bounding surface A,W
p
the work due to pressure and "
vis
the
viscous dissipation. We dene a dissipation function " as
" =
Z
V
1

ij
@u
i
@x
j
dV
=
Z
V

@u
i
@x
j
@u
i
@x
j
+
@u
j
@x
i
@u
i
@x
j

2
3
@u
k
@x
k
@u
l
@x
l

dV; (7.12)
and the equation
" ="
vis
; (7.13)
denes the averaged kinematic viscosity for the volume V .
7.2.2 Discretized form
If we assume a FV spatial discretization and a generic discretization in time then the
discretized equation (7.11) is contaminated by the truncation errors and a local residual
"
n
(m)
=
E
kin
(m)
t
+F
kin
(m)
+F
ac
(m)
F
vis
(m)
W
p
(m)
+"
vis
(m)
(7.14)
78
does not vanish. In the above equation the subscript [ ]
(m)
refers to the m
th
control
volume. We call the residual"
n
the numerical dissipation rate because it has been shown
that if integrated over a suciently large control volume it has a predominantly dissipative
character [58, 23]. Equation (7.13) allows for the denition of a numerical kinematic
viscosity for the m
th
control volume as

n
(m)
=
"
n
(m)
"
(m)
: (7.15)
The above denitions can be extended to a larger sub-domain or to the entire computa-
tional domain
"
n
sub
=
M
X
m
"
n
(m)
;
n
sub
=
"
n
sub
"
sub
; (7.16)
"
n
tot
=
N
X
m
"
n
(m)
;
n
tot
=
"
n
tot
"
tot
; (7.17)
where M is the total number of adjacent cells of a given sub-domain and N is the total
number of cells of the simulation domain.
7.3 Procedure
This is a posteriori procedure where all quantities are computed using the numerical

ow eld. In principle this procedure does not require any information about the solver, all
that is needed are temporal snapshots of the conserved variables (the density, components
of the momentum/velocity, and the energy) and the corresponding mesh. In the present
work, except for the time discretization of the kinetic energy, all quantities appearing
in equation (7.14) have been calculated on-the-
y using the post-processing capabilities
79
of the code itself. We can rely on the low-order post-processing tools embedded in the
software since Schranner et al. [58] showed that even with a low-order reconstruction
one can obtain fair estimates of the terms in the kinetic energy balance. For the time
discretization a second-order three-points nite dierence formula is used, i.e.
E
kin
(m)
t
=aE
kin
(t t

)
(m)
+bE
kin
(t)
(m)
+cE
kin
(t + t
+
)
(m)
; (7.18)
where the weights a, b and c are
a =
t
+
t

1
t

+ t
+
; (7.19)
b =ac; (7.20)
c =
t

t
+
1
t

+ t
+
: (7.21)
The kinetic energy, acoustic and viscous
uxes are calculated as follows
F
kin
(m)
=

X
r
(e
kin
u
j
n
j
)
(r)
A
(r)
!
(m)
; (7.22)
F
ac
(m)
=

X
r
(pu
j
n
j
)
(r)
A
(r)
!
(m)
; (7.23)
F
vis
(m)
=

X
r
(u
i

ij
n
j
)
(r)
A
(r)
!
(m)
; (7.24)
80
where subscript [ ]
(r)
represents a sum over the r
th
face of the m
th
control volume and
A
(r)
is the area of ther
th
face. Volume terms such as kinetic energy, pressure work and
viscous dissipation rate are calculated in the following way
E
kin
(m)
=

1
2
u
i
u
i

(m)
V
(m)
; (7.25)
W
p
(m)
=

p
@u
i
@x
i

(m)
V
(m)
; (7.26)
"
vis
(m)
=
(m)

@u
i
@x
j
@u
i
@x
j
+
@u
j
@x
i
@u
i
@x
j

2
3
@u
k
@x
k
@u
l
@x
l

(m)
V
(m)
; (7.27)
where V
(m)
is the volume of the m
th
control volume. These steps have been rst vali-
dated for Star-CCM+ by repeating the same test cases as in [58] for both incompressible
and compressible solvers. Even though in this study a specic commercial software has
been used, it is important to stress that the solver is treated as a black box and that all
that is required to perform this analysis are temporal snapshots of the
ow elds.
7.4 Taylor Green vortex
ow
7.4.1 Numerical method
Star-CCM+ solves either the incompressible or compressible three-dimensional inte-
gral Navier-Stokes equations in conservative form. The equations are solved in a precon-
ditioned dimensional form on an unstructured grid [17].The time is advanced through a
81
dual time-stepping implicit scheme. The incompressible solver uses a Rhie-Chow pressure-
velocity coupling and a SIMPLE algorithm, while the for the compressible solver the in-
viscid
uxes are evaluated by using the Weiss-Smith preconditioned Roe's
ux-dierence
splitting scheme. Both schemes are formally at best second order accurate. The viscous

uxes are evaluated by a standard central dierence scheme.
7.4.2 Problem setup
This is the same case which was analyzed previously to validate the procedure applied
to INCA with the ALDM method as an implicit LES model [58, 36]. The Taylor Green
Vortex (TGV) problem is originally described by Taylor et al. [64] The initial conditions,
following [59], are
u =u
0
[sin(x)cos(y)cos(z)]; (7.28)
v =u
0
[cos(x)sin(y)cos(z)]; (7.29)
w = 0; (7.30)
=
0
; (7.31)
p =p
0
+

0
16
[(cos(2x) +cos(2y))(cos(2z) + 2) 2]; (7.32)
where subscript [ ]
0
indicate reference quantities. The initial
ow eld is let to evolve
for 10 non-dimensional time units. The time unit is T = l
0
=u
0
where l
0
= 1 is the
reference length. The time step is kept constant as t = T=50. Star-CCM+ solves the
dimensional form of the Navier-Stokes. Under ideal gas assumption, to obtain the same
initial condition as in Schranner et al. [58] and Hickel et al. [36] we need to pick the
82
(a) T = 0 (b) T = 10
Figure 7.1: Isosurface of Q-criterion (Q = 0) colored with velocity magnitude at initial and nal
time, inviscid case.
following reference values: u
0
=Ma
p

p
0
=
0
and
0
=p
0
=(RT
0
). The reference pressure
is p
0
= 101325 Pa, the Mach number Ma = 0:0845, the heat capacity ratio
= 1:4, the
reference temperature T
0
= 300 K and the specic gas constant R = 287:02 J=(kg K).
Even though the code solves the dimensional Navier-Stokes equations, all results will
be presented in non-dimensional form. Unless stated otherwise the resolution for all
simulations is of 64 cells per spatial direction for a cube with sides that are 2 long
and that has periodic boundary conditions. Initially all simulations are run with the
incompressible solver for a direct comparison with the work of Schranner et al. [58],
afterwards selected cases are run with the compressible solver.
83
0 2 4 6 8 10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
t
ε
n

64
128
INCA
× 10
−2
Figure 7.2: Time-evolution of numerical dissipation rate "
n
tot
, Re =1. Star-CCM+ with 64
3
cells (blue line), Star-CCM+ with 128
3
cells (red dashed line) and INCA from [58] (black circles).
7.4.3 Full periodic domain
For the rst test case we use as a control volume the whole domain so that we can
calculate the total numerical dissipation rate without
ux terms. For Re =1 we show
"
n
tot
plotted against time in Figure 7.2. The numerical dissipation rate is quite high,
it is around three times as much as in Schranner et al. [58] throughout the temporal
domain. To make sure that the "
n
tot
is calculated properly this case has been repeated
with 128 cells per spatial direction. This case can also be seen in Figure 7.2. We can
observe that for the laminar part of the simulation ( t< 3:5 ) there is indeed fairly good
agreement with the reference results. Enstrophy (
= 1=2
R
V
!
i
!
i
dV , where !
i
are the
component of the vorticity vector) is often used in TGV simulations as a metric to judge
the degree of resolution. It is hypothesized that for the inviscid TGV problem enstrophy
should increase to innity at the non-dimensional time t 5 [59]. The resolution of
this singularity is a common metric to evaluate a given scheme/resolution for a TGV
simulation. In Figure 7.3 we can see that for our current simulation we fall in between
the behavior of a spectral scheme, which has very little numerical dissipation and a
84
0 1 2 3 4 5 6 7
5
10
15
20
t
Ω/Ω(0)

64
128
WENO
spectral
Figure 7.3: Time-evolution of normalized enstrophy
(t)=
(0), Re =1. Star-CCM+ with 64
3
cells (blue line), with 128
3
cells (red dashed line), WENO (black squares) and spectral (black
circles) from [59].
fth-order WENO scheme which has substantial numerical dissipation, with Star-CCM+
closer to the latter. This is expected since the current scheme is formally at best second
order, therefore quite dissipative. The result for the 128
3
case is also included. While
this is not a fair comparison neither with the spectral nor the WENO schemes used
with the 64
3
resolution, it is still a useful way to show that the numerical dissipation
of the current scheme approaches the spectral numerical dissipation when the resolution
increases. As further proof of the signicant amount of numerical dissipation, if we plot
the time-dependent normalized kinetic energy in log-log scale (see Figure 7.4) we can
observe that for late times an energy cascade is established and we recover the kinetic
energy decay law E
kin
t
1:2
even if we are solving the inviscid equation.
For the same conguration, as in the reference paper, a sweep of Reynolds number
has been ran (Re =1, 3000, 1600, 800, 400, 200, and 100). We can see "
n
tot
plotted
against time in Figure 7.5 for the nite Reynolds number cases. The general trend is the
same as in the innite Reynolds number case, and the peak is about 2:5 times the peak
in the reference study. In Figure 7.6 we show a full sweep of Reynolds numbers to better
85
10
−1
10
0
10
1
10
−0.2
10
−0.1
10
0
t
E/E(0)

64
128
∝ t
−1.2
Figure7.4: Time-evolution of normalized kinetic energyE
kin
(t)=E
kin
(0),Re =1. Star-CCM+
with 64
3
cells (blue line), with 128
3
cells (red dashed line) and t
1:2
law (black dots).
appreciate the behavior of the numerical dissipation rate as a function of the Reynolds
number.
7.4.4 Sub-domains
As test cases we choose some of the sub-domains used in [58]. The rst cubic sub-
domain used in the reference was an octant of the domain. Due to symmetries of the
TGV problem this case is only useful to check that all
uxes cancel out. The second test
case simulated is the cubic sub-domain one (CSD1) which consists of cells 10 through
56 in all spatial directions. The results for the innite Reynolds number case are shown
in Figure 7.7. Here and in the following Figures terms of the kinetic energy balance are
plotted with a negative sign in order to sum up to "
n
. All the current results are plotted
against the results from the reference paper and also plotted are the various terms that
contribute to the numerical dissipation rate. We can see that we recover the same trend
observed in the innite Reynolds number case with the full periodic box, i.e. for the
laminar part of the simulation the results agree with the spectral reference while the
peak in the turbulent part of the simulation is about 3 times as high. In Figure 7.7
86
0 2 4 6 8 10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
t
ε
n

current
INCA
× 10
−2
(a) Re = 3000
0 2 4 6 8 10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
t
ε
n

current
INCA
× 10
−2
(b) Re = 1600
0 2 4 6 8 10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
t
ε
n

current
INCA
× 10
−2
(c) Re = 800
0 2 4 6 8 10
0.0
0.5
1.0
1.5
t
ε
n

current
INCA
× 10
−2
(d) Re = 400
0 2 4 6 8 10
0.0
0.5
1.0
1.5
t
ε
n

current
INCA
× 10
−2
(e) Re = 200
0 2 4 6 8 10
0.0
0.5
1.0
1.5
t
ε
n

current
INCA
× 10
−2
(f) Re = 100
Figure 7.5: Time-evolution of numerical dissipation rate "
n
tot
, niteRe. Star-CCM+ (blue line)
and INCA from [58] (black circles).
87
0 2 4 6 8 10
0.0
0.5
1.0
1.5
2.0
2.5
3.0
t
ε
n

∞
3000
1600
800
400
200
100
× 10
−2
Figure 7.6: Time-evolution of numerical dissipation rate "
n
tot
. Re =1 (blue line), Re = 3000
(red dashed line),Re = 1600 (black dash-dotted line),Re = 800 (magenta dotted line),Re = 400
(solid blue line with markers), Re = 200 (red dashed line with markers), and Re = 100 (black
dash-dotted line with markers) .
we can see a noisy oscillatory behavior in the acoustic
ux term and the work due to
pressure. These oscillations are out of phase and cancel each other out. We speculate
that this is an unphysical energy exchange between the two terms due to the fact that
incompressibility is not strictly enforced. For the nite Re cases we see the same trend
observed for the periodic box cases. A sweep of Re has been performed and the results
can be seen in Figure 7.8. To better visualize the additional dissipation provided by the
scheme with respect to the physical dissipation we have plotted the molecular viscosity
together with the total eective viscosity +
n
, i.e. a sum of the numerical and the
physical viscosity, see Figure 7.9.
The only non-cubic sub-domain considered here is the NCSD2 from reference [58],
i.e. for the x-direction cells 8 to 23 and for the y and z-direction cells 1 to 31. This
sub-domain is not chosen randomly, it was selected because this is the region within the
octant with non-negligible numerical dissipation within an octant. The results of the
inviscid case for the sub-domain NCSD2 are shown in Figure 7.10.
88
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ǫ
n
ΔE
kin
/Δt
F
kin
F
ac
W
p
INCA
× 10
−2
Figure 7.7: Time-evolution of numerical dissipation rate "
n
sub
, Re =1, sub-volume CSD1.
Numerical dissipation rate "
n
sub
(blue line), rate of change of kinetic energyE
kin
sub
=t (black
dot-dashed line), kinetic energy
uxF
kin
sub
(blue dot-dashed line with markers), acoustic
ux
F
ac
sub
(magenta dotted line) and INCA from [58] (black circles).
7.4.5 Surface versus Volume Integrals
For certain incompressible solvers pressure can be dened up to a constant without
changing the dynamics of the
ow eld since only the gradient of pressure appears directly
in the momentum equations. For this reason simply adding W
p
sub
to the energy balance
can potentially lead to wrong estimates of the numerical dissipation rate. One way to
check for consistency of the results is to recast the pressure terms asF
ac
(m)
W
p
(m)
=pt
(m)
, where pt
(m)
is the numerical form of the pressure transport term dened as pt
(m)
=

u
j
@p
@x
j

(m)
V
(m)
. As can be seen comparing Figures 7.7 and 7.11 the two methods are
equivalent. Additionally, other surface terms have been calculated as volume terms with
similar conclusions; there seem to be negligible dierences in calculating the quantities
in the kinetic energy balance as either volume or surface integrals.
89
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(a) Re = 3000
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(b) Re = 1600
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(c) Re = 800
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(d) Re = 400
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(e) Re = 200
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
F
vis
ε
vis
W
p
INCA
× 10
−2
(f) Re = 100
Figure 7.8: Time-evolution of numerical dissipation rate "
n
sub
, nite Re, sub-volume CSD1.
Numerical dissipation rate "
n
sub
(blue line), rate of change of kinetic energyE
kin
sub
=t (black
dot-dashed line), kinetic energy
uxF
kin
sub
(blue dot-dashed line with markers), acoustic
ux
F
ac
sub
(magenta solid line with markers), viscous
ux F
vis
sub
(black dashed line with markers)
viscous dissipation rate"
vis
sub
(red dashed line) and INCA from [58] (black circles).
90
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(a) Re = 3000
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(b) Re = 1600
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(c) Re = 800
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(d) Re = 400
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(e) Re = 200
0 2 4 6 8 10
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
t
ν

ν + ν
n
ν
× 10
−2
(f) Re = 100
Figure 7.9: Time-evolution of physical and numerical viscosity, nite Re, sub-volume CSD1.
Total viscosity +
n
(blue line), physical viscosity (red dashed line).
91
0 2 4 6 8 10
−2.0
−1.0
0.0
1.0
2.0
3.0
4.0
t
ε
n

ε
n
ΔE
kin
/Δt
F
kin
F
ac
W
p
INCA
× 10
−3
Figure 7.10: Time-evolution of numerical dissipation rate "
n
sub
, Re =1, sub-volume NCSD2,
including work due to pressure. Numerical dissipation rate "
n
sub
(blue line), rate of change of
kinetic energyE
kin
sub
=t (black dot-dashed line), kinetic energy
uxF
kin
sub
(blue dot-dashed
line with markers), acoustic
uxF
ac
sub
(magenta solid line with markers), W
p
sub
(red dashed line
with markers) and INCA from [58] (black circles).
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ǫ
n
ΔE
kin
/Δt
F
kin
pt
INCA
× 10
−2
Figure 7.11: Time-evolution of numerical dissipation rate"
n
sub
,Re =1, sub-volume CSD1 with
pressure gradient transport term. Numerical dissipation rate "
n
sub
(blue line), rate of change of
kinetic energyE
kin
sub
=t (black dot-dashed line), kinetic energy
uxF
kin
sub
(red dashed line),
pressure transportpt
sub
(magenta solid line with markers) and INCA from [58] (black circles).
92
7.4.6 Eects of Compressibility
All simulations so far have been performed with an incompressible solver to have a
one to one comparison with the results available in literature. Due to convergence issues
for the compressible case the time step has been reduced to t =T=500. For the CSD1
sub-domain all Reynolds number cases have been run for a complete comparison with
the incompressible simulations. Results for the nite Re cases can be seen in Figure 7.12
where the physical and the numerical dissipation rates are plotted. We nd a good
agreement with the incompressible solver at all Reynolds numbers (see Figure 7.8).
7.5 Discussion of the validation case
Even lacking a direct validation through an external spectral code as in the reference
paper of [58] the methods applied to a low-order incompressible and compressible schemes
compares well to the available reference data and it is self consistent. It should be noted
that for low-order schemes the work due to pressure, even though should be zero, might
not be negligible. The extension to a compressible solver, at least at low Mach numbers,
is straightforward. From our results it seems to not matter if terms in the discrete kinetic
energy balance, equation (7.14), are calculated as
uxes or as volume terms. The method
seems self consistent and robust as all our results are in qualitative agreement with the
results of [58]. Therefore we conclude that the method works well for low-order schemes
and that it is ready to be applied to more realistic
ow geometries.
93
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(a) Re = 3000
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(b) Re = 1600
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(c) Re = 800
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(d) Re = 400
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(e) Re = 200
0 2 4 6 8 10
−1.0
−0.5
0.0
0.5
1.0
1.5
t
ε
n

ε
n
ε
vis
INCA
× 10
−2
(f) Re = 100
Figure 7.12: Time-evolution of numerical dissipation rate "
n
sub
for a compressible solver,
Ma = 0:0845, nite Re, sub-volume CSD1. Numerical dissipation rate "
n
sub
(blue line), viscous
dissipation rate"
vis
sub
(red dashed line) and INCA from [58] (black circles).
94
SD1a SD3
SD1
SD2 SD2a SD3a
Figure 7.13: Sub-domains location for NACA-0012.
7.6 NACA-0012
In Chapter 6 it was suggested that there was signicant numerical dissipation in
simulations performed with Star-ccm+. Furthermore, we observed that adding an explicit
LES model have negligible in
uence on the results. This can be explained by having the
numerical dissipation being dominant over the SGS dissipation provided by the explicit
models. In the rst half of this Chapter we have presented and tested a method to
quantify the numerical dissipation, applying this method to the NACA-0012 simulations
should allow for the quantication of the above statement.
7.6.1 Results
In the present Section the procedure to quantify the numerical dissipation rate and
viscosity is applied only to the UDNS case from Section 6.6.2. Initially three sub-domains
have been selected to be analyzed, one in the laminar region (SD1), one in the transitional
region (SD2) and one in the fully turbulent region (SD3). The location and size of the
sub-domains can be seen in Figure 7.13. They all have approximatively the same number
95
of cells ( 20; 000) and they all extend for the entire domain length in the spanwise
direction. Starting from a converged statistically steady state solution the simulation has
been run for an additional time unit. The results for the numerical and viscous dissipation
rate are shown in Figures 7.14(a)(c)(e) while the physical and total viscosity are shown
in Figure 7.15(a)(c)(e). All three sub-domains behave as expected. The sub-domain
SD1 encompasses a region where the boundary layer is fully laminar, therefore the
ow
is well resolved and the added numerical dissipation is very small, almost negligible (see
Figure 7.15(a) and from now on refer to Table 7.1 for a summary of time averaged values).
In SD1 the numerical viscosity is around 4% of the physical viscosity. Sub-domain SD2
is in the transitional region of the
ow where there is a signicant amount of turbulence.
Here consequently the numerical dissipation is substantial (see Figure 7.15(c)), with the
numerical dissipation being more than 1:5 times the physical one. Lastly sub-domain
SD3 is in a fully turbulent region and we can observe a similar behavior as SD2 with an
even greater level of numerical dissipation (see Figure 7.15(e)). It is interesting to note
that "
vis
sub
in the sub-domains SD1 and SD2 is considerably greater than in SD3. This
behavior is due to the presence of the separated shear layer in the rst two sub-domains
leading to a higher viscous dissipation rate.
To date the methodology proposed by Schranner et al. [58] has been applied only to
a Taylor-Green vortex
ow which is a
ow that evolves in a box with periodic boundary
conditions. In order to both conrm that the higher viscous dissipation is eectively
coming from the shear layer and to investigate the eects of the presence of a no-slip
boundary condition, another set of three sub-domains have been analyzed with approx-
imatively the same sub-domain volumes. The new sub-domains that will be referred to
96
as SD1a, SD2a and SD3a have the same position along the cord as the previous set and
they are shifted away from the wall in the normal direction, the exact locations can be
seen in Figure 7.13. The number of cells per sub-domain is around 12; 000. The lower
count of cells is due to the cell size increasing in the wall-normal direction moving away
from the airfoil. It seems that indeed the high values of "
vis
sub
can be attributed to the
presence of a shear layer, for the value of the viscous dissipation rate in SD2a is similar
to the values in SD3 and SD3a, see Figure 7.15. The procedure appears pretty robust as
for the new sub-domains it predicts values of numerical viscosity that are slightly higher
for SD2a and SD3a with respect to the original sub-domains. This is consistent with the
fact that at the wall turbulence is damped.
A separate more careful analysis is needed for the sub-domain SD1a. As can been
seen from Figures 7.14(b) and 7.15(b) and Table 7.1, in this sub-domain the numerical
dissipation rate is negative, furthermore the calculated numerical viscosity in absolute
value is one order of magnitude higher than the physical viscosity. The method seems to
fail in this region of the
ow. The reason for this failure is that the region SD1a is laminar,
away from the walls, where the velocity gradients (@u
i
=@x
j
) are negligible. Since from
equation (7.12) the dissipation function " is a function of velocity gradients, in region
SD1a " (and the viscous dissipation) becomes negligible because of the absence of free-
stream turbulence. Additionally, we note that the computed numerical dissipation rate in
SD1a is two orders of magnitude smaller than in other regions considered. Therefore using
equation (7.15) involves a division of two small quantities which can lead to an unreliable
result, in this case an excessive value of
n
. What seems to be more troubling is that
the computed numerical dissipation is negative. This can be explained by looking at the
97
Sub-domain "
vis
10
5
"
n
10
5
" 10
4

n
10
4

tot
10
4

tot
=
SD1 1.42 0.05 0.62 2.29 0.08 2.37 1.04
SD2 1.76 2.82 0.79 2.24 3.55 5.79 2.59
SD3 0.58 1.60 0.28 2.05 5.66 7.71 3.75
SD1a 0.00 -0.02 0.00 2.14 -30.2 -28.1 -13.1
SD2a 0.67 1.28 0.30 2.21 4.19 6.40 2.89
SD3a 0.34 1.22 0.17 2.05 7.24 9.29 4.53
Table 7.1: Summary of averaged in time dissipation rates and viscosities, where
tot
= +
n
.
individual terms of the kinetic energy balance (see Figures 7.16 and 7.17). First we recall
that the terms in the balance are approximated with at best second-order reconstruction.
With this in mind we can observe that the numerical dissipation rate order of magnitude
is 10
7
while the kinetic energy
ux, the acoustic
ux, the work due to pressure and the
viscous dissipation rate are on the order of 10
5
, thus a 1% error in any of the latter
quantities could easily cause a sign switch in the numerical dissipation rate. Based on
these observations, in the regions of laminar
ows that do not contain strong velocity
gradients (for instance in our simulations whenever the sub-domain of interest is outside
the boundary layer) results produced by the method may be unreliable. However, this
is because the numerical dissipation in such regions is very small, within error bounds of
contributing terms, and thus of little signicance in the simulations. In other words, in
those regions the mesh resolution is adequate and the simulation results are not aected
by the truncation errors.
As already mentioned the UDNS simulation performs surprisingly well in the predic-
tion of the pressure and skin friction coecient (see Figure 6.11). This seems to suggest
that the numerical dissipation provided by the scheme is about the correct amount, i.e.
it is about the amount that an explicit sub-grid scale model would provide. To quantify
98
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(a) SD1
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(b) SD1a
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(c) SD2
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(d) SD2a
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(e) SD3
0 0.2 0.4 0.6 0.8 1
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
x 10
−5
t
ε

ε
n
ε
vis
(f) SD3a
Figure 7.14: Time-evolution of numerical and viscous dissipation rates. "
n
sub
(blue solid line)
and "
vis
sub
(red dashed line).
99
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(a) SD1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
−5
−4
−3
−2
−1
0
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(b) SD1a
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(c) SD2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(d) SD2a
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(e) SD3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν
(f) SD3a
Figure 7.15: Time-evolution of physical and numerical viscosities. Total viscosity +
n
(blue
solid line), physical viscosity (red dashed line) and average total viscosity +
n
(dot-dashed
black line).
100
0 0.2 0.4 0.6 0.8 1
−15
−10
−5
0
5
x 10
−6
t
ΔE
kin
/Δt

SD1
SD1a
(a) E
kin
=t
0 0.2 0.4 0.6 0.8 1
1
2
3
4
5
6
7
x 10
−5
t
F
kin

SD1
SD1a
(b) F
kin
0 0.2 0.4 0.6 0.8 1
−3
−2.5
−2
−1.5
−1
−0.5
0
x 10
−5
t
F
ac

SD1
SD1a
(c) F
ac
0 0.2 0.4 0.6 0.8 1
−2.6
−2.4
−2.2
−2
−1.8
−1.6
−1.4
−1.2
x 10
−7
t
F
vis

SD1
SD1a
(d) F
vis
Figure 7.16: Time-evolution of single terms of the kinetic energy balance, Eq. (7.14). Plotted
are the values for sub-domain SD1 (solid blue line) against SD1a (red dashed line). (a) change in
time of kinetic energy, (b) kinetic energy
ux, (c) acoustic
ux, (d) viscous
ux.
101
0 0.2 0.4 0.6 0.8 1
−4
−3
−2
−1
0
1
x 10
−5
t
W
p

SD1
SD1a
(e) W
p
0 0.2 0.4 0.6 0.8 1
−1.5
−1
−0.5
0
x 10
−5
t
ε
vis

SD1
SD1a
(f) "
vis
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
t
ε

SD1
SD1a
(g) "
0 0.2 0.4 0.6 0.8 1
−4
−2
0
2
4
6
8
x 10
−7
t
ε
n

SD1
SD1a
(h) "
n
Figure 7.17: Time-evolution of single terms of the kinetic energy balance, Eq. (7.14). Plotted
are the values for sub-domain SD1 (solid blue line) against SD1a (red dashed line). (e) pressure
work, (f) viscous dissipation rate, (g) dissipation function and (h) numerical dissipation rate.
102
this hypothesis a sub-grid scale viscosity has been estimated via the classical Smagorinsky
model

SGS
(m)
=C
2
s
4V
2=3
(m)

jS
ij
j

(m)
; (7.33)
whereC
s
is a constant andjS
ij
j = (2S
ij
S
ij
)
1=2
whereS
ij
is the strain tensor. In literature
the value ofC
s
can vary from 0:2 for fully isotropic turbulence to 0:1 for channel or shear

ows, the reduction of C
s
for the latter cases is justied by the fact the large scale non-
zero eld gradients add a contribution tojS
ij
j that has nothing to do with turbulence
[56].
SGS
(m)
is calculated for every cell in the sub-domain and volume averaged to obtain

SGS
. In Figure 7.18 the total viscosity provided by the code ( +
n
) and the total
viscosity that should be predicted by the SGS model ( +
SGS
) are compared. For all
sub-domains the SGS viscosity is provided for two C
s
values, 0:1 as a reference and a
value that allows the average SGS viscosity to match the averaged numerical viscosity,
0:06, 0:24, and 0:36 for sub-domains SD1, SD2 and SD3 respectively. Sub-domain SD1 is
in a laminar region thereforeC
s
should be zero. As expected we need a small value of the
Smagorinsky constant in order to match the numerical viscosity with the sub-grid scale
viscosity. For the sub-domains SD2 and SD3 the numerical dissipation provides a larger
value of viscosity than the Smagorinsky model would predict. For the sub-domain SD3,
in a region of fully developed turbulent boundary layer, a good match is achieved with
a value of C
s
= 0:36 which is 80% higher than what is typically found in literature for
isotropic turbulence and 260% higher than for wall bounded
ows. This suggests that the
numerical dissipation does not provide the correct amount of sub-grid-scale dissipation.
On the other hand, for this particular
ow conguration, it seems that the excess of
103
numerical dissipation, with respect to the sub-grid-scale dissipation, is not enough to be
appreciable in the chosen metrics (Pressure and skin friction coecients, separation and
reattachment points). These ndings also explain why there is a negligible dierence in
the results when an explicit LES model is added (see Figure 6.11) [14]. The numerical
dissipation is dominant over the SGS dissipation, therefore adding an explicit model does
not improve the results.
104
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(a) SD1, Cs = 0:10
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(b) SD1, Cs = 0:06
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(c) SD2, Cs = 0:10
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(d) SD2, Cs = 0:24
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(e) SD3, Cs = 0:10
0 0.2 0.4 0.6 0.8 1
0
0.5
1
1.5
x 10
−4
t
ν

ν + ν
n
ν + ν
n
ν + ν
SGS
ν
(f) SD3, Cs = 0:36
Figure 7.18: Time-evolution of physical, numerical and sub-grid-scale viscosities with dierent
values of C
s
. Total viscosity +
n
(blue line), physical viscosity (red dashed line), average
total viscosity +
n
(dot-dashed black line) and physical viscosity plus sub-grid scale viscosity

SGS
(blue line with markers) .
105
Chapter 8
Conclusions
8.1 Simulations
We have conducted a study to apply and evaluate LES techniques and an IB method
to simulate separated
ows. For LES we adopted the resolution target of the order of
O(1%) of DNS resolution, typically achievable for wall resolved LES of fully turbulent

ows.
We have used the numerical code INCA developed at TUM that uses an immersed
boundary method to simulate a
ow around a NACA-0012 airfoil for which a detailed
DNS database is available [40]. Two-dimensional and three-dimensional simulations at
low resolution around the reduction target have been performed, and at this resolution our
simulations reproduce qualitative features of the
ow such as the structure of the vorticity
eld and the presence of a separation. We also obtain a fairly good quantitative prediction
for the pressure coecient and the location of the primary separation point. However, we
have also determined that for the IB method a better quantitative agreement with the
benchmark data requires a DNS-like resolution; any computational savings with respect
106
to actual DNS can only be obtained by decreasing the resolution in the homogeneous
(spanwise) direction and in regions away from the wall. The reattachment point in
the three-dimensional case is quite sensitive to the resolution used. We speculate that
this behavior for the ALDM cases may be attributed to the noise coming from the IB
procedure at the wall rather than from not having the correct SGS dissipation. This claim
is supported by the results of two-dimensional simulations, where a dissipative code (5th
order WENO scheme) provided an excellent agreement with the DNS benchmark. A
DMD analysis revealed that the
ow dynamics predicted by the dissipative method is
fully captured with a single dynamic mode and its harmonics. Hence, vortex shedding
takes place at a constant streamwise coordinate, leading to pronounced secondary peaks
of C
f
and C
p
in agreement with the DNS. However, with a less dissipative method, the
vortex shedding is essentially dominated by two dierent low-frequency dynamic modes,
which leads to a slight upstream-downstream oscillation of the vortex shedding location
and a smeared C
f
and C
p
distribution in the temporal average. A separate conclusion
regarding LES for such problems in the framework of IB methods is that because of
diculties in resolving geometry of a curved wall with a Cartesian grid, the IB methods
are more appropriate for a high resolution/delity DNS and LES rather than for very
coarse LES which can better take advantage of high-aspect-ratio body tted grids to
resolve the wall region.
The incorrect reattachment point prediction in the three-dimensional case with the
5th order WENO scheme can be explained by the insucient dissipation given by the
lack of spanwise resolution. This is demonstrated by the good agreement between our
DNS and the reference DNS results of Jones et al.
107
The results obtained using Star-CCM+ help to give a broader answer to the feasibility
of very low resolution LES for separated
ows. The good performance of Star-CCM+
shows that the target of 1% of DNS resolution can be reached, conrming what has already
been observed for a dierent geometrical conguration [10]. However it is important to
note that the good performance achieved with the commercial code cannot guarantee
good predictions for new congurations. This is because the numerical dissipation of
the scheme dominates the SGS dissipation provided by explicit models. While in this
particular case, at this specic resolution, the numerical dissipation mimics the physical
SGS dissipation correctly, there is no mechanism in the numerical method that would
guarantee the same for other problems.
Since for this particular case the agreement with the reference DNS results was very
good a sweep of angles of attack for the two-dimensional case has been done. The sweep
of angles of attack allow for the prediction of the jump in the lift and drag coecients;
the polar plot has some similarity to the experimental results of Spedding and McArthur
[61]. In our case the LSB appears almost immediately and the jump in lift and drag
coecients is located between 1:5

and 3

.
8.2 Numerical dissipation
For the rst time the method proposed by Schranner et al. [58] has been applied
to a compressible low-order Navier-Stokes solver with a realistic
ow conguration. The
method seems quite robust and self-consistent and it seems to be easily applicable to
wall-bounded
ows. A previously unknown limitation to the method has been found,
108
specically the method will fail to compute accurately the numerical viscosity for lam-
inar
ows that do not contain strong velocity gradients, for instance in our simulations
whenever the sub-domain of interest is chosen outside the boundary layer. A posteriori
this seems a quite obvious nding given that the numerical viscosity is computed from
equation (7.15) and the dissipation function tends to zero in absence of velocity gradients.
The above described limitation of the method was not observed before due to the fact
that in TGV
ows, even though at initial times the
ow is fully laminar, strong gradients
are always present in the
ow due to initial conditions. This limitation does not dimin-
ish the utility of this tool since the method is designed for the estimation of numerical
dissipation in turbulent
ows.
The application of the method to the NACA-0012 simulations performed with Star-
CCM+ enabled the quantication of the initially qualitative ndings reported in Chap-
ter 6. Two important observations have been quantied. First we have quantied the
amount of numerical dissipation in dierent regions of the
ow eld in a coarse UDNS
with a low-order scheme. Second we have been able to quantify the previously qualitative
observation that the numerical dissipation dominates over the sub-grid-scale dissipation
when an explicit model is added. The analyses conrmed that the addition of an ex-
plicit LES model to a coarse UDNS with a low-order scheme is at best useless and can
potentially even deteriorate the quality of a UDNS as previously reported [43, 9].
109
8.3 Outlook
This study points out the limitations of both IB codes and low-order methods codes.
On the other hand it also shows that very coarse LES with realistic geometries is feasible.
The missing link to enable predictive capabilities is an unstructured high-order code which
would allow for the computation of complex geometries and little numerical dissipation.
The method presented in Chapter 7 can be used to evaluate a posteriori the quality of a
code-grid resolution combination. Furthermore, it can be used to control the numerical
dissipation even in low-order schemes by increasing the resolution where needed. This
can be implemented in an iterative procedure and would allow for a predictive tool, albeit
computationally expensive.
110
References
[1] M. Alam and N.D. Sandham. Direct numerical simulation of `short' laminar separa-
tion bubbles with turbulent reattachment. J. Fluid Mech., 410:1{28, 2000.
[2] J. H. Almutairi, L. E. Jones, and N. D Sandham. Intermittent bursting of a laminar
separation bubble on an airfoil. AIAA J., 48(2):414{426, 2010.
[3] A.E. Alving and H.H. Fernholz. Turbulence measurements around a mild separation
bubble and downstream of reattachment. J. Fluid Mech., 322(1):297{328, 1996.
[4] J. D. Anderson. Modern compressible
ow: with historical perspective. McGraw-Hill,
3rd edition, 1990.
[5] J. D. Anderson. Computational
uid dynamics. McGraw-Hill, 1st edition, 1995.
[6] J. D. Anderson. Ludwig Prandtls boundary layer. Phys. Today, 58(12):42{48, 2005.
[7] J.P. Boris, F.F. Grinstein, E.S. Oran, and R.L. Kolbe. New insights into large eddy
simulation. Fluid Dyn. Res., 10:199{228, 1992.
[8] F. Cadieux, J. A. Domaradzki, T. Sayadi, S. Bose, and F. Duchaine. DNS and LES of
separated
ows at moderate Reynolds numbers. In Preceedings of the 2012 Summer
Program, pages 77{86. Center for Turbulence Research, Stanford University, 2012.
[9] F. Cadieux, J. A. Domaradzki, T. Sayadi, and T. Bose. DNS and LES of laminar sep-
aration bubbles at moderate Reynolds numbers. ASME J. Fluids Eng., 136:061102,
2014.
[10] F. Cadiuex, G. Castiglioni, J.A. Domaradzki, T. Sayadi, S. Bose, M. Grilli, and
S. Hickel. LES of separated
ows at moderate Reynolds numbers appropriate for
turbine blades and unmanned aero vehicles. In Proceedings of 8th International
Symposium on Turbulence and Shear Flow Phenomena (TSFP8), 2013.
[11] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods,
fundamentals in single domains. Springer, 2006.
[12] G. Castiglioni and J.A. Domaradzki. A numerical dissipation rate and viscosity in

ow simulations with realistic geometry using low-order compressible Navier-Stokes
solvers. Comp. & Fluids, submitted.
111
[13] G. Castiglioni, J.A. Domaradzki, M. Grilli, and S. Hickel. Numerical modeling
of separated
ows at moderate Reynolds numbers appropriate for turbine blades
and unmanned aero vehicles. In Proceedings of the rst TRR40 Summer Program,
Technische Universit at M unchen, pages 67{76, 2011.
[14] G. Castiglioni, J.A. Domaradzki, V. Pasquariello, S. Hickel, and M. Grilli. Numerical
simulations of separated
ows at moderate Reynolds numbers appropriate for turbine
blades and unmanned aero vehicles. Int. J. Heat and Fluid Flow, 49:91{99, 2014.
[15] G. Castiglioni, Domaradzki J.A., V. Pasquariello, S. Hickel, and M. Grilli. Numerical
modeling of 3-D separated
ows at Reynolds numbers appropriate for turbine blades
and unmanned aero vehicles. In Proceedings of the second TRR40 Summer Program,
Technische Universit at M unchen, pages 177{190, 2013.
[16] P. Catalano and R. Tognaccini. Turbulence modeling for low-Reynolds-number
ows.
AIAA J., 48(8):1673{1685, 2010.
[17] CD-Adapco. Star-CCM+ manual, 2013. Version 8.02.
[18] H. Choi and P. Moin. Eects of the computational time step on numerical solutions
of turbulent
ow. J. Comp. Phys., 113(1):1 { 4, 1994.
[19] H. Choi and P. Moin. Grid-point requirements for large eddy simulation: Chapmans
estimates revisited. Phys. Fluids, 24(1):011702, 2012.
[20] Joshua N.N. Counsil and K. Goni Boulama. Low-Reynolds-number aerodynamic per-
formances of the NACA 0012 and Selig{Donovan 7003 airfoils. J. Aircr., 50(1):204{
217, 2013.
[21] L. Davidson, D. Cokljat, J. Fr ohlich, M. A. Leschziner, C. Mellen, and W. Rodi.
LESFOIL: Large Eddy Simulation of Flow Around a High Lift Airfoil: Results of the
Project LESFOIL Supported by the European Union 1998-2001, volume 83. Springer,
2003.
[22] J. A. Domaradzki, Z. Xiao, and P. K. Smolarkiewicz. Eective eddy viscosities in
implicit large eddy simulations of turbulent
ows. Phys. Fluids, 15(12):3890{3893,
2003.
[23] J.A. Domaradzki and S. Radhakrishnan. Eective eddy viscosities in implicit large
eddy simulations of decaying high Reynolds number turbulence with and without
rotation. Fluid Dyn. Res., 36:385{406, 2005.
[24] Y. Dubief and F. Delcayre. On coherent-vortex identication in turbulence. J.
Turbul., 1(1):011{011, 2000.
[25] S. Eisenbach and R. Friedrich. Large-eddy simulation of
ow separation on an airfoil
at a high angle of attack andRe = 10
5
using Cartesian grids. Theor. Comput. Fluid
Dyn., 22:213{225, 2008.
112
[26] C. Fletcher. Computational techniques for
uid dynamics 1, volume 1. Springer,
1991.
[27] M. Galbraith and M. Visbal. Implicit large eddy simulation of low reynolds number

ow past the SD7003 airfoil. In 46th AIAA Aerospace Sciences Meeting and Exhibit,
pages 2008{225, 2008.
[28] E. Garnier, N. Adams, and P. Sagaut. Large eddy simulation for compressible
ows.
Springer, 2009.
[29] E. Garnier, M. Mossi, P. Sagaut, P. Comte, and M. Deville. On the use of shock-
capturing schemes for large-eddy simulation. J. Comp. Phys., 153(2):273{311, 1999.
[30] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgrid-scale eddy
viscosity model. Phys. Fluids A, 3:1760{1765, 1991.
[31] S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes. Math.
Comp., 67(221):73{85, 1998.
[32] M. Grilli, S. Hickel, X.Y. Hu, and N.A. Adams. Conservative immersed bound-
ary method for compressible
ows. Academy Colloquium on Immersed boundary
methods: current status and future reaserch directions, 2009. Amsterdam, The
Netherlands.
[33] F. Grinstein, L. Margolin, and W. Rider. Implicit Large Eddy Simulation: Comput-
ing Turbulent Flow. Cambridge University Press, 2007.
[34] C. H aggmark. Investigations of disturbances developing in a laminar separation
bubble
ow. PhD thesis, Kungl Tekniska H ogskolan, 2000.
[35] R. Hain, C.J. Kaehler, and R. Radespiel. Dynamics of laminar separation bubbles
at low-Reynolds-number aerofoils. J. Fluid Mech., 630:129{153, 2009.
[36] S. Hickel, N.A. Adams, and J.A. Domaradzki. An adaptive local deconvolution
method for implicit LES. J. Comp. Phys., 213:413{436, 2006.
[37] S. Hickel and J. Larsson. An adaptive local deconvolution model for compressible
turbulence. In Proceedings of the 2008 Summer Program, pages 85{96. Center for
Turbulence Research, Stanford University, 2008.
[38] H. Hu, Z. Yang, and H. Igarashi. Aerodynamic hysteresis of a low-Reynolds-number
airfoil. J. Aircr., 44:2083{2086, 2007.
[39] L.E. Jones. Numerical Studies of the Flow around an Airfoil at Low Reynolds Num-
ber. PhD thesis, University of Southampton, 2008.
[40] L.E. Jones, R.D. Sandberg, and N.D. Sandham. Direct numerical simulations of
forced and unforced separation bubbles on an airfoil of incidence. J. Fluid Mech.,
602:175{207, 2008.
113
[41] L.E. Jones, R.D. Sandberg, and N.D. Sandham. Stability and receptivity charac-
teristics of a laminar separation bubble on an airfoil. J. Fluid Mech., 648:257{296,
2010.
[42] R. Kojima, T. Nonomura, A. Oyama, and K. Fujii. Large-eddy simulation of low-
Reynolds-number
ow over thick and thin NACA airfoils. J. Aircr., 50(1):187{196,
2013.
[43] A.G. Kravchenko and P. Moin. On the eect of numerical errors in large eddy
simulations of turbulent
ows. J. Comp. Phys., 131(2):310{322, 1997.
[44] A. Leonard. Energy cascade in large-eddy simulations of turbulent
uid
ows. Adv.
in geophys., 18:237{248, 1975.
[45] D.K. Lilly. A proposed modication of the Germano subgrid-scale closure method.
Phys. Fluids A, 4:633{635, 1992.
[46] J. C. M. Lin and L. L. Pauley. Low-Reynolds-number separation on an airfoil. AIAA
J., 34:1570{1577, 1996.
[47] G. Lodato, P. Domingo, and L. Vervisch. Three-dimensional boundary conditions
for direct and large-eddy simulation of compressible viscous
ows. J. Comp. Phys.,
227(10):5105{5143, 2008.
[48] M. P. Martin, U. Piomelli, and G. V. Candler. Subgrid-scale models for compressible
large-eddy simulations. Theoretical and Computational Fluid Dynamics, 13(5):361{
376, 2000.
[49] M. Meyer, A. Devesa, S. Hickel, X.Y. Hu, and N.A. Adams. A conservative immersed
interface method for Large-Eddy Simulation of incompressible
ows. J. Comp. Phys.,
229(18):6300{6317, 2010.
[50] F. Nicoud and F. Ducros. Subgrid-scale stress modeling based on the square of the
velocity gradient tensor. Flow, Turbul. Combust., 62(3):183{200, 1999.
[51] R. L. Panton. Incompressible
ow. John Wiley & Sons, 3rd edition, 2006.
[52] T. J. Poinsot and S.K. Lele. Boundary conditions for direct simulations of compress-
ible viscous
ows. J. Comp. Phys., 101(1):104{129, 1992.
[53] S. B. Pope. Turbulent
ows. Cambridge University Press, 2000.
[54] M. Popov, K. Mikhaylyukov, A. Ryabov, F. Mendonca, and C. Tzanos. Large eddy
simulation of turbulent channel
ow with heat transfer at low Reynolds and Mach
numbers. ICHMT Digital Libr. Online, 10, 2006.
[55] D. H. Rudy and J. C. Strikwerda. A nonre
ecting out
ow boundary condition for
subsonic Navier-Stokes calculations. J. Comp. Phys., 36(1):55{70, 1980.
[56] P. Sagaut. Large eddy simulation for incompressible
ows: an introduction. Springer,
3rd edition, 2006.
114
[57] P. J. Schmid. Dynamic mode decomposition of numerical and experimental data. J.
Fluid Mech., 656:5{28, 2010.
[58] F.S. Schranner, J.A. Domaradzki, S. Hickel, and N.A. Adams. Assessing the numer-
ical dissipation rate and viscosity in numerical simulations of
uid
ows. Comp. &
Fluids, 114:84{97, 2015.
[59] C.-W. Shu, W.-S. Don, D. Gottlieb, O. Schilling, and L. Jameson. Numerical con-
vergence study of nearly incompressible, inviscid TaylorGreen vortex
ow. J. Sci.
Comp., 24(1), 2005.
[60] P.R. Spalart and M.K. Strelets. Mechanisms of transition and heat transfer in a
separation bubble. J. Fluid Mech., 403:329{349, 2000.
[61] G.R. Spedding and J. McArthur. Span eciencies of wings at low Reynolds numbers.
J. Aircr., 47:120{128, 2010.
[62] D. Taieb, G. Ribert, and A. Hadjadj. Numerical simulations of shock focusing over
concave surfaces. AIAA J., 48:1739{1747, 2010.
[63] T. Tantikul and J. A. Domaradzki. Large eddy simulations using truncated Navier-
Stokes equations with the automatic ltering criterion. J. Turbul., 11(21):1{24, 2010.
[64] G.I. Taylor and A.E. Green. Mechanism of the production of small eddies from large
ones. Proc. R. Soc. Lond. A, 158(895):499{521, 1937.
[65] K. W. Thompson. Time dependent boundary conditions for hyperbolic systems. J.
Comp. Phys., 68(1):1{24, 1987.
[66] B. Vreman. Direct and large-eddy simulation of the compressible turbulent mixing
layer. PhD thesis, University of Twente, 1995.
[67] B. Vreman, B. Geurts, and H. Kuerten. A priori tests of large eddy simulation of
the compressible plane mixing layer. J. Eng. Math., 29(4):299{327, 1995.
[68] D. C. Wilcox. Turbulence modeling for CFD. DCW industries La Ca~ nada, 3rd
edition, 2006.
[69] P.G. Wilson and L.L. Pauley. Two- and three-dimensional large-eddy simulations of
a transitional separation bubble. Phys. Fluids, 10:2932{2940, 1998.
[70] Z. Yang and P.R. Voke. Large-eddy simulation of boundary-layer separation and
transition at a change of surface curvature. J. Fluid Mech., 439:305{333, 2001.
115

Abstract (if available)

Abstract Flows over airfoils and blades in rotating machinery, for unmanned and micro-aerial vehicles, wind turbines, and propellers consist of a laminar boundary layer near the leading edge that is often followed by a laminar separation bubble and transition to turbulence further downstream. Typical Reynolds averaged Navier-Stokes turbulence models are inadequate for such flows. Direct numerical simulation is the most reliable, but is also the most computationally expensive alternative. This work assesses the capability of immersed boundary methods and large eddy simulations to reduce the computational requirements for such flows and still provide high quality results. Two-dimensional and three-dimensional simulations of a laminar separation bubble on a NACA-0012 airfoil at Rec = 5×10⁴ and at 5° of incidence have been performed with an immersed boundary code and a commercial code using body fitted grids. Several sub-grid scale models have been implemented in both codes and their performance evaluated. For the two-dimensional simulations with the immersed boundary method the results show good agreement with the direct numerical simulation benchmark data for the pressure coefficient Cp and the friction coefficient Cf, but only when using dissipative numerical schemes. There is evidence that this behavior can be attributed to the ability of dissipative schemes to damp numerical noise coming from the immersed boundary. For the three-dimensional simulations the results show a good prediction of the separation point, but an inaccurate prediction of the reattachment point unless full direct numerical simulation resolution is used. The commercial code shows good agreement with the direct numerical simulation benchmark data in both two and three-dimensional simulations, but the presence of significant, unquantified numerical dissipation prevents a conclusive assessment of the actual prediction capabilities of very coarse large eddy simulations with low order schemes in general cases. Additionally, a two-dimensional sweep of angles of attack from 0° to 5° is performed showing a qualitative prediction of the jump in lift and drag coefficients due to the appearance of the laminar separation bubble. The numerical dissipation inhibits the predictive capabilities of large eddy simulations whenever it is of the same order of magnitude or larger than the sub-grid scale dissipation. The need to estimate the numerical dissipation is most pressing for low-order methods employed by commercial computational fluid dynamics codes. Following the recent work of Schranner et al., the equations and procedure for estimating the numerical dissipation rate and the numerical viscosity in a commercial code are presented. The method allows for the computation of the numerical dissipation rate and numerical viscosity in the physical space for arbitrary sub-domains in a self-consistent way, using only information provided by the code in question. The method is first tested for a three-dimensional Taylor-Green vortex flow in a simple cubic domain and compared with benchmark results obtained using an accurate, incompressible spectral solver. Afterwards the same procedure is applied for the first time to a realistic flow configuration, specifically to the above discussed laminar separation bubble flow over a NACA 0012 airfoil. The method appears to be quite robust and its application reveals that for the code and the flow in question the numerical dissipation can be significantly larger than the viscous dissipation or the dissipation of the classical Smagorinsky sub-grid scale model, confirming the previously qualitative finding.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Large eddy simulations of laminar separation bubble flows

PDF

On the simulation of stratified turbulent flows

PDF

Large eddy simulations of turbulent flows without use of the eddy viscosity concept

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Boundary layer and separation control on wings at low Reynolds numbers

PDF

Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations

PDF

Numerical simulations of linearly stratified flow around a sphere

PDF

Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders

PDF

RANS simulations for flow-control study synthetic jet cavity on NACA0012 and NACA65(1)412 airfoils.

PDF

Direct numerical simulation of mixing and combustion under canonical shock turbulence interaction

PDF

Biomimetics and bio-inspiration for moderate Reynolds number airfoils and aircraft

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Three-dimensional kinetic simulations of whistler turbulence in solar wind on parallel supercomputers

PDF

A subgrid-scale model for large-eddy simulation based on the physics of interscale energy transfer in turbulence

PDF

Dynamic modeling and simulation of flapping-wing micro air vehicles

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

PDF

Studies of combustion characteristics of heavy hydrocarbons in simple and complex flows

PDF

Numerical study of focusing effects generated by the coalescence of multiple shock waves

PDF

Grid-based Vlasov method for kinetic plasma simulations

PDF

Modeling and simulation of complex recovery processes

Asset Metadata

Creator Castiglioni, Giacomo (author)

Core Title Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 06/22/2015

Defense Date 05/20/2015

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag airfoil,CFD,laminar separation bubble,large eddy simulations,OAI-PMH Harvest,turbulence

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Domaradzki, Julian A. (committee chair), Eliasson, Veronica (committee member), Nakano, Aiichiro (committee member), Spedding, Geoffrey R. (committee member)

Creator Email castigli@usc.edu,giacomocastiglioni@hotmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-575375

Unique identifier UC11301640

Identifier etd-Castiglion-3500.pdf (filename),usctheses-c3-575375 (legacy record id)

Legacy Identifier etd-Castiglion-3500.pdf

Dmrecord 575375

Document Type Dissertation

Format application/pdf (imt)

Rights Castiglioni, Giacomo

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA