
Full Citation 
Material Information 

Title: 
10.2.5  Numerical study on the motion characteristics of a freely falling twodimensional circular cylinder in a channel Particle Bubble and Drop Dynamics 

Series Title: 
7th International Conference on Multiphase Flow  ICMF 2010 Proceedings 

Physical Description: 
Conference Papers 

Creator: 
Son, S.W. Jeong, H.K. Yoon, H.S. Ha, M.Y. 

Publisher: 
International Conference on Multiphase Flow (ICMF) 

Publication Date: 
June 4, 2010 
Subjects 

Subject: 
freely falling circular cylinder gap ratio density ratio transverse motion 
Notes 

Abstract: 
A twodimensional circular cylinder freely falling in a channel has been simulated by using Direct forcing/Fictitious domain – lattice Boltzmann method (DF/FDLBM) in order to analyze the characteristics of motion originated by the interaction between the fluid flow and the cylinder. The wide range of the solid/fluid density ratio has been considered to identify the effect of the solid/fluid density ratio on the motion characteristics such as the falling time, the transverse force and the trajectory in the streamwise and transverse directions. In addition, the effect of the gap between the cylinder and the wall on the motion of a twodimensional freely falling circular cylinder has been revealed by taking into account a various range of the gap size. As the cylinder is close to the wall at the initial dropping position, vortex shedding in the wake occurs early since the shear flow formed in the spacing between the cylinder and the wall drives flow instabilities from the initial stage of freely falling. In order to consider the characteristics of transverse motion of the cylinder in the initial stage of freely falling, quantitative information about the cylinder motion variables such as the transverse force, trajectory and settling time has been investigate. 

General Note: 
The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: BioFluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and NanoScale Multiphase Flows; Microgravity in TwoPhase Flow; Multiphase Flows with Heat and Mass Transfer; NonNewtonian Multiphase Flows; ParticleLaden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows 
Record Information 

Bibliographic ID: 
UF00102023 

Volume ID: 
VID00249 

Source Institution: 
University of Florida 

Holding Location: 
University of Florida 

Rights Management: 
All rights reserved by the source institution and holding location. 

Resource Identifier: 
1025SonICMF2010.pdf 

