7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
An Immersed Boundary Method for Interacting Particles
Marian Zastawny* Berend van Wachem* and Jos6 Oliveira *
Department of Mechanical Engineering, Imperial College London, Exhibition Road, London, SW7 2AZ
t Sygma Motors, Sao Jose dos Campos, SP, Brazil
B.vanWachem@imperial.ac.uk
Keywords: Immersed Boundary Method, Direct Numerical Simulation, Particulate Flows
This work describes a novel implicit mirroring immersed boundary method (MIB) for simulating the detailed flow
around arbitrary moving bodies. An immersed boundary method uses a nonboundary conforming grid in the whole
domain, to enable efficient solution of the flow equations, and a Lagrangian representation of the boundary of an
arbitrary particle (the immersed boundary) to determine the behaviour of the particle. The boundary of the particle is
triangulated and at the intersections of the triangles with the grid, socalled immersed boundary conditions are applied.
The implicit immersed boundary condition used in this work, mirrors the velocity field along the normal of the local
triangulated immersed boundary (IB) segments ensuring that the fluid exactly follows the IB. As a result, a fictitious
velocity field inside the immersed boundary is developed, which is excluded from the discretised NavierStokes
equations to preserve mass. The numerical method is based on a finite volume approach on a collocated variable
arrangement with a fully coupled pressurevelocity solver. In this paper, the method is employed on the uniform flow
past a rotating sphere at various Re numbers. Because of the rotation, next to the drag force a lift force is developed
as well. Both drag force and lift force are compared to theoretical and experimental findings. The method is also used
to simulate a sedimenting sphere and is compared to literature data. Both cases show that the force on the particles is
predicted with second order accuracy.
Introduction
Multiphase flows are a common phenomena. They
exist widely in nature and can be observed in such pro
cesses as formation of clouds, erosion, sediment trans
portation in rivers etc. Multiphase flows are also dom
inant in industrial processes. Particle transport, drying,
combustion, flows through filters are just a few exam
ples. So far, however, only a limited understanding has
been achieved by the researchers working in this field.
Despite the rapid development in the area of the com
putational fluid dynamics (CFD) the fundamental prin
ciples describing multiphase flows are still not fully de
scribed.
The main challenge in understanding the multiphase
flows lies in their complexity. Presence of multiple par
ticles immersed in a fluid implies that to predict the be
haviour of the flow, it is not enough to find the forces ex
erted by the fluid on particles, but also particleparticle
interactions have to be considered. Additional chal
lenges arise when turbulence is involved, making exact
evaluation of the flow so expensive (in a computational
sense), that it becomes unreasonable to perform it.
Resolving highly turbulent flows requires fine mesh
ing to compute all the flow details along with the effects
of interaction between the fluid and particles. Direct
numerical simulations (DNS), the approach attempting
to compute all the features of the flow is therefore lim
ited only to simple cases with small number of bodies
present in the fluid and small turbulence intensity. For
more practical cases different methods, such as Large
Eddy Simulations (LES) or Reynolds Averaged Navier
Stokes equations (RANS) are employed. Instead of per
forming detailed calculations at the smallest scale, they
solve only the part of the problem approximating the fine
details through the turbulence models. Still, in order to
properly simulate the complex coupling between parti
cles and fluid additional terms have to be introduced in
the NavierStokes equations. Due to its ability to capture
fine details of the flow, DNS can be employed to study
the fundamental phenomena and to derive models that
can be applied in more complicated cases.
Throughout the years, various methods coupling the
particles with surrounding fluid have been developed.
Among the oldest ones is the arbitraryLagrangian
Eulerian method (4), where an unstructured grid is cre
ated around the bodies and the mesh is adapted as they
move. Although this method works very well when the
deformations are small, like in aerodynamic problems,
the need to create a new mesh every time step if the de
formations are large, is time consuming and may limit
accuracy.
The immersed boundary (IB) method tackles this
problem by using Cartesian grid for the fluid, while the
presence of particles is accounted by modifying flow
variables where the IB crosses the Eulerian mesh cells.
Peskin (10) was among the first who used the IB
method. In his implementation the coupling between
the phases is achieved by calculating point forces rep
resenting the influence of the boundaries. The forces are
distributed over the mesh using a distribution function.
This method is first order accurate is space and time. It
also leaves a blurred representation of the boundary.
Another way to implement the IB method, sometimes
called immersed interface method, is to apply the bound
ary conditions directly at the interface of the particle. It
is introduced by Mohd & Yusof in (8), where a smooth
velocity gradient is applied over the IB. Despite being
able to achieve some promising results, the method ex
hibits problems with mass conservation in the boundary
cells.
This paper presents an implicit mirroring immersed
boundary method for moving bodies. The first ideas on
this method are introduced in (7), although a number of
improvements have been made. The particle surface is
triangulated and the intersections of the surface with the
fluid grid are found. The fluid velocities are mirrored
through the surface to guarantee that no slip condition is
satisfied at the surface. The fictitious field generated in
this way is excluded from the pressure equation to pre
serve mass. The method has shown its ability to predict
drag on a sphere with second order accuracy. It is imple
mented in a fully coupled flow solver, to allow one outer
iteration per timestep.
The capability of the method to deliver second order
accurate results of the forces acting on a particle will be
researched in this paper based on two test cases. The
flow around rotating sphere in the low Reynolds num
ber regime up to Re 0.5 along with simulations at
Re 100 are computed and compared with theoretical,
experimental and computational results. Subsequently,
the final velocity of a spherical particle sedimenting in a
quiescent fluid is calculated and compared with experi
mental data.
Implicit Immersed Boundary Method
The mirroring IB (MIB) method is introduced and val
idated by Mark and van Wachem in (7). Due to the fact
that it uses a nonboundary conforming grid it is capable
to resolve detailed flow around an arbitrary body. The
surface of the particle is triangulated and an immersed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
boundary condition is applied at the cells where the body
surface triangles intersect the fluid Eulerian grid. In the
MIB method, the velocity field around the body is mir
rored through each triangle adding a constraint to the
NavierStokes equations and creating a fictitious flow
inside the particle. In order to preserve mass, the fic
titious flow is discarded in the pressure correction equa
tion. The method has been proven to be second order
accurate in predicting the drag force on a particle. The
Figure 1: Mirroring Immersed Boundary Method im
plementation.
general algorithm of the method can be seen in Fig. 1.
At the beginning of the simulation, the Eulerian grid
is loaded along with the triangulated immersed bound
aries. Then the grid cells are tagged according to their
position relatively to the IB triangles as presented in Fig.
2. Three types of cells are distinguished: cells inside the
body (IB), cells next to the at least one outside cell (IIB)
and cells outside the IB (tagged with 0). Finally, for
each cell inside and near the IB (IIB) the centre is mir
rored through the closest surface triangle to create the r'
point.
The coefficients resulting from the linearised Navier
Stokes equations are modified by the IB condition. For
the external cells (tag 0), the coefficients are determined
from central discretisation of the momentum equations.
Since the velocity at points inside the IB has to be equal
to the velocity of the body, coefficients at the cells far in
Figure 2: The markings of the cell centred u veloci
ties. The exterior points, are marked by 0, the IIB ve
locity points, are marked by 1 while the interior velocity
points, are marked by 2. r' is the cell centre mirrored
through the IB.
side the IB are set accordingly. In case of the IIB points
the procedure is slightly more complicated. Firstly the
velocity at the 'I point is determined by trilinear inter
polation from the neighboring external cells. Then it
is mirrored through the interface onto the IIB cell, and
used to set the coefficients at the IIB point so that the no
slip condition is satisfied at the surface.
The pressure equation is linearised excluding the IIB
velocities. This ensures that there is no mass flow
through the surface, significantly enhancing conver
gence and accuracy of the simulation. When all the
coefficients are known the matrix is assembled and the
NavierStokes equations are solved. At this stage the
flow field is fully determined.
In order to calculate forces acting on the IB, the pres
sure and velocity need to be extrapolated onto the trian
gulated surface of the body. For each surface triangle,
a central point is found. This point is projected into the
flow domain in the direction perpendicular to the sur
face. The pressure and velocity at the point are obtained
using trilinear interpolation from neighboring cells and
subsequently extrapolated onto the corresponding trian
gle. The total force acting on the particle is found by
integrating the sum of pressure and viscous forces act
ing on each triangle of the surface.
In the case of a moving particle the motion of the cen
tre and rotation of the body are resolved using the New
ton second law and the new position of the particle is
found. After this stage the simulation moves to the next
time step. The method is combined with a fully cou
pled incompressible pressurevelocity solver, allowing
one outer iteration per timestep. The time discretisation
Outside the 1B
1 1 0O
2 2 1 0
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
is second order backward Euler.
Test cases
Uniform flow past a rotating sphere
This test case describes the flow past a rigid sphere
rotating with angular velocity c. Since in most multi
phase flow applications the particles are approximated
as spheres, solution of the problem has a large impor
tance. Depending on the particle Reynolds number the
character of their interaction with the flow can vary. In
dilute flows where the distances between particles are
large and the motion of the body is controlled mainly by
the surrounding fluid, the particle is unlikely to achieve
high spinning velocity, therefore forces induced by the
rotation are small and rarely taken into account.
However, if a dense flow or a flow in proximity of the
wall is considered, the resulting collisions may lead to
high transverse rotational velocities achieved by the par
ticles. This in turn can change the behaviour of flow past
a particle significantly introducing additional lift force
and causing transition of the flow character at smaller
Reynolds number as shown in (5). Hence it is important
to understand how the rotation of the particle influences
the particleflow interaction.
One of the first studies in this topic was performed
by Rubinow and Keller (12). This work derives a for
mula relating the lift force to the rotational speed of the
particle for low Reynolds numbers by combining the
solutions near the sphere and at infinity. The resulting
compact formula for lift coefficient experienced by the
sphere is given by C1 2c*, where w* R= is the
dimensionless angular speed, c is the rotational velocity
of the sphere, R is the particle diameter and Uo is the
free stream velocity. It has also been shown that at small
Reynolds numbers (Re < 1) the rotation does not affect
the drag coefficient acting on the sphere.
Uniform flow past a sphere rotating transverse di
rection is also analysed numerically by other authors
in (3), (11), (6), (5). This set of papers provides ex
tensive data for the force coefficients along with be
haviour of the wake for the flow past a rotating sphere
in range of Reynolds numbers from 1 up to 500 with
nondimensional angular velocities of 0 < c* < 1.2.
Additionally, experimental investigation of the topic
can be found in an article by Oesterl6 (1), where for
mula for approximated lift coefficient is obtained for the
Reynolds numbers from 10 to 140 and for dimensionless
angular velocities c* from 1 to 4.
CL 0.45 + 2c* 0.45)
exp(0.075w 04Re0'7) (1)
This paper presents the results for the simulations of
the flow around a rotating sphere. Two different flow
regimes, defined by the particle Reynolds number, are
analysed to show ability of the MIB method to success
fully capture the forces acting on a sphere. Firstly, the
low Reynolds number range where the lift coefficients
can be found from the solution introduced by Rubinow
and Keller (12). Then simulations at Re 100 are per
formed and compared with other numerical studies of
the problem.
The considered flow domain is illustrated in Fig. 3.
The sphere is placed in the centre of a cube which
stretches 5 sphere diameters in each direction. A uni
form flow with velocity U, 1 m/s enters the domain
from the left hand side. The Reynolds number is ad
justed by changing the viscosity p of the fluid, while the
dimensionless angular velocity is determined by specif
ing angular speed of the sphere c. A full slip, (i.e. no
velocity gradient) boundary condition is applied on the
walls of the cube.
U
(i 
Re= U ,w*
p 2U,
Figure 3: The analysed domain for a flow past a rotat
ing sphere.
Sedimenting Sphere
The second case concerns the study of a glass sphere
settling in quiescent water, as described in (9). Grav
ity acting on particles can have a major influence espe
cially when relatively large bodies are considered. This
includes, but is not limited to, dune flows or slug flows
observed in pipes transporting solid particles.
The main issue related to analysing the motion of set
tling spheres is the fact that the Reynolds number, de
fined as Re = D, is not known a priori, as the settling
velocity of the particle Up is determined from the drag
coefficient which in turn depends on the density ratio
and the Reynolds number. Hence, even though analyt
ical solution is available for low Reynolds numbers, it
has limited applicability.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
The problem has been investigated experimentally
by several authors such as Mordant (9). Additionally,
Cheng in his work (2) reviews available results and sug
gests an explicit semiempirical formula for calculat
ing the settling velocity of the sphere based on dimen
sionless diameter D* (Ag/v2)1/3D, where A
(pp p)/p is the dimensionless density difference be
tween particle (pp) and the fluid (p):
Up = D(Ag)1/3 (2)
3CD
where the settling drag coefficient CD is found from:
432
CD *3(1 + 0.022D*3)054
+ 0.47[1 ep(0.15D*045)] (3)
The performed simulation shows the capability of the
MIB method to capture gravity effects and body forces
acting on the particle and therefore successfully predict
the settling velocity of the sphere. The accuracy of the
prediction of the final velocity of the particle will be
in investigated through grid refinement. The simulation
setup can be found in the Fig. 4.
20 D
Figure 4:
case.
P
Up
10 D
Flow geometry for the sedimenting sphere
The computational domain is ten sphere diameters
wide and is at least 20 diameters long. The particle is
placed close to the top boundary in the centre of the of
the domain with zero initial velocity. It is allowed to ac
celerate due to gravity until a uniform settling velocity
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Up is achieved. The vertical size of the domain is of suf
ficient length to achieve a steady state. The simulation
parameters (viscosity and density of the fluid, particle
diameter and density) are chosen to match those of the
experiment (9), resulting in Re 41. The pressure at
the walls is set by extrapolation while a full slip bound
ary condition is applied to the velocities.
Results
The presented observations are obtained by solv
ing the linearised NavierStokes equations including the
mirrored IB conditions in a fully coupled framework. A
second order accurate central scheme is applied in space
while second order accurate twostep backward Euler
scheme is used for time discretisation. The fluids are
assumed to be incompressible and eulerian, while the
particles described by the IB are rigid, i.e. no surface
deformation is considered.
Uniform Flow past a rotating sphere
The analysis for the uniform flow past a rotating
sphere are performed with Re = 0.5 and a dimension
less rotation velocity of c* = 0.5. The number of com
putational cells is varied to establish the accuracy of the
method.
Fig. 5 presents the development of predicted drag and
lift coefficients for different grid refinements as a func
tion of time. After short transition period the force coef
ficients achieve constant values, which converge asymp
totically as the grid is refined. For the finest mesh spac
ing, the drag coefficient is CD 46 while the lift coeffi
cient approaches the value of CL 1.05. These values
are in very good agreement with the results predicted
from theory of Stokes Flow, (i.e. CD = j = 48) and
results found by Rubinow and Keller (12) (i.e. CL
2wc* 1).
The second order convergence of the lift coefficient
can be clearly identified when considering Fig. 6. One
can achieve reasonable estimation of the forces acting
on the sphere by using a mesh as dense as 80 cells in
each dimension. Further increase of the number of cells
will result in more accurate results, but simulations with
more then 100 mesh cells in each direction are much
more expensive while the precision no longer increases
significantly.
In addition to aforementioned accuracy study, sepa
rate analysis are realized and compared with the theoret
ical predictions. The considered flow domain is discre
tised into 80 uniform mesh cells in each direction. Five
Reynolds number settings, from Re = 0.1 to Re =0.5
are investigated, along with 3 values of w*. The results
are presented in Fig. 7 and Fig. 8. As expected and
[ff 40
* 60
* 80
100
125
0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 A
Time [(s
Figure 5: The lift coefficients development in time for
flow at Re = 0.5 and c* = 0.5.
predicted by Rubinow and Keller, the rotation does not
have any influence on the drag coefficient exhibited by
the sphere, and the change of Reynolds number has no
effect on the value of the lift force, hence the correspond
ing plots are not shown. On the other hand the drag
coefficient changes according to the theoretical formula
CD = = 48, while the lift is directly proportional
to the nondimensional angular velocity as predicted by
(12), CL 2W*.
The small difference between the numerical results
and theoretical values can be explained by the influence
of boundary conditions applied to the flow. Theoretical
considerations assume an infinite domain surrounding
the sphere, while in the case of simulations the compu
tational box size is limited. In this work the domain it is
set to reach 5 sphere diameters in each direction, a com
promise between computational speed and accuracy. As
the Reynolds number is decreased to very small values
the influence of the box size becomes more significant
as it can be observed in Fig. 7.
Second order accuracy of the mirroring IB method
in predicting the lift coefficient is also shown for the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Number of grid cells in each direction I]
Figure 6: The variation of lift coefficient for flow at
Re = 0.5 and w* = 0.5 as a function of grid refinement
level in each direction.
8.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50
Reynolds Number []
Figure 7: Variation of the drag coefficient as a func
tion of Reynolds number in the limit of Stoke's flow
with w* 0.5. The dashed line represents theoreti
cal prediction for Stokes flow. Solid line represents the
numerical results.
Reynolds number Re 100 and w* 0.25. Smilarly
as above the change in Re is obtained by varying the vis
cosity of the fluid. As expected, after a transient period
the force coefficients achieve a constant value, what can
be observed in Fig. 9, where the drag and lift coefficients
are plotted as a function of time. Second order conver
gence can be also seen in the case of the lift coefficient
as the values approach a value of CL = 0.24, what is
presented in Fig. 10.
The comparison with the other works is shown in Ta
ble 1. The values obtained from the simulation are sim
ilar to those found in literature. It is worth mentioning
that the current results are achieved with much coarser
grid spacing then used in cited computational work.
1.4
D.2 0.3 0.4 0.5 0.6
Dimensionless angular velocity []
0.7 0.8
Figure 8: Variation of the lift coefficient as a function
of angular velocity w* for a flow at Re = 0.5. dashed
line represents theoretical prediction as in (12).
Table 1:
literature
Comparison of the results with the data from
In addition, a separate simulation with Re = 100,
* = 1 on a 100 cubed grid is performed and compared
with experimental data from (1). The lift coefficient de
velopment in time is presented in Fig. 11 in combina
tion with value predicted by formula (1) introduced by
Oesterl6.
The equation 1 gives a quite reasonable approxima
tion of the expected values. However, it has to be men
tioned that the experimental results are highly scattered
and experiment for these specific conditions has not
been performed. Even though, the simulations results
are in good agreement with empirical predictions.
Sedimenting Sphere
In order to analyse the velocity of a sphere settling in
quiescent fluid, the motion of a glass sphere with diam
eter D = 0.5mm sedimenting in water tank is resolved.
Fig. 12 presents the velocity of the particle as a function
of time for different grid refinement levels in comparison
with experimental results of Mordant (9). As the particle
accelerates from zero velocity, the drag starts to increase
as well. When the gravity force becomes balanced by the
hydrodynamic forces, an equilibrium is achieved and a
constant settling velocity is observed.
 Rubinow & Keller
 Simulations
 Stokes Flow
' Simulations
Reference CD CL
Giacobello (6) 1.15 0.29
Kim (5) 1.1 0.25
Niazmand (3) 1.12 0.25
Kurose (11) 1.1 0.23
Current work 0.9 0.24
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0.20
D 50 60 70 80 90 100 110 120 1
Number of grid cells in each direction []
Figure 9: The drag and lift coefficients development in
time for Re = 100.
The obtained value of Up = 0.073m/s is also in a
good agreement with the velocity predicted by the equa
tion 2: Up 0.077m/s.
For coarser grids, the velocity tends to be overesti
mated. However, as the grid is refined, the results are
getting closer to the experimental values. When the
number of grid cells is very low (40 grid cells in horizon
tal direction), wiggles in velocity are observed. They are
result of the difference in the levels of grid refinement of
the fluid and the triangulated sphere surface. The effect
of changing the amount of computational cells can be
seen in Fig. 13. Again, the second order accuracy in
predicting the settling velocity is observed.
Conclusions
The article presented ability of the mirroring im
mersed boundary MIB method to capture various fea
tures of the multiphase flows. In the described method,
the interface is implemented in an implicit way, by in
troducing an artificial flow field inside of the particle,
causing an exact noslip condition on the surface of the
Figure 10: The change of lift coefficient with grid re
finement for Re = 100.
0.5 1.0 1.5 2.0 2.5 3.0 3.5
Time [s]
Figure 11: The lift coefficients development in time for
Re 100 and w* = 1.
immersed boundary. This artificial flow field is excluded
from the continuity equation to improve mass conser
vation and the stability of the method. The method
is shown to be second order accurate in predicting the
forces acting on a body moving in a fluid for two dif
ferent cases: calculating the lift coefficient of a rotating
sphere in a uniform flow and finding the final velocity of
a sphere settling in a fluid. In both cases, the obtained
results were in good agreement with the values found in
literature.
References
[1] T Bui Dinh B. Oesterle. Experiments on the lift
of a spinning sphere in a range of intermediate
reynolds numbers. Experiments in Fluids, 25:16
22, 1998.
[2] N.S. Cheng. Comparison of formulas for drag co
0.00
E 0.05
8o.1o
g 0.15
Time s]
Figure 12: The velocity profile of a glass sphere sedi
menting in water tank.
simulation
0.100  experimental
0.095
S0.090
0.085
0.080 _
0.075  
0,070  
40 50 60 70 80 90
Number of grid cells in horizontal directions I]
100
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
[7] A. Mark and B.G.M. van Wachem. Derivation and
validation of a novel implicit secondorder accu
rate immersed boundary method. Journal of Com
putational Physics, 227:66606680, 2008.
[8] J. MohdYusof. Combined immersedboundary/b
spline methods for simulations of flow in complex
geometries. Technical report, Center for Turbu
lence Research, Annual Research Briefs, 1997.
[9] J.F. Pinton N. Mordant. Velocity measurement of
a settling sphere. The European Physical Journal
B, 18:343352, 2000.
[10] C.S. Peskin. Numerical analysis of blood flow
in the heart. Journal of Computational Physics,
25:220252,1977.
[11] S. Komori R. Kurose. Drag and lift forces on a
rotating sphere in a linear shear flow. Journal of
Fluid Mechanics, 384:183206, 1999.
[12] S.I. Rubinow and J.B. Keller. The transverse force
on a spinning sphere moving in a viscous fluid. J
FluidMech., 11:447459, 1961.
Figure 13: The lift coefficients development in time for
Re 100 and w* = 1.
efficient and settling velocity of spherical particles.
Powder Technology, 189:395398, 2009.
[3] M. Renksizbulut H. Niazmand. Surface effects on
transient threedimensional flows around rotating
sphere at moderate reynolds numbers. Computers
& Fluids, 32:14051433,2003.
[4] H.H. Hu. Direct simulation of flows of solidliquid
mixtures. Int. J Multiphase Flow, 22(2):335352,
1996.
[5] Kim. Laminar flow past a sphere rotating in the
transverse direction. Journal of Mechanical Sci
ence and Technology, 23:578589, 2009.
[6] S. BalachandarM. Giacobello, A. Ooi. Wake struc
ture of a transversly rotating sphere at moderate
reynolds numbers. Journal of Fluid Mechanics,
621:103130,2009.
