ABSTRACT
Title of dissertation: AN EMBEDDED-BOUNDARY FORMULATION FOR
LARGE-EDDY SIMULATION OF TURBULENT FLOWS
INTERACTING WITH MOVING BOUNDARIES
Jianming Yang, Doctor of Philosophy, 2005
Directed by: Assistant Professor Elias Balaras
Department of Mechanical Engineering
A non-boundary-conforming formulation for simulating transitional and tur-
bulent flows with complex geometries and dynamically moving boundaries on fixed
orthogonal grids is developed. The underlying finite-difference solver for the filtered
incompressible Navier-Stokes equations in both Cartesian and cylindrical coordi-
nates is based on a second-order fractional step method on staggered grid. To sat-
isfy the boundary conditions on an arbitrary immersed interface, the velocity field
at the grid points near the interface is reconstructed locally without smearing the
sharp interface. The complications caused by the Eulerian grid points emerging from
a moving solid body into the fluid phase are treated with a novel ?field-extension?
strategy. To treat the two-way interactions between the fluid and structure, a strong
coupling scheme based on Hamming?s fourth-order predictor-corrector method has
been developed. The fluid and the structure are treated as elements of a single
dynamical system, and all of the governing equations are integrated simultaneously,
and iteratively in the time-domain.
A variety of two and three-dimensional fluid-structure interaction problems of
increasing complexity have been considered to demonstrate the accuracy and the
range of applicability of the method. In particular, forced vibrations of a rigid
circular cylinder including the harmonic in-line vibrations in a quiescent fluid and
the transverse vibrations in a free-stream, and the vortex-induced vibrations of an
elastic cylinder with one and two degrees of freedom in a free-stream are presented
and compared with reference simulations and experiments. Three-dimensional DNS
and LES of fluid flows involving stationary complex geometries include the flow
past a sphere at Re = 50 ? 1,000, the transitional flow past an airfoil with a 10?
attack angle at Re = 10,000. Then, the turbulent flow over a traveling wavy wall
at Re = 10,170 are simulated are compared with the detailed DNS using body-
fitted grid in the literature. Finally, the simulation of the transitional flow past
a prosthetic mechanical heart valve with moving leaflets at Re = 4,000 has been
performed. All results are in good agreement with the available reference data.
AN EMBEDDED BOUNDARY FORMULATION FOR
LARGE-EDDY SIMULATION OF TURBULENT FLOWS
INTERACTING WITH MOVING BOUNDARIES
by
Jianming Yang
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2005
Advisory Committee:
Assistant Professor Elias Balaras, Chair/Advisor
Professor Richard V. Calabrese
Professor Ugo Piomelli
Associate Professor Arnaud Trouv?e
Professor James M. Wallace
c? Copyright by
Jianming Yang
2005
Dedicated to my parents,
to my wife, Guangming,
and to my son, Reece
ii
ACKNOWLEDGMENTS
I would like to express my sincere gratitude to my advisor, Prof. Elias Balaras,
for his great patience, frequent encouragement, and generous support throughout
this work. Elias not only introduced me to the large-eddy simulation of turbulence,
but also guided me into the depth and breadth of my research project. I have learned
a great deal from his invaluable insight into many of the numerical and physical
problems. This manuscript has also benefited from his countless corrections.
I am very grateful to the members of my advisory committee, Prof. Richard V.
Calabrese, Prof. Ugo Piomelli, Prof. Arnaud Trouv?e, and Prof. James M. Wallace,
for reading this dissertation and offering many constructive comments and valuable
suggestions. To Prof. Piomelli in particular, thanks are due for teaching me three
courses including almost everything in Fluid mechanics and numerical simulations.
His lectures were among the best ones I have ever attended. To Prof. Wallace,
my gratitude extends to his very careful reading and many grammar and spelling
corrections to this manuscript. I also would like to thank Prof. Trouv?e and Prof.
Piomelli for their constant input to my work in our weekly CFD group meeting.
My thanks are due to Prof. Sergio Preidikman for providing me the ODE
solver for rigid body motion and many valuable discussions. I am also grateful to
Mr. Antonio Cristallo for his help on the heart valve computation.
My life and study in Maryland has been enriched greatly by my friends and
iii
colleagues here. I would like to thank all of them for their help and support, espe-
cially, Dr. Ning Li for many useful discussions and helps during the course of this
work.
During my PhD study, my family have been very supportive. Their uncondi-
tional support has given me the strength to finish my study. I am very grateful to
them. I am also much indebted to my wife Guangming and my son Reece for all
their love and support.
This work was supported in part by the University of Maryland, College Park,
the National Institutes of Health under Grant No. R01-HL-07262, and the National
Science Foundation under Grant No. CTS-0347011. These sponsors are gratefully
acknowledged.
iv
NOMENCLATURE
Roman Symbols
a,A,Ax,Ay Oscillation amplitude.
ai,bi,ci Coefficient in quadratic polynomial.
ami,alj,ank Coefficients in the discretizated Poisson equation in cross-
stream, spanwise, and streamwise directions, respectively.
A Coefficient matrix.
Ai Operator in the ith momentum equation containing terms
treated explicitly.
Bi Operator in the ith momentum equation containing terms
treated implicitly.
c Damping coefficient; Chord length of an airfoil; Phase speed of
wave.
C Model coefficient in the Smagorinsky model.
C Damping matrix.
Cf Skin friction coefficient.
Cp Pressure coefficient.
v
CD Drag coefficient.
CL Lift coefficient.
CS Coefficient of side force.
Cx Coefficient of force in x direction.
Cy Coefficient of force in y direction.
D Diameter of cylinder/sphere.
e Local truncation error.
f Scalar function; Frequency.
f Momentum forcing vector; Hydrodynamic force vector.
f0 Shedding frequency of stationary cylinder.
fe Excitation frequency.
fi i = 1,2,3 Momentum forcing in i direction.
fn Natural vibrating frequency of a structure.
Fi i = 1,2,3 i direction force on a structure.
G Filtering operation using the grid-filter.
?G Filtering operation using the test-filter.
i,j,k Unit normal vectors.
k Time step index in time-advancement scheme; Spring constant;
Wave number.
kprimel Modified wave number in spanwise direction.
vi
K Stiffness matrix.
KC Keulegan-Carpenter number.
l Wave number in spanwise direction.
L Reference length scale.
Li Length of computational domain in i direction.
Lij Resolved turbulent stress tensor.
m Structural mass.
M Mass matrix.
n Time step index in time-advancement scheme; Normal direc-
tion of the boundary; Unit normal vector; Degrees of freedom;
Mass ratio.
Nx,Ny,Nz Number of grid points in the transverse, spanwise, and stream-
wise directions, respectively.
Nr,N?,Nz Number of grid points in the radial, azimuthal, and axial di-
rections, respectively.
p Pressure normalized by ?U2.
r A spatial ray.
R Radius of a pipe.
Q The second invariant of velocity gradient tensor.
Re Reynolds number.
s Arclength.
vii
Sij Strain rate tensor based on the resolved velocity field, ui.
St Strouhal number.
t Time.
Tij Subgrid scale stress tensor at the test-filter level.
ui Resolved velocity components: velocity components filtered at
the grid-filter level.
?ui Components of the predicted velocity field.
?uij Velocity components filtered at both the grid-filter and test-
filter level.
u? Friction velocity.
ui General notation for velocity components:
Cartesian coordinates: u1 = u,u2 = v,u3 = w.
Cylindrical coordinates: u1 = ur,u2 = u?,u3 = uz
u,v,w Cartesian coordinates: cross-stream, spanwise, and streamwise
velocity components, respectively.
Cylindrical coordinates: radial, azimuthal, and axial velocity
components, respectively.
ux,uy,uz Cartesian coordinates: cross-stream, spanwise, and streamwise
velocity components, respectively.
ur,u?,uz Cylindrical coordinates: radial, azimuthal, and axial velocity
components, respectively.
viii
ub,vb,wb Normal velocity components at the boundary of the computa-
tional domain.
ut,vt,wt Tangential velocity components at the boundary of the com-
putational domain.
U Reference velocity.
Ured Reduced velocity.
Umax Maximum velocity.
U? Free-stream velocity.
Ubulk Bulk velocity.
Ucl Centerline velocity in a pipe.
Uconv Convective velocity used in the convective boundary condition.
W Weighting function.
xi General notation for (non-dimensional) spatial coordinates:
Cartesian coordinates: x1 = x,x2 = y,x3 = z.
Cylindrical coordinates: x1 = r,x2 = ?,x3 = z.
X0,Y0,Z0 Structural displacement.
x Spatial position vector.
Xi Arclength coordinate.
x,y,z Cartesian coordinates: Spatial coordinates in the cross-stream,
spanwise, and streamwise directions, respectively.
z Cylindrical coordinates: axial coordinate.
ix
Greek Symbols
?i Factor in RK time-advancement scheme.
?ij Kronecker?s delta.
? Discrete form of the partial operator.
? Filter width.
? Filter width at the grid-filter level.
?? Filter width at the test-filter level.
?x,?y,?z Cartesian coordinates. Grid-spacing in the cross-stream, span-
wise, and streamwise directions, respectively.
?r,??,?z Cylindrical coordinates. Grid-spacing in the radial, azimuthal,
and axial directions, respectively.
?t Time step.
?i Factor in RK time-advancement scheme.
? Wavelength.
? Dynamic viscosity of a fluid.
? Total viscosity: ? = 1/Re+?t.
?m Kinematic viscosity of a fluid.
?t Turbulent eddy viscosity.
x
? Angular frequency oscillation; Vorticity.
? Scalar function.
? Scalar to project the predicted velocity field into a divergence-
free space.
? Immersed interface.
? Density of a fluid.
?i Factor in RK time-advancement scheme.
? Time
?ij Subgrid scale stress tensor at the grid-filter level.
?w Shear stress at the wall.
? Cylindrical coordinates: Azimuthal coordinate.
? Structural rotation.
?,?,? Cartesian coordinates: spatial coordinates in the computa-
tional space in the cross-stream, spanwise, and streamwise di-
rections, respectively.
Cylindrical coordinates: spatial coordinates in the computa-
tional space in the the radial, azimuthal, and axial directions,
respectively.
? Damping ratio.
xi
Ohter Symbols
?() First derivative.
?() Second derivative.
<> Averaging operator.
()+ Variable given in wall units: ()+ = ?()u?/?.
Abbreviations
AB Adams-Bashforth.
CFL Courant-Friedrichs-Levy.
CN Crank-Nicholson.
CPU central processing unit.
DNS direct numerical simulation.
FFT fast Fourier transform.
LES large-eddy simulation.
MPI message passing interface.
ODE ordinary differential equation.
RHS right-hand side.
RK Runge-Kutta.
SGS subgrid scale.
xii
TABLE OF CONTENTS
List of Tables xv
List of Figures xvi
1 Introduction 1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Literature Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Present Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2 Basic Navier-Stokes Solver 12
2.1 Formulation in Cartesian Coordinates . . . . . . . . . . . . . . . . . . 12
2.1.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 12
2.1.2 Subgrid Scale Stress Models . . . . . . . . . . . . . . . . . . . 14
2.1.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.4 Time Advancement Scheme . . . . . . . . . . . . . . . . . . . 21
2.2 Transformation in Cylindrical Coordinates . . . . . . . . . . . . . . . 22
2.2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.2 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . 26
2.2.3 Time Advancement Schemes . . . . . . . . . . . . . . . . . . . 32
2.3 Poisson Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.4.1 Dirichlet and Neumann Boundary Conditions . . . . . . . . . 37
2.4.2 Convective Boundary Condition . . . . . . . . . . . . . . . . . 39
2.4.3 Periodic Boundary Conditions . . . . . . . . . . . . . . . . . . 39
2.5 Validation of Basic Solver . . . . . . . . . . . . . . . . . . . . . . . . 40
3 Embedded Boundary Method 49
3.1 Interface-grid Relation . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.1 Interface Description . . . . . . . . . . . . . . . . . . . . . . . 49
3.1.2 Tagging of Points on the Eulerian Grid . . . . . . . . . . . . . 51
3.1.3 Establishment of Interface-Normal Intersections . . . . . . . . 53
3.2 Treatment of Stationary Immersed Boundaries . . . . . . . . . . . . 56
3.2.1 Momentum Forcing . . . . . . . . . . . . . . . . . . . . . . . . 56
3.2.2 Generalized Interpolation Scheme . . . . . . . . . . . . . . . . 58
3.2.3 Treatment of Eddy Viscosity . . . . . . . . . . . . . . . . . . . 63
3.3 Treatment of Moving Immersed Boundaries . . . . . . . . . . . . . . 65
3.3.1 Complications in Time Advancement . . . . . . . . . . . . . . 65
3.3.2 Field Extension . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.3.3 Numerical Procedure . . . . . . . . . . . . . . . . . . . . . . . 69
3.4 Surface Force Calculation . . . . . . . . . . . . . . . . . . . . . . . . 71
3.5 Fluid-Structure Interaction . . . . . . . . . . . . . . . . . . . . . . . . 75
3.6 Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
xiii
4 Two-Dimensional Validation Cases 83
4.1 Formal Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.2 Forced Vibrations of a Cylinder . . . . . . . . . . . . . . . . . . . . . 88
4.2.1 In-Line Oscillating Cylinder in a Fluid at Rest . . . . . . . . . 88
4.2.2 Transversely Oscillating Cylinder in a Free-Stream . . . . . . . 96
4.3 Vortex-Induced Vibrations of an Elastic Cylinder . . . . . . . . . . . 107
4.3.1 One Degree of Freedom . . . . . . . . . . . . . . . . . . . . . . 107
4.3.2 Two Degrees of Freedom . . . . . . . . . . . . . . . . . . . . . 117
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5 Three-Dimensional Applications 130
5.1 Flow Past a Sphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.1.1 Axisymmetric Flow . . . . . . . . . . . . . . . . . . . . . . . . 130
5.1.2 Planar-Symmetric Flow . . . . . . . . . . . . . . . . . . . . . 134
5.1.3 Transitional Flow . . . . . . . . . . . . . . . . . . . . . . . . . 140
5.2 Transitional Flow Past an Airfoil . . . . . . . . . . . . . . . . . . . . 150
5.3 Turbulent Flow Over a Traveling Wavy wall . . . . . . . . . . . . . . 154
5.4 Pulsatile Flow Past a Prosthetic Bileaflet Heart Valve . . . . . . . . . 165
6 Conclusions and Future Directions 175
6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
6.2 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178
Bibliography 181
xiv
LIST OF TABLES
2.1 Grid resolution in wall units compared with data in [1] . . . . . . . . 43
2.2 Mean flow properties compared with data in [1] . . . . . . . . . . . . 43
5.1 Prediction of CD for axisymmetric flow past the sphere . . . . . . . . 132
xv
LIST OF FIGURES
2.1 A sketch of grid cell and the variable arrangement in Cartesian coor-
dinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Staggered grid and variable arrangement in x?z coordinates. . . . . 20
2.3 A sketch of grid cell and the variable arrangement in cylindrical co-
ordinates. a) A typical grid cell; b) a grid cell near the centerline. . . 27
2.4 Staggered grid and the variable arrangement in r?? plane. ? p, uz;
square ur; triangle u?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5 Centerline treatment of the cylindrical coordinates. a) The required
boundary conditions and variable collacation on the ghost cell; b)
Variables defined across the centerline on r?? plane. . . . . . . . . . 30
2.6 The implementation of Dirichlet and Neumann boundary condition. . 37
2.7 The implementation of periodic boundary condition. Arrows show
the directions of data duplication. . . . . . . . . . . . . . . . . . . . . 40
2.8 Typical computational grid. Shaded area is the core-region; un-
shaded area is the outer-region. . . . . . . . . . . . . . . . . . . . . . 42
2.9 Mean streamwise velocity profile in wall unit. DNS data from [16] . . 44
2.10 Radial turbulence intensity. DNS data from [16] . . . . . . . . . . . . 45
2.11 Azimuthal turbulence intensity. DNS data from [16]. . . . . . . . . . 45
2.12 Axial turbulence intensity. DNS data from [16]. . . . . . . . . . . . . 46
2.13 Total stress balance. DNS data from [16]. . . . . . . . . . . . . . . . . 46
2.14 Ratio of eddy viscosity to molecular viscosity. . . . . . . . . . . . . . 47
2.15 Instantaneous field in a polar plane. (a) Axial velocity; (b) Axial
vorticity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
2.16 Instantaneous vortical structure by the Q-criterion (colored by axial
vorticity). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.1 The parametrized description of interfaces of arbitrary shapes using
marker particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
xvi
3.2 Grid?interface relation. (a) Parametrized interface immersed in the
underlying Cartesian grid; (b) Zoom in the vicinity of interface where
the inside/outside status of the Eulerian grid points are shown. square
Fluid points; squaresolid Solid points. . . . . . . . . . . . . . . . . . . . . . . . 52
3.3 Identification of boundary points. (a) triangle Forcing points, square fluid
points, and squaresolid solid points for momentum forcing procedure; (b) trianglesolid
Pseudo-fluid points, square fluid points, and squaresolid solid points for field ex-
tension procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.4 Schematicofthesolutionprocedureforinterface-normalintersections.
triangle Forcing points, trianglesolid Pseudo-fluid points, square fluid points, and squaresolid solid
points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.5 Previous Interpolation schemes. (a) One-dimensional scheme in [15];
(b) Two-dimensional scheme in [4]. Cases (1) and (3) illustrate two
possible interpolation stencils depending on the interface topology
and local grid size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.6 Extra ghost cells required for the interpolation scheme in [4] when
domain decomposition technique used for parallelization. . . . . . . . 60
3.7 Generalized Interpolation stencil. . . . . . . . . . . . . . . . . . . . . 61
3.8 Treatment of turbulent eddy viscosity for LES. (a) Test-grid cells
for the boundary points involving solid points from the interior; (b)
Interpolation stencils for ?t at the boundary points based the well-
definedvaluesfromthefluidpointsawayfromtheimmersedboundary
and the boundary condition for ?t at the interface. . . . . . . . . . . . 63
3.9 A simplified one-dimensional model problem that demonstrates the
possible changes of flags of the Eulerian grid nodes near a moving
interface. (a) Solid phase encroaches upon the fluid phase; (b) Solid
phase withdraws from the fluid phase. . . . . . . . . . . . . . . . . . . 66
3.10 The field-extension procedure. Shaded triangles are possible extrap-
olation stencils for the pseudo-fluid points at substep tk?1 where the
solution will be extended. . . . . . . . . . . . . . . . . . . . . . . . . 67
3.11 Extension of pressure field. . . . . . . . . . . . . . . . . . . . . . . . . 69
xvii
3.12 Surface force calculation. (a) Arrangement of flow variables on the
staggered grid: triangleright u velocity, triangle v velocity, square pressure, ? marker
particles. The corresponding filled symbols on the interface represent
the location where normal intersects the boundary. (b) detail in the
vicinity of the boundary for the grid of u velocity component. (c)
detail in the vicinity of the boundary for the the grid of pressure.
Interpolation stencil is the shaded triangle and line with arrow is the
normal to the interface. . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.13 Schematic of rigid-body motion. (a) One degree of freedom of transla-
tion in vertical direction; (b) One degree of freedom of rotation about
the origin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.14 Domain decomposition parallelization. (a) Slab decomposition in
streamwise direction for momentum equations; (b) Slab decompo-
sition in wave-number space for the Poisson equation. . . . . . . . . . 80
3.15 Inter-slab communication via layers of ghost cells. . . . . . . . . . . . 81
3.16 Parallel Speedup for a 96?64?384 Grid. . . . . . . . . . . . . . . . 82
4.1 Flow around a periodic array of cylinders in a plane channel: Sketch
of the computational domain. . . . . . . . . . . . . . . . . . . . . . . 84
4.2 Streamlines on a 150?150 uniform grid colored by u velocity compo-
nent (a reference frame moving with the body is used for the moving
cylinder case). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 Instantaneous pressure contours at two different locations. . . . . . . 86
4.4 a50 L1 norm of the error; ? L2 norm of the error. Open symbols are
for w and filled for u. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.5 In-line oscillating cylinder in a fluid at rest (Re = 100 and KC = 5).
Pressure and vorticity contours at four different phase-angles. ?1.1 <
P < 0.6 with intervals of 0.09 and ?26 < ?y < 26 with intervals of
0.95. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.6 Measured (left) and computed (right) velocity vectors and streamlines
in the vicinity of the cylinder at Re = 100 and KC = 5. . . . . . . . . 90
4.7 In-line oscillating cylinder in a fluid at rest (Re = 100 and KC = 5).
Profiles of u velocity component at: (a) 180?; (b) 210?; (c) 330?.
Symbols: experiment; Lines: Computation. . . . . . . . . . . . . . . 91
xviii
4.8 In-line oscillating cylinder in a fluid at rest (Re = 100 and KC = 5).
Profiles of w velocity component at: (a) 180?; (b) 210?; (c) 330?.
Symbols: experiment; Lines: Computation. . . . . . . . . . . . . . . . 92
4.9 Convergence of the u velocity component for different grid levels at
Re = 100 and KC = 5. . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.10 Convergence of the w velocity component for different grid levels at
Re = 100 and KC = 5. . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.11 Time history of in-line force . . . . . . . . . . . . . . . . . . . . . . . 96
4.12 Cylinder oscillating in a free-stream: instantaneous streamline pat-
terns when cylinder is located at its extreme upper position. Stream-
lines are colored with transverse velocity component values. . . . . . . 97
4.13 Cylinder oscillating in a free-stream: instantaneous vorticity contours
when cylinder is located at its extreme upper position. ?8 < ?y < 8
with intervals of 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.14 Cylinder oscillating in a free-stream: instantaneous pressure contours
when cylinder is located at its extreme upper position. . . . . . . . . 100
4.15 Evolution of the drag and lift coefficients in time . . . . . . . . . . . . 102
4.16 Comparison of force coefficients and phase angle . . . . . . . . . . . . 104
4.17 Distribution of the pressure coefficients on the cylinder?s surface,
when located at the extreme upper position. grid 800 ? 640;
grid 400?320; grid 200?160; ? body-fitted reference com-
putation in [23]. (a) fe/f0 = 0.80; (b) fe/f0 = 0.90; (c) fe/f0 = 1.00;
(d) fe/f0 = 1.10; (e) fe/f0 = 1.12; (f) fe/f0 = 1.20. . . . . . . . . . . 105
4.18 Distribution of the skin friction coefficients on the cylinder?s surface,
when located at the extreme upper position. grid 800 ? 640;
grid 400?320; grid 200?160; ? body-fitted reference com-
putation in [23]. (a) fe/f0 = 0.80; (b) fe/f0 = 0.90; (c) fe/f0 = 1.00;
(d) fe/f0 = 1.10; (e) fe/f0 = 1.12; (f) fe/f0 = 1.20. . . . . . . . . . . 106
4.19 Time history of the displacements for the freely vibrating cylinder in
the cross-stream direction: start-up phase. . . . . . . . . . . . . . . . 109
4.20 Time history of the displacements for the freely vibrating cylinder in
the cross-stream direction: stable periodic phase. . . . . . . . . . . . 111
4.21 Freely vibrating cylinder in the cross-stream direction: instantaneous
vorticity contours. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
xix
4.22 Freely vibrating cylinder in the cross-stream direction: instantaneous
pressure contours. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
4.23 Comparison of the force coefficients for the cases of fixed cylinder and
the freely oscillating cylinder. ? fixed cylinder; a50 freely oscillating
cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
4.24 Comparison of the oscillation frequency and the maximum oscillation
amplitude. ? present computations; diamondmath experiment in [2]; a50 computa-
tion in [53]; ? computation in [31] . . . . . . . . . . . . . . . . . . . . 116
4.25 Vortex-induced vibration of an elastic cylinder with two degrees of
freedom: instantaneous vorticity contours. . . . . . . . . . . . . . . . 119
4.26 Vortex-induced vibration of an elastic cylinder with two degrees of
freedom: instantaneous pressure contours. . . . . . . . . . . . . . . . 120
4.27 X-Y phase plots. X and Y are in different scales. . . . . . . . . . . . . 122
4.28 Time history of the displacements and lift coefficients. . . . . . . . . . 123
4.29 Phase plots of transverse displacement, lift coefficient. . . . . . . . . . 124
4.30 Phase plots of drag, lift coefficients. . . . . . . . . . . . . . . . . . . . 125
4.31 Theoscillationfrequency,transverseoscillationamplitudeandstream-
wise oscillation amplitude as functions of reduced velocity. diamondmath two de-
grees of freedom cases; ? one degree of freedom cases from Sec. 4.3.1;
a50 one degree of freedom cases from [53] . . . . . . . . . . . . . . . . . 127
4.32 Comparison of force coefficients of the freely oscillating cylinders with
one degree of freedom and two degrees of freedom. a50 two degrees of
freedom cases; ? one degree of freedom cases from Sec. 4.3.1 . . . . . 128
5.1 Part of grid 40 ? 40 ? 100 for axisymmetric flow past the sphere.
Uniform grid in the red box. . . . . . . . . . . . . . . . . . . . . . . . 131
5.2 Error in the drag coefficient for the axisymmetric flow past the sphere
as a function of grid resolution. . . . . . . . . . . . . . . . . . . . . . 133
5.3 Axisymmetric streamlines past the sphere (colored by the streamwise
velocity): (a) Re = 50; (b) Re = 100; (c) Re = 150; (d) Re = 200. . . 135
5.4 Vorticitycontoursforaxisymmetricflowpastthesphere: (a)Re = 50;
(b) Re = 100; (c) Re = 150; (d) Re = 200. . . . . . . . . . . . . . . . 136
xx
5.5 Streamlines of projected velocity field for flow past the sphere at
Re = 250 (colored by the streamwise velocity): (a) x?z plane; (b)
y ?z plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.6 Surface pressure distribution on the sphere at Re = 250: (a) front
view; (b) rear view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
5.7 Vorticity contours for flow past the sphere at Re = 250: (a) Spanwise
vorticity in x?z plane; (b) Spanwise vorticity in y ?z plane. . . . . 139
5.8 Time history of drag, x, and y side force coefficients at Re = 300. . . 141
5.9 Spanwise vorticity contours for flow past the sphere at Re = 300 at
the four quarter periods in the x?z plane. . . . . . . . . . . . . . . . 142
5.10 Spanwise vorticity contours for flow past the sphere at Re = 300 at
the four quarter periods in the y ?z plane. . . . . . . . . . . . . . . . 143
5.11 Instantaneous vortical structures for flow past the sphere at Re = 300
(colored with the streamwise vorticity) at the four quarter periods:
x?z view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
5.12 Instantaneous vortical structures for flow past the sphere at Re = 300
(colored with the streamwise vorticity) at the four quarter periods:
y ?z view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.13 Instantaneous azimuthal vorticity contours for flow past the sphere
at Re = 1000: (a) x?z plane; (b) y ?z plane. . . . . . . . . . . . . . 146
5.14 Instantaneous vortical structures for flow past the sphere at Re =
1000 (colored with the streamwise vorticity): (a) x?z view; (b) y?z
view; (c) oblique view. . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.15 Time history of drag, x, and y side force coefficients and frequency
spectrum of CD at Re = 1000. . . . . . . . . . . . . . . . . . . . . . . 148
5.16 Time history and frequency spectrum of ur at r = 0.3D, z = 2.0D. . . 149
5.17 Time history and frequency spectrum of uz at r = 0.3D, z = 2.0D. . 149
5.18 Part of computational grid in x?z plane around the airfoil. Every
three grid lines are shown. . . . . . . . . . . . . . . . . . . . . . . . . 151
5.19 Mean streamlines of the flow past an E387 airfoil. Streamlines are
colored with streamwise velocity component. . . . . . . . . . . . . . . 151
5.20 Mean pressure coefficient on the airfoil surface. . . . . . . . . . . . . . 152
xxi
5.21 Time history of drag and lift coefficients for the flow past an airfoil. . 153
5.22 Frequency spectra of the drag and lift coefficients for the flow past
an airfoil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
5.23 Turbulent intensities, turbulent shear stress, and mean eddy viscosity
for the flow past an airfoil. . . . . . . . . . . . . . . . . . . . . . . . 155
5.24 Close-up views of the streamwise turbulent intensity, turbulent shear
stress < uw >, and mean eddy viscosity < ?t > for the flow past an
airfoil. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
5.25 Instantaneous contours of the spanwise and streamwise vorticity com-
ponents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
5.26 Instantaneous vortical structures using isosurfaces of Q = 10 colored
by the streamwise vorticity for the flow past an airfoil. . . . . . . . . 157
5.27 A sketch of the geometry for the case of turbulent flow over a traveling
wavy wall. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
5.28 Meanstreamlinepatterninaframeofreferencefollowingthetraveling
wave. (a) c/U = 0.0; (b) c/U = 0.4; (c) c/U = 1.2. ? indicates the
center of the corresponding recirculation region in the reference DNS
[54]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
5.29 Phase-averaged statistics for c/U = 0.4. . . . . . . . . . . . . . . . . . 161
5.30 Phase-averaged statistics for c/U = 1.2. . . . . . . . . . . . . . . . . . 163
5.31 Mean force acting on the traveling wavy boundary as a function of
c/U. , ? is total friction force, Ff; , a50 is the total pressure
force, Fp. Lines are the DNS results in [54], and symbols are the
present results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
5.32 Isosurfaces of Q = 8.0 colored with the streamwise vorticity at an
instant in time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
5.33 Geometry of the simplified bileaflet prosthetic heart valve placed in
a rigid wall model of the aortic root. The leaflets are shown in their
fully open position, and part of the housing wall has been removed
for clarity. Two planes of the Cartesian grid are also shown. . . . . . 167
5.34 Imposed flow rate and leaflet kinematics for the heart valve problem:
(a) Bulk velocity at the inlet; (b) opening angle of the leaflets, as
functions of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
xxii
5.35 Instantaneous flow fields at an x?z plane cutting through the middle
of the leaflets. Left: streamwise velocity (?1.5 < u/Umaxb < 3, in-
tervals are 0.05); Right: spanwise vorticity (?20 < ?yD/Umaxb < 20,
intervals are 0.2). The phase reference, t/T, is from Fig. 5.34. . . . . 170
5.36 Instantaneous vortical structures visualized using isosurfaces of Q =
600 colored by the streamwise vorticity at: (a) t/T = 0.3; (b) t/T =
0.4; (c) t/T = 0.8; (d) t/T = 0.9. The phase reference, t/T, is from
Fig. 5.34. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172
xxiii
Chapter 1
Introduction
1.1 Background
Numerical simulations of turbulent and transitional flows with dynamically
moving boundaries are amongst the most challenging problems in computational
mechanics. The main difficulty comes from the fact that the spatial domain occu-
pied by the fluid changes with time, and the location of the boundary is usually
an unknown itself that depends on the fluid flow and the solution of additional
equations describing the motion/deformation of the body. There is only a limited
number of special cases where established boundary-conforming formulations can
be directly applied to such problems with a relatively small overhead. The use of
moving reference frames [31], or coordinate transformations [40, 54] are characteris-
tic examples. In more complex configurations formulations utilizing moving and/or
deforming grids that continuously adapt to the changing location of the body have
to be adopted (see for example [53, 24] ). For problems that involve multiple bodies
undergoing large motions and/or deformations, these algorithms are fairly compli-
cated and have an adverse impact on the accuracy and efficiency of the fluid solvers.
Therefore, although a variety of fluid-structure interaction algorithms has been
developed over the years, relatively few applications in turbulent and transitional
flows have been reported. In most cases, this is due to prohibitively high computa-
1
tional cost, or dissipative discretizations that limit the applicability of such methods
to classical turbulence modeling strategies. Further advancements in this field can be
achieved by coupling state-of-the-art tools to model turbulence and transition (i.e.
large-eddy simulations (LES) or hybrid formulations) with cost-effective numerical
methods applicable to problems with large boundary motions and deformations.
An alternative class of methods that has the potential to overcome some of the
above limitations are non-boundary-conforming formulations. In such methods the
requirement that the grid conforms to a solid boundary is relaxed, and the effect of
a complex object on the flow is introduced through proper treatment of the solution
variables at the grid cells in the vicinity the body. The basic advantage of these
formulations is the simplification of grid generation, especially in cases of moving
boundaries where the need for regeneration or deformation of the grid is eliminated.
In addition, efficient, Cartesian solvers can be directly applied to complex flow prob-
lems. Both the above features are particularly attractive for Direct Numerical Sim-
ulations (DNS) or LES of turbulent and transitional flows, where the use of highly
efficient, energy conserving solvers is imperative for accurate computations. It is
therefore conceivable that successful integration of non-boundary-conforming strate-
gies with robust Cartesian or cylindrical coordinate solvers developed for DNS/LES
will open a wide new area of applications for these tools. Example applications of
such a tool include a variety of low and moderate Reynolds number turbulent flow
problems from engineering, biology, and medicine, where fluid-structure interactions
are central to the dynamics of the flow.
The objective of this work is the development of an embedded-boundary
2
method for DNS and LES of turbulent and transitional flows interacting with dy-
namically moving boundaries. The target application fields are the a variety of
low and moderate Reynolds number turbulent flows with complex geometries and
moving boundaries in engineering, biology, and medicine.
1.2 Literature Survey
Over the past decades a variety of non-boundary-conforming methods with
various degrees of accuracy and complexity have been proposed. The so-called
immersed-boundary formulation pioneered by Peskin [43] represents a family of
methods where a set of body forces is used to represent the effect of an object to the
flow. Initially the method was used to study fluid-structure interaction problems in
the cardiovascular circulation [44, 45]. In these computations the vascular boundary
was modeled as a set of elements linked by springs. As a result the forces required
to enforce boundary conditions could be evaluated in a straightforward manner (i.e.
Hooke?s law). In the case of rigid immersed bodies, however, the corresponding
forces are not known a priori and must be calculated by some feedback algorithm.
Lai and Peskin [30] suggested a formulation where the body is allowed to move a lit-
tle -rather than being fixed- by connecting it to a very stiff spring. They tested this
approach for the flow around a cylinder with satisfactory results, although the spec-
ification of the stiffness constant is somehow ad hoc. Goldstein et al. [19] introduced
an alternative approach where a feedback-forcing scheme is used to asymptotically
enforce the desired boundary conditions on a solid boundary. Application of the
3
method to three-dimensional computations of turbulent flow in plane and ribbed
channels yielded results in good agreement with reference data [20, 21]. An advan-
tage of immersed-boundary formulations is their straightforward implementation in
existing solvers (modifications are confined to the RHS of the equations of motion).
On the other hand, the need for a smooth transition between the fluid and the solid
body spreads the forcing function over several grid cells and introduces some blur-
ring between the two regions. This feature can decrease the order of accuracy of the
scheme near the body or increase the resolution requirements making their use in
turbulent flows problematic.
Another class of methods which does not suffer from the ?blurring? men-
tioned above are the so-called Cartesian or cut-cell formulations. In this case the
solid boundary is tracked as a sharp interface and the grid cells at the body in-
terface are modified according to their intersections with the underlying Cartesian
grid. The discrete operators at these cells are then modified to reflect the desired
boundary conditions. Successful applications of such methods in two dimensional
flows with stationary boundaries can be found in [73, 7, 52, 64, 71]. However, due
the large number of possible intersections between the grid and the boundary a
variety of interface-cells is generated leading to an equally large number of special
treatments. Also in complex configurations the unavoidable generation of irregularly
shaped cells with very small size can have an adverse impact of the conservation and
stability properties of the solver. Recently Ye et al. [71] suggested a cell merging
scheme to address this problem. This formulation was also extended to treat moving
boundaries with good results for a variety of two-dimensional problems [63]. The
4
extension of the methodology in complex three-dimensional configurations remains
to be investigated.
Recently, Fadlun et al. [15] introduced an embedded-boundary method, which
shares a number of features with both approaches discussed above. As in immersed-
boundary formulations the method still utilizes a force-field to enforce boundary
conditions. In this case, however, the forces are not specified in the continuous space
by means of some physical arguments, but rather in the discrete space by directly
requiring the solution to respect the desired boundary conditions. This process is
equivalent to a local reconstruction of the solution near the interface and enforces the
desired boundary conditions ?exactly?, as in cut-cell formulations. The encouraging
results reported in [65, 15] together with the straightforward implementation of the
method in existing Navier-Stokes solvers, motivated a number of recent studies,
where alternative embedded-boundary formulations based on the principals outlined
above have been proposed. The main difference between them is the way the solution
is reconstructed near the interface. In Fadlun et al. [15] and Balaras [4], for example,
the solution is reconstructed at the fluid nodes closest to the immersed boundary
(fluid points with at least one neighbor in the solid phase). In the former study a one-
dimensionalinterpolationschemealonganarbitrarygridlineisusedforthispurpose,
while in the latter the reconstruction is performed along the well defined line normal
to the interface. In Kim et al. [29], Majumdar et al. [35], or Tseng and Ferziger [61]
on the other hand, the solution is reconstructed at ?ghost-cells?, which are points
inside the solid phase with at least one neighbor in the fluid phase. Despite the
different strategies, however, both approaches have the velocity boundary conditions
5
?implicitly? build into the reconstruction stencil, which results in a sharply defined
interface.
A related issue, which also has implications on the overall accuracy of the
approach, is the way the direct-forcing is incorporated into the time advancement
scheme (fractional-step methods are almost exclusively used in the all the above
studies). In [15], for example, the use of a one-dimensional interpolation strat-
egy facilitates the imposition of boundary conditions to the predicted velocity by
directly modifying the coefficients of the standard linear system. In [4] where a
multidimensional reconstruction is used in the framework of an explicit time ad-
vancement scheme the same result can be simply achieved by directly modifying the
predicted velocity field. As a result in both cases the actual value of the forcing
function is not explicitly computed and the boundary conditions on the immersed
interface are satisfied by a direct modification of the discrete operators. When a
semi-implicit scheme (i.e. third-order Runge-Kutta scheme for advective terms, and
Crank-Nikolson for the viscous terms) is used, however, the multi-dimensional re-
construction stencil involves values from the yet to be determined predicted velocity.
To implicitly satisfy the boundary conditions, as in [15, 4], one needs to modify the
standard linear system of equations in the predictor step, which is now going to
become a sparse system with a costly solution. To alleviate this problem, Kim et
al. [29], proposed an approach where the forcing is estimated by taking a prelim-
inary fully explicit ?predictor? step (a forward Euler step for the viscous terms in
their case). The estimated forcing function is then incorporated into the RHS of the
regular linear system to enforce the desired boundary conditions. Balaras and Yang
6
[6] used a similar approach in the framework of a semi-implicit, finite-difference,
fractional step method in cylindrical coordinates with good results.
The boundary motion/deformation over a fixed grid is usually the source of
additional complications in most of the non-boundary-conforming strategies that
were discussed in the previous paragraphs. In classical immersed-boundary formu-
lations [43, 30], such problems are usually minimal due to the smooth transition
from the solid to the fluid. In cut-cell [73, 7, 52, 64, 71] or embedded-boundary
formulations [15, 29, 61, 4], however, complications are encountered due the fact
that the role of the Eulerian grid points near the interface changes from timestep
to timestep (for example a point in the solid body can became a point in the fluid
at the next timestep and vise-versa), as the body moves through the fixed grid. As
a result the velocity and pressure for some points in the flow will get non-physical
values due to their previous association with the solid phase. In the framework
of cut-cell formulations Udaykumar et al. [63] proposed a cell merging scheme to
handle the cells that emerge from the body which they call ?freshly-cleared? cells.
In embedded-boundary formulations on the other hand, no systematic study that
identifies and addresses problems associated with large boundary motions has been
reported. Due to the explicit dependence of such methods on the details of the
adopted numerical method it is difficult to formulate a general approach. In the
case of ?ghost-cell? methods, for example, as the body moves through the fixed grid
some of the ?ghost-cells? will emerge into the fluid and will became fluid cells. Since
they were previously in the solid they have no history in the fluid phase and no
physically realizable value for the velocity and pressure at the previous timestep(s).
7
As a consequence their treatment to some degree would have to be ad-hoc. In
formulations where the solution is reconstructed at the fluid nodes closest to the
boundary, (i.e. [15]), the points that emerge from the solid become the boundary
points that are central to the reconstruction procedure, and therefore their history
in the fluid phase is irrelevant. The points that require special treatment in this
case are the boundary points that move further into the fluid.
The later approach is the one that will be discussed in detail in the present
study, where a variance of the embedded-boundary formulation proposed in [4] will
be extended to moving boundary problems. To address the problems due to bound-
ary motion we will propose a field-extension strategy that practically extends the
pressure and velocity fields into the solid body to implicitly treat problematic grid
cells in a robust manner. We will also discuss a straightforward method to compute
the local force distribution on a complex body that is based on the same recon-
struction stencil used for the velocity field. As it will be demonstrated in the results
section the reconstructed local forces are in good agreement with reference results
obtained with boundary-conforming formulations.
1.3 Present Contributions
The major contributions of this work are summarized as follows:
? Extension of the Cartesian, finite-difference, fractional-step, Navier-Stokes
solver [3] in cylindrical coordinates, including dynamic SGS models for LES.
? Development of a generalized local reconstruction scheme for immersed bound-
8
ary treatment in the embedded-boundary formulation.
? Development of a novel field-extension strategy for treating grid points emerg-
ing from the solid into the fluid phase.
? Development of a strong coupling scheme for the simulation of Fluid-Structure
Interaction problems.
? First comprehensive validation of non-boundary-conforming methods on sur-
face force distribution for moving boundary problems.
? First simulation of transitional flow past prosthetic mechanical heart valve
with moving leaflets in a realistic geometry.
1.4 Outline
The outline of this dissertation is as follows:
? In Chap. 2, the basic fluid solver is described in detail. The governing equa-
tions, the subgrid-scale stress model, spatial discretization and time advance-
ment schemes in both Cartesian and cylindrical coordinates are all documented
there. The discussion of boundary conditions is included. A validation of the
basic fluid solver is performed for the LES of fully developed turbulent pipe
flow.
? Chap. 3 presents the embedded-boundary formulation, the strong coupling
scheme for fluid-structure interaction problems, and the parallelization of the
9
whole algorithm. For the embedded-boundary formulation, the establishment
of the interface-grid relation, the treatment of stationary immersed bound-
aries and moving immersed boundaries are given in detail. A novel surface
force calculation procedure, which computes the hydrodynamic loads to the
structure, is also covered.
? A series of two-dimensional cases with increasing complexity are presented in
Chap. 4. First, the two dimensional flow around one element of a moving
cylinder-array in a planar channel is computed to demonstrate the formal
accuracy of the method. Then, laminar flow problems involving prescribed
motions of two-dimensional bodies are simulated: the flow induced by the
harmonic in-line oscillation of a circular cylinder in a quiescent fluid and the
flow from a transversely oscillating cylinder in a free-stream. Finally, the
vortex-induced vibrations of circular cylinders are used in the present study to
evaluate the accuracy and efficiency of the proposed fluid-structure interaction
scheme. Two different configurations are considered: one degree-of-freedom
oscillations in the cross-stream direction and two degree-of-freedom oscillations
in both the streamwise and cross-stream directions.
? In Chap. 5, four three-dimensional large-scale simulations are presented. The
first case is the flow past a sphere at Reynolds numbers from 50, to 300,
and 1000. The transitional flow past an airfoil at Reynold number 10,000
is also simulated. Then, LES of flow over a traveling wave are presented in
comparison to reference DNS results. To demonstrate the robustness and
10
applicability of the method in highly unsteady turbulent flows that involve
multiple moving complex boundaries, the flow around a bileaflet prosthetic
heart valve is also presented.
? Finally, in Chap. 6 conclusions and future directions are given.
11
Chapter 2
Basic Navier-Stokes Solver
The aim of this study is to develop a computational tool for DNS and LES of
turbulent flows interacting with complex moving boundaries. Rigid bodies under-
going motion that is either prescribed or governed by additional ODEs that have to
be solved as a coupled system together with the Navier-Stokes equations are consid-
ered. The overall solver will be given in two parts. In this chapter the basic solver
in Cartesian and cylindrical coordinates is discussed in detail. A brief validation
for fully developed turbulent pipe flow is also included. The treatment of immersed
boundary and fluid-structure coupling scheme will be discussed in the next chapter.
2.1 Formulation in Cartesian Coordinates
2.1.1 Governing Equations
The governing equations for unsteady, incompressible, viscous flow of a Newto-
nian fluid with constant density can be written in Cartesian coordinates as follows:
?ui
?t +
?(uiuj)
?xj = ?
?p
?xi +
1
Re
?2ui
?xj?xj, (2.1)
?ui
?xi = 0, (2.2)
12
where xi and xj (i,j = 1,2,3) are the Cartesian coordinates, ui and uj are the
velocity components in the corresponding directions, normalized by a reference ve-
locity U, t is the time normalized by L/U with L the reference length scale, p is the
pressure normalized by ?U2 with ? the density of the fluid, and Re is the Reynolds
number defined as Re = ?UL/? with ? the dynamic viscosity of the fluid.
For the laminar flow cases or DNS of turbulent flow cases considered in this
study, the above equations are directly integrated in space and time without in-
troducing any model. In the LES approach, however, a spatial filtering operation
is applied and it separates the large, energy carrying eddies, which are directly re-
solved, from the small scales, which are modeled. In the present finite-difference
implementation a top-hat filter is implicitly applied by the discrete operators. The
resulting equations governing the evolution of the large scales have the following
form:
?ui
?t +
?(uiuj)
?xj = ?
?p
?xi +
1
Re
?2ui
?xj?xj ?
??ij
?xj , (2.3)
?ui
?xi = 0, (2.4)
where the overbar denotes a filtered variable and
?ij = uiuj ?uiuj, (2.5)
are the subgrid scale (SGS) stresses.
13
2.1.2 Subgrid Scale Stress Models
The spatial filtering operation on the Navier-Stokes equations produces the
subgrid scale stresses, which have to be parametrized in order to close the system
of equations (2.3) and (2.4).
The most widely used model for the subgrid scale stresses is the Smagorinsky
eddy viscosity model [56]. In this model, the subgrid scale stresses are related to
the resolved strain rate tensor as follows:
?ij ? 13?ij?kk = ?2?tSij = ?2C?2|S|Sij (2.6)
where ?ij is Kronecker?s delta, Sij is the resolved strain rate tensor,
Sij = 12
parenleftbigg?u
i
?xj +
?uj
?xi
parenrightbigg
. (2.7)
with magnitude:
|S| =
radicalBig
2SijSij, (2.8)
?t is the turbulent eddy viscosity, ? is the filter size, defined as (?x?y?z)1/3 with
?x, ?y, and ?z the local grid size in x, y, and z coordinate direction, respectively.
In the model C is a user-specified constant which usually varies between 0.1 ?
0.2 depending on the flow [32, 12]. This, together with the fact that the model
does not vanish in the laminar flow regions and does not have the proper limiting
behavior near solid boundaries make its use in complex flows problematic.
The dynamic modeling procedure [17] addresses most of the above deficiencies.
The basic idea behind the family of dynamic models is to use information available
at the resolved scales to predict the unresolved subgrid scales by the resolved scales.
14
An additional spatial filter, or test filter, is applied to the Navier-Stokes equations
for this purpose, which are filtered by the grid filter. The filter width of the test-
filter, hatwide?, is larger (usually double size) than that of the grid-filter, ?. Similar to
the subgrid scale stresses, the subtest scale stresses will appear in the resulting
equations:
Tij = hatwideruiuj ?hatwideuihatwideuj. (2.9)
Assume the subtest scale stresses can be modeled using the same model as the
subgrid scale stresses in Eq. (2.6):
Tij ? 13?ijTkk = ?2Chatwide?
2
|hatwideS|hatwideSij, (2.10)
where the filtered strain rate tensor is
hatwideS
ij =
1
2
parenleftBigg
?hatwideui
?xj +
?hatwideuj
?xi
parenrightBigg
, |hatwideS| =
radicalBig
2hatwideSijhatwideSij. (2.11)
The subgrid scale stresses and the subtest scale stresses can be related to the smallest
resolved turbulent stresses, Lij, and the Germano identity [18] as follows
Lij = Tij ?hatwide?ij = hatwidestuiuj ?hatwideuihatwideuj. (2.12)
with the substitution of Eqs. (2.6) and (2.10) into Eq. (2.12), an equation with the
only unknown parameter C is obtained
Lij ? 13?ijLkk = ?2?2
bracketleftBig
(hatwide?/?)2C|hatwideS|hatwideSij ? hatwiderC|S|Sij
bracketrightBig
. (2.13)
In this study, the trapezoidal rule is used to evaluate the test-filtering operation (On
a uniform grid, the trapezoidal rule gives a sequence of weights (W?1,W0,W1) =
15
(1/4,1/2,1/4)), and consistently, the ratio between the test filter and grid filter
hatwide?/? is chosen to be ?6 [34].
Various approaches were used to solve the above equation. In general, C is
assumed to be a weak function in space and could be taken outside the test-filter
operation. One commonly used approach is the least-square method proposed by
Lilly [33]:
C(x,t) = ?12 LijMijM
ijMij
, (2.14)
where Mij is:
Mij = hatwide?
2
|hatwideS|hatwideSij ??2hatwidest|S|Sij. (2.15)
The model parameter C given by (2.14) may become locally negative, resulting in
negative eddy viscosity that are known to cause numerical instabilities. Averag-
ing C in the homogeneous directions can help alleviate this problem, but in the
present study the target applications are biological flows with complex geometries
and moving boundaries, in which a homogeneous direction is rarely available in
practical problems. For this reason, the Lagrangian averaging procedure proposed
by Meneveau et al. [36] is adopted. In this approach, the averaging is done over the
fluid particle pathlines rather than the homogeneous directions and the Lagrangian
averaging operator is defined as:
< ? >=
integraldisplay t
??
?(tprime)W(t?tprime)dtprime, (2.16)
where W(t) is an exponential weighting function chosen to give more weight to
recent times in flow history. Details of the Lagrangian dynamic model can be found
in [36].
16
If the SGS stresses in Eq. (2.3) are substituted by the model in Eq. (2.6) and
p is replaced by p+ 13?kk, then the governing equations for LES can be written as:
?ui
?t = ?
?(uiuj)
?xj ?
?p
?xi +
?
?xj
parenleftbigg
??ui?x
j
parenrightbigg
+ ??x
j
parenleftbigg
?t?uj?x
i
parenrightbigg
, (2.17)
?ui
?xi = 0, (2.18)
where ? = ?t + 1/Re is the total viscosity. Note that in the cross diffusive terms,
the constant molecular viscosity 1/Re is removed due to the continuity equation.
2.1.3 Spatial Discretization
A standard second-order central-difference scheme on a staggered grid is used
in the present study. A typical grid cell and the staggered variable arrangement is
shown in Fig. 2.1, with the pressure and scalar variables located at the center of the
grid cell and velocity components located at the cell face centers. Below the half-cell
nomenclature (such as i? 12, j ? 12, k ? 12) for velocity components is discarded to
be consistent with the numerical implementation in the code.
The finite difference schemes can be either constructed in physical space or
computational space. Usually, the construction is more convenient when the non-
uniform grids in physical space are converted to the uniform grids in computational
space. Here, the mapping between the physical space and computational space is
x = x(?), y = y(?), z = z(?), (2.19)
and
? = ?(x), ? = ?(y), ? = ?(z), (2.20)
17
where ?, ?, and ? are the coordinates in the computational space, which correspond
to x, y, and z in physical space, respectively. Therefore, the derivatives in physical
space are transformed to computational space as follows:
?
?x =
?
???x,
?
?y =
?
???y,
?
?z =
?
???z, (2.21)
with
?x = 2??x
i+1 ?xi?1
, ?y = 2??y
j+1 ?yj?1
, ?z = 2??z
k+1 ?zk?1
, (2.22)
where xi?1, xi+1, yj?1, yj+1, zk?1, and zk+1 are discrete coordinates of the Carte-
sian grid in physical space. A sketch of two-dimensional uniform grid and variable
arrangement in x ? z coordinates is shown in Fig. 2.2, in which the definition of
x, z grid lines and relative imposition of variables are given. In the computational
space, the cell sizes are set to unity, i.e., ?? = ?? = ?? = 1. By introducing this
Figure 2.1: A sketch of grid cell and the variable arrangement in Cartesian coordi-
nates.
18
transformation, the Eqs. (2.17) and (2.18) now can be written as
?ux
?t = ?
bracketleftbigg
?x?(uxux)?? +?y?(uyux)?? +?z?(uzux)??
bracketrightbigg
??x?p??
+
bracketleftbigg
?x ???
parenleftbigg
??x?ux??
parenrightbigg
+?y ???
parenleftbigg
??y?ux??
parenrightbigg
+?z ???
parenleftbigg
??z?ux??
parenrightbiggbracketrightbigg
+
bracketleftbigg
?x ???
parenleftbigg
?t?x?ux??
parenrightbigg
+?y ???
parenleftbigg
?t?x?uy??
parenrightbigg
+?z ???
parenleftbigg
?t?x?uz??
parenrightbiggbracketrightbigg
,(2.23)
?uy
?t = ?
bracketleftbigg
?x?(uxuy)?? +?y?(uyuy)?? +?z?(uzuy)??
bracketrightbigg
??y?p??
+
bracketleftbigg
?x ???
parenleftbigg
??x?uy??
parenrightbigg
+?y ???
parenleftbigg
??y?uy??
parenrightbigg
+?z ???
parenleftbigg
??z?uy??
parenrightbiggbracketrightbigg
+
bracketleftbigg
?x ???
parenleftbigg
?t?y?ux??
parenrightbigg
+?y ???
parenleftbigg
?t?y?uy??
parenrightbigg
+?z ???
parenleftbigg
?t?y?uz??
parenrightbiggbracketrightbigg
,(2.24)
?uz
?t = ?
bracketleftbigg
?x?(uxuz)?? +?y?(uyuz)?? +?z?(uzuz)??
bracketrightbigg
??z?p??
+
bracketleftbigg
?x ???
parenleftbigg
??x?uz??
parenrightbigg
+?y ???
parenleftbigg
??y?uz??
parenrightbigg
+?z ???
parenleftbigg
??z?uz??
parenrightbiggbracketrightbigg
+
bracketleftbigg
?x ???
parenleftbigg
?t?z?ux??
parenrightbigg
+?y ???
parenleftbigg
?t?z?uy??
parenrightbigg
+?z ???
parenleftbigg
?t?z?uz??
parenrightbiggbracketrightbigg
(2.25)
?x?ux?? +?y?uy?? +?z?uz?? = 0. (2.26)
The discretization of above equations on a staggered grid in computational
space is described given in detail in [3]. Here the discretization of the diagonal
convective term in the w-momentum equation is given as an example. Note that an
arithmetic average is used to obtain variables at grid locations between the points
where the variables are defined (see Fig. 2.2):
?(ww)
?z
vextendsinglevextendsingle
vextendsingle
w
i,j,k
? ?zvextendsinglevextendsinglewk ?(ww)??
vextendsinglevextendsingle
vextendsingle
w
i,j,k
= ?zvextendsinglevextendsinglewk 1??
bracketleftbiggw
i,j,k +wi,j,k+1
2 ?
wi,j,k?1 +wi,j,k
2
bracketrightbigg
. (2.27)
19
The diagonal diffusive term can be evaluated as
?
?z
parenleftbigg
??w?z
parenrightbiggvextendsinglevextendsingle
vextendsingle
w
i,j,k
? ?zvextendsinglevextendsinglewk ???
parenleftbigg
??z?w??
parenrightbiggvextendsinglevextendsingle
vextendsingle
w
i,j,k
= ?zvextendsinglevextendsinglewk 1??
bracketleftbiggparenleftbigg
??z?w??
parenrightbiggvextendsinglevextendsingle
vextendsingle
p
i,j,k+1
?
parenleftbigg
??z?w??
parenrightbiggvextendsinglevextendsingle
vextendsingle
p
i,j,k
bracketrightbigg
(2.28)
where,
parenleftbigg
??z?w??
parenrightbiggvextendsinglevextendsingle
vextendsingle
p
i,j,k+1
= ?|i,j,k +?|i,j,k+12 ?zvextendsinglevextendsinglepk+1w|i,j,k+1 ?w|i,j,k?? ,
parenleftbigg
??z?w??
parenrightbiggvextendsinglevextendsingle
vextendsingle
p
i,j,k
= ?|i,j,k?1 +?|i,j,k2 ?zvextendsinglevextendsinglepk w|i,j,k ?w|i,j,k?1?? . (2.29)
In above equations, the symbols |pi,j,k and |wi,j,k refer to quantities evaluated at p
points and w points, respectively. Hereinafter, ux ? u ? ur, uy ? v ? u?,
and uz ? w ? uz will be used interchangeably for both Cartesian and cylindrical
coordinates.
Figure 2.2: Staggered grid and variable arrangement in x?z coordinates.
20
2.1.4 Time Advancement Scheme
The fractional step method is used to integrate the governing equations in
time. In the current work, there are several choices of time advancement schemes.
For solving the equations in Cartesian coordinates, only explicit schemes are used in
all simulations, although the implicit scheme, which is discussed in the next section,
canalsobeadopted. Inparticular, thesecond-orderAdams-Bashforth(AB2)scheme
or the low-storage third-order Runge-Kutta (RK3) scheme is used:
?uki ?uk?1i
?t = ?kA
parenleftbiguk?1
i
parenrightbig+?
kA
parenleftbiguk?2
i
parenrightbig??
k
?pk?1
?xi , (2.30)
?2?k
?xi?xi =
1
?k?t
??uki
?xi , (2.31)
uki = ?uki ??k?t??
k
?xi , (2.32)
pk = pk?1 +?k, (2.33)
where k is the substep index, which ranges from 1 to 3 for Runge-Kutta scheme
and equals 1 for Adams-Bashforth scheme, ?uki is the intermediate velocity and ? is
the scalar used to project ?uki into a divergence-free space, A is the spatial operator
containing the convective, viscous and SGS terms, ?t is the time step, the RK3
coefficients are
?1 = 8/15, ?1 = 8/15, ?1 = 0;
?2 = 2/15, ?2 = 5/12, ?2 = ?17/60;
?3 = 1/3, ?3 = 3/4, ?3 = ?5/12,
(2.34)
with
3summationdisplay
k=1
?k =
3summationdisplay
k=1
(?k +?k) = 1, (2.35)
21
and the AB2 coefficients are
?1 = 1, ?1 = 3/2, ?1 = ?1/2. (2.36)
with
?1 = ?1 +?1 = 1. (2.37)
The following stability criterion (or the generalized CFL number including
the time step constraint from the viscous terms) is adopted [1]:
CFL = ?t
bracketleftbigg|u|
?x +
|v|
?y +
|w|
?z +2?
parenleftbigg 1
?x2 +
1
?y2 +
1
?z2
parenrightbiggbracketrightbigg
. (2.38)
The theoretical stability limit for the third-order Runge-Kutta scheme is ?3. Since
the cross-terms are not included in above equation, the actual CFL number in
the simulations is lower. CFL = 1.2 is used in most of the simulations in this
work. On the other hand, the theoretical stability limit for the second-order Adams-
Bashforth scheme is CFL < 1. However, the actual CFL number in the simulations
has to be much lower, e.g., CFL = 0.6 to prevent the code from blowing up.
In this study, the Adams-Bashforth scheme is used together with the fourth-order
Hamming?s predictor-corrector for fluid-structure interaction problems, in which a
constant timestep is used to maintain the accuracy of the schemes.
2.2 Transformation in Cylindrical Coordinates
Cylindrical coordinates are convenient for the description of many fluid flows
of practical interest. Also, they provide an easy way of grid adaptation for a variety
22
of external flows, e.g., insect and bird flight, fish swimming, etc., when combined
with immersed boundary formulations.
2.2.1 Governing Equations
The corresponding governing equations in cylindrical coordinates can be writ-
ten as follows:
?ur
?t +
1
r
?(rurur)
?r +
1
r
?(u?ur)
?? +
?(uzur)
?z ?
u2?
r
= ??p?r + 1Re
bracketleftbigg ?
?r
parenleftbigg1
r
?
?r(rur)
parenrightbigg
+ 1r2 ?
2ur
??2 +
?2ur
?z2 ?
2
r2
?u?
??
bracketrightbigg
, (2.39)
?u?
?t +
1
r
?(ruru?)
?r +
1
r
?(u?u?)
?? +
?(uzu?)
?z +
uru?
r
= ?1r ?p?? + 1Re
bracketleftbigg1
r
?
?r
parenleftbigg
r?u??r
parenrightbigg
? u?r2 + 1r2 ?
2u?
??2 +
?2u?
?z2 +
2
r2
?ur
??
bracketrightbigg
, (2.40)
?uz
?t +
1
r
?(ruruz)
?r +
1
r
?(u?uz)
?? +
?(uzuz)
?z
= ??p?z + 1Re
bracketleftbigg1
r
?
?r
parenleftbigg
r?uz?r
parenrightbigg
+ 1r2 ?
2uz
??2 +
?2uz
?z2
bracketrightbigg
, (2.41)
1
r
?(rur)
?r +
1
r
?u?
?? +
?uz
?z = 0, (2.42)
where the subscript r represents the radial direction, ? represents the azimuthal
direction, and z represents the axial direction (or streamwise direction in this study).
As with the case of Cartesian coordinates a filtering operation can be applied
to the above equations to define the resolved and unresolved scales. The result-
ing equations governing the evolution of the large scales for incompressible flow in
23
cylindrical coordinates can be written as follows:
?ur
?t = ?
bracketleftbigg1
r
?(rurur)
?r +
1
r
?(u?ur)
?? +
?(uzur)
?z
bracketrightbigg
? ?p?r
+ ??r
parenleftbigg
?1r ??r(rur)
parenrightbigg
+ 1r2 ???
parenleftbigg
??ur??
parenrightbigg
+ ??z
parenleftbigg
??ur?z
parenrightbigg
+ ??r
parenleftbigg
?t1r ??r(rur)
parenrightbigg
+ 1r ???
parenleftbigg
?t?u??r
parenrightbigg
+ ??z
parenleftbigg
?t?uz?r
parenrightbigg
+ u
2
?
r
? ? 2r2 ?u??? ??t 1r2 ?u??? ? u?r2 ??t?? ? 2urr ??t?r , (2.43)
?u?
?t = ?
bracketleftbigg1
r
?(ruru?)
?r +
1
r
?(u?u?)
?? +
?(uzu?)
?z
bracketrightbigg
? 1r ?p??
+ 1r ??r
parenleftbigg
?r?u??r
parenrightbigg
+ 1r2 ???
parenleftbigg
??u???
parenrightbigg
+ ??z
parenleftbigg
??u??z
parenrightbigg
+ 1r ??r
parenleftbigg
?t?ur??
parenrightbigg
+ 1r2 ???
parenleftbigg
?t?u???
parenrightbigg
+ 1r ??z
parenleftbigg
?t?uz??
parenrightbigg
? uru?r
+ ? 2r2 ?ur?? +?t 1r2 ?ur?? ??u?r2 + 2urr2 ??t?? ? u?r ??t?r , (2.44)
?uz
?t = ?
bracketleftbigg1
r
?(ruruz)
?r +
1
r
?(u?uz)
?? +
?(uzuz)
?z
bracketrightbigg
? ?p?z
+ 1r ??r
parenleftbigg
?r?uz?r
parenrightbigg
+ 1r2 ???
parenleftbigg
??uz??
parenrightbigg
+ ??z
parenleftbigg
??uz?z
parenrightbigg
+ 1r ??r
parenleftbigg
?tr?ur?z
parenrightbigg
+ 1r ???
parenleftbigg
?t?u??z
parenrightbigg
+ ??z
parenleftbigg
?t?uz?z
parenrightbigg
, (2.45)
1
r
?(rur)
?r +
1
r
?u?
?? +
?uz
?z = 0, (2.46)
In the above equations ?t is the SGS eddy viscosity as in the Cartesian formulation.
In the Smagorinsky model for example:
?ij ? 13?ij?kk = ?2?tSij = ?2C?2|S|Sij (2.47)
24
with i,j,k = r,?,z and ? = (r?r???z)1/3. ?r, ??, and ?z the local grid size in
r, ?, and z coordinate direction, respectively.
The strain rate tensor in cylindrical coordinates can be written as follows:
Srr = ?ur?r ,
Sr? = 12
parenleftbigg?u
?
?r +
?ur
?? ?
u?
r
parenrightbigg
,
S?? =
parenleftbigg1
r
?u?
?? +
ur
r
parenrightbigg
,
S?z = 12
parenleftbigg?u
?
?z +
1
r
?uz
??
parenrightbigg
,
Szr = 12
parenleftbigg?u
r
?z +
?uz
?r
parenrightbigg
,
Szz = ?uz?z . (2.48)
Note that the equations above in cylindrical coordinates and the equations in
Cartesian coordinates have many common terms. This fact was utilized to formu-
late a generalized code for both cylindrical and Cartesian coordinates. The major
differences between the equations in two coordinates systems are: a) r and 1/r
terms in cylindrical coordinates; and b) additional terms due to coordinates curva-
ture in cylindrical coordinates. In the current study, a flag is used in the input file
to identify the reference frames. When this flag is switched on then all terms are
calculated and r is used as defined in cylindrical coordinates. However, when this
flag is switched off, then r is set to constant 1 and the additional terms (lines 4 and
5 in Eqs. (2.43) and (2.44)) are ignored in the calculation.
25
2.2.2 Spatial Discretization
A staggered grid is also used in the case of cylindrical coordinates. The variable
arrangement for a typical grid cell and grid cell near the axis is shown in Fig. 2.3.
Fig. 2.4 shows the staggered grid in r ?? plane. It is evident that only the radial
velocity component is located at the centerline. The grid is uniform in the azimuthal
direction, but non-uniform in the radial and axial (streamwise) directions. Here, the
mapping between the physical space and computational space is derived as:
r = r(?), ? = ?(?), z = z(?), (2.49)
and
? = ?(r), ? = ?(?), ? = ?(z), (2.50)
where ?, ?, and ? are the coordinates in the computational space, which correspond
to r, ?, and z in physical space, respectively. Therefore, the derivatives in physical
space are transformed to computational space as follows:
?
?r =
?
???r,
?
?? =
?
????,
?
?z =
?
???z, (2.51)
where
?r = 2??r
i+1 ?ri?1
, ?? = 2???
j+1 ??j?1
, ?z = 2??z
k+1 ?zk?1
. (2.52)
The cell sizes in the computational space are set to unity, i.e., ?? = ?? = ?? = 1.
By introducing this transformation, Eqs. (2.17) and (2.18) now can be written
26
(a) (b)
Figure 2.3: A sketch of grid cell and the variable arrangement in cylindrical coordi-
nates. a) A typical grid cell; b) a grid cell near the centerline.
Figure 2.4: Staggered grid and the variable arrangement in r?? plane. ? p, uz; square
ur; triangle u?.
27
as:
?ur
?t = ?
bracketleftbigg1
r?r
?(rurur)
?? +
1
r??
?(u?ur)
?? +?z
?(uzur)
??
bracketrightbigg
??r?p??
+
bracketleftbigg
?r ???
parenleftbigg1
r??r
?(rur)
??
parenrightbigg
+ 1r2?? ???
parenleftbigg
????ur??
parenrightbigg
+?z ???
parenleftbigg
??z?ur??
parenrightbiggbracketrightbigg
+
bracketleftbigg
?r ???
parenleftbigg1
r?t?r
?(rur)
??
parenrightbigg
+ 1r?? ???
parenleftbigg
?t?r?u???
parenrightbigg
+?z ???
parenleftbigg
?t?r?uz??
parenrightbiggbracketrightbigg
+ u
2
?
r
? (2? +?t) 1r2???u??? ? u?r2????t?? ? 2urr ?r??t?? , (2.53)
?u?
?t = ?
bracketleftbigg1
r?r
?(ruru?)
?? +
1
r??
?(u?u?)
?? +?z
?(uzu?)
??
bracketrightbigg
? 1r???p??
+
bracketleftbigg1
r?r
?
??
parenleftbigg
?r?r?u???
parenrightbigg
+ 1r2?? ???
parenleftbigg
????u???
parenrightbigg
+?z ???
parenleftbigg
??z?u???
parenrightbiggbracketrightbigg
+
bracketleftbigg1
r?r
?
??
parenleftbigg
?t???ur??
parenrightbigg
+ 1r2?? ???
parenleftbigg
?t???u???
parenrightbigg
+ 1r?z ???
parenleftbigg
?t???uz??
parenrightbiggbracketrightbigg
? uru?r
+ (2? +?t) 1r2???ur?? ??u?r2 + 2urr2 ????t?? ? u?r ?r??t?? , (2.54)
?uz
?t = ?
bracketleftbigg1
r?r
?(ruruz)
?? +
1
r??
?(u?uz)
?? +?z
?(uzuz)
??
bracketrightbigg
??z?p??
+
bracketleftbigg1
r?r
?
??
parenleftbigg
?r?r?uz??
parenrightbigg
+ 1r2?? ???
parenleftbigg
????uz??
parenrightbigg
+?z ???
parenleftbigg
??z?uz??
parenrightbiggbracketrightbigg
+
bracketleftbigg1
r?r
?
??
parenleftbigg
?t?z?ur??
parenrightbigg
+ 1r?? ???
parenleftbigg
?t?z?u???
parenrightbigg
+?z ???
parenleftbigg
?t?z?uz??
parenrightbiggbracketrightbigg
,(2.55)
1
r?r
?(rur)
?? +
1
r??
?u?
?? +?z
?uz
?? = 0. (2.56)
The discretization of all derivatives on the staggered grid in the computational space
is similar to that in Cartesian coordinates that has described in th previous section.
Also in this case, the arithmetic average is used to obtain variables at grid locations
between the points where the variables are defined.
28
A major complication of the cylindrical coordinate formulation is the mathe-
matical singularity at the centerline, which is not present in Cartesian coordinates.
Due the singularity at the axis, both the radial and the azimuthal velocity
components cannot be defined at the axis even when the flow field itself does not
have a singularity problem at the center line. Over the past years ways of addressing
the singularity issue at the centerline with various degrees of complexity have been
proposed[14,1,66,38]. Heretheapproachproposedin[46]isadopted. Basically, the
radial and the azimuthal velocity components are assumed to have multiple values
at the axis and a linear averaging of two symmetrical points over the axis provides
the value at the center line. Notice that for the radial velocity component this is an
average over two grid cells and has a decreased order of accuracy. Therefore, finer
grid spacing near the axis is required to obtain good resolution of the flow field near
the axis.
Fig. 2.5 shows the current centerline treatment. In particular, ghost cells are
used to implement the centerline boundary conditions. In Fig. 2.5a the required
boundary conditions on the ghost cell for solving the Navier-Stokes equations are
given. It is obvious from Fig. 2.5b, which shows the variable collocation on the grid
cell across the centerline, the variables on the ghost cell in Fig. 2.5a can be obtained
29
(a) (b)
Figure 2.5: Centerline treatment of the cylindrical coordinates. a) The required
boundary conditions and variable collacation on the ghost cell; b) Variables defined
across the centerline on r?? plane.
30
from the grid cell across the centerline in Fig. 2.5b:
u?(??r/2,?,z) = ?u?(?r/2,? +pi,z) or v1,j,k = ?v2,Ny/2+j,k,
uz(??r/2,?,z) = uz(?r/2,? +pi,z) or w1,j,k = w2,Ny/2+j,k,
p(??r/2,?,z) = p(?r/2,? +pi,z) or p1,j,k = p2,Ny/2+j,k,
?t(??r/2,?,z) = ?t(?r/2,? +pi,z) or ?t 1,j,k = ?t 2,Ny/2+j,k. (2.57)
where the boundary conditions for p are not directly utilized in the current formu-
lation.
Inspection of Eqs. (2.53)?(2.56) to locate the potential problematic terms
due to the 1/r factors shows that many terms turn out to be free of the singularity
problem because of the r factors in front of the fluxes or derivatives (with the
assumption that those fluxes and derivatives are finite values, which are usually
sound in physical sense). However, a centerline value for ur is still required for those
terms in red color. Notice that we write the diffusive terms in r direction in the
forms appeared in Eq. (2.53) in order to be able to use Crank-Nicolson scheme in
radial direction [66]. Here, the centerline bounary condition for ur is obtained by
averaging values of the opposing ur across the centerline:
ur(r = 0,?,z) = 12 [ur(?r,?,z)?ur(?r,? +pi,z)] (2.58)
or
u1,j,k = 12parenleftbigu2,j,k ?u2,Ny/2+j,kparenrightbig. (2.59)
The above approximation is a linear interpolation over a distance of 2?r. Therefore,
the grid needs to maintain good quality near the centerline to avoid large discretiza-
31
tion error. For all simulations in this work, this boundary condition is used and no
obvious problem in the solution was found due to this centerline treatment.
2.2.3 Time Advancement Schemes
The fractional step method, which is discussed in the previous section, is
used to integrate the governing equations in cylindrical coordinates in time. The
extremely small grid size in azimuthal direction near the centerline imposes severe
timestep constraints if all terms are treated explicitly. To remove this limitation,
the diffusive terms in the azimuthal direction are treated implicitly. In particular, a
second-order Adams-Bashforth (AB2) scheme or the low-storage third-order Runge-
Kutta(RK3)schemeisusedfortermstreatedexplicitly, andthesecond-orderCrank-
Nicholson (CN2) scheme is used for terms treated implicitly:
bracketleftbigg
1? ?k?t2 1r2 ???
parenleftbigg
? ???
parenrightbiggbracketrightbigg
?ukr = RHSk?1r
= uk?1r +?k?tAk?1r +?k?tAk?2r + ?k?t2 Bk?1r ??k?t?p
k?1
?r , (2.60)bracketleftbigg
1? ?k?t2 1r2 ???
parenleftbigg
(? +?t) ???
parenrightbiggbracketrightbigg
?uk? = RHSk?1?
= uk?1? +?k?tAk?1? +?k?tAk?2? + ?k?t2 Bk?1? ??k?t1r ?p
k?1
?? , (2.61)bracketleftbigg
1? ?k?t2 1r2 ???
parenleftbigg
? ???
parenrightbiggbracketrightbigg
?ukz = RHSk?1z
= uk?1z +?k?tAk?1z +?k?tAk?2z + ?k?t2 Bk?1z ??k?t?p
k?1
?z . (2.62)
32
where the coefficients are the same as in the previous section, and the operators A
and B are defined as:
Ar = ?
bracketleftbigg1
r
?(rurur)
?r +
1
r
?(u?ur)
?? +
?(uzur)
?z
bracketrightbigg
+ ??r
bracketleftbigg
(? +?t)1r ?(rur)?r
bracketrightbigg
+ ??z
parenleftbigg
??ur?z
parenrightbigg
+ 1r ???
parenleftbigg
?t?u??r
parenrightbigg
+ ??z
parenleftbigg
?t?uz?r
parenrightbigg
+ u
2
?
r ?(2? +?t)
1
r2
?u?
?? ?
u?
r2
??t
?? ?
2ur
r
??t
?r , (2.63)
A? = ?
bracketleftbigg1
r
?(ruru?)
?r +
1
r
?(u?u?)
?? +
?(uzu?)
?z
bracketrightbigg
+ 1r ??r
parenleftbigg
?r?u??r
parenrightbigg
+ ??z
parenleftbigg
??u??z
parenrightbigg
+ 1r ??r
parenleftbigg
?t?ur??
parenrightbigg
+ 1r ??z
parenleftbigg
?t?uz??
parenrightbigg
? uru?r +(2? +?t) 1r2 ?ur?? ??u?r2 + 2urr2 ??t?? ? u?r ??t?r , (2.64)
Az = ?
bracketleftbigg1
r
?(ruruz)
?r +
1
r
?(u?uz)
?? +
?(uzuz)
?z
bracketrightbigg
+ 1r ??r
parenleftbigg
?r?uz?r
parenrightbigg
+ ??z
bracketleftbigg
(? +?t)?uz?z
bracketrightbigg
+ 1r ??r
parenleftbigg
?tr?ur?z
parenrightbigg
+ 1r ???
parenleftbigg
?t?u??z
parenrightbigg
, (2.65)
Br = 1r2 ???
parenleftbigg
??ur??
parenrightbigg
, (2.66)
B? = 1r2 ???
bracketleftbigg
(? +?t) ?u???
bracketrightbigg
, (2.67)
Bz = 1r2 ???
parenleftbigg
??uz??
parenrightbigg
. (2.68)
The spatial discretization of above equations results in a series of cyclic tridi-
agonal equations, which are solved using the solver for cyclic tridiagonal systems
from Numerical Recipes [51]. No boundary condition is required for the variables at
the intermediate sub-steps due to the periodicity in azimuthal direction.
The following stability criterion from [1] is adopted here
CFL = ?t
bracketleftbigg|u
r|
?r +
|u?|
r?? +
|uz|
?z +4?
parenleftbigg 1
?r2 +
1
?z2
parenrightbiggbracketrightbigg
. (2.69)
33
In the simulations of flow past a sphere, CFL = 1.2 is used.
2.3 Poisson Equation
The Poisson equation Eq. (2.31) can be rewritten here in vector form as
?2? = f = 1?
n?t
?? ?uni , where ? = pn ?pn?1, (2.70)
where n is sub step index for RK3 scheme or time step index for AB2 scheme. The
discrete form of the above equation in Cartesian coordinates is
parenleftbigg ?2
?2x +
?2
?2y +
?2
?2z
parenrightbigg
?i,j,k = fi,j,k, (2.71)
and in cylindrical coordinates is
bracketleftbigg1
r
?
?r
parenleftbigg
r ??r
parenrightbigg
+ 1r2 ?
2
?2? +
?2
?2z
bracketrightbigg
?i,j,k = fi,j,k. (2.72)
The discrete operators above have to be consistent with the ones used in the momen-
tum equations. Otherwise, mass conservation cannot be guaranteed. Here only the
solution procedure in cylindrical coordinates will be discussed because it is evident
that the equation for Cartesian coordinates will be recovered when r in Eq. (2.72)
is set to unity. The discretized form of Eq. (2.72) can be written as:
1
rp|i?x|
p
i
1
??2
braceleftbig[r
u|i?x|ui (?i+1,j,k ??i,j,k)]?
bracketleftbigr
u|i?1?x|ui?1 (?i,j,k ??i?1,j,k)
bracketrightbigbracerightbig
+ 1r2
p|i
?y|pj 1??2 braceleftbigbracketleftbig?y|vj (?i,j+1,k ??i,j,k)bracketrightbig?bracketleftbig?y|vj?1 (?i,j,k ??i,j?1,k)bracketrightbigbracerightbig
+ ?z|pk 1??2 braceleftbig[?z|wk (?i,j,k+1 ??i,j,k)]?bracketleftbig?z|wk?1 (?i,j,k ??i,j,k?1)bracketrightbigbracerightbig
= fi,j,k, (2.73)
34
or
ami ?i?1,j,k +bmi ?i,j,k +cmi ?i+1,j,k
+ alj ?i,j?1,k +blj ?i,j,k +clj ?i,j+1,k
+ ank ?i,j,k?1 +bnk ?i,j,k +cnk ?i,j,k+1 = fi,j,k, (2.74)
with the coefficients
ami = 1??2 1r
p|i
?x|piru|i?1?x|ui?1,
cmi = 1??2 1r
p|i
?x|piru|i?x|ui ,
bmi = ?ami ?cmi,
alj = 1??2 1r2
p|i
?y|pj?y|vj?1,
clj = 1??2 1r2
p|i
?y|pj?y|vj,
blj = ?alj ?clj,
ank = 1??2?z|pk?z|wk?1,
cnk = 1??2?z|pk?z|wk ,
bnk = ?ank ?cnk. (2.75)
The above equation is solved using a combination of a fast Fourier transform
(FFT) method from FFTPACK [58] and a direct solution procedure from FISH-
PACK [57]. The major advantage of this direct solver is that it can converge the
solution to machine accuracy with only one iteration, which is a great saving com-
paring with the iterative solvers. However, to be able to utilize Fourier transforms,
the computational grid must be uniform in the direction in which the FFT is per-
formed. In this work, the grid is uniform in the spanwise direction, i.e., y or ?
35
direction. Therefore, the fast Fourier transform in the spanwise direction can trans-
form Eq. 2.70 into a set of two-dimensional Helmholtz equations in the uncoupled
wave number space:
bracketleftbigg1
r
?
?r
parenleftbigg
r ??r
parenrightbigg
+ ?
2
?2z +
1
r2k
prime
l
bracketrightbigg
??i,j,l = ?fi,l,k, (2.76)
and the modified wave number kprimel is defined as
kprimel = 2??2
bracketleftbigg
1?cos
parenleftbigg2pil
N?
parenrightbiggbracketrightbigg
(2.77)
where l is the wave number, N? is the number of grid cells (not including ghost cells)
and ?? is the cell size in the spanwise direction. Eq. (2.74) can be rewritten as
ami ??i?1,l,k +
parenleftbigg
bmi ? k
prime
l
r2p|i
parenrightbigg
??i,l,k +cmi ??i+1,l,k
+ ank ??i,l,k?1 +bnk ??i,l,k +cnk ??i,l,k+1 = ?fi,l,k. (2.78)
The above equations can be solved separately for each wave number. Both the
real and imaginary part of each wave number are solved using the ?BLKTRI? rou-
tine, in which a generalized cyclic reduction algorithm [57] is implemented, in the
FISHPACK library.
In cases where periodic boundary conditions are used in both streamwise and
spanwise directions, an FFT is also applied in the streamwise direction, and the
resulting series of (cyclic) tridiagonal equations are solved using the (cyclic) tridi-
agonal system solver from Numerical Recipes [51].
36
Figure 2.6: The implementation of Dirichlet and Neumann boundary condition.
2.4 Boundary Conditions
In this work, ghost cells are used to implement all boundary conditions. The
advantages of the ghost cell approach are that the grid spacing near the bound-
ary is continuous and the code can have a generalized form for all grid points at
which the equations are solved. In addition, it provides a natural way to implement
parallelization via domain decomposition technique.
2.4.1 Dirichlet and Neumann Boundary Conditions
TheimplementationofDirichletandNeumannboundaryconditionsisisshown
in Fig. 2.6, in which the lower-left corner of a ? ? ? plane of the computational
domain and the collocation of the variables is given. For velocity component normal
37
to the wall, e.g., u to the lower boundary and w to the left boundary in the figure,
the Dirichlet condition can be directly enforced as
u1,j,k = ub, and wi,j,1 = wb, (2.79)
where ub and wb are the prescribed normal velocity components on the lower wall
and left wall, respectively.
The non-slip conditions for the wall-tangential components are implemented
through the use of ghost cells:
vt = 12 (v1,j,k +v2,j,k) or v1,j,k = 2vt ?v2,j,k,
wt = 12 (w1,j,k +w2,j,k) or w1,j,k = 2wt ?w2,j,k, (2.80)
where vt and wt are prescribed tangential components in y and z direction, respec-
tively. For stationary wall without tangential movement, vt = 0 and wt = 0.
The Neumann boundary condition for an arbitrary variable ? can be written
as
??
?n = f, (2.81)
where n is the normal direction of the boundary, and f is a known function. Refer-
ring to the lower boundary in Fig. 2.6, this condition can be implemented with the
help of ghost cells as follows:
?1,j,k = ?2,j,k ?f?x, (2.82)
or
?1,j,k = ?2,j,k (2.83)
for homogeneous boundary condition.
38
2.4.2 Convective Boundary Condition
The convective boundary condition proposed by Orlanski [41] is used for out-
flow boundaries in this study. This condition has been found very successful in
convecting structures out of the domain without distorting the flow in the compu-
tational domain. In such a case the boundary velocity can be obtained from the
following equation:
?ui
?t +Uconv
?ui
?z = 0, (2.84)
where ui is any velocity component, and Uconv is the convective velocity, which is
set to the mean streamwise velocity at the exit plane. Here, the outflow boundary is
always normal to the streamwise direction, which is z direction for both Cartesian
and cylindrical coordinates. However, the corresponding term from the continuity
equation has to be used instead of ?ui/?z when the exit plane is normal to the
radial or azimuthal direction in cylindrical coordinates.
Eq. (2.84) is discretized using explicit Euler scheme in time, one-sided dif-
ference formulas for the streamwise velocity component and central difference for
the other two velocity components in space. Note that the predicted streamwise
velocity component is adjusted every timestep to globally conserve mass.
2.4.3 Periodic Boundary Conditions
The implementation of periodic boundary condition through ghost cells is also
very simple. As shown in Fig. 2.7, the solution values on the left side are directly
39
Figure 2.7: The implementation of periodic boundary condition. Arrows show the
directions of data duplication.
copied to the ghost cells on the right side and vise versa:
ui,j,Nz+2 = ui,j,2, vi,j,Nz+2 = vi,j,2, wi,j,Nz+2 = wi,j,2, pi,j,Nz+2 = pi,j,2;
ui,j,1 = ui,j,Nz+1, vi,j,1 = vi,j,Nz+1, wi,j,1 = wi,j,Nz+1, pi,j,1 = pi,j,Nz+1. (2.85)
This way periodicity of the function and its first derivative is imposed.
2.5 Validation of Basic Solver
The basic solver for the Navier-Stokes equations in Cartesian coordinates has
been validated extensively [3, 4, 5]. Here the basic solver in cylindrical coordinates
is validated for the case of fully developed turbulent pipe flow. This problem has not
been studied numerically as extensively as its Cartesian equivalent, the turbulent
channel flow. Although the geometry is equally simple, the presence of a singularity
at the axis gives rise to numerical difficulties, making it a more challenging problem.
40
In the following we will present LES of this problem and will compare our results
with the DNS data reported in [16].
As already discussed, a semi-implicit scheme (third-order Runge-Kutta scheme
for explicit terms and second-order Crank-Nicolson scheme for implicit terms) is
used for the time advancement. To maximize the timesteps the following procedure
proposed in [1] is adapted in this study. Basically, the computational domain is
divided into two separate regions: the core-region, which includes the axis and
extends to a given radius (usually in the middle of the radius R), and the outer-
region, which includes all the rest. Fig. 2.8 shows how this dividing is performed
for a typical computational grid. Then all the convection and diffusion terms in the
azimuthal direction are treated implicitly in the core-region, while all other terms
are treated explicitly. The azimuthal momentum equation is linearized and solved
first to provide a predicted velocity component for the other two equations. On the
other hand, all the convection and diffusion terms in the radial direction are treated
implicitly in the outer-region with all other terms treated explicitly. Again the
radial momentum equation is linearized and solved before the other two equations
as the predicted radial velocity is needed by those two equations. By doing this, the
time step can be twenty times larger than that with only the diffusion terms in the
azimuthal direction treated implicitly. Both the plane-averaged and the Lagrangian
dynamic SGS model are implemented in cylindrical coordinates in this study.
The computational domain is the same as in [14] and [1]: 0 ? z ? Lz = 10R,
0 ? r ? R, and 0 ? ? ? 2pi with R the radius of the pipe. Periodic boundary
conditions are applied in the streamwise and azimuthal directions. Computations
41
Figure 2.8: Typical computational grid. Shaded area is the core-region; un-shaded
area is the outer-region.
on three grids, 32?64?48, 48?64?64, and 64?96?64 (Nr?N??Nz), have been
performed. Here only the results on grid 1 and 3 will be reported as the results on
grid 2 are always in the middle of the other two grids. The grid spacing is uniform
in the streamwise and azimuthal directions, and stretched in the radial direction to
cluster points near the wall. The Reynolds number based on the friction velocity u?
and the pipe radius R is 180.
Table 2.1 summarizes the grid resolution in wall units for the three different
grids used in this works. The corresponding resolution for the reference DNS and
LES has been added for comparison. Some mean flow properties together with the
data from [1] are listed in Table 2.2. All data have been normalized using the friction
velocity u? and the radius R. In the table Ucl is the center line velocity and Ubulk
is the bulk velocity. Cf is the friction coefficient based on the bulk velocity and
42
Table 2.1: Grid resolution in wall units compared with data in [1]
.
Nr ?N? ?Nz ?z+ (R??)+ (?r??)+ ?r+wall ?r+cl ?r+max
Grid 1: 32?64?48 37.5 17.7 0.85 1.54 8.62 8.62
Grid 3: 64?96?64 28.13 11.8 0.29 0.39 4.39 4.39
68?128?256 DNS [1] 7.03 8.84 0.13 0.17 2.61 5.94
38?64?32 LES [1] 56.3 17.7 0.42 0.69 4.28 8.87
Table 2.2: Mean flow properties compared with data in [1]
.
Nr ?N? ?Nz Ucl Ubulk Ucl/Ubulk Cf ?103
32?64?48 LES 19.27 15.15 1.27 8.50
64?96?64 LES 19.00 14.76 1.29 9.05
68?128?256 DNS [1] 19.32 14.70 1.31 9.25
38?64?32 LES [1] 19.15 15.11 1.27 8.76
calculated according to Cf = ?w/parenleftbig12?U2bulkparenrightbig. From Table 2.2 it is evident that our
flows match closely those in [1].
The mean streamwise velocity profile in wall units is shown in Fig. 2.9. Results
from grid 3 agree well with the DNS data. As expected the LES on grid 1 results
in a high intercept of the log-law.
Radial, azimuthal, and axial turbulence intensities are shown in Figs. 2.10, 2.11,
and 2.12, respectively (only the resolved part of the intensities plotted here). The
43
Figure 2.9: Mean streamwise velocity profile in wall unit. DNS data from [16]
overall agreement between our simulations the DNS data is very good. However, the
intensities near the center line are under-predicted due to the coarse grid spacing
there. Figure 2.13 shows the resolved turbulent shear stress and the total stress
balance in the pipe. The viscous shear stress, resolved shear stress, and the SGS
shear stress are summed up to obtained the total shear stress. The results from
both grids shown here give total stress very close to the DNS data. The finer grid
performs better.
Figure 2.14 shows the ratio of the eddy viscosity to molecular viscosity as a
function of radius. As expected, fine grid gives lower eddy viscosity. Due to the
coarse grid spacing near the center line, there is a bump of the eddy viscosity from
the plane-averaged model near the center line. Also, it shows a sharp decrease to
44
Figure 2.10: Radial turbulence intensity. DNS data from [16]
Figure 2.11: Azimuthal turbulence intensity. DNS data from [16].
45
Figure 2.12: Axial turbulence intensity. DNS data from [16].
Figure 2.13: Total stress balance. DNS data from [16].
46
Figure 2.14: Ratio of eddy viscosity to molecular viscosity.
(a) (b)
Figure 2.15: Instantaneous field in a polar plane. (a) Axial velocity; (b) Axial
vorticity.
47
Figure 2.16: Instantaneous vortical structure by the Q-criterion (colored by axial
vorticity).
the center line due to the decreased strain-rate as we average the axial velocity
component over the center line.
Fig. 2.15a shows an instantaneous snapshot of the streamwise velocity in
a polar plane and Fig. 2.15b shows an instantaneous snapshot of the streamwise
vorticity in the same plane. The latter gives a very good illustration of the near-wall
flow structure. A snapshot of the instantaneous vortical structure visualized by the
Q-criterion (For details on the identification of coherent structures using the ?Q-
criterion? the reader is referred to [25]) is shown in Fig. 2.16. The quasi-streamwise
structures that are responsible for most turbulent production can be observed.
48
Chapter 3
Embedded Boundary Method
In this chapter, the embedded-boundary formulation is introduced to extend
the Cartesian/cylindrical grid solver to cases with complex moving boundaries. In
Sec. 3.1, establishment of the interface-grid relation is discussed. Then the treatment
of stationary immersed boundaries and moving immersed boundaries is given in
Sections 3.2 and 3.3, respectively. The surface force calculation procedure, which
computes the hydrodynamic loads to the structure, will be introduced in Sec. 3.4.
Then the fluid-structure coupling scheme is discussed in Sec. 3.5. To utilize the
current high-performance parallel computing platforms, the above algorithms have
been parallelized. The details are given in Sec. 3.6.
3.1 Interface-grid Relation
3.1.1 Interface Description
There are two approaches for interface description: Eulerian and Lagrangian
methods. In Eulerian methods, such as level set [42], the interface is implicitly
given by a field function, namely, the signed shortest distance to the interface.
In the Lagrangian methods on the other hand, such as front tracking and marker
particles, the interface is explicitly described independent of the underlying grid. In
the present study a front tracking scheme originally proposed for multi-phase flow
49
problems is used for the description of the fluid/solid interface [62].
Figure 3.1: The parametrized description of interfaces of arbitrary shapes using
marker particles.
As shown in Fig. 3.1, a two-dimensional immersed interface, ?, can be rep-
resented by a series of interfacial marker particles, which are defined by arc length
coordinates vectorX(s,t). The immersed interface can have arbitrary shapes and it can
be open or closed. The marker particles are evenly attached to the interface with
a spacing approximating the local grid size, and the beginning of the arclength co-
ordinate is defined such that the fluid (or interested side of the interface) is always
to the left of the observer as one moves along the interface toward increasing s.
For each marker particle with arclength coordinates vectorXi, the functions defining the
coordinates can be written as
x(s,t) = axs2 +bxs+cx and y(s,t) = ays2 +bys+cy. (3.1)
These functions are generated at each sub-step of the splitting scheme for moving
interfaces. The coefficients ax,y, bx,y, cx,y can be obtained by fitting quadratic
50
polynomials to particle (i) and its two neighbors (i?1) and (i+1).
The normal from any location on the interface to the fluid can be calculated
by the following equations:
nx = ?ysradicalbig(x2
s +y2s)
, ny = xsradicalbig(x2
s +y2s)
, (3.2)
where the derivatives, xs, ys can be evaluated from the functions in Eq. (3.1) above
as follows,
xs(s,t) = 2axs+bx, ys(s,t) = 2ays+by. (3.3)
The coordinates vectorX(s,t) and the coefficients are stored for each marker particle on
the interface. For three-dimensional interfaces Bi-spline fitting can be used.
3.1.2 Tagging of Points on the Eulerian Grid
Having defined the immersed interface as a series of marker particles, one can
now establish the relationship between these particles and the underlying Eulerian
grid. Fig. 3.2a shows the parametrized interface immersed in a Cartesian grid. Here
an approach similar to that presented in [4] is adopted for the tagging process and
is summarized as follows:
1. Determine the part of the grid occupied by the interface as the coordinates of
each marker particle are known from the previous Section.
2. For those grid points, a search for the closest marker particle, sb, is performed.
A ray, r, is shot from this marker particle to the grid point, and the dot product
of r and nb is calculated.
51
(a) (b)
Figure 3.2: Grid?interface relation. (a) Parametrized interface immersed in the un-
derlyingCartesiangrid; (b)Zoominthevicinityofinterfacewheretheinside/outside
status of the Eulerian grid points are shown. square Fluid points; squaresolid Solid points.
3. If r?nb < 0, then this grid point is inside the interface and assigned a tag of
?1; otherwise, it is outside the interface and maintains its initial tag of 1.
Fig. 3.2bshowsazoom-inofthevicinityoftheinterfacewiththeinside/outside
status of the Eulerian grid points shown. In the next step the boundary points are
identified, which are points in the fluid phase with at least one neighboring point
in the solid. These points will be central to the reconstruction procedure described
in the next section. An example of the result of the flagging process is shown in
Fig. 3.3a where all Eulerian grid points are split into three different categories: (a)
forcing points, which are grid points in the fluid phase that have one or more neigh-
boring points in the solid phase; (b) solid points, which are all the points in the
52
solid phase; (c) fluid points, which are all the remaining points in the fluid phase.
In the solution procedure, the fluid points are the unknowns, the forcing points are
boundary points, while the solid points do not influence the rest of the computation.
For a stationary boundary the above tagging and flagging process is done only once
at the beginning of the computation.
For a moving body the process is repeated at each timestep. In addition,
another set of flags is used for the field extension treatment, which is shown in
Fig. 3.3b. Again, all Eulerian grid points are split into three different categories:
(a) pseudo-fluid points, which are grid points in the solid phase that have one or
more neighboring points in the fluid phase; (b) solid points, which are all other
remaining points in the solid phase; (c) fluid points, which are all the grid points
in the fluid phase. In the field extension procedure, the solution at those pseudo-
fluid points are extrapolated from the known solution of the fluid points and the
interface. It is obvious that the forcing points and the pseudo-fluid points are the
points closest to the interface inside the fluid and inside the solid, respectively. The
proper manipulation of those boundary points is very important to the success of
the embedded boundary formulation.
3.1.3 Establishment of Interface-Normal Intersections
With the boundary points identified, the next task is to establish the infor-
mation required for the reconstruction procedure. Central to this algorithm is the
normal from the boundary points to the interface. This normal passes through
53
(a) Momentum Forcing (b) Field Extension
Figure 3.3: Identification of boundary points. (a) triangle Forcing points, square fluid points,
and squaresolid solid points for momentum forcing procedure; (b) trianglesolid Pseudo-fluid points, square
fluid points, and squaresolid solid points for field extension procedure.
boundary point (xi,yj) and intersects the interface at sn, or (xn,yn). Therefore, the
unit normal vector of this line is identical to the unit normal vector of the interface
on point (xn,yn) and can be written as,
xi ?xnradicalbig
(xi ?xn)2 +(yj ?yn)2 = nx =
?ysradicalbig
(x2s +y2s),
yj ?ynradicalbig
(xi ?xn)2 +(yj ?yn)2 = ny =
xsradicalbig
(x2s +y2s). (3.4)
Obviously, from the above equations it can be deduced that
(xi ?xn) xs +(yj ?yn) ys = 0, (3.5)
and, substituting Eq. (3.3) into Eq. (3.5), it gives
(xi ?xn) (2axsn +bx)+(yj ?yn) (2aysn +by) = 0, (3.6)
54
Figure 3.4: Schematic of the solution procedure for interface-normal intersections.
triangle Forcing points, trianglesolid Pseudo-fluid points, square fluid points, and squaresolid solid points
or
sn = ?bx(xi ?xn)?by(yj ?yn)2a
x(xi ?xn)+2ay(yj ?yn)
. (3.7)
With the substitution of Eq. (3.1) into Eq. (3.7) above, an equation with only sn
unknown can be obtained:
(2a2x +2a2y) s3n
+(3axbx +3ayby) s2n
+(2axcx +2aycy +b2x +b2y ?2axxi ?2ayyj) sn
+(bxcx +bycy ?bxxi ?byyj) = 0. (3.8)
55
The equation is solved iteratively using the Newton-Raphson method. The initial
solution to this equation is the closest interfacial marker particle, sb. Fig. 3.4 shows
the schematic of the solution procedure from sb to sn. After obtaining sn, the
interface-normal intersection coordinates (xn,yn) and the unit normal vector n at
the intersection can be calculated from Eq. (3.1) and Eq. (3.2).
3.2 Treatment of Stationary Immersed Boundaries
3.2.1 Momentum Forcing
To demonstrate the basic philosophy of the embedded-boundary formulation
let us examine the special case where an Eulerian grid point coincides with the
immersed interface, ?, on which a Dirichlet boundary condition, u?, has to be
imposed (see Case 2 in Fig. 3.5). The magnitude of external forcing function that
will enforce the above boundary condition can be obtained from Eq. (2.30) by simply
setting, ?uki = u?, and solving for fki :
fki = u? ?u
k?1
i
?t ?RHS
k
i, (3.9)
In the framework of the explicit fractional-step method discussed in the pre-
vious chapter, the use of fi, given by Eq. (3.9) will enforce the proper boundary
conditions on the predicted velocity, ?uki , (substitution of ( 3.9) into (2.30) gives
?uki = u?). This choice does not compromise the overall temporal accuracy of the
splitting scheme, and uki will also satisfy the same boundary condition to the order
of ?t2 [29].
56
A specific problem, which is related to the evaluation of fk for the embedded
boundary formulation, is raised when the implicit scheme is used, i.e., the momen-
tum forcing fk is no longer able to be obtained from the predicted field ?u, which
is obtained using the explicit scheme in the previous section. That is why the mo-
mentum forcing terms are omitted in the above equations. Kim et al. [29] proposed
a provisional explicit step for the evaluation of fk, which uses RK3 for Ai and the
forward Euler scheme for Bi in the above equations:
uk ?uk?1
?t = ?kA
parenleftbiguk?1parenrightbig+?
kA
parenleftbiguk?2parenrightbig+?
kB
parenleftbiguk?1parenrightbig
? ?k?pk +fk. (3.10)
This approach is adopted in the code and it works well for various cases. Notice that,
for the flow past a sphere reported in Chapter 5, the geometry is two-dimensional in
computational space without any variation in the azimuthal direction. Therefore,
the implicit treatment in the azimuthal direction will not be affected by the re-
quirement of fulfilling boundary conditions on embedded boundaries. This fact can
be utilized to simplify the time advancement procedure and the above provisional
explicit step can be omitted. On the other hand, the simplification is not applicable
any more when the diffusive terms in the radial direction are also treated implicitly.
The detailed implementation has been discussed in Sec. 2.2.3. As mentioned
before, for the flow past a sphere reported in this work, direct forcing of the predicted
velocity field without the preliminary Euler step still can be applied due to the two-
dimensionality of the geometry in cylindrical coordinates. Nevertheless, simulations
of the flow past a sphere have also been conducted with the preliminary step for
57
(a) (b)
Figure 3.5: Previous Interpolation schemes. (a) One-dimensional scheme in [15]; (b)
Two-dimensional scheme in [4]. Cases (1) and (3) illustrate two possible interpola-
tion stencils depending on the interface topology and local grid size.
the momentum forcing and the results have been found to be the same as the other
approach used.
3.2.2 Generalized Interpolation Scheme
The above idealized case, although it demonstrates the basic approach, rarely
appears in practical applications, where the Eulerian grid nodes almost never co-
incide with the immersed boundary. In such cases, fi, has to be computed at grid
points near ? and not exactly on ? the interface. For reasons that will became
apparent in the next section, we will compute fi at all fluid points that have at least
one neighbor in the solid phase (these points were labeled forcing points in the pre-
58
vious section). In this case, however, u? is not known and has to be reconstructed
using information from the interface and the surrounding velocity field.
In Fig. 3.5a, the interpolation stencil used by [15] is shown. They employed a
second order linear interpolation stencil along the grid lines. But this method has
ambiguity as shown in case (3) in the figure that there is no unique direction like
case (1) over which the interpolation can be performed. This ambiguity might not
decrease the accuracy of the method when the boundaries are stationary. However,
it might be a problem, as the flow has a determined direction near the immersed
interface when the interface is moving.
To eliminate this ambiguity, Balaras [4] proposed to perform the interpolation
along the well defined line normal to the boundary as illustrated in Fig. 3.5b: ini-
tially a virtual point is located along the normal; then, the virtual point together
with the point on the interface is used to perform the linear interpolation to find
u? in Eq. (3.9) at the location of the forcing point. The velocity at the virtual
point is computed from the surrounding fluid nodes using bi-linear interpolation.
In this last step, one also has to impose the constraint that the stencil does not
involve other forcing points, which can be easily achieved by gradually moving the
virtual point further away from the boundary (see for example Case 1 in Fig. 3.5b).
Implementation details and a series of applications in laminar and turbulent flows
are given in [4].
The interpolation stencil in this method, however, involves a search procedure
for the suitable fluid points and those points in the stencil may contain points that
do not belong to the neighboring grid lines of the forcing point considered as shown
59
Figure 3.6: Extra ghost cells required for the interpolation scheme in [4] when
domain decomposition technique used for parallelization.
in Fig. 3.6. When domain decomposition is employed in parallel computing, the
neighboring block should provide boundary conditions to the current block. In a
second order algorithm on Cartesian grids, only one layer of ghost cells are required
for the boundary conditions in the same way as a ?real? boundary treated in a
finite difference method. But the example case requires extra information from
the neighboring block. Of course, the information can be provided through extra
communication, provided that the structure of the variable arrays will be modified.
The reconstruction scheme that is used in the present study is based on the
work reported in [4], where u? at the forcing points is computed by means of linear
interpolation along the well defined line normal to the boundary. In this study a
60
Figure 3.7: Generalized Interpolation stencil.
variation of the above method is adopted that is better suited to moving boundary
problems. As will be demonstrated in the results section, the proposed method is as
accurate and robust as the one in [4]. It utilizes, however, a more compact stencil
and allows for the computation of all components of the strain rate tensor on the
interface in a straightforward manner. The later feature is particularly attractive in
fluid/structure interaction problems, while the former simplifies parallelization. The
similarities between the two approaches became apparent when comparing Fig. 3.5b
and Fig. 3.7. In the present approach we basically omit the virtual point and replace
it with two points on an x-grid line, y-grid line, or along the diagonal. Consequently,
the interpolation procedure is now a single-step process that involves two points
from the fluid and one on the interface (shaded area in Fig. 3.7). The interpolation
61
coefficients for the case of linear interpolation can be computed as follows: let us
assume that any variable ? in the two-dimensional space can be written in the
following form,
? = b1 +b2x+b3y. (3.11)
The coefficients b1, b2, and b3 in Eq. (3.11) can be found by solving the following
system: ?
??
??
??
?
b1
b2
b3
?
??
??
??
?
= A?1
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
?1
?2
?3
?
??
??
??
?
=
?
??
??
??
?
1 x1 y1
1 x2 y2
1 x3 y3
?
??
??
??
?
?1?
??
??
??
?
?1
?2
?3
?
??
??
??
?
, (3.12)
where (x1,y1), (x2,y2), and (x3,y3) in the 3 ? 3 matrix A are the coordinates of
the three points in the interpolation stencil shown in Fig. 3.7. For the case of a
stationary body, the inversion of matrix, A, at every forcing point is performed in
the beginning of the simulation and then stored in memory. For moving boundaries,
however, the inversion is performed at every time-step, since new forcing points
emerge as a consequence of the boundary motion. The extension of the above
interpolation procedure to three dimensions is straightforward (one needs to add
a term of the form, b4z, to the polynomial in Eq. (3.11) to reflect the dependence
of the solution near the interface on an additional spatial dimension. Higher-order
reconstructions can also be readily achieved using the same overall procedure (see
for example [35, 61]).
62
(a) (b)
Figure 3.8: Treatment of turbulent eddy viscosity for LES. (a) Test-grid cells for the
boundary points involving solid points from the interior; (b) Interpolation stencils
for ?t at the boundary points based the well-defined values from the fluid points away
from the immersed boundary and the boundary condition for ?t at the interface.
3.2.3 Treatment of Eddy Viscosity
For the simulation of laminar flows or direct numerical simulation of turbu-
lence, the embedded boundary formulation discussed above can satisfy the boundary
conditions on the immersed interfaces ?exactly? to the second-order accuracy of the
spatial discretization scheme for the Navier-Stokes equations. Also, the pressure
boundary conditions are implicitly satisfied as the results of linearization of the
velocity near the interface. However, for the SGS models the treatment of the
turbulent viscosity near the immersed boundary has to be considered.
The eddy viscosity ?t at the boundary points is required in the computation
63
of diffusion fluxes near the interface, and its calculation is no longer straightforward
due to the presence of solid body on the fixed grid. In the case of a dynamic SGS
model for example, the evaluation of all test-filtered quantities in the vicinity of the
immersed boundary, involves points from the interior of the solid body. A schematic
of the test-grid cells on a Cartesian grid with an immersed boundary is shown in
Fig. 3.8a. A direct computation ignoring the interface can be problematic and the
resulting eddy viscosity may be contaminated by the non-physical solution at the
solid points. Explicit modification of the filtering operator or redefining the test-grid
cells as in a cut-cell formulation can be considered, but a large number of cases will
be introduced and the algorithm becomes complicated.
Toavoidcomplexmodificationsofthefilteringoperatoratthesepoints, alinear
reconstruction procedure similar to that used for the velocity field is also applied
on ?t, which follows the same method proposed in [4]. The interpolation stencil for
?t is shown in Fig. 3.8. The linearized eddy viscosity is only an approximation of
the realistic value, which usually does not follow a linear distribution near the solid
boundary. In [4] extensive testing of this approach in LES over immersed boundaries
has been performed and the error in the results due to the linearized eddy viscosity
near the solid boundary has been found to be small.
64
3.3 Treatment of Moving Immersed Boundaries
3.3.1 Complications in Time Advancement
The general algorithm outlined in sections 3.1 and 3.2 above, is directly ap-
plicable to moving boundary problems provided that the grid-interface relation and
the interpolation coefficients are re-evaluated every time the location of the inter-
face is updated. In moving boundary problems, however, complications that are
usually related to the time advancement scheme can also arise. This is due to
the fact that the role of the Eulerian grid points near the interface changes from
timestep to timestep (for example a forcing point can became a fluid point at the
next timestep and vise-versa), as the solid body moves through the fixed grid. For
the time advancement scheme used in the present study, the evaluation of the RHS
of the momentum equation at substep k, requires physical values of the velocity
vector and pressure, as well as their derivatives from substep k?1 at all fluid points
(see Eq. 2.30). Due to the fact that the interface changes locations, it is possible
that some of the required values from substep k ?1 are not physical.
To identify the specific cases that are potential sources of error during the time
advancement procedure, let us consider a simple, one-dimensional problem where
the solid-phase: 1) encroaches upon the fluid phase, and 2) withdraws from the fluid
phase as shown in Figs. 3.9a and 3.9b respectively. Due to the CFL restriction of the
present scheme the interface cannot move by more than one gird cell in each substep,
resulting to two possible changes in the flags of the points near the interface. For
case 1: (a) forcing points =? solid points, and (b) fluid points =? forcing points.
65
(a) (b)
Figure 3.9: A simplified one-dimensional model problem that demonstrates the
possible changes of flags of the Eulerian grid nodes near a moving interface. (a)
Solid phase encroaches upon the fluid phase; (b) Solid phase withdraws from the
fluid phase.
Both of the above do not cause problems to the time advancement scheme. For case
1a, the solution at tk in all forcing points will be reconstructed anyway, and does not
depend on the history of the velocity or pressure field. Case 1b, where Eulerian grid
points move into the solid, also does not cause problems since the interior treatment
of the body does not influence the solution.
In case 2, where the solid phase withdraws from the fluid (see Fig. 3.9b), the
possible flag changes are the following: a) solid points =? forcing points; (b) forcing
points =? fluid points. The former case does not generate any problem for the
same reason that was discussed above (case 1a). In case 2b, however, it is apparent
from Eq. (2.30) that for all the newly defined fluid points the RHS which involves
66
Figure 3.10: The field-extension procedure. Shaded triangles are possible extrapo-
lation stencils for the pseudo-fluid points at substep tk?1 where the solution will be
extended.
derivatives from previous sub-steps will not be correct. This is because these points
at substep k ?1 were forcing points and although they have the correct values for
the velocity and pressure, their derivatives will not have physical values since they
involve points from the solid phase. These unphysical values can introduce spurious
vorticity near the boundary, leading to large errors.
3.3.2 Field Extension
To avoid forbiddingly complex special treatments of the computation of these
derivatives every time such a change of flag is detected during the time advancement
67
procedure, we propose a field-extension procedure, in which the velocity and pressure
fields are ?extended? in the solid phase at the end of each substep. Practically, with
this procedure the velocity fields are extrapolated at the Eulerian nodes defined as
pseudo-fluid points in Sec. 3.1 in the solid as illustrated in Fig. 3.10. The pressure
fields are also extended into the domain where the pressure solution is non-physical.
This way, not only the velocity and pressure at the forcing points at substep k?1,
but also their derivatives will have physical values, eliminating problems in the
computation of the RHS in the next substep for points in Case 2b. The values
of the velocity at the pseudo-fluid points is reconstructed using a procedure which
is similar to the one used for the forcing points that is described in section 3.2.
Central to this procedure is again the point on the interface where the normal from
the corresponding pseudo-fluid point intersects it. In addition to this point the
stencil involves two more points from the fluid phase as shown in Fig. 3.10. The
interpolation coefficients are computed by solving a system of equations identical
to that given by Eq. (3.12). The overall field-extension procedure is similar to the
generalized ghost-cell approach discussed in [26, 61]. In both these cases, however,
extrapolation is used to enforce boundary conditions on a stationary boundary.
We should also note that due to the staggering of the mesh the selection of
points on the pressure grid where the ?field-extension? is performed is not only
based on their relative location with the interface but also on their association to
velocity values that are non-physical. As a result the pressure is also ?extended?
to some nodes on the pressure grid that would have been classified as ?boundary?
nodes on the velocity grid as shown in Fig. 3.11
68
(a) (b)
Figure 3.11: Extension of pressure field. (a) Definition of pseudo-fluid points for
pressure; (b) Extrapolation stencils for pressure pseudo-fluid points. a3 Forcing
points for u velocity component, triangle Forcing points for v velocity component, trianglerightsld
Solid points for u velocity component, trianglesolid Solid points for v velocity component, square
fluid points for pressure, and squaresolid pseudo-fluid points for pressure.
3.3.3 Numerical Procedure
Having established the treatment of all Eulerian points in the vicinity of a
stationary or moving interface, we can summarize the overall algorithm as follows:
? Given the location of the interface at step k, identify fluid, forcing, solid, and
pseudo-fluid points on the Eulerian grid. This procedure needs to be performed
only once in the beginning of the computation for problems with stationary
boundaries.
69
? Calculate the predicted velocity field ?uki from Eq. (2.30).
? Reconstruct the predicted velocity ?uki at the forcing points as discussed in
Section 3.2.
? Solve the pressure Poisson equation (Eq. 2.31).
? Update the velocity field to uki (Eq. 2.32) and pressure field to pk (Eq. 2.33).
? Perform the field-extension procedure described above for the pseudo-fluid
points. This step is omitted for stationary boundary problems.
As it will be demonstrated in the results section the above algorithm is very
robust and satisfies the boundary conditions on an immersed interface ?exactly?,
within the overall accuracy of the discretization scheme. Also, contrary to alterna-
tive reconstruction procedures (i.e. generalized ghost cell approaches) it does not
require special treatment of the grid cells emerging from the solid, which is usually
ad-hoc (see [63] for a more detailed discussion on this issue). We should also note
that in the above algorithm, boundary conditions for the pressure near the inter-
face are not imposed directly, but they are essentially implicit into the RHS of the
Poisson equation (2.31). To clarify this statement one can consider the projection
of the momentum equation of a fluid particle on the solid boundary in the direction
normal to the interface:
Du
Dt ?n = ??p?n+
1
Re?
2u?n. (3.13)
With the present reconstruction procedure the viscous term turns out to be zero as
all velocity components are linearized in the vicinity of ? (see Sections 3.2 and 3.3),
70
and the above reduces to:
?p
?n = ?
Du
Dt ?n, (3.14)
Eq. (3.14) represents the boundary condition for the pressure which is usually en-
forced ?directly? in boundary-conforming or cut-cell formulations involving moving
boundaries (see for example [63]). In these cases, however, it is obtained from
the projection of the inviscid momentum equation in the direction normal to the
boundary. A discussion on the behavior of the pressure in embedded-boundary for-
mulations is also given in [15]. In the results section, we will show that the above
argument is sound and the pressure behaves correctly in the vicinity of the immersed
boundaries.
3.4 Surface Force Calculation
In the case of non-boundary-conforming formulations the computation of the
local forces on the surface of the body is a non-trivial problem. A good discussion
of this issue can be found in the study by Lai and Peskin [30], where the accuracy
of several different force calculation algorithms is discussed in the framework of
their immersed-boundary formulation. In the present work we propose a method to
compute the local force distribution on a complex body that is based on the same
reconstruction stencil used for the velocity field. In general, the force f per unit area
on a surface element with an outward normal n can be written as:
fi = ?jinj =
bracketleftbigg
?p?ij +?
parenleftbigg?u
i
?xj +
?uj
?xi
parenrightbiggbracketrightbigg
nj, (3.15)
where fi is the surface force component in the xi direction, ?ji is stress tensor,
71
Figure 3.12: Surface force calculation. (a) Arrangement of flow variables on the
staggered grid: triangleright u velocity, triangle v velocity, square pressure, ? marker particles. The
corresponding filled symbols on the interface represent the location where normal
intersects the boundary. (b) detail in the vicinity of the boundary for the grid of
u velocity component. (c) detail in the vicinity of the boundary for the the grid
of pressure. Interpolation stencil is the shaded triangle and line with arrow is the
normal to the interface.
72
? is the dynamic viscosity of the fluid, and nj is the direction cosine of n in xj
direction. To compute f from Eq. (3.15) one would need to find p and ?ui/?xj on
the body surface. In the framework of the present scheme, the later term can be
computed in a straightforward manner using the stencil and interpolation coefficients
that were utilized to reconstruct the velocity field near the interface (Eqs. 3.11 and
3.12). Due to grid staggering, however, the derivatives for each velocity component
are computed on different points on the interface. Fig. 3.12a, for example, shows
the staggered grid arrangement of all flow variables in the vicinity of a solid body
and the corresponding points on the interface where ?ui/?xj is computed for each
velocity component. In a detailed view around a u velocity point (see Fig. 3.12b),
one can identify the interpolation stencil that was used in the reconstruction step. It
involves points (2) and (3) from the flow, and point (1) on the interface. Considering
Eq. (3.11) at point (1) and differentiating with respect to xj, we can compute the
desired velocity derivatives as follows:
?u
?x = b2 and
?u
?y = b3, (3.16)
where the coefficients b2 and b3 have been already computed form Eq. (3.12). In a
similar manner, the derivatives for v and w can be found at the interface-normal
intersection points of each component using the corresponding interpolation stencils.
The computation of the surface pressure, p, is a little more complicated if
one is to avoid simple extrapolation from the interior, which can result in large
errors. In the present study we have tested several different strategies and found
that an accurate and cost/efficient approach to compute the surface pressure is
73
the one that utilizes the momentum equation normal to the boundary ( Eq. 3.14),
in addition to information from the interior of the flow. In particular, Eq. (3.14)
can be used to obtain an estimate of ?p/?n on any point on the interface, since
the material derivative in the RHS of can be computed directly from the known
boundary velocities. The normal derivative of pressure can then be expressed as:
?p
?n =
?p
?xnx +
?p
?yny, (3.17)
where nx and ny are the components of normal unit vector, n, in the x and y
directions, respectively. Referring to the stencil shown in Fig. 3.12c, differentiation
of Eq. (3.11) at point (1) yields:
?p
?x = b2 and
?p
?y = b3. (3.18)
This stencil is built as if the pressure near the interface will be reconstructed in
a manner similar to the velocity field, even though it is only used for the surface
pressure computation.
Using equations (3.17) and (3.18) we can rewrite Eq. (3.12) as follows:
?
??
??
??
?
b1
b2
b3
?
??
??
??
?
= A?1
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglevextendsingle
?p
?n
p2
p3
?
??
??
??
?
=
?
??
??
??
?
0 nx ny
1 x2 y2
1 x3 y3
?
??
??
??
?
?1?
??
??
??
?
?p
?n
p2
p3
?
??
??
??
?
. (3.19)
We can now use Eq. (3.19) to compute b1, b2 and b3, and then compute the surface
pressure in a straightforward manner from Eq. (3.11). As also mentioned above,
due to the grid staggering the velocity derivatives and surface pressure have to be
computed on different points on the interface. We then use cubic spline interpolation
74
to compute the corresponding values at the locations of the marker particles. Finally
the local forces can be computed from Eq. (3.15) at the marker particle location.
Numerical integration is used to calculate integral quantities like drag and lift forces.
3.5 Fluid-Structure Interaction
Up to this point the imposition of boundary conditions on a body moving
through a fixed orthogonal grid has been discussed in detail. The boundary con-
ditions were assumed to be known functions of space and time. This information
would be sufficient to solve problems with prescribed motion. However, in cases
where the location of the body is determined by the fluid force, additional ODEs
governing the motion of the body have to be considered. Also the coupling between
the Navier-Stokes equations and ODEs becomes an important aspect.
In the present study the fluid and the structure will be treated as elements of
a single dynamical system, and all governing equations will be integrated simulta-
neously, and iteratively in the time-domain. A fundamental complication with the
application of a time-domain approach to fluid-structure interaction problems like
the one described above is that the prediction of the hydrodynamic loads requires
knowledge of the motion of the structure and vice-versa. To overcome this com-
plication an iterative predictor-corrector scheme that accounts for the interaction
between the hydrodynamic loads and the motion of the structure has been devel-
oped. The procedure is based on Hamming?s 4th-order, predictor-corrector method
[9]. Hamming?s scheme was adapted to the present problem to avoid evaluating the
75
(a) (b)
Figure 3.13: Schematic of rigid-body motion. (a) One degree of freedom of transla-
tion in vertical direction; (b) One degree of freedom of rotation about the origin.
hydrodynamic loads at fractional time steps and to accommodate hydrodynamic
loads that are proportional to the acceleration (the so called added-mass effect).
In the following paragraphs the structural solver and the coupling scheme will be
discussed in detail. A schematic of rigid-body motion is shown in Fig. 3.13. The
structure is modeled as a mass-spring-damper system. Fig. 3.13a shows one degree
of freedom of translation in the vertical direction, and Fig. 3.13b shows the one
degree of rotation about the origin. The motions are governed by the same ODE.
Although in Chap. 5 the rigid-body rotation of the leaflets will be considered, the
motion is prescribed instead of governed by a separated equation. Therefore, in
this work only the problem of vortex-induced vibrations for a circular cylinder is
76
considered. Here we introduce the equation of motion that generally governs the
motion of the cylinder oscillating in the (X ?Y) plane:
M?x(t)+C?x(t)+Kx(t) = f [x(?), ?x(?),?x(t)], ?? ? 0 ? ? ? t (3.20)
where M is the mass matrix, C is the damping matrix, K is the stiffness matrix
and f is the vector of generalized hydrodynamic forces per unit span, and x(t) =
X0 (t)i+Y0 (t)j+?(t)k is the vector of generalized structural displacements. Here,
X0 (t) and Y0 (t) are the displacements of the center of mass of the cylinder in the
x and y directions, respectively, and ?(t) is the rotation of the cylinder around an
axis perpendicular to the (X ?Y) plane. Equations (2.1) and (2.2) governing the
dynamics of the fluid and equation (3.20) governing the dynamics of the body need
to be solved in a coupled manner. The overall algorithm will be given later.
Equation (3.20) can be rewritten in non-dimensional form, as a system of 2n
first-order, non-linear, ordinary-differential equations (n is the number of degrees-
of-freedom of the structure) as follows:
?y(t) = F[y(?), ?y(?)], ?? ? 0 ? ? ? t, (3.21)
where half of the vector F represents generalized velocities and the other half repre-
sents generalized forces divided by the corresponding inertias. In general, the loads
depend explicitly on y and implicitly on the history of the motion and the accel-
eration of the structure. For this reason Hamming?s 4th-order, predictor-corrector
method [9] is used to integrate equation (3.21) in the time domain. The details of
the basic numerical procedure used to determine the current value of the vector y
are given next:
77
1. Let tj = j?t denote the time at the j-th time step, where ?t is the time-step
size used to obtain the numerical solution, and
yj = y(tj), ?yj = ?y(tj), Fj = F[y(tj)] (3.22)
2. Compute the predicted solution, pyj, and modify it, 1yj, using the local trun-
cation error, ej?1, from the previous time step:
pyj = yj?4 + 4
3?t
parenleftbig2Fj?1 ?Fj?2 +2Fj?3parenrightbig, 1yj =p yj + 112
9 e
j?1 (3.23)
3. Correct the modified, predicted solution by:
k+1yj = 1
8
bracketleftbig9 yj?1 ?yj?3 +3?t parenleftbig kFj +2 Fj?1 ?Fj?2parenrightbigbracketrightbig (3.24)
where kFj = Fparenleftbigkyjparenrightbig and 1yj = pyj. k is the iteration index. Convergence
is achieved when the iteration error, ej = vextenddoublevextenddouble k+1yj ? kyjvextenddoublevextenddouble? is less than a
prescribed tolerance epsilon1.
4. Compute the local truncation error, ej, and the final solution, yj, at step j
and advance to the next timestep:
ej = 9121parenleftbig k+1yj ? pyjparenrightbig, yj = k+1yj ?ej (3.25)
Note that in the above procedure information from four previous timesteps is re-
quired. Therefore, starting from the initial conditions, Euler, Adams-Moulton, and
Adams-Bashforth methods were used.
The fluid solver described in previous chapter is integrated into the above
procedure as follows:
78
1. Find the predicted location and velocity of the body using equation (3.23).
2. Find the predicted fluid velocity and pressure fields using equations (2.30)?
(2.33) with the boundary conditions provided by step 1. Then, compute the
resulting loads on the structure.
3. Compute the new location and velocity of the body using equation (3.24).
4. Check for convergence. If ej is greater than the prescribed tolerance, epsilon1, repeat
steps 2 to 4. If convergence is achieved go to step 5.
5. Find final position and velocity of the body using equation (3.25), and the
final fluid pressure and velocity fields from equations (2.30)?(2.33).
In all computations reported in this study a tolerance of epsilon1 = 10?8 was used.
The number of iterations required for convergence at each time step varied from 2
to 6 depending from the stiffness of the problem.
3.6 Parallelization
The parallelization is implemented via slab decomposition and the parallel
communication between among blocks is done using the MPI library. To minimize
the communications at the block (slab) interface, the decomposition is carried out
in the direction with the largest number of points (usually the streamwise direc-
tion). Fig. 3.14a shows the slab decomposition in the streamwise direction for the
momentum equations. Here, each block is equally sized, which greatly simplifies
the mesh partition and provides optimal load balance. Minor modifications are
79
(a) (b)
Figure 3.14: Domain decomposition parallelization. (a) Slab decomposition in
streamwise direction for momentum equations; (b) Slab decomposition in wave-
number space for the Poisson equation.
introduced to the original serial solver, as ghost cells at the block interfaces are
utilized to provide information from neighboring blocks. A schematic of message
communication among slabs is shown in Fig. 3.15. In the current implementation,
only one layer of ghost cells from the neighboring block are required as the stan-
dard second-order central difference scheme is used for the spatial discretization.
Usually the solution of the pressure Poisson equation is the most expensive part
of a Navier-Stokes solver, and the parallelization via domain decomposition adds
further complications. In the present study, the computational grid is always uni-
form in the spanwise direction; thus the fast Fourier transformation (FFT) can be
applied to the Poisson equation in the spanwise direction and a series of decou-
pled two-dimensional problems are obtained and solved using the direct solver from
80
Figure 3.15: Inter-slab communication via layers of ghost cells.
FISHPACK [57]. The FFT is performed in each block and no information from
neighboring blocks is required. However, global communication is required to swap
the slabs in th estreamwise direction into slabs in wavenumber space, and swap
back after the direct solution of the two-dimensional problem for each wavenumber
(Fig. 3.14b shows a slab decomposition in wave number space). Nevertheless, the
overhead for the global communication is small and, as shown in the results which
follow, the speedup for the problem considered is linear.
The issue of the parallelization of the embedded boundary formulation has
been addressed in Sec. 3.2.2 and will be discussed here.
The parallel speedup for a simulation including the embedded boundary treat-
ment is shown in Fig. 3.16. Linear speedup is obtained up to 16 processors, which
81
Figure 3.16: Parallel Speedup for a 96?64?384 Grid.
was the limit of the cluster where the preliminary testes were carried out.
82
Chapter 4
Two-Dimensional Validation Cases
To evaluate the accuracy of the proposed methodology several cases with in-
creasing complexity have been considered. The first case in Sec. 4.1, which simulates
the two dimensional flow around one element of a moving cylinder-array in a planar
channel, is a demonstration of the formal accuracy of the method. Then, laminar
flow problems involving prescribed motions of two-dimensional bodies are computed
in Sec. 4.2: the flow induced by the harmonic in-line oscillation of a circular cylinder
in a quiescent fluid (Sec. 4.2.1), and the flow from a transversely oscillating cylinder
in a free-stream (Sec. 4.2.2). Both cases are characterized by complex vortex-vortex
and vortex-structure interactions and present a challenge to non-boundary conform-
ing formulations. In addition, detailed experimental data and numerical results from
boundary-fitted methods are available in the literature, allowing for a detailed val-
idation. Then, in Sec. 4.3 vortex-induced vibrations of circular cylinders will be
used in the present study to evaluate the accuracy and efficiency of the proposed
fluid-structure interaction scheme. Two different configurations are considered: one
degree-of-freedom oscillations in the cross-stream direction (Sec. 4.3.1); and two
degree-of-freedom oscillations in both the streamwise and cross-stream directions
(Sec. 4.3.2). In both cases the selected parametric space was as close as possible to
well documented laminar two-dimensional flow problems in the literature.
83
4.1 Formal Accuracy
(a) Stationary cylinder (b) Moving cylinder
Figure 4.1: Flow around a periodic array of cylinders in a plane channel: Sketch of
the computational domain.
There are no cases in fluid dynamics that involve complex moving boundaries
and have analytical solutions, to serve as reference for the formal accuracy study
of the proposed method. For this reason we have designed a test problem which
involves a body moving through a fixed grid, and performed a detailed numerical
accuracy study. In particular, we have considered the flow around a periodic array
of cylinders moving with constant velocity W = ?1 in a planar channel. A sketch
of the computational domain and boundary condition setup is shown Fig. 4.1a.
Non-slip conditions are used on the walls and periodic boundary conditions in the
direction of motion. When the reference frame moves with the cylinder, the case
84
(a) Stationary cylinder (b) Moving cylinder
Figure 4.2: Streamlines on a 150?150 uniform grid colored by u velocity component
(a reference frame moving with the body is used for the moving cylinder case).
shown in Fig. 4.1b is recovered, where the cylinder is stationary and the walls have a
velocity of W = 1. The Reynolds number based on the diameter and velocity of the
cylinder (or the velocity of the moving wall) is 100. The corresponding flow fields
are visualized in Fig. 4.2 where streamlines are shown. A recirculation bubble is
formed and restricted between two cylinders, which keeps the flow two-dimensional
and steady. Note that for the moving cylinder a reference frame following the
cylinder is used. In both cases, the flow is identical even though we have two
distinct computational configurations.
Fig. 4.3 shows the instantaneous pressure contours when the cylinder is at
two different locations: the center and the boundary of the computational domain.
The almost identical pressure distributions in Fig. 4.3a and Fig. 4.3b also verify
85
(a) t/T = 0.5 (b) t/T = 1.0
Figure 4.3: Instantaneous pressure contours at two different locations.
the stableness and accuracy of the current method.
Calculations on four levels of uniform grids were carried out: 30?30, 90?90,
150 ? 150, and 450 ? 450. The above choice -there is not the usual doubling of
grid points in the next finer grid level- is due to the staggered variable arrange-
ment. We have selected the specific grids in a way that every node on coarser grids
can be directly located on the reference, finest, grid to avoid interpolations when
comparing the two solutions. The L1 and L2 norms of the error, which measure the
difference between the solutions from coarser grids and the reference grid, versus the
corresponding grid resolution are shown in Fig. 4.4. The results from the stationary
cylinder case and the moving cylinder case are given in Fig. 4.4a and Fig. 4.4b,
respectively. Note that we use the solution on the finest grid from the stationary
cylinder case for both plots presented. For the moving cylinder case, the L1 and L2
86
(a) Stationary cylinder (b) Moving cylinder
Figure 4.4: a50 L1 norm of the error; ? L2 norm of the error. Open symbols are for
w and filled for u.
norms of the error computed using the solution on the 450?450 grid are almost the
same as those in Fig. 4.4b. In the figure, lines representing the first- and second-
order accuracy slopes are included. It is clearly shown that the error decreases with
?x2, verifying the second-order accuracy of the present method. This means that,
within the order of the accuracy of the discretization scheme, the moving boundaries
are represented ?exactly?. By comparing Fig. 4.4b with Fig. 4.4a, we can see that
the error of the moving cylinder case is of the same order as the stationary cylinder
case.
87
4.2 Forced Vibrations of a Cylinder
4.2.1 In-Line Oscillating Cylinder in a Fluid at Rest
The interaction of an oscillating circular cylinder with a fluid at rest is a well
documented model-problem in the literature because of its relevance to a variety
of applications. The two key parameters in this flow are the Reynolds number,
Re = ?UmaxD/?, and the Keulegan-Carpenter number, KC = Umax/fD (Umax is
the maximum velocity of the cylinder, D is the diameter of the cylinder, ? is the
density of the fluid, ? is the dynamic viscosity of the fluid, and f is the character-
istic frequency of the oscillation). The parametric space we selected in the present
work corresponds to the LDA experiments and numerical simulations reported by
D?utsch et al. [13]. In particular, the cylinder?s translational motion is described by
a simple harmonic oscillation, x(t) = ?Asin(2pift), where A is the amplitude of the
oscillation. The Reynolds and Keulegan-Carpenter numbers were set to Re = 100
and KC = 5 respectively. For this setup the flow remains two-dimensional with
periodic vortex shedding. The size of the computational domain is 50D ? 30D in
x and z, respectively, with the cylinder located at the center. Radiative boundary
conditions are applied to all the far-field boundaries. Three different grid levels were
used (240?120, 480?240, and 960?480) to investigate the sensitivity of the results
to numerical resolution. Near the cylinder the grid is approximately uniform in both
directions with local spacing of 0.05D, 0.025D, and 0.0125D respectively.
All computations were started from a quiescent field and the integration in
time was performed until periodic vortex shedding was established. In Fig. 4.5
88
Pressure contours Vorticity contours
(a) 0?
(b) 96?
(c) 192?
(d) 288?
Figure 4.5: In-line oscillating cylinder in a fluid at rest (Re = 100 and KC = 5).
Pressure and vorticity contours at four different phase-angles. ?1.1 < P < 0.6 with
intervals of 0.09 and ?26 < ?y < 26 with intervals of 0.95.
89
Measured Computed
(a) 180?
(b) 210?
(c) 330?
Figure 4.6: Measured (left) and computed (right) velocity vectors and streamlines
in the vicinity of the cylinder at Re = 100 and KC = 5.
90
x = ?0.6D x = 0.0D x = 0.6D x = 1.2D
Figure 4.7: In-line oscillating cylinder in a fluid at rest (Re = 100 and KC =
5). Profiles of u velocity component at: (a) 180?; (b) 210?; (c) 330?. Symbols:
experiment; Lines: Computation.
91
x = ?0.6D x = 0.0D x = 0.6D x = 1.2D
Figure 4.8: In-line oscillating cylinder in a fluid at rest (Re = 100 and KC =
5). Profiles of w velocity component at: (a) 180?; (b) 210?; (c) 330?. Symbols:
experiment; Lines: Computation.
92
x = ?0.6D x = 0.0D x = 0.6D x = 1.2D
Figure 4.9: Convergence of the u velocity component for different grid levels at
Re = 100 and KC = 5. (a) 180?; (b) 210?; (c) 330?. : grid 240?120; :
grid 480?240 ; : grid 960?480.
93
x = ?0.6D x = 0.0D x = 0.6D x = 1.2D
Figure 4.10: Convergence of the w velocity component for different grid levels at
Re = 100 and KC = 5. (a) 180?; (b) 210?; (c) 330?. : grid 240?120; :
grid 480?240 ; : grid 960?480.
94
the pressure and vorticity contours at four different phase-angles are shown. One
can observe the two fixed stagnation points at both ends of the cylinder that also
define the axis of symmetry of the flow. As the cylinder moves to the left, two thin
boundary layers develop on the upper and lower sides (see Fig. 4.5a), that eventually
separate giving rise to two counter-rotating vortexes of exactly the same strength.
The vortex generation procedure stops as the body reaches its extreme left location,
as shown in Fig. 4.5b. Then the cylinder moves backwards, and the same process
takes place on its right side. During this part of the cycle, the vortex pair generated
earlier is split and pushed away the cylinder (see Fig. 4.5c). These results are in
very good agreement with the corresponding results reported in [13], demonstrating
that the present method can properly capture the dynamics of the vorticity field.
Fig. 4.6 shows the comparisons of measured and computed velocity vector fields at
three different phase angles of the cylinder motion. Very good agreement between
the experimental and computed results is emphasized by the coincidence of the
measured and computed streamlines in the vicinity of the cylinder.
To further examine the accuracy of the proposed method, quantitative com-
parisons with the experimental results are considered. Fig. 4.7 and Fig. 4.8 show the
computed velocity profiles on the intermediate grid, at four x locations and three
different phase-angles, in comparison with the experiments. The agreement is very
good. We also investigate the grid independence of our solution by comparing the
velocity profiles obtained from the three grids in Fig. 4.9 and Fig. 4.10. It is very
clear that even the coarse grid 120 ? 240 gives very good results. As the grid is
refined, there is a obvious convergence of the solution. In Fig. 4.11 the evolution in
95
Figure 4.11: Evolution in time of the in-line force acting on a cylinder oscillating
in a fluid at rest at Re = 100 and KC = 5. ? Boundary-conforming simulation
[13]; present computation 240?120; present computation 480?240;
present computation 960?480.
time of the in-line force, Fx(t), acting on the cylinder is shown for the two finest grids
in comparison with the reference simulation results. Details on the actual definition
and normalization of Fx(t) can be found in [13]. Here it is important to point out
that Fx(t) is practically calculated by integrating the surface shear stress and pres-
sure using the algorithm described in Section 3.4. On both fine grids the agreement
with the reference computations is very good (results are within 2%) demonstrating
that the forces on the cylinder?s surface (or at least their integral at an instant in
time) can also be predicted very accurately with the proposed approach.
4.2.2 Transversely Oscillating Cylinder in a Free-Stream
The fluid-structure interaction problem becomes more complicated when the
oscillating cylinder is in a free-stream. In the present study we only consider cases of
96
(a) fe/f0 = 0.80 (b) fe/f0 = 0.90
(c) fe/f0 = 1.00 (d) fe/f0 = 1.10
(e) fe/f0 = 1.12 (f) fe/f0 = 1.20
Figure 4.12: Cylinder oscillating in a free-stream: instantaneous streamline patterns
when cylinder is located at its extreme upper position. Streamlines are colored with
transverse velocity component values.
97
transverse oscillation. The relevant parameters for such a problem are the Reynolds
number Re = ?U?D/? (U? is the free-stream velocity), the excitation frequency,
fe, and the oscillation amplitude, A. When A is larger than a threshold value the
synchronization condition can appear -manifested by a jump in the phase-angle, ?,
between the lift force and the body motion- when the frequency of the oscillation
is varied around the natural shedding frequency, f0, of the stationary cylinder. The
variations in this phase angle are associated with the change in sign of the energy
transfer between the fluid and the cylinder. Capturing the physics of these flows,
even at low Reynolds numbers, requires a very accurate prediction of the vorticity
production and dynamics on the surface of the body, which makes their simulation
a very challenging task, especially using non-boundary-conforming formulations.
The choice of parameters in this flow reproduces the conditions in the ex-
periments by Gu et al. [22] and the subsequent simulations by Guilmineau &
Queutey [23], where an accurate boundary-conforming numerical method is used.
To conduct quantitative comparisons, however, we mainly use the numerical results
since the amount of quantitative information reported in the experiments is lim-
ited. To match the above conditions the prescribed cylinder motion is chosen to
be z(t) = Acos(2pifet), where z(t) is the cross-stream location of the cylinder as a
function of time. Six different cases were considered at Re = 185, A = 0.2D and
0.8 ? fe/f0 ? 1.2. This way we cover a wide spectrum of the conditions in [22, 23].
The computational domain is 50D?30D, in the streamwise and cross-stream
directions respectively. The cylinder is located at a distance of 10D from the the
inlet boundary, and 15D from the upper and lower boundaries. At the outflow
98
(a) fe/f0 = 0.80 (b) fe/f0 = 0.90
(c) fe/f0 = 1.00 (d) fe/f0 = 1.10
(e) fe/f0 = 1.12 (f) fe/f0 = 1.20
Figure 4.13: Cylinder oscillating in a free-stream: instantaneous vorticity contours
when cylinder is located at its extreme upper position. ?8 < ?y < 8 with intervals
of 0.5.
99
(a) fe/f0 = 0.80 (b) fe/f0 = 0.90
(c) fe/f0 = 1.00 (d) fe/f0 = 1.10
(e) fe/f0 = 1.12 (f) fe/f0 = 1.20
Figure 4.14: Cylinder oscillating in a free-stream: instantaneous pressure contours
when cylinder is located at its extreme upper position.
100
boundary a convective boundary condition is used, and radiative conditions are
imposed at the freestream boundaries. The grid is stretched in both directions to
cluster points in the vicinity of the body. To investigate the dependency of our
results on grid resolution three different meshes with an increasing number of nodes
are considered: 200?160, 400?320, and 800?640 points in the streamwise and
cross-stream directions respectively. For all three grids the flow past a fixed cylinder
at Re = 185 is initially computed, to obtain the natural shedding frequency, f0. As
for the previous test problem the results converge with grid refinement. For example,
CD = 1.366, CDrms = 0.029, and CLrms = 0.461 on the finest grid, which are within
1% of the values on the previous grid level. The Strouhal number is 0.197 for all
the three grids.
Fig. 4.12 shows the instantaneous streamline patterns, Fig. 4.13 shows the
vorticity contours on the finest grid, and Fig. 4.14 shown the instantaneous pressure
contours when the cylinder is at the extreme upper position. Six cases are shown,
where fe/f0 is gradually increased from 0.8 to 1.2. A sudden change in the topology
of the near wake can be observed at fe/f0 = 1.10, (see Fig. 4.12d, 4.13d, and 4.14d),
as manifested by the appearance of concentric, closed streamlines, and remains the
same up to the highest value of excitation frequency. This topology suggests the
existence of concentrated vorticity in the near wake, which can be verified from the
vorticity contours shown in the same figure. The different wake structure is also
reflected in the lift and drag coefficients shown in Fig. 4.15. For values of fe/f0
greater than 1.0, both the drag and lift exhibit signs of the influence of a higher
harmonic. This behavior can be attributed to the strengthening of the opposite-
101
Figure 4.15: Evolution of the drag and lift coefficients in time for the case of
a cylinder oscillating in a free-stream: CD; CL. (a) fe/f0 = 0.8; (b)
fe/f0 = 0.9; (c) fe/f0 = 1.0; (d) fe/f0 = 1.1; (e) fe/f0 = 1.12; (f) fe/f0 = 1.2.
102
sign vorticity formed at the base of the cylinder as fe/f0 increases, which alters the
vorticity in the top and bottom shear layers. The above results are in excellent
agreement with reference data in the literature [22, 23].
Quantitative comparisons of the predicted drag and lift coefficients for all
frequencies with the corresponding values from the simulations by Guilmineau &
Queutey [23] are shown in Fig. 4.16a. The variation of the average drag coefficient,
CD, with fe/f0 is in agreement with the reference data, although it is a little higher
in our simulations (approximately 5%) for all frequency ratios. In the reference
simulations, however, CD was found to be very sensitive to grid resolution and they
reported higher values when the grid was refined for the case of fe/f0 = 1.10. The
r.m.s values of the drag and lift coefficients are less sensitive to grid resolution and
are in excellent agreement with the reference data. The same applies to the phase
angle, ?, betweenthe liftcoefficient and theverticaldisplacement of thecylinderthat
is shown in Fig. 4.16b. Both ? and CLrms exhibit a sudden change at fe/f0 = 1.10
when the vortex switching occurs as shown in [22, 23].
Next, we will examine how accurately the local forces on the cylinder?s surface
are predicted. The fact that CD, CDrms, and CLrms above are in good agreement
with the reference simulations, does not guaranty the same agreement on the local
skin friction distribution, for example, which contributes a small percentage in the
above quantities. Guilmineau & Queutey [23] reported the local distribution of
surface pressure and vorticity, when the cylinder is located at the extreme upper
position. A comparison with our results on all grids for a range of frequencies is
shown in Fig. 4.17 and 4.18. The angle ? on the x axis, is measured clockwise from
103
Figure 4.16: (a) Comparison of force coefficients for the case of the cylinder oscil-
lating in a free-stream: & ? CD; & triangle CDrms; & a50 CLrms; lines are
present results on the 800?640 grid and symbols are the corresponding data from
[23]. (b) Comparison of the phase angle between the lift coefficient and the vertical
displacement: present results; a50 [23]
104
Figure 4.17: Distribution of the pressure coefficients on the cylinder?s surface, when
located at the extreme upper position. grid 800 ? 640; grid 400 ? 320;
grid 200?160; ?body-fittedreference computationin [23]. (a)fe/f0 = 0.80; (b)
fe/f0 = 0.90; (c) fe/f0 = 1.00; (d) fe/f0 = 1.10; (e) fe/f0 = 1.12; (f) fe/f0 = 1.20.
105
Figure 4.18: Distribution of the skin friction coefficients on the cylinder?s surface,
when located at the extreme upper position. grid 800?640; grid 400?320;
grid 200?160; ?body-fittedreference computationin [23]. (a)fe/f0 = 0.80; (b)
fe/f0 = 0.90; (c) fe/f0 = 1.00; (d) fe/f0 = 1.10; (e) fe/f0 = 1.12; (f) fe/f0 = 1.20.
106
the stagnation point. It can be seen that our results follow the reference data very
well especially on the finest grid. As expected the skin friction is more sensitive
to grid resolution and is under-predicted near the peaks on the coarser grids. The
pressure coefficient on the other hand is fairly insensitive and almost identical on
all grids. We should also note that the body-fitted grid that is used in [23] is much
finer near the cylinder?s surface compared to our finest grid. The local grid spacing
normal to the boundary is 0.001D in [23] and 0.005D for our 800?640 grid, which
further demonstrates the accuracy of the overall approach.
4.3 Vortex-Induced Vibrations of an Elastic Cylinder
Vortex-induced vibration of a circular cylinder is a problem that has been ex-
tensively studied both experimentally and numerically (see [67] for a recent review),
and will be used in the present study to evaluate the accuracy and efficiency of the
proposed methodology. Two different configurations are considered: 1. Free oscil-
lations in the cross-stream direction (one degree-of-freedom); 2. Free oscillations in
both the streamwise and cross-stream directions (two degree-of-freedom). In both
cases the selected parametric space was as close as possible to well documented
laminar two-dimensional flow problems in the literature.
4.3.1 One Degree of Freedom
For a cylinder freely oscillating in the cross-stream direction, which is mod-
eled as mass-damper-spring system, the non-dimensional form of equation (3.20)
107
becomes:
?y +2?
parenleftbigg 2pi
Ured
parenrightbigg
?y +
parenleftbigg 2pi
Ured
parenrightbigg2
y = 12nCL, (4.1)
where y = Yo/D is the dimensionless vertical displacement, ? = c2?km is the damping
ratio, CL = fy/(12?DLU2) is the lift force coefficient, n = m/?D2 is the mass
ratio with ? the fluid density, and Ured = U/fnD is the reduced velocity with
fn = 12piradicalbigk/m the natural vibrating frequency of the structure. In the above m
is system mass, c the damping coefficient, k the spring constant, Yo the transverse
displacement of the system centroid, and fy the instantaneous lift force on cylinder.
Note that the same reference scales as in the Navier-Stokes equations (the cylinder
diameter, D, and the freesream velocity, U), were used to obtain the dimensionless
form of the equations.
The parametric space was selected to match an experiment by Anagnostopou-
los and Bearman [2], which has been used for validation in several other numerical
simulations (see for example [53, 31]). The mass ratio was set to n = 117.10 and
the damping ratio to ? = 0.0012. The computational domain is 40D ? 20D, in
the streamwise and cross-stream directions, respectively. The cylinder is located at
a distance of 10D from the the inlet boundary, where a uniform velocity is speci-
fied, and 10D from the upper and lower boundaries, where radiative conditions are
imposed. At the outflow boundary a convective boundary condition is used [41].
The number of grid points in all computations was 320 ? 240 in the streamwise
and transverse directions, respectively. The grid was stretched in both directions
and the resulting resolution near the cylinder was approximately 0.02D?0.02D. A
108
Displacements Force coefficients
(a) Re = 90
(b) Re = 100
(c) Re = 110
(d) Re = 120
Figure 4.19: Time history of the displacements for the freely vibrating cylinder in
the cross-stream direction: start-up phase.
109
constant time step of 0.005D/U was used for all cases.
For all cases the flow over a stationary cylinder for the same Reynolds number
was initially computed. Then, the cylinder was allowed to move freely in the trans-
verse direction. Integration in time was performed until a periodic state of constant
maximum amplitude was reached. The large mass ratio of this system makes the
problem fairly stiff and very long integration times were necessary for most cases.
It is interesting to investigate how the cylinder behaves after the release. Fig. 4.19
shows the start-up phase and Fig. 4.20 shows the stable phase of the cylinder dis-
placement and corresponding drag and lift force coefficients. Outside the lock-in
regime the vibration amplitude is modulated and has a very small magnitude, e.g.,
as shown in Fig. 4.19a for Re = 90. Correspondingly, the drag coefficient is oscil-
lating in a very narrow range. Also, the time histories of cylinder displacement and
the force coefficients do not change much even after reaching the stable phase as
shown in Fig. 4.20a. In the lock-in regime, for most cases the amplitudes are also
modulated (see Fig. 4.19c,d). However, for Re = 100, one can observe a monotoni-
cally growing oscillation amplitude, as shown in Fig. 4.19b. Nevertheless, all cases
in the lock-in regime approach the stable periodic states. Similar phenomena had
been observed by [31], but with a small shift of the range of Reynolds numbers.
The vorticity and pressure contours for one of the cases in the lock-in regime
is shown in Fig. 4.21. The expected 2S mode pattern in Fig. 4.21b, c(two single
vortexes per cycle of motion) can be seen in the wake [68]. Also, the smoothness of
the vorticity contours and pressure countours near the cylinder is a farther indication
that the method properly captures the complex dynamics of the flow near near the
110
Displacements Force coefficients
(a) Re = 90
(a) Re = 95
(b) Re = 100
(c) Re = 110
(d) Re = 120
Figure 4.20: Time history of the displacements for the freely vibrating cylinder in
the cross-stream direction: stable periodic phase.
111
immersed interface.
The comparison of CD, CrmsD , and CrmsL between the fixed cylinder cases and
the freely oscillating cylinder cases is shown in Fig. 4.23. As the Reynolds number
increases, the rms values of the drag and lift coefficients increase slowly, and the
mean values of the drag coefficient decrease slowly. However, for the freely oscillating
cylinder cases, in the lock-in regime, a significant different plot is presented, while
outside the lock-in regime, there is little difference between the fixed and the freely
oscillating cases.
The maximum oscillation amplitudes and frequencies are shown in Fig. 4.24.
The corresponding experimental results [2] and previous computations [53, 31] have
been added for comparison. For all simulations the maximum vibration amplitude
in the lock-in regime is lower that one in the experiment, where also the critical
Reynolds number for lock-in is higher (approximately Re = 103, while in most com-
putations it is below 95). Schulz & Kallinderis [53] speculated that this difference
could be due to three-dimensional effects induced by the absence of end-plates in
the experimental apparatus. Discrepancies can also be observed between the simu-
lations. Our results are in good agreement with the computations reported in [31],
where a boundary conforming methodology utilizing a moving frame of reference
was used. This indicates the accuracy of the proposed methodology. An extensive
grid refinement study, however, remains to be performed to clarify this issue.
112
(a) Re = 90
(b) Re = 95
(c) Re = 100
(d) Re = 110
(e) Re = 120
Figure 4.21: Freely vibrating cylinder in the cross-stream direction: instantaneous
vorticity contours.
113
(a) Re = 90
(b) Re = 95
(c) Re = 100
(d) Re = 110
(e) Re = 120
Figure 4.22: Freely vibrating cylinder in the cross-stream direction: instantaneous
pressure contours.
114
Figure 4.23: Comparison of the force coefficients for the cases of fixed cylinder and
the freely oscillating cylinder. ? fixed cylinder; a50 freely oscillating cylinder
115
Figure 4.24: Comparison of the oscillation frequency and the maximum oscillation
amplitude. ? present computations; diamondmath experiment in [2]; a50 computation in [53]; ?
computation in [31]
116
4.3.2 Two Degrees of Freedom
Although the one degree of freedom case has shed light on the physics of the
vortex-induced vibration problems, there are questions to be answered regarding
to what extent the freedom of oscillating in streamwise direction can change the
response of the system. Recently, a paper by Jauvtis & Williamson [27] reported
a substantial increase of the oscillation amplitude for systems with two degrees of
freedom when the mass ratio is smaller than a critical value in their water channel
experiments with Reynolds numbers ranged from 1,000 to 15,000. In the laminar
flow regime, whether such a dramatical change exists or not is still in question since
there are no experiment results on systems with two degrees of freedom reported
in the literature for low Reynolds numbers. Zhou et al. [74] carried out numerical
simulations using the discrete vortex method on vortex-induced vibrations of a cir-
cular cylinder with two degrees of freedom at Reynolds number 200. They compared
these results to those from one degree of freedom cases, and found no substantial
changes.
Here we shall not aim at reproducing the results in [27] as the moderate
Reynolds number flows in their experiments require fully three dimensional numeri-
cal simulations. By contrast, we shall stay in the low Reynolds number, laminar flow
regime, and select a low mass ratio n = 2.04 and a low damping ratio ? = 0.00425
for the system, which are in the range used in [27]. The Reynolds number is set
to Re = 200. The reduced velocity Ured, however, was varied from 1 to 11. The
computational domain is 36D?20D, in the streamwise and cross-stream directions,
117
respectively. The cylinder is located at a distance of 6D from the the inlet boundary,
where a uniform velocity is specified, and 10D from the upper and lower bound-
aries, where radiative conditions are imposed. At the outflow boundary a convective
boundary condition is used [41]. The number of grid points in all computations was
400 ? 360 in the streamwise and transverse directions, respectively. The grid was
stretched in both directions and the resulting resolution near the cylinder was ap-
proximately 0.02D ? 0.02D, which is the same as the one degree of freedom case
above. The time step was set to ?t = 0.0025D/U (forUred = 4.08, ?t = 0.002D/U).
As for the previous case, in all computations the corresponding stationary
cylinder problem was first solved and then free vibrations in both directions were
allowed. Due to the low mass ratio in this case, a stable periodic state was achieved
much faster compared to the previous problem.
Fig. 4.25 and Fig. 4.26 show the instantaneous vorticity and pressure contours
for several cases with increasing Ured. Outside the lock-in regime, e.g., Ured = 1.96
and Ured = 10.71, the contour patterns are very similar to those of the fixed cylinder
cases. Whithin the lock-in regime, the vibrations of the cylinder modify the vortex
shedding patterns substantially. However, for the current low Reynolds number
cases, only the ?2S? mode is captured as in the one degree of freedom cases.
The X,Y trajectories of several cases are shown in Fig. 4.27. All these figures
illustrate highly periodical behaviors and symmetries in the transverse direction.
They all exhibit the well-known ?figure eight? type of motion, while the different
shapes reflect the differences of phase angles between the X and Y displacements. In
the lock-in regime, very high amplitudes are observed with a maximum value of 0.597
118
(a) Ured = 1.96
(b) Ured = 4.08
(c) Ured = 4.76
(d) Ured = 10.71
Figure 4.25: Vortex-induced vibration of an elastic cylinder with two degrees of
freedom: instantaneous vorticity contours.
119
(a) Ured = 1.96
(b) Ured = 4.08
(c) Ured = 4.76
(d) Ured = 10.71
Figure 4.26: Vortex-induced vibration of an elastic cylinder with two degrees of
freedom: instantaneous pressure contours.
120
at Ured = 4.08, which agrees with the saturation maximum amplitude for laminar
flows (A ? 0.6) reported in [27] very well. They reported a crescent trajectory for
the ?super-upper? branch for their moderate Reynolds number experiments, which
is not presented in our low Reynolds number simulations.
Fig. 4.28 shows the time history of the displacements and lift coefficients for
several cases with increasing Ured, and Fig. 4.29 shows the corresponding phase
plots of displacement and lift coefficient. Outside the lock-in regime, both the force
and the displacement are very close to sinusoidal; while in the lock-in regime, the
shape of displacement doesn?t vary much from sinusoidal, but the force is quite
different from it, which illustrates a superposed oscillating frequency along with the
fundamental oscillating frequency.
The phase plots of drag and lift coefficients are presented in Fig. 4.30. Outside
the lock-in regime, the phase plots give the similar ?figure eight? as phase plots of
the streamwise and transverse displacement. However, the phase plots for the cases
in the lock-in regime are much more complex and resemble the corresponding plots
of the ?super-upper? branch discussed in [27] for much higher Reynolds numbers.
The intricate shapes indicate the development of oscillating frequencies different
from the fundamental oscillating frequencies in both the streamwise and transverse
directions.
Fig. 4.31 shows the oscillation frequency, transverse oscillation amplitude, and
the streamwise oscillation amplitude as functions of the reduced velocity. A lock-
in regime ranging from Ured ? 4 ? 6 is evident as the transverse and streamwise
oscillation amplitudes are prominent in this range. Also, the oscillation frequencies
121
(a) Ured = 1.96 (b) Ured = 4.08
(c) Ured = 4.76 (d) Ured = 10.71
Figure 4.27: X-Y phase plots. X and Y are in different scales.
122
Transverse displacements Transverse force coefficients
(a) Ured = 1.96
(b) Ured = 4.08
(c) Ured = 4.76
(d) Ured = 10.71
Figure 4.28: Time history of the displacements and lift coefficients.
123
(a) Ured = 1.96 (b) Ured = 4.08
(c) Ured = 4.76 (d) Ured = 10.71
Figure 4.29: Phase plots of transverse displacement, lift coefficient.
124
(a) Ured = 1.96 (b) Ured = 4.08
(c) Ured = 4.76 (d) Ured = 10.71
Figure 4.30: Phase plots of drag, lift coefficients.
125
are about the same values as the natural oscillation frequencies of the structure. We
also put the results in terms of reduced velocity from the one degree of freedom case
in the previous section. There is an apparent overlapping lock-in regime. However,
the two degrees of freedom case has a wider lock-in range, which is probably a result
of the much lower mass ratio.
It is also interesting to compare CD, CrmsD , and CrmsL between the one degree-
of-freedom and the two degree-of-freedom cases, as shown in Fig. 4.32. Again, there
is significant similarity between these two cases. Just like the oscillation frequency
and amplitude data, the two degree-of-freedom case has wider lock-in regime than
that of the one degree-of-freedom case.
4.4 Summary
In this chapter, the formal second-order accuracy of the embedded boundary
method presented in the previous chapter has been demonstrated via a specially-
designed case: a periodic array of moving cylinders in a two-dimensional planar
channel. This case is characterized by rigid bodies moving through the computa-
tional domain, which presents an extra difficulty to body-fitted methods.
After the establishment of the formal accuracy, the study has been concen-
trated on the vortex-induced vibrations of a circular cylinder. First, the forced
vibrations are studied, i.e., the flow induced by the harmonic in-line oscillation of
a circular cylinder in a quiescent fluid, and the flow from a transversely oscillat-
ing cylinder in a free-stream. Both cases are well documented in the literature,
126
Figure 4.31: The oscillation frequency,transverse oscillation amplitude and stream-
wise oscillation amplitude as functions of reduced velocity. diamondmath two degrees of freedom
cases; ? one degree of freedom cases from Sec. 4.3.1; a50 one degree of freedom cases
from [53]
127
Figure 4.32: Comparison of force coefficients of the freely oscillating cylinders with
one degree of freedom and two degrees of freedom. a50 two degrees of freedom cases;
? one degree of freedom cases from Sec. 4.3.1
128
and experimental and numerical results are available for comparison. The pro-
posed method was found to reproduce all features of the flow, including the detailed
force distribution on the body with the same degree of accuracy as the reference
boundary-conforming computations.
Then the proposed strong coupling scheme for fluid-structure interaction prob-
lems are verified through the studies of vortex-induced vibrations of an elastic cir-
cular cylinder with one degree of freedom and two degrees of freedom. For the
one degree-of-freedom case, the lock-in regime is accurately captured. The compu-
tational results have been compared with experimental data and other simulations
using body-fitted grids. The current results are in good agreement with other simula-
tions and reasonable agreement with the experiment. For the two degree-of-freedom
case at Reynolds number 200, the results are similar to the one degree-of-freedom
case but with a wider lock-in regime in reduced velocity space and higher oscillation
amplitude in the lock-in regime.
129
Chapter 5
Three-Dimensional Applications
5.1 Flow Past a Sphere
The flow past the immersed bluff-body is of interest in many engineering ap-
plications. Here, the three-dimensional case of the flow past a sphere for a Reynolds
numbers ranging from 50 to 1000 will be studied using the method developed in this
work. The flow is steady and axisymmetric for Reynolds numbers up to 200. While
at Re = 250 the flow is still steady, the axisymmetry breaks down. At Re = 300
the flow is no longer steady and is dominated by periodic vortex shedding. For
Re = 1000, the flow becomes chaotic and transitions to turbulence.
5.1.1 Axisymmetric Flow
Cylindrical coordinates are used for a more efficient distribution of the grid
pointsforthesimulationsoftheflowpastasphere. ForthecasesofRe = 50,100,150,
and 200, (Re = ?UD/?, where U is the freestream velocity, D the diameter of the
sphere), the size of the computational box in the streamwise direction is 30D with
the sphere located in the middle, and 15D in the radial direction. To investigate the
influence of grid resolution on the results three different grids have been considered
for these low Reynolds number cases. Grid 1 involves 40 ? 40 ? 100 computa-
tional points in the radial, azimuthal, and streamwise directions, respectively. In
130
Figure 5.1: Part of grid 40?40?100 for axisymmetric flow past the sphere. Uniform
grid in the red box.
the other two cases the grid is refined in the streamwise and radial directions (Grid
2: 80?40?200 and Grid 3: 160?40?400, respectively). A part of the computa-
tional grid 1 for the flow around the sphere is shown in Fig. 5.1, in which a uniform
grid distribution is adopted around the sphere. All three grids are stretched in the
streamwise and radial directions to cluster points near the surface of the sphere.
The resulting average grid spacing near the sphere for the three different resolutions
is approximately 0.1D, 0.05D and 0.025D, respectively. In all cases a uniform ve-
locity field is specified at the inflow plane. A convective boundary condition is used
at the outflow boundary [41], and radiative boundary conditions are applied at the
freestream boundary.
Table 5.1 shows the present results in comparison with the experimental results
131
Table 5.1: Prediction of CD for axisymmetric flow past the sphere
.
Re 50 100 150 200
Grid 1: 40?40?100 1.708 1.230 1.045 0.944
Grid 2: 80?40?200 1.610 1.118 0.920 0.807
Grid 3: 160?40?400 1.586 1.095 0.894 0.776
Experiment [10] 1.574 1.087 0.889 0.776
Reference Simulations [28] 1.575 1.100 0.900 0.775
in [10], and the well-resolved simulations by Johnson and Patel [28] where body-
fittedgridsareused. ForallReynoldsnumbersandgridresolutionsthemainfeatures
of the flow are properly captured. The drag coefficient on the coarsest grid, is a
little higher (approximately 8%) in comparison with the reference experimental and
numerical data. As the grid is refined, however, the agreement is very good. The
error in drag coefficient, which is computed using the experimental results in [10]
as the reference values, as a function of grid spacing near the sphere is shown in
Fig. 5.2. It can be seen that the error reduces with a second order slope, which
is consistent with the order of accuracy of the method. Streamlines for the these
Reynolds numbers are shown in Fig. 5.3. The flow direction is from left to right. A
steady separation bubble is formed behind the sphere for all four Reynolds numbers.
The flow is axisymmetric and has the same topological structure for all four cases.
The only differences are the separation angle, vortex center, and the length of the
recirculation bubble.
132
Figure 5.2: Error in the drag coefficient for the axisymmetric flow past the sphere
as a function of grid resolution.
133
The only non-zero vorticity component is in the azimuthal direction and its
contours are shown in Fig. 5.4. As the Reynolds number increases, the boundary
layer becomes thinner and the vorticity extends to further downstream. Also, it is
evident that the vorticity changes the sign across the separation point at the sphere
surface.
5.1.2 Planar-Symmetric Flow
For Re = 250 and 300 cases the computational domain is expanded in stream-
wise direction, becoming ?10D ? 40D with the sphere located at the origin, to
capture the vortical structure in the wake. And the grid is refined to contain
112?64?420 nodes in radial, azimuthal, and streamwise directions, respectively.
Compared to the grid used for lower Reynolds number calculations, this grid is
stretched more to cluster enough points near the sphere surface. With this res-
olution approximately 10 grid points are located in the boundary layers near the
stagnation point.
At Re = 250 and 300, the flow no longer retains its axial symmetry. The flow
past a sphere transits from aixsymmetric flow to three dimensional flow at about
Re = 211 through a linear instability [60]. However, a plane of symmetry does still
exist. In the current study, this plane of symmetry is almost (one grid cell difference)
perpendicular to the boundary plane of the azimuthal direction, which is the y ?z
plane in this work. Due to the use of ghost cells in the simulations, the plane of
symmetry doesn?t coincide with the y ? z plane exactly. Note that the results in
134
(a)
(b)
(c)
(d)
Figure 5.3: Axisymmetric streamlines past the sphere (colored by the streamwise
velocity): (a) Re = 50; (b) Re = 100; (c) Re = 150; (d) Re = 200.
135
(a)
(b)
(c)
(d)
Figure 5.4: Vorticity contours for axisymmetric flow past the sphere: (a) Re = 50;
(b) Re = 100; (c) Re = 150; (d) Re = 200.
136
(a)
(b)
Figure 5.5: Streamlines of projected velocity field for flow past the sphere at Re =
250 (colored by the streamwise velocity): (a) x?z plane; (b) y ?z plane.
this part have been rotated to have the symmetry plane coincide with the y ? z
plane. In [28], the plane of symmetry is near the azimuthal boundary plane. This
difference may be due to the different initial disturbances in two simulations.
For Re = 250, the streamlines obtained from the projected velocity field in
both x?z and y?z planes are shown in Fig. 5.5. The symmetry of the streamlines
in the x ? z plane is obvious with the upper part mirroring the lower part of the
figure. Note that in the y ?z plane the velocity component perpendicular to this
137
(a) (b)
Figure 5.6: Surface pressure distribution on the sphere at Re = 250: (a) front view;
(b) rear view.
plane is non-zero and the streamlines are different from the real three-dimensional
case due to the projection operation.
The contour plot of the surface pressure is shown in Fig. 5.6. The pressure
distribution near the stagnation points is still axisymmetric. However, on the rear
part of the sphere the pressure losses its axial symmetry. Note that the surface
pressure is obtained using the surface force calculation procedure in Chap. 3.
The contours of the azimuthal vorticity in the x?z and y?z planes are shown
in Fig. 5.7. It is evident that the vorticity in the x ? z plane are symmetric and
very similar to those of the lower Reynolds number cases. However, this symmetry
breaks down in the y ?z plane and one side of the vorticity contours is extended
farther downstream than the other side.
138
(a)
(b)
Figure 5.7: Vorticity contours for flow past the sphere at Re = 250: (a) Spanwise
vorticity in x?z plane; (b) Spanwise vorticity in y ?z plane.
The drag coefficient for Re = 250 is CD = 0.70 from the present simulation,
which is the same as in [28]. Due to the loss of axisymmetry, a non-zero lift force is
generated. And the lift coefficient is CL = 0.062, which is also the same as in [28].
For Re = 300 the flow past the sphere becomes unsteady and the vortices
are shed, although the planar-symmetry is retained. The Strouhal number in the
current study is St = 0.133, which is in good agreement with the Strouhal number
of 0.137 in [28] (within 3%). Fig. 5.8 shows the force coefficients as a function
of time. The mean drag coefficient is CD = 0.655, which is the same as that in
[28]. Note that in the figure the side force coefficients Cx and Cy are given directly
without a rotation to have y?z plane as the plane of symmetry. The combination
139
of these two components gives a lift coefficient CL = 0.064, which is again in very
good agreement with other reference data (CL = 0.065 in [28]).
The contours of the azimuthal vorticity component in the x?z plane is shown
in Fig. 5.9 for the four quarter periods. The definition of the period of shedding is
chosen to be the same as in [28]. From the first to the fourth panel, the pinch-off
and the convection downstream of a segment of vorticity is evident. A further check
of Fig. 5.10, which shows the azimuthal vorticity component in the y ? z plane,
reveals that the above phenomenon corresponds to a small vortex shed from the
vortex attached to the lower part of the sphere. Also, the vortex shed from the
upper part of the sphere has an interesting finger shape.
The vortical structures in the wake are shown in Fig. 5.11 and 5.12, where iso-
surfaces of the second invariant of the velocity gradient tensor (or the Q-criterion)
are used for their visualization. Hairpin like vortices originating from the surface
of the sphere can be observed, with new vortical structures developing around their
legs as they are convected downstream. This behavior is nearly the same as the one
shown in [28].
5.1.3 Transitional Flow
The flow past the sphere loses its planar-symmetry as the Reynolds number
further increases. A DNS of the flow at Re = 1000 has been conducted to estab-
lish the applicability of the current method in transitional and turbulent flows in
cylindrical coordinates. The computational domain is chosen to the same as in [60]:
140
Figure 5.8: Time history of drag, x, and y side force coefficients at Re = 300.
?4.5D ? 25D in the streamwise direction with the sphere located at the origin,
and 0 ? 4.5D in radial directions. The computational grid is 300 ? 64 ? 768 in
the radial, azimuthal, and streamwise directions, respectively. The total number
of points is nearly 15 million. The simulation is started from the interpolated flow
field of Re = 300. Due to the limited computing resources, this simulation has not
been run for enough time to obtain statistically stationary results: the simulation
has been carried out for a total of about 100 time units with the first 60 time units
discarded. For the current grid, a mean separation angle 100? is obtained, while in
[60] the averaged separation angle was 102?. Although our result is close to that
in [60], further refinement near the sphere surface may be necessary to capture the
thin boundary layers.
141
(a) t = 0T
(b) t = 14T
(c) t = 24T
(d) t = 34T
Figure 5.9: Spanwise vorticity contours for flow past the sphere at Re = 300 at the
four quarter periods in the x?z plane.
142
(a) t = 0T
(b) t = 14T
(c) t = 24T
(d) t = 34T
Figure 5.10: Spanwise vorticity contours for flow past the sphere at Re = 300 at the
four quarter periods in the y ?z plane.
143
(a) t = 0T
(b) t = 14T
(c) t = 24T
(d) t = 34T
Figure 5.11: Instantaneous vortical structures for flow past the sphere at Re = 300
(colored with the streamwise vorticity) at the four quarter periods: x?z view.
144
(a) t = 0T
(b) t = 14T
(c) t = 24T
(d) t = 34T
Figure 5.12: Instantaneous vortical structures for flow past the sphere at Re = 300
(colored with the streamwise vorticity) at the four quarter periods: y ?z view.
145
(a)
(b)
Figure 5.13: Instantaneous azimuthal vorticity contours for flow past the sphere at
Re = 1000: (a) x?z plane; (b) y ?z plane.
At this Reynolds number, small scales are present in the wake and the flow
becomes chaotic further downstream. Those small scale structures are generated
by the breakdown of the large cylindrical vortical structures that are formed just
downstream of the sphere. The mechanism of this breakdown is due to a Kelvin-
Helmholtz-like instability of the shear layer developed from the boundary layer sep-
aration of the sphere. The instantaneous azimuthal vorticity contours in the x?z
and y?z plane and the instantaneous vortical structures are shown in Fig. 5.13 and
5.14, respectively. The roll-up of the shear layer, the development of large hairpin
structures, and the breakdown of larger structures into smaller scales can be clearly
observed. The flow becomes turbulent in the far wake. Those results are in very
good agreement with the DNS using a high-order spectral method in [60].
The time history of the drag, and side force coefficients in x and y directions
are shown in Fig. 5.15a. Unlike the Re = 300 case, which has a predominant side
146
(a)
(b)
(c)
Figure 5.14: Instantaneous vortical structures for flow past the sphere at Re = 1000
(colored with the streamwise vorticity): (a) x?z view; (b) y ?z view; (c) oblique
view.
147
(a) Force coefficients (b) Frequency spectrum of CD
Figure 5.15: Time history of drag, x, and y side force coefficients and frequency
spectrum of CD at Re = 1000.
force which gives rise to a net lift force, for Re = 1000 both side forces in the x
and y directions are oscillating around zero. A mean drag coefficient CD = 0.46
is obtained for the limited sample size. The frequency spectrum of CD is given in
5.15b. The Strouhal number is St1 = 0.183, which is smaller than that in [60]. The
reason for this difference may be the earlier separation in our simulation and the
limited sample.
Fig. 5.16 and 5.17 show the time series of the axial and streamwise velocity
components at r = 0.3D and z = 2.0D and their frequency spectra, respectively.
Both components show a major Strouhal number of St1 = 0.183; they also present
a second higher frequency at St2 = 0.336, which is again smaller than the corre-
sponding value of St2 = 0.35 in [60].
148
(a) (b)
Figure 5.16: Time history and frequency spectrum of ur at r = 0.3D, z = 2.0D.
(a) (b)
Figure 5.17: Time history and frequency spectrum of uz at r = 0.3D, z = 2.0D.
149
5.2 Transitional Flow Past an Airfoil
A more challenging test-case is the computation of transitional flow around
an Eppler E387 airfoil. The angle of attack is 10? and the chord Reynolds number
is Re = ?Uc/? = 10,000 with U the free-stream velocity and c the chord length of
the airfoil. The Reynolds number regime is characteristic of bird flight.
The computational domain used is ?10D ? 21D and ?10D ? 10D in the
streamwise and transverse direction, respectively. The airfoil is arranged such that
there is a 10D distance to the inflow, the upper and the lower boundaries, and
a 20D distance to the outflow boundary. This simulation is very demanding in
terms of computing resources. A LES using the Lagrangian dynamic model has
been performed on a grid with approximately 18?106 nodes (560?48?672 in the
transverse, spanwise, and streamwise directions, respectively). The length scale in
the spanwise direction is chosen to be 4c. A part of the computational grid is shown
in Fig. 5.18. The computation has been run at a constant CFL = 1.2, which gives
a timestep varying in 4 ? 8?10?4c/U.
The flow field is characterized by the massive separation as shown in Fig.
5.19 where the mean streamlines are shown. A primary recirculation bubble, which
starts from the laminar separation point and extends until the trailing edge can be
observed. A second recirculation zone induced by the flow below the airfoil develops
at the trailing edge. Interestingly, a third separation bubble, which is attached to
the airfoil, beneath the large one can be observed. A careful examination of Fig.
5.20, which gives the pressure coefficient Cp on the upper and lower surface of the
150
Figure 5.18: Part of computational grid in x ? z plane around the airfoil. Every
three grid lines are shown.
Figure 5.19: Mean streamlines of the flow past an E387 airfoil. Streamlines are
colored with streamwise velocity component.
151
Figure 5.20: Mean pressure coefficient on the airfoil surface.
airfoil, shows that an adverse pressure gradient developes on the upper surface as
the flow approaching this second separation point. In Fig. 5.20, a small bump
of Cp is evident as a consequence of reattachment. The time history of the drag
and lift force coefficients is given in Fig. 5.21. A mean drag coefficient CD = 0.13
and a mean lift coefficient CL = 0.86 can be derived from the limited time sample.
The frequency sepctra of the drag and lift coefficients are shown in Fig. 5.22. The
lower Strouhal number is St1 = 0.05 and the averaged higher Strouhal number is
St2 = 1.3. These two different Strouhal numbers may be associated with the two
different vortex sheddings: one resulting from the separation on top of the airfoil
and the other from the trailing edge separation. However, a detailed analysis of the
time series of the variables at different locations is necessary to obtain the timing of
152
Figure 5.21: Time history of drag and lift coefficients for the flow past an airfoil.
the vortex shedding.
The turbulent intensities, turbulent shear stress, and the eddy viscosity are
shown in Fig. 5.23 (Note that the turbulent intensities and shear stress only contain
the resolved parts). The turbulence activities of the transverse and spanwise com-
ponents exist mainly within the shear layer formed at the trailing edge, where the
laminar flow beneath the airfoil interacts with the transitional flow resulting from
the massive separation. The close-up views of streamwise turbulent intensity, the
turbulent shear stress, and the mean eddy viscosity normalized using the molecu-
lar viscosity ?m are given in Fig. 5.24. It is apparent that a weak shear layer is
generated as the flow separates from the upper surface of the airfoil.
Fig. 5.25 shows the instantaneous spanwise and vorticity components, and Fig.
153
(a) (b)
Figure 5.22: Frequency spectra of the drag and lift coefficients for the flow past an
airfoil.
5.26 shows the instantaneous vortical structures. At the early stage of separation,
the vortices are well-organized and large rollers are generated. The flow becomes
turbulent very soon after it separates from the trailing edge.
5.3 Turbulent Flow Over a Traveling Wavy wall
In this section we will investigate the accuracy and efficiency of the proposed
approach in the framework of LES of turbulent flows. In particular, we will perform
LES for the case of turbulent flow over a flexible wall undergoing streamwise trav-
eling wave motion. A key parameter in this flow is the ratio of the wave speed, c,
to the external mean flow, U, which has a substantial effect on turbulence produc-
tion near the surface. Depending the value of c/U turbulence can be enhanced or
fully suppressed leading to relaminarization. The parametric space in our compu-
154
(a) < uu > /U2
(b) < vv > /U2
(c) < ww > /U2
(d) < uw > /U2
(e) < ?t > /?m
Figure 5.23: Turbulent intensities, turbulent shear stress, and mean eddy viscosity
for the flow past an airfoil.
155
(a) < ww > /U2
(b) < uw > /U2
(c) < ?t > /?m
Figure 5.24: Close-up views of the streamwise turbulent intensity, turbulent shear
stress < uw >, and mean eddy viscosity < ?t > for the flow past an airfoil.
156
(a) ?y
(b) ?z
Figure 5.25: Instantaneous contours of the spanwise and streamwise vorticity com-
ponents.
Figure 5.26: Instantaneous vortical structures using isosurfaces of Q = 10 colored
by the streamwise vorticity for the flow past an airfoil.
157
Figure 5.27: A sketch of the geometry for the case of turbulent flow over a traveling
wavy wall.
tations is selected to match the conditions in the DNS by Shen et al. [54]. These
simulations were performed using a pseudospectral method in the horizontal di-
rections, and 2nd order finite-differences in the wall-normal direction. They also
used a coordinate transformation to map the physical space into a computational
space that eliminates the need to deform the grid with the wave motion. A sketch
of the computational domain is shown in Fig. 5.27. The wavy wall is undergoing
a vertical oscillation in the form of a streamwise-traveling wave. The location of
the wall boundary as a function of time is: yw(t) = asink(x?ct), where a is the
magnitude of the oscillation, k = 2pi/? is the wavenumber with ? the wavelength,
and c is the phase-speed of the traveling wave. For such a configuration the only
non-zero component of the velocity vector on the boundary is the vertical one given
by: vw(t) = ??acosk(x ? ct), where ? = kc is the frequency of oscillation. The
components of the velocity vector are denoted by u, v, and w in the x (streamwise),
y (wall-normal), and z (spanwise) directions respectively.
158
(a) c/U = 0.0
(b) c/U = 0.4
(c) c/U = 1.2
Figure 5.28: Mean streamline pattern in a frame of reference following the traveling
wave. (a) c/U = 0.0; (b) c/U = 0.4; (c) c/U = 1.2. ? indicates the center of the
corresponding recirculation region in the reference DNS [54].
159
Wehaveconsideredthesamewavesteepness, ka = 0.25, andReynoldsnumber,
Re = ?U?/? = 10,170 (U is the mean freestream velocity), as in the DNS by Shen et
al. [54]. Computations at three different c/U ratios are presented: c/U = 0.0,
c/U = 0.4 and c/U = 1.2. The first case, where c/U = 0.0, corresponds to a sta-
tionary wavy wall, and it is included to highlight the effect of the wave motion on the
flow physics. The size of the computational domain is 2??2/pi??2?, in the stream-
wise, wall-normal and spanwise directions respectively. The corresponding number
of points in each direction is 288 ? 88 ? 64 for all simulations. The grid is uni-
form in the streamwise and spanwise directions, and is stretched in the wall-normal
direction. Grid refinement studies for the case of a stationary wavy boundary un-
der similar flow conditions reported in [4] guided the generation of the above grid.
Following Shen et al. [54], periodic boundary conditions are used in the stream-
wise and spanwise directions and a slip-wall at the top boundary. Phase-averaged
statistics with respect to the motion of the wavy boundary were accumulated for
approximately 40?/U, according to the relation [54]: x?ct = xprime + n?, where n is
an integer and 0 ? xprime < ?. In Fig. 5.28 the mean streamline pattern in a frame
of reference moving with the phase speed c is shown for all c/U ratios. For the
stationary wall case (see Fig. 5.28a) the flow separates just after the crest and forms
a large recirculating bubble on the downslope side of the wavy boundary. As the
c/U increases, the recirculating region becomes more elongated and extends over
the crest (see Fig. 5.28b). For larger values flow separation is altogether suppressed
and the streamlines follow the wavy terrain as shown in Fig. 5.28c for c/U = 1.2.
On the same figure the centers of the corresponding recirculation regions from the
160
Figure 5.29: Phase-averaged statistics for c/U = 0.4.
161
reference DNS [54] are shown for comparison. The agreement with our results is
very good indicating that the mean flow dynamics are properly captured with the
present method.
In Figs. 5.29 and 5.30, detailed phase-averaged statistics in a fixed frame of
reference for c/U = 0.4 and c/U = 1.2 are shown respectively. The differences
compared to the stationary wavy wall case are more evident in this reference frame,
with streamlines that begin from, and end on, the boundary as shown in part (a)
of both figures. From the isolines of the wall-normal velocity component, < v >,
it is also evident that as c/U increases from 0.4 to 1.2, the vertical flow induced
by the moving wall increases substantially (see Figs. 5.29c and 5.30c). This has an
effect on the shape of the streamlines which is concave for c/U = 0 (see Fig. 5.28a),
and becomes flat for c/U = 0.4 and convex for c/U = 1.2. As expected, c/U,
also affects the turbulent fluctuations. The turbulent kinetic energy is dramatically
decreased and the maximum is relocated further downstream when c/U increases
from 0.4 to 1.2, as shown in part (d) of both figures. A comparison of Figs. 5.29 and
5.30 to the corresponding ones in [54] -the scale and contour intervals are selected
appropriately to facilitate direct comparisons- shows very good agreement with the
reference results.
To further investigate the accuracy of the proposed approach the forces on the
wavy boundary are computed, and then compared to the corresponding reference
data. These quantities are very sensitive to the adopted numerical methodology and
are a stringent test for the accuracy of the proposed formulation. Following Shen et
al. [54], the local friction force, ffx, and local pressure force, fpx, on an element, ds,
162
Figure 5.30: Phase-averaged statistics for c/U = 1.2.
163
Figure 5.31: Mean force acting on the traveling wavy boundary as a function of
c/U. , ? is total friction force, Ff; , a50 is the total pressure force, Fp. Lines
are the DNS results in [54], and symbols are the present results.
on the solid boundary can be written as:
braceleftBigg ff
x = ?
bracketleftbigg
?2?u?xdywdx +
parenleftbigg?u
?y +
?v
?x
parenrightbiggbracketrightbigg 1
ds,
fpx = pdywdx 1ds.
(5.1)
The velocity derivatives and surface pressure in Eq. (5.1) are computed using the
method discussed in section 3.4. Then, one can integrate ffx and fpx over the wavy
surfacetoobtainthetotalfrictionforce, Ff, andtotalpressureforce, Fp, respectively.
A comparison of these quantities with the reference data is shown in Fig. 5.31. Note
that we also include a simulation at c/U = 0.8, where all other parameters are the
same as in the cases discussed above. Our results are in excellent agreement with
DNS results. As in the reference simulations the variation of the friction force is
relatively small and the value is always above zero, while the pressure force decreases
164
as c/U increases and becomes negative around c/U = 1.0.
Finally, a series of instantaneous realizations from all cases were carefully
examined to verify the smoothness of the velocity, vorticity and pressure fields near
the moving boundary, and visualize the characteristic coherent vortical structures
present in this flow. An example is shown in Fig. 5.32 where isosurfaces of the
second invariant of the velocity gradient tensor, Q, are used to identify the near
wall vortices. Three cases are included: c/U = 0.0, c/U = 0.4, and c/U = 1.2.
For all cases the isosurfaces are colored with the streamwise vorticity values. For
the stationary wavy boundary (c/U = 0.0) the presence of strong quasi-streamwise
vortices that begin in the upslope area of the wave and extend over the trough is
evident. As c/U increases, there are fewer and more elongated vortices. For the
maximum ratio considered in the present study (c/U = 1.2) these structures are
fully suppressed. This trend is in agreement with the phase-averaged turbulent
statistics shown in Figs. 5.29 and 5.30, and the results reported in [54].
5.4 Pulsatile Flow Past a Prosthetic Bileaflet Heart Valve
To further demonstrate the capabilities of the present method in handling real-
isticturbulentandtransitionalflowproblemsthatinvolvecomplexthree-dimensional
boundaries consisting of multiple moving parts, we have computed the flow past a
bileaflet, mechanical heart valve in the aortic position. The geometry of the valve
is shown in Fig. 5.33. The shape and size of the leaflets roughly mimics the St.
Jude Medical (SJM) Standard bi-leaflet mechanical heart valve, which is commonly
165
(a) c/U = 0.0
(b) c/U = 0.4
(c) c/U = 1.2
Figure 5.32: Isosurfaces of Q = 8.0 colored with the streamwise vorticity at an
instant in time.
166
Figure 5.33: Geometry of the simplified bileaflet prosthetic heart valve placed in
a rigid wall model of the aortic root. The leaflets are shown in their fully open
position, and part of the housing wall has been removed for clarity. Two planes of
the Cartesian grid are also shown.
used in clinical practice. In our case, however, the hinge region has been simplified
to only allow one degree of freedom for each leaflet (rotation around a fixed axis
in the y direction). Also, the leaflets contact the housing wall tightly and there is
no gap between the hinge and the housing wall. The housing is also shown in Fig.
5.33, where a straight pipe with rigid walls expands and then contracts to mimic the
geometry of the aortic root. The overall set-up resembles the ones commonly used
in in-vitro experiments to test the hemodynamic performance of prosthetic valves
in the aortic position (see for example [72] for a recent review).
167
In the actual case the motion of the leaflets is a product of their interaction
with the pulsatile blood flow: the valve opens at the beginning of the systole and
closes before the start of the diastole. When the leaflets are fully open the blood
flow rapidly accelerates through the valve and reaches its peak velocity. Since the
objective of the present simulations, however, is to demonstrate the applicability of
themethodtocomplexmovingboundaries, themovementoftheleafletsisprescribed
according to a simplified law that resembles the real movement of the leaflets as
determined by their interaction with the blood flow. The fully open position of the
leaflet forms a 85? angle with the y ?z plane, while a rotation of 53? towards the
housing wall closes the valves completely. In Fig. 5.34 the variation of the bulk
velocity and the corresponding angle of the leaflets are shown for a full pulsatile
cycle. The Reynolds number based on the pipe diameter and the maximum bulk
velocity was set to 4,000, which is a realistic value for this application.
The above configuration which involves the interaction of the flow with a
stationary (the aortic chamber) and two moving (the leaflets) boundaries is a great
challenge for any structured or unstructured boundary-conforming method. In the
present computation the overall geometry is immersed into a Cartesian grid and the
boundary conditions on the stationary and moving boundaries are imposed using
the method described in the previous sections. The computational grid involves
420 ? 200 ? 200 nodes in streamwise, spanwise (direction parallel to the rotation
axis of the leaflets), and transverse directions, respectively. The mesh is stretched in
the x?y plane and is kept uniform in the z direction (see Fig. 5.33), in a way that
the leaflets are moving in an approximately isotropic grid with average dimension of
168
Figure 5.34: Imposed flow rate and leaflet kinematics for the heart valve problem:
(a) Bulk velocity at the inlet; (b) opening angle of the leaflets, as functions of time.
0.005D. At the inflow boundary, located 4D upstream of the valve, fully developed
flow corresponding to the flowrate variation shown in Fig. 5.34 is imposed. At the
outflow boundary (8D downstream of the valve) a convective boundary condition is
used. The total number of nodes in this computation is 16.8 million.
The flow in the proximal area of the leaflets is very complicated, and it is
dominated by intricate vortex-leaflet and vortex-vortex interactions. To illuminate
basic flow patterns and demonstrate that the present method can accurately capture
the thin shear layers that form on both moving and stationary immersed boundaries,
instantaneous contours of the streamwise velocity and spanwise vorticity are shown
in Fig. 5.35 for a few characteristic phase angles during the pulsatile cycle. When the
leaflets are fully open, a strong jet is formed from the central orifice which interacts
with the wake of the leaflets as can be seen from the velocity and vorticity contours
169
(a) t/T = 0.2
(b) t/T = 0.4
(c) t/T = 0.6
(d) t/T = 0.8
Figure 5.35: Instantaneous flow fields at an x?z plane cutting through the middle
of the leaflets. Left: streamwise velocity (?1.5 < u/Umaxb < 3, intervals are 0.05);
Right: spanwise vorticity (?20 < ?yD/Umaxb < 20, intervals are 0.2). The phase
reference, t/T, is from Fig. 5.34.
170
in Fig. 5.35a. Two shear layers are generated and roll-up into large vortices which
interact with other vortices inside the chamber. A similar flow pattern has been
observed downstream of a bi-leaflet valve in a recent highly resolved PIV experiment
where a much more detailed model of the left ventricle was used [47]. There are
several other locales where weaker shear layers are present: at the edge of the leaflet
housing and outer surfaces of the leaflets, for example. In general, at this stage in
the pulsatile cycle the flow has similar characteristics with what has been observed
in steady experiments with the leaflets at the fully open position [8].
As the flow rate is decreasing and the leaflets begin to move towards the fully
closed position ( see Fig. 5.35b), the central orifice shrinks and the jet becomes
thinner. The interaction of this jet and the ambient decelerating flow generates
stronger shear layers than before. The flow rate keeps decreasing until it reaches
zero, and the leaflets reach their maximum closing angle. The central jet during
this stage becomes weaker and is finally dissipated as there is no flow through the
central orifice. When the leaflets begin to open again and the flow rate starts to
increase, the ring vortex at the edge of leaflet housing begins its development (see
Fig. 5.35c). The flow accelerates until it reaches the maximum flow rate, and the
leaflets reach their fully open position again as shown in Fig. 5.35d.
To better illuminated the highly three-dimensional nature of the complex vor-
tex interactions just downstream of the leaflets, isosurfaces of Q are shown in
Fig 5.36. In (a) and (b) of Fig. 5.36, two consecutive instantaneous realizations
(t/T = 0.3 and t/T = 0.4) correspond to instances in the process of the closing
of the leaflets. the motion of the leaflets generate strong vortices (structures B) at
171
(a)
(b)
Figure 5.36: (a) and (b). For caption see next page.
172
(c)
(d)
Figure 5.36: Instantaneous vortical structures visualized using isosurfaces of Q =
600 colored by the streamwise vorticity at: (a) t/T = 0.3; (b) t/T = 0.4; (c)
t/T = 0.8; (d) t/T = 0.9. The phase reference, t/T, is from Fig. 5.34.
173
their tips. On the other hand, after the leaflets return to their fully open position,
as the flow rate increases a ring vortex is generated at the edge of the leaflet housing
Fig 5.36c (t/T = 0.8). At a later time (Fig. 5.36d, t/T = 0.9) this ring vortex is
shed from the edge and interact with the surrounding flow, which again gives rise
to complex vortex-vortex interactions.
174
Chapter 6
Conclusions and Future Directions
The objective of the present study is the development of a highly accurate,
cost-effective strategy for the numerical simulations of a wide range of challenging
problems in biological flows and other fields. Brief conclusions and further possible
developement of current work are given in the chapter.
6.1 Conclusions
First, the original Cartesian solver [3] is extended to a generalized basic solver
for the Navier-Stokes equations in both Cartesian and cylindrical coordinates. The
singularity at the centerline is addressed and a simple approach is adopted. To
remove the severe timestep constraint in the azimuthal direction near the centerline,
the diffusive terms in the azimuthal direction in the momentum equations are treated
implicitly. The solver is validated through a LES of fully developed turbulent pipe
flow. Good agreement with the reference DNS data has been obtained.
Despite the increased attention that a variety of immersed boundary methods
have been receiving in recent years, systematic studies examining their accuracy and
applicability in cases of moving boundaries have not been reported. In the present
study we initially identified problems that are unique to these formulations, and are
a consequence of the boundary motion on a fixed grid. In particular, we have shown
175
that in a typical fractional-step scheme , when a point near the interface (forcing
point)movesintotheflow(fluidpoint), carrieswithitunphysicalinformationrelated
to the pressure and velocity derivatives. As a result, the evaluation of the RHS of the
momentum equations at this splitting cycle is problematic. To address this issue
we proposed a robust field-extension procedure, where the velocity and pressure
fields are ?extended? in the solid phase at the end of each substep. This way the
pressure and velocity derivatives at the problematic points have values that reflect
the proper boundary conditions. The overall extension procedure, although it is
tailored to our reconstruction scheme that utilizes the well defined normal to the
interface [4], is applicable in a straightforward manner to any embedded-boundary
formulation that reconstructs the solution around points in the fluid phase [15]. In
cases where some type of a generalized ghost-cell approach is used [29, 61], then the
grid points emerging from the solid need to be treated in a manner similar to that
suggested, for example, in [63].
There are a variety of fluid-structure interaction problems that arise in engi-
neering, biology and medicine. The numerical simulations of these problems require
the solution of the Navier-Stokes equations for the fluid flows in a changing domain
and the coupled dynamical equations for the structures that may be undergoing large
deformation and/or displacement. The coupling between the flow and the structure
leads to a highly nonlinear system, which presents a great challenge to numerical
algorithms on accuracy, robustness, and efficiency. We have proposed a strong cou-
pling scheme, where the fluid and the structure are treated as elements of a single
dynamical system, and all of the governing equations are integrated simultaneously
176
and interactively in the time domain. This algorithm utilizes an iterative predictor-
corrector scheme that accounts for the interaction between the hydrodynamic loads
and the motion of the structure.
In terms of validating the proposed methodologies, a variety of problems with
increasing complexity were considered. First, the formal accuracy of the method
(2nd order) was established for the problem of an array of cylinders moving in a
plane channel. Then, computations for two cases of laminar flow interacting with
moving boundaries were conducted. The flow induced by the harmonic in-line os-
cillation of a circular cylinder in a quiescent fluid, and the flow from a transversely
oscillating cylinder in a free-stream. Both cases are well documented in the liter-
ature, and experimental and numerical results are available for comparison. The
proposed method was found to reproduce all features of the flow, including the
detailed force distribution on the body with the same degree of accuracy as the ref-
erence boundary-conforming computations. In addition, vortex-induced vibrations
of an elastic cylinder with one and two degrees of freedom in a free-stream are simu-
lated. The results are compared with reference experiments and simulations. Good
agreement was observed.
Several three-dimensional cases have been considered. Flows involving station-
ary bodies were considered. The laminar flow past a sphere was examined carefully
and the results are in excellent agreement with the well-resolved DNS using a body-
fitted grid. Transitional flows past a stationary sphere and an airfoil were simulated.
Then, computations of the flow over a traveling wavy wall provided a validation of
the method in LES of turbulent flows with moving boundaries. Finally LES of tran-
177
sitional flow around a prosthetic heart valve with moving leaflets demonstrated the
robustness and applicability of the method in turbulent/transitional flows with mul-
tiple moving boundaries. Despite the simplification on the geometry the obtained
results revealed flow patterns that have been reported in the literature.
6.2 Future Directions
In general the results in this study show that embedded-boundary formulations
can be cost-efficient strategies especially for moving boundary problems without sac-
rificing accuracy. A limiting factor in all these methods, however, is the inflexibility
in clustering grid points near a complex body. As the Reynolds number increases
and more points need to be clustered near the body to resolve the thinner boundary
layers, this can result in very large grids. To investigate if they are still cost efficient
compared to boundary conforming formulations, one needs to look at the rate of
increase of the total number of grid points for both categories of methods and the
cost per node for each method as a function of the Reynolds number. This is a fairly
difficult task since it depends on a variety of problem specific parameters. In gen-
eral for the case of a single body, when using boundary-fitted grids, an increase in
the Reynolds number requires refinement of the mesh in the wall-normal direction
only; in non-boundary-conforming formulations like the present one, however, to
achieve refinement normal to the body the grid usually needs to be refined in two or
three directions. As a result the required number of grid points as a function of the
Reynolds number increases faster for non-boundary-conforming formulations com-
178
pared to boundary conforming ones. Actually, for laminar boundary layers it can be
shown that the ratio of the boundary-conforming to the non-boundary-conforming
grid size increases approximately as Re1.5 [37]. However, given that methods like
the present one are substantially less expensive per node than boundary-conforming
ones (especially in case of boundaries undergoing large motions and/or deformations
where grid re-generation is a major obstacle that compromises the accuracy and ef-
ficiency of boundary fitted approaches) they can be cost efficient for a fairly wide
range of Reynolds numbers. The computation of the flow around the heart valve,
for example, has been performed on a Linux desktop workstation equipped with
four AMD Opteron 846 2GHz processors (each with 1M cache, and 8GB memory
in total) and takes about 6 CPU hours for one pulsatile cycle. A 10 cycle simula-
tion with 17 million points can be completed in less than 3 days. For a variety of
problems, especially from biology where only low/moderate Reynolds numbers are
encountered, methods such the present one have a lot of advantages.
In the near future, there are three major directions in the further extension of
the current work: a) description of the geometries/evolution of immersed interfaces;
b) dynamical grid adaptation; and c) introduction of new physics phenomena and
models.
a) The immersed interfaces can be two and/or three dimensional; their motion
can be prescribed or governed by the interaction between different phases separated
by the interfaces; the interfaces can collapse or break up. For the description of
evolving interfaces, both explicit and implicit methods can be applied, while their
appropriateness in the framework of our method and in the context of the physical
179
problems studied needs further investigation.
b) As the immersed body/interface interacting with the flow field, the location
of interface becomes unpredictable and the dynamical grid adaptation will be im-
perative to a successful simulation. One key issue here is to extend the applicability
of our embedded boundary method and at the same time to maintain the efficiency
of Cartesian solvers.
c) Complex geometries and moving interfaces are encountered in many fields
where fluid flow plays a major role. Some can be solved readily in the current
framework with the addition of new boundary conditions at the interface, such as
free surface and interfacial flows, etc. However, new physical models are required
for most situations, e.g. fluid-structure interaction problems in biomimetic and
physiological flows, reacting flows, and aero-acoustics among others.
180
BIBLIOGRAPHY
[1] K. Akselvoll and P. Moin, ?Large-eddy simulation of turbulent confined coan-
nular jets and turbulent flow over a backword-facing step.? Report TF-63, Ther-
mosciences Div., Dept. Mech. Eng., Stanford University, 1995.
[2] P. Anagnostopoulos and P.W. Bearman, ?Response characteristics of a vortex-
excited cylinder at low Reynolds numbers.? J. Fluids Struct., 6:39?50, 1992.
[3] E. Balaras, ?Large eddy simulation of high Reynolds number wall bounded
shear layers.? Ph.D. Dissertation, Swiss Federal Institute of Technology, 1995.
[4] E. Balaras, ?Modeling complex boundaries using an external force field on fixed
Cartesian grids in large-eddy simulations.? Comput. Fluids, 33:375?404, 2004.
[5] E. Balaras, U. Piomelli, and J.M. Wallace, ?Self-similar states in turbulent
mixing layers.? J. Fluid Mech., 446:1?24, 2001.
[6] E. Balaras and J. Yang, ?Nonboundary conforming methods for large-eddy
simulations of biological flows.? ASME J. Fluids Eng., 127:851-857, 2005.
[7] S.A. Bayyuk, K.G. Powell, and B. van Leer, ?A simulation technique for 2-
D unsteady inviscid flows around arbitrary moving and deforming bodies of
arbitrary geometry.? AIAA Paper 93-3391-CP, 1993.
[8] P. Browne, A. Ramuzat, R. Saxena, and A.P. Yganathan, ?Experimental inves-
tigation of the steady flow downstream of the St. Jude bileaflet heart valve: a
181
comparison between LDV and PIV techniques.? Annals Biomed. Eng., 28:39?
47, 2000.
[9] B. Carnahan, H.A. Luther, and J.O. Wilkes, Applied Numerical Methods, John
Wiley & Sons, Inc., New York, 1969.
[10] R. Clift, J.R. Grace, and M.E. Weber, Bubbles, Drops and Particles, Academic
Press, San Diego, CA, 1978.
[11] G.S. Constantinescu and S.K. Lele, ?A highly Accurate Technique for the Treat-
ment of Flow Equations at the Polar Axis in Cylindrical Coordinates Using
Series Expansions.? J. Comput. Phys., 183:165?186, 2002.
[12] J.W. Deardorff, ?A numerical study of three-dimensional turbulent channel flow
at alrge Reynolds numbers.? J. Fluid Mech.. 41:453?480, 1970.
[13] H. D?utsch, F. Durst, S. Becker, and H. Lienhart, ?Low-Reynolds-number flow
around an oscillating circular cylinder at low Keulegan-Carpenter numbers.? J.
Fluid Mech.. 360:249?271, 1998.
[14] J.G. Eggels, F. Unger, M.H. Weiss, J. Westerweel, R.J. Adrian, R. Friedrich,
and F.T.M. Nieuwstadt, ?Fully developed turbulent pipe flow: a comparison
between direct numerical simulation and experiment.? J. Fluid Mech.. 268:175?
209, 1994
[15] E.A. Fadlun, R. Verzicco, P. Orlandi, and J. Mohd-Yusof, ?Combined
immersed-boundary finite-difference methods for three-dimensional complex
flow simulations.? J. Comput. Phys., 161:35?60, 2000.
182
[16] K. Fukagata, and N. Kasagi, ?Highly Energy-Conservative Finite Difference
Method for the Cylindrical Coordinate System.? J. Comput. Phys., 181:478?
498, 2002.
[17] M. Germano, U. Piomelli, P. Moin, and W.H. Cabot, ?A dynamic subgrid-scale
eddy viscosity model.? Phys. Fluids A, 3:1760?1765, 1991.
[18] M. Germano, ?Turbulence ? The filtering approach.? J. Fluid Mech., 238: 325?
336, 1992.
[19] D. Goldstein, R. Handler, and L. Sirovich, ?Modeling a no-slip flow boundary
with an external force field.? J. Comput. Phys., 105:354?366, 1993.
[20] D. Goldstein, R. Handler, and L. Sirovich, ?Direct numerical simulation of
turbulent flow over a modeled riblet covered surface.? J. Fluid Mech., 302:333?
, 1995.
[21] D. Goldstein, and T.C. Tuan, ?Secondary flow induced by riblets.? J. Fluid
Mech., 363:115?, 1998.
[22] W. Gu, C. Chyu, and D. Rockwell, ?Timing of vortex formation from an oscil-
lating cylinder.? Phys. Fluids, 6:3677?3682, 1994.
[23] E. Guilmineau, and P. Queutey, ?A Numerical simulation of vortex shedding
from an oscillating circular cylinder.? J. Fluids Struct., 16:773?794, 2002.
183
[24] H. H. Hu, N. A. Patankar, and M. Y. Zhu. ?Direct numerical simulations
of fluid-solid systems using the arbitrary Largrangian-Eulerian Technique.? J.
Comput. Phys., 169:427?462, 2001.
[25] J.C.R. Hunt, A.A. Wray, and P. Moin, ?Eddies, Streams and Convergence
Zones in Turbulent Flows.? , Proceedings of the Summer Program 1988, Center
for Turbulence Research, 193, 1988.
[26] G. Iaccarino, and R. Verzicco, ?Immersed boundary technique for turbulent
flow simulations.? Applied Mechanics Reviews, 56:331?347, 2003.
[27] N. Jauvtis, and C.H.K. Williamson, ?The effect of two degrees of freedom on
vortex-induced vibration at low mass and amping.? J. Fluid Mech., 509:23?62,
2004.
[28] T. A. Johnson and V. C. Patel, ?Flow past a sphere up to a Reynolds number
of 300.? J. Fluid Mech., 378:353?385, 1999.
[29] J. Kim, D.Kim, and H. Choi, ?An immersed-boundary finite-volume method
for simulations of flow in complex geometries.? J. Comput. Phys., 171:132?150,
2001.
[30] M.C. Lai and C.S. Peskin, ?An immersed boundary method with formal second-
order accuracy and reduced numerical viscosity.? J. Comput. Phys., 160:705?
719, 2000.
184
[31] L. Li, S.J. Sherwin, and P.W. Bearman. ?A moving frame of reference algo-
rithm for fluid/structure interaction of rotating and translating bodies.? Int. J.
Numer. Meth. Fluids, 38:187?206, 2002.
[32] D.K. Lilly. ?The representation of small-scale turbulence in numerical simula-
tion experiments.? Proceddings of the IBM Scientific Computing Symposium on
Environmental Sciences. Yorktown Heights, N.Y., 1967.
[33] D.K. Lilly, ?A proposed modification of the Germano subgrid-scale closure
method.? Phys. Fluids A, 4:633? , 1992
[34] T.S. Lund, ?On the use of discrete filters for large eddy simulation.? Annual
Research Briefs, 1997, Center for Turbulence Research, 83?95, 1997.
[35] S. Majumdar, G. Iaccarino, and P. Durbin, ?RANS solvers with adaptive struc-
tured boundary non-conforming grids.? Annual Research Briefs, 2001, Center
for Turbulence Research, 353?366, 2001.
[36] C. Meneveau, C.S. Lund, and W.H. Cabot, ?A Lagrangian dynamic subgrid-
scale model of turbulence.? J. Fluid Mech., 319:353?385, 1996.
[37] R. Mittal and J. Iaccarino, ?Immersed boundary methods.? Annu. Rev. Fluid
Mech., 37:239?261, 2005.
[38] Y. Morinishi, O.V. Vasilyev, and T. Ogi, ?Fully conservative finite difference
scheme in cylindrical coordinates for incompressible flow simulations.? J. Com-
put. Phys., 197:686?710, 2004.
185
[39] A.H. Nayfeh, S. Emam, S. Preidikman, and D.T. Mook. ?An Exact Solution
for the Natural Frequencies of Flexible Beams Undergoing Overall Motion.? J.
Vibration and Control, 9:1221?1229, 2003.
[40] D.J. Newman and G.E. Karniadakis. ?A direct numerical simulation study of
flow past a freely vibrating cable.? J. Fluid Mech., 344:95?136, 1997.
[41] I. Orlanski. ?Simple boundary-condition for unbounded hyperbolic flows.? J.
Comput. Phys., 21:251?269, 1976.
[42] S. Osher and R. Fedkiw. Level set methods and dynamic implicit surfaces.
Springer, New York, 2003.
[43] C.S. Peskin, ?Flow patterns around heart valves: a numerical method.? J.
Comput. Phys., 10:252?271, 1972.
[44] C.S. Peskin and D.M. McQueen, ?Numerical analysis of blood flow in the
heart.? J. Comput. Phys., 25:220?252, 1980.
[45] C.S. Peskin and D.M. McQueen, ?Cardiac fluid dynamics.? Critical Reviews in
Biomedical Engineering, 20:451?459, 1992.
[46] C.D. Pierce, ?Progress-variable approach for large eddy simulation of turbulent
combustion.? Ph.D. Dissertation, Stanford University, 2001.
[47] O. Pierrakos, P.P. Vlachos, and D.P. Telionis, ?Time-resolved DPIV analysis
of vortex dynamics in a left ventricular model through bileaflet mechanical and
porcine heart valve prostheses.? J. Biomech. Eng., 126:714?726, 2004.
186
[48] U. Piomelli, E. Balaras, and A. Pascarelli, ?Turbulent structures in accelerating
boundary layers.? J. Turbulence, 1:1?19, 2000.
[49] P. Ploumhans, G.S. Winckelmans, J.K. Salmon, A. Leonard, and M.S. War-
ren, ?Vortex methods for direct numerical simulation of three-dimensional bluff
body flows: Application to the sphere at Re=300, 500, and 1000.? J. Comput.
Phys., 178:427?463, 2002.
[50] S. Preidikman and D.T. Mook. ?Time-Domain Simulations of Linear and Non-
Linear Aeroelastic Behavior.? J. Vibration and Control, 6:1135?1176, 2000.
[51] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical
Recipes in Fortran 77: The Art of Scientific Computing. Cambridge University
Press, N.Y., 1992.
[52] J.J. Quirk, ?An alternative to unstructured grids for computing gas dynamic
flows around arbitrary two-dimensional bodies.? Comput. Fluids 23:125?, 1994.
[53] K.W. Schulz and Y. Kallinderis. ?Unsteady flow structure interaction for incom-
pressible flows using deformable hybrid grids.? J. Comput. Phys., 143:569?597,
1998.
[54] L. Shen, X. Zhang, D.K.P. Yue, and M.S. Triantafyllou. ?Turbulent flow over a
flexible wall undergoing a streamwise travelling wave motion.? J. Fluid Mech.,
484:197?221, 2003.
[55] D. Shiels, A. Leonard, and A. Roshko. ?Flow-induced vibration of a circular
cylinder at limiting structural parameters.? J. Fluids Struct., 15:3?21, 2001.
187
[56] J. Smagorinsky. ?General circulation experiments with the primitive equations.
I. The basic experiment.? Monthly Weather Review, 91:99?164, 1963.
[57] P.N. Swarztrauber, ?Direct method for the discrete solution of separable elliptic
equations.? SIAM J. Num. Anal., 11:1136?, 1974.
[58] P.N. Swarztrauber, ?Fast Fourier Transforms Algorithms for Vector Comput-
ers.? Parallel Computing. 45?63, 1984.
[59] F. Tessicini, G. Iaccarino, M. Fatica, M. Wang, and R. Verzicco, ?Wall mod-
eling for large-eddy simulation using an immersed boundary method.? Annual
Research Briefs 2002, Center for Turbulence Research, 181?378, 2002.
[60] A.G. Tomboulides and S.A. Orszag, ?Numerical investigation of transitional
and weak turbulent flow past a sphere.? J. Fluid Mech., 416:45?73, 2000.
[61] Y.-H. Tseng and J.H. Ferziger, ?A ghost-cell immersed boundary method for
flow in complex geometry.? J. Comput. Phys., 192:593?623, 2003.
[62] H.S. Udaykumar, R. Mittal, and W. Shyy, ?Computation of solid-liquid phase
fronts in the sharp interface limit on fixed grids.? J. Comput. Phys., 153:535?
574, 1999.
[63] H.S. Udaykumar, R. Mittal, P. Rampunggoon, and A. Khanna, ?A sharp inter-
face Cartesian grid method for simulating flows with complex moving bound-
aries.? J. Comput. Phys., 174:345?380, 2001.
188
[64] H.S. Udaykumar, W. Shyy, and M.M. Rao, ?Elafint: A mixed Eulerian-
Lagrangian method for fluid flows with complex and moving boundaries.? Int.
J. Numer. Methods Fluids 22:691?, 1996.
[65] R. Verzicco, J. Mohd-Yusof, P. Orlandi, and D. Haworth, ?Large-eddy simula-
tion in complex geometric configurations using body forces.? AIAA J., 38:427?,
2000.
[66] R. Verzicco and P. Orlandi, ?A finite-difference scheme for three-dimensional
incompressible flows in cylindrical coordinates.? J. Comput. Phys., 123:402?
414. 1996.
[67] C.H.K. Williamson and R. Govardhan. ?Vortex-induced vibrations.? Annu.
Rev. Fluid Mech., 36, 413?455, 2004.
[68] C.H.K. Williamson and A. Roshko. ?Vortex formation in the wake of an oscil-
lating cylinder.? J. Fluids Struct., 2:355?381, 1988.
[69] C.T. Yamamoto, J.R. Meneghini, F. Saltara, R.A. Fregonesi, and J.A. Ferrari
Jr., ?Numerical simulations of vortex-induced vibration on flexible cylinders.?
J. Fluids Struct., 19:467?489, 2004.
[70] J. Yang and E. Balaras. ?An embedded-boundary formulation for large-eddy
simulation of turbulent flows interacting with moving boundaries.? J. Comput.
Phys., 2006.
189
[71] T. Ye, R. Mittal, H.S. Udaykumar, and W. Shyy, ?An accurate Cartesian grid
method for viscous incompressible flows with complex immersed boundaries.?
J. Comput. Phys., 156:209?240, 1999.
[72] A.P. Yoganathan, Z. He, and C.S. Jones, ?Fluid mechanics of heart valves.?
Annu. Rev. Biomed. Eng., 6:331?362, 2004.
[73] D. De Zeeuw, and K.G. Powell, ?An adaptively refined Cartesian mesh solver
for the Euler equations.? AIAA Paper 91-1542-CP, 1991.
[74] C.Y. Zhou, R.M.C. So, and K. Lam. ?Vortex-induced vibrations of an elastic
circular cylinder.? J. Fluids Struct., 13:165?189, 1999.
190