Full Text 
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Numerical study on the motion characteristics of a freely falling twodimensional circular
cylinder in a channel
Sung Wan Son1, Hea Kwon Jeong2, Hyun SikYoon3 and Man Yeong Ha*1
School of mechanical engineering, Pusan Nation University, Korea
2Technical Research Laboratories, Posco, Korea
3Advanced Ship Engineering Research Center, Pusan National University, Pusan, Korea
*corresponding author, Email: myha@pusan.ac.kr
Keywords: freely falling circular cylinder, gap ratio, density ratio, transverse motion
Abstract
A twodimensional circular cylinder freely falling in a channel has been simulated by using Direct forcing/Fictitious domain 
lattice Boltzmann method (DF/FDLBM) in order to analyze the characteristics of motion originated by the interaction
between the fluid flow and the cylinder. The wide range of the solid/fluid density ratio has been considered to identify the
effect of the solid/fluid density ratio on the motion characteristics such as the falling time, the transverse force and the
trajectory in the streamwise and transverse directions. In addition, the effect of the gap between the cylinder and the wall on
the motion of a twodimensional freely falling circular cylinder has been revealed by taking into account a various range of the
gap size. As the cylinder is close to the wall at the initial dropping position, vortex shedding in the wake occurs early since the
shear flow formed in the spacing between the cylinder and the wall drives flow instabilities from the initial stage of freely
falling. In order to consider the characteristics of transverse motion of the cylinder in the initial stage of freely falling,
quantitative information about the cylinder motion variables such as the transverse force, trajectory and settling time has been
investigate.
1. Introduction
The particlefluid interaction has been widely applied in
the fields such as the chemical, civil and aerospace
engineering as well as biological science. Some examples
are the transport of radionuclides, plasma spray coating,
fluidized bed reactor, blood flow and droplet formation.
Moreover, the particlefluidstructure interaction has
influence on the response, stability and life of the structure.
It's important to understand the behavior of particle
suspensions or sedimentations, which has been attracting
lots of researcher's attention for the past decades, both
experimentally [14] and numerically [2,4,517].
For the representative studies about the twodimensional
motion of the cylinder, Hu et al.[5], Feng et al.[6] and
Hu[9] analyzed the sedimentation of the cylinder in a
narrow cavity by solving NavierStokes equations. Hu et
al.[5] showed that a circular cylinder sedimenting in a
narrow channel at small Reynolds number (Re) drifts to the
centre of the channel. At large Re, however, the cylinder
drifts off the centre of the channel and the rotation of the
cylinder is due to uneven shear and oscillation of the
cylinder is due to vortex shedding. Feng et al.[6]
demonstrated the characteristics of the motion of a circular
cylinder sedimenting based on the region of Re dependent
on the cylinder diameter and terminal velocity of the
circular cylinder in a narrow channel. Hu[9] studied the
rotation of a circular cylinder setting close to a solid wall
using the lubrication theory and found that the direction of
rotation reverses at large Re while the rotation is in the
direction as if rolling up the nearer wall at small Re. On
other matters, Namkoong et al.[15] reported the
twodimensional motion of a circular cylinder freely falling
in an infinite fluid for the range of Re<188. They deduced
correlations of the relationship between St (Strouhal
number) and Re from their numerical results.
According to our survey, many researchers numerically
studied the motion of various shape particles such as
circular cylinder, elliptic cylinder and rectangular cylinder
as well as sphere. Based on the results reported by Joseph's
group, however, many researchers including Joseph's
group performed the numerical calculation about the
motion of the cylinder to validate numerical scheme
proposed by them[7,1214,16,17]. Moreover, most
simulations except for some results reported by Feng et
al.[6], which deal with the characteristics of the motion of
the cylinder, mentioned above have been studied for the
motion of the cylinder at nonzero small Re (O
below 200. For a small Re, the previous results by Feng et
al[6] and Hu[9] show the transient trajectory of the cylinder
directs toward the centre of the channel without oscillation.
For a large Re (the density of the cylinder is much larger
than one of the small Re), the trajectory of the cylinder may
be different from a small Re cases. Also, the studies for
characteristics of the motion of the cylinder sedimenting
according to the variation of the density ratio (the ratio
between a cylinder and fluid) and the gap between the
initial position of the cylinder and the channel wall at large
Re are not enough.
In this study the motion of a circular cylinder freely
falling in a narrow channel is studied numerically. The
primary objective of the present study is to investigate the
characteristics of a single cylinder motion freely falling
(including the rotation of the cylinder) in a narrow channel
according to the variation of density ratio and gap ratio at
large Re (Re>200). To evaluate the translational motion of
the cylinder, the transient trajectory, transient transverse
force and settling time of the cylinder have been considered.
The characteristics of the rotation of the cylinder according
to the various gap ratios were also reported. Moreover, the
effect of the rotation of a cylinder on the transverse motion
of one was described.
2. Numerical Method and Validation
2.1 Numerical Method
f/ D
G
h
L t
0 0x
Figure 1: Schematics of system
Figure 1 shows the fluid domain (Q0) represented by
the Eulerian coordinate and the circular cylinder arranged
by the Lagrangian points with M Lagrangian points
uniformly distributed. The Lagrangian points are arranged
as shown in Figure 2.
y
'/z
A AV
I N, 4
I
r12
X
(a) (b)
Figure 2: (a) Definitions for a circular cylinder to arrange
Lagrangian forcing points (r, is the actual cylinder radius).
(b) Arrangement of Lagrangian points in case of a circular
cylinder.
To solve the incompressible NavierStokes equation for
these two different coordinate, the DF/FD LBE was used.
The interaction between the fluid and the cylinder is
calculated by the direct forcing and fictitious domain
method proposed by Yu & Shao[15]. The result of this
interaction is introduced to the governing equation as a
form of external force by the equilibrium velocity approach
Buick & Greated[19].
To investigate the motion of a freely falling circular
cylinder in a closed channel, as shown in figure 1, the
incompressible continuity and momentum equation are
used and given by
V u 0,
+u Vu
+u.Vu
1
Vp+ vV2u +f ,
Pf
where u, P, P v and f are flow velocity, pre
fluid density, kinematic viscosity and external f
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
respectively.
To find the solutions of Eq. (1) and Eq. (2) in the present
study, the LBE is used. Eq. (3) is so called lattice BGK
equation based on the Boltzmann equation [Chen & Doolen
(1998)] and given as
f,(x+c At,t+At) f(x,t)= [f(x,t) f"(x,t)], (3)
r
where f x, c r and f," mean the density
distribution function, position vector, lattice velocity vector,
single relaxation time and equilibrium density distribution
function, respectively. The subscript a is the direction of
particle and depends on the lattice model. In all the
simulations, a twodimensional 9bit model was used. The
equilibrium density distribution function is obtained as
follows:
eq 3ce u 9(ca U)2 3U2
f = co p, 1 + 2 4 (4)
S2c4 2c2
where c = Ax/ At and Ax and At are the
magnitude and time spacing of the particle, respectively.
The weighting coefficient co is coo=4/9 and
co =1/9 for a =1,2,3,4 and co =1/36 for
a = 5, 6,7, 8. Each particle velocity vector is represented
by Eq. (5).
0, a=0
c= (cos[(al)nr/2], sin[(al)nr/2])c, a=1,2,3,4 (5)
c2(cos[(a5)n /2+ n/4], sin[(a5) /2+ / 4])c, a=5,6,7,8.
The macroscopic variables such as the fluid density ( p)
and flow velocity (u) in the computational domain are
defined as
Pf =If, (6)
Pfu = fac (7)
The single relaxation time is related to the kinematic
viscosity by
(Ax)
v = (2r 1) (8)
Thb hncir irlda nf th nFI/Fn mthnrl ic tn vtenrd n
problem on a geometrically complex domain to a larger
(1) simpler domain. To obtain the volume force at a location of
the lattice node, a discrete quantity the Lagrangian force
(F,) at points X, and the Eulerian force (f ) at each
lattice site x, are transferred and represented as follows:
sure,
brce,
F (X,)= f(x,)5,(x,XI), (9)
M
where AV, is a finite volume (surface area for 2D) of
each Lagrangian point and calculated by the same manner
as Uhlmann (2005) and Yu & Shao (2007). 6, is a
smoothed approximation of the discrete Dirac Delta
function reported by Roma Peskin & Berger (1999). For
the case of 2D, the discrete delta function 8h is defined
by
1 xXO y Y
h,(x Xo) =h2 ( X ) h )
h h h
1(53 r 3(1 r)2+1 ,0.5
(r) =I1+
0
3
hrl 0.5
, otherwise
where h is the lattice size.[Yu and Shao (2007)]
When the equilibrium velocity with momentum body
force from immersed boundary method is substituted in the
equilibrium density distribution function, the change in the
momentum at each lattice is obtained. The external force
from Eq. (10) is introduced to Eq. (3) using an equilibrium
velocity approach as follow:
PfU, =pu, + rf (13)
, (x, t)= f (pf,u,) (14)
where u, is the equilibrium velocity at each lattice site
and fe* is the modified equilibrium density distribution
function from the equilibrium velocity approach and
introduced to Eq. (3) instead of fjq.
2.2 Validation
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the free fall of cylinder in the vertical channel given by
Glowinski et al.[13], Wang et al.[15], Wan & Turek[19]
where the starting position of cylinder free fall is located at
the vertical centerline. The specific data for the geometry
and physical properties used in this validation study are
given in Table 1. The cylinder is released at t= and
accelerated by the gravity.
Table 1 Geometry, initial and boundary conditions and
physical parameters
Case 1 2
Computational domain 2 (cm) x (cm) 2 x 6 2 x 6
Diameter of cylinder D, (cm) 0.25 0.25
Density of cylinder p (g/cm3) 1.25 1.5
Density of fluid p (g/cm3) 1.00 1.00
Fluid viscosity v (g/cm3) 0.1 0.01
Initial cylinder location (x, y) (cm) (1,4) (1,4)
Figure 3 shows the comparison of the present
computational results for the time histories of translational
trajectory and translational velocity of the center of a
circular cylinder for Pr =1.5 and v = 0.01 g/cm s (case 2)
with previous results by Glowinski et al.[13] and Wan &
Turek[19]. As shown in this figure, the present
computational results agree very well with computational
results obtained by Glowinski et al.[13] and Wan &
Turek[19]. In this validation test, we considered two
different size of lattice with Ax =1/150 and 1/250 and the
results using Ax =1/150 are almost the same as that using
Ax =1/250. We also calculated the maximum particulate
Reynolds number, Rep which is defined as
Rep,m =max[LpD Uc2 +VT2 /v] where p D,, Uc,
V, and v represent the cylinder density, cylinder
diameter, transverse velocity, translational velocity and
fluid viscosity, respectively. Table 1 shows the comparison
of the present computational results for the maximum
particulate Reynolds number (Repm.) during a cylinder
sedimentation with previous results by Glowinski et al.[13],
Wang et al.[15] and Wan & Turek[19].
To validate the computer code and methodology
developed in the present study, we considered the case of
Table 2 Comparisons of the maximum particulate Reynolds number (Repm ) during free fall of cylinder
Present Glowinski et al.[13] Wang et al.[15] Wan and Turek [19]
Ax 1/50 1/100 1/192 1/256 1/72 1/144 1/256 1/48 1/96
Case 1
Rep,ma 16.87 17.24 17.27 17.31 16.962 17.216 17.307 17.42 17.15
S Ax 1/150 1/250 1/256 1/384 1/72 1/144 1/256 1/48 1/96
Case 2
Rep,,ma 489.1 478.69 450.7 466 502.37 503.26 503.38 442.19 465.52
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
UI.
. 100
Present, Ax=l/150
......... Present, Ax=l/250
10 o Wan and Turek (2007)
S Glowinski et at (2001)
0
5
10
150.1 0.2 0.3 0.4
0 0.1 0.2 0.3 0.4
Figure 3: Time histories of (a) translational trajectory and
(b) translational velocity of the center of a circular cylinder
for pr =1.5 and = 0.01g/(cms).
As shown in this Table 2, our present results for Rem
represents well the computational results given by
Glowinski et al.[13], Wang et al.[15] and Wan & Turek[19]
even though the size of lattice used in the present study is
smaller than that used by previous studies.
In order to validate the present computational results for
the rotation of cylinder due to the interaction between the
fluid and the cylinder, we calculated the rotation rate ( )
of the suspended cylinder in simple shear flow and
compared our computational results with the experimental
results given by Zettner & Yoda[20] and computational
results by Ding & Aidun[21]. Here the height of channel
considered in the present comparison is 2 D and the
width of channel is 4 D where D is the diameter of a
circular cylinder (D ) as 0.25cm. The circular cylinder is
located at the center of channel and the upper and lower
walls of channel moves in the opposite direction in order to
make a shear flow. The viscosity and the density ratio
between the fluid and cylinder used in this comparison are
0.05 g/cm s and 1, respectively. Figure 4 shows the
comparison of the present computational results for Oc /
of the circular
101
Re
Figure 4: Comparison of the simulation results on the
~c / X of rotation of a circular cylinder in shear flow with
present and previous results.
cylinder rotation in the shear flow with previous
experimental and computational results by Zettner & Yoda
[20] and Ding & Aidun[21]. Here the Reynolds number
(Rex ) used in this validation is defined as
Rez = D /v, (15)
where y = du / dy presents the shear rate. Rex s
considered in this validation are 1, 10, 76.8 and 100,
respectively. The computational domain used in this
validation test is 145 x 73 lattice units for all Re .
When a cylinder suspends freely in between two parallel
plates moving in opposite directions, our computational
results are in good agreement with previous experimental
results by Zettner & Yoda[20] and numerical results by
Ding & Aidun[21].
3. Results and Discussion
A Study for motion characteristics of a freely falling
circular cylinder in the long channel was performed. For
simulation of free fall of a single circular cylinder for 2D,
the specific information for the geometry and physical
properties used in this study are given in Table 3. The
gravitational acceleration is g=981 cm/s2. The cylinder is
released at t= and accelerated by the gravity. For an
evaluation of the transient characteristics of cylinder by the
initial location of cylinder, the gap ratio is defined as
G, = G/D,. The gap ratio is varied from 0.1 to 3.5 for all
density ratios considered in the present study. G =0.1 is
the case that the cylinder is closest to the wall and G, =3.5
is the case that the cylinder is located at the centre of the
channel. Both the cylinder and the fluid are at rest at time
t=0 and the cylinder starts abruptly freefall motion due to
gravity in channel for t > 0.
0 Zetner & Yoda (2001, Ep.)
+ Ding & Aidun (2000, Num.)
O Present
+I
,. I .. I ,
nr I
Table 3 Geometry, initial/boundary conditions and
physical/numerical parameters
values
Computational domain 2 (cm) (cm)
Diameter of cylinder D, (cm)
Density of cylinder Pp (g/cm3)
Density of fluid Pf (g/cm3)
Fluid viscosity v (g/cm3)
Initial cylinder location (x, y) (cm)
Density ratio P, = P, P
Gapratio G =G /D
Lattice size Ax (cm)
Time step At (sec)
2 15
0.25
1.25, 1.5, 1.75
1.00
0.01
(1, 14)
1.25, 1.5, 1.75
0.1, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5
1/150
5 105
3.1 Flow field
In this section, we present some typical results on the
vortical structures and the cylinder behaviors for a freely
falling cylinder in the narrow vertical channel.
Figure 5 shows the time evolution of vorticity field around
a cylinder for different G, values of 0.1, 1.0, 2.5 and 3.5
at p =1.5 while the cylinder falls freely in the narrow
vertical channel. When G, =3.5, the starting position of the
cylinder is at the horizontal centre of channel (x =0.0) and
as a result the transverse force acting on the cylinder by the
fluid in the horizontal direction is balanced. As a cylinder
falls freely, a pair of symmetric vortex is formed around the
cylinder and grows in its size by maintaining its symmetric
shape, so that we cannot observe any vortex shedding in the
wake of cylinder as shown in Figure 5(d).
When G =2.5, the starting position of cylinder moves
slightly to the right wall of channel and is located at x =1.0.
Because the distance from the right and left walls to the
cylinder center is not the same at G =2.5, the transverse
force acting on the cylinder caused by the fluid is not
balanced. As a result, a pair of symmetric vortices formed
around the cylinder at the initial stage of free fall of
cylinder does not maintain its symmetric shape as the
cylinder falls freely and the counterclockwise vortex
formed in the left side of near wake around the cylinder
start to be inclined to right side at t =0.25. As the cylinder
falls further freely, we can observe the generation of the
periodic Karman vortex shedding in the wake of cylinder as
time goes by. We can also observe the formation of weak
shear layer on the right wall of channel and its development
as a function of time, because of the interaction between
the right wall of channel and the vortex adjacent to the right
wall in the wake region of the cylinder.
As the distance between the initial position of cylinder and
the right wall of channel decreases with decreasing G,
values from 3.5 to 0.1, the transverse force imbalance
acting on the cylinder caused by the fluid increases. As a
result, as G, decreases, the vortex shedding in the wake
occurs in the earlier time, the extent of vortex shedding
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
formed in the wake region increases, and the interaction
between the vortex in the wake and the right wall of
channel becomes stronger. We can also observe the
transverse motion of cylinder in addition to the free falling
motion of cylinder more clearly with decreasing G .
Especially, when G,=0.1 in which the distance between
the cylinder and the right wall of channel at the starting
point of cylinder free fall is very narrow, the interaction
between the fluid flow around the falling cylinder and the
right wall of channel becomes very strong.
3.2 Force in the transverse directions
Figures 6(a) and 6(b) show the time history of transverse
force and the time history of transverse trajectory, velocity
and force acting on the cylinder by the fluid, respectively,
for different density ratios of p,.=1.25, 1.5 and 1.75 at
G =0.1. After the starting transient time period shown in
figures 6(a) and 6(b), the regularly periodic transverse
force acts on the falling cylinder. When the cylinder density
becomes larger than the fluid density with increasing
density ratio, the amplitude and frequency of oscillating
transverse force acting on the cylinder decrease.
Figure 6(c) shows the distribution of instantaneous
vorticity contours at five different instants denoted by A, B,
C, D and E in figure 6(b), respectively, at the initial stage
just after the cylinder falls freely, when p, =1.5 and
G =0.1. In figure 6(c), the positive and negative values of
instantaneous vorticities are denoted by the solid and dotted
lines, respectively, in the contour range from 300 to 300
with 14 levels. When t =0.04 corresponding to A in figure
6(b) just after the cylinder falls freely in the channel, the
negative strong repulsive transverse force from the right
wall of channel acts on the cylinder because the gap
between the cylinder and the right wall of channel is very
narrow at the value of G =0.1 and as a result the cylinder
migrates to the left direction. At this instant, the negative
and positive vorticities start to be developed from the left
side of cylinder and on the right wall of channel,
respectively, as shown in figure 6(c). When t =0.07
corresponding to B in figure 6(b), the negative repulsive
transverse force between the cylinder and the right wall of
channel reaches a peak value. Because the gap between the
cylinder and the right wall of cylinder at t =0.07 becomes
larger than that at t =0.04 due to the continuous movement
of cylinder to the left direction in the presence of negative
repulsive transverse force, the positive vorticity starts to be
developed from the right side of cylinder whereas the
negative vorticity in the left side of cylinder develops
further and rotates slightly in the clockwise direction. The
positive vorticity on the right wall of channel is elongated.
As the cylinder keeps moving to the left direction due to its
inertia, the magnitude of negative transverse force
decreases and reaches a zero value at t =0.1 corresponding
to C in figure 6(b). At this instant of time, the negative
vorticity in the left side of cylinder rotates further to the
clockwise direction and becomes very close to the positive
vorticity rotating to the counterclockwise direction from
the right side of cylinder.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
10.200
10.200
6'
0.200
i
10.200
0.250
0.300
(a) Gr =0.1
t=0.250
1=0.300
(b) G( =1.0
t0.250
t0.300
m
(c) G = 2.5
t=0.250
1=0.300
(d) G =3.5
10.350
0.350
0.350
0.350
t0.400
t0.400
10.400
t0.400
A
10.450
i
04
"4:
t10.450
J
10.450
t10.450
6
10.500
OD
t0.500
*S
t0.500
t=0.500
I
Figure 5: Evolution process of vorticity field and a cylinder position during sedimentation for p, =1.5.
(vorticity range: 50 50)
The negative vorticity starts to be developed on the right
wall of channel due to the interaction with the positive
vorticity formed in the right side of cylinder whereas the
positive vorticity on the right wall of channel keeps being
elongated in the upward direction. During the time interval
between t =0.1 and t =0.175, the transverse force acting
on the cylinder becomes positive, keeps increasing and
reaches a peak value at t =0.14 corresponding to D in
figure 6(b). The cylinder keeps moving to the left direction
during the time interval between t =0 and t =0.14. The
distance between the cylinder and the right wall of channel
at t =0.14 has a maximum value during the starting time
period of t =00.2 seconds. At t =0.14, the negative
vorticity rotating clockwisely in the left side of cylinder
starts to be torn by the positive vorticity rotating
counterclockwisely from the right side of cylinder. During
the time intervalbetween t =0.14 and t =0.2, the
transverse force acting on the cylinder decreases with a
positive value and the transverse motion of cylinder
changes its direction from left to the right direction in the
presence of positive transverse force.
10.050
4
t0.050
o
t0.050
0
t0.050
*
t0.150
t0.150
b
t0.150
6
t0.150
10.100
I
t=0.100
I
10.100
6
t=0.100
6
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
20
 p,=1.25
S p= 1. 5
10  p 1.75
A C BE
(c)
Figure 6: (a) Time histories of transverse force for various
density ratios at G,=0.1, (b) Time history of transverse
force, trajectory and velocity for p,=1.5 at G,=0.1 at
initial stage of freely falling and (c) instantaneous vorticity
contours (Contour values range from 300 to 300 with 14
levels; Positive solid, Negative dashed).
At t =0.175, the negative vorticity rotating clockwisely
from the left side of cylinder is almost torn by the positive
vorticity rotating clockwisely from the right side of cylinder.
So, when G =0.1, we can observe the strong interaction
between the fluid flow around the falling cylinder and the
right wall of channel, giving the strong influence on the
motion of falling cylinder, as shown in figure 6(c).
As G, increases from 0.1 to 1.0, the gap between the
falling cylinder and the right wall of channel increases and
as a result the influence of the presence of channel right
wall on the transverse force acting on the cylinder
decreases. Figures 7(a) and 7(b) show the time history of
transverse force and the time history of transverse trajectory,
velocity and force acting on the cylinder by the fluid,
respectively, for different density ratios of p,=1.25, 1.5
and 1.75 at G =1.0. The variation of transverse force
acting on the cylinder as a function of time at G =1.0 is
generally similar to that at G =0.1. As a result the
regularly oscillating period of transverse force is followed
by the starting transient oscillating period and the
amplitude and frequency of oscillating transverse force
decreases as the density ratio increases. However, when
G, =1.0, the starting transient time period of oscillating
A,
B
I
C
8
Figure 7: (a) Time histories of transverse force for various
density ratios at G,=1.0, (b) Time history of transverse
force, trajectory and velocity for p,=1.5 at G,=1.0 at
initial stage of freely falling and (c) instantaneous vorticity
contours (Contour values range from 300 to 300 with 14
levels; Positive solid, Negative dashed).
transverse force is longer than that when G, =0.1.
Figure 7(c) shows the distribution of instantaneous
vorticity contours at five different instants denoted by A, B,
C, D and E in figure 7(b), respectively, at the initial stage
just after the cylinder falls freely, when p, =1.5 and
G =1.0. When t =0.04 corresponding to A in figure 7(b)
just after the cylinder falls freely in the channel, the
transverse force acting on the cylinder by the fluid at
G=1.0 is still almost zero unlike to the case of G,=0.1
with the negative transverse force and as a result the
distribution of vorticity formed in the wake of cylinder at
G, =1.0 is almost symmetric. When t =0.12 corresponding
to B in figure 7(b), the small negative peak transverse force
at G =1.0 starts to act on the cylinder and as a result the
symmetric shape of vorticity formed in the wake of
cylinder starts to be broken. Unlike to the case of G =0.1,
when G, =1.0, the negative vorticity is formed on the right
wall of channel and the size of positive vorticity in the right
side of cylinder is larger than the size of negative vorticity
Fx/110
 ........ U'/10
'..   x
D E
0.2 0.3.. U/10.
. .. ...3 0.4 0.
0.1 0.2 0.3 0.4 0.5
 F ,/
D
A C ....
0 0.
0 0.1 0.2 0.3 0.4 0.5
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
in the left side of cylinder. The positive vorticity in the right
side of channel interacts with the negative vorticity formed
on the right wall of channel and is elongated in the
streamwise direction. At t =0.15 corresponding to C in
figure 7(b), the positive transverse force starts to act on the
cylinder while the cylinder still keeps moving to the left
direction. The negative vorticity in the left side of cylinder
rotates in the clockwise direction whereas both the positive
vorticity in the right side of cylinder and the negative
vorticity on the right wall of channel are elongated in the
streamwise direction. At t =0.19 corresponding to D in
figure 7(b), the positive transverse force acting on the
cylinder has a positive peak value and the cylinder changes
its direction to the right. The positive vorticity in the right
side of cylinder starts to be torn by the negative vorticity
rotating clockwisely from the left side of cylinder. The
negative vorticity on the right wall of channel is elongated
further due to the interaction with the positive vortex in
right side of cylinder. At t =0.22 corresponding to E in
figure 7(b), the transverse force acting on the cylinder
becomes almost zero and the distance between the cylinder
and the right wall of channel is smallest during the time
period of t =00.25. The negative vorticity rotating
clockwisely from the left side of cylinder starts to be torn
by the positive vorticity rotating counterclockwisely from
the right side of cylinder. The positive vorticity separated at
t =0.19 is convected to the streamwise direction. The
negative and positive vorticies are formed in series on the
right wall of channel which matches with corresponding
positive and negative vorticies around the cylinder.
When G, is 2.5, the influence of the presence of
channel right wall on the transverse force acting on the
cylinder is much smaller than that when G, is 0.1 and 1.0,
because the starting position of cylinder free fall is close to
the channel center line and the distance between the falling
cylinder and the right wall of channel at G, = 2.5 becomes
larger than that at G,= 0.1 and 1.0. Figures 8(a) and 8(b)
show the time history of transverse force and the time
history of transverse trajectory, velocity and force acting on
the cylinder by the fluid, respectively, for different density
ratios of p, =1.25, 1.5 and 1.75 at G =2.5. Similar to the
cases of G, =0.1 and 1.0, the transverse force acting on the
cylinder at G =2.5 oscillates regularly as a function of
time after the starting transient oscillating period and its
amplitude and frequency increase with increasing density
ratio. However, when G =2.5, the starting transient time
period of oscillating transverse force is much longer than
that when G, =0.1 and 1.0.
Figure 8(c) shows the distribution of instantaneous
vorticity contours at five different instants denoted by A, B,
C, D and E in figure 8(b), respectively, at the initial stage
just after the cylinder falls freely, when p, =1.5 and
G =2.5. At t =0.06 and 0.13 corresponding to A and B in
figure 8(b) during the starting transient period of time after
the cylinder falls freely in the channel, a pair of symmetric
vorticities in the wake of cylinder and the negative
vorticity on the right wall of channel are formed, because
A
I
S F, /10
..     /
0 A B C D
0 0.1 0.2 0.3 0.4 (
B
i6
C
SD
Figure 8: (a) Time histories of transverse force for various
density ratios at G,=2.5, (b) Time history of transverse
force, trajectory and velocity for p,=1.5 at G,=2.5 at
initial stage of freely falling and (c) instantaneous vorticity
contours (Contour values range from 300 to 300 with 14
levels; Positive solid, Negative dashed; Dashed dot:
horizontal centre line of channel).
the transverse force acting on the cylinder by the fluid at
G =2.5 is almost zero for the longer period of time than
that at G =0.1 and 1.0. At t =0.18 and 0.215
corresponding to C and D in figure 8(b), because the small
positive transverse force starts to act on the cylinder, the
symmetric shape of a pair of vortices starts to be broken.
As a result the positive vorticity in the right side of cylinder
is slightly longer in the streamwise direction than the
negative vorticity in the left side of cylinder due to the
interaction between the positive vorticity in the right of
cylinder and the right wall of channel. During this period of
time, a pair of vortices in the wake of cylinder and the
negative vorticity on the right wall of channel is elongated
in the streamwise direction as time goes on. At t =0.26
corresponding to E in figure 8(b)), a pair of vortices are
elongated in the streamwise direction and starts to oscillate
in the presence of negative transverse force.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
3.3 Trajectories in the transverse directions
1.8
IP
1.6
1.4
1.2
1
(b)
Figure 9: (a) Transverse trajectories of a circular cylinder
at various gap ratio and density ratio on xy plane, and (b)
time evolution of transverse trajectories of a circular
cylinder at various gap ratios and density ratios at initial
stage of freely falling.
Figures 9(a) and 9(b) show the trajectory of a cylinder in
the x y plane and the initial evolution of transverse
trajectory of a cylinder, respectively, as a function of time
for different gap ratios and density ratios. The interaction
between the presence of channel wall and vortex shedding
formed in the wake of cylinder determines the trajectory of
cylinder which falls freely in the channel. At the initial
stage of free fall of cylinder, the effect of channel wall
presence on the cylinder trajectory is larger than the effect
of vortex shedding in the wake. As a result, as the gap ratio
decreases, the transverse repulsive force increases and the
distance which the cylinder moves to the left direction
from the starting inlet position increases as shown in
figures 9(a) and 9(b). During this starting period of time, a
pair of vortex formed in the wake of cylinder is almost
symmetric and as a result the trajectory of cylinder does
not depend on the density ratio. In the following transient
period of time after the initial stage of free fall of cylinder,
a symmetric shape of vortices formed in the wake of
cylinder is broken and starts to oscillate due to the vortex
shedding. As a result the trajectory of cylinder oscillates
due to this vortex shedding with increasing amplitude and
frequency as the density ratio increases, in addition to the
cylinder movement to the left caused by the transverse
Figure 10: Settling time of a circular cylinder at various
gap ratios and density ratios.
repulsive force in the presence of channel wall.
As the gap ratio increases, the cylinder moves the longer
distances from the starting point of cylinder free fall before
it starts the oscillatory motion. Especially, when G = 3.5,
the cylinder falls a quite long distance vertically from the
starting position ( y=56 ) to y 32 without any
transverse motion in the transverse direction as shown in
Figures. 9(a) and 9(b), unlike to different cases of G, =0.1,
1.0 and 2.5 with some transverse motion to the left.
However, the distance for the cylinder to fall freely without
any oscillatory motion in the transverse direction does not
depend on the density ratio, meaning that the trajectory of
cylinder starts to oscillates after the cylinder moves the
same distance from the starting position of cylinder free
fall for different density ratios at the same gap ratio. After
the transient periodic trajectory of cylinder, the
quasisteady periodic trajectory of cylinder is followed
with different amplitudes and frequencies depending on the
density ratio until the cylinder arrives at the bottom of
channel. When G =0.1, 1.0 and 2.5, the cylinder hits the
bottom of channel at the similar position because the
cylinder movement to the left increases due to increasing
repulsive force in the presence of channel wall as the gap
ratio decreases. However, when G = 3.5, the extent in
which the cylinder trajectory deviates from the vertical
centerline is very small because the effect of the channel
wall presence on the cylinder trajectory is negligible.
3.4 The effects of rotation of the cylinder
Figure 10 shows the settling time (t,) of free fall cylinder
as a function of gap ratio for different density ratios. Here
the settling time represents the time required for a free fall
cylinder to arrive at the bottom of the channel. As the gap
ratio increases for the specified density ratio, the distance,
which the trajectory of cylinder travels, decreases with
decreasing transverse motion of cylinder before the
cylinder hits the bottom of channel and as a result the
settling time decreases linearly except the case of G, =3.5.
When G =3.5, the settling time does not follow the linear
p,=1.25
A p,=.5
p,=1.75
. U .
*
A A A A A A A
0. '' 0.5 1 .5 2 2.5 3 .5.
0 0.5 1 1.5 2 2.5 3 3.5
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Figure 11: Time histories of rotational velocity(co ) and
angle of rotation ( 0) of a circular cylinder for various gap
ratios at p, =1.5.( G, =0.1(dashdot), 1.0(dashed),
2.5(solid).
variation as a function of time and is less than the value
when it follows the linear variation. As the density ratio
increases for the specified gap ratio, the settling time
decreases because the weight of cylinder with higher
density ratio is larger than that with lower density ratio.
When the cylinder falls freely, the cylinder has the
rotational motion in addition to the linear motion. Figure 11
shows the time histories of the angular velocity (co) and
the angular position of the cylinder relative to its initial
angle (O ) when the cylinder starts to fall freely for
different values of G,=0.1, 1.0 and 2.5 at p, =1.5. When
G =0.1, the shear force difference acting on the left and
right sides of cylinder is very large due to a small gap
between the cylinder and the right wall of channel, and as a
result the cylinder starts to rotate in the clockwise direction
with a large value of angular velocity as soon as it starts to
fall freely from its starting position. The rotational direction
obtained from the present computation under the narrow
gap condition between the cylinder and the right wall of
channel is qualitatively similar to that obtained from the
steady state simulation results based on the lubrication
theory by Hu (1995). At the starting time when the cylinder
starts to fall freely, the value of the clockwise angular
velocity increases abruptly and has a maximum value at
t 0.1. In the following time, the angular velocity shows
the transient periodic shape with decreasing absolute
magnitude until it reaches the steady periodic state, which
matches well the history of cylinder trajectory as shown in
figure 9. If the cylinder trajectory and the angular velocity
reach the steady periodic state, the angular velocity
oscillates regularly around a value of zero with a small
amplitude, meaning that the cylinder rotation becomes
much smaller than the linear motion of cylinder at the
steady periodic state. While the cylinder falls freely at
G =0.1 by following the trajectory shown in figure 9, the
angular position of the cylinder relative to the initial angle
( O) keeps rotating in the clockwise direction, with a rapid
decrease in its magnitude at the initial stage of free fall of
cylinder which is followed by the gradual decrease in its
Figure 12: Comparison of time histories of transverse
trajectory and velocity from 'thought experiments'.
magnitude. When the cylinder hits the bottom of channel at
Gr=0.1, 6 134.70.
When the gap ratio increases from 0.1 to 1.0 and 2.5, the
gap between the cylinder and the right wall of channel
increases and as a result the shear force difference acting
on the left and right sides of cylinder at G =1.0 and 2.5 is
much smaller than that at G =0.1. Thus the angular
velocity at G =1.0 and 2.5 does not decrease abruptly at
the initial stage of cylinder free fall unlike to the case of
G =0.1. The angular velocity at G =1.0 and 2.5 oscillates
regularly around a value of zero in both transient and
steady periodic states. The angular position of the cylinder
relative to the initial angle at G,=1.0 and 2.5 also keeps
rotating in the clockwise direction according to the history
of cylinder trajectory, but the absolute magnitude of O6 at
G =1.0 and 2.5 is much smaller than G =0.1. When the
cylinder hits the bottom of channel at G,=1.0 and 2.5,
q z 24.2 and 4.25 respectively. Thus, when
G, =1.0 and 2.5, the effect of cylinder rotation is relatively
small, compared to case of G =0.1.
In order to consider the effect of cylinder rotation on the
trajectory of cylinder and the transverse force acting on the
cylinder, we calculated the transverse movement of
cylinder ( xe ) and the transverse velocity (Uc ) of cylinder
by considering both cases with and without the cylinder
rotation. Figure 12 shows the time history of the transverse
trajectory and the transverse velocity (Uc) of cylinder for
both cases with and without the cylinder rotation for
p,=1.5 and G,=0.1.
4. Conclusions
A freely falling circular cylinder in a long channel has
been simulated using the DF/FDLBM method for 2D,
which combines the desired features of the
Directforcing/Fictitious domain method and lattice
Boltzmann method.
The gap ratio between the wall and the cylinder is varied
from 0.1 to 3.5 for three density ratios (p, =1.25, 1.5 and
1.75) considered in the present study.
For all p, s, when G, is decreased, the interaction
between the fluid flow around the falling cylinder and the
right wall of channel becomes very strong. Therefore the
vortex shedding in the wake region of the cylinder is
formed earlier and influence on the transverse trajectory as
well as translational trajectory of the cylinder, because the
transverse force imbalance acting on the cylinder caused by
the fluid increases.
For all G, s, when p, is increased, the transverse force
acting on the cylinder by the fluid is larger than small p,.
However the period of the oscillation with transverse
direction is decreased and the amplitude is increased.
Acknowledgements
This work was supported by the Korea Foundation for
International Cooperation of Science &
Technology(KICOS) through a grant provided by the Korea
Ministry of Education, Science & Technology(MEST) in
2009 (No. K2070200001307E020001310).
References
1. Fortes, A. F., Joshep, D. D. & Lundgren, T.S., Nonlinear
mechanics of fluidization of beds of spherical particles, J.
Fluid Mech., 177, 467483 (1987)
2. Cate, A. ten, Nieuwstad,C. H., Derksen, J. J. & Van den
Akker, H. E. A., Particle imaging velocimetry experiments
and latticeBoltzmann simulations on a single sphere
settling under gravity, Phys. Fluids, 14, 11, pp.40124025
(2002)
3. Pan, T W, Joseph, D. D., Bai, R., Glowinski, R. &
Sarin, V, Fluidization of 1204 spheres: simulation and
experiment, J. Fluid Mech., 451, 169191 (2002)
4. Jenny, M., Dusek J. & Bouchet, G., Instability and
transition of a sphere falling or ascending freely in a
Newtonian fluid, J. FluidMech., 508, 201239 (2111 1).
5. Hu, H. H., Joseph, D. D. & Crochet, M. J., Direct
simulation of fluid particle motions, Theoret. Comput. Fluid
Dynamics, 3, 285306 (1992)
6. Feng, J., Hu, H. H., & Joseph, D. D., Direct simulation of
initial value problems for the motion of solid bodies in a
Newtonian fluid Part 1. Sedimentation, J. Fluid Mech., 261,
95134 (1994)
7. Ladd, A. J. C., Numerical simulations of particulate
suspensions via a discretized Boltzmann equation. Part 1.
Theoretical foundation, J. FluidMech., 271, 285309 (1994)
8. Ladd, A. J. C., Numerical simulations of particulate
suspensions via a discretized Boltzmann equation. Part 2.
Numerical results, J. FluidMech., 271, 311339 (1994)
9. Hu, H. H., Motion of a circular cylinder in a viscous
liquid between parallel pates, Theoret. Comput. Fluid
Dynamics, 7, 441 (1995)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
10. Hu, H. H., Direct simulation of flows of solidliquid
mixtures, It. J. Multiphase Flow, 22, 2, 335352 (1996)
11. Qi, D., LatticeBoltzmann simulations of particles in
nonzero Reynoldsnumber flows, J. Fluid Mech., 385,
4162 (1999)
12. Feng, Z. G. & Michaelides, E. E., The immersed
boundarylattice Boltzmann method for solving
fluidparticles interaction problems, J. Comput. Phys., 195,
602628 (2" 114)
13. Glowinski, R., Pan, T. W., Hesla, T. I., Joseph, D. D. &
Periaux, J., A fictitious domain approach to the direct
numerical simulation of incompressible viscous flow past
moving rigid bodies: Application to particulate flow, J.
Comput. Phys., 169, 363426 (2001)
14. Yu, Z., & Shao, Z., A directforcing fictitious domain
method for particulate flows, J. Comput. Phys., 227,
292314 (2007)
15. Namkoong, K., Yoo J. Y. & Choi, H.G., Numerical
analysis of twodimensional motion of a freely falling
circular cylinder in an infinite fluid, J. Fluid Mech., 604,
3353 (2008)
16. Wang, Z., Fan, J., & Luo, K., Combined multidirect
forcing and immersed boundary method for simulating
flows with moving particles, It. J. Multiphase Flow, 34,
283302 (2008)
17. Uhlmann, M., An immersed boundary method with
direct forcing for the simulation of particulate flows, J.
Comput. Phys., 209, 448476 (2005)
18. Buick, J. M. & Greated, C. A., Graivity in a lattice
Boltzmann model, Phys Rev. E, 61, 5307 (2000)
19. Wan, D., & Turek, S., An efficient multigridFEM
method for the simulation of solidliquid two phase flows, J.
Comp. Appl. Math., 203, 561580 ( 2007)
20. Zettner, C. M. and Yoda, M. The circular cylinder in
simple shear at moderate Reynolds numbers: An
experimental study, Experiments in Fluids, 30, 346353
(2001)
21. Ding, E. and Aidun, C. K., The dynamics and scaling
law for particles suspended in shear flow with inertia, J.
FluidMech., 423, 317344 (2000)

