Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders
(USC Thesis Other)
Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
i
Numerical Study of Flow Characteristics of Controlled
Vortex Induced Vibrations in Cylinders
By
Amirhossein Eftekharian
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CIVIL ENGINEERING)
August 2018
© AMIRHOSSEIN EFTEKHARIAN 2018
ii
Dedicated to my beloved family especially my late father, Manouchehr
iii
Acknowledgement
First and foremost, I am grateful to my PhD advisor, Professor Jiin-Jen Lee, for his
invaluable advice, understanding, and assistance throughout my PhD program as well as
providing the opportunity for me to pursue my studies at the University of Southern California. It
has been a special honor for me working as his very last PhD student at USC as well. I would
also like to express my gratitude to my dissertation committee, Professors Sami Masri, Carter
Wellford, Julian Domaradzki, and Patrick Lynnett for accepting to review my work as well as
their time, profound comments, guidance, and mentorship.
I would like to thank the USC High Performance Computing and Communications staff for
showing support and enthusiasm regarding my simulations on the USC cluster. I would also like
to thank the CD-adapco engineers and customer service for their academic support of the CFD
software.
I would like to express my appreciation to my friends Nima Jabbari, Iman Yadegaran, Cyrus
Ashayeri, Amir Mani, Azadeh Rouhani, Mohammad Mehryar, Sepehr Dara, Babak Zaeriyan,
Hossein Yousefpour, Reza Jafarkhani, Ali Bolourchi, and Farrokh Jazizadeh and my former
group members especially Mehrdad Bozorgnia for their valuable support, insight, and friendship.
I would also like to thank the USC Civil and Environmental Engineering and Aerospace and
Mechanical Engineering Departments faculty and staff more specifically Professors Najmedin
Meshkati, Gregg Brandow, Henry Koffman, Iraj Nasseri, and David Wilcox. I feel also thankful
to work with the USC Foundation for Cross-Control Connection and Hydraulic Research staff
iv
and engineers in the course of my PhD. I would like to take the opportunity to thank Dr Frank
Abdi and the AlphaSTAR Corporation team as well.
Last but not least, my special thanks goes to my lovely family in Iran, my mother Pari, my
late father, Manouchehr who passed away few months before the completion of my PhD, my
sisters Leila and Sara, and my brother Hossein, for their never-ending love, support, and
encouragement from day one of my education. I would have never made it through this journey
without them.
v
Abstract
Vortex Induced Vibration (VIV) is the vibration of a solid object due to its interaction with
the surrounding fluid. VIV in cylinders has been getting a lot of attention because it is a
fundamental fluid mechanics problem which deals with complicated flow structures such as
separations-reattachments, laminar-to-turbulent transitions, and streamlines curvatures.
Additionally, it is a real engineering challenge with applications in offshore risers, bridges, and
power plant towers. Studying VIV by imposing a prescribed oscillation to the cylinder, which is
called controlled VIV, can lead to a better insight of flow characteristics because of the reduction
in the complexity of the problem. Moreover, experimental studies have shown that controlled
VIV, for transient Reynold numbers, is able to predict the free VIV if a fine resolution is used in
selection of the experimental parameters.
The current study has been focused on numerical simulation of vortex shedding behind a
stationary and moving cylinder for Re = 4,000. Numerous simulations with different mesh
configurations, time step set-ups, and turbulence models as well as the mesh motion techniques
have been used to study different parameters’ effects on the statistical values of fluid excitation.
2D simulations with RANS turbulence models show that, for neither of cases, this type of
simulation is capable of reproducing the experimental data accurately. 2D simulations with
RANS predict specifically much higher lift forces and their resulted time history of forces does
not demonstrate any amplitude modulations for the stationary cylinder case. The vortex shedding
pattern behind the moving cylinder does not agree with the experiment either and 2D simulations
seem to be unable to predict the switches in the shedding pattern observed in the experiment.
vi
Solving the stationary problem in a 3D domain with RANS seems to mitigate some of the 2D
simulation drawbacks. Furthermore, 3D simulations reveal the significance of the damping
introduced by some of the RANS models near the cylinder wall. Damping in the RANS closure
seems to be the main obstacle to the development of 3D flow structures such as breakage of
planar shear layers in the spanwise direction and existence of inline vorticity. Anisotropy and
nonlinearity in turbulence models seem to neutralize the damping effects in some cases, but it
seems that models with no damping functions are capable of triggering three dimensionality. In
addition, damping functions show their trace in the result of 2D simulations as well: Models with
damping functions lead to lower lift values in general, and damping seems to eliminate or
mitigate amplitude modulations and moderate the vorticity field in the moving case.
Migrating to scale resolving turbulence models, i.e. DES and LES, immediately shows its
effects on the statistical values of flow for the stationary cylinder simulation. CFD simulations
show that using a detached eddy simulation (DES) model with the same 3D mesh count drops
the statistical values of lift and drag dramatically. A detailed adjustment of radial,
circumferential, and spanwise discretization in an O-type domain can lead to an acceptable set of
results using DES based on the Spalart-Allmaras model with a relatively low mesh count and
high time step. This is significant due to the fact that this CFD study requires less computational
resources compared to the majority of available data in the literature. Moreover, the mesh
configuration derived here can later be used as a baseline for the cases where the cylinder is
under a prescribed motion.
3D simulations using the scale resolving turbulence models for Re = 4,000 illustrates that
CFD seems to be capable of predicting the vortex shedding mode and the transferred energy
between the cylinder and fluid for the studied range of oscillation amplitude and frequency with
vii
an acceptable accuracy. However, CFD results do now demonstrate any switches between vortex
shedding patterns in time observed in the experimental data for the studied range of parameters.
On the other hand, the fluid excitation predicted by CFD has a distinguished match with the test
data implying that CFD seems to be able to reproduce the experimental data however the newly
discovered shedding mode is not seen in CFD results in this range of parameters.
Implementing more accurate and at the same time more computationally expensive
simulations in terms of mesh count, temporal and spatial resolutions, turbulence model physical
fidelity, and mesh motion technique seems to add more details to the predictions of flow.
However, the overall behavior of predictions does not change and the added value of the more
demanding simulations is insignificant for CFD simulation of flow past an oscillating cylinder
for the studied range of parameters.
viii
Table of Contents
Acknowledgement ........................................................................................................................ iii
Abstract .......................................................................................................................................... v
Chapter 1- Background and Literature Review ........................................................................ 1
1-1- Introduction ........................................................................................................................... 1
1-2- Literature Review ................................................................................................................. 2
1-2-1- Experimental Studies ...................................................................................................... 4
1-2-2- Numerical Studies ......................................................................................................... 16
Chapter 2- Numerical Formulations ......................................................................................... 21
2-1- Introduction ......................................................................................................................... 21
2-2- Governing Equations .......................................................................................................... 22
2-3- Discretization ...................................................................................................................... 23
2-3-1- Momentum Equation in Discrete Form......................................................................... 24
2-3-2- Continuity Equation in Discrete Form .......................................................................... 25
2-4- SIMPLE Solver Algorithm ................................................................................................. 28
2-5- Numerical Method for Solving Algebraic Equations ......................................................... 29
2-6- Turbulence Modeling .......................................................................................................... 31
2-6-1- Direct Numerical Simulation (DNS) ............................................................................. 31
2-6-2- Reynolds Averaged Navier-Stokes (RANS) ................................................................. 31
2-6-3- Large Eddy Simulation (LES) ....................................................................................... 35
2-6-4- Detached Eddy Simulation (DES) ................................................................................ 40
Chapter 3- Analysis of Flow past a Cylinder by 2D Simulations ........................................... 42
3-1- Introduction ......................................................................................................................... 42
3-2- Two Dimensional Simulation of Flow past a Stationary Cylinder ..................................... 44
3-2-1- Flow Domain and Boundary Conditions....................................................................... 44
3-2-2- Two Dimensional Domain Meshing ............................................................................. 46
3-2-3- Results for 2D Simulation of a Stationary Cylinder ..................................................... 57
3-2-4- Discussion on Results for the Stationary Cylinder ....................................................... 62
3-3- Two Dimensional Simulation of Flow past a Moving Cylinder ......................................... 66
ix
3-3-1- Parameters of the Model ............................................................................................... 67
3-3-2- Apply the Motion to the Cylinder ................................................................................. 68
3-3-3- Results for 2D Simulation of the Moving Cylinder ...................................................... 70
3-4- Concluding Remarks ........................................................................................................... 74
3-4-1- Two Dimensional Simulation of Flow past a Stationary Cylinder ............................... 74
3-4-2- Two Dimensional Simulation of Flow past a Moving Cylinder ................................... 75
Chapter 4- Analysis of Flow past a Stationary Cylinder by 3D Simulations and RANS
Models .......................................................................................................................................... 77
4-1- Introduction ......................................................................................................................... 77
4-2- Model Parameters ............................................................................................................... 79
4-2-1- Planar Mesh Configuration ........................................................................................... 79
4-2-2- Mesh Discretization and Domain Size in the z Direction ............................................. 81
4-2-3- Turbulence Model ......................................................................................................... 83
4-2-4- Boundary Conditions on Top and Bottom .................................................................... 83
4-2-5- Controlling Parameters and Running the Simulations .................................................. 83
4-3- Results of 3D Simulation for Flow past a Stationary Cylinder .......................................... 84
4-4- Discussion ........................................................................................................................... 91
4-4-1- Isotropy ......................................................................................................................... 91
4-4-2- Wall Treatment.............................................................................................................. 94
4-4-3- Low-Re Models ............................................................................................................. 95
4-5- Justifications ....................................................................................................................... 96
4-5-1- Why Do the Standard and Realizable k – ε Models with all y
+
Wall Treatment Show no
Three Dimensionality? ............................................................................................................. 97
4-5-2- Why Do the AKN Low-Re Low y
+
, v
2
f – Low y
+
, and EB – Low y
+
k – ε Models
Produce Three Dimensionality? ............................................................................................... 97
4-5-3- Why Does the Reynolds Stress Model with Two Layer Treatment Have no Three
Dimensionality? ....................................................................................................................... 97
4-5-4- Why Does the Nonlinear Constitutive Relation Help to Develop Three Dimensionality?
.................................................................................................................................................. 98
4-5-5- Why Do all Versions of the k – ω and Spalart-Allmaras Models Show Some Level of
Three Dimensionality? ............................................................................................................. 98
4-5-6- Why Is Three Dimensionality Observed in Some of the Models with High y
+
Treatments? .............................................................................................................................. 99
x
4-6- Back to 2D Simulation of Flows past a Cylinder ............................................................. 100
4-6-1- Effects of Turbulence Models on the Stationary Cylinder Results ............................. 100
4-6-2- Effects of Turbulence Models on the Moving Cylinder Results................................. 101
4-7- Concluding Remarks ......................................................................................................... 101
Chapter 5- Analysis of Flow past a Stationary Cylinder by 3D Simulations and DES
Models ........................................................................................................................................ 104
5-1- Introduction ....................................................................................................................... 104
5-2- Three Dimensional Simulation of Flow past a Stationary Cylinder ................................. 105
5-2-1- Flow Domain and Boundary Conditions..................................................................... 105
5-3- Results of 3D DES Simulation ......................................................................................... 117
5-4- Concluding Remarks ......................................................................................................... 122
Chapter 6- Analysis of Flow past an Oscillating Cylinder with Three Dimensional
Simulations ................................................................................................................................ 123
6-1- Introduction ....................................................................................................................... 123
6-2- Model Parameters ............................................................................................................. 124
6-2-1- Applying the Cylinder Oscillation in Modeling ......................................................... 126
6-2-2- Turbulence Modeling .................................................................................................. 129
6-2-3- Simulation Initiation and Inner Iteration ..................................................................... 129
6-2-4- Added Mass ................................................................................................................. 130
6-2-5- High Performance Computing .................................................................................... 131
6-3- Results of Flow past an Oscillating Cylinder ................................................................... 133
6-3-1- Calculating the Phase Lag ........................................................................................... 136
6-3-2- Results for Total and Vortex Lifts and Their Phase Lags ........................................... 137
6-4- Discussion on Results ....................................................................................................... 149
6-5- Comparison between CFD and Test data ......................................................................... 150
6-5-1- Total and Vortex Lift and Their Corresponding Phase Lags in Time......................... 150
6-5-2- Fluid Excitation for Selected Data Points ................................................................... 155
6-6- Comparing CFD Results for Different Data Points .......................................................... 158
6-6-1- Changes in Total Lift for a Cut of A* = 0.75 .............................................................. 159
6-6-2- Changes in Total Lift for a Cut of λ* = 6.0 ................................................................. 160
6-6-3- Sensitivity of Total Lift to Turbulence Modeling ....................................................... 160
xi
6-6-4- Sensitivity of Total Lift to Spatial Resolution ............................................................ 163
6-6-5- Sensitivity of Total Lift to Temporal Resolution ........................................................ 167
6-6-6- Sensitivity of Total Lift to Mesh Motion Technique .................................................. 167
6-6-7- Effects of Convection Terms Resolution Schemes on Total Lift ............................... 169
6-6-8- The Last Comparison .................................................................................................. 170
6-6-9- Discussion on Numerical Dissipation ......................................................................... 172
6-7- Concluding Remarks ......................................................................................................... 173
Conclusions ................................................................................................................................ 175
References .................................................................................................................................. 178
xii
List of Figures
Figure 1-1- Tacoma Narrow Bridge collapse in 1940 (picture from www.exit133.com) .............. 2
Figure 1-2 – Regimes of fluid across smooth circular cylinders (Lienhard 1966) – Figure by MIT
OCW ............................................................................................................................................... 5
Figure 1-3 – Vortex pattern for different vortex regimes (Williamson, 1996) ............................... 6
Figure 1-4 – Components of fluid excitation on a cylinder ............................................................ 7
Figure 1-5 – Response amplitude and frequency for low mass damping. (a) Amplitude response
for low (■), and high (◊) mass damping, (b) Oscillation frequency through the excitation regime
in low mass damping (Khalak & Williamson, 1997) ..................................................................... 8
Figure 1-6 – Williamson and Roshko’s regimes of vortex shedding............................................ 10
Figure 1-7 – Map of vortex shedding regimes in Morse and Williamson’s studies ..................... 11
Figure 1-8 – Vortex shedding modes observed in Morse and Williamson’s experiments ........... 12
Figure 1-9 – Measured (○) and predicted (●) amplitude response for free VIV at (a) low mass-
damping and (b) high mass-damping (Re = 4,000) ...................................................................... 15
Figure 1-10 – Map of the zone where free vibration can occur in the amplitude – velocity plane
for m
*
=10 ...................................................................................................................................... 16
Figure 2-1 – Control Volume for cell (i,j) at time level n, with flux densities H
n
on sides, side
vectors S , and center data u (from Causon et al., 2011) .............................................................. 24
Figure 2-2 – Scales of turbulence that are resolved or modeled in different methods (Sodja 2007)
....................................................................................................................................................... 35
Figure 2-3 – Turbulence energy spectrum for LES model. Wavenumbers greater than k2 need to
be modeled (from Apte) ................................................................................................................ 36
Figure 3-1 – RMS values of lift coefficient in cylinders for different Reynolds numbers (Norberg
2003) ............................................................................................................................................. 43
Figure 3-2 – Strouhal number of cylinders for different Reynolds numbers (Norberg 2003) ...... 44
Figure 3-3 – Domain size for 2D simulations ............................................................................... 45
Figure 3-4 – Boundary conditions for the 2D domain .................................................................. 46
Figure 3-5 – Using prism layers near walls to resolve turbulent flows (CD-adapco 2016) ......... 47
Figure 3-6 – Plot of u+ vs y+ for a typical boundary layer (Wilcox 2006) .................................. 49
Figure 3-7 – A big jump in transition from prism layers to polygonal mesh (poor transition) .... 50
Figure 3-8 – An example of a good transition from prism layers to polygonal mesh .................. 50
xiii
Figure 3-9 – Mesh configuration for 2D simulation of flow past a stationary cylinder ............... 51
Figure 3-10 – Details of mesh for flow past a stationary cylinder ................................................ 52
Figure 3-11 – Lift coefficient for cases with or without steady run initiation .............................. 56
Figure 3-12 – A typical plot of drag coefficient in time for 2D simulation of a stationary cylinder
....................................................................................................................................................... 57
Figure 3-13 – A typical plot of lift coefficient in time for 2D simulation of a stationary cylinder
....................................................................................................................................................... 57
Figure 3-14 – Comparison between 2D RANS and 3D LES; on left, 2D RKE with 15,778 cells,
on right, 3D LES with 3.78 million cells (Thingbø 2013) ............................................................ 62
Figure 3-15 – Oscillating lift and drag for cylinders at a subcritical Reynolds number
(Humphreys 1960) ........................................................................................................................ 63
Figure 3-16 – Lift (…) and drag (―) coefficients time history for Re = 3,900 using 2D DNS
(Singh and Mittal 2005) ................................................................................................................ 64
Figure 3-17 – Force and energy-transfer time traces for A
*
= 0.8 and λ
*
= 5.6 (2P - 2P
O
region)
from Morse and Williamson (2009) ............................................................................................. 66
Figure 3-18 – 2P
O
shedding mode for A
*
= 0.8 and λ
*
= 5.6 from Morse and Williamson (2009) 67
Figure 3-19 – 2D domain boundary conditions and velocity field for simulating the moving
cylinder ......................................................................................................................................... 69
Figure 3-20 – Lift coefficient for Moving_12_00_00 (left) and Moving_26_00_00 (right) ........ 71
Figure 3-21 – Vortex shedding pattern for Moving_12_00_00 (left) and Moving_26_00_00
(right) ............................................................................................................................................ 71
Figure 3-22 – Lift coefficient and vortex pattern for Moving_27_RKE case (RKE with Curvature
Corrections)................................................................................................................................... 72
Figure 3-23 – Lift coefficient and vortex pattern for Moving_27_AKNKE case (AKNKE Low
y
+
) .................................................................................................................................................. 72
Figure 3-24 – Lift coefficient and vortex pattern for Moving_27_KW_00 case (KWSST All y
+
)
....................................................................................................................................................... 72
Figure 3-25 – Lift coefficient and vortex pattern for Moving_27_KW_01 case (SKW All y
+
) .. 73
Figure 3-26 – Lift coefficient and vortex pattern for Moving_27_KWSST case (KWSST Low y
+
with Curvature Correction) ........................................................................................................... 73
Figure 3-27 – Lift coefficient and vortex pattern for Moving_27_RSM case (RSM All y
+
with
Linear Pressure Strain Correlation) .............................................................................................. 73
Figure 4-1 – Mode A (left) and B (right) of transformation to 3D vortex shedding (Williamson,
1996) ............................................................................................................................................. 77
xiv
Figure 4-2 – Stages of transition from the primary von Karman vortex shedding to a 3D flow
(Williamson, 1996) ....................................................................................................................... 78
Figure 4-3 – Planar mesh configuration used for 3D simulation of flow past a cylinder ............. 80
Figure 4-4 – Domains with thickness of 3D (36 layers) on left and with thickness of 1D (8
layers) on right .............................................................................................................................. 82
Figure 4-5 – Drag (left) and lift (right) coefficients for the case of 3D simulation by the Standard
k – ε with the thickness of 3D and using 36 elements in z ........................................................... 84
Figure 4-6 – Vorticity (left) and velocity (right) magnitudes in the depth of domain for the case
of 3D simulation by the Standard k – ε with the thickness of 3D and using 36 elements in z ..... 85
Figure 4-7 –Stream-wise vorticity for the case of 3D simulation by the Standard k – ε with the
thickness of 3D and using 36 elements in z .................................................................................. 85
Figure 4-8 – The lift coefficient resulted from the 3D (on left) and 2D (on right) simulations by
the Standard k – ε .......................................................................................................................... 85
Figure 4-9 – Drag (left) and lift (right) coefficients for the case of 3D simulation by the k – ω
SST Low y
+
treatment with the thickness of 3D and using 36 elements in z ............................... 86
Figure 4-10 – Vorticity (left) and velocity (right) magnitudes in the depth of domain for the case
of 3D simulation by the k – ω SST Low y
+
treatment with the thickness of 3D and using 36
elements in z ................................................................................................................................. 86
Figure 4-11 –Stream-wise vorticity for the case of 3D simulation by the k – ω SST Low y
+
treatment with the thickness of 3D and using 36 elements in z .................................................... 87
Figure 4-12 – The lift coefficient resulted from the 3D (on left) and 2D (on right) simulations by
the k – ω SST Low y
+
treatment ................................................................................................... 87
Figure 5-1 – Circular domain size used for simulation............................................................... 106
Figure 5-2 – Boundary conditions of the simulation domain – top view ................................... 106
Figure 5-3 – Boundary conditions of the simulation domain – 3D view.................................... 107
Figure 5-4 – Structured mesh in the r-z plane with a geometric growth in r .............................. 108
Figure 5-5 – Planar mesh in the r-θ plane ................................................................................... 110
Figure 5-6 – Details of mesh for flow past a stationary cylinder with a structured 3D mesh .... 111
Figure 5-7 – Three dimensional structured mesh for simulating flow past a cylinder ............... 111
Figure 5-8 – Different in lift coefficient time trace resulted from DES between a case starting
with a steady RANS run and a case starting from still conditions.............................................. 113
Figure 5-9 – The lift coefficient resulted from the DES S-A with 351,000 cells (on left) and the
3D S-A URANS Low y
+
treatment with 358,410 cells (on right) .............................................. 114
Figure 5-10 – Lift and drag coefficient resulted from DES with the Standard Spalart – Allmaras
..................................................................................................................................................... 115
xv
Figure 5-11 – A section of velocity and vorticity magnitudes resulted from DES with the
Standard Spalart – Allmaras ....................................................................................................... 115
Figure 5-12 – Spanwise vorticity field resulted from DES with the Standard Spalart – Allmaras
..................................................................................................................................................... 115
Figure 5-13 – Matlab script to calculate the statistical values of flow past a stationary cylinder
..................................................................................................................................................... 116
Figure 5-14 – Frequency spectrum of lift force for flow past a stationary cylinder ................... 117
Figure 5-15 – Plot of u
+
vs y
+
for a typical boundary layer (Wilcox 2006) ............................... 121
Figure 6-1 – Map of vortex shedding regimes in Morse and Williamson’s studies ................... 124
Figure 6-2 – Plot of u
+
vs y
+
for a typical boundary layer (Wilcox 2006) .................................. 125
Figure 6-3 – Boundary conditions and the velocity field to model the cylinder oscillation ....... 127
Figure 6-4 – Vector of velocity for flow past a moving cylinder using MRF ............................ 127
Figure 6-5 – Time trace of simulated motion by MRF for a case with an oscillation amplitude of
0.028 m and oscillation period of 2.4 seconds ............................................................................ 128
Figure 6-6 – An example of a PBS file to run a Star-CCM+ file in the batch mode .................. 131
Figure 6-7 – Black x shows the location of simulation points in the λ
*
vs A
*
map .................... 133
Figure 6-8 – Force and energy-transfer time traces for A
*
= 0.8 and λ
*
= 5.6 (2P - 2P
O
region)
from Morse and Williamson (2009) ........................................................................................... 135
Figure 6-9 – Normalized fluid excitation map and selected data points for simulation ............. 135
Figure 6-10 – Dominant frequency of lift for flow past an oscillating cylinder ......................... 137
Figure 6-11 – Total lift coefficient and its phase lag for A* = 0.52 and λ* = 5.9 ...................... 138
Figure 6-12 – Vortex lift coefficient and its phase lag for A* = 0.52 and λ* = 5.9.................... 138
Figure 6-13 – Total lift coefficient and its phase lag for A* = 0.55 and λ* = 6.1 ...................... 139
Figure 6-14 – Vortex lift coefficient and its phase lag for A* = 0.55 and λ* = 6.1.................... 139
Figure 6-15 – Total lift coefficient and its phase lag for A* = 0.6 and λ* = 6.0 ........................ 140
Figure 6-16 – Vortex lift coefficient and its phase lag for A* = 0.6 and λ* = 6.0...................... 140
Figure 6-17 – Total lift coefficient and its phase lag for A* = 0.6 and λ* = 6.25 ...................... 141
Figure 6-18 – Vortex lift coefficient and its phase lag for A* = 0.6 and λ* = 6.25.................... 141
Figure 6-19 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 5.0 ...................... 142
Figure 6-20 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 5.0.................... 142
Figure 6-21 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 5.5 ...................... 143
Figure 6-22 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 5.5.................... 143
xvi
Figure 6-23 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.0 ...................... 144
Figure 6-24 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.0.................... 144
Figure 6-25 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.25 .................... 145
Figure 6-26 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.25.................. 145
Figure 6-27 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.5 ...................... 146
Figure 6-28 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.5.................... 146
Figure 6-29 – Total lift coefficient and its phase lag for A* = 0.78 and λ* = 6.0 ...................... 147
Figure 6-30 – Vortex lift coefficient and its phase lag for A* = 0.78 and λ* = 6.0.................... 147
Figure 6-31 – Total lift coefficient and its phase lag for A* = 0.8 and λ* = 5.0 ........................ 148
Figure 6-32 – Vortex lift coefficient and its phase lag for A* = 0.8 and λ* = 5.0...................... 148
Figure 6-33 – Total lift coefficient and its phase lag for A* = 0.8 and λ* = 5.6 ........................ 149
Figure 6-34 – Vortex lift coefficient and its phase lag for A* = 0.8 and λ* = 5.6...................... 149
Figure 6-35 – Experimental data for total lift and its phase lag for A* = 0.8 and λ* = 5.6 ........ 151
Figure 6-36 – Experimental data for total lift and its phase lag for A* = 0.78 and λ* = 6.0 from
Morse and Williamson (2009a) .................................................................................................. 151
Figure 6-37 – Experimental data for total and vortex lifts and their corresponding phase lags for
A* = 0.8 and λ* = 5.0 Morse and Williamson (2009a) .............................................................. 152
Figure 6-38 – Comparison between CFD simulations and test data of total lift coefficient for A*
= 0.8 and λ* = 5.6 ....................................................................................................................... 152
Figure 6-39 – Comparison between CFD simulations and test data of total lift phase lag for A* =
0.8 and λ* = 5.6 ........................................................................................................................... 153
Figure 6-40 – Comparison between CFD simulations and test data of total lift coefficient for A*
= 0.78 and λ* = 6.0 ..................................................................................................................... 153
Figure 6-41 – Comparison between CFD simulations and test data of total lift phase lag for A* =
0.78 and λ* = 6.0 ......................................................................................................................... 153
Figure 6-42 – Comparison between CFD simulations and test data of vortex lift coefficient for
A* = 0.78 and λ* = 6.0 ............................................................................................................... 154
Figure 6-43 – Comparison between CFD simulations and test data of vortex lift phase lag for A*
= 0.78 and λ* = 6.0 ..................................................................................................................... 154
Figure 6-44 – Comparison between CFD simulations and test data of total lift coefficient for A*
= 0.8 and λ* = 5.0 ....................................................................................................................... 154
Figure 6-45 – Comparison between CFD simulations and test data of vortex lift coefficient for
A* = 0.8 and λ* = 5.0 ................................................................................................................. 155
Figure 6-46 – Normalized fluid excitation map and selected data points for simulation ........... 156
xvii
Figure 6-47 – Time trace of normalized fluid excitation for some of the selected data points .. 157
Figure 6-48 – CFD results of normalized fluid excitation (C
L
sinϕ) for the selected data points 158
Figure 6-49 – Total lift for a cut of A* = 0.75 ............................................................................ 159
Figure 6-50 – Total lift for a cut of λ* = 6.0 ............................................................................... 160
Figure 6-51 – Total lift for A* = 0.78 and λ* = 6.0 resulted from DES SA, DES EB k-ε, and LES
WALE models on different mesh configurations ....................................................................... 162
Figure 6-52 – Effects of 1
st
cell thickness on total lift for A* = 0.78 and λ* = 6.0 resulted by LES
..................................................................................................................................................... 164
Figure 6-53 – Effects of 25% increase in z discretization on total lift for A* = 0.78 and λ* = 6.0
resulted by LES ........................................................................................................................... 165
Figure 6-54 – Total lift for A* = 0.78 and λ* = 6.0 resulted by 2.6 million and 4 million cells 166
Figure 6-55 – Total lift for A* = 0.78 and λ* = 6.0 resulted by different mesh configurations . 166
Figure 6-56 – Temporal resolution effects on the total lift for A* = 0.78 and λ* = 6.0 ............. 167
Figure 6-57 – Total lift resulted MRF vs morphing for A* = 0.7 and λ* = 6.0 with
LES_14_WALE .......................................................................................................................... 168
Figure 6-58 – Effects of morphing field type on total lift for A* = 0.7 and λ* = 6.0 with
LES_14_WALE .......................................................................................................................... 169
Figure 6-59 – Effects of convective scheme on total lift for A* = 0.75 and λ* = 5.5 by LES
WALE and LES Smagorinsky with 4 million cells .................................................................... 170
Figure 6-60 – Total lift for A* = 0.8 and λ* = 5.6 resulted from a DES simulation with 2.6
million cells and MRF vs an LES simulation with 4 million cells and morphing ...................... 171
xviii
List of Tables
Table 1-1 - Different parameters in vortex induced vibration ........................................................ 3
Table 3-1 – Experimental values for drag, lift, and Strouhal number – Re = 4,000 ..................... 43
Table 3-2 – Properties of model for a stationary cylinder 2D simulation .................................... 54
Table 3-3 – Results for 2D simulation of the stationary cylinder with RANS ............................. 59
Table 3-4 – Results for flow past a cylinder at Re ≈ 4,000 by different CFD studies .................. 65
Table 3-5 – Results for flow past a moving cylinder in a 2D domain with RANS models .......... 70
Table 4-1 – Properties of the planar mesh for 3D simulation with RANS ................................... 80
Table 4-2 – Results of 3D and corresponding 2D simulations for different turbulence models .. 88
Table 4-3 – Summary of models with or without triggering 3 dimensionality in a 3D domain ... 90
Table 5-1 – An example of calculations for the 3D mesh .......................................................... 109
Table 5-2 – Details of mesh parameters for flow past a stationary cylinder .............................. 118
Table 5-3 – Parameters of simulations and final results for flow past a stationary cylinder ...... 119
Table 5-4 – Comparison between the current study results and model parameters with most
recently CFD studies ................................................................................................................... 122
Table 6-1 – Comparison of computational speed for different MPI flags .................................. 132
Table 6-2 – Data points used in simulations for wavelength and amplitude of oscillation ........ 134
Table 6-3 – Mean value of C
L
sinϕ for selected data points ........................................................ 156
Table 6-4 – Mesh details for flow past a moving cylinder with A* = 0.78 and λ* = 6.0 ........... 163
Table 6-5 – Details of the coarsest and most refined mesh configurations used for A* = 0. 8 and
λ* = 5.6 ....................................................................................................................................... 171
1
Chapter 1- Background and
Literature Review
1-1- Introduction
Vortex Induced Vibration (VIV) is the phenomenon of vibration of a structure due to its
interaction with a surrounding fluid. This vibration is caused by vortices formed behind the
structure which are shed with some specific frequency. When the vortex generation frequency
matches the structure’s natural frequency, the structure starts to vibrate. In most cases, for
structures with a symmetric shape (like cylinders), VIV is self-excited and self-limiting (Blevins
1990), meaning that the oscillation eventually tends to a stable vibration with some finite
amplitude. This is because of the inherent damping in the system that prevents the vibration
amplitude to grow up.
Even though the amplitude of vibration is limited and finite, VIV is still a big concern in
many structural systems such as offshore structures, suspended bridges, and power plant towers.
VIV damages to structures are usually in form of fatigue in structural components and its effects
are investigated carefully in design and proper mitigation methods may be implemented if the
expected behavior from the structure is not fulfilled.
Vortex shedding from solid objects, as discussed, has the potential to be a hazard to structural
systems. Probably, the collapse of Tacoma Narrow Bridge in 1940 is the most famous incident
which was caused by excessive motions of the deck due to wind and structure interaction and
vortex shedding.
2
Figure 1-1- Tacoma Narrow Bridge collapse in 1940 (picture from www.exit133.com)
Because of its importance and complexity, VIV has been getting a lot of attention in research
and industry. VIV research literature has focused more on circular cylinders as it is a real case in
many structures like offshore risers and it is still a challenging problem. Depending on flow
characteristics, the fluid separates the cylinder at different locations and this alters the flow
pattern behind the cylinder and formation of vortices. Behind the cylinder, there exist a complex
system of detachments – reattachments, interactions between vortices with different scales, and
flow structures in span-wise direction. These all add to the complexity of the system.
1-2- Literature Review
In this research, the focus is on the case where there is only one cylinder which is free to
move only in the direction perpendicular to the flow (cross-flow). This is a fundamental case for
studying vortex shedding, and a lot of experimental, analytical, and computational studies have
been done on it. For studying cylinders vortex induced vibration, several parameters should be
considered. Probably the most important parameter is the Reynolds number which shows the
ratio of inertial forces to viscous ones. For circular cylinders, Reynolds number is written in the
form of:
3
μ
D U ρ
e R ( 1-1)
Where, U is the flow velocity, ρ is the fluid density, D is the cylinder diameter, and µ is the
fluid dynamic velocity.
Other parameters include reduced velocity, frequency ratio, amplitude ratio, mass ratio,
damping ratio, Strouhal number, and fluid force coefficients which will be introduced and
discussed later.
Table 1-1 - Different parameters in vortex induced vibration
Reduced velocity
*
U
D f
U
U
N
*
Frequency ratio
*
f
N
v
f
f
f
*
Amplitude ratio
*
A
D
A
A
*
Mass ratio
*
m
4 /
2
*
L D
m
m
Damping ratio
A
m m k
c
2
Strouhal number St
U
D f
St
v
Force coefficient C
DL U
F
C
2
2
1
In the table, f
N
is the natural frequency of the cylinder, f
v
is the vortex shedding frequency, A
is the oscillation amplitude, m is the cylinder mass, L is the cylinder length, c is the structural
damping, m
A
is the added mass, k is the structural stiffness, and F is the fluid forces on the
cylinder.
Since vortex shedding is a complex fluid mechanics problem and it usually involves high
turbulent flow structures, analytical approaches are limited to very simplified cases. For
4
example, Iwan and Blevins (1974) proposed an analytical formulation for predicting the response
of fixed and forced elastically mounted elements where they introduced a hidden variable for the
fluid dynamic effects. Their model uses a form of Van der Pol equation and its parameters need
to be determined from experimental data. Their model predicts the peak amplitude in free
vibration relatively well in the case of a cylinder, but at a different frequency ratio. It also does
not simulate lock-in, which is observed in experiments.
1-2-1- Experimental Studies
Several experimental studies have been conducted on different aspects of flow past a
cylinder. Experiments show that the flow and vortex pattern behind the cylinder, amplitude and
frequency of vibration, and forces acting on the cylinder mainly depend on Reynolds number.
Considering a stationary cylinder in a uniform flow field, vortices start to generate for Re > 5
(Blevins 1990). Williamson (1996) explains different processes of vortex shedding as the
Reynolds number increases: First, two stationary and stable vortices are formed in the range 5 <
Re < 49. Vortices start to move downstream for Re > 49 and up to Re < 260 vortex generation is
mainly 2 dimensional and the flow is laminar. The regime for 190 < Re < 260 is called wake
transition where, first, vortex loops start to develop and after that fine scale stream-wise vortices
are formed. As the Reynolds number gets bigger, the vortex wake becomes turbulent and this
causes formation of fine-scale vortex structures which have stream-wise components. At around
Re = 1,000, secondary (Kelvin-Helmholtz) vortices start to grow and vortex shedding is
completely 3 dimensional.
5
Figure 1-2 – Regimes of fluid across smooth circular cylinders (Lienhard 1966) – Figure by MIT OCW
Up to Re = 2×10
5
, the boundary layer around the cylinder is still laminar and vortex wake is
fully turbulent. This range is called subcritical. At Re = 2×10
5
, the boundary layer becomes
unstable and on one side of the cylinder becomes turbulent while the other side is still laminar.
This is the case up to Re = 3.5×10
6
(Lienhard 1966), and this range is called the drag crisis as a
significant reduction occurs in fluid forces and the base pressure. In the range of Re > 3.5×10
6
,
which is the supercritical regime, the boundary layer is completely turbulent and fluid forces
start to grow with velocity.
6
Figure 1-3 – Vortex pattern for different vortex regimes (Williamson, 1996)
Forces acting on the cylinder are usually described by their components; drag which is the
component in the direction of the flow, and lift which is the cross flow force component. To have
dimensionless quantities, force coefficients can also be used instead of the force values.
DL U
F
C
Drag
D
2
2
1
( 1-2)
DL U
F
C
Lift
L
2
2
1
( 1-3)
7
Where, F
Drag
and F
Lift
are the in-flow and cross-flow forces acting on the cylinder,
respectively.
Figure 1-4 – Components of fluid excitation on a cylinder
If the cylinder is free to oscillate, at lower velocities, the oscillation amplitude is small and
negligible and the vortex shedding frequency increases with velocity similar to the case of a
stationary cylinder. However, when the shedding frequency gets close to the cylinder natural
frequency, the amplitude grows rapidly. In fact, experiments show that, depending on Reynolds
number and mass-damping ratio, for frequency ratio around unity, there can exist a sudden jump
in amplitude which corresponds to some changes in shedding patterns (depending on mass-
damping ratio). Interestingly, with the further increase in velocity, the frequency ratio does not
change for some range of velocity. This phenomenon is called lock-in meaning that vortex
generation synchronizes itself with the cylinder oscillation. But unlike the frequency, the
amplitude drops with further increase in velocity.
Studying these changes in flow and amplitude in more details reveals some interesting
phenomena. These changes are corresponding to different branches in the amplitude vs velocity
response curves. Khalak and Williamson (1999) observed that if the mass-damping ratio is low,
there are 3 branches for free vibration, which they call them initial, upper, and lower branches.
Bishop and Hassan (1964) and Feng (1968) showed that if the mass-damping ratio is high, the
8
free vibration response has only 2 branches (initial and lower branch). Switches between these
branches have different characteristics, for example Khalak and Williamson (1999) showed that
the switch between the initial and upper branches has a hysteresis behavior, while they observed
an intermittent switch between the upper and lower branch.
Figure 1-5 – Response amplitude and frequency for low mass damping. (a) Amplitude response for low
(■), and high (◊) mass damping, (b) Oscillation frequency through the excitation regime in low mass
damping (Khalak & Williamson, 1997)
9
Different branches can have different vortex generation patterns and lift characteristics as
shown by experiments (Govardhan & Williamson 2000). It has been shown that, for low mass
damping, φ
vortex
which is the phase between the cylinder oscillation and the force due to the
vortex shedding, is continuous during the transition between upper and lower branches while
φ
vortex
experiences a jump at the initial – upper branch transition which should be associated with
change in the vortex shedding mode.
A part of the literature is on study of vortex shedding when the cylinder is under a controlled
motion with a prescribed amplitude and frequency (Controlled Vortex Induced Vibration).
Unlike, the freely oscillating case, here, the oscillation is imposed to the cylinder and the vortex
structures and forces on the cylinder are measured. Although this is not the case in most
applications but there are several reasons why studying controlled VIV is beneficial. Controlled
vortex shedding is obviously an easier path since it has less controlling parameters and studies
have shown that if it is conducted well, it is capable of predicting the freely VIV result. As
Bearman (2009) explained, on one hand, in the freely oscillating problem, a small change in flow
velocity can lead to big differences in amplitude of motion and also change in shedding pattern.
In fact, the peak amplitude of upper branch depends on the Reynolds number (Govardhan &
Williamson 2006). These fluctuations are hard to study under conditions of free vibration. On the
other hand, the importance of flow history and the role of harmonics in the displacement trace of
a freely oscillating cylinder should be studied well. Therefore, if controlled VIV is opted for
studying the freely oscillating case, a large number of test cases are required in order to predict
the conditions under which fluid excites the cylinder and transfers energy to the cylinder.
Morse and Williamson (2009) at Cornell University performed an extensive series of tests on
vortex shedding under controlled motion. About 6,000 test cases were done in a continuously
10
flowing tank with a refined resolution and because of their resolution they have come up with
some interesting results. First, they were able to predict the freely oscillating case by their
controlled VIV data with an acceptable accuracy. They also justified different phenomena seen
in free oscillation such as hysteresis and intermittent switches between different response
branches.
In 1988, Williamson and Roshko did a similar study with a much less number of test cases
and they showed that, depending on prescribed amplitude and frequency of motion, there exist
different vortex shedding patterns (or modes) behind the cylinder. Basically, these modes are 2S,
2P, C(2S), and P+S. In 2S, two (single) vortices are shed in one oscillation, while 2 pairs of
vortices are formed per oscillation if the 2P mode is observed. In P+S, one single and one pair
are shed, and finally, C(2S) is just the coalescing 2S mode. These modes have been confirmed in
free vibration.
Figure 1-6 – Williamson and Roshko’s regimes of vortex shedding
Morse and Williamson (2009) worked on controlled vibration vortex shedding with a higher
Reynolds number compared to the Williamson and Roshko’s studies (4,000 vs 1,000) and
11
conducting about 6,000 test cases enabled them to have a remarkable resolution which leads to a
precise study of different phenomena observed or expected in free oscillation.
First of all, they observed a new stable shedding mode which they call it 2P
O
(or 2P
Overlap
). In
this mode, 2 pairs of vortices are shed per oscillation (like 2P), but unlike 2P, one of the vortices
is smaller and dissipates faster. In studying the possible modes for free oscillation, they saw that
2P
O
has the highest motion amplitude. They also showed that it is possible to have this mode for
upper branch of free oscillation response. Thus, it is an important mode.
Figure 1-7 – Map of vortex shedding regimes in Morse and Williamson’s studies
Morse and Williamson have used the normalized wave length (λ
*
) in their studies instead of
reduced velocity (U
*
).
*
*
*
f
U
D f
U
( 1-4)
In controlled vortex shedding f is the frequency of oscillation.
12
It should be noted that it may not be possible to discover an overlap region in free vibration
experiments as 2P and 2P
O
modes yields different values of fluid excitation. This mode was
already observed in Govardhan and Williamson (2000) free oscillation diagram.
Figure 1-8 – Vortex shedding modes observed in Morse and Williamson’s experiments
According to the equation of motion for the cylinder, if the motion is sinusoidal, the resulting
force (lift) will also have a sinusoidal force with the same frequency as the motion, but there will
be a phase lag in the response.
The equation of motion for the cylinder is as follows
) (t F ky y c y m
( 1-5)
When the cylinder motion is synchronized with vortex shedding, the cylinder motion and the
cross-flow force acting on it can be approximated with sinusoidal forms Morse and Williamson
(2009) :
) 2 ( sin t f A y
( 1-6)
13
) 2 ( sin ) (
1
t f F t F ( 1-7)
The phase lag, φ, is shown to be an important parameter in free vibration as when the initial
branch switches to the upper branch, there will be a significant change in the phase lag.
Khalak and Williamson (1999) showed that, for free oscillation, the steady state response has
a form like
,
sin
4
1
*
2
*
*
* 3
*
f
f
U
C m
C
A
A
L
( 1-8)
EA
A
C m
C m
f
*
*
*
( 1-9)
Where C
A
is the added mass coefficient and is equal to 1.0 for circular cylinders. C
EA
is the
effective added mass coefficient due to the transverse force in phase with the cylinder
acceleration.
2
*
*
* 3
cos
4
1
f
U
A
C
C
L
EA
( 1-10)
The energy transferred from the fluid (E
IN
) and the energy dissipated by the cylinder damping
(E
OUT
) per one oscillation are defined as
sin
1
AF E
IN
( 1-11)
f cA E
OUT
2 2
2 ( 1-12)
If the oscillation is steady, the input and output energy are equal. By normalizing these two
equations we will have
*
2
*
*
* * 3
4
sin
f
f
U
C m A
C
A
L
( 1-13)
14
The phase lag, ϕ, as discussed, has an important role in defining the free vibration response.
According to Equation 1-13, sin
L
C is the fluid excitation or normalized input energy over an
oscillation to the system. So, if the phase lag lies between 0 and 180 degrees, this quantity is
positive and there should exist the possibility of free vibration.
Mass damping is also an important parameter in free vibration. As Brika & Laneville (1993)
found, for high mass damping, two free vibration branches, namely initial and lower, are
available, while according to Khalak and Williamson (1999), in the case of low mass damping,
the cylinder can also experience a higher amplitude mode (upper branch) during free vibration.
Comparing the prediction of free vibration by controlled VIV data with the experiment
shows a good agreement for both low and high mass damping. This is due to the fine resolution
of test cases and careful matching of all conditions (especially the Reynolds number) in free and
controlled vibration studies.
15
Figure 1-9 – Measured (○) and predicted (●) amplitude response for free VIV at (a) low mass-damping
and (b) high mass-damping (Re = 4,000)
By introducing the concept of ‘energy portrait’, which is the plot of energy excitation (E
*
in
)
and the energy dissipated by the structural damping (E
*
out
) as a function of A
*
while U
*
is kept
constant, Morse and Williamson were able to explain different phenomena observed in free
vibration.
A stable condition for equilibrium is where E
*
out
intersects with E
*
in
curve and the slope of
E
*
in
curve is negative. So, there can be different paths for getting stable conditions as U
*
increases or decreases which can lead to a hysteresis between the initial and the upper branches.
There can be situations where two modes are available for the same U
*
(2P
O
overlaps 2P and
2S regions), so if the shedding mode changes there should a change in amplitude (increase or
16
decrease) to get to a new stable point. This is why the amplitude changes intermittently between
the upper and the lower branches.
Finally, Morse and Williamson proposed a region where free vibration is possible for Re =
4,000. It should be noted that the overlap region between 2P
O
and 2S is, although leads to a
positive response, excluded because the switch between these two modes is not stable and does
not produce a steady state.
Figure 1-10 – Map of the zone where free vibration can occur in the amplitude – velocity plane for
m
*
=10
1-2-2- Numerical Studies
With developing computational facilities, there has been an increase in utilizing numerical
simulations to study flow past a cylinder. Most of these studies have been aimed to validate the
new numerical techniques or codes for the case of a challenging flow that involves curved
streamlines, separation, detachment and reattachment, and unsteadiness. But there has been
numerical research for a better understanding of vortex structures and the transition from 2D to
3D as well.
17
Since flow past a cylinder involves many complex flow structures, it is very crucial how to
resolve the near wall behavior numerically. There are many approaches for near wall modeling
and depending on the flow Reynolds number one or more methods will be suitable. Because of
the turbulent nature of flow past a cylinder for Re > 300, turbulence modeling and how to
simulate or model the turbulence scales will be important. There are several models, methods,
and techniques for turbulence modeling which will be discussed later.
In this section, some previous numerical research studies and their results for the stationary,
freely oscillating, and controlled vortex shedding cases are addressed.
For the stationary cylinder case, up to date (Kara et al., 2015), due to high computational
costs, Re = 10,000 is the highest Reynolds number that Direct Numerical Simulation (DNS) has
been used for,. The DNS studies of Dong and Karniadakis (2005) had good agreement with the
experimental tests regarding the statistics of flow parameters and vortex structures in the near
wake. Breuer (2000) conducted Large Eddy Simulation (LES) for Re = 140,000 and found
encouraging agreement for this high velocity. Norberg (2003) listed several numerical studies’
result and methodology done by different researchers at different Reynolds number with a focus
on the lift coefficient. There have been also several studies following the Williamson’s
experimental work on onset of three dimensionality in vortex shedding published in 1988.
Pereira et al. (2015) studied flow past a cylinder for Re = 3,900 with different Reynolds
Averaged Navier-Stokes (RANS), LES, and hybrid LES/RANS models in 2D and 3D domains
and concluded that 3D simulations with LES or hybrid models lead to a better agreement.
According to their simulations, when modeling with RANS, 3D simulations can improve the
results up to 50%. Prsic and Myrhaug (2012) did a similar study for Re = 3,900 and 13,100 and
proposed that near wake mesh resolution is the most influential parameter in predicting the flow
18
forces on the cylinder correctly, and refinement along the cylinder does not affect the results
significantly.
CFD simulation of Vortex Induced Vibration (VIV) in circular cylinders has been limited to
low Reynolds numbers due to the limitations of computational resources. Yang et al. (2008)
investigated 2D cylinder VIV for Re = 200 and observed some type of branching behavior that
had been seen experimentally in higher Re. Singh and Mittal (2005) studied hysteresis response
of the kind observed in experiments. They also investigated a range of Re from 50 to 500 with
2D simulations and observed P+S mode in the upper branch for Re > 300 which was the first
time reported in simulations or experiments. Willden (2006) did 2D LES for Re = 10,000 and
observed 2P and P+S modes for the lower branch and 2S for the initial branch. Since the
agreement between their results and the experiment was poor, they concluded that 2D modeling
is not able to correctly simulate the upper branch which had already been known that involved
three dimensionality (Hover et al., 2004).
Navrose and Mittal (2013) used 3D DNS to simulate free vibration for Re = 1,000 and they
observed 2S mode along the upper branch for the first time. They also studied hysteresis and
intermittent switches between branches and illustrated that for this range of Re, the onset of lock-
in/synchronization is at the point where initial branch switches to the upper branch. They also
showed that the phase lag between the lift and the cylinder oscillation is constant through the
initial and the upper branch and there is a 180 degree jump in the upper-lower transition.
Meneghini and Bearman (1995) conducted controlled VIV with high amplitudes for Re =
200 and they particularly focused on lock-in region. Udaykumar et al. (2001) replicated the same
19
study to show their sharp-interface IBM method’s capabilities. In both of these studies, P+S
mode was observed for high amplitudes. Dong and Karniadakis (2005) performed 2D DNS
which had a very good agreement with the experiment, but their study was limited in terms of
number of cases and time of simulation because of computational expenses. Nguyen and
Temarel (2014) used 2D RANS models for simulating flow past a cylinder at Re = 10,000 for a
stationary case and also forced and free oscillations. They finally concluded that 2D RANS
simulations are capable of predicting the overall behavior of phenomena observed in different
cases and the agreement with the experimental data is only qualitative. These simulations
particularly fail to reproduce the experimental results regarding the lift coefficient.
Kim et al. (2014) studied forced oscillation of a cylinder at Re = 5,500 with 3D LES and
focused on a narrow range of frequency ratio in order to investigate the transitions between
modes and vortex patterns. They observed significant changes in the lift coefficient properties
after passing a certain frequency ratio; the lift amplitude got bigger and it became in-line with the
cylinder oscillation. They also showed that for Re = 10,000 , at a frequency close to the von
Karman frequency, flow properties change and lift may experience a drop or an increase
depending the oscillation amplitude.
Kara et al. (2015) used the immersed boundary method (IBM) with a fluid solver which was
strongly coupled with the equation of motion through Hamming’s 4th order predictor-corrector
method. They studied flow past a stationary cylinder and also free and controlled VIV for
different Reynolds numbers in the laminar region (80 < Re < 185) and got good agreement with
experimental data.
20
Chaplin et al. (2005) conducted a series of studies for a very large length to diameter ratio of
469 by 2D simulation based on strip theory in which it was assumed that during lock-in, the flow
is locally 2D. Strip theory is currently the most practical option for long cylinders (Kara et al.
2015). But, it is shown that in general, CFD overestimates the in-line motion and underestimates
the cross-flow response. The error in the cylinder curvature, which is very important in
predicting the fatigue damages, computed by CFD can be very significant. Newman and
Karniadakis (1997) also studied the span-wise effects in vortex shedding by DNS for different
low Reynolds number both for free and controlled VIV and observed standing and travelling
wave responses along the cylinder.
In Chapter 2, the mathematical formulation implemented in Star-CCM+ will be presented.
Chapter 3 discusses the 2D simulation of a stationary and moving cylinder, their modeling
challenges, results of simulations, and justifications. Chapter 4 covers the 3D simulation of a
stationary cylinder with RANS model and explains the differences between 2D and 3D
simulations.
21
Chapter 2- Numerical Formulations
2-1- Introduction
The focus of this research is on numerical simulation of flow past a controlled moving
cylinder. For this purpose, a commercial fluid mechanics solver called Star-CCM+ is used.
Star-CCM+ is a CFD solver developed by CD-adapco which has been successfully
implemented in many engineering applications. The USC’s water resources engineering group
has used this software for more than 8 years. In fact, Bozorgnia (2012) conducted numerical
studies on wave/deck interaction in highway bridges when they are exposed to hurricanes. The
simulations were based on multi-phase flows in 2D and 3D domains, and he got good agreement
with experimental data. The other recent CFD study of this group was on CFD simulation of
breakwater overtopping by Bozorgnia, Eftekharian, and Lee (2014). They implemented 2D
simulations with RANS turbulence models for analyzing vortex structures when a solitary wave
overtops the breakwater. Comparing the CFD results with the experiments by Zhuang and Lee
(1996) shows the CFD simulations do a good job even for the zone right after the breakwater and
near the seabed where a big rotational flow is developed.
In this section, features and capabilities of the software are addressed. More details can be
found in the Star-CCM+ user manual (CD-adapco 2016) and Bozorgnia (2012).
22
2-2- Governing Equations
Star-CCM+ solves Navier-Stokes equations. According to CD-adapco (2016) the integral
form of these equations are as follows. For the control volume V, bounded by closed surface A,
the continuity equation will be
0
A
g
V
d dV
t
a v v
( 2-1)
And the momentum equation is
V
L u p g r
A A
A
g
V
dV d d p
d dV
t
) ( f f f f f f a T a I
a v v v v
( 2-2)
In these equations, ρ is the fluid density, v the velocity with components of
i
u ,
g
v the grid
velocity, p the pressure, I the identity matrix, T the stress tensor, and f ’s represent different
body forces such as rotation and gravity.
The stress tensor, T , has two components for turbulent flows, laminar and turbulent
t l
T T T ( 2-3)
where,
I v v v T ) (
3
2
) (
T
l
( 2-4)
µ is the fluid dynamic viscosity. If the Boussinesq approximation is used for turbulence
modeling,
t
T will have a form similar to
l
T :
I v v v T ) (
3
2
) (
T
t t
( 2-5)
Here,
t
is the turbulent viscosity of the flow.
23
The governing equations are solved using a segregated flow model. Segregated flow models,
as opposed to coupled flow ones, solve for the velocity components and pressure in an uncoupled
way. These models are suitable for incompressible flows with constant density (CD-adapco
2016). A predictor-corrector approach is utilized to obtain a linkage between the continuity and
the momentum equations. The entire process can be defined as using a co-located variable
arrangement and a Rhie-and-Chow-type pressure-velocity coupling with algorithms like
SIMPLE (CD-adapco 2016). The software has two more solvers for segregated flow model, one
for velocity, which solves the momentum equations to get the intermediate velocity field, and the
other one for pressure which solves the discrete equation for pressure correction and updated the
pressure field (CD-adapco 2016).
2-3- Discretization
The governing equations discussed in the previous section need to be discretized and solved
numerically. Star-CCM+ implements Finite Volume Method (FVM) for this purpose. In Finite
Volume Method, the solution domain is further subdivided into a finite number of small
elements or cells which each one is called a control volume (CV). The Navier-Stokes equations
are applied to each control volume as well as the entire domain. Since surface integrals for inner
faces cancel out, by summing all the equations for all control volumes, a global conservation
equation will be obtained. This leads to a set of linear algebraic equations that has a total number
of unknowns equal to the number of cells.
24
Figure 2-1 – Control Volume for cell (i,j) at time level n, with flux densities H
n
on sides, side vectors S ,
and center data u (from Causon et al., 2011)
It is most convenient to calculate the volume and surface integrals by using the midpoint rule
approximations. In this way, simple algebraic expressions, which are second order, will be
obtained. In the case of control volumes, since the values of the midpoint are calculated, it is
needed to only multiply these midpoint values by the volume of the CV’s to get the values for
the whole system. For the surface integrals, as the value of variables are now known at the center
of the faces, more approximations are required.
2-3-1- Momentum Equation in Discrete Form
Using the momentum equation, Equation 2-2, for the cell-centered CV of cell-0, the
following discretized equation will be obtained (CD-adapco 2016):
f f
f
f
f
g
p V
dt
d
a T a I a v v v v
0
( 2-6)
in which V is the volume of cell-0 and a is the area of face f.
25
For each velocity component an implicit, first or second order in time, linear system of
equations can be used in the discrete form. Here, the transient terms, body forces, and convective
flux are discretized as a transport of a simple scalar quantity. For the stress tensor, T, the velocity
tensor gradient of faces,
f
v , will be computed in terms of cell velocities according to the
following equation (CD-adapco 2016)
ds v v v v
_____ _____
f f f
( 2-7)
where
0 1
v v v ,
2
1 0
_____
v v
v
f
, and
ds a
a
.
0
v and
1
v are the explicitly
computed velocity gradient tensors in cell-0 and -1. ds is the vector between cells centroids. For
boundary faces, the software used available information at the boundary depending on the flow
regime and the type of boundary condition (CD-adapco 2016).
2-3-2- Continuity Equation in Discrete Form
The discrete form of continuity equation is as follows (CD-adapco 2016):
0
*
f
f
f
f
f
m m m
( 2-8)
f m
*
, the uncorrected mass flow rate, is computed after the discretized momentum equations
are solved with the following equation (CD-adapco 2016):
f f f
f G m
2
*
*
1
*
0
v v
a
( 2-
9)
where
*
0
v and
*
v
1
are the cell-0 and -1 velocities obtained after solving the discretized
momentum equation.
f
g f
G v a is the grid flux.
f
is the Rhie-and-Chow-type dissipation at
the face:
26
ds
______
*
f
L R f f
p p p Q
( 2-10)
in which,
a α
1 0
1 0
a a
V V
Q
f f f
( 2-11)
V
0
and V
1
are the volume of cell-0 and -1 respectively.
0
a and
1
a are the average of
momentum coefficients for all components of momentum for cell-0 and -1 respectively.
L
p and
R
p are the pressure values for cell-0 and -1 respectively.
______
*
f p is the volume-weighted average
of the cell gradients of pressure (
2
1
*
0
* ______
*
p p
p f
).
The mass flow correction,
f
m
, is given by the following equation (CD-adapco 2016):
upwind
T
f
f
f f
p
p
m
p p Q m
*
1 0
( 2-12)
Here,
0
p and
1
p
are the cell pressure corrections and
upwind
p
will be:
0 for
0 for
*
1
*
0
f
f
upwind
m p
m p
p
( 2-13)
The discrete pressure correction equation is obtained from Equations 2-8 and 2-13 and can be
written in coefficient form as:
r p a p
n
n n p
( 2-14)
Where r is just simply the residual of mass flow rate into the cell:
f
f m r
*
( 2-15)
27
If the velocity values are specified on the boundaries, like walls, inlets, or symmetry
boundaries, the value of f m
*
is calculated directly with known values of f
*
v on boundaries
(CD-adapco 2016):
f
f
f
f G m
* *
v a
( 2-16)
For this type of boundaries, the Neuman condition is used for pressure
0
p p
f
( 2-17)
and the mass flux corrections are zero.
For specified-pressure boundaries, like pressure outlets, the pressure corrections are not zero
and the uncorrected boundary mass flux is computed as following (CD-adapco 2016)
f f f f
f G m v a
*
( 2-18)
f
v is the boundary velocity and
f
is the dissipation coefficient given as:
___
______
*
0
*
0
*
ds p p p Q
f f f
( 2-19)
where:
a α
__
__
0
0
a
v
Q
f
f
( 2-20)
For subsonic outflows, 0
f
p ,
0
p p
upwind
and Equation 2-12 will become (CD-adapco
2016):
0
p
p
Q m
T
f f f
a v ( 2-21)
28
For subsonic inflows, 0
upwind
p and
f
m
can be written in the following form (CD-adapco
2016)
2
f f f
f f
f
Q
Q
m
v a v
a v
( 2-22)
and
0 2
2
p
Q
Q
p
f f f
f f
f
a v v
v
( 2-23)
2-4- SIMPLE Solver Algorithm
SIMPLE stands for Semi-Implicit Method for Pressure Linked Equation. Many CFD
software packages utilize this algorithm to solve the Navier-Stokes Equations. Star-CCM+
follows the following steps to apply SIMPLE solver for the Navier-Stokes Equations (CD-
adapco 2016):
1) Set the boundary conditions,
2) Compute reconstruction gradient of velocity and pressure,
3) Compute velocity and pressure gradients,
4) Solve the discretized momentum equation to create the intermediate velocity field v*,
5) Compute uncorrected mass fluxes at faces
*
f
m ,
6) Solve pressure correction equation to produce cell values of the pressure correction pʹ,
7) Update pressure field p
n+1
= p
n
+ ωpʹ where ω is the under relaxation factor for
pressure,
29
8) Update boundary pressure correction
b
p ,
9) Correct face mass fluxes m m m
f
n
f
* 1
,
10) Correct cell velocities:
V
p
n
p V
a
v v
* 1
( 2-24)
where, p is gradient of the pressure corrections,
V
p
a
is the vector of central coefficients for
discretized linear system representing the velocity equation and V is the cell volume,
11) Update density due to pressure changes,
12) Free all temporary storage.
2-5- Numerical Method for Solving Algebraic Equations
Applying the finite volume formulation to the Navier-Stokes Equations will lead to the
coefficients of linear equation system that should be solved implicitly. For the transported
variable , the algebraic system for a typical point p at iteration k+1 can be written as (CD-
adapco 2016):
b a a
n
k
n n
k
p p
1 1
( 2-25)
Where the summation is over all neighbor cells of cell p. b is obtained from the previous
iteration step and the coefficients a
p
and a
n
are evaluated directly from the discretized equations.
Introducing a new variable, ω, to the above equation will help to cut out steep oscillations and
avoid big changes in
1 k
p
over an iteration. So, a new equation will be generated (CD-adapco
2016):
30
k
p
p
n
k
n n
k
p
p
a
b a
a
1
1 1
( 2-26)
Here the right hand side is evaluated at the previous time step, k, and acts as a source term. It
would be easier to define a delta form like
k
p
k
p p
1
rather than solving this equation. So,
the final system that needs to be solved becomes (CD-adapco 2016):
n
k
n n
k
p p
n
n n p
p
a a b a
a
( 2-27)
The right hand side of the above equation is the residual term which represents the
discretized form of the original transport equation for the k-th time step.
The system of equations for CFD problems are usually very huge and contains millions of
equations. For solving these systems of linear algebraic equations there are generally two
approaches: Direct methods and indirect or iterative methods. The iterative methods, like Jacobi
or Gauss-Seidal, are based on a repeated application of relatively simple algorithm. They are
usually far more efficient for large sets of equation. However, these methods usually lack fast
convergence rates which can be a big issue for large sets of equations.
Recently, multigrid acceleration techniques have been developed for improving the
convergence rate deficiency of iterative algorithms and currently are implemented in CFD codes.
The concept of multigrid is based on fact that coarser mesh sizes have a more rate of
convergence due to shorter wave lengths for error. So those methods that use the multigrid
concept, first, form a coarse mesh by the finer cells, and transfer the residual from a fine scale to
a coarse level, and finally, transfer back the correction from the coarse level to the original finer
level. More information can be found in CD-adapco (2016), Bozorgnia (2012), and Versteeg and
Malalasekera (2007).
31
2-6- Turbulence Modeling
One of the challenges in simulating any fluid mechanics problem is how to deal with the flow
turbulence. Many engineering applications in fluid mechanics involve some sort of turbulent
flows and it is required to treat turbulence as accurate as possible in order to predict the
phenomena well. The issue is, depending on flow properties and structures, turbulence can be
very costly and even with super computers that have parallel processing capabilities, it is
impossible to model turbulence accurately in most applications. However, there have been
models and techniques developed since 1960s to have an acceptable accuracy depending on
applications and purposes. Here, the methods and technique that Star-CCM+ utilizes to resolve
turbulence are presented.
2-6-1- Direct Numerical Simulation (DNS)
In DNS all turbulence scales are solved explicitly. Although this is the best approach to
simulate the flow characteristics, the computational cost of DNS is extraordinarily high and its
applications have been limited for low velocity cases. As discussed in the previous chapter for
flow past a cylinder Re = 10,000 is the highest value so far that DNS has been used for
simulation. Star-CCM+ is capable of doing DNS in a sense that choosing LES with a very small
sub-grid scale will result in DNS simulation.
2-6-2- Reynolds Averaged Navier-Stokes (RANS)
The idea behind RANS models is decomposing the flow parameters into two components;
one that is the averaged over time, and the fluctuating part:
x x t , ( 2-28)
32
Where
T
dt t
T
0
,
1
x x x
Applying Equation 2-28 into the Navier-Stokes equations and taking average of all equations
will result in a new set of equations that are called Reynolds Averaged Navier-Stokes equations.
The Reynolds averaged continuity equation looks pretty much the same as the original continuity
equation, however, the Reynolds averaged momentum equation has now an extra term:
_____
j i ij
u u R
( 2-29)
R
ij
is called Reynolds stress. In RANS models this term is expressed in terms of averaged
quantities of the flow. There are different types of RANS models and they are usually
categorized by the number of added equations to the system. Zero-equation models relate
Reynolds stress terms to velocity via empirical expressions. Star-CCM+ does not use this type of
turbulence model.
In one- and two-equation models, the Reynolds stress is related to velocity changes by the
Boussinesq eddy viscosity approximation. This approximation is based on the analogy to laminar
flows (Wilcox 2006):
ij ij T j i
k S u u
3
2
2
_____
( 2-30)
In which, S
ij
is the mean-strain rate tensor,
T
is the turbulent (eddy) viscosity of flow, and k
is the turbulence kinetic energy. This equation is sometimes called the (linear) constitutive
relation.
i
j
j
i
ij
x
U
x
U
S
2
1
( 2-31)
_____
2
1
i i
u u k ( 2-32)
33
Different models have different approaches for modeling
T
. Star-CCM+ has 3 different
models for this purpose:
Spalart-Allmaras Models: This is a one-equation model that solves for
T
directly. It was
developed for aerospace engineering applications and since it solves only an extra equation,
computationally it is fast.
k - ε Models: k-ε is the most popular RANS model and has been used in different
applications since 1970s. It solves two transport equations, one for k, the turbulence kinetic
energy, and the other one for ε, the turbulent dissipation. According to Wilcox (2006)
__________
k
i
k
i
x
u
x
u
( 2-33)
υ is the fluid viscosity. Relating
T
with k and ε, each k-ε model variant may use a different
expression. For example, in the Standard version of k - ε model,
T
is computed as following
(Wilcox 2006)
2
k
C
T
( 2-34)
here,
C is a constant and equal to 0.09.
Star-CCM+ has several versions for k - ε model and depending on the application, one or
more options are recommended.
k - ω Models: Developed originally by Wilcox in 1988, it is another two-equation model that
instead of ε, a new variable called the specific dissipation rate, ω, is modeled and solved. Again,
34
different versions have different definitions for ω, but in general ω is proportional to the ratio of
ε and k. The standard k - ω models uses the following expression for ω (Wilcox 2006)
k
( 2-35)
In Star-CCM+ there are two options for k - ω models; the Standard (Wilcox) k - ω and k - ω
SST (Mentor) model (CD-adapco 2016). The k - ω SST is basically a model that blends k - ω
and k - ε models; it is a k - ω model very close to wall boundaries and gradually transforms to a k
- ε model for areas farther away from the wall.
Reynolds Stress Models (RSM): In RSM’s, instead of using the Boussinesq approximation
and turbulent viscosity, all components of the Reynolds stress tensor in Equation 2-29 are solved
individually by corresponding transport equations. Since the Reynolds stress matrix is
symmetric, this means that three and six equations for Reynolds stress components are required
in 2D and 3D respectively. However, modeling only the Reynolds stresses is not a complete
solution and one more equation is necessary to close the system. Star-CCM+ uses the ε equation
from the Standard k - ε model for this purpose.
Reynolds Stress Models have been developed to improve the shortcomings of one- and two-
equation models more specifically the Boussinesq approximation deficiencies. Some of the areas
that the Boussinesq eddy viscosity approximation falls back include flow with sudden changes in
mean strain rate, flow over curved surfaces, flow in duct with secondary motions, flow in
rotating fluids, and three-dimensional flows (Wilcox 2006). It should be noted that, using
Reynolds Stress Models is not the only solution to these problems as there have been a lot of
efforts for improving the behavior of two equation models and some are implemented in Star-
35
CCM+ such as nonlinear constitutive relation, flow curvature correction, and considering
secondary flows.
2-6-3- Large Eddy Simulation (LES)
The idea of using large eddy simulation is coming from the fact that small scales of
turbulence have sort of universal characteristics, so they can be modeled, but the rest of scales,
which they are now called large, should be simulated. For this purpose, a cut-off scale is defined
as a sub-grid scale (SGS) and any scales smaller than SGS will be modeled. This works like a
filter but here, unlike RANS where the filtering is in time, the turbulence scales are filtered.
Figure 2-2 – Scales of turbulence that are resolved or modeled in different methods (Sodja 2007)
Even with using the sub-grid scale which eliminates the need of so many very small-size
cells near walls in DNS and allows to use a bigger time step, LES is still computationally very
expensive and is not a very common tool in industry. Moreover, LES requires more details and
caution regarding the boundary and initial conditions compared to RANS models.
36
Figure 2-3 – Turbulence energy spectrum for LES model. Wavenumbers greater than k2 need to be
modeled (from Apte)
As discussed above, in LES, the governing equations are filtered in space instead of time.
This means that any flow quantity can be decomposed into two terms (
~
and )
~
( 2-36)
Where
~
is the filtered quantity and is the sub-filtered (or sub-grid) quantity.
Using this decomposition in the Navier Stokes equations and applying a filter into the
equations will result in a set of equations much like RANS, but here, instead of Reynolds
stresses, there exist new terms that are called sub-grid scale stresses. Sub-grid scale stresses
represent the interaction between large eddies (which are resolved) and small eddies (unresolved
ones). Using a Boussinesq approximation,
t
T can be written as:
I v S T k
t t t
~
3
2
2 ( 2-37)
Where S is the mean strain rate tensor for the resolved velocity field, v
~
is the filtered
velocity, ρ is the density, and k is the sub-grid scale turbulent kinetic energy.
37
S has a form similar to RANS formulations:
T
v v S
~ ~
2
1
( 2-38)
The sub-grid scale turbulent viscosity μ
t
is required to close the governing equations. It
should be noted that μ
t
is described by the Sub-grid Scale Model (SGS) used in the model and it
represents the effects of small eddies on large resolved structures. Star-CCM+ has three different
sub-grid scale models that are briefly described here:
Smagorinsky Sub-grid Scale
The Smagorinsky Sub-grid Scale model was the very first sub-grid scale model proposed by
Smagorinsky (1963), and it uses a mixing length hypothesis to model the sub-grid scale stresses:
S
2 2
S t
C ( 2-39)
The model is based on the local equilibrium balance between shear production and
dissipation in the SGS turbulent kinetic energy (k). The constant
S
C is not universal and this is
one of shortages of this model. Nonetheless, the value of
S
C is usually recommended to set to
0.1. The other feature the Smagorinsky’s model is that near walls it requires a damping function
to reduce the value of
S
C near wall boundaries. A very typical damping function is the original
the Van Driest damping function. According to CD-adapco the Van Driest damping function is a
non-local operation and it needs the use of a KD-tree data structure to communicate wall
information to the interior. More importantly for flows with complex structures, this KD-tree
must be broadcast to each parallel computing node, resulting in a memory and communication
overhead (CD-adapco). This would mean that for complex flows, a very large computational
resource is required if the Smagorinsky model with the Van Driest damping function is
implemented.
38
Dynamic Smagorinsky Subgrid Scale
In the Dynamic Smagorinsky Sub-grid Scale method, unlike the original Smagorinsky
model, the constant
S
C is replaced by a local coefficient which changes in time and it is
computed by test-filtering the flow quantities for any length scale greater than the grid length
scale (CD-adapco). This model eliminates the need of damping functions for wall bounded
flows.
WALE Subgrid Scale
The WALE (Wall-Adapting Local-Eddy Viscosity) Sub-grid Scale model is the most recent
development in sub-grid scale model. The model uses a form of velocity gradient tensor in its
formulations and similar to the Smagorinsky model, the WALE model coefficient is not
universal. However, according to CD-adapco, test validations have shown that the WALE model
is less sensitive to the value of the model coefficient and like the Dynamic Smagorinsky Sub-
grid Scale model it does not require any damping functions for adjustments near wall boundaries
(CD-adapco).
Resolving Convection Terms
One important point in using LES models is how to discretize convective terms. Any
convective term can be discretized at a face and written as:
f f f f
m m a v ( 2-40)
39
Here
f
m
is the mass flow rate at the face, a the area of the face, and
f
is the fluid quantity at
the face. Using the cell value of ϕ for calculating the property at the face can potentially affect
the accuracy and stability of the numerical scheme. There are several approaches to compute the
value at the face based on cell value such as first-order upwind, second-order upwind, central
differencing, bounded central differencing, and combinations of these method. Among these
techniques, for LES simulations, CD-adapco recommends bounded central differencing and
central differencing.
Basically central differencing scheme (CDS) is a second order accurate scheme and it uses a
linear interpolation between the two nearest neighboring cell center values for the face value. It
should be noted that the central differencing scheme is very sensitive to dispersive errors
especially for positive-definite quantities such as temperature or energy and this can create
instabilities. On the other hand, the CDS does not introduce any extra decaying agents (like
turbulent viscosity) to the model and unlike the second order upwind scheme, the CDS scheme is
advantageous for discretizing quantities like velocity.
CD-adapco recommends that, for complex flows, the bounded central-differencing scheme is
used. The bounded central differencing scheme is essentially a central differencing technique but
it switches to a first order upwind scheme when the boundedness criterion is not met. (Thingbø
(2013) and CD-adapco). This feature makes the bounded central differencing scheme a more
robust approach. However, the scheme inherits the properties of upwind techniques such as fast
decay of turbulent kinetic energy.
40
2-6-4- Detached Eddy Simulation (DES)
There have been a lot of efforts on mixing LES and RANS models and take advantage of
both approaches. These models are called hybrid LES/RANS and DES is one of them. In DES,
boundary layers and irrotational flow regions are resolved by a base RANS closure, but the
model is modified so that, if the grid is refined enough, it behaves like an LES model for
detached flow regions (CD-adapco 2016).
Star-CCM+ has three options for the RANS portion of its DES model, and the user can either
choose a Spalart-Allmaras based DES, a k - ω SST, or a k – ε EB ones.
In the original Detached Eddy Simulation model, the Standard Spalart - Allmaras one
equation model was rewritten by Spalart et al. (1997) in a way that the model switches to an LES
model based on the distance function of the Standard RANS model:
DES
C d d , min
~
( 2-41)
Where C
DES
is the model constant, ∆ is the maximum grid size at the specific location, and d
is the distance from wall. d
~
is the length scale distance also known as the dissipation length
scale. This was an effort to resolve the entire boundary condition by the RANS closure while the
model can take advantage of LES features where it is required.
There have been some modifications to the original model as an effort to improve the
model’s abilities to determine where to switch from RANS to LES and vice versa. This type of
model is called Delayed Detached Eddy Simulation (DDES). The most recent efforts here are by
introducing Improved Delayed Detached Eddy Simulation (IDDES). In IDDES, the sub-grid
scale model depends on the wall distance (d) which makes the model more capable for very thin
41
near-wall regions (CD-adapco 2016). It should be noted that IDDES is the recommended model
by CD-adapco.
Another point that should be noted about the DES family is that all these models work based
on the distance from wall. This means regardless of the flow structure right at the boundary, a
RANS closure will be used. Therefore the DES model carries all characteristics of RANS models
at the boundary. Using wall distance also comes in some models with damping functions. This is
the case of Sapalart - Allmaras based DES models (CD-adapco (2016) and Spalart et al. (1997)).
The final note is that Spalart (2000) explains and compares different turbulence models and
discusses the challenges in modeling. He predicts that, for an airfoil application, LES will be
ready for everyday usage in 2045, while it takes longer for computational facilities to get to a
point that DNS will be a usable method (2080)
42
Chapter 3- Analysis of Flow past a
Cylinder by 2D Simulations
3-1- Introduction
The ultimate purpose of this research is to numerically study flow past a cylinder when it is
under a controlled and prescribed motion for Re = 4,000 and evaluate the numerical tool
capabilities against the experimental studies of Morse and Williamson (2009). In this chapter,
study of flow past a cylinder in a 2D domain is discussed.
2D simulation is chosen because, compared to 3D, 2D simulations are more economic as less
resources are required for meshing (pre-processing), solving the governing equation
(processing), and even post-processing. 3D domains have more cells and faces and the governing
equation includes more equations.
It should be noted that, as discussed before, flow past a cylinder for this range of Reynolds
number is 3 dimensional, meaning that 2D simulations may fall back regarding predicting and
reproducing the experimental results. However, since 2D modeling is much faster, it gives the
luxury of doing a large number of simulations for studying different parameters’ effects.
To find a proper configuration regarding the numerical characteristics of the simulation such
as mesh size and pattern, time step, turbulence model, and initial and boundary conditions for the
main problem, which is the vortex shedding of a controlled vibration, it is better to first simulate
the case where the cylinder is fixed (stationary cylinder).
43
Thus, in this chapter, the simulation and results for a stationary cylinder is presented first,
and then the moving cylinder case will be discussed.
In order to get to the best possible configuration of mesh and time step, the simulation results
should be checked against the experimental tests. The focus here will be on the values of drag
and lift coefficients and also Strouhal number. The drag coefficient for a cylinder is usually
expressed as the averaged value through time. For lift coefficient, however, as the mean value of
lift is zero, the root mean square (RMS) value of lift coefficient is used in the literature. For
Strouhal number, the shedding frequency is averaged for a certain number of cycles, e.g. 20
cycles. Experimental values of these quantities for Re = 4,000 are as follows:
Table 3-1 – Experimental values for drag, lift, and Strouhal number – Re = 4,000
Mean drag coefficient (Wieselsberger 1922) 0.99±0.05
RMS of lift coefficient (Norberg 2003) 0.09
Strouhal number (Norberg 2003) 0.21
Figure 3-1 – RMS values of lift coefficient in cylinders for different Reynolds numbers (Norberg 2003)
44
Figure 3-2 – Strouhal number of cylinders for different Reynolds numbers (Norberg 2003)
3-2- Two Dimensional Simulation of Flow past a Stationary
Cylinder
Star-CCM+ is basically a 3D fluid solver, but it is able to simulate 2 dimensional domain as
well. In 2D, however, there are some limitations on the features of the software and the most
important one is turbulence modeling options. Star-CCM+ does not allow the user to use LES or
DES in 2D as these models are meaningful for 3D flows.
Before discussing the results, the numerical properties of flow past a cylinder in 2D domains
are presented.
3-2-1- Flow Domain and Boundary Conditions
A rectangular domain is chosen for 2D simulation. The domain should be big enough that the
boundaries do not affect the flow. The following figure shows the flow domain for 2D
simulations.
45
Figure 3-3 – Domain size for 2D simulations
For boundary conditions, the left side is a velocity inlet, the right side is defined as a pressure
outlet. The top and bottom boundaries are usually defined as symmetry planes, meaning that the
shear stress on these boundaries is set to zero. In this research, however, these boundaries are
chosen as velocity inlets with velocity components similar to the left side boundary. There are
two reasons for this, one, these boundaries need to be modeled as velocity inlets for the moving
cylinder case which will be discussed later. Two, the domain is shown to be wide enough that the
type of boundary condition on top and bottom does not affect the result. So, for consistency
between the stationary and moving cases, the top and bottom boundaries are also defined as
velocity inlet.
The cylinder wall is modeled as a wall with no slip boundary condition.
46
Figure 3-4 – Boundary conditions for the 2D domain
3-2-2- Two Dimensional Domain Meshing
Star-CCM+ has a very robust built-in mesher. This is a very practical feature for, especially,
complex geometries. Different mesh type are available in Star-CCM+. For 3D modeling, the
following options can be used to mesh the domain (CD-adapco 2016):
Polyhedral mesh
Tetrahedral mesh
Trimmer mesh
Thin layer mesher
The area near wall (the boundary layer) can be meshed with prism layers. Prism layers
basically consist of cells with high aspect ratio. According to CD-adapco, prism layers allow the
solver to resolve the near wall region more accurately. This is important especially for flow
separation and determining the gradient of velocity normal to walls (CD-adapco 2016). Near
walls, the gradient of flow quantities is very high for turbulent flows and resolving it requires
47
more efforts and prism layers are designed for this purpose. Using prism layers also helps to
reduce number of cells in these regions. As they have a very high aspect ratio, this means that
with less total number of cells, an acceptable accuracy is obtained for resolving the near wall
flow regime.
Figure 3-5 – Using prism layers near walls to resolve turbulent flows (CD-adapco 2016)
For 2D simulation of flow past a cylinder, the domain is meshed by polygons and the area
near the cylinder is meshed by prism layers. Polygons have more sides so they are able to capture
gradients of flow properties better. However, polygonal mesher produces more cells compared to
trimmer (rectangular), at least for this case, plus that having more sides means more
computations. Thus, using polygonal mesh results in more computational time and cost.
However, since flow past a cylinder is a sensitive problem and includes flow separations,
reattachments, and high gradients, it would be good to see how the polygonal mesh configuration
performs for this case.
48
In order to reduce computational costs for meshing and solving the governing equations, it
has been tried to mesh different regions of the domain with different sizes. For this purpose, 3
ring-like zones right after the prism layer, and 3 conical zones beyond the rings are defined in a
way that, the cell size increases as getting distance from the cylinder. Obviously, the flow
regimes near the inlet and top and bottom boundaries are simple and do not need a refined mesh.
The transition between different zones should be in a way that a sudden big change in size of
cells is avoided,
Prism layers, as discussed before, need more details and attention for discretization. The
parameters that are considered in meshing the near wall region include:
Number of circumferential elements,
Total thickness of prism (boundary) layer,
Number of prism layers,
Prism layers stretching factor,
1
st
prism layer thickness and aspect ratio
Last prism layer aspect ratio
For Re = 4,000 the boundary layer is still laminar, so there should be no need to have a very
thick prism layer or a lot of prism layers to resolve the boundary layer. The thickness of the 1
st
prism layer is important as it defines the value of y
+
.
y
+
is a dimensionless parameter for distance from the wall. Its definition is as follows
(Wilcox 2006)
u y
y
( 3-1)
49
where, y is the distance from the wall,
u is the friction velocity, and υ is the fluid viscosity.
w
u ( 3-2)
w
is the wall shear stress.
If the velocity is normalized by
u and draw against y
+
, the velocity profile for a turbulent
boundary layer will be obtained
u
u
u
( 3-3)
Figure 3-6 – Plot of u+ vs y+ for a typical boundary layer (Wilcox 2006)
The region where the relationship between u
+
and y
+
is logarithmic is called the log layer and
for the y
+
values smaller than the log layer (0< y
+
< 30), the flow is defined as the viscous
sublayer (Wilcox 2006).
Since for Re = 4,000 the boundary layer is laminar, it is better to mesh the near wall region in
a way that y
+
values will be as small as possible. The thickness of the 1
st
prism layer defines the
y
+
value and in this study has been determined in a way that y
+
will always be less that unity.
50
Although using prism layers allows the user to have elements with high aspect ratios, the
aspect ratio should not be too big in order to resolve the boundary layer more accurately.
Because, the prism layer thickness value is stretched, the 1
st
prism layer has the highest aspect
ratio. For having a better resolution of boundary layer, the aspect ratio of the 1
st
prism layer is
bounded to 50.
The last prism layer aspect ratio is also important in order to have a good transition to
polygonal mesh region. It has been tried to have an aspect ratio of about 2 for the last prism
layer.
Figure 3-7 – A big jump in transition from prism layers to polygonal mesh (poor transition)
Figure 3-8 – An example of a good transition from prism layers to polygonal mesh
The size of the 1
st
layer of polygonal mesh right after the boundary layer follows the prism
layer proportions and the mesh size gradually increases (or decreases) to the specified value for
the 1
st
ring. The mesh size as discussed above is controlled by rings and cones to have a smooth
51
and at the same time fined mesh for the regions near the cylinder. The maximum mesh size of
the entire domain is set to be less than 1.0D.
The final configuration of mesh is shown below
Figure 3-9 – Mesh configuration for 2D simulation of flow past a stationary cylinder
It should be noted that, this mesh will be used for the moving cylinder case as well. So, the
refinements by rings and cones should be in a way that they will be able to resolve the flow
structures behind the moving cylinder properly.
52
Figure 3-10 – Details of mesh for flow past a stationary cylinder
All the simulations are run with an implicit 2nd order temporal discretization. Two
parameters are chosen to control the simulations: maximum value of wall y+ and maximum local
convective CFL number. CFL or the Courant number is the ratio of the physical time-step to the
mesh convection time scale (CD-adapco 2016)
U
dx
t
CFL
( 3-4)
where, Δt is the time step, U the fluid velocity relative to the mesh, and dx is the distance that
fluid moves in one time step.
For this research study, the following conditions are tried to be held for simulations:
1) Maximum Convective CFL is less or around unity for the entire domain,
2) For the cylinder wall, Max (y
+
) < 1.0
53
As discussed before having y
+
< 1.0 helps to resolve the boundary layer more properly. It
should also be noted that, although all simulations are implicitly run, limiting CFL to about 1.0
will lead to a better solution. CFL = 1.0 means that the fluid moves by about one cell per time
step.
The time step is kept the same throughout a simulation, this means setting the time step in a
way that CFL < 1.0 is not possible before hand and needs a trial and error procedure. For this, a
field function is defined to find the minimum cell size in the domain.
Cell
Area Surface Size Cell _ ( 3-5)
To estimate the time step properly, the far field velocity with some safety factor is used.
Cell
U
Size Cell CFL
t
min
max
( 3-6)
This equation shows that, if the minimum cell size decreases (a more refined mesh), the
maximum time step for keeping CFL < 1.0 will be lower and this implies that for a more refined
mesh, a smaller time step is required. Therefore, for a finer mesh, more computational resources
and time will be needed.
For the turbulence modeling, as discussed before, the only option for 2D domain is using
RANS models. At the beginning, according to CD-adapco’s suggestion, the two-layer Realizable
k – ε model was selected. Star-CCM+ have three wall treatment options for most of the
turbulence models; high y
+
, low y
+
, and all y
+
or two layer. According to CD-adapco (2016), the
high y
+
wall treatment is basically the classic wall-function approach, in which all quantities are
calculated based on equilibrium turbulent boundary layer theory. Low y
+
wall treatments are
fundamentally the low-Reynolds number modifications for turbulence models, where no wall
54
laws are required because the viscous sublayer is resolved. The all y
+
approach mixes the high
and low y
+
treatments and depending on the mesh resolution chooses one of them (CD-adapco
2016).
Different turbulence models for some of the mesh configurations have been used to
investigate how these models affect the result. It is worth noting that different turbulence models
need different mesh considerations; for example k – ω models are well-known to be more
sensitive to mesh size, and they need more refinements. However, having y
+
values less than one
supposedly provides a good refinement at least normal to the wall for all turbulence models.
To reduce the computational time and to be able to run more simulations, all simulations are
initiated by a steady run. Here, this approach is explained by an example for one of the cases:
A domain with the computational quantities presented in the following table is
Table 3-2 – Properties of model for a stationary cylinder 2D simulation
Total number of cells 15,778
Time step (sec) 0.004
ΔtU/D 0.011
Number of circumferential elements 126
Number of prism layers 6
Boundary layer thickness 0.025D
Cell size right after the boundary layer 0.03D
Turbulence model Realizable k – ε all y
+
- Initiating the simulation from rest
The simulation should normally get initiated from rest and the velocity field increases to the
ambient velocity gradually. So, in this example, the velocity at the inlet boundaries is increased
from [0,0,0] to [U,0,0] in 10 seconds. Considering that the Strouhal number is 0.21, 10 seconds
is about 5.81 times the natural shedding frequency, so should be large enough. For each time
55
step, 25 iterations are run to decrease the residuals as much as possible. The forces tend to a
quasi-static response eventually and they look like sinusoidal signals with a constant amplitude
and frequency.
- Initiating the simulation by a steady run
If the simulation starts with the boundary velocity of [U,0,0], the cylinder drag and lift
experience very large and unphysical values for the first couple of iterations. They damp to
smaller values in time and eventually a quasi-static response is obtained for both. To mitigate
these huge jumps in drag and lift, a steady run will be helpful. Therefore, here, the flow is run
with the boundary velocity of [U,0,0] under a steady simulation for 1,000 iterations. So, when
the unsteady implicit simulation starts, the flow field has already somehow been formed behind
the cylinder.
To reduce the effects of this unphysical steady run and decrease the values of residual for
simulation parameters, the implicit simulation is first run for 1.0 second with 100 iterations for
each time step. After one second, the iteration number is changed to 25. In this way, the quasi-
static result is achieved in a much smaller time (in more than 10 seconds less) which is quite
significant considering the size of domain and the time step. The final response for drag and lift
is almost identical to the previous case.
56
Figure 3-11 – Lift coefficient for cases with or without steady run initiation
The comparison shows the maximum lift amplitude is very close for both cases. The same
comparison for the drag and Strouhal number is valid as well. It is known that turbulent flows are
very sensitive to initial and boundary conditions and clearly these two cases have different initial
conditions, but the final result is almost the same. The reason could lie in using RANS models
for turbulence; the nature of averaging in RANS mitigate the dependency on initial conditions
and this is one the features of RANS models (Wilcox 2006).
The last important subject for numerical simulation of flow past a cylinder is the number of
iterations per time step and how or when to stop iterating for each time step. In this research
study, calculations proceed to the next time step after certain number of iterations. Number of
iterations depends on time step, coarseness of the mesh, and turbulence model, it is set between
15 and 50.
57
3-2-3- Results for 2D Simulation of a Stationary Cylinder
The model is run with different mesh configurations and turbulence models. Then the values
of lift, drag, and Strouhal number are evaluated for each case after reaching to a quasi-steady
solution. The response for lift looks like a sinusoidal and drag eventually fluctuates between two
limits. For the before-mentioned case, the results are shown below
Figure 3-12 – A typical plot of drag coefficient in time for 2D simulation of a stationary cylinder
Figure 3-13 – A typical plot of lift coefficient in time for 2D simulation of a stationary cylinder
For CD
mean
, the average of the maximum and minimum of drag is calculated. To find the
RMS of lift, since the response is almost periodic, the RMS value will be the maximum over √ 2
.
58
The following table summarizes the results for several cases with their mesh configurations,
time step, and turbulence model.
59
Table 3-3 – Results for 2D simulation of the stationary cylinder with RANS
Sim
Name
# Cells ΔtU/D
Turbulence
Model
Max
y+
# Inner
Iteration
Max
CFL
CL,
rms
CD,
mean
St #
# Cells
on
Cylinder
BL
thickness
# Prism
Layers
Cell Size
right
after the
B.L.
Cell size
in the
wake
Cell size
in far
wake
12_00_
00
508,555 2.77E-05
RKE-Two
Layer
0.21 20 0.007 - - - 158 0.035D 10 0.02D 0.025D 0.025D
13_00_
00
208,786 2.77E-04
RKE-Two
Layer
0.065 20 0.12 0.329 0.95 0.235 450 0.0112D 10 0.0075D 0.025D 0.2D
15_00_
00
213,532 2.77E-04
RKE-Two
Layer
0.067 25 0.19 0.346 0.96 0.235 628 0.0112D 10 0.005D 0.025D 0.3D
16_00_
00
163,094 2.77E-04
RKE-Two
Layer
0.067 20 0.193 0.354 0.96 0.234 628 0.0112D 10 0.005D 0.03D 0.3D
17_00_
00
177,540 2.77E-04
RKE-Two
Layer
0.42 25 0.056 0.368 0.985 0.231 210 0.075D 12 0.015D 0.025D 0.3D
18_00_
00
106,367 2.77E-04
RKE-Two
Layer
0.41 25 0.06 0.318 0.948 0.234 218 0.02D 6 0.015D 0.025D 0.3D
19_00_
00
98,153 2.07E-03
RKE-Two
Layer
0.4 25 0.91 0.318 0.948 0.234 419 0.02D 6 0.0075D 0.03D 0.3D
19_00_
01
98,153 2.07E-03
SST KW
Low y+
0.454 25 1.09 1.379 1.75 0.252 419 0.02D 6 0.0075D 0.03D 0.3D
19_00_
02
98,153 2.07E-03
St Spallart
Alimaras -
Low y+
0.307 25 1.05 1.138 1.553 0.245 419 0.02D 6 0.0075D 0.03D 0.3D
20_00_
00
102,447 1.66E-03
RKE-Two
Layer
0.241 25 0.986 0.322 0.95 0.235 483 0.02D 8 0.0065D 0.03D 0.3D
60
Sim
Name
# Cells ΔtU/D
Turbulence
Model
Max
y+
# Inner
Iteration
Max
CFL
CL,
rms
CD,
mean
St #
# Cells
on
Cylinder
BL
thickness
# Prism
Layers
Cell Size
right
after the
B.L.
Cell size
in the
wake
Cell size
in far
wake
21_00_
00
73,161 2.77E-03
RKE-Two
Layer
0.27 25 0.68 0.339 0.966 0.233 232 0.04D 9 0.014D 0.035D 0.3D
21_00_
01
73,161 2.77E-03
SST KW
low y+
0.315 25 0.79 1.379 2.1 0.253 232 0.04D 9 0.014D 0.035D 0.3D
21_00_
02
73,161 2.77E-03
St Spallart
Alimaras -
Low y+
0.21 25 0.77 1.138 1.585 0.246 232 0.04D 9 0.014D 0.035D 0.3D
22_00_
00
125,711 1.38E-03
RKE - Two
Layer
0.38 20 0.9 0.325 0.952 0.235 583 0.005D 6 0.015D 0.03D 0.3D
23_00_
00
63,482 3.60E-03
RKE - Two
Layer
0.495 30 0.78 0.322 0.952 0.233 203 0.025D 6 0.0155D 0.035D 0.3D
23_00_
01
63,482 3.60E-03
SKE - Two
Layer
0.495 30 0.72 0.233 0.92 0.201 203 0.025D 6 0.0155D 0.035D 0.3D
23_00_
02
63,482 3.60E-03
V2F KE -
low y+
0.54 30 0.87 1.089 1.45 0.244 203 0.025D 6 0.0155D 0.035D 0.3D
23_00_
03
63,482 3.60E-03
EB KE -
low y+
0.54 30 0.87 1.138 1.51 0.248 203 0.025D 6 0.0155D 0.035D 0.3D
24_00_
00
54,519 3.32E-03
RKE Two
layer
0.574 40 0.96 0.327 0.95 0.233 251 0.02D 6 0.0125D 0.04D 0.4D
24_00_
01
54,519 3.32E-03
SKE Two
layer
0.572 40 0.9 0.226 0.913 0.202 251 0.02D 6 0.0125D 0.04D 0.4D
Sim
# Cells ΔtU/D
Turbulence Max # Inner Max CL, CD,
St #
# Cells
on
BL # Prism
Cell Size
right
Cell size
in the
Cell size
in far
61
Name Model y+ Iteration CFL rms mean Cylinder thickness Layers after the
B.L.
wake wake
24_00_
02
54,519 3.32E-03
EB KE All
y+
0.623 40 1.18 1.124 1.515 0.251 251 0.02D 6 0.0125D 0.04D 0.4D
24_00_
03
54,519 3.32E-03
V2F KE
All y+
0.622 40 1.15 1.068 1.45 0.244 251 0.02D 6 0.0125D 0.04D 0.4D
25_00_
00
129,222 2.77E-03
RKE Tow
layer
0.382 25 0.936 0.318 0.94 0.235 309 0.015D 6 0.01D 0.03D 0.2D
25_00_
01
129,222 2.77E-03
SKE Tow
layer
0.381 20 0.931 0.227 0.908 0.203 309 0.015D 6 0.01D 0.03D 0.2D
25_00_
02
129,222 2.77E-03 LAMINAR - 25 1.815 1.393 1.65 0.262 309 0.015D 6 0.01D 0.03D 0.2D
26_00_
00
84,694 3.87E-03
SKE Two
Layer
0.496 25 0.732 0.235 0.918 0.203 209 0.025D 6 0.015D 0.035D 0.2D
27_K
W_01_
01
45,632 3.87E-03
SKW All
y+
0.644 25 0.697 1.492 1.65 0.22 157 0.025D 6 0.02D 0.04D 0.4D
27_K
W_02
45,632 4.15E-03
Rey Stress
Two Layer
0.519 30 0.672 0.655 1.212 0.251 157 0.025D 6 0.02D 0.04D 0.4D
62
3-2-4- Discussion on Results for the Stationary Cylinder
Comparing the result of 2D simulation with the experimental data shows that the values of
lift coefficient are higher in simulation. RANS models, especially the k – ε family, do predict the
value of mean drag and the Strouhal number with an acceptable accuracy. For example, the
Realizable k – ε model, in some cases, has the value of mean drag in the range of experimental
data (0.95 ~ 1.00), and the Strouhal number predicted by the Standard k – ε is very close to the
experiment (with only 3 % error). The best prediction for the lift coefficient is coming from the
Standard k – ε with about 150% error. In general, the values of drag and lift resulted from any
RANS model except the Realizable and Standard k – ε are much higher and these models
perform very poorly in this regard.
Comparing the result with other researchers’ numerical efforts is also notable. Thingbø
(2013) conducted 3D LES in Star-CCM+ for Re = 3,900 and the statistical values of lift, drag,
and Strouhal number are in a very remarkable agreement with the experiment. In her study,
Thingbø (2013) used about 4 million structured cells in a 3 dimensional domain with thickness
of 4D and very small time step ( Δt·U/D = 5.13×10
-5
). Comparing the current 2D RANS results
with Thingbø’s results also reveals another characteristics of 2D RANS simulation:
Figure 3-14 – Comparison between 2D RANS and 3D LES; on left, 2D RKE with 15,778 cells, on right,
3D LES with 3.78 million cells (Thingbø 2013)
63
The noticeable difference between the result of 2D RANS and 3D LES, beside the high
values of lift, is the shape of response. In 3D LES, the coefficients fluctuate over time with some
significant amplitude modulations, while in 2D RANS, there is almost no or insignificant
amplitude modulations in the response. The lift response looks more like a pure sinusoidal signal
and drag has some very small amplitude modulations as shown before.
These amplitude modulations are observed in experiments as well.
Figure 3-15 – Oscillating lift and drag for cylinders at a subcritical Reynolds number (Humphreys 1960)
So, the discrepancies of 2D simulation using RANS models are not coming from mesh
refinements or size of time increment. As shown in Table 3-3, using more a refined mesh or a
smaller time step, after some point does not help with getting a better result. The differences
between the 2D RANS simulation results and the 3D LES are possibly coming from using
RANS models and simulating in a 2D domain.
Vortex shedding is a three dimensional phenomenon for Re = 4,000 and 2D planar shear
layers are broken in the span-wise direction. Formation of 3D structures and disruption of strong
2D shear layers shed from the cylinder reduce the lift value. This can be observed in Norberg's
(2003) experiments too (Figure 3-1); the value of lift drops in the range of 300 < Re < 1,000
where it is shown by Williamson (1988) that laminar vortex shedding transforms to a fully 3D
64
flow field. It is worth noting that even using DNS in a 2D domain can lead to a similar result as
shown by Singh and Mittal (2005). They used a finite element technique for solving flow past a
cylinder in a 2D domain with no turbulence model for a large range of Reynolds number with a
focus on drag crisis. Their results of lift and drag time history for Re = 3,900 is as following:
Figure 3-16 – Lift (…) and drag (―) coefficients time history for Re = 3,900 using 2D DNS (Singh and
Mittal 2005)
Using RANS models even in 3D domains results in higher values for lift and drag. Young
and Ooi (2010) compared RANS and LES for Re = 3,900 where they used the same planar mesh
configuration with different discretization in the thickness of the domain. They showed that the
Standard k – ω predicts higher values for lift and drag. They also showed that RANS models in
3D can simulate the amplitude modulations in the response.
RANS models are notorious for generating high and artificial (turbulent) viscosity to the
system (Salvetti et al., 2010) which could reach up to 1,000 times the fluid viscosity. RANS
models, especially one- and two-equation ones, are also known to fail against flows with high
separations and high flow stream curvatures (Wilcox 2006). One- and two-equation models were
originally devised for 2D fully turbulent flows and one should expect some level of inaccuracy
65
when modeling a highly turbulent flow with transitions, separations, and reattachments. On the
other hand, they are easy to use, well suited for coding, and computationally inexpensive.
Comparing the current study result with other numerical efforts show that predicting the
statistical values of forces and the Strouhal number numerically is not an easy work and needs a
lot of mesh and time step refinements-adjustments and use of complicated turbulence models.
This is especially valid for the lift coefficient and there are a lot of studies in the literature that
the authors did not report their model’s lift prediction. The following table shows how the
current study performs compared with other researchers’:
Table 3-4 – Results for flow past a cylinder at Re ≈ 4,000 by different CFD studies
Author Method CD
mean
CL
rms
St #
Current Study 2D RANS 0.99 0.226 0.203
Thingbø (2013) 3D LES 1.01 0.09 0.21
Breuer (1998) 3D LES 1.08 - -
Lysenko et al (2012) 3D LES 0.97 0.09 0.209
Wornom et al (2011) 3D VMS LES 0.99 0.108 0.21
Ouvrard et al (2008) 3D VMS hybrid RANS/LES 1.00 - 0.212
Prsic & Myrhaug (2012) 3D LES 1.1 0.228 0.205
Wu (2011) 2D RANS 1.518 1.058 0.211
Tuan & Temarel (2014) 2D RANS (RKE) 0.98 0.27 -
Experiment - 0.99±0.05 0.09 0.21
In the next chapter, some efforts on 3D simulation of flow past a stationary cylinder will be
presented. RANS models will be utilized for 3D simulations to see if going to 3D from 2D
exhibits any changes and improvements in results. However, the 2D simulation of moving
cylinder case will be discussed in the coming section.
66
3-3- Two Dimensional Simulation of Flow past a Moving
Cylinder
In this section, the flow past a moving cylinder is discussed. Morse and Williamson (2009)
did an extensive experimental study on vortex shedding behind a cylinder when it is under a
controlled prescribed oscillation. They observed a new vortex shedding pattern which borrows
some characteristics from well-known 2P shedding pattern, but one of the vortices is weaker and
dissipates much faster through the wake. They called this new pattern (mode) 2P
O
, which has
overlap with 2P and 2S modes, and showed that it is a stable mode and associated with the
highest possible amplitude in freely vortex induced vibration for Re = 4,000.
They put the result for one case in their publication that lies in the 2P - 2P
O
region.
Figure 3-17 – Force and energy-transfer time traces for A
*
= 0.8 and λ
*
= 5.6 (2P - 2P
O
region) from
Morse and Williamson (2009)
67
The RMS value of the lift coefficient is 0.678 according to Figure 3-17. In this figure, y/D is
the imposed sinusoidal oscillation to the cylinder, C
Y
is the lift coefficient, φ is the phase lag
between C
Y
and y/D and É
*
is the instantaneous rate of energy transfer. The thick solid line in the
energy diagram shows the cycle averages. According to this diagram, 2P
O
mode yields a net
positive fluid excitation while the 2P mode has a net negative energy. The corresponding vortex
shedding pattern to this state (A
*
= 0.8 and λ
*
= 5.6) is shown below:
Figure 3-18 – 2P
O
shedding mode for A
*
= 0.8 and λ
*
= 5.6 from Morse and Williamson (2009)
The goal in this section is to see how 2D simulation with RANS models performs when the
cylinder is under a controlled motion.
3-3-1- Parameters of the Model
The same domain with the same type of turbulence models will be used. The mesh
configurations with the best results for the stationary case are selected. Since the cylinder will be
under motion, the time step chosen should be lower than the stationary case in order to keep the
68
Courant number less than unity. Maximum y
+
values for most cases were under 0.5 in the
stationary cylinder analysis, so the value of y
+
should most likely be less than one for the moving
case too.
3-3-2- Apply the Motion to the Cylinder
The only new issue for the moving case is how to simulate the cylinder motion. There are
three options in Star-CCM+ for applying motions to the model:
1) Mesh morphing: In mesh morphing the mesh needs to be regenerated for each time step.
So, for example, the sinusoidal motion is applied to the no slip boundary of the cylinder
and the code regenerates the entire mesh at the beginning of each time step. The problem
with this method is that the mesh quality can be affected throughout the simulation. It has
been already shown that flow past a cylinder is very sensitive to mesh and distorting the
mesh will definitely result in a poor prediction especially in the wall region. The other
issue is the cost of meshing for each time step which can be significant especially for
unstructured meshes considering the cell size refinements near the cylinder.
2) Mesh overset: One independent mesh is generated and moves on top of a fixed mesh and
the fluid solver needs to couple the solution of both meshes in the overlapped zones.
There is no need to regenerate a new mesh for each time step and the moving mesh
moves as a rigid body. This method also can give a poor mesh quality for flow past a
cylinder and is not recommended.
3) Moving Reference Frame (MRF): Instead of moving the boundaries, the associated
velocity field can be applied and added to the domain. In this way, the quality of mesh is
preserved since it looks like the original mesh moves like a rigid body. It is also
69
computationally cheaper because the mesh moves incrementally, and there is no need for
the mesh regeneration or connecting two different meshes. This method is suitable for
flow past a cylinder because of the simple geometry and is used in this research. The only
concern here is that the domain should be wide enough so that the flow properties will be
independent of the side boundaries as long as these boundaries are modeled as velocity
inputs.
More details of mesh motion techniques will be presented in Chapter 6.
It is required to define an imposed velocity field to the domain. So if the cylinder motion is a
sinusoidal function, then velocity field for MRF will be
t f A f
D
y
t f A
D
y
2 cos 2 2 sin
* *
( 3-7)
In Star-CCM+, the motion can be applied to the domain as a field function. So, a vector of
velocity is defined and the reference motion of the domain is changed from static to the defined
motion.
Figure 3-19 – 2D domain boundary conditions and velocity field for simulating the moving cylinder
70
The mesh configuration obtained from the stationary cylinder simulation is used for
analyzing the effects of cylinder oscillation on shedding pattern and lift coefficient. At the
beginning the Standard and Realizable k – ε models are selected since they had the best
performance for the stationary case. But, for one case, other turbulence models are also tested.
3-3-3- Results for 2D Simulation of the Moving Cylinder
Results show that, surprisingly, the Standard and Realizable versions of k – ε lead to, again, a
response with no or insignificant amplitude modulations and the RSM value of the lift coefficient
is less than what Morse and Williamson’s data shows.
The values of maximum y
+
, maximum CFL number, and the RMS of lift coefficient for
different simulations are presented in the following table
Table 3-5 – Results for flow past a moving cylinder in a 2D domain with RANS models
Sim Name
Time Step
(s)
Δt.U/D Turbulence Model Max y+ Max CFL CL,rms
_12_00_00 0.001 2.77E-03 RKE all y+ 0.282 0.808 0.287
_26_00_00 0.00135 3.73E-03 SKE all y+ 0.68 1.46 0.215
_27_AKNKE 0.0015 4.15E-03 AKNKE Low y+ 0.8 1.03 1.6
_27_KW_00 0.001 2.77E-03 KW SST All y+ 0.82 0.78 1.73
_27_KW_01 0.00135 3.73E-03 SKW All y+ 0.923 1.219 1.723
_27_KWSST 0.00128 3.54E-03
KWSST Low y+ with
Curvature Correction
0.8 0.95 1.494
_27_RKE 0.0013 3.60E-03
RKE All y+ w/ Curvature
Correction
0.67 0.84 0.651
_27_RSM 0.0015 4.15E-03 RSM All y+ Linear 0.74 1.01 0.921
It can be seen that using other models results in a higher lift coefficient. The statistical value
of lift for the case where the Realizable k – ε is used with the curvature correction (_27_RKE) is
71
very close to the experiment, but as discussed before, the shape of response is not what the
experiment shows.
Figure 3-20 – Lift coefficient for Moving_12_00_00 (left) and Moving_26_00_00 (right)
Figure 3-20 shows the lift for cases with the Realizable (_12_00_00) and Standard
(_26_00_00) k – ε models. It is clear that the lift has no amplitude modulations. However, vortex
pattern for these two cases is very different; The Realizable version leads to a mode that looks
like 2P, but the shedding mode resulted by the Standard k – ε looks more like a stretched 2S
mode.
Figure 3-21 – Vortex shedding pattern for Moving_12_00_00 (left) and Moving_26_00_00 (right)
Vortices dissipate much faster when the Standard k – ε is used while the 2P mode in
_12_00_00 case is preserved for a long distance. The vorticity generated by the Stnadard k – ε is
smaller according to this figure. The mesh is much more refined in the near and far wake for
_12_00_00 and this can be a reason but the fast dissipation of vortices is most probably because
of the turbulence model.
72
Figure 3-22 – Lift coefficient and vortex pattern for Moving_27_RKE case (RKE with Curvature
Corrections)
The _27_RKE case has a mesh configuration very similar to _26_00_00, but it can be
observed that the 2P mode is more stable than the stretched 2S in _26_00_00.
Using other turbulence models rather than the Realizable and Standard k – ε will produce
some level of amplitude modulations in lift. However, the RMS value of the lift coefficient is
higher than the experiment and the shedding mode does not look like 2P, 2P
O
, or 2S in most
cases:
Figure 3-23 – Lift coefficient and vortex pattern for Moving_27_AKNKE case (AKNKE Low y
+
)
Figure 3-24 – Lift coefficient and vortex pattern for Moving_27_KW_00 case (KWSST All y
+
)
73
Figure 3-25 – Lift coefficient and vortex pattern for Moving_27_KW_01 case (SKW All y
+
)
Figure 3-26 – Lift coefficient and vortex pattern for Moving_27_KWSST case (KWSST Low y
+
with
Curvature Correction)
Figure 3-27 – Lift coefficient and vortex pattern for Moving_27_RSM case (RSM All y
+
with Linear
Pressure Strain Correlation)
The figure of vortex shedding for all cases shown above is the vorticity field between -15 and
15 (1/s). Comparing the vorticity field for all cases shows that the magnitude generated by the
Standard and Realizable k – ε is lower, so is the lift coefficient. However using these two models
results in a symmetric lift and a distinguishable shedding pattern. The difference between the
Standard and Realizable k – ε models with other RANS models is more pronounced for 3D
simulation of a stationary cylinder and will be discussed in more details in the next chapter.
74
Thus, 2D simulation of flow past a moving cylinder with RANS does not produce a good
result and agreement with the reported experimental data is very poor. In case of the Standard
and Realizable k – ε models, the lift coefficient shows no amplitude modulations and it looks like
that the phase lag between the lift and the cylinder oscillation is constant through time. The other
turbulence models produce some level of amplitude modulations but in most cases the response
is not symmetric and the vortex shedding has no pattern. 2D simulation seems like to not be able
to produce the 2P
O
mode either.
Flow past a moving cylinder has even more complex structures compared to the stationary
cylinder. So, one should not expect that 2D simulation with RANS turbulence models performs
well for the moving case. This proves the necessity of using 3D models with more complicated
turbulence models to capture more details behind the cylinder.
In the next chapter, the flow past a stationary cylinder in a 3D domain will be studied to see
how using a 3D domain can affect the result. The RANS models will be still used and this helps
to compare the result of the same 2D planar mesh with the corresponding 3D mesh. Different
mesh configurations, domain thicknesses, and turbulence models will be used to study their
effects in more details.
3-4- Concluding Remarks
3-4-1- Two Dimensional Simulation of Flow past a Stationary Cylinder
Simulating the stationary cylinder problem is useful for achieving a proper mesh
configuration and time discretization, and figuring out which turbulence model performs better
for flow past a cylinder. These properties will later be implemented for analyzing the moving
cylinder case.
75
Comparing the results of 2D simulations with the experimental data shows that some of the
simulations are able to predict the mean drag and the Strouhal number with an acceptable
accuracy. Refining the near-wall mesh after some point does not produce a better result and
results seem to be independent of the mesh size (at least for the Standard and Realizable k – ε
models). However, there are two significant differences between the experimental data and the
2D simulation results: 1) Higher values of lift coefficient, and 2) The time histories do not show
any amplitude modulations. The value of lift and drag for most of turbulence models are much
higher than the experiments.
These two discrepancies are coming from the nature of 2D domains and RANS models when
dealing with a 3D phenomenon with laminar-turbulent transitions, flow separations, and
reattachments. Strong 2D shear layers produced by vortex shedding, in reality, are broken in the
3
rd
dimension which is not feasible in 2D domains. On the other hand, RANS models introduce
unnecessary viscosity to the system which affects the fluid excitation and increases the fluid
forces.
3-4-2- Two Dimensional Simulation of Flow past a Moving Cylinder
Using the mesh configuration from the stationary case, simulating the flow when the cylinder
is under a sinusoidal motion requires a smaller time step regarding accuracy and in general
simulations need more computational resources. The statistical values of lift resulted by 2D
simulation by the Standard and Realizable k – ε are lower than the experiment while the other
RANS models have a higher lift. However, these two k – ε models lead to a lift time history with
no significant amplitude modulations.
The shedding pattern of the 2D simulation does not look like what was observed in the
experiment. Among all turbulence models, the Realizable k – ε develops a 2P shedding mode
76
which is the nearest result to the experiment and in most of cases RANS models do not produce a
distinguishable shedding pattern. Using the Realizable k – ε with the curvature corrections has
the lift value very close to the experiment, however, the time history has no amplitude
modulation and the phase lag between the motion and lift seems be to constant through time.
In general, simulating the flow past a moving cylinder with 2D simulations and RANS
models has less encouraging results compared to the stationary case because shortcomings of 2D
simulations and RANS are more noticeable in this case. 3D simulations with more complicated
turbulence models like DES could help with getting a better result and will be investigated later
for this study.
77
Chapter 4- Analysis of Flow past a
Stationary Cylinder by 3D Simulations
and RANS Models
4-1- Introduction
Williamson (1988) explains how the flow past a cylinder transforms from a 2D to a 3D
structure. For 170 < Re < 180, the first stage of transformation occurs, which he calls it “Mode
A”. In this range, stream-wise vortex loops are established and by further increase of the
Reynolds number, finer-scale stream-wise vortex pairs start to form in the range of 230 < Re <
260 (“Mode B”)
Figure 4-1 – Mode A (left) and B (right) of transformation to 3D vortex shedding (Williamson, 1996)
78
However, for Re > 1,000 the flow regime is called the shear-layer transition (Williamson,
1996), where secondary (or Kelvin-Helmholtz) vortices in the shear layer begin to form and the
flow is completely three dimensional.
Figure 4-2 – Stages of transition from the primary von Karman vortex shedding to a 3D flow
(Williamson, 1996)
79
Numerical simulation of flow past a cylinder in a 2D domain with RANS models shows that
the result has two significant differences compared to the experimental data; 1) high values of
force coefficients especially the lift, and 2) no or little amplitude modulations in time history.
One of the reasons is resolving the 3D turbulent wake in a 2D domain, and this chapter focuses
on the effects of 3D modeling on the result. Since, RANS models are the only option for 2D
simulations in Star-CCM+, the 3D simulations are also run by RANS models so in this way, it
will be possible to make comparison between 2D and 3D simulations with focus on the domain
dimension.
In this chapter, first the model parameters are addressed first and then the results will be
discussed.
4-2- Model Parameters
For meshing, first a planar (x-y plane) mesh configuration is constructed with polygons and
prism layers first, and it is then extruded in z direction to create a 3D domain. For comparing 2D
and 3D simulations it is essential to have the same planar mesh.
4-2-1- Planar Mesh Configuration
In order to be able to conduct many simulations to study different parameters effects, a
relatively coarse planar mesh in the far wake is used. The near wake and the zone close to the
cylinder is still meshed with more care and refined enough to have an acceptable result.
The same boundary conditions and mesh configuration (ring-like near the cylinder and
conical in the wake) are used. The planar mesh parameters are shown below
80
Table 4-1 – Properties of the planar mesh for 3D simulation with RANS
Total number of cells 15,778
Number of circumferential elements 126
Number of prism layers 6
Boundary layer thickness 0.025D
Cell size right after the prism layer 0.03D
Cell size in the near wake 0.1D
Cell size in the far wake 0.4D
The planar mesh is shown in the following figure
Figure 4-3 – Planar mesh configuration used for 3D simulation of flow past a cylinder
Simulations show that this mesh configuration results in CL
RMS
= 0.247 with the Standard k –
ε model, which is comparable with the best result derived from more refined meshes shown in
the previous chapter.
81
4-2-2- Mesh Discretization and Domain Size in the z Direction
When doing 3D simulation of flow past a cylinder, two important parameters are the
thickness of the domain (size in z direction) and the discretization in z meaning that how and
how many elements should be used in z direction.
For the thickness, two sizes are used, 1D and 3D. After reviewing the literature, it was
noticed that the most used size for the z direction is πD but there has been research studies which
show that this is not an enough depth for capturing the whole phenomenon and the depth should
be much bigger. However, for example, Thingbø (2013) worked on a domain with a depth of 4D
and was able to get a very good result with LES for Re = 3,900. Since this chapter focuses on the
effects of migrating from 2D to 3D domain and getting a good result is not important here, using
the depth of 3D should be enough to observe different 3 dimensional phenomena.
Solving a 3D fluid mechanics needs more computational resources; there are more cells and
faces in 3D and obviously there are more equations to solve and more boundary conditions to
satisfy. To reduce the computational costs of meshing and solving the governing equations,
instead of using polyhedral elements, prismatic elements will be used. Generating a polyhedral
mesh will create so many cells because of refinements near the cylinder. Having very small cells
needs a smaller time step if the goal is still limiting the Courant number to 1, which even adds
more computational costs. Moreover, the planar polygonal mesh is not necessarily the exact
projection of the polyhedral mesh on the x-y plane.
So, the planar mesh in extruded in z direction. Star-CCM+ has a meshing technique for
extruding the planar mesh called thin mesher. In thin mesher, a thickness is set so that the
software realizes the domain is thin and then discretizes the domain according to the desired
number of elements in z direction. This helps a lot to reduce the cell count, but most of the
82
domain cells will have poor dimensional proportions because near the cylinder the cells look like
too stretched in z direction while near the boundaries cells are too thin.
Refinement in z direction is also important. In this chapter, three values per 1D are used: 8,
12, and 16 discretization per 1D of thickness. This helps to compare the effects on z refinements
on the results as well. Other researchers have used different values for z discretization; for
example Thingbø (2013) used 10.5 elements per 1D. Prsic and Myrhaug (2012) conducted LES
for Re = 3,900 and 13,100 with different domain thicknesses and z refinements where they used
25 and 37.5 elements per 1D and concluded that the planar mesh refinement is much more
important and the refinement along the cylinder does not affect the result significantly.
Figure 4-4 – Domains with thickness of 3D (36 layers) on left and with thickness of 1D (8 layers) on right
83
4-2-3- Turbulence Model
Sticking to RANS models, almost all the options for turbulence available in Star-CCM+ has
been used and some interesting features of these models are found comparing the 2D and 3D
simulation results.
4-2-4- Boundary Conditions on Top and Bottom
The type of boundary condition on the top and bottom of domain can have significant effects.
Pereira et al. (2015) used cyclic and symmetric boundaries with RANS, LES and hybrid
LES/RANS models for Re = 3,900 and according to their simulations using symmetry planes for
the top and bottom results in a better result. In this chapter, both types are used, but since the
focus is on how 3D simulations perform against 2D ones, no conclusion has been drawn
regarding the type of boundary.
4-2-5- Controlling Parameters and Running the Simulations
The same controls implemented for 2D simulations are used for 3D modeling as well; 1)
keeping the Courant number less than unity, and 2) limit the maximum value of y
+
to 1.0. All
simulations are run with an initial steady state and then a 2
nd
order implicit unsteady with 100
iterations per time step for the 1
st
second and 50 iterations per time step for the rest of time. To
choose the proper time step, a similar procedure to what was done for 2D simulation is followed
by finding the smallest cell size for a better estimation.
84
4-3- Results of 3D Simulation for Flow past a Stationary
Cylinder
As discussed before, two changes should be expected comparing 3D and 2D simulations: 1)
drop in the values of force coefficients especially the lift, and 2) existence of amplitude
modulations in time histories. Starting with the best option resulted from 2D simulation for
turbulence, the Standard and Realizable k – ε all y
+
models, in a domain with the thickness of 1D
leads to no significant difference between 2D and 3D simulations.
It was thought that the thickness of 1D is not sufficient, so the depth was increased to 3D and
no improvement was observed yet. No significant amplitude modulations were seen in time
histories either. Change of the top and bottom boundary conditions type did not help either.
Figure 4-5 – Drag (left) and lift (right) coefficients for the case of 3D simulation by the Standard k – ε
with the thickness of 3D and using 36 elements in z
Studying different quantities in the z direction shows that there is no variation in z axis and it
looks like modeling the third dimension is redundant. The other important aspect of 3D
simulations using the Standard and Realizable k – ε all y
+
models is the absence of stream-wise
vortices which should exist in the wake for Re = 4,000 according to Williamson (1996).
85
Figure 4-6 – Vorticity (left) and velocity (right) magnitudes in the depth of domain for the case of 3D
simulation by the Standard k – ε with the thickness of 3D and using 36 elements in z
Figure 4-7 –Stream-wise vorticity for the case of 3D simulation by the Standard k – ε with the thickness of
3D and using 36 elements in z
Comparing the 2D and 3D simulation results shows that no improvements or changes are
achieved by doing the computationally heavy 3D simulation:
Figure 4-8 – The lift coefficient resulted from the 3D (on left) and 2D (on right) simulations by the
Standard k – ε
Therefore, according to these simulations, the Standard and Realizable k – ε models are not
able to produce 3 dimensionality for flow past a cylinder. Three dimensionality, here, simply
86
means the existence of amplitude modulations in force time histories, variations of flow
parameters in z direction, significant stream-wise vorticity, and a drop in statistical values of the
lift and drag coefficient compared to 2D simulation.
However, switching to k – ω models and keeping all other parameters like the boundary
conditions and the z discretization the same, leads to three dimensionality. The lift and drag
amplitudes modulate significantly through time as shown below.
Figure 4-9 – Drag (left) and lift (right) coefficients for the case of 3D simulation by the k – ω SST Low y
+
treatment with the thickness of 3D and using 36 elements in z
Unlike the cases of using the Standard and Realizable k – ε models, flow properties vary in z
direction and stream-wise vortices have a significant magnitude compared to the total vorticity in
the wake.
Figure 4-10 – Vorticity (left) and velocity (right) magnitudes in the depth of domain for the case of 3D
simulation by the k – ω SST Low y
+
treatment with the thickness of 3D and using 36 elements in z
87
Figure 4-11 –Stream-wise vorticity for the case of 3D simulation by the k – ω SST Low y
+
treatment with
the thickness of 3D and using 36 elements in z
Finally, the lift and drag’s statistical properties drop in the 3D simulation.
Figure 4-12 – The lift coefficient resulted from the 3D (on left) and 2D (on right) simulations by the k – ω
SST Low y
+
treatment
Thus, observing three dimensionality for flow past a cylinder should depend on the
turbulence model. Star-CCM+ has several options for RANS models and this helps to run
different simulations with different models and see which one leads to three dimensionality and
which one does not.
Here, the results of 3D simulations with different turbulence model, depth, and z
discretization and comparison with the corresponding 2D simulation are presented in Table 4-2.
In the table L/D is the ratio of domain depth to the cylinder diameter and using a prismatic
mesh is denoted by Layered Mesh.
88
Table 4-2 – Results of 3D and corresponding 2D simulations for different turbulence models
Mesh Type L/D # of element in Z Turbulence Model CD,mean CL,rms Modulated?
2D -
SKE All y+ 0.8715 0.247 -
Polyhedral 1 - SKE All y+ 0.881 0.291 No
Layered 1 16 SKE All y+ 0.913 0.222 No
Layered 3 24 SKE All y+ 0.916 0.245 No
Layered 3 36 SKE All y+ 0.916 0.245 No
2D -
RKE All y+ 0.997 0.339 -
Layered 1 16 RKE All y+ 0.952 0.344 No
2D -
RKE All y+ w/ Curvature Correction 1.033 0.295 -
3D 3 36 RKE All y+ w/ Curvature Correction 0.993 0.212 No
2D -
KWSST Low y+ 1.65 1.345 -
Layered 3 24 KWSST Low y+ 1.31 0.679 Yes
Layered 3 36 KWSST Low y+ 1.259 0.643 Yes
2D -
KWSST Low y+ w/ Curvature Correction 1.695 1.393 -
Layered 3 24 KWSST Low y+ w/ Curvature Correction 1.257 0.654 Yes
89
Mesh Type L/D # of element in Z Turbulence Model CD,mean CL,rms Modulated?
2D - - RSM - Linear Pressure Strain - 2 Layer 1.23 0.644 -
Layered 3 36 RSM - Linear Pressure Strain - 2 Layer 1.243 0.661 No
2D
RSM Quadratic Pressure Strain - High y+ 1.425 1.1 -
Layered 3 36 RSM Quadratic Pressure Strain - High y+ 1.21 0.634 Yes
Layered 3 24 RSM - Linear Pressure Strain - High y+ 1.204 0.582 Yes
Layered 3 24
RSM - Linear Pressure Strain - 2 Layer w/
Daly_Harlow Gradient Diffusion
1.248 0.691 No
2D -
SKE - Low Re - Low y+ 1.085 0.459 -
Layered 1 16 SKE - Low Re - Low y+ 1.0375 0.45 No
Layered 3 24
SKE - Low Re - Low y+ w/ Cubic
Constitutive Relation
1.254 0.681 Yes
2D
StKW Low y+ 1.673 1.45 -
Layered 3 24 StKW Low y+ 1.34 0.8 Yes
2D
St Spalart Allmaras Low y+ 1.558 1.17 -
Layered 3 24 St Spalart Allmaras Low y+ 1.324 0.764 Yes
3D 3 24 AKNKE Low Re 1.025 0.179 Yes
3D 3 24 V2F KE Low y+ 1.312 0.828 Yes
3D 3 24 EB KE Low y+ 1.332 0.796 Yes
90
As it can be seen from Table 4-2, some of the 3D simulations show signs of three
dimensionality and some do not. Below is the summary of the 3D simulations with and without
three dimensionality.
Table 4-3 – Summary of models with or without triggering 3 dimensionality in a 3D domain
Models with 3 Dimensionality
Models without 3 Dimensionality
All St KW models
RKE - High y+ w/ and w/o Curvature
Correction
All SST KW models
SKE - High y+ w/ Linear Constitutive
Relation
All Spalart-Allmaras models
SKE - Low y+ w/ Linear Constitutive
Relation
RSM Linear Pressure Strain - High y+
RMS Linear Pressure Strain All y+
RSM Quadratic Pressure Strain - High y+
RSM Linear Pressure Strain - All y+ w/
Daly_Harlow Gradient Diffusion
AKN KE - Low Re - Low y+
RKE All y+
V2F KE - Low y+
SKE - All y+ w/ Linear Constitutive
EB KE - Low y+
SKE - Low y+ w/ Cubic Constitutive Relation
SKE - High y+ w/ Cubic Constitutive
Relation
Interestingly, all the models without any three dimensionality have an ε equation in common.
As mentioned before, the Reynolds Stress models need an extra equation to close the system of
equations and Star-CCM+ uses the models with the ε equations. In the following section, these
results will be discussed and the turbulence models will be examined in more details.
91
4-4- Discussion
There are some subjects that need to be addressed first to have a better insight of why some
models result in three dimensionality and some do not. But before that, it is worthwhile to
mention that the tabulated result suggests some statements regarding the discretization in z.
First of all, using polyhedral instead of prismatic mesh does not improve the result although
the cell size especially near the boundary layer is much more refined and the time step is smaller.
Comparing the polyhedral mesher and prismatic mesher with 16 layers in a domain with the
thickness of 1D shows that the minimum cell size is almost half in the polyhedral mesh meaning
that the time step should also be half in the case of polyhedral. However, the final result using
the Standard k – ε model with all y
+
treatment is worse for the polyhedral mesh.
Having a better refinement in z direction can help to get a better result. This is not the case of
the Standard k – ε model, however, because as shown before, the Standard k – ε has no three
dimensionality. But, using the k – ω SST model with 36 layers shows some improvements
compared to 24 layers with the same thickness and the statistical properties of the lift and drag
drop slightly. This improvement is not very significant considering the increase in cell count (1.5
times) as Prsic and Myrhaug (2012) showed for the case LES as well.
RANS models have some characteristics that can affect the result and in this case may or
may not trigger three dimensionality. Here, these characteristics are presented:
4-4-1- Isotropy
Isotropy comes into RANS models in some ways:
1) The turbulence kinetic energy (k): Using k in one- and two-equation models introduces
isotropy. Basically k can be computed by taking the trace of the Reynolds stress tensor:
92
ii i i
w v u u u k
2
1
2
1
2
1
2 2 2
( 4-1)
The k equation in k – ω and k – ε models makes these models isotropic which means that
these models may not capture the anisotropic features of flow. In fact, one of the reasons for
developing Reynolds Stress models was to introduce anisotropy by solving for all individual
Reynold stress components instead of k only.
2) Eddy viscosity model: In the Boussinesq approximation, Reynolds stresses are related to
the mean strain rate tensor, S
ij
, by the turbulent eddy viscosity,
T
, in one- and two-
equation models. The eddy viscosity is originally an isotropic parameter meaning that it
has the same value for all the stress components. Sometimes the term linear constitutive
relation is used for the original Boussinesq approximation.
Even using the turbulent eddy viscosity
T
in these models involves isotropy since there is
no reason to assume that the eddy viscosity is the same for all directions.
Some efforts have been done on introducing anisotropy to eddy viscosity model which are
called nonlinear constitutive relations. For example, Star-CCM+ has the following quadratic
constitutive relation as an option for k - ε models (CD-adapco 2016):
T T
t
T
t t Linear t Quad t
k
C
k
C
k
C
W W I W W
W S S W S S I S S T T
:
3
1
4
4 :
3
1
4
3
2 1 , ,
( 4-2)
where,
I v S T k
t t Linear t
3
2
2
,
( 4-3)
S and W are the mean strain rate and vorticity tensors respectively. C
1
, C
2
, and C
3
are the
model coefficients.
93
The cubic constitutive relation is also available and used in some of the simulations to
investigate the effects of this type of nonlinearity in the results.
3) Dissipation: Dissipation occurs at the smallest scales of turbulence, so according the
Kolmogorov (1941) hypothesis, dissipation is locally isotropic (Wilcox 2006). This is the
basis of many ε based models. Even in Reynolds Stress Models the same assumption is
used in most of the ε based models. According to Wilcox (2006) in Reynolds Stress
models:
ij ij
3
2
( 4-4)
in which,
__________
k
i
k
i
x
u
x
u
( 4-5)
However, in reality dissipation is anisotropic especially near the wall and some researchers
have tried to include some sort of anisotropy in ε. For example for generalizing a low Reynolds
number model it is proposed that (Wilcox 2006)
ij s ij ij
b f 2
3
2
( 4-6)
where b
ij
is the dimensionless Reynolds stress anisotropy tensor and f
s
is a low Reynolds
number damping function.
k
k u u
b
ij j i
ij
2
3
2
( 4-7)
Therefore, there are many sources of isotropy in RANS model which can affect the solution
for simulating a basically anisotropic problems, like many 3D flow regimes.
94
4-4-2- Wall Treatment
Originally, turbulence models, i.e. k – ε, were devised for high Reynolds number and were
valid and accurate only outside the viscous-dominated regions, i.e. near the wall. So, a set of
mathematical expressions are required for near the wall behavior. There are three categories for
wall treatment that vary depending on the type of turbulence model and application. These
categories in Star-CCM+ and many other CFD packages are the high y+, low y+, and all y+
treatments.
- High y
+
treatment: Star-CCM+ uses a wall-type approach for high y+ treatment. Here it
is assumed that the turbulence model is valid for only the fully turbulent regions and the
viscous affected region (the viscous sublayer) is not resolved. For this purpose, the
domain should be meshed in a way that the centroid of the 1
st
cell near the wall is located
in the logarithmic region of the boundary layer (CD-adapco 2016). This implies that y
+
values should lie between 30 and 50 but the values up to 200 are still good. Since the
viscous sublayer is not modeled, using this type of treatment results in a faster solution
and a coarser mesh which usually converges faster. Obviously, high y
+
treatments should
be used when resolving the viscous sublayer is not necessary or important.
- Low y
+
treatment: This treatment is suitable only for low Reynolds number models and
assumes that the entire boundary layer including the viscous sublayer is properly resolved
(CD-adapco 2016). Therefore, the near wall mesh needs to be well refined and the value
of y
+
must be as low as possible (preferably less than one but values between 3 and 5 are
still acceptable). Because the flow quantities are integrated down to the wall, this requires
a lot computational efforts which makes this approach more expensive. On the other
hand, this is the desirable method when dealing with viscous-turbulent transitions. It
95
should be noted that since the viscous sublayer is very thin using a low y
+
treatment for
high Reynolds number requires a very fine mesh and is not practical.
- All y
+
treatment: The all y
+
treatment in Star-CCM+ is a blending treatment that uses a
low y
+
treatment when y
+
values are small (y
+
→ 0) and switches to a high y
+
treatment
for higher y
+
(y
+
> 30). It is claimed that for the values of y
+
between the viscous sublayer
and the logarithmic layer, i.e. the buffer layer, the method has an acceptable accuracy
according to CD-adapco. In most cases, the term “two layer” can be interpreted as using
an all y
+
treatment in Star-CCM+.
4-4-3- Low-Re Models
Low Reynolds number models are needed for resolving the viscous sublayer. The k – ω
family and Spalart-Allmaras models are designed such that they can be applied throughout the
boundary layer. In contrast, the k – ε models require more modification for low Reynolds
number. This is also the case for Reynolds Stress models.
There are several low-Re k – ε models and four are available in Star-CCM+ namely the
Standard Low-Re k – ε, the Abe-Kondoh-Nagano k – ε Low-Re, v
2
f - k – ε and EB k – ε. The
low Reynolds number version of Standard k – ε and the Abe-Kondoh-Nagano k – ε Low-Re
models use damping functions to modify the model’s coefficient (CD-adapco 2016). These
damping functions are basically functions of a turbulence Reynolds number, Re
T
≡ k
2
/(ευ)
(Wilcox 2006), and include the distance from the wall (CD-adapco 2016). However, the v
2
f - k –
ε and EB k – ε models do not need any damping functions and do not use wall distance.
Badically, using a wall function that is a function of wall distance introduces some level of
anisotropy to the model because wall distance functions act directionally normal to the wall. The
96
Abe-Kondoh-Nagano k – ε Low-Re was devised by Abe, Kondoh, and Nagano (1995) for heat
transfer and reattachment at low Reynolds number. It implements a wall function for the (linear)
eddy viscosity term and as shown before, the results for flow past a cylinder are the best among
all the simulation.
The v
2
f – k – ε , developed by Durbin (1991), is a k – ε model with 2 additional equations;
one for
2
v
and the other one for an elliptic relaxation function, f. It is originally designed for low
Reynolds number and unlike the Standard k – ε model does not require any wall functions and is
capable of integrating throughout the viscous sublayer (CD-adapco 2016). Solving for
2
v
independently means that this model is anisotropic near the wall. The EB k – ε model (Billard &
Laurence 2012) is basically a version of v
2
f – k – ε with elliptic blending functions and
technically anisotropic (CD-adapco 2016).
Unlike k – ε models, k – ω and Spalart-Allmaras models do not need damping functions and
they are able to resolve the viscous sublayer without wall distance functions. There is a low
Reynolds number k – ω model based on the Standard version which uses some modified
coefficients and is available is Star-CCM+ (CD-adapco 2016).
4-5- Justifications
By comparing the results from different turbulence models it seems that two factors play
significant roles in triggering three dimensionality: Wall y
+
treatment, and isotropy. The results
here are tried to be justified by addressing the following questions:
97
4-5-1- Why Do the Standard and Realizable k – ε Models with all y
+
Wall
Treatment Show no Three Dimensionality?
The all y
+
treatment for the range of y
+
< 1.0 acts like the low y
+
treatment. This means that
both models will use damping functions. Most of the k – ε models are notorious for producing so
much damping, especially the Standard version (Wilcox 2006). So, the models seem to be so stiff
that they do not allow any three dimensionality development or damp the three dimensionality.
Beside this, both models are isotropic, so isotropy may also prevent three dimensionality.
The Standard k – ε low-Re low y
+
treatment seems to perform with the same characteristics.
Although the model is revised for low Reynolds numbers and low y
+
values, the damping
produced by damping functions seems to eliminate three dimensionality.
4-5-2- Why Do the AKN Low-Re Low y
+
, v
2
f – Low y
+
, and EB – Low y
+
k – ε
Models Produce Three Dimensionality?
As discussed, the v
2
f and EB k – ε models are anisotropic and do not use any damping
functions for low y
+
treatment. So it seems like they are capable of triggering three
dimensionality. AKN low-Re k – ε, on the other hand, has damping functions. However, this
model is designed for such problems with low Reynolds number and strong separations and
reattachments. Thus, there seems to be no wonder why using AKN k – ε results in three
dimensionality.
4-5-3- Why Does the Reynolds Stress Model with Two Layer Treatment Have
no Three Dimensionality?
The Reynolds Stress models are anisotropic in a sense that they solve for all Reynolds stress
components and do not use any constitutive relations. However, they usually have many
common features with the Standard k – ε model including the ε equation and wall treatments.
98
Therefore, it seems that using the Reynolds Stress model for y
+
< 1.0 produces unnecessary
damping which seems to stop developing three dimensionality. Even, using the Daly-Harlow’s
generalized gradient diffusion model instead of the isotropic model for the turbulent diffusion
term does not help to initiate three dimensionality. The pressure-strain correlation is an important
term in Reynolds Stress model, but Star-CCM+ does not have any options for using the quadratic
correlation for the low y
+
treatment and the user can select it with a high y
+
treatment only.
4-5-4- Why Does the Nonlinear Constitutive Relation Help to Develop Three
Dimensionality?
Choosing a cubic constitutive relation instead of the Boussinesq approximation triggers three
dimensionality in the original and low-Re Standard k – ε models. Although damping functions
are still in use, but apparently the nonlinearity in calculating the stress tensor seems to dominate
and turn three dimensionality on. This could show the importance of anisotropy and nonlinearity
when dealing with flows with significant 3D effects.
4-5-5- Why Do all Versions of the k – ω and Spalart-Allmaras Models Show
Some Level of Three Dimensionality?
It was observed in simulations that the results of all versions of the k – ω SST, Standard k –
ω, and Spalart-Allmaras models have three dimensionality. All of these models are isotropic and
no nonlinear terms are included. However, being able to resolve the whole boundary layer
without any wall functions or damping functions seems to allow three dimensionality. This
feature is very important considering that k – ε models require two more equations in order to be
able to integrate throughout the viscous sublayer without damping functions and they are more
expensive.
99
All unsteady 3D simulations are initiated after a steady run for reducing the computing time.
To make sure that the disturbance induced by this steady state run does not introduce unreal 3D
structures, for one case, a k – ω SST simulation is initiated from zero in which the velocity field
is increased to the ambient velocity in 10.8 seconds. This simulation does trigger three
dimensionality as well, but the lift and drag coefficients are higher. Thus, observing three
dimensionality seems to be depending on the turbulence model not initial or boundary
conditions.
4-5-6- Why Is Three Dimensionality Observed in Some of the Models with
High y
+
Treatments?
Interestingly, in some cases using high y
+
wall treatments helps to have three dimensionality
like the case of RSM with or without a quadratic pressure-strain correlation, or the Standard k – ε
with a cubic constitutive relation. However, it should be noted that when the y
+
values are low
and always less than 1.0 in all simulations, using a high y
+
treatment is physically invalid. High
y
+
treatments assume that the centroid of the 1
st
prism layer lies in the log layer (30 < y
+
< 50)
and they use an analytical solution at the wall. So, using a high y
+
treatment for y
+
< 1.0 imposes
a logarithmic velocity profile in a region where experimental measurements show that the
velocity changes linearly with the wall distance (u
+
= y
+
). Therefore, the results of any high y
+
wall treatment should not be included in any conclusions.
So, apparently, damping functions in the majority of k – ε models and also the Reynolds
Stress models seem to be a significant obstacle for the model to trigger three dimensionality.
Using higher order approximations like nonlinearity in the constitutive relation seems to help to
overcome the deficiencies caused by damping functions. It seems like these nonlinearities
100
neutralize the effects of damping functions and release three dimensionality. Anisotropy can also
facilitate the development of three dimensionality if damping functions are not implemented.
In summary, it should be noted that for problems with severe streamline curvature,
substantial separations and reattachments, and laminar to turbulent flow transient, using RANS
models, especially those originally devised for fully turbulent flows and employ wall functions,
is not the best practice and other turbulence resolving approaches should be implemented.
4-6- Back to 2D Simulation of Flows past a Cylinder
After investigating the performance of different turbulence models in 3D simulation of flow
past a stationary cylinder, some analogies could be drawn for 2D simulations of flow past a
stationary and moving cylinder.
4-6-1- Effects of Turbulence Models on the Stationary Cylinder Results
By comparing the results for different model a trend can be drawn. The values of lift and
drag are constantly increasing with this order: The Standard k – ε all y
+
, then the Realizable k – ε
all y
+
, the Standard low-Re k – ε low y
+
, the RSM two layer. All these models use damping
functions for resolving the viscous sublayer. Interestingly, the error for all other models, those
that have no damping functions, is much higher than this category and the value of lift is
considerably high.
Therefore it seems that in 2D simulation of flow past a stationary cylinder, damping
functions play a significant role as well.
101
4-6-2- Effects of Turbulence Models on the Moving Cylinder Results
Here, the effects of damping are much more pronounced; the Standard k – ε all y
+
model
results in a shedding pattern that looks like a stretched 2S mode. This means that the damping of
this model seems to be so high that it does not allow to generate vortices in pair unlike other
models. The Realizable k – ε all y
+
and RSM lead to a better shedding mode compared to other
models and as discussed before the curvature correction helps to improve the value of lift
although it does not improve the shedding pattern.
The time history of lift for different model also indicates the effects of damping; the lift for
the Standard and Realizable k – ε all y
+
models shows no amplitude modulations and in the case
of the Realizable model with curvature corrections the signal looks completely damped. It seems
like damping functions here prevent the lift amplitude from modulating. Even though, the result
of RSM has some modulations but these modulations are less compared to the k – ω models.
Again the same trend in the values of lift can be observed for the moving case.
4-7- Concluding Remarks
2D simulations are not able to produce any 3 dimensionality, thus moving to a 3D domain
should bring about 4 differences which all are covered as a concept called ‘three dimensionality’.
Three dimensionality in vortex shedding implies 1) a drop in the values of drag and lift, 2)
significant amplitude modulations in time traces, 3) variation of different quantities in z
direction, and 4) existence of stream-wise vortices in the wake.
Running the simulations with the same planar mesh and different discretization methods in z
allows to investigate the effect of each parameter in results. Simulations show that the z
refinement does not alter the final result significantly considering the extra efforts resulted by
102
refinements. However, the domain thickness (size in z direction) and the type of boundary
condition on top and bottom can be influential according to the literature, and the current study
confirms that the thickness of 1D is not sufficient to capture the 3D effects of flow.
Studying the result of different RANS models indicates that three dimensionality can and
cannot be triggered depending on the turbulence model. Two major features of RANS seem to
play a role here, one, the isotropy and nonlinearity of the model, and two, wall treatments and
using damping functions to resolve the y
+
< 1 regions. All those models with no three
dimensionality, the Realizable and Standard k – ε all y
+
treatment, the low-Re Standard k – ε
low y
+
, and RSM two-layer with linear pressure-strain correlation, have damping functions for y
+
< 1.0 . Damping seems to be so much in these models that does not allow any 3D structure to
form or develop.
However, the cases of V2F k – ε low y
+
, and the Standard k – ε all y
+
with a cubic
constitutive relation show that anisotropy/nonlinearity seems to be able to deactivate the
damping functions impacts and three dimensionality can be released. AKN k – ε low y
+
model
produces the closest result to the experiment and has three dimensionality. Damping functions
are tweaked for this model in a way that the model is capable of treating 3D flows and
separations well.
Damping in low y
+
treatments of k – ε models seems to be responsible for getting lower lift
coefficients for 2D simulations of the stationary cylinder, as well as the absence of amplitude
modulations in lift for the moving cylinder vortex shedding. Simulations show that the result of
lift is much lower for the Realizable and Standard k – ε all y
+
models and this could be a result
of the damping functions. This trend is also evident for the moving case. Moreover, the shedding
pattern and lift response of these two models are much more damped compared to other models.
103
In the end, 3D simulation can have some improvements over the 2D simulation considering
the features and capabilities of the turbulence model. RANS models in general are not suitable
for 3D flows with considerable streamline curvatures, large separations and reattachments, and
transitions. LES and DNS are proven to be able to capture these phenomena with a very high
cost. The hybrid LES/RANS like DES seems to an intermediate solution for the current era,
however, this study shows that using RANS models with damping functions may not be a good
choice for hybrid methods.
104
Chapter 5- Analysis of Flow past a
Stationary Cylinder by 3D Simulations
and DES Models
5-1- Introduction
Flow past a moving cylinder under a controlled oscillation is the main subject of interest in
this study. However, keeping in mind that the phenomenon is a very complex fluid mechanics
challenge and requires a lot of computational resources, it would be beneficial to reduce the
complexity and first try to simulate a flow passing a stationary cylinder with the same Reynolds
number. This would help to come up with a configuration for flow domain, mesh structure,
spatial resolution (mesh refinement), temporal resolution, turbulence modeling, boundary
conditions, and other details for the solver.
In previous chapters, the same approach was pursued as well. In chapter 3, flow past a
cylinder with a 2D mesh was studied and it was shown that 2D simulations have some physical
and numerical limitations for flow past a stationary and moving cylinder. Chapter 4 was
specified to 3D simulation of a stationary cylinder using RANS closures. As discussed in chapter
4, RANS models have some drawbacks when dealing with a complex turbulent flow that consists
of transitions and massive separations and reattachments near boundaries.
105
In this chapter there is going to be an effort to demonstrate the capabilities of detached eddy
simulation (DES) in its simplest form for capturing the characteristics of flow past a cylinder.
5-2- Three Dimensional Simulation of Flow past a Stationary
Cylinder
DES is the main subject in this section. Before presenting and discussing the results, first the
simulation properties such as domain size, meshing, and time step along with DES model details
are illustrated.
5-2-1- Flow Domain and Boundary Conditions
In previous chapters a rectangular domain with prism layers and planar polygonal cells was
used for simulations. It was decided for the rest of study, a circular shape with a structured mesh
is a better option especially for the case of moving cylinder.
A polygonal/polyhedral mesh, because of having more cellular faces, is claimed to be able to
capture the flow characteristics better. However, keeping in mind that an oscillation will be
imposed to the mesh later, this type of mesh can add more complexity to the problem, so it
would be better to go with a more structured mesh. Beside this, structured meshes use and need
less memory and computational resources which is important for a moving mesh problem.
After reviewing the literature a size of 30D×πD was chosen for simulation. This is a very
typical size for such a CFD problem and has been used by other researchers with a lot of success.
For the boundary conditions, the left half of domain is selected as a velocity inlet with
components of [U,0,0], the right half as a pressure outlet that applies the pressure outlet
106
boundary condition only in the direction of flow (x direction). The top and bottom of cylinder
(top and bottom boundaries) are defined as a periodic boundary, and the cylinder is set to behave
like a wall with no slip boundary condition.
Figure 5-1 – Circular domain size used for simulation
Figure 5-2 – Boundary conditions of the simulation domain – top view
107
Figure 5-3 – Boundary conditions of the simulation domain – 3D view
As discussed above, a structured (direct) mesher is used to mesh the domain. Instead of x and
y for the planar mesh to have a better control on mesh generation, r and θ are used. In order to
create a structured mesh as simple as possible, discretization in z and θ are uniformly distributed,
but for the radial discretization 3 regions are built and in each the size in r increases with a one-
sided geometric growth. This helps to produce less discretization in r with a desired refinement
near the cylinder wall where it is already known that radial discretization plays an important role
in capturing flow characteristics. This will be discussed in more details later.
Therefore, to get most benefits of the direct mesher, the r-z plane is meshed with the desired
discretization first (uniform in z and three geometric-growing zones in r), and then this 2D mesh
is revolved in θ to create a 3D mesh.
108
Figure 5-4 – Structured mesh in the r-z plane with a geometric growth in r
There are two factors that are considered for the radial discretization. First and most
important is the 1st layer size (thickness). The first cell radial size, as discussed in previous
chapters, directly controls the wall y+ value. It should be noted once again that since for the
studied Reynolds number the boundary layer is still laminar, it is tried to keep the value of wall
y+ as minimum as possible and very close to 1.0. This will make sure that the solver is able to
capture the entire boundary layer without any further numerical assumptions.
The other parameter for the radial discretization is that it is tried to keep the cell aspect ratio
in the r-θ plane between 0.5 and 2 and for the majority of plane close to 1. Having a bounded
aspect ratio can be used as a guarantee that the vortical structures are captured with a better
accuracy. If defined as θ/r, considering the 1
st
thin layer, the aspect ratio is highest next to the
wall, and lowest at the inlet and pressure boundaries. However, three mesh zones in r are tuned
in a way to have aspect ratios close to 1 for the most important zones of the mesh.
Starting from Zone 1 and knowing the 1
st
cell size and the zone size, a series of repetitive
calculations can be done for the geometric growth rate (r
geom
).
N
geom geom
r r
/ 1
1
( 5-1)
109
Where δ is the size of zone, α is the first cell size, and N is the number of discretization in
that zone. For Zone 2, the 1
st
cell size can be computed from the Zone 1 last cell size and its
growth rate. Similar calculations are valid for Zone 3.
In Table 5-1 an example of mesh calculations are presented for the 3D o-shape domain.
Table 5-1 – An example of calculations for the 3D mesh
Considering the flow pattern, it would be thought why not using a non-uniform discretization
in θ with a more concentration around a line parallel to the flow passing through the center of the
cylinder (θ = 0). Although there are more complex vortical flow structures due to vortex
generation around θ = 0 than other regions for example θ = 180°, but having a uniform mesh in
θ will be more beneficial for the moving cylinder case where changes in mesh and cell quantities
Domain Radius (*D) 15
Spanwise Size (*D) π
Circumferential Discretization 314
Number of Zones 3
Spanwise Discretization 60
Spanwise Size (*D) 0.052
Planar Cell Count 43,960
Radial Discretization 140
3D Cell Count 2,637,600
Zone Size (*D) 0.25
1st Cell Radial Size (*D) 0.005
Radial Discretization 25
Geometric Growth Rate 1.053
1st Cell Planar Aspect Ratio 2.00
Last Cell Planar Aspect Ratio 0.87
Zone Size (*D) 4.75
Radial Discretization 85
Geometric Growth Rate 1.023
1st Cell Planar Aspect Ratio 0.83
Last Cell Planar Aspect Ratio 0.90
Zone Size (*D) 9.5
Radial Discretization 30
Geometric Growth Rate 1.057
1st Cell Planar Aspect Ratio 0.88
Last Cell Planar Aspect Ratio 0.48
Name: DES_05
Zone 1
Zone 2
Zone 3
110
for each time step and more importantly generating a new mesh at each time step can be captured
more efficiently with less computational efforts.
Figure 5-5 – Planar mesh in the r-θ plane
The discretization in θ is chosen in a way to control the aspect ratio of cells closest to the
wall. Prism layers are not used anymore so it is not feasible to have high aspect ratios near the
wall. It should be noted that having prism layers provides a mesh with less circumferential
discretization. However, prism layers essentially introduce approximations to the model (i.e.
changes in normal direction to the wall are much drastic than parallel to it) and for flow past a
cylinder where the wall is curved and high gradients of velocity and pressure along with severe
separations exist at the wall, it would be more accurate to drop prism layers and use a direct
mesh instead.
111
Figure 5-6 – Details of mesh for flow past a stationary cylinder with a structured 3D mesh
The final 3D mesh will look like Figure 5-7
Figure 5-7 – Three dimensional structured mesh for simulating flow past a cylinder
Once again there are two controls for each simulation, 1) the wall y+ values are kept close to
one for the entire simulation time, 3) the convective CFL number should be limited to one. It
112
should be noted that since a directed mesh with much smaller cell sizes are used in this new
mesh and flow patterns are captured by a turbulence model with a higher fidelity, the time step to
keep CFL close to 1.0 will be much smaller.
Having much more cell count, a smaller time step, and a more accurate turbulence model will
all result in high computational costs. In fact, for solving such a CFD problem, regular computers
are not useful anymore and supercomputers with large memories are more favorable. All
simulations are run on the USC supercomputer. More details on high performance computing
will be discussed in the next chapter.
For the turbulence model, it is tried to use the simplest RANS closure of detached eddy
simulation (DES) which is the Spalart - Allmaras IDDES. For some cases, the SST k-ω model
has been as well. There have been some other studies in the literature that LES models are
implemented in their simulations with quite acceptable results, however, the main goal in this
section is how to improve the literature by using less cell count and a turbulence model with a
lower fidelity (DES vs LES). Considering the high cost of such simulations, coming up with a
simulation with lower cost is a substantial step. It should be noted that RANS models alone are
shown, in Chapter 4, to be unable to capture the flow patterns and predicting the statistical
parameters of flow (i.e. mean drag coefficient, RMS value of lift, and Strouhal number)
accurately.
Employing a second order solver, simulation proceeds to next time step after a certain
number of iterations. It is observed that the residuals of model variables converge to constant
values after some iterations meaning that more iterations do not improve the calculations.
Therefore, after running the simulation for few time steps with high inner iterations, a fix number
is set for the rest of time. This fix number depends on the temporal resolution (time step), but, for
113
the stationary case, it is set to a number between 20 and 50. It should be noted that for a turbulent
transient flow getting to a constant residual value is not always a case, thus, it is a rule of thumb
that a threshold is usually picked so if the residual value drops to a certain level, e.g. %5 of the
initial residual, iterations stop.
All simulations, like the case of RANS models, are first initiated by a steady state run. Since
for steady runs, only RANS models are in access, the same RANS closure as in the DES model
is selected and the simulation is run for about 2,000 iterations to get to a relatively developed
wake condition where some vertical structures are formed after the cylinder. This method helps
to save a lot of CPU time, because running DES from a still condition requires some shedding
time to achieve a quasi-steady response. For the case of 2D RANS model it was shown that the
error corresponding to this approach is negligible, for the 3D DES simulations, one case was run
from still conditions, and while the simulation takes about 50% more time to get to a quasi-
steady state, the error for statistical values is not that significant. In fact as it is shown in
Figure 5-8 the lift resulted from both approaches behaves very similar. The Strouhal number
error is very negligible is less than 2%.
Figure 5-8 – Different in lift coefficient time trace resulted from DES between a case starting with a
steady RANS run and a case starting from still conditions
114
As it is expected, DES is a much better tool compared to 3D URANS and 2D URANS for
capturing flow patterns behind the cylinder. In fact, a pilot run on the same simulation domain
with only 350,000 cells, which is quite coarse for a DES, with the Standard Spalart-Allmaras
simulation indicated that DES reduces the lift by more than 80% compared to a simulation with
almost the same cell count using URANS with the Standard Spalart-Allmaras model (CL
RMS
~
0.135 vs 0.764 from DES and URANS respectively). This reduction is shown in Figure 5-9
where the lift coming from both simulations are plotted on almost the same scale. The same DES
simulation also results in a reduction of 25% of the lift compared to the best 3D URANS model
presented in Chapter 4.
Figure 5-9 – The lift coefficient resulted from the DES S-A with 351,000 cells (on left) and the 3D S-A
URANS Low y
+
treatment with 358,410 cells (on right)
DES simulations also have all features of flow discussed in Chapter 4, where it was shown
that flow past a cylinder is a 3D phenomenon and flow quantities should have fluctuations in the
spanwise (z) direction. Moreover, lift and drag forces should experience modulations in time.
115
Figure 5-10 – Lift and drag coefficient resulted from DES with the Standard Spalart – Allmaras
Figure 5-11 – A section of velocity and vorticity magnitudes resulted from DES with the Standard Spalart
– Allmaras
Figure 5-12 – Spanwise vorticity field resulted from DES with the Standard Spalart – Allmaras
116
Each simulation is run to get to a quasi-steady state in lift and drag. Then the simulation is
stopped after some number of vortex shedding (between 20 and 50). To compare the results with
the test data, time histories are truncated and a Matlab script is used to calculate the statistical
value of drag and lift. Part of the same script operates an fft operation on the lift time trace to
compute the dominant frequency of shedding. This frequency is later used to calculate the
Strouhal number.
U
D f
St
dominant
( 5-2)
Figure 5-13 – Matlab script to calculate the statistical values of flow past a stationary cylinder
The frequency spectrum of lift shows that the flow has one distinct dominant frequency. This
frequency corresponds to vortex generations behind the cylinder.
117
Figure 5-14 – Frequency spectrum of lift force for flow past a stationary cylinder
5-3- Results of 3D DES Simulation
Having an extensive background in 2D and 3D simulations with URANS models, getting to a
proper mesh to match the experimental data should be easier. A combination of 1st cell size, z-
and θ discretization along with efforts on keeping the planar aspect ratio bounded lead to a
relatively low cell count and large time step which both are major factors for solving such a
challenging fluid mechanics problem.
118
Table 5-2 – Details of mesh parameters for flow past a stationary cylinder
119
In Table 5-2 details of calculations for different mesh setups used to get to a proper
configuration are presented. It is tried to see different parameters’ effects, such as θ and z
discretization, number of discretization in r, and growth rates in r, on the final lift and drag
values. Once again, it should be noted that the main focus in this section is using the Standard
Spalart Allmaras as the RANS closure, however, some simulations are run with the SST k-ω as
well to only see the effects on the final answer as well as control parameters (the CFL number
and wall y
+
values). The top and bottom boundaries are set to a period boundary, but symmetry
planes are also used for some pilot runs to see the influence. Since having a low wall y
+
has been
one of the controls, for the majority of simulations a Low y
+
model for the RANS closure is
used.
Table 5-3 – Parameters of simulations and final results for flow past a stationary cylinder
According to Table 5-3, the Strouhal number has a very low error for all simulations. This
indicates that DES can be used as a reliable approach for solving flow past a stationary cylinder
for this range of Reynolds number. Looking at lift and drag, unlike the 2D and 3D RANS
120
simulations the values are lower than the experimental data which may imply that DES
formulations damp the turbulent viscosity built up in the flow field. Comparing the results from
different mesh configurations, it seems that a more mesh cells does not necessarily makes the
simulation more accurate in terms of predicting the statistical values of flow parameters. This is
the case of DES_04 which has the highest cell count. The z discretization seems to also play an
important role to improve the result. This can be seen as the difference between DES_05 and
DES_07 which both have the same planar mesh but in DES_07, z discretization is lower. Since it
has been tried to keep the planar mesh aspect ratio bounded between 0.5 and 2, it would be hard
to investigate the impacts of r and θ discretization. However, it should be noted that for having a
better normal-to-wall resolution the growth rate in r should be kept as small as possible
particularly near the wall. Therefore, it would be beneficial to keep a balance between r and θ
discretization. The other important factor that should be emphasized again is the role of 1
st
cell
thickness which, as discussed before, determines the value of y
+
at the wall.
The influence of the SST k-ω model in DES was investigated for some cases. The most
important impact of this model on results is the change of wall y
+
values which are consistently
higher compared to the Standard Spalart-Allmaras model. This means that, if wall y
+
values and
planar aspect ratio values are assumed to be important, using the SST k-ω model will need
higher resolution in r and consequently in θ. For the case of wall y
+
value, it should be noted that,
as shown in Figure 5-15, the difference between value of 1 and 2, especially in RANS models, is
insignificant. In fact, it is a rule of thumb in the community that for resolving the boundary layer,
it is certainly great to have wall y
+
values limited to 1 but values up to 5 can be still acceptable.
However, it should be kept in mind that the mesh will later will be used in a moving cylinder
case, where wall y
+
values will be definitely higher because of the cylinder motion. Therefore,
121
even though, results of the SST k-ω model, in the case of DES_00_SST, are in a good agreement
with the experimental data, it was decided not to use this closure for the rest of research as the
main model.
Figure 5-15 – Plot of u
+
vs y
+
for a typical boundary layer (Wilcox 2006)
Comparison of the current CFD results with some previous numerical studies on the same
range of Reynold number shows that the current DES simulation (DES_05) works with relatively
low computational costs in terms of the turbulence modeling, mesh count, and even the time
step. For comparison the CFD results of Kim et al. (2014), Thingbø (2013), Lysenko et al
(2012), Prsic & Myrhaug (2012), Pereira et al (2015), and Wornom et al (2011) are listed in
Table 5-4. As Table 5-4 shows the current DES simulation is among the lowest for the cell count
and the time step. More importantly, the other studies have implemented LES models which are
more expensive than DES. This can imply that, if DES is used in conjunction with a tuned mesh
it can be quite capable of predicting the statistical values of forces for flow past a stationary
cylinder in the transient Reynold number range.
122
Table 5-4 – Comparison between the current study results and model parameters with most recently CFD studies
5-4- Concluding Remarks
In conclusion for this part, it has been shown that, using DES over URANS has a clear
superiority for flow past a stationary cylinder at Re = 4,000. It has been also shown that, if mesh
parameters, meaning 1
st
cell thickness, discretization in θ and z, and the planar aspect ratio, are
selected properly, DES could be a great solution for accurately predicting flow characteristics for
this range of Reynolds number when the computational costs of simulations are considered and
compared to LES which is already proven to be a reliable solution for such a CFD challenge.
Results of this chapter helps a lot with coming up with a proper set-up for the model to
simulate the flow past a moving cylinder which will be covered in the next chapter.
123
Chapter 6- Analysis of Flow past an
Oscillating Cylinder with Three
Dimensional Simulations
6-1- Introduction
This chapter focuses on flow past a cylinder when the cylinder is moving with a prescribed
sinusoidal form. The goal is to evaluate capabilities of Star-CCM+ in simulating such a velocity
field with more sophisticated turbulence models. For this purpose, experimental data from Morse
and Williamson (2009) will be used to verify the CFD simulations.
As mentioned before, Morse and Williamson (2009) came up with a new distinct vortex
shedding pattern, which they called 2P
O
, which possesses some interesting characteristics. They
showed that this new shedding mode is an overlap between the 2P and 2S modes, and depending
on the amplitude and frequency of oscillation, it can be stable and sustain for some time. Their
experimental data also shows that, for a constant amplitude and frequency, flow characteristics
switch between 2P
O
mode and 2S or 2P modes in time. More importantly, the 2P
O
mode is
associated with the highest possible amplitude for having a freely oscillating cylinder at Re =
4,000.
In this chapter simulation results based on the domain and mesh used for the stationary
cylinder will be discussed. For applying the motion to the cylinder two techniques are used and
124
several DES and LES based models are employed to simulate the flow past an oscillation
cylinder in the 2P
O
region of Morse and Williamson’s study for Re = 4,000.
Figure 6-1 – Map of vortex shedding regimes in Morse and Williamson’s studies
6-2- Model Parameters
Simulating flow past a moving cylinder requires many details and parameters to adjust
properly. However, the stationary cylinder case can help to resolve the majority of these
parameters. In fact, setting two controls for the model (wall y
+
~ 1.0 and the convective CFL
number < 1.0) along with the statistical flow characteristics at Re =4,000 (meaning CL
rms
,
CD
mean
, and Strouhal number) sort of guarantee that the parameters can be applicable to the
moving case as well.
The stationary cylinder simulations showed that having around 2.6 million cells in an o-shape
domain of 30D in diameter and thickness of πD with first layer thickness of 0.005D provided a
good spatial resolution for this problem. Regarding the temporal resolution, ∆tU/D = 0.005
guaranteed that CFL < 1 for the stationary case. Keeping in mind that applying an oscillation to
125
the cylinder will increase the magnitude of velocity in the system, the time step for the moving
case should be smaller to maintain CFL < 1.0.
It is also worth noting that higher velocities due to the cylinder oscillation will potentially
increase the wall shear stress and consequently the wall y
+
value. However, wall y
+
values
around 2 to 3 can still be acceptable for a proper resolution of the boundary layer near the
cylinder wall. This is illustrated in the plot of u
+
vs y
+
where as it can be seen, the deviation from
the line of u
+
= y
+
for values of y
+
up to 5 is insignificant.
Figure 6-2 – Plot of u
+
vs y
+
for a typical boundary layer (Wilcox 2006)
For turbulence modeling, a DES model based on the Standard Spalart-Allmaras closure was
shown to lead to acceptable results in terms of statistical values of flow properties. Therefore,
this model can be a good starting point for the moving cylinder problem. Boundary conditions
were also tested by the stationary cylinder simulation. The only difference for simulation of a
moving cylinder is how to incorporate the motion of the cylinder in the model.
126
6-2-1- Applying the Cylinder Oscillation in Modeling
As discussed before, for applying a motion to the domain there are three techniques available
in the software. These techniques are:
1) Moving Reference Frame (MRF):
This technique was used in the 2D simulation presented in Chapter 3 before. According to
the Star-CCM+ documentation, in MRF, it is required to have a reference frame, which can be
either fixed or moving itself, so the mesh motion is then defined with respect to that reference
frame. If it exists a moving reference frame, the velocity of a material point P with respect to the
moving reference frame will be:
t MRF r ,
v v v ( 6-1)
Where
r
v is the relative velocity with respect to the moving frame, v is the velocity in the
fixed reference (lab) frame, and
t MRF ,
v is the moving reference frame velocity. All governing
equations then need to be rewritten with relative velocities.
Therefore, as mentioned before, a velocity field function can be defined and by defining a
motion based on this velocity field and setting the reference frame to moving reference frame,
the cylinder oscillation can be incorporated to the model. The cylinder oscillation has a sine
form, so the velocity field will be cosine.
t f A f
D
y
t f A
D
y
2 cos 2 2 sin
* *
( 6-2)
127
Figure 6-3 – Boundary conditions and the velocity field to model the cylinder oscillation
It should be noted that using such a technique does not apply the velocity field in the entire
domain. In fact, MRF conveys momentum onto the flow where there are solid wall boundaries
only and there will be no additional momentum introduced by MRF to other boundaries or
fluid/fluid interfaces. Simulations also confirm that MRF affects only the region near the
cylinder wall. Figure 6-4 shows the top view of the velocity vector for a simulation using MRF
for a section of the domain cutting the y axis. It is clear from the figure that MRF introduces the
velocity field associated with the cylinder oscillation only at the cylinder boundary and the
locations far from the cylinder do not get affected by this velocity field.
Figure 6-4 – Vector of velocity for flow past a moving cylinder using MRF
128
Tracking a point located at the cylinder wall illustrates that the modeled cylinder motion by
MRF is exactly the same as what it is expected. This is shown in Figure 6-5.
Figure 6-5 – Time trace of simulated motion by MRF for a case with an oscillation amplitude of 0.028 m
and oscillation period of 2.4 seconds
It should be noted that for the first 2.4 second only flow past a stationary cylinder is
simulated as an approach to initialize the model and there is no motion introduced to the system.
Morphing:
When morphing is activated, the software deforms the mesh with the specified morphing
velocity which is denoted by
g
v at the beginning of each time step. Then, as discussed before,
g
v will appear in all governing equations. Morphing requires mesh regeneration which can be
costly for high values of mesh count. More importantly, as the mesh is compressed and stretched,
in this problem near the cylinder, the mesh quality will decrease and this can have negative
effects on the boundary layer resolution.
To minimize the bad effects of morphing on the mesh near the wall, there are several
techniques and one way can be solving a Laplace equation for the morpher velocity and
introducing a diffusive term in it. This diffusive term can hold different forms and Kim (2014)
129
showed that the quadratic distance-based formula led to the minimum distortion. A similar
technique is used in Star-CCM+ as well. More information regarding mesh morphing can be
found in the Star-CCM+ documentation as well as Kim (2014) and Jasak & Tukovic (2004).
2) Mesh overset:
In mesh overset, one mesh is moving on top of another fixed mesh to simulate any mesh
motion. Overset is a faster option compared to morphing, however, according to CD-adapco, it
can also degrade the mesh quality near the cylinder.
In this research, the first two techniques are used to model the cylinder motion. First the
MRF approach is employed as it is faster and does not disturb the mesh quality. Later on,
morphing will be applied to study the differences between morphing and MRF.
6-2-2- Turbulence Modeling
Simulations for the stationary cylinder showed that DES with a Spalart-Allmaras closure
works well in predicting the statistical characteristics of the flow. It was decided to start with this
model for the oscillating cylinder as well. Later on, other DES models, i.e. k-ω SST and k-ε EB,
as well as LES models will be employed to compare the results.
6-2-3- Simulation Initiation and Inner Iteration
Like the simulations for flow past a stationary cylinder, for a moving cylinder, the simulation
is initiated by a steady state run with high iterations. For some of simulations a transient run for
the fixed cylinder is performed and then the cylinder oscillation is applied. The rest of
simulations are continued with the cylinder oscillation right after the steady run. The final results
are shown to be the same for both scenarios. The solution proceeds to next time step, like
stationary cylinder simulations, after a certain number of iteration. This number lies between 10
130
to 40 depending on the time step and temporal resolution but it is kept fixed for the entire
simulation.
6-2-4- Added Mass
One more important factor that should be included in modeling the flow past an oscillating
cylinder is the added mass. Added mass is the inertia force due to the displaced fluid and it is in
the direction of fluid acceleration.
a F
fluid displaced mass added
m
( 6-3)
Where a is the acceleration of cylinder motion. Considering the oscillation form for the
cylinder, a simple series of calculations shows that the coefficient of lift due to added mass can
be written as:
D
y
CL
mass added 2
*
3
2
( 6-4)
Where y is the cylinder motion and λ
*
is the normalized wavelength of oscillation.
In simulations, since the displaced fluid is not modeled, the drag and lift forces calculated
from pressure and shear stresses distributions on the cylinder wall do not include any added mass
effects. Therefore, the coefficient using these forces are noted as vortex forces coefficients. The
coefficient calculated from added mass is also denoted as the potential coefficient. The total
force which what is measured in the experiment will be the sum of vortex and potential forces.
p v t
CL CL CL
( 6-5)
Here CL
t
is the total lift, CL
v
the vortex lift, and CL
p
is the added mass effect. As shown
above, CL
p
is always in line with the cylinder oscillation. This decomposition of a hydrodynamic
force (in this case lift) into the potential and vortex components has been used by several
131
research studies including Lighthill (1986), Govardhan & Williamson (2000), and Morse &
Williamson (2009).
6-2-5- High Performance Computing
Solving CFD simulations like flow past a cylinder requires a lot of computational resources
due to high cell counts and small time steps. In this research, the USC supercomputer is used for
all simulations.
Simulation files are submitted to the supercomputer by PBS scripts. A PBS script is basically
a file that contains a series of commands and tells the system what computational resources and
for how long are required for a specific file. PBS files can be used to open the software and run
in the batch mode.
Figure 6-6 – An example of a PBS file to run a Star-CCM+ file in the batch mode
Figure 6-6 is an example of a PBS file where 300 CPUs (processes) with at least a 6.2 GB
memory for each process are requested for 24 hours.
132
For parallel processing, Star-CCM+ uses only InfiniBand (IB) computational nodes and it is
not compatible with MyriNet nodes any more. This is requested by “feature=octocore” in the
example above. It was previously shown by Bozorgnia (2012) that the number of CPUs per
computational node (ppn) can have some influence on computational time. This is because of
communication time between each CPU in one node and Bozorgnia (2012) showed that for his
particular CFD problem the optimum number for ppn in running Star-CCM+ on MyriNet nodes
is 2. However, the IB nodes are newer technology and ppn does not affect the computational
time that much, therefore it is decided not to limit the scheduler for choosing nodes with specific
ppn’s.
The other important factor in parallel processing is how different computational nodes are
communicating with each other. For this purpose a message passing interface (MPI) standard can
be used. An MPI is a library that allows the supercomputer to pass information between its
nodes. There are several MPI standards which are usually open source and the USC
supercomputer has some of them in its library.
A series of simple tests is performed to select the fastest MPI for this research study. For this
purpose, a 50 CPU job is requested and three cases are run: 1) specifying the MPI flag to
OpenMPI, 2) specifying the MPI flag to Platform, and 3) no flag for MPI. The Platform MPI by
IBM is the recommended MPI by CD-Adapco and it is confirmed that not using a flag for MPI
will call Platform. Results of these tests are listed in Table 6-1
Table 6-1 – Comparison of computational speed for different MPI flags
MPI Flag OpenMPI Platform No flag
Time to open the file (sec) 30 60 147
Number of iterations after 5 minutes 603 578 65
133
According to Table 6-1, OpenMPI has superiority over Platform and no flag (presumably
using Platform), therefore it is the choice for running Star-CCM+ simulations on the USC
supercomputer.
Finally it should be noted that, simulations depending on cell count, temporal resolution,
mesh motion, and turbulence model as well as number of CPUs requested take from 65 to 160
clock hours.
6-3- Results of Flow past an Oscillating Cylinder
Because of the importance of the new 2P
O
vortex shedding mode in the experimental studies
of Morse and Williamson (2009), 13 data points in the 2P
O
region are picked to simulate. These
points are shown in Figure 6-7.
Figure 6-7 – Black x shows the location of simulation points in the λ
*
vs A
*
map
The values of λ
*
and A
*
for data points are listed in Table 6-2
134
Table 6-2 – Data points used in simulations for wavelength and amplitude of oscillation
Morse and Williamson showed that in the lower part of the 2P
O
region, energy transferred to
the fluid switches to positive. This means that the oscillation characteristics for this part of the
map can lead to a freely oscillating vortex shedding. The energy transferred to the fluid
corresponds to C
L
sinϕ where ϕ is the phase lag between the lift and cylinder oscillation. More
importantly, according to the Morse and Williamson’s map, a portion of 2P
O
region with positive
excitation overlaps with the 2P region with negative excitation. If the oscillation properties lie in
this part of the map, the fluid excitation switches signs which correspond to switches between
shedding modes. This is illustrated for one point in Figure 6-8.
135
Figure 6-8 – Force and energy-transfer time traces for A
*
= 0.8 and λ
*
= 5.6 (2P - 2P
O
region) from Morse
and Williamson (2009)
Those 13 data points are selected in way to study switches between modes and positive
excitation to the system.
Figure 6-9 – Normalized fluid excitation map and selected data points for simulation
136
6-3-1- Calculating the Phase Lag
The phase lag of lift with respect to the cylinder oscillation (ϕ) is proven to be an important
feature of energy transfer. In fact, Morse and Williamson showed that if the phase lag lies
between 0 and 180 degrees, free oscillation is possible, otherwise fluid is not capable of exciting
the cylinder and free oscillation cannot take place.
t f F F t f D A y
Lift
2 sin 2 sin
1
*
( 6-6)
It was mentioned before that the vortex lift can be extracted by integrating the pressure and
shear stresses on the cylinder wall and the potential lift can be computed from the cylinder
oscillation. Thus, the total lift time trace will be available. The challenge, now, is how to find the
phase between the total lift and the cylinder oscillation (which has a phase lag of zero). For this
purpose, the Hilbert transform is used. Hilbert transform creates an analytical signal and provides
instantaneous properties of any signal meaning the magnitude and the phase, therefore, it can be
used as a tool to find the phase lag and later on lift magnitude in real time. Keeping in mind that
the total lift has one dominant frequency equal or very close to the frequency of oscillation (f),
the phase lag of lift can be calculated by subtracting the term 2πft from the instantaneous phase.
t f t C angle
L
2 Hilbert ( 6-7)
Since Hilbert transform gives the magnitude of the signal as well, it can also be used for
computing the term C
L
sin(ϕ):
t f t C angle t C Mag C
L L L
2 Hilbert sin * Hilbert sin ( 6-8)
It should be noted that all simulations are conducted on the part of A*-λ* map that has a
dominant frequency, so assuming that the lift will follow the same frequency is quite acceptable.
In fact, operating an fft function on the total lift shows that the signal has only one dominant
frequency equals to the oscillation frequency. This is illustrated in Figure 6-10 for two cases.
137
Figure 6-10 – Dominant frequency of lift for flow past an oscillating cylinder
To discuss different characteristics of flow, the vortex lift phase lag is also calculated by
Hilbert transform. The results for the selected data points are discussed in the following section.
It showed be noted that the phase lag of total lift is denoted by ϕ
t
and the phase lag of vortex lift
is denoted by ϕ
v
.
6-3-2- Results for Total and Vortex Lifts and Their Phase Lags
In this section results are shown and discussed. In each case, the potential lift coefficient is
also plotted to illustrate the differences in magnitude and phase for the total and lift coefficients.
For some of the data points more than one simulation either with different meshes, or turbulence
models, or temporal resolutions, or even mesh motion techniques has been performed. The
effects of each of these changes will be discussed later.
In all plots, time is normalized by T which is the period of oscillation for the cylinder.
U
D
T
*
( 6-9)
138
Figure 6-11 – Total lift coefficient and its phase lag for A* = 0.52 and λ* = 5.9
Figure 6-12 – Vortex lift coefficient and its phase lag for A* = 0.52 and λ* = 5.9
Figure 6-11 and Figure 6-12 show the results for A* = 0.52 and λ* = 5.9 which has the
lowest amplitude and lies at the bottom of the 2P
O
region. As it is clear, the vortex lift has a
lower magnitude and its phase is 180 degrees ahead compared to the potential lift. Having a low
vortex lift will lead to a total lift almost equal to the potential lift as it can be seen in Figure 6-11.
139
Figure 6-13 – Total lift coefficient and its phase lag for A* = 0.55 and λ* = 6.1
Figure 6-14 – Vortex lift coefficient and its phase lag for A* = 0.55 and λ* = 6.1
For A* = 0.55 and λ* = 6.1, the same trend can be seen; the vortex lift has a -180 degree lag
and its magnitude is still low, therefore, the total lift will be almost in line with the cylinder
oscillation. The only difference in this case compared to the first data point is that the total lift
phase lag almost diminishes and is very close to zero.
140
Figure 6-15 – Total lift coefficient and its phase lag for A* = 0.6 and λ* = 6.0
Figure 6-16 – Vortex lift coefficient and its phase lag for A* = 0.6 and λ* = 6.0
By increasing the oscillation amplitude the trend starts to change. As it is shown in
Figure 6-15 and Figure 6-16, while the vortex list still has a negative phase lag, this time the
phase lag of total lift is slightly negative and fluctuates about zero.
141
Figure 6-17 – Total lift coefficient and its phase lag for A* = 0.6 and λ* = 6.25
Figure 6-18 – Vortex lift coefficient and its phase lag for A* = 0.6 and λ* = 6.25
The point A* = 0.6 and λ* = 6.25 lies right at the border of the 2P
O
region. It carries almost
the same characteristics of the previous data point (A* = 0.6 and λ* = 6.0). The only difference
here is that because of higher wavelength the potential lift and consequently the vortex lift are
slightly lower compared to the previous point.
142
Figure 6-19 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 5.0
Figure 6-20 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 5.0
By further increase in the cylinder oscillation amplitude, the vortex lift changes
characteristics. As shown in Figure 6-20, the vortex lift magnitude is very close to the potential
lift. More importantly, its phase lag is slightly more positive. This causes a negative phase lag for
the total lift which is illustrated in Figure 6-19. It should be noted that the data point A* = 0.75
and λ* = 5.0 is at the border of 2P
O
region and it overlaps with the 2S region as well.
143
By keeping the oscillation amplitude constant and increasing the wavelength (λ*) phase lag
starts showing some changes. These changes are illustrated in the following figures.
Figure 6-21 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 5.5
Figure 6-22 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 5.5
144
Figure 6-23 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.0
Figure 6-24 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.0
Comparing Figure 6-20, Figure 6-22, Figure 6-24, Figure 6-26, and Figure 6-28 shows that
while the vortex lift phase lag remains almost the same with increasing λ*, the magnitude of
vortex lift drops compared to the potential lift.
145
Figure 6-25 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.25
Figure 6-26 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.25
146
Figure 6-27 – Total lift coefficient and its phase lag for A* = 0.75 and λ* = 6.5
Figure 6-28 – Vortex lift coefficient and its phase lag for A* = 0.75 and λ* = 6.5
The drop in vortex lift seems to be the main reason for changes in the total lift phase lag. As
it is clear from Figure 6-19, Figure 6-21, Figure 6-23, Figure 6-25, and Figure 6-27, the total lift
phase lag (ϕ
t
) moves toward zero and becomes more positive with increase in λ*. It also seems
like by increasing λ*, the magnitude of total lift increases as well.
147
Figure 6-29 – Total lift coefficient and its phase lag for A* = 0.78 and λ* = 6.0
Figure 6-30 – Vortex lift coefficient and its phase lag for A* = 0.78 and λ* = 6.0
148
Figure 6-31 – Total lift coefficient and its phase lag for A* = 0.8 and λ* = 5.0
Figure 6-32 – Vortex lift coefficient and its phase lag for A* = 0.8 and λ* = 5.0
Results for A* = 0.78 and λ* = 6.0, A* = 0.8 and λ* = 5.0, and A* = 0.8 and λ* = 5.6
(Figure 6-29 to Figure 6-34) show the same trend as discussed for “A* = 0.75” data points; the
vortex lift behaves the same way and it is relatively out of phase. Moreover, its magnitude is on
par with the potential lift. This pulls the total lift to a negative phase lag zone.
149
Figure 6-33 – Total lift coefficient and its phase lag for A* = 0.8 and λ* = 5.6
Figure 6-34 – Vortex lift coefficient and its phase lag for A* = 0.8 and λ* = 5.6
6-4- Discussion on Results
Total lift as explained before, has two components; vortex and potential. The potential lift
has a constant amplitude and it is always in line with the cylinder oscillation. Therefore, if there
are going to be any modulations in magnitude and switches in phase lag, the vortex lift should
take care of them. By a careful examination of the vortex lift, it seems that high oscillation
150
amplitude cases behave differently compared to low oscillation amplitude ones. In low cylinder
amplitudes, the vortex lift is out of phase from the very beginning and the phase lag is almost
close to -180 degrees throughout the time. This is not the case for high amplitude simulations,
where the vortex lift follows the cylinder motion for some time and there is no phase lag for most
of the data points. At some point, the vortex lift stops following the potential lift and shows an
abrupt drop in the magnitude. After this point the magnitude of vortex lift goes back up,
however, these changes cause a phase lag which this time is not quite out of phase and it
fluctuates between -90 and -180 degrees.
The other difference between vortex lift of low and high amplitudes is that as mentioned
before, the magnitude of vortex lift is much lower compared to the potential lift if the cylinder
oscillation amplitude is relatively low. On the other hand, for high oscillation amplitudes, the
magnitude of vortex lift is quite comparable with the potential lift and this is very distinct for low
wavelength cases.
6-5- Comparison between CFD and Test data
6-5-1- Total and Vortex Lift and Their Corresponding Phase Lags in Time
In the previous section, the CFD results were discussed and it was shown that for high
cylinder oscillation amplitudes, the total lift after some time follows a trend and exhibits
moderate magnitude modulations and moderate phase lag changes. Experimental studies by
Morse and Williamson (2009) shows that in the 2P
O
-2P overlap, the total lift demonstrates
magnitude modulations and phase lag switches in time. This is illustrated in Figure 6-35 and
Figure 6-36.
151
Figure 6-35 – Experimental data for total lift and its phase lag for A* = 0.8 and λ* = 5.6
Figure 6-36 – Experimental data for total lift and its phase lag for A* = 0.78 and λ* = 6.0 from Morse
and Williamson (2009a)
This is not the case in the 2P
O
-2S overlap, where Morse and Williamson (2009a) and (2009b)
showed that total lift does not have any significant magnitude modulations or phase lag switches.
In this region, however, the vortex lift displays these characteristics as shown in Figure 6-37.
152
Figure 6-37 – Experimental data for total and vortex lifts and their corresponding phase lags for A* =
0.8 and λ* = 5.0 Morse and Williamson (2009a)
Comparing the CFD results with experimental data for two points in the 2P
O
-2P overlap (A*
= 0.8 - λ* = 5.6 and A* = 0.78 - λ* = 6.0) shows that the total lift from CFD simulations matches
in overall magnitude. However, CFD results do not show any significant magnitude modulations
and phase lag switches observed in the experiment. The reason could lie in the behavior of
vortex lift as explained before. Phase lag switches in total lift seem to be affected by magnitude
modulations of vortex lift.
Figure 6-38 – Comparison between CFD simulations and test data of total lift coefficient for A* = 0.8
and λ* = 5.6
153
Figure 6-39 – Comparison between CFD simulations and test data of total lift phase lag for A* = 0.8 and
λ* = 5.6
This means that when vortex lift experiences any abrupt changes in magnitude, since it is out
of phase with the cylinder oscillation, the total lift will be affected and its magnitude will change.
For few first periods, the vortex lift follows the cylinder oscillation as shown in Figure 6-34 and
Figure 6-30. After some time, it seems that the vortex lift cannot keep up with the cylinder
oscillation and its magnitude suddenly drops.
Figure 6-40 – Comparison between CFD simulations and test data of total lift coefficient for A* = 0.78
and λ* = 6.0
Figure 6-41 – Comparison between CFD simulations and test data of total lift phase lag for A* = 0.78
and λ* = 6.0
Figure 6-42 and Figure 6-43 illustrate the vortex lift and its corresponding phase lag of A* =
0.78 and λ* = 6.0 for CFD and test. The magnitude of vortex lift coming from CFD is lower than
the test, but it shows modulations like the test. In the case of vortex lift phase lag both CFD and
154
test do not have any significant switches, however, while the test shows that the vortex lift is
perfectly out of phase, the phase lag from CFD is slightly toward zero.
Figure 6-42 – Comparison between CFD simulations and test data of vortex lift coefficient for A* = 0.78
and λ* = 6.0
Figure 6-43 – Comparison between CFD simulations and test data of vortex lift phase lag for A* = 0.78
and λ* = 6.0
For the point A* = 0. 8 - λ* = 5.0 (2S-2P
O
overlap), again, the predicted total lift coefficient
by CFD is slightly lower than the test but the general behavior of total lift resembles the test.
However, as it is clear in Figure 6-45, the vortex lift coming from CFD is much bigger than what
the experiment reported.
Figure 6-44 – Comparison between CFD simulations and test data of total lift coefficient for A* = 0.8
and λ* = 5.0
155
Figure 6-45 – Comparison between CFD simulations and test data of vortex lift coefficient for A* = 0.8
and λ* = 5.0
It is worth mentioning that examining the experimental results for A* = 0.8 - λ* = 5.0 and A*
= 0.78 - λ* = 6.0 shows that, even though both data points have a high oscillation amplitude, the
vortex lift of λ* = 5.0 is almost half of λ* = 6.0 which is a huge reduction. This is most likely
due to a shift in the shedding mode as λ* = 5.0 lies in the 2S region while λ* = 6.0 is in 2P. It
seems that CFD does not see such a shift in modes, because the vortex lift of λ* = 5.0 is higher
than λ* = 6.0. This increase in the vortex lift can also be seen when λ* is reduced for A* = 0.75
and A* = 0.6.
6-5-2- Fluid Excitation for Selected Data Points
As mentioned before C
L
sinϕ plays an important role in predicting the freely oscillating
cylinder phenomenon and if it is positive this can be an indicator that the flow excites the
cylinder and cylinder can move freely. By comparing the results of CFD simulations with the
experiment for the selected data points some interesting features can be concluded.
Using the Hilbert transform helps with find the instantaneous magnitude and phase of the
total lift coefficient. For each data point, the quantity of C
L
sinϕ is calculated and then the average
of it can be compared with the test data as shown in Figure 6-46.
156
Figure 6-46 – Normalized fluid excitation map and selected data points for simulation
Table 6-3 summarizes the result of C
L
sinϕ for the selected data points.
Table 6-3 – Mean value of C
L
sinϕ for selected data points
157
`
Figure 6-47 – Time trace of normalized fluid excitation for some of the selected data points
From Table 6-3, it can be seen that CFD results for most of the data points are pretty good if
they are compared with the value of C
L
sinϕ reported for the 2P region. In fact, the values listed in
Table 6-3 are very close to the test data for 2P in Figure 6-46. This is clear in Figure 6-48 where
values of CFD prediction are shown in orange boxes.
158
Figure 6-48 – CFD results of normalized fluid excitation (C
L
sinϕ) for the selected data points
Since CFD simulations do not predict any significant phase lag switches and the value of
C
L
sinϕ is in agreement with the test, it can be implied that CFD seems to overlook the 2P
O
region. Even values of points A* = 0.8 - λ* = 5.0 and A* = 0.75 - λ* = 5.0, which according to
the experiment lie in the 2S-2P
O
region, seem to be an extension to the 2P region. Therefore,
CFD seems to be able to capture the 2P region, as shown for the point A* = 0.75 - λ* = 6.5 and
somewhat for A* = 0.52 - λ* = 5.9, A* = 0.6 - λ* = 6.25, A* = 0.75 - λ* = 6.25 and its results are
in agreement with the test. However, CFD predictions do not display any changes in lift behavior
and consequently the shedding mode unlike what the experiment suggests.
6-6- Comparing CFD Results for Different Data Points
Results of CFD simulations for 13 data points in the A*-λ* map can be further used to
compare and discuss the effects of different parameters on lift. It should be noted that for some
159
of the selected data points more than one simulation has been conducted (in some cases 17
simulations with different parameters). Having several simulations for each data point will help
to compare different modeling parameters on the final result which will be presented in
following sections.
6-6-1- Changes in Total Lift for a Cut of A* = 0.75
For A* = 0.75, 5 different data points (λ* = 5, 5.5, 6.0, 6.25, 6.5) are conducted and
comparison of the total lift for these points are shown in Figure 6-49.
Figure 6-49 – Total lift for a cut of A* = 0.75
As the figure suggests, the magnitude of lift drops with increasing λ*. Previously, it was
discussed that for a constant cylinder amplitude, increase in λ*, decreases the vortex lift. Because
λ* appears in the potential lift formulation with a power of -2, the comparison that Figure 6-49
demonstrates is meaningful.
160
6-6-2- Changes in Total Lift for a Cut of λ* = 6.0
There are 4 data points with λ* = 6.0 (A* = 0.6, 0.7, 0.75, and 0.78). The total lift of three of
these points are illustrated in one plot in Figure 6-50.
Figure 6-50 – Total lift for a cut of λ* = 6.0
According to Figure 6-50, as the oscillation amplitude goes up, the total lift experience an
increase as well. The same trend can be drawn for a cut of λ* as well: The potential lift changes
linearly with the oscillation amplitude and the vortex lift was previously shown to increase from
low to high oscillation amplitude.
In the following sections, the effects of numerical parameters including mesh discretization,
turbulence modeling, mesh motion, wall region resolution, and temporal resolution on total lift
will be presented.
6-6-3- Sensitivity of Total Lift to Turbulence Modeling
Detached Eddy Simulation (DES) with the Standard Spalart-Allmaras closure (DES S-A)
showed a great performance for flow past a stationary cylinder so it was the first choice for the
161
moving cylinder case as well. Later on it was decided to use different turbulence models to test
each model capabilities for flow past a cylinder.
One of the characteristics of any DES model is the use of wall distance functions. Basically
any DES model switches between RANS and LES based on a tuned operator which is a function
of the distance to the wall. Using a wall distance function implies that the turbulence model,
regardless of the flow regime, will always resolve the near wall region with a RANS closure.
This can be a great feature especially if a laminar sublayer exists where LES resolutions either
struggle or are very expensive. However, if the regime consists of separations, like flow past a
cylinder, using a RANS closure comes with all deficiencies of RANS models that discussed
before.
Other DES models (the SST k-ω and the EB k-ε) as well as LES models (WALE,
Smagorinsky, and Dynamic Smagorinsky) are tested for some of the data points and comparison
between different models on the same mesh is presented in this section.
162
Figure 6-51 – Total lift for A* = 0.78 and λ* = 6.0 resulted from DES SA, DES EB k-ε, and LES WALE
models on different mesh configurations
Simulations on 3 different mesh configurations with different turbulence models for one data
point shows that, for the same mesh configuration, none of the turbulence models have any
superiority over others in terms of overall prediction of total lift. LES WALE seems to be able to
capture more details and magnitude modulations compared to DES models but the frequency,
phase lag, and magnitude for all models are very similar. It should be noted that it is predictable
163
that LES captures more details of the flow as it does not use any time averaging for flow
quantities.
This comparison implies that for the moving cylinder, DES is still a reliable technique, when
all costs of LES are considered.
6-6-4- Sensitivity of Total Lift to Spatial Resolution
To study the changes due to different mesh configurations, a series of simulations with the
LES WALE is performed for A* = 0.78 and λ* = 6.0. The mesh details are presented in the
following table.
Table 6-4 – Mesh details for flow past a moving cylinder with A* = 0.78 and λ* = 6.0
LES_05 is the baseline for this analysis as it is the mesh configuration selected as the best for
flow past a stationary cylinder. LES_08 is basically the same mesh and the only difference is in
the first cell size which is reduced by half to see the effects of wall y
+
value. LES_09 and
LES_10 are based on LES_05 and LES_08 respectively with a higher z discretization. The
spanwise discretization is shown to have some impacts on spatial resolution for LES simulation
by some researchers in the literature and in LES_13 and LES_14 it is tried to mesh the domain
with a high z discretization value.
164
According to Figure 6-52, having a smaller 1
st
cell size will help to capture more details of
the flow. However, for the majority of the time history both meshes predictions are very close.
This seems to show that having wall y
+
limited to 2 (LES_05) and limited to 1 (LES_08) does
not make a drastic change in resolving the boundary layer.
Figure 6-52 – Effects of 1
st
cell thickness on total lift for A* = 0.78 and λ* = 6.0 resulted by LES
Figure 6-53 illustrates the impact of a 25% increase in the spanwise discretization. As it can
be seen, having a higher resolution in z will provide more details, however, especially when
values of wall y
+
are limited to 1 (LES_08 and LES_10), the high resolution mesh does not result
in a significant change in the trend of total lift. It should be noted that a 25% increase in z
discretization means that the cell count is up by 25% which consequently means more expensive
simulations.
165
Figure 6-53 – Effects of 25% increase in z discretization on total lift for A* = 0.78 and λ* = 6.0 resulted
by LES
Just to compare the coarsest and most refined meshes, results of LES_05 (2.6 million cells)
and LES_14 (4 million cells) are shown in Figure 6-54. It should be noted that that the
computational time of LES_14 is more than twice higher than LES_05.
166
Figure 6-54 – Total lift for A* = 0.78 and λ* = 6.0 resulted by 2.6 million and 4 million cells
Figure 6-55 compares all these mesh configurations for A* = 0.78 and λ* = 6.0. As discussed
earlier, having higher resolutions in r (more specifically the 1
st
cell thickness) and z can provide
more details of lift changes in time in LES simulations. However, all meshes display the same
trend and final results are quite comparable from the coarsest to the most refined mesh.
Figure 6-55 – Total lift for A* = 0.78 and λ* = 6.0 resulted by different mesh configurations
167
6-6-5- Sensitivity of Total Lift to Temporal Resolution
All simulations are conducted in a way that the convective CFL number is limited to 1. For
A* = 0.78 and λ* = 6.0, an extra LES WALE simulation is done by decreasing the time step. For
the original simulation CFL number is bounded to 0.7 meaning that the new simulation has a
maximum 0.07 for CFL which is pretty small.
Figure 6-56 – Temporal resolution effects on the total lift for A* = 0.78 and λ* = 6.0
Plots in Figure 6-56 would suggest that using a much lower time step (10 times smaller) does
not necessarily lead to a drastic change in total lift. Even though the model is 10 times slower for
the new simulation, the general trend of lift seems to stay the same. Of course, the simulation is
capable of catching more fluctuations and details of the flow, but, the expense of simulation may
not justify the final result.
6-6-6- Sensitivity of Total Lift to Mesh Motion Technique
For some cases, in order to see the effects of mesh motion techniques, the same simulation is
performed by Moving Reference Frame (MRF) and morphing. It should be noted that morphing
168
requires more computational time as the model regenerates the mesh at the beginning of each
time step following the motion of cylinder.
Moreover, morphing deforms the mesh near the cylinder wall so the quality of mesh changes
from time to time. However, since the time step is very small and the mesh motion velocity is not
very high, these distortions do not affect the near wall spatial resolution (more specifically wall
y
+
) much. In fact, analyses show that for the same simulation, the value of wall y
+
is not very
different when morphing is used instead of MRF. This is also true for the maximum CFL number
for the entire domain.
Figure 6-57 – Total lift resulted MRF vs morphing for A* = 0.7 and λ* = 6.0 with LES_14_WALE
Morphing can be applied in the model based on a displacement or a velocity field (
t A y sin vs t A y cos ). As both Figure 6-57 and Figure 6-58 demonstrate either
technique results in relatively the same total lift for A* = 0.7 and λ* = 6.0 with an LES
simulation on a very refined mesh. Morphing based on a displacement field seems to predict
more details and fluctuations (especially compared to MRF), but considering the higher
computational costs of morphing, MRF would be still a good choice for mimicking the cylinder
oscillation.
169
Figure 6-58 – Effects of morphing field type on total lift for A* = 0.7 and λ* = 6.0 with LES_14_WALE
6-6-7- Effects of Convection Terms Resolution Schemes on Total Lift
Using a central differencing scheme (CDS) for calculating the velocities at the cell faces, as
discussed in Chapter 2, introduces no extra numerical (artificial) viscosity into the flow.
However, it is a sensitive scheme and can cause numerical instabilities. On the other hand, the
bounded central differencing is the recommended scheme for convective terms because it
switches to an upwind scheme if quantities pass a bounded limit, therefore it is stable. It should
be noted that upwind schemes are susceptible to adding viscosity to the system and this can have
impacts for sensitive problems like flow past a cylinder.
The majority of simulations are performed by the bounded central differencing, however, a
few of them are redone using the central differencing. Here, results of both schemes for one data
point is presented to investigate the effects of central differencing on total lift.
170
Figure 6-59 – Effects of convective scheme on total lift for A* = 0.75 and λ* = 5.5 by LES WALE and
LES Smagorinsky with 4 million cells
It should be noted that the time step is set lower for the simulation with the central
differencing to have the CFL number limited to 1. Figure 6-59 shows that central differencing
seems to resolve more details of the total lift. However, even though the time step is lower and
CFL is still less than one for this case, the model becomes unstable at t/T = 43.28, which
indicates the sensitivity of the scheme.
6-6-8- The Last Comparison
At the end of this section, it would be noteworthy to compare the very first simulation to the
very last one. For this purpose, total lift for A* = 0.8 and λ* = 5.6 resulted from the coarsest
mesh with MRF using DES S-A and the bounded central differencing scheme are plotted against
the most refined mesh with morphing using LES WALE and the central differencing scheme.
Details of each mesh are in Table 6-5.
171
Table 6-5 – Details of the coarsest and most refined mesh configurations used for A* = 0. 8 and λ* = 5.6
It is worth noting that a simple series of calculations shows that the computational cost of the
second simulation for each time step is roughly 4.35 times the first simulation with the coarsest
mesh. Considering the temporal resolution, this would mean that to run both simulations for the
same real time (seconds), the simulation using morphing and LES requires 8.7 times more
computational resources.
Figure 6-60 – Total lift for A* = 0.8 and λ* = 5.6 resulted from a DES simulation with 2.6 million cells
and MRF vs an LES simulation with 4 million cells and morphing
It is shown in Figure 6-60 that, considering the cost difference between the simulations, the
DES simulation seems to do a very good job compared to a more refined LES simulation. This
172
comparison would indicate that, as long as the temporal and spatial resolutions are studied well,
MRF can be used to model the cylinder oscillation in a simulation with a relatively coarse mesh
using the simplest form of DES.
6-6-9- Discussion on Numerical Dissipation
Numerical dissipation is introduced to the model by many parameters including but not
limited to finite volume (any spatial) discretization, iterative solution, and approximations and
assumptions of the turbulence model. Numerical dissipation would become more important for
CFD challenges are more sensitive to the viscosity because numerical dissipation shows its
effects with added (artificial) viscosity. In fact, there are several studies that focus on the
inherent numerical dissipation of LES simulations, its magnitude compared to the SGS
dissipation, and techniques for understanding and mitigating its effects. One example is the
worked done by Castiglioni and Domaradzki (2015) where they showed that, for a realistic CFD
problem with a laminar sublayer as well as transitions and separations, the generated viscosity by
the numerical dissipation can be an order of magnitude higher than the viscosity of the sub-grid
scale model.
In this study, even though there has been no direct attempt to calculate and minimize the
effects of numerical dissipation, but using a structured mesh, having low wall y+ values, and a
proper temporal resolution would potentially help with reducing the numerical viscosity effects
on the model. In fact, keeping in mind that the lift value is almost lowest for the studied range of
Reynolds number (Figure 3-1), testing the model configuration, i.e. the mesh, time step, and
turbulence model, for the stationary case would work as an indicator of the numerical dissipation
presents in the system.
173
Moreover, the sensitivity analysis of simulations to different parameters can also been used
as a measure to check the numerical dissipation effects on the lift prediction. As discussed in this
section, the lift coefficient seems to be not so sensitive to the numerical parameters refinements.
This would mean that the inherent dissipation coming from numerics, although changing and
most likely decreasing from simulations with a coarse mesh using DES to simulations with a
finer mesh using LES and CDS, does not seem to be effective on the final answer.
6-7- Concluding Remarks
This chapter focused on simulating flow past an oscillating cylinder with a 3D domain. Using
the model parameters from the stationary cylinder simulations including the mesh configuration,
temporal resolution, and turbulence modeling helped as a starting point for the oscillating case.
To investigate the capabilities of CFD simulations, experimental data of Morse and Williamson
(2009) is used and 13 data points in their amplitude vs frequency (A* - λ*) map were selected.
Decomposing the total lift into the vortex and potential components was shown to be useful to
examine the lift behavior for different data points. Moreover, the instantaneous properties of lift
time history was computed by the Hilbert transform for each simulation case.
CFD simulations showed that in general an increase in the cylinder oscillation amplitude
leads to higher total lift values. CFD also illustrated that higher oscillation frequencies (or
lowering the wavelength λ*) generally causes the lift to grow up. Comparing the CFD results
with the available test data for three points showed that the magnitude of lift predicted by CFD is
in relatively a good agreement, however, unlike the test data for the overlapped 2P
O
region, CFD
results do not display any significant amplitude modulations in lift. More importantly, CFD
results for the 2P
O
region do not predict any notable phase lag switches for the lift time history.
174
The reason for these predictions would be that the vortex lift is unable to follow the cylinder
oscillation and the potential lift after some time and it drops suddenly and seems to reach to a
stable condition. This observation seems to be more noticeable for high amplitude and high
frequency cases.
On the other hand, in terms of reproducing the values of normalized fluid excitation, CFD
seems to do a very good job. In fact, CFD results match with the experimental data for most
cases. It was seen that CFD appears to be able to predict the 2P region with quite an acceptable
accuracy. The absence of switches in the phase lag seems to prevent the CFD results from giving
any predictions for the 2P
O
region. Finally, after examining the CFD results, it seems that CFD
somehow overlooks the 2P
O
region suggested by the experiment.
Investigating the sensitivity of total lift to different CFD simulation properties including the
mesh configurations and the spatial resolution, the turbulence modeling, the temporal resolution,
the mesh motion technique, and the convective terms resolution scheme indicated that adding
more refinements to the solution by implementing more accurate approaches leads to more
details in lift predictions, but the general trend, i.e. the magnitude and the phase lag, does not
change much. Therefore, the final result of these refinements seems to be unable to justify the
additional complexity and computational cost.
175
Conclusions
This research has been focusing on the numerical simulation of flow past a cylinder for Re =
4,000. For this purpose, Star-CCM+, a CFD solver developed by CD-adapco, is used. The study
starts by simulating flow past a stationary cylinder with RANS turbulence model in a 2D
domain. This step helps to figure out a proper pattern and size for the domain discretization and
choosing an appropriate time step for the problem of flow past a cylinder with a prescribed
motion. A rectangular domain meshed with polygonal cells and prism layers at the cylinder wall
is used for 2D simulations.
The results indicate that, regarding the statistical characteristics of the flow (the mean drag,
the RMS of lift, and the Strouhal number), some simulations, depending on the turbulence
model, are able to reproduce the drag coefficient and Strouhal number with an acceptable
accuracy. However, the lift predicted by all cases is much higher than the experiment. The other
discrepancy of 2D simulations is the absence of amplitude modulations in the time history of
forces unlike what experiments show. These two differences are the result of the nature of RANS
models and 2D simulation in attempting to resolve a 3D flow pattern with separations,
reattachments, and transitions.
The next part of the study focuses on the 2D simulation of flow past a moving cylinder. The
mesh configurations used in the stationary cylinder case are subjected to a cross-flow sinusoidal
oscillation. The results for the lift show that, again, depending on the turbulence model, time
traces do not demonstrate any amplitude modulations in some cases and unlike the experiment,
the phase lag between lift and the cylinder motion seems to be constant in time. Interestingly,
176
those simulations with amplitude modulations do not lead to a distinguishable vortex shedding
pattern. In none of the cases, the vortex shedding pattern looks like what the experiment reveals.
The study continues with investigation of 3D simulations for a stationary cylinder. 3D
simulations should have some improvements over 2D ones including amplitude modulations,
disturbance of 2D flow structures in the 3
rd
dimension, and drop in statistical values of lift and
drag. However, by using the same planar mesh configurations extruded in the spanwise direction,
results show that some turbulence models do not lead to any improvements, i.e. three
dimensionality, and some do.
Studying different turbulence models in more details indicates that the models without three
dimensionality employ damping functions for resolving the near cylinder region. These damping
functions seem to eliminate any three dimensionality and immobilize the 3
rd
dimension.
Introducing anisotropy or nonlinearity seems to help with counterbalancing the effects of
damping functions. On the other hand, the turbulence models that do not require any damping
functions for the wall region trigger three dimensionality. The damping effects are detectable in
2D simulations as well and by analyzing the 2D results it seems that, in general, models with
damping functions demonstrate lower lift and in the moving case and there are no amplitude
modulations in their time histories.
After examining the RANS models characteristics for simulating flow structures due to
vortex generation behind the cylinder, scale resolving turbulence models, i.e. DES and LES, are
employed. A new o-shaped domain with a structured mesh is employed for this section and by
simulating flow past a stationary cylinder a proper set-up for the mesh details, the turbulence
model, and temporal resolution is derived. Results show that if the spatial and temporal
resolutions are adjusted with appropriate controls, the statistical properties of flow can be
177
predicted correctly by a relatively coarse mesh and higher time step using a turbulence model
with less fidelity compared to other similar numerical studies.
The final step of this study centers on 3D simulating of flow past an oscillating cylinder.
CFD results show that, for a specific range of the oscillation amplitude and frequency, the lift
coefficient goes up by increasing the amplitude and frequency. Unlike the experimental study
conducted by Morse and Williamson (2009) where it is shown that the vortex shedding pattern
can alter in time for a constant amplitude and frequency, the current CFD simulations do not
predict any switches in shedding modes in time. These switches are supposed to be manifested
by noticeable lift magnitude modulations along with significant lift force phase lag changes
compared to the cylinder oscillation.
On the other hand, the normalized fluid excitation produced by CFD has a noteworthy
agreement with the available test data. In fact, comparing the test data with the CFD results
indicates that CFD seems to be capable of predicting the vortex shedding mode and the energy
transfer to the fluid. However, it seems that CFD does not exhibit any formation of new shedding
modes for the studied range of amplitude and frequency.
Conducting more accurate yet more computationally expensive simulations shows that
having more refinements in terms of spatial and temporal resolutions as well as the choice of
turbulence model and mesh motion technique seems to add more details in predictions for the lift
coefficient. However, the overall pattern of results does not change from simpler to more
sophisticated approaches and using more computationally demanding methods does not seem to
be justified for this type of CFD challenge. The final conclusion is that CFD seems to be a
reliable tool for predicting the details of flow regime for flow past a stationary and moving
cylinder for Re = 4,000.
178
References
Abe, K., Kondoh, T. & Nagano, Y., 1995. A new turbulence model for predicting fluid flow and
heat transfer in separating and reattaching flows—II. Thermal field calculations.
International Journal of Heat and Mass Transfer, 38(8), pp.1467–1481.
Apte, S., Introduction to Turbulence Modeling; Lecture Notes of ME569 - Oregon State
University,
Bearman, P.W., 2009. Understanding and predicting vortex-induced vibrations. Journal of Fluid
Mechanics, 634, p.1.
Billard, F. & Laurence, D., 2012. A robust k − ε − v 2 / k elliptic blending turbulence model with
improved predictions in near-wall , separated and buoyant flows. International Journal of
Heat and Fluid Flow, 33(1), pp.45–58.
Bishop, R.E.D. & Hassan, A.Y., 1964. The lift and drag forces on a circular cylinder oscillating
in a flowing fluid. Proceedings of Royal Society of London.
Blevins, R.D., 1990. Flow-induced vibration 2nd ed., Van Nostrand Reinhold.
Bozorgnia, M., 2012. Computational Fluid Dynamic Analysis of Highway Bridge
Superstructures Exposed to Hurricane Waves (PhD Thesis). University of Southern
California.
Bozorgnia, M., Eftekharian, A. & Lee, J.-J., 2014. CFD Modeling of a Solitary Wave
Overtopping Breakwater of Varying Submergence. In Coastal Engineering Proceedings. p.
7.
Breuer, M., 2000. A challenging test case for large eddy simulation: High Reynolds number
circular cylinder flow. International Journal of Heat and Fluid Flow, 21(5), pp.648–654.
Breuer, M., 1998. Large eddy simulation of the subcritical flow past a circular cylinder:
Numerical and modeling aspects. International Journal for Numerical Methods in Fluids,
28(9), pp.1281–1302.
Brika, D. & Laneville, A., 1993. Vortex-induced vibrations of a long flexible circular cylinder.
Journal of Fluid Mechanics, 250, pp.481–508.
Castiglioni, G. & Domaradzki, J. a., 2015. A numerical dissipation rate and viscosity in flow
simulations with realistic geometry using low-order compressible Navier-Stokes solvers.
179
Computers and Fluids, 119, pp.37–46. Available at:
http://dx.doi.org/10.1016/j.compfluid.2015.07.004.
Causon, M.D., Mingham, C.G. & Qian, L., 2011. Introductory Finite Volume Methods for PDEs,
bookboon.com.
CD-adapco, 2016. STAR-CCM+ Documentation.
Chaplin, J.R. et al., 2005. Blind predictions of laboratory measurements of vortex-induced
vibrations of a tension riser. Journal of Fluids and Structures, 21(1 SPEC. ISS.), pp.25–40.
Dong, S. & Karniadakis, G.E., 2005. DNS of flow past a stationary and oscillating cylinder at Re
= 10000. Journal of Fluids and Structures, 20(4 SPEC. ISS.), pp.519–531.
Durbin, P.A., 1991. Near-wall turbulence closure modeling without “damping functions.”
Theoretical and Computational Fluid Dynamics, 3(1), pp.1–13.
Feng, C.C., 1968. The measurement of vortex idunced effects in flow past stationary and
oscillating circular and D-section cylinders (MSc Thesis). The University of British
Columbia, Canada.
Govardhan, R. & Williamson, C.H.K., 2006. Defining the “modified Griffin plot” in vortex-
induced vibration: revealing the effect of Reynolds number using controlled damping.
Journal of Fluid Mechanics, 561, p.147.
Govardhan, R. & Williamson, C.H.K., 2000. Modes of vortex formation and frequency response
of a freely vibrating cylinder. Jourmal of Fluid Mechanics, 420, pp.85–130.
Hover, F.S., Davis, J.T. & Triantafyllou, M.S., 2004. Three-dimensionality of mode transition in
vortex-induced vibrations of a circular cylinder. European Journal of Mechanics - B/Fluids,
23(1), pp.29–40.
Humphreys, J.S., 1960. On a circular cylinder in a steady wind at transition Reynolds numbers.
Journal of Fluid Mechanics, 9(04), pp.603–612.
Iwan, W.D. & Blevins, R.D., 1974. A Model for Vortex Induced Oscillation of Structures.
Journal of Apllied Mechanics.
Jasak, H. & Tukovic, Z., 2004. Automatic mesh motion for the unstructured Finite Volume
Method. Transactions of Famena, 30(2), pp.1–20.
Kara, M.C., Stoesser, T. & McSherry, R., 2015. Calculation of fluid–structure interaction:
methods, refinements, applications. In Proceedings of the Institution of Civil Engineers.
180
Khalak, A. & Williamson, C.H.., 1999. Motions, Forces and Mode Transitions in Vortex-
Induced Vibrations At Low Mass-Damping. Journal of Fluids and Structures, 13(7-8),
pp.813–851.
Khalak, A. & Williamson, C.H.K., 1997. Fluid forces and dynamics of a hydroelastic structure
with very low mass and damping. Journal of Fluids and Structures, 11, pp.973–982.
Kim, S., 2014. LARGE EDDY SIMULATION OF FORCES AND WAKE MODES OF AN
OSCILLATING CYLINDER. University of Southampton.
Kim, S., Wilson, P. a. & Chen, Z.M., 2014. Numerical simulation of force and wake mode of an
oscillating cylinder. Journal of Fluids and Structures, 44, pp.216–225.
Lienhard, J.H., 1966. Synopsis of lift, drag, and vortex frequency data for rigid circular
cylinders,
Lighthill, J., 1986. Fundamentals concerning wave loading on offshore structures. Journal of
Fluid Mechanics, 173, pp.667–681.
Lysenko, D. a., Ertesvag, I.S. & Rian, K.E., 2012. Large-eddy simulation of the flow over a
circular cylinder at Reynolds number 3900. Flow Turbulence Combust, pp.1–29.
Meneghini, J.R. & Bearman, P.W., 1995. Numerical Simulation High Amplitude Oscillatory
Flow About a Circular Cylinder. Journal of Fluids and Structures, 9, pp.435–455.
Morse, T.L. & Williamson, C.H.K., 2009. Fluid forcing, wake modes, and transitions for a
cylinder undergoing controlled oscillations. Journal of Fluids and Structures, 25(4),
pp.697–712. Available at: http://dx.doi.org/10.1016/j.jfluidstructs.2008.12.003.
Morse, T.L. & Williamson, C.H.K., 2009. Prediction of vortex-induced vibration response by
employing controlled motion. Journal of Fluid Mechanics, 634, p.5.
Navrose & Mittal, S., 2013. Free vibrations of a cylinder: 3-D computations at Re = 1000.
Journal of Fluids and Structures.
Newman, D.J. & Karniadakis, G.E., 1997. A direct numerical simulation study of flow past a
freely vibrating cable. Journal of Fluid Mechanics, 344(1), pp.95–136.
Norberg, C., 2003. Fluctuating lift on a circular cylinder: Review and new measurements.
Journal of Fluids and Structures, 17(1), pp.57–96.
Ouvrard, H. et al., 2008. Variational Multiscale LES and Hybrid RANS / LES Parallel
Simulation of Complex Unsteady Flows. In VECPAR’08 — 8th International Meeting High
Performance Computing for Computational Science. pp. 465–478.
181
Pereira, F., Vaz, G. & Eca, L., 2015. Flow Past a Circular Cylinder: A Comparison Between
RANS and Hybrid Turbulence Models for a Low Reynolds Number. In Proceedings of
ASME 34th International Conference on Ocean, Offshore and Arctic Engineering.
Prsic, M. & Myrhaug, D., 2012. LARGE EDDY SIMULATIONS OF THREE-DIMENSIONAL
FLOW AROUND A PIPELINE IN A UNIFORM CURRENT. In Proceedings of the ASME
2012 31st International Conference on Ocean, Offshore and Arctic Engineering. pp. 1–10.
Salvetti, M. V et al., 2010. Quality and Reliability of Large-Eddy Simulations II, Springer
Science & Business Media.
Singh, S.P. & Mittal, S., 2005. Flow past a cylinder: Shear layer stability and drag crisis.
Internation Journal for Numerical Methods in Fluids, 47, pp.75–98.
Singh, S.P. & Mittal, S., 2005. Vortex-induced oscillations at low reynolds numbers: Hysteresis
and vortex-shedding modes. Journal of Fluids and Structures, 20(8), pp.1085–1104.
Smagorinsky, J., 1963. General Circulation Experiments With the Primitive Equations. Monthly
Weather Review, 91(3), pp.99–164. Available at:
http://journals.ametsoc.org/doi/abs/10.1175/1520-
0493(1963)091<0099:GCEWTP>2.3.CO;2.
Sodja, J., 2007. Turbulence models in CFD,
Spalart, P.., 2000. Strategies for turbulence modelling and simulations. International Journal of
Heat and Fluid Flow, 21(3), pp.252–263.
Spalart, P.R. et al., 1997. Comments on the feasibility of LES for wings, and on a hybrid
RANS/LES approach. In 1st AFOSR Int. Conf. on DNS/LES Aug.
Thingbø, S.S., 2013. Simulation of viscous Flow around a circular Cylinder with STAR-CCM +
(MSc Thesis). Norwegian University of Science and Technology.
Tuan, L. & Temarel, P., 2014. Numerical Simulation of an Oscillating Cylinder in Cross-Flow
At a reynolds Number of 10,000: Forced and Free Oscillations. In Proceedings of the ASME
2014 33rd International Conference on Ocean, Offshore and Arctic Engineering. pp. 1–10.
Udaykumar, H.S. et al., 2001. A Sharp Interface Cartesian Grid Method for Simulating Flows
with Complex Moving Boundaries. Journal of Computational Physics, 174(1), pp.345–380.
Versteeg, H.K. & Malalasekera, W., 2007. An Introduction to Computational Fluid Dynamics:
The Finite Volume Method, Pearson Education Limited.
Wieselsberger, C., 1922. New Data on the Laws of Fluid Resistance - Technical Notes,
Wilcox, D.C., 2006. Turbulence Modeling for CFD, DCW Industries.
182
Willden, R.H.J., 2006. Numerical Simulation of the Multiple Branch Transverse Response. In
2006 ASME Pressure Vessels and Piping Division Conference.
Williamson, C.H.. & Roshko, A., 1988. Vortex Formation in the Wake of an Oscillating
Cylinder. Journal of Fluids and Structures.
Williamson, C.H.K., 1988. The existence of two stages in the transition to three ‐ dimensionality
of a cylinder wake. Physics of Fluids, 31(11), pp.3165–3168.
Williamson, C.H.K., 1996. Vortex dynamics in the cylinder wake. Annu. Rev. Fluid Mech., 28,
pp.477–539.
Wornom, S. et al., 2011. Variational multiscale large-eddy simulations of the flow past a circular
cylinder: Reynolds number effects. Computers & Fluids, 47(1), pp.44–50.
Wu, W., 2011. Two-Dimensional RANS Simulation of Flow Induced Motion of Circular Cylinder
With Passive Turbulence Control (PhD Thesis). The University of Michigan.
Yang, J., Preidikman, S. & Balaras, E., 2008. A strongly coupled, embedded-boundary method
for fluid-structure interactions of elastically mounted rigid bodies. Journal of Fluids and
Structures, 24(2), pp.167–182.
Young, M. & Ooi, A., 2010. Comparative Assessment of LES and URANS for Flow Over a
Cylinder at a Reynolds Number of 3900. In 16th Australasian Fluid Mechanics Conference
(AFMC). pp. 1063–1070.
Zhuang, F. & Lee, J., 1996. A viscous rotational model for wave overtopping over marine
structure. In Coastal Engineering Proceedings. pp. 2178–2191.
Abstract (if available)
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
A Boltzmann model for tracking aerosols in turbulent flows
PDF
Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves
PDF
Numerical and experimental study on dynamics of unsteady pipe flow involving backflow prevention assemblies
PDF
Numerical modeling of separated flows at moderate Reynolds numbers appropriate for turbine blades and unmanned aero vehicles
PDF
RANS simulations for flow-control study synthetic jet cavity on NACA0012 and NACA65(1)412 airfoils.
PDF
Kinetic studies of collisionless mesothermal plasma flow dynamics
PDF
Flow and thermal transport at porous interfaces
PDF
Reconstruction and estimation of shear flows using physics-based and data-driven models
Asset Metadata
Creator
Eftekharian, Amirhossein
(author)
Core Title
Numerical study of flow characteristics of controlled vortex induced vibrations in cylinders
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Civil Engineering
Publication Date
08/08/2018
Defense Date
05/08/2018
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
computational fluid dynamics (CFD),controlled VIV,DES,flow past a cylinder,fluid excitation,fluid structure interaction,LES,lift and drag coefficients,numerical simulations,OAI-PMH Harvest,RANS,turbulence modeling,vortex induced vibrations (VIV),vortex shedding,vortex shedding modes
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Lee, Jiin-Jen (
committee chair
), Domaradzki, Julian (
committee member
), Masri, Sami (
committee member
), Wellford, Carter L. (
committee member
)
Creator Email
amirhose@usc.edu,eftekharian.amirhossein@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-64694
Unique identifier
UC11671354
Identifier
etd-Eftekharia-6706.pdf (filename),usctheses-c89-64694 (legacy record id)
Legacy Identifier
etd-Eftekharia-6706.pdf
Dmrecord
64694
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Eftekharian, Amirhossein
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
Tags
computational fluid dynamics (CFD)
controlled VIV
DES
flow past a cylinder
fluid excitation
fluid structure interaction
lift and drag coefficients
numerical simulations
RANS
turbulence modeling
vortex induced vibrations (VIV)
vortex shedding
vortex shedding modes