NUMERICAL SIMULATION OF PARTICLE-LADEN FLOWS
IN CURVED PIPES
By
Pusheng Zhang

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2012

ABSTRACT
NUMERICAL SIMULATION OF PARTICLE-LADEN FLOWS IN CURVED PIPES
By
Pusheng Zhang
Particulate flows in curved pipes are commonplace in a variety of settings and involve aerosol
deposition in airways as well as solids or droplets impacting pipe walls. The presence of bends in
such flows is well-known to be associated with large pressure drops and phase separation due to
the centrifugal force. The consideration of bends is thus important in designing multiphase flow
equipment and pipes. For example, natural gas transported in pipes is often wet and water
droplets, as they pass through the bend, may impact the bend wall where they form a film; this
may well influence the performance of separation equipment located downstream. Phenomena
associated with flows through curved pipes, such as the formation of secondary flow patterns, are
thus studied in this work in order to evaluate their effects on the disperse phase. In addition, the
accuracy of the numerical simulations for such flows is often in question when compared with
experiments. This may be due to a number of factors which include the selection of an
appropriate multiphase model and turbulence closure, or numerical aspects such as the quality of
the treatment of the near wall behavior.
In the first part of this work, dilute particle flows in curved pipes are modeled using one-way
and two-way coupled models. In one-way coupling simulations, the influence of the particles on
the carrier phase is ignored. The influence of turbulence closures on the flow and particle
trajectories are investigated. The “standard” k-ε model and the Reynolds Stress Model (RSM)
based on the Reynolds-Averaged Navier-Stokes (RANS) equation are employed with different
near-wall treatments. For two-way coupling simulations, the drift flux model based on the

mixture theory is used to consider the interaction between the phases. A realizable k-ε model is
employed to close the RANS equation and the Enhance Wall Treatment (EWT) is applied for the
flow in the near-wall region. Results show that the pressure drop of a single phase flow along the
curved pipe is well predicted by the turbulent closures studied.
For one-way coupled simulations, RSM with EWT is accurate in estimating grade efficiency
curves. Compared to other possible combinations, using RSM with EWT can improve the
accuracy by as much as 19% in a 90°bend and up to 30% in a 180°bend. The results of these
simulations have allowed the development of an improved correlation for predicting grade
efficiency curves. For two-way coupled simulations, results show that the pressure drop is
significantly affected by the disperse phase. The computed pressure has a good agreement with
the empirical correlation of Paliwado. Bend design using the mixture model shows that 90°and
180°bends with the curvature ratios equal to 5 and 7 respectively can be used to achieve a high
deposition efficiency with a relatively low pressure drop.
Following the above studies, a modification of the Immersed Boundary (IB) method, an
approach for the simulation of particles moving in a fluid, is introduced to perform a preliminary
validation of the closures used in multiphase flow modeling. An algorithm is implemented to
combine the MacCormack scheme with the IB method. The technique is applied to particulateladen flow simulations. This data is used to verify the performance of the mixture model on
estimating particle behaviors in Couette flows and Poiseuille flows. Results show that IB method
based on the MacCormack scheme is promising in dealing with fluid-structure interaction
especially for particulate-laden flows. Comparison between the DNS data and the mixture model
indicates that the mixture model is not able to capture particle migration. To improve the
performance, the lift force needs to be considered in the model used to close the slip velocity.

Dedicate this work to my parents and my wife, Shaowen.

iv

ACKNOWLEDGMENTS

I would like to give my deepest gratitude to my advisor, Dr. André Bé
nard for his
invaluable guidance, understanding, and constant encouragement during the entire period of 4
years of the present research at MSU. His patience and support helped me overcome many
difficulties and finish this dissertation.
I am deeply grateful to Dr. Charles Petty for his countless assistance and guidance to
always keep my research on the right track. I am also thankful to him for the many valuable
discussions that helped me understand my research area better.
I would like to gratefully and sincerely thank Dr. Randy Roberts for his insightful
comments and constructive criticisms at different stages of my research.
I would like to acknowledge Dr. Farhad Jaberi for setting high standards in his teaching
and helping his students to achieve those standards in the class.
I also want to thank Dr. Neil Wright’s help for his serving my dissertation defense
committee and reviewing the manuscript of my dissertation.
I want to thank all my instructors at MSU, who taught and encouraged me to reach this
point. I want to thank department staffs for their secretarial support on my study and life at MSU.
I would like to thank my family for their generous support, and in particular, to express
my deep appreciation to my wife, Shaowen, for her patience, assistance, continuous
encouragement, and support.
Finally, I gratefully acknowledge the support for this work from the NSF-I/UCRC grant
IIP-0934374 on Multiphase Transport Phenomena and the Chevron Energy Technology
Company.
v

TABLE OF CONTENTS

LIST OF TABLES ..................................................................................................................................ix
LIST OF FIGURES ................................................................................................................................. x
1 INTRODUCTION ................................................................................................................................ 1
1.1 Background .......................................................................................................................... 1
1.2 Objectives of this work ....................................................................................................... 4
1.3 Outlines of this dissertation ................................................................................................ 7
2 SIMULATION OF DILUTE PARTICULATE-LADEN TURBULENT FLOWS IN
CURVED PIPES USING A EULERIAN-LAGRANGIAN METHOD.................................... 9
2.1 Introduction ......................................................................................................................... 9
2.2 Modeling the Continuous Phase ...................................................................................... 11
2.2.1 Governing Equations.................................................................................................. 13
2.2.2 Near-wall Treatment .................................................................................................. 14
2.3 Modeling the Discrete Phase ............................................................................................ 14
2.3.1 Stochastic Turbulent Model ...................................................................................... 16
2.3.2 Empirical Models for Grade Efficiency.................................................................... 17
2.4 Numerical Approaches ..................................................................................................... 19
2.5 Test Cases ........................................................................................................................... 20
2.6 Results and Discussions .................................................................................................... 26
2.6.1 Flow Patterns and Secondary Flow Intensity .......................................................... 26
2.6.2 Pressure Drop along Curved Pipes ........................................................................... 30
2.6.3 Modeling Particle Deposition .................................................................................... 35
2.7 Summary ............................................................................................................................ 48
3 SIMULATION OF DILUTE TURBULENT MIST FLOW IN CURVED PIPES
USING AN EULERIAN-EULERIAN METHOD ....................................................................... 51
3.1 Introduction ....................................................................................................................... 51
3.2 Multiphase flow modeling ................................................................................................ 53
3.2.1 Multiphase turbulent flow modeling......................................................................... 54
3.2.2 Turbulence closure ..................................................................................................... 57
vi

3.2.3 Near-wall treatment ................................................................................................... 57
3.3 Numerical methodologies ................................................................................................. 58
3.4 Geometry and meshing ..................................................................................................... 58
3.5 Results and discussion ....................................................................................................... 59
3.5.1 Model assessment ........................................................................................................ 59
3.5.2 Applications................................................................................................................. 66
3.6 Summary ............................................................................................................................ 77
4 SIMULATION OF FLUID-STRUCTURE INTERACTION USING AN ELASTIC
FORCING METHOD BASED ON AN IMMERSED BOUNDARY METHOD................. 78
4.1 Introduction ....................................................................................................................... 78
4.2 The MacCormack/IB method .......................................................................................... 82
4.2.1 Elastic forcing method................................................................................................ 84
4.2.2 Construction of the Dirac delta function δh ............................................................. 87
4.3 Algorithm ........................................................................................................................... 90
4.4 Results and discussion ....................................................................................................... 93
4.4.1 Poiseuille flow .............................................................................................................. 93
4.4.2 Flow passing a stationary cylinder ............................................................................ 98
4.4.3 Particle impinging on a wall .................................................................................... 103
4.5 Summary .......................................................................................................................... 107
5 DIRECT NUMERICAL SIMULATION OF PARTICULATE-LADEN FLOW
USING AN FICFICIOUS DOMAIN METHOD BASED ON THE IMMERSED
BOUNDARY METHOD ................................................................................................................ 109
5.1 The fictitious domain method ........................................................................................ 109
5.2 Collision scheme .............................................................................................................. 111
5.3 Test cases .......................................................................................................................... 114
5.3.1 Flow passing a stationary cylinder .......................................................................... 114
5.3.2 Particle sedimentation .............................................................................................. 115
5.3.3 Particle behavior ....................................................................................................... 127
5.4 Validation of the mixture model in FLUENT ............................................................... 136
5.4.1 Transport of a non-uniform shear flow .................................................................. 137
5.4.2 Transport of a uniform unidirectional flow ........................................................... 141
5.5 Summary .......................................................................................................................... 144
6 SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS .......................................... 146
vii

6.1 Summary and Conclusions ............................................................................................. 146
6.2 Recommendations for future work ................................................................................ 149
APPENDICES ....................................................................................................................................... 152
A. Realizable k-ε equation for two-phase flows .................................................................. 152
B. Analysis of force acting on a particle .............................................................................. 154
BIBLIOGRAPHYS .............................................................................................................................. 156

viii

LIST OF TABLES

Table 2.1: Summary of the parameters used in the test case studies ............................................ 20
Table 2.2: Boundary conditions used in all simulations ............................................................... 21
Table 2.3: Mesh grid nodes used for the solution independent study........................................... 24
Table 3.1: Boundary conditions used in all simulations ............................................................... 59
Table 3.2: Summary of the parameters used in the test cases ....................................................... 66
Table 3.3: Pressure drop and deposition efficiency in a U-bend with different incline angles….73
Table 4.1: Differencing sequence for the space derivative discretization .................................... 92
Table 4.2: Mesh dependence study in space and time for flow with Re=100 ............................ 103
Table 5.1: Comparison of results from the current schemes and other methods (Re=100)........ 115

ix

LIST OF FIGURES

Figure 1.1: Flow chart shows the structure of this dissertation. ………………………………….8
Figure 2.1: Schematic display of 3D geometric parameters of a bend with a bend angle θ. ........ 22
Figure 2.2: O-type structured mesh used for a cross section; this mesh allow easy
refinement close to the wall and prevents a singularity at the center of the pipe. ........................ 23
Figure 2.3: Independent study of particle deposition using RSM with (a) SWF and (b)
EWT. ............................................................................................................................................. 25
Figure 2.4: The streamlines show the development of flow patterns in a tight bend (δ=1.5);
the secondary flows are also shown at different cross-sections. ................................................... 27
Figure 2.5: The streamlines associated with secondary flow patterns are studied at a 45°
,
90°and 135°deflection and shown in (a), (b) and (c) by the stand k-ε model and in (d), (e)
and (f) by the RSM. ...................................................................................................................... 29
Figure 2.6: Plots of the streamlines associated with the in–plane flow patterns at 135°in a
U-bend with δ=5.6 to demonstrate the effect of an increasing flow Re number. ......................... 30
Figure 2.7: A comparison between the computed (markers) and estimated (lines) pressure
drop values by using different closure models combined with EWT. .......................................... 33
Figure 2.8: A comparison between the computed (markers) and estimated (lines) pressure
coefficient (Cp) from near-wall treatments is made at the outer bend, inner bend, and the
sides............................................................................................................................................... 34
Figure 2.9: A comparison between computed and experimental data for the grade
efficiency shows a discrepancy occurs at small value of St for SWF and NEWF and great
improvement by using EWT in: (a) the 90-degree bend, (b) the 180-degree bend. ..................... 36
Figure 2.10: Comparison of the models (makers) and the empirical models (lines) shows
that the k-ε model over-estimates particle deposition when particles experience a longer
residence time in (b) the 180° bend (τ=4.4ms) than in (a) the 90° bend (τ=2.2ms) ..................... 38
Figure 2.11: Comparison of grade efficiency shows that the particle cutsize is 30% smaller
in the 180°bend due to the longer particle residence time. .......................................................... 39
Figure 2.12: Comparison of bends with a varying curvature ratio for a fixed flow Reynolds
number in a 90° bend shows the grade efficiency is increased with the increase of the
curvature ratio. .............................................................................................................................. 41
x

Figure 2.13: Comparison of results for curved pipes with different curvature ratio and fixed
particle residence time shows the grade efficiency is increase with the decrease of bend
curvature ratio. .............................................................................................................................. 42
Figure 2.14: (a) Comparison shows that the models are not appropriate for varying δ in a
90° bend (b) Comparison shows the new model match the results relatively well for δ
between 3 and 8 but has reduced accuracy at extreme curvature ratios. ...................................... 43
Figure 2.15: Comparison of grade efficiency shows a deviation of 9% when the Re number
increases by a factor of 10. ........................................................................................................... 45
Figure 2.16: Comparison results show the particle cutsize is increased as the increase of the
pipe dimension. ............................................................................................................................. 46
Figure 2.17: Estimation of deposition patterns of particles changes with Stokes number. .......... 48
Figure 3.1: The air and water streamlines show the development of the flow patterns and
the interaction between two phases. The contour of αp shows the liquid film formation
along the 180°bend. “G” and “L” denotes gas and liquid, respective. ........................................ 61
Figure 3.2: Water film forms at the inner bend wall affected by the secondary flow pattern. ..... 63
Figure 3.3: Volume fraction distribution of the water phase shows the location of liquid
film at the pipe wall of the outlet cross section. ........................................................................... 64
Figure 3.4: Comparison shows good agreement of the numerical results with the
Paliwoda’s empirical model .......................................................................................................... 65
Figure 3.5: Comparison between the single-phase air flow (b) and the air-water flows (c-f)
on the secondary flow patterns at the cross section of 90° deflection (a) changing with
different droplet size ..................................................................................................................... 68
Figure 3.6: Volume fraction distribution shows the liquid film location at the cross section
of 90°deflection changing with different droplet size ................................................................. 70
Figure 3.7: Secondary flow patterns and flow intensity at the cross section of 90°
deflection for flow with a different αp .......................................................................................... 70
Figure 3.8: Pressure drop of the bend section changing with (a) the droplet size, (b) volume
fraction of the water phase ............................................................................................................ 71
Figure 3.9: Droplet deposition efficiency changing with (a) the droplet size, (b) volume
fraction of the water phase ............................................................................................................ 72

xi

Figure 3.10: Schematic of a U-bend changing from horizontal position to vertical position
with an increment of 15° For interpretation of the references to color in this and all other
.
figures, the reader is referred to the electronic version of this dissertation
…………………...........................................................................................................................74
Figure 3.11: Tendency of bend pressure drop and droplet deposition efficiency in a
horizontal 90°bend and a 180°bend shows the bend geometry is optimal at curvature ratio
equal to 5 and 7, respectively. ....................................................................................................... 76
Figure 4.1: Schematic configuration shows a Eulerian mesh for the fluid domain Ωf and a
Lagrangian mesh around the fluid-object interface S . ................................................................ 80
Figure 4.2: Schematic of spring construction for a movable rigid structure using a
Lagrangian forcing method. .......................................................................................................... 87
Figure 4.3: Schematic configuration shows two mesh systems and the Eulerian points
involved in a spreading procedure and an interpolation procedure for one Lagrangian point. .... 89
Figure 4.4: Distribution of (r) shows how the Dirac delta function works. ............................... 90
Figure 4.5: Schematic of a flow passing two parallel planes (a) without an immersed
cylinder and (b) with an immersed cylinder ................................................................................. 94
Figure 4.6: Comparison of the velocity profiles between the analytical solution and the
numerical results cross the channel shows a very good agreement. ............................................. 95
Figure 4.7: Comparison of the vector and the vorticity between (a) a Poiseuille flow (b) a
Poiseuille flow interacted with a stationary structure. .................................................................. 96
Figure 4.8: Comparison of the velocity profiles along Line AB between two cases .................... 97
Figure 4.9: Comparison of the pressure distribution along Line CD between two cases ............. 97
Figure 4.10: Schematic of the computational domain used for flow passing a stationary
cylinder ....................................................................................................................................... 101
Figure 4.11: The developing vorticity patterns with a changing time shows the flow
eventually breaks into a periodic Ká n Vortex Street. ........................................................... 102
rmá
Figure 4.12: Schematic of computational domain for a disc settling under gravity and a
stationary wall ............................................................................................................................. 104
Figure 4.13: The vertical velocity and center location of the disc changing with time shows
the settling disc is (1) accelerated, (2) with constant velocity, (3) decelerated, and (4) at rest.
..................................................................................................................................................... 105

xii

Figure 4.14: Schematic of velocity contour shows a disc is settling under gravity and
impinging on a stationary wall. ................................................................................................... 107
Figure 5.1: Schematic of the collision scheme used to prevent collision between (a) two
particles (b) a particle and a stationary wall ............................................................................... 113
Figure 5.2: A curve shows the repulsive force acting on a rigid particle is decreased
dramatically as the increase of particle distance within the safety zone and is zero out of
the safety zone............................................................................................................................. 114
Figure 5.3: Comparison of (a) disc position and (b) disc velocity between the elastic
forcing method and the fictitious domain method shows that the fictitious domain method
(short dash) outperforms the elastic forcing method (solid line) in estimating the terminal
velocity of a disc falling in a flow with Re=100. ........................................................................ 118
Figure 5.4: Comparison of (a) disc position and (b) disc velocity between the elastic
forcing method and the fictitious domain method shows that both methods are numerically
symmetric in estimating the terminal velocity of a disc falling in a flow with small
Reynolds number (Re=10 for the fictitious domain method). .................................................... 120
Figure 5.5: Schematic demonstration of two particles (a) drafting, (b) kissing, (c) tumbling
and (d) being apart after tumbling. ............................................................................................. 123
Figure 5.6: Schematic demonstration of 240 particles doing sedimentation under gravity in
a sealed container at (a) t=0s, (b) t=2.5s, (c) t=5s, (d) t=7s, and (e) t=20s. ................................ 124
Figure 5.7: Schematic of computational domain for a neutrally buoyant elliptical disc
moving in a Stokes shear flow. ................................................................................................... 128
Figure 5.8: Comparison of results from the DNS method for different Rep and Jeffery’s
solution shows that the angular velocity and the period are well predicted by the present
DNS scheme for Rep=0.4. .......................................................................................................... 129
Figure 5.9: Schematic of the computational domain dimension [L1 L2] and the initial
location of an elliptical disc located at (1/4L1, 1/5L2) in a Couette flow with Re=700. ............ 132
Figure 5.10: Comparison of the particle location changing with time in the vertical
direction shows that particles attain the centerline in a Couette flow regardless of the size
and shape of a particle................................................................................................................. 132
Figure 5.11: Schematic of the computational domain dimension [L1 L2] and the initial
location of an elliptical disc located at (1/4L1, 7/16L2) in a Poiseuille flow with Re=700. ....... 133

xiii

Figure 5.12: Comparison of the particle location changing with time in the vertical
direction shows that particles attain approximately the location 60% from the centerline to
the wall in a Poiseuille flow regardless of the size and shape of a particle. ............................... 133
Figure 5.13: Schematic demonstration of the velocity contour and the concentration
distribution of a uniform particulate-laden flow with αp=10% in a Couette flow at t=0s,
5.5s, 100s, 125s, 135s and 160s. ................................................................................................. 135
Figure 5.14: Schematic demonstration of the velocity contour and the concentration
distribution of a uniform particulate-laden flow with αp=10% in a Poiseuille flow at t=0s,
1s, 30s, 60s, 100s and 416.0 s ..................................................................................................... 136
Figure 5.15: Initial condition shows the concentration of the dispersed phase in the
computational domain for (a) the present scheme (b) the mixture model .................................. 138
Figure 5.16: Comparison of the concentration of the dispersed phase at t=10s for (a) the
present scheme (b) the mixture model ........................................................................................ 139
Figure 5.17: Comparison of the concentration of the dispersed phase at t=25s for (a) the
present scheme (b) the mixture model ........................................................................................ 140
Figure 5.18: Comparison of the concentration of the dispersed phase at t=60s for (a) the
present scheme (b) the mixture model ........................................................................................ 141
Figure 5.19: Schematic of the computational domain for a uniform flow with αp=10%. .......... 143
Figure 5.20: Comparison of volume fraction distribution along a vertical cross-section
between two methods at t=50s and t=100s shows that the mixture model fails to predict a
higher concentration near the center for a uniform particulate-laden shear flow. ...................... 143
Figure 5.21: Comparison of volume fraction distribution along a vertical cross-section
between two methods at t=20s and t=40s shows that the mixture model fails to predict a
higher concentration close to the wall for a uniform particulate-laden pressure-driven flow
(Re=700). .................................................................................................................................... 144

xiv

CHAPTER 1
INTRODUCTION

1.1 Background
A particulate-laden flow is a two-phase flow in which one phase is continuously connected
called carrier phase and the other phase composes of discrete immiscible particles. Particulateladen flows in curved channels are commonplace in nature and often seen in a variety of setting
associated with biological, chemical, nuclear, pharmaceutical and food industries. They involve
phenomena as varied as aerosol deposition in airways or pollution control systems, as well as
solids or droplets impacting pipeline walls. Flows in most industrial applications are turbulent.
Therefore understanding of the behavior of particle turbulent flows in curved pipes is of
tremendous importance to help industries in cost saving and productivity increase.
Compared to single-phase flows, particulate-laden flows are significantly more complex. New
parameters are introduced due to the presence of particles such as the particle volume fraction αp,
mass loading φ, Stokes numbers St, and particle Reynolds number Rep. The discrete phase can
impact the pressure drop and the turbulent intensity of the carrier phase in a significant way. In
general, the pressure drop of a particulate-laden flow is increased with the increase of the particle
volume fraction (Hoang & Davis 1984). However, Marcus et al. (1990) observed that in some
-5

cases involving a small mass loading (φ<4) and transport of fine particles (in an order of 10 m),
the increasing αp can actually decrease the pressure drop due to a reduction in the gas-phase
stress. The particle diameter dp, volume fraction αp, and Reynolds number Rep have been shown
to be responsible for turbulent intensity modulation. Theofanous & Sullivan (1982) showed

1

theoretically that the turbulence intensity is a function of the particle volume fraction. Gore &
Growe (1989) provided a critical parameter d p /   0.1 , where  is the fluid integral length
scale, to predict turbulence attenuation and augment. Turbulence intensity is suppressed for

d p /   0.1 and enhanced for d p /   0.1 . Hetsroni (1989) stated that particles with Rep larger
than 400 would enhance the turbulence due to particle vortex shedding.
The Stokes number provides a measure of response time of the particles to the flow and is
defined as St  p / f where p is the particle relaxation time and f is the fluid characteristic
time scale. Particles with a small Stokes number ( St

1 ) would follow the carrier flow closely

and particles with a large St ( St  1 ) would go across the flow streamlines. Starkey (1956)
experimentally observed the non-uniform distribution of particles in pipe flows and
demonstrated that neutrally buoyant spherical particles under certain conditions have lateral
migration in a Poiseuille flow. Particle behavior such as particle clustering is not fully
understood to date.
To study particulate-laden flows, it is helpful to classify the flows into two distinct categories:
dilute flows and dense flows. Elghobashi (1991) used the mean distance S between the centers
of two neighboring particles to define a suspension flow. He stated that the flow is considered to
be dilute if S  10d p and dense if S  10d p . Crowe et al. (1998) defines a particulate-laden flow
based on the forces dominated to control the particle motion. A flow is dilute if the particle
motion is determined by hydrodynamic forces, such as drag and lift and is dense if controlled by
particle collision. A simple classification of dilute and dense flows is based on the volume
fraction of the dispersed phase αp. For industrial processes, flows are considered dilute if αp <
10%. Depending on the mass loading (or mass fraction) used for the dispersed phase, dilute
2

flows can be further classified into one-way coupling and two-way coupling. If the particle mass
fraction is small, it is reasonable to neglect the influence of the dynamics of the dispersed phase
on the carrier phase called one-way coupling. However, a dilute flow has to be considered twoway coupling if the mass fraction is so large that the influence of the dispersed phase on the
carrier phase cannot be neglected.
Recent advances on computational fluid dynamics (CFD) software and computer simulations
provide an efficient approach for studying dilute flows. Models describing particulate-laden flow
by convention can be divided into two groups: Eulerian-Lagrangian methods (Maxey and Riley
1983) and Eulerian-Eulerian methods (Manninen et al. 1996). Although all of the approaches
treat the carrier phase as a continuum, the particle phase is treated differently in these two groups.
In the Eulerian-Lagrangian approaches, the particles are marked as discrete objects in the fluid
flow. The motion of each particle is computed through a force balance equation. The EulerianLagrangian methods seem quite efficient in dealing with flows of one-way coupling. However,
the methods are shown to be computationally expensive when a large number of particles are
involved and their impact on the carrier phase becomes significant due to a high mass loading. A
more sophisticated way of dealing with two-way coupling in industrial applications is to employ
a Eulerian-Eulerian method or a two-fluid model which treats both phases as continua. The
interaction between the phases is considered in a momentum term which requires to be closed by
models.
An alternative approach of studying particulate-laden flow is direct numerical simulation
(DNS). To tackle the interaction between the phases, the no-slip boundary condition has to be
imposed on the surface of the moving particles. A popular approach in dealing with such
problems is to enforce the no-slip boundary condition directly on the structure surface which

3

requires an adaptive mesh to be built on the particle interface and extended to the fluid domain
(called body-conformal mesh). However, building a body-conformal mesh is cumbersome even
for a simple geometry like particles. This can be overcome by using a non-body conformal
Cartesian grid. Because the structure surface does not line up with the Cartesian grid,
modification to the fluid equation is required in the vicinity of the boundary. Clarke et al. (1986)
proposed a cut-cell approach which imposes the boundary conditions through a cut-cell
procedure. Majumdar et al. (2001) introduced a ghost-cell method where the boundary
conditions are imposed by fixing suitable values of the solution on the ghost cells outside the
computational domain. The most efficient and flexible method in this category is called the
immersed boundary (IB) method which considers impact of the immersed structure in the force
density term of the fluid equation. The IB method is therefore employed in this dissertation for
DNS of particulate-laden flow. Although solving a turbulent flow using DNS is still challenging
nowadays due to the limited computational power, DNS data for flows with a moderate
Reynolds number can be used as experimental data to validate multiphase flow models.

1.2 Objectives of this work
Dilute particulate-laden flow systems constitute one of the most widely used conveying
systems in industrial applications. This dissertation considers dilute flows conveying in curved
pipes. A Eulerian-Lagrangian particle tracking method and a Eulerian-Eulerian method in the
commercial software ANSYS FLUENT are employed to model one-way coupling and two-way
coupling dilute flows, respectively. Direct numerical simulations of dilute suspension flows in 2D Couette and Poiseuille flows are conducted using an immersed boundary method for model
validation.
The objectives of this work are:
4

(1) Provide computational guidelines for modeling pressure drop and particle deposition of
one-way coupling dilute gas/solid flow in the turbulent regime in curved pipes and
develop an empirical model for estimating particle grade efficiency including the impacts
of particle Stokes number, bend angle and curvature ratio.
The accuracy of using the Discrete Phase Model (DPM) (a Eulerian-Lagrangian
method) in FLUENT in predicting particle deposition of one-way coupling dilute
turbulent flows in curved pipes is often in question when compared with experiments.
This may be due to a number of factors which include the selection of an appropriate
multiphase model and turbulence closure or numerical aspects such as the quality of the
treatment of the near wall behavior. The accuracy of different closure models and the
performance of using different near-wall treatments are thus studied in order to evaluate
their effects on the pressure drop and the deposition efficiency of a particulate flow.

(2) Systematically study and discuss the feasibility of using a drift flux model based on the
mixture theory (a Eulerian-Eulerian method) in estimating pressure drop and liquid film of
two-way coupling dilute gas/liquid (mist) flows in the turbulent regime in curved pipes.
Design a curved pipe can promote film formation without causing a large pressure drop.
This work is motivated by the possibility of using computational fluid dynamics (CFD)
as a design tool applied to curved pipes feeding a gas/liquid separator. The question is to
identify the curvature of such pipes that can promote film formation upstream of the
separator and thus pre-condition the flow without creating a large pressure drop. As
regards to liquid film modeling, two-way coupling needs to be considered for a
particulate-laden flow. The one-way coupling DPM is shown accurate in estimating

5

particle deposition in Chapter 2. However, this method is not efficient and only suitable to
predict a very thin film with the thickness limited to 500 μm (ANSYS FLUENT 12.1,
2009). The mixture theory with a new drift flux model in FLUENT is considered to be an
efficient version of the Eulerian-Eulerian method in dealing with two-way coupling
particulate-laden flows. However, there is a lack of literature on the application of the drift
flux model in estimating pressure drop and liquid film formation in curved pipes.

(3) Developed an algorithm to combine the immersed boundary methods with the explicit
MacCormack solver for direct numerical simulation of particulate-laden flow. DNS data is
used to validate the performance of the mixture model in FLUENT on estimating the
particle behavior of suspension flows.
Particle migrating cross streamline in a unidirectional flow, i.e. Couette flow and
Poiseuille flow, is caused due to the lateral force induced by fluid inertia (Ho and Leal
1973). This phenomenon is known to be responsible for non-uniform concentration
distribution of particles in pipe flow. The mixture model in FLUENT has been broadly
used in industrial applications, such as sedimentation, cyclone separators and, particleladen flows. To accurately predict the particle concentration distribution, the model must
possess the capability of capturing the phenomenon of inertia-induced cross-stream
migration of small suspended particles. Nevertheless, study on the performance of the
mixture model on estimating particle migration cannot be found in literature.

6

1.3 Outlines of this dissertation
The structure of this work is shown in the flow chart in Fig 1.1. Chapter 2 of this dissertation
investigates the performance of various Reynolds-Averaged Navier-Stokes (RANS) models and
near-wall treatments with the discrete phase model on modeling one-way coupling dilute
particulate-laden flows particle deposition in curved pipes. The study focus on the flow patterns,
secondary flow intensity, pressure drop, and particle deposition. An empirical model is
developed to estimate the grade efficiency of particle depositing on the curved pipe wall
according to the particle Stokes number, the bend angel and curvature ratio. In Chapter 3, twoway coupling dilute mist flows in curved pipes are simulated using a Eulerian-Eulerian technique
referred to as a drift flux model based on the mixture theory. Bend design related the pressure
drop and liquid film formation is conducted. Chapter 4 of this work introduces an immersed
boundary method based on the elastic forcing method to tackle fluid-structure interaction. An
algorithm is developed to combine the MacCormack scheme with the immersed boundary
method for DNS of particulate-laden flow simulation. In Chapter 5, a direct forcing method
proposed by Uhlmann (2005) is used to improve the stability and the accuracy of the immersed
boundary method. The results are used to validate the performance of the mixture model in
FLUENT on estimating particle migration in a Couette flow and a Poiseuille flow.

7

Figure 1.1: Flow chart shows the structure of this dissertation.

8

CHAPTER 2
SIMULATION OF DILUTE PARTICULATE-LADEN TURBULENT FLOWS IN
CURVED PIPES USING A EULERIAN-LAGRANGIAN METHOD

2.1 Introduction
Curved pipes form essential components of piping systems of the oil and gas industry. The
presence of bends is well-known to be associated with complex flow patterns as well as large
pressure drops and to affect the performance of downstream equipment. This is especially
important for multiphase flows as the significant pressure gradient and secondary flow around
the bends would affect the phase distribution. A variety of phenomena associated with flows
through curved pipes, such as the formation of secondary flow patterns, are thus studied in this
work in order to evaluate their effects on the disperse phase.
Much work has been done on studying pressure drop, flow patterns, and particle deposition in
curved pipes. Thomson (1876) first observed the curvature effects of bends on flows. Eustice
(1910) also observed the existence of secondary flows by injecting ink into water passing
through a coiled pipe. Wilson et al. (1922) observed that the pressure drop is dependent on the
flow Reynolds number and Dean (1928) studied theoretically curved pipe flows and identified
the condition for the onset of secondary vortices. Ito (1959) found that secondary flows can
cause a rapid rise in friction and lead to a much increased pressure drop. Tunstall and Harvey
(1968) observed the presence of a main (or primary) flow recirculation at the inner wall for tight
bends (δ < 3). The literature on flow through curved pipes is vast and Berger et al. (1983)
provided a comprehensive review of this subject. The intensity of secondary flows in bends
depends on the combination of the main flow Reynolds number (Re=Ud/ν) and the curvature
ratio (δ=Rb /Rt). It can be characterized by a dimensionless number called the Dean number

9

which is defined here as De=Re/δ

1/2

. Ito (1987) demonstrated that the size of the secondary flow

patterns matches the size of the duct radius.
Based on theoretical calculations, Cheng and Wang (1975) derived a formula to estimate
particle deposition in bends but it does not account for secondary flow patterns. The empirical
model is accurate for Reynolds numbers in the range of 1000<Re<5000. Pui et al. (1987)
observed that complex flow patterns as well as turbulent fluctuations play a major role in particle
deposition patterns. They proposed an empirical model for the particle deposition efficiency for
turbulent flows in a 90°bend that is a function of only the Stokes number (St). Peters and Leith
(2003) measured the deposition efficiency of high Re number flows in large industrial curved
pipes with different bend angles. Brockman (1993) extended the empirical model of Pui et al.
(1987) by accounting for the bend angle in addition to the Stokes number. McFarland et al.
(1997), using numerical results, developed an empirical model that accounts for the bend angle,
the Stokes number, and the bend curvature ratio to estimate the deposition efficiency. In spite of
the vast amount of work performed, there is still a need to identify conditions under which
computer simulations can provide relatively accurate results.
Computer simulations provide an efficient approach for studying flows through curved pipes
under various conditions. Practical simulations can be performed by solving the filtered NavierStokes equation using a Large-Eddy Simulation (LES) for an instantaneous solution or by
solving the Reynolds Average Navier-Stokes (RANS) equation for an ensemble averaged
solution (with appropriate closure models for the Reynolds stress or the sub-grid stress tensors).
For example, Breuer et al. (2006) and Berrouk et al. (2008) used LES to study particle deposition
and their computed results are in good agreement with the work of Pui et al. (1987) except for
small Stokes number St<0.2. A RANS approach was selected for this work because of the large

10

number of simulations expected. In addition, a Reynolds stress model was selected to complete
the formulation of the RANS equation due to the anticipated presence of strong streamline
curvatures.
In this Chapter, the performance of RANS models with different near wall treatments is
evaluated on estimating pressure drop and deposition efficiency of a dilute turbulent suspension
flow through curved pipes. The quality of the results is assessed in various ways. A comparison
of the flow patterns at different pipe cross sections is made between the standard k-ε and
Reynolds stress models. The intensity distribution of secondary flows is investigated at different
cross sections of the U-bend with a fixed curvature ratio. Numerical calculations from different
near-wall treatments based on the standard k-ε model and RSM of the pressure drop are
compared to the experimental work of Sudo et al. (1998 and 2000). The performance of the k-ε
model and RSM, in combination with different near-wall treatments in estimating particle
deposition, is also examined, along with the effect of turbulent fluctuations on particle deposition.
Grade efficiency curves are compared to the experimental work of Pui et al. (1987), Brockmann
(1993) and McFarland et al. (1997). In addition, the effects of particle residence time, bend
curvature ratio, and flow Reynolds number on the grade efficiency are investigated.

2.2 Modeling the Continuous Phase
A dilute particle-laden flow with a low mass loading is considered in this Chapter. For this
kind of flows, the impact of the particles on the carrier fluid is ignorable (one-way coupling) and
no particle-particle interaction is involved. Particles move with the local gas velocity, but no
momentum and mass of the particles are transferred to their surrounding gas flow. Under this
assumption, the single-phase gas flow and the particle trajectories can be solved independently.

11

The continuous phase is solved using a Eulerian approach. The Reynolds-Averaged NavierStokes (RANS) equation governs the transport of all mean-flow properties of turbulent flows
with the range of all scales being modeled. This allows computing costs to be reduced provided
that an accurate closure model is used. The commercial software ANSYS FLUENT 12.1 was
used for the simulations presented below. The software provides several models for closing the
RANS equation and these include models based on eddy-viscosity closure approaches such as
the “standard” k-ε and k-ω models. Various formulations of the Reynolds Stress Model (RSM)
are also available; these involve solving transport equations for approximating the Reynolds
stress and an equation for the dissipation rate. To estimate the grade efficiency, Matida et al.
(2004) used the Spalart-Allmaras model and the k-ω model in a mouth-throat geometry and
found that the grade efficiency curves computed from those models are highly overestimated for
the entire range of particle diameters. The discrepancy from the results of Pui et al. (1987)
reaches the order of 15%-25%. Furthermore, the agreement to the experiment of the standard k-ε
model is observed to be worse than that of the k-ω model. To avoid using the Boussinesq
approximation (or eddy-viscosity models), McFarland et al. (1997) employed the RSM and
calculated the turbulent fluctuations using one-step correction approach of Abuzeid et al. (1991).
However, deviations from the results of Pui et al. (1987) were still as high as 17% at Stokes
number St=0.2. Although much work has been done using different closure models, the effect of
the near-wall treatments on pressure drop along curved pipes and particle deposition on the pipe
wall is rarely discussed. An accurate representation of the flow in the near-wall region is needed
to model particle deposition as turbulent fluctuations strongly affect particle trajectories at the
near-wall region. ANSYS FLUENT 12.1 allows to use the standard wall function (SWF), the
non-equilibrium wall function (NEWF) or the enhanced wall treatment (EWT) as the near-wall

12

treatments (ANSYS FLUENT 12.1, 2009). Zhang et al. (2012) used the RSM and the EWT for
near-wall treatment. Results are compared to the work of Pui et al. (1987) and show that the
EWT based on the RSM closure is promising in predicting particle deposition in a 90°bend. An
error of only about 3% is observed at small Stokes number.

2.2.1 Governing Equations
Since the flow is one-way coupling, a transient solution to the continuous flow field is not
necessary and the flow can be solved in a steady state. The governing equations for the
continuous phase include the averaged continuity equation

Ui
0
x i

(2.1)

and the RANS equation

uj

Ui
1 P

 gi 

2S ji  uju
i
x j
 x i x j





(2.2)

where Sji is the mean rate-of-strain tensor, p is the mean pressure, U represents the average
velocity (or mean velocity) profile, u' is the fluctuating velocity, ρ is the continuous phase
density, ν denotes kinematic viscosity, and g is the gravity. A model for uuj is needed to close
i
Eq. (2.2).
The standard k-ε model and the Reynolds stress model are used in this study and are not
presented here for the sake of brevity. The standard k-ε model, based on an eddy-viscosity
closure, is still commonly used and requires a transport equation for the kinetic energy k and for
the dissipation rate ε (Launder and Sharma, 1972). The RSM is more complicated as it requires

13

the solution of six transport equations for individual Reynolds stresses plus one equation for the
dissipation rate (Launder et al., 1975; Gibson and Launder, 1978).

2.2.2 Near-wall Treatment
In wall bounded flows, modeling the near-wall region is problematic as the closure models are
valid only for the core flow. The quality of the approximated near-wall region significantly
impacts the numerical solutions for the whole domain as it is a major source of vortices and
turbulence. In dealing with wall-bounded flows, a common approach for modeling the near-wall
region flow is to use semi-empirical wall functions. The standard wall function (SWF) and the
non-equilibrium wall function (NEWF) are used in this work. The SWF is a semi-empirical
formula proposed by Launder and Spalding (1974) based on the assumption that the production
of kinetic energy and its dissipation rate in the wall-adjacent cell are equal. Kim and Choudhury
(1995) provided a correction to the SWF which resulted in a two-layer-based NEWF. An
alternative method of dealing with near-wall region is to capture the flow in that region by
solving the one-equation model of Wolfshtein (1969); such an approach is called the enhanced
wall treatment (EWT). The EWT is especially recommended for flow exhibiting strong
streamline curvatures, boundary separation, or reattachment. However, extremely fine meshes
near the wall are needed which leads to more demanding computations.

2.3 Modeling the Discrete Phase
The one-way coupling Discrete Phase Model (DPM), or Lagrangian particle tracking method,
proposed by Maxey and Riley (1983) is used to track particles moving in the continuous phase
through the curved pipes. In this method, fictitious particles are released in the continuous phase

14

and treated as a discrete phase. The concentration of particles is assumed to be dilute i.e. the
particles do not interact and they have no effects on the continuous phase. The trajectory for a
particle is obtained by integrating the force balance, i.e.
g (p  )
 FD (u i  u p,i )  i
 Fs
dt
p

du p,i

(2.3)

and

dx p,i
dt

 u p,i

(2.4)

and where

FD 

1 CD Rep
p 24

(2.5)

In the above equations, u, up are the instantaneous gas and particle velocities; xp,i is the position
vector of particle i; dp is the particle diameter; ρp is the particle density; τp is the particle
2

relaxation time obtained from τp=ρpdp /18µ; Rep is the relative Reynolds number from
Rep=ρdp|up,i-ui|/ µ; CD is the particle drag coefficient (particles are assumed to be spherical and
the spherical drag law is used to compute CD (Morsi and Alexander, 1972)); Fs is an additional
acceleration term that could be important under certain circumstance (e.g., severe pressure
gradient, or great density difference).
Equation (2.3) and (2.4) are solved using a Runge-Kutta scheme. The maximum number of
time step is set to 100,000 in all cases for this work, which is sufficient to ensure that all the
particles either ended on the wall or at the pipe exit (there were no particles not accounted for in
*

the simulations). The discretized time step ∆t is obtained from ∆t=∆t / λ where λ is the step

15

length factor equal to the maximum number of time step divided by number of cross-sections
*

along the main-flow direction; ∆t is the transit time required for a particle to travel through the
current mesh cell. This value is estimated according to the current particle velocity and the mesh
size in the main flow direction following guidelines in ANSYS FLUENT 12.1 Documentation
(2009).

2.3.1 Stochastic Turbulent Model
The effects of turbulent fluctuations on the particle deposition are modeled by using a discrete
random walk model (DRWM), also called the stochastic turbulent model (Jang and Acharya,
1988). The instantaneous gas velocity ui in Eq. (2.3) is a combination of the mean gas velocity ui
and fluctuating components ui´ If turbulent fluctuations are not considered, trajectories are
.
computed by substituting the instant gas velocity with the mean gas velocity obtained from the
continuous phase calculation. To include the turbulent influence, the fluctuating components are
approximated using an ad-hoc random distribution given by:

u  1 uu; v   2 vv; w  3 ww

(2.6)

where ζ1, ζ2 and ζ3 are normally distributed random numbers with zero mean and a unit standard
deviation. The r.m.s. of fluctuation in Eq. (2.6) is not obtained from the transportation equations
for the Reynolds stresses but is assumed isotropic and equal to 2k / 3 so that the above
approximation can be used.

16

The fluctuating components are discrete piecewise parameters that are function of time. Their
random value is kept constant over an interval of time τe estimated by the lifetime of the eddy
defined as
k


e  -0.3 log(r)

(2.7)

where r is a uniform random number between 0 and 1. When this time is reached, a new value of
the instantaneous velocity is obtained by applying a new value of ζ. The performance of the
stochastic model on particle deposition is studied later in this work.

2.3.2 Empirical Models for Grade Efficiency
Grade efficiency curves can be obtained from the results of particle tracking simulations in
order to evaluate the performance of a device in separating particles. The grade efficiency is
defined as
(St) 

Mc Nc

M
N

(2.8)

where Mc is the mass of particles separated from the total particle mass M released. Since
monodisperse particles having the same density are tracked, Mc and M are replaced respectively
by Nc which is the number of particles trapped in the bend plus the straight pipe after the bend
and by the total number of particles tracked N. The grade efficiency, or deposition efficiency, is
presented in this paper as a function of the Stokes number St. The Stokes number provides a
measure of response time of the particles to the flow and is defined as the ratio of the particle
relaxation time p and the fluid characteristic time scale the system response time f . It is defined
as
17

St 

p
f



2
Cp d p /18

(d / 2) / U0

2
Cpd p U0



(2.9)

18d / 2

where C is the Cunningham correction factor (Cunningham,1910) used to consider
noncontinuum effects when computing the drag on small particles, U0 is the mean flow velocity
[m/s] and ρP denotes the particle density. For given values of U0, d, ρp and μ, the particle size is
then proportional to St

1/2

.

The computed grade efficiency curves are compared to empirical models. Pui et al. (1987)
conducted experiments for particulate flows passing through a 90° bend with an 8.51mm
diameter. Based on the experiments, they provided an empirical model used to estimate the grade
efficiency in 90-degree bends, given as

  (1 100.963St ) 100%

(2.10)

Brockmann (1993) extended Pui et al.’s work to include the effect of bend angles. The
expression of the model is





  1  e1.412St 100%

(2.11)

where θ is bend angle in radians. Both empirical models are stated with a 3% uncertainty. Later,
McFarland et al. (1997) developed a model based on their numerical results to account for the
effect of the St number, the bend angle and the bend curvature ratio. The model is given as

4.61  aSt


2
2 
 1  bSt  cSt  d St 

  100%  exp 

where a =  0.9526  0.05686δ, b 

c  0.306 

1.895

2
 ,
0.5 

d

(2.12)

0.297  0.0174
1  0.07  0.01712

,

0.131  0.0132  0.0003832
1  0.129  0.01362
18

In this work, particles mimicking droplets are released at an upstream position of 3d prior to
the entrance of the bend (z´
=0). At this location, the flow can be treated as fully developed and
the bend effect on particles can be neglected. 25,000 monodisperse spherical droplets are
3

released at the mean flow velocity U0 with the density 895 kg/m . The density used here is
consistent with that used in Pui et al.’s experiment. The particles are released at 3D upstream
from the bend inlet and injected normal to the cross section in a square region (the sides of the
square is 0.5d). This location keeps the particle away from the wall. Meanwhile, it avoids the
influence of the bend and allows the particle to well mix with the continuous phase before the
flow enters a bend.

2.4 Numerical Approaches
The governing equations for the continuous phase are solved assuming an adiabatic, steady
flow in FLUENT using the pressure-based segregated algorithm, Semi-Implicit Method for
Pressure-Linked Equations (SIMPLE) put forward by Patankar and Spalding (1972). The
transport equations for k, ε and the Reynolds stress equations (for the RSM) are spatially
discretized using a second-order upwind scheme. This is combined with a correction equation
derived from the continuity equation and the momentum equation for pressure is discretized
using the PRESTO! scheme (PRESsure STaggering Option) (ANSYS FLUENT 12.1, 2009)
which is well-suited for flows with strong streamline curvatures. The stability of the numerical
computations is improved by reducing the under-relaxation factor to 0.2 for the pressure and 0.5
for the momentum. The residual criterion for convergence (the “monitor”) is set to 10

-5

for all

variables. Two monitors for the difference of the mass flow rate between the inlet and the outlet
of a pipe and the area-average weighted pressure at the bend exit are used to ensure a converged
19

solution. Finally, to solve Eq. (2.3) for the discrete phase, a 5

th

order Runge-Kutta scheme is

used (Cash and Karp, 1990).

2.5 Test Cases
The parameters used in three test studies are summarized in Table 2.1. The geometric
configuration and the Reynolds number are selected based on the experimental setups available
for comparison. A tight bend with a curvature ratio δ=1.5 and a larger bend with δ=5.6 are used
to study the curvature effect on the flow patterns. In addition, the effect of the flow Re number
on the flow patterns is investigated by changing the flow velocities. The local pressure is
estimated numerically from different near-wall treatments based on the standard k-ε model and
RSM in a 90°bend and a 180°bend and compared to the experimental work of Sudo et al. (1998
& 2000). The performance of the near-wall treatments in estimating particle deposition is
evaluated by comparing the grade efficiency to the empirical work from the literature. Varying
the bend curvature ratio and flow residence time is achieved by changing the bend radius and
bend angles. Finally, the effect of the pipe dimension on the particle deposition is studied by
increasing the pipe diameter from d=8.51mm to d=100mm while maintaining the Re number, the
curvature ratio and the bend angle.

Table 2.1: Summary of the parameters used in the test case studies

Diameter, d [mm]
Curvature Ratio, δ
Bend Angle, θ [o]
Re Number Range
Near-Wall Treatments
Closure Models

Qualitative Study
on Flow Patterns
8.51
1.5,5.6
180
10,000~ 50,000
EWT
k-ε model, RSM

Study on Pressure
Drop
104
4.0
90,180
60,000
SWF,NEWF,EWT
k-ε model, RSM
20

Study on Particle
Deposition
8.51
1.5~20
60~180
10,000~100,000
SWF,NEWF,EWT
k-ε model, RSM

Table 2.2: Boundary conditions used in all simulations

Outlet

turbulent intensity [%]
turbulent viscosity ratio
pressure outlet [Pa]

Calculated from the
known Re and d
0.05
1.0
Pstatic=0

Wall

no-slip boundary [m/s]

uwall=0

velocity inlet [m/s]
Inlet

In Fig. 2.1, the geometry of a curved pipe is shown with an unspecified bend angle θ. The
outer and inner bends of the pipes are labeled as “O” and “I”, respectively. The boundary
conditions in all simulations are specified to be uniform velocity at the pipe entrance. To ensure
that the turbulent flow is fully developed before it enters the bend, a long entrance length l is
roughly estimated using Eq. (2.13) (Spurk, 1997).
1/6

l=4.4Re

(2.13)

21

Figure 2.1: Schematic display of 3D geometric parameters of a bend with a bend angle θ.

At the pipe wall, a no-slip boundary condition is applied. The wall is “sticky” so that particles
adhere to it on contact, which represent the conditions used in Pui et al.’s experiments (1987).
Numerically, the trajectory calculations are terminated when particles enter the cell adjacent to
the boundary. The boundary condition for the outflow is a pressure outlet with the static pressure
equal to zero. The applied boundary conditions are summarized in Table 2.2. No heat transfer is
considered in the process. Since the direction of the inlet flow is normal to the inlet cross section
and its velocity is uniform, the turbulent intensity and turbulent viscosity ratio at the inlet are set
to be small as recommended in ANSYS FLUENT 12.1 Documentation (2009).
22

Figure 2.2: O-type structured mesh used for a cross section; this mesh allow easy refinement
close to the wall and prevents a singularity at the center of the pipe.

The mesh construction for each cross section along the curved pipe is shown in Fig. 2.2. The
+

density of the mesh is a function of the near-wall treatment used. The dimensionless distance y

= yuτ / ν is used to describe the distance from the wall, where y is the real distance of a location
in the calculation domain from the wall; uτ denotes the friction velocity and ν is the kinematic
viscosity of air. For SWF and NEWF, simulations using the distance of the first cell centroid
+

from a wall has to be laid above the buffer layer at y between 30 and 300 so that those semiempirical equations are relatively accurate. For EWT, simulations using the centroid of wall-

23

+

+

adjacent cells is made at y ≈ 1 and at least 3 grid nodes are put inside the viscous sublayer (y <5)
to ensure the accuracy of the numerical solution. A higher Re could cause a larger friction
+

velocity and hence increase the y value at the wall-adjacent cells. Therefore, a refinement to the
+

wall-adjacent mesh size is needed to meet the y ≈ 1 requirement according to the flow Re used.
Along the mainflow direction of the pipe, the number of cross sections is changed according to
the length of the pipes. To study the simulation convergence of the near-wall treatments in
estimating particle deposition, different mesh density for pipe cross-sections and for the mainflow direction is used for SWF and EWT based on RSM. To reduce the computational effort,
solution from the entrance length is used as the inlet boundary condition and only the segment
after the entrance length (shown in Fig. 2.1) is refined for investigation. The mesh used for the
grid independence study is summarized in Table 2.3. The mesh size is doubled at the crosssections as well as the main-flow direction. For example, SWF-Mesh1 uses 896 grid nodes for
one pipe cross-section and 156 in the main-flow direction. For SWF-Mesh2, it possesses 1716
and 325 for one pipe cross-section and the main-flow direction, respectively. Meanwhile, the
+

refinement has to maintain a mesh that meets the y requirement of the utilization of SWF and
EWT.
Table 2.3: Mesh grid nodes used for the solution independent study

SWF-Mesh1
SWF-Mesh2
SWF-Mesh3
EWT-Mesh1
EWT-Mesh2
EWT-Mesh3

Number of nodes for Number of cross-sections along
Total mesh grid nodes
one cross-section
the main-flow direction
896
156
139,776
1716
325
557,700
3441
650
2,236,650
1044
156
162,864
2093
325
678,132
4096
650
2,662,400

24

(a)

(b)

Figure 2.3: Independent study of particle deposition using RSM with (a) SWF and (b) EWT.

25

Before performing computations on pressure drop, flow patterns and particle deposition, a
mesh independence study was performed. The impact of the mesh on the particle deposition,
which appears to most sensitive to the mesh, are reported below. The meshes shown in Table 2.3
are used to construct a 90° bend with δ=5.6 and d=8.51mm. The RSM with SWF and EWT is
employed to compute the flow field. The resulting grade efficiency curves are compared to the
empirical model of Pui et al. (1987) in Fig. 2.3. It can be seen that the SWF-Mesh2 and EWTMesh2 are appropriate i.e the solution does not change with a refined mesh. A comparison of the
empirical model of Pui et al (1987) with the results obtained from SWF-Mesh2 and EWT-Mesh2
show a 22% and 3% discrepancies respectively at St=0.001. In the following, Mesh2 is used in
the computations.

2.6 Results and Discussions
The complexity of the flow patterns observed in curved pipes and secondary flow intensity
are first discussed below, followed by a comparison between the computation and the
measurement of pressure drops and grade efficiency.

2.6.1 Flow Patterns and Secondary Flow Intensity
The flow patterns observed from the numerical simulations are presented in Fig 2.4-2.6. The
streamlines are used to represent the flow pattern only. No information of flow intensity is
indicated by the streamline density. The streamlines at the middle plane of a tight bend with
δ=1.5 are shown in Fig. 2.4, in which a strong recirculation of the primary flow can be observed
on the inner bend as the fluid is passing through the tight bend. Secondary flow patterns can be
seen at various cross sections of the bend (at 0° 45° 90° 135° 180°of deflection and after the
,
,
,
,
recirculation, respectively). No secondary flows are observed at the 0°deflection cross-section.
26

A pair of counter-rotating vortices is found at the 45°deflection and becomes distorted at 90°
.
The pair then evolves or split into two pairs of vortex cells at the 135°deflection and becomes
more distorted at the 180°deflection. The additional pair of vortices at the inner bend is caused
by the recirculation which generates a low velocity zone at the inner bend. The recirculation
leads to a high pressure drop at the bend. Moreover, the recirculation affects the downstream
flow patterns and partially formed Dean vortices are observed.

Figure 2.4: The streamlines show the development of flow patterns in a tight bend (δ=1.5); the
secondary flows are also shown at different cross-sections.
27

No primary recirculation (along the main flow direction) is observed in a bend with δ=5.6
from the simulations performed with the standard k-ε model and RSM. Numerical results of flow
patterns as well as secondary flow intensity are presented in Fig. 2.5. Flow patterns situated at a
45° 90°and 135°deflection angles estimated by the standard k-ε model combined with EWT are
,
shown in Fig 2.5a, 2.5b and 2.5c and compared with those simulated by the RSM with the same
near-wall treatment at the corresponding locations (shown in Fig 2.5d, 2.5e and 2.5f). The flow is
turbulent with a Re=10,000. The streamlines at all cross sections show pronounced secondary
flows with faster flowing central part of the flow. These results are consistent with the
observations of Berger et al. (1983), who found that the pressure gradient caused by the
centrifugal force around bends drives the flow near the wall to move inward, while the fastermoving fluid in the core is swept outward.
The graphs at the 45°deflection show that very different patterns are observed for RSM and
the standard k-ε model. The RSM estimates two pairs of counter-rotating vortices induced by the
Coriolis and centrifugal forces, but the standard k-ε model only estimates one pair. At the 90°
deflection, the standard k-ε model estimates an extra pair of vortices at the inner bend, which is
very weak. In addition, vortices simulated by the RSM are more distorted than those estimated
by the other model. At the 135°deflection, the vortices by RSM are located very close to the
wall and highly distorted, which is not the case for the secondary flow estimated by the standard
k-ε model. Experimental measurements performed in a 180˚ bend show that curvature effects are
responsible for producing anisotropic turbulence in the core of the turbulent flow (Chang et al.,
1983). Ito (1987) observed that anisotropic turbulence cannot be reproduced by calculations
performed with a k-ε model using a standard Boussinesq model for closing the Reynolds stress.

28

The selection of the closure model is probably the reason for the differences in the flow patterns
estimated by the two closure models.

Figure 2.5: The streamlines associated with secondary flow patterns are studied at a 45° 90°and
,
135°deflection and shown in (a), (b) and (c) by the stand k-ε model and in (d), (e) and (f) by the
RSM.

In Fig. 2.6, flows with a Reynolds number equal to 10,000, 30,000 and 50,000 are used and
the flow patterns at a 135°deflection of the same bend are compared. Instead of two cells of
vortices for flow with Re=10,000, four and six cells are observed at the higher Reynolds number.

29

Figure 2.6: Plots of the streamlines associated with the in–plane flow patterns at 135°in a Ubend with δ=5.6 to demonstrate the effect of an increasing flow Re number.

2.6.2 Pressure Drop along Curved Pipes
Calculations of the pressure coefficient (Cp) along the 90° bend and the 180°bend are
compared to the experiments of Sudo et al. (1998 and 2000) in which the pressure coefficient is
computed with
1
Cp  (P  Pref ) / ( U02 )
2

(2.14)

where P denotes the local static pressure, Pref represents a reference pressure and ρ is air density
3

(1.225 kg/m ). The flow has a mean velocity of 8.7 m/s and is passing through a bend with a
radius ratio of δ=4. The flow is turbulent with Re=60,000. The local static pressure is estimated
by using the standard k-ε model and RSM with different near-wall treatments.
Comparison between the numerical results and the experimental pressure along the outer bend,
inner bend, and at the center are shown in Fig 2.7 and 2.8. The numerical static pressure is
extracted from the locations slightly inside the outer bend wall (1 μm away from the wall), the
inner bend wall and right at the center, respectively. The variables z' and z are the pipe
longitudinal coordinates used to locate a cross section along the pipe direction and d is the pipe
30

diameter. The cross sections at the beginning and the end of the bend are located at z'=0 and z=0,
respectively. The reference pressure is the pressure located at the outer bend at the cross section
z' =-17.6d from the beginning of the bend. Comparison of the local static pressure estimated by
the standard k-ε model and RSM combined with EWT in the 90-degree bend and the 180-degree
bend is provided in Fig. 2.7. In Fig. 2.7a, the local pressure through the bend estimated from the
standard k-ε model and the RSM is plotted. The lines obtained from the standard k-ε model and
the RSM are hardly distinguishable. For the 90°bend, the maximum discrepancy between the
experimental Cp and the numerical Cp estimated by both models is located at the 45°deflection
of the inner bend wall, the numerical error is 7%. In the 180°bend, the performance of both
models in estimating the local pressure worsens. On the outer bend wall, the maximum error
occurs at 60°deflection with discrepancies of 21% and 26% for the standard k-ε model and the
RSM respectively. On the inner bend wall, the maximum discrepancy is observed at the 75°
deflection for both models with the same 26% error. Simulations results from both models with
EWT match the empirical model curve well in the 90-degree bend, but underestimate the local
pressure at the inner wall and the outer wall of the bend portion in the 180-degree bend, as
shown in Fig. 2.7b. This seems to indicate that both models are performing well in short bends
(small residence time), but for flows with a longer residence time the local pressure is not
estimated accurately.
In Fig. 2.8a, all three near-wall treatments are observed to be capable of modeling the local
pressure along the outer wall and the side wall of the 90°bend. At the inner bend, however,
results computed from the SWF and NEWF overestimate the local pressure at the inner wall of
the straight pipe after the bend. A large discrepancy of 34.5% is observed at the inner pipe wall
located at the downstream 0.5d. The plot shows that results from EWT have a better agreement

31

at those locations. For the 180°bend in Fig. 2.8b, SWF and NEWF seem to improve the results
along the inner wall and the outer wall of the U-bend portion. The maximum errors at the 60°
deflection of the outer bend and at the 75°deflection of the inner bend are reduced to 18% and
5%, respectively. Nevertheless, all methods estimate an accurate pressure drop when compared
to the experiment measurement. Hence, they all appear suitable for bend design in estimating
pressure drop.

32

(a)

(b)

Figure 2.7: A comparison between the computed (markers) and estimated (lines) pressure drop
values by using different closure models combined with EWT.
33

(a)

(b)

Figure 2.8: A comparison between the computed (markers) and estimated (lines) pressure
coefficient (Cp) from near-wall treatments is made at the outer bend, inner bend, and the sides.

34

2.6.3 Modeling Particle Deposition
Grade efficiency curves for the 90°and the 180°bends computed from the RSM, combined
with different wall functions and the stochastic turbulent model, are plotted against St and shown
in Fig. 2.9 along with the results computed from the empirical models. The results from SWF
and NEWF are in relatively good agreement for large St. However, when the value of St
becomes small, a large discrepancy from the experiments can be observed. Fig. 2.9a shows that
the maximal discrepancy (of 22%) occurs for a very small value of St (St=0.001) for both SWF
and NEWF in the 90°bend. In the 180°bend, Fig. 2.9b shows a 40% discrepancy for SWF and
43% for NEWF. Numerical error of SWF and NEWF increases when a flow is under a high
streamline curvature for a longer travelling time of the particles in bends. Small particles are very
sensitive to turbulent fluctuation because of their fast response to the system. The failure of the
SWF and NEWF to adequately capture the turbulent fluctuations near the wall may explain this
deviation. On the other hand, results computed using an EWT agree very well with the
experimental data and the maximum deviation of the comparison is only 3% at St=0.001.
Comparison in Fig. 2.9b of the grade efficiency between numerical simulation and the models in
the U-bend shows a higher deviation (of 10%) than that in the 90°bend at St=0.001. But a great
improvement by using EWT can be observed in the Figures. The EWT appears to capture more
accurately the actual velocity profile in the boundary layer especially if separation and strong
streamline curvatures occur. The SWF and the NEWF were formulated from empirical data
obtained from a straight bend and are probably not appropriate for particle deposition in flows
with strong streamline curvatures. In addition, since SWF and NEWF use thicker wall-adjacent
cells, the odds increase for small particles to hit the wall-adjacent cells and be counted as

35

(a)

(b)

Figure 2.9: A comparison between computed and experimental data for the grade efficiency
shows a discrepancy occurs at small value of St for SWF and NEWF and great improvement by
using EWT in: (a) the 90-degree bend, (b) the 180-degree bend.
36

“deposited”. Therefore, a more accurate estimation of small particle deposition is expected when
compared to SWF and NEWF. All the near-wall treatments perform well on estimating the grade
efficiency for large particles. With the particles become larger, the influence of the continuous
flow becomes minor due to their large inertia.
To see whether using the stochastic model is necessary, particle tracking is also run without
considering the effect of turbulent fluctuation by turning off the stochastic model. Comparison in
Fig. 2.9a shows a very poor agreement in the 90°bend without using the stochastic model. Small
particles are especially sensitive to the fluctuating velocity because of their low inertia, which
account for the great discrepancy of the model when not considering particle dispersion.
In Fig. 2.10, a comparison is presented between the grade efficiency curves obtained from the
numerical results and the models. The standard k-ε model and RSM are used in the computations
in combination with EWT and the stochastic model. In Fig. 2.10a, both grade efficiency curves
obtained from the models agree with the experiment well in the 90° bend. It is difficult to
identify which model is more accurate in estimating the grade efficiency. However, the k-ε
model estimates a maximum deviation of 10% at St=0.05 while the RSM obtains 3% at St=0.001.
Moreover, the error averaged over the range of St is 4.5% for the RSM, which is slightly lower
than that for the k-ε model. Apparently, one can see from the Fig. 10.b that the performance of
the standard k-ε model worsens in estimating the particle deposition in the 180° bend. In this
case, the standard k-ε model and the RSM estimate the maximum discrepancy of 13% at St=0.1
and 8% at St=0.001, respectively. The error averaged over the range of St is 9.5% and 6% for the
standard k-ε model and the RSM, respectively.

37

(a)

(b)

Figure 2.10: Comparison of the models (makers) and the empirical models (lines) shows that the
k-ε model over-estimates particle deposition when particles experience a longer residence time in
(b) the 180° bend (τ=4.4ms) than in (a) the 90° bend (τ=2.2ms)
38

The results from the RSM closure appear promising and more accurate than models with the
Boussinesq eddy-viscosity approximation in simulating flow encountering strong streamline
curvatures. This is attributed to the ability to capture the anisotropy of the turbulence fluctuating
motion and thus lead to improved results of particle deposition. The results discussed below are
thus all from simulations performed using an RSM with EWT due to the accuracy of this
combination.

Figure 2.11: Comparison of grade efficiency shows that the particle cutsize is 30% smaller in the
180°bend due to the longer particle residence time.

In Fig. 2.11 computed grade efficiency curves are shown for curved pipes for a changing
particle residence time. The flow Re number is 10,000 in the 90°bend and the 180°bend with

39

the curvature ratio equal to 5.6. The time needed for the flow to the 90°bend is 2.2 ms and it is
twice for the 180°bend. The grade efficiency for the U-bend is obviously larger than that for the
90°bend for the whole range of St, which indicates that particles experiencing a longer time
inside a bend have larger chance to deposit on a pipe wall and thus higher grade efficiency can
be obtained. The cutsize of a separation equipment is the size associated with a 50% probability
that it will leave the equipment (the bend) with the fluid. The particle cutsize for the 180°
-bend is
3/5 of that for the 90°bend.
The grade efficiency curves for 90-degree bends with different curvature ratios are calculated
and plotted against the Stokes number in Fig. 2.12. The bend curvature ratio is increased from
1.5 to 20. The pipe diameter and the flow Reynolds number are maintained the same as 8.51mm
and 30,000, respectively. Results show that the grade efficiency increases with an increase of the
bend curvature ratio. The empirical model of Pui et al. (1987) appears to be appropriate to a bend
with certain curvature ratio and should account for, or be modified to include, the bend curvature
ratio. For a 90-degree bend, an increase of the bend curvature ratio actually increases the particle
residence time in a bend, particles thus experience a centrifugal force for a longer duration.
The trend in Fig. 2.12 is likely due to the combined effects of the curvature ratio and the
particle residence time. In order to investigate the effect of the bend curvature ratio, the particle
residence time was fixed by using the same bend length in the simulations. The bend angles are
thus different due to different curvature ratios used. Bends with curvature ratios of 3, 5 and 9
corresponding to different bend angles of 180° 108°and 60°are studied to match the total travel
,
time. A comparison of the computed grade efficiency with Brockmann’s model (Eq. (2.11)) is
shown in Fig. 2.13. The plot shows that the increasing curvature ratio reduces the grade
efficiency of a curved pipe given the same total travel time. The centrifugal force is inversely

40

proportional to the bend radius, thus the increasing grade efficiency is due to a larger centrifugal
force. This trend is consistent with the observation of the experiment done by McFarland et al.
(1997). In Fig. 2.13, Brockmann’s model appears to provide accurate results except for cases
where the bend with curvature ratio equal to 9, which probably indicates that the model should
also account for the effect of a bend curvature ratio.

Figure 2.12: Comparison of bends with a varying curvature ratio for a fixed flow Reynolds
number in a 90°bend shows the grade efficiency is increased with the increase of the curvature
ratio.

41

Figure 2.13: Comparison of results for curved pipes with different curvature ratio and fixed
particle residence time shows the grade efficiency is increase with the decrease of bend curvature
ratio.

Fig. 2.14a shows the comparison between the different numerical results and the empirical
models from literature. Since the effect of δ is not included in the models of Pui et al. and
Brockman, the models are not accurate in estimating η in a bend with a changing δ. Although
MarFarland et al.’s model includes the effect of δ, a comparison with the grade efficiency curves

42

(a)

(b)

Figure 2.14: (a) Comparison shows that the models are not appropriate for varying δ in a 90°
bend (b) Comparison shows the new model match the results relatively well for δ between 3 and
8 but has reduced accuracy at extreme curvature ratios.
43

estimated from the models of Pui et al. and Brockmann are very different and it does not match
the numerical results for the cases studied. Based on the numerical results, a new model is
proposed to account for the effect of the St number, the bend angle, and the bend curvature ratio
on the grade efficiency. To fit the numerical results, the model is given as







  1  exp  0.528St

21/ 0.5    100%



(2.15)

This model is derived by curve fitting, without accounting for the numerous physical processes
encountered over the range of parameters considered. The comparison is shown in Fig. 2.14b.
The curves show that the new model fits the numerical results well for the cases when the
curvature ratio is from 3 to 8. For δ=3, the curve obtains a maximum discrepancy of 5% at
St=0.5. The maximum discrepancies for δ=5.6 and δ=8 are 7% and 7.5% at small Stokes number
St=0.1 and St=0.05, respectively. According to the previous study in this work, some errors are
caused by the overestimation of the grade efficiency of the numerical model at a small St number.
This indicates that the maximum deviation at small St is lower than the plot shows. Therefore,
the new model is actually very accurate for the cases discussed above. Extraordinary cases like
δ=1.5 and δ=20 are not fitted by this model. For δ=1.5, a recirculation flow is generated in the
tight bend. For δ=20, an extremely long bend may need to consider other factors for particle
deposition such as entrainment. Therefore, those cases should be excluded when using this
model.
Fig. 2.15 shows the grade efficiency of the 90-degree bend under different flow Re numbers.
The pipe diameter and the bend curvature ratio are fixed at 8.51mm and 5.6, respectively. The Re
number is increased from 10,000 to 100,000 and the results are plotted against St. The
comparison shows that the grade efficiency has a deviation of 9% when the Re number increases
by a factor of 10, which indicates that the effect of flow Re number is insignificant. This
44

conclusion is consistent with the statement of Pui et al. (1987). The effect of changing Re
number is actually accounted by using the St number. Since the particle size dp is inversely
proportional to the mean axial flow velocity U0, higher Re number will give a smaller particle
cutsize.

Figure 2.15: Comparison of grade efficiency shows a deviation of 9% when the Re number
increases by a factor of 10.

45

Since one is often interested in studying particle deposition in bends of various dimensions,
the effect of changing the diameter is shown in Fig. 2.16. The Reynolds number is 30,000 for all
the simulated cases and the bend curvature ratio is held at 5.6 in a 90°bend. The plot shows that
the particle cutsize is significantly affected by the pipe dimension. The results indicate that the
grade efficiency of the particles is decreased with the increase of the pipe size.

Figure 2.16: Comparison results show the particle cutsize is increased as the increase of the pipe
dimension.

46

Particle deposition patterns in the 90°bend with d=8.51mm and δ=5.6 are plotted with an
increasing value of the St number in Fig. 2.17. The flow Reynolds number for this case is 10,000.
All results are computed using EWT with a mesh that satisfies the required conditions. When
values of the St number are very small and equal to 0.001, approximately 3% of particles were
captured by the pipe wall. Deposition of small-size droplet due to turbulent fluctuations can
account for this small amount of deposition. From the definition of Stokes number (see Eq. (2.9)),
it can be seen that low values of the St number indicate a fast particle response time, and thus
correspond to particles that will follow closely the continuous flow. In this case, most particles
pass through the bend and only a small number of particles are trapped on the pipe wall because
of the flow fluctuating velocity. When the particle size increases (higher St number), particles
with high inertia tend to cross streamlines and deposit on the outer part of the bends because of
the centrifugal force (see Fig. 2.17e-g). Similar deposition patterns can be observed in a 180degree bend.
To decrease the cutsize and increase the pressure drop, either decrease the curvature ratio,
increase the residence time, or decrease the duct diameter, holding the other two parameters
constant. If pressure drop is not a concern, high particle deposition efficiency can be expected for
situations where a small pipe diameter d, a small bend curvature ratio τ, as well as a long
residence time τ are adopted. To accurately estimate the particle deposition, a turbulence closure
appropriate for flows with streamline curvatures should be used. In this work, the RSM is
employed with the linear pressure-strain correlation. The selection of an EWT as the near-wall
treatment is crucial in correctly estimating particle deposition along with constructing a highquality mesh.

47

Figure 2.17: Estimation of deposition patterns of particles changes with Stokes number.

2.7 Summary
A dilute turbulent flow passing through circular curved pipes and particle deposition are
studied in this work using computational fluid dynamic simulations. The simulations are

48

performed by solving the RANS equation and a force balance equation on the particles. Since the
accuracy of such numerical simulations is often in question when compared with experiments,
numerous case studies were performed to identify the potential source of the discrepancies.
Empirical models were used to validate the computational results.
Results of these computations show that turbulent bend flows possess, as expected,
complicated patterns influenced by the flow Re number and the bend configuration. Computed
results for the pressure drop through a bend are obtained from different near-wall treatments
based on the k-ε model and the RSM; all are close to experimental measurements provided that
the mesh meets the y+ requirement for the wall treatment selected except the local pressure at the
outer and inner pipe of the U-bend.
The grade efficiency for particle deposition in curved pipes is influenced by the bend diameter
d, the curvature ratio δ and the particle residence time inside the bend τ for a fixed flow rate. To
decrease the cutsize and increase the pressure drop, either decrease the curvature ratio, increase
the residence time, or decrease the duct diameter, holding the other two parameters constant. The
selection of the turbulence model and the wall treatment were found to be crucial in these studies.
Using RSM with EWT provides a much improved accuracy over the standard k-ε model in
estimating particle deposition patterns. The inaccuracy of SWF and NEWF for particle
deposition, when compared to EWT, indicates that these near-wall treatments are not appropriate
for modeling particle deposition on curved surfaces; this is likely due to the importance of
modeling the near wall fluctuations. In addition, a new model is developed to estimate the grade
efficiency of particulate flow in bends which include the effect of the curvature ratio. The given
model is an improvement in the empirical models of Pui et al. (1987) and Brockmann (1993).

49

The study of this work also provides the range of application of the empirical models discussed.
Broad guidelines for designing a bend with a high deposition efficiency are finally provided.

50

CHAPTER 3
SIMULATION OF DILUTE TURBULENT MIST FLOW IN CURVED PIPES USING AN
EULERIAN-EULERIAN METHOD

3.1 Introduction
In parts of the study presented in in Chapter 2, it is found that bends are associated with a
large pressure drop for single-phase flows. Particle trajectories obtained from the Lagrangian
particle tracking method also indicate that phase separation may occur when particulate-laden
flows pass through curved pipes. For gas/liquid mist flows with a large mass loading, Hoang and
Davis (1984) observed that the pressure drop associated with particle-laden flows is significantly
greater than that of a single-phase flow. Crane (1957) summarized that the total pressure drop for
an incompressible, adiabatic flow should include the contribution of wall friction, flow
detachment, momentum exchange between phases and downstream flow velocity profile
recovery. When mist flows pass through curved pipes, droplets may impinge on the pipe wall
where they accumulate and form a wall film. In general, the structure of this film is stratified
along the outer bend wall (Hoang and Davis 1984). However, Banerjee et al. (1967) observed
that within certain ranges of liquid flow rate and gas flow rate this film can travel inversely from
the outer band wall to the inner bend wall. This phenomenon is referred to as “film inversion”.
The consideration of bends is important in designing multiphase flows equipment and pipes. The
liquid film may be mitigated by varying the pipe geometry to prevent pipe erosion and reduce
pressure drop or, alternatively, promoted to improve the performance of a downstream separator.
The liquid film formation for gas/liquid flows passing through curved pipes has been vastly
studied in the literature, mostly from an experimental standpoint [5, 27, 34, 36, 38, 66, 90].
Computational fluid dynamics (CFD) provides a means to study such complex flows. In Chapter

51

2, the one-way coupling DPM has been proven accurate and efficient for flows under the
assumption of one-way coupling. To study the impact of the disperse phase and liquid film
formation, the flow has to be considered two-way coupling. Two-way coupling mist flows can be
modeled using the two-way coupling DPM where particle impact is treated as a momentum
source in the transient momentum equation. However, using the two-way coupling DPM is very
limited in film thickness for liquid film modeling (ANSYS FLUENT 12.1, 2009). Another
drawback is that the method is computationally expensive in tracking a large amount since
particles for fully coupled converged solutions are difficult to achieve. An alternative approach
of modeling a two-way coupling flow is to use the drift flux model based on the mixture theory
proposed by Manninen et al. (1996). The drift flux model has been broadly used in simulating
aerosol deposition in human airways and fine particles (typically smaller than 1µm) in laminar
flows were usually studied in this field [62, 88, 110]. Hossain and Naser (2004) studied a
turbulent bend flow with larger particles (2 µm - 20 µm) and showed that the method is
promising in predicting the concentration of a gas-solid granular flow in a pipe loop with four 90°
bends. However, the application of the drift flux model in predicting pressure drop and liquid
film formation in bends cannot be found in literature.
In this Chapter, a study is performed for evaluating the performance of the drift flux model in
predicting pressure drop and film formation of a turbulent mist flow when passing through a
curved pipe. In addition, the numerical method is used as a tool for bend design related to
pressure drop and deposition efficiency. The Realizable k-ε model was used to close the RANS
equation. It would have been preferable to close the RANS equation with the Reynolds Stress
Model but converged solutions for most multiphase flow problems studied here could not be
obtained. The flow at the near-wall region is resolved by the enhanced wall treatment (EWT).

52

Numerical results allow verifying if the model is able to capture phenomena such as secondary
flows and film inversion. Furthermore, the accuracy of the model was validated by comparing
the numerical results of pressure drop with the empirical models from literature. Subsequently,
the influence of parameters such droplet size dp and the volume fraction of the water phase αp on
the flow patterns is studied, as well as on pressure drop, secondary flow intensity and the liquid
film thickness. Finally, a bend allows for a low pressure drop and high droplet deposition
efficiency is identified for a given flow condition.

3.2 Multiphase flow modeling
Multiphase flows can be categorized as gas-liquid, gas-solid, solid-liquid and liquid-liquid
flows. Widely different processes are encountered in each flow and this makes it nearly
impossible to find a universal model appropriate for simulating all possible situations. The
selection of appropriate closures for a given theory will thus vary according to the nature of the
flow studied. Eulerian-Eulerian methods found in the literature include the Volume of Fluid
(VOF) model, the drift flux model, and the two-fluid model. The application of a VOF model
requires significant computer resources to track in detail the movement of a liquid/gas interface.
The two-fluid model can tackle a wide variety of flows but also requires significant computing
effort due to the number of differential equations to solve. The drift flux model in the mixture
theory appears to be efficient and can be accurate, it requires a smaller set of governing equation
and provide the mixture properties of the multiphase flow studied [Manninen et al., 1996]; it is a
statistical approach that does not provide the same information as a two-equation model. The
mixture theory with a drift flux model of [Manninen et al., 1996] was selected in this work since
it involves design and numerous simulations.

53

3.2.1 Multiphase turbulent flow modeling
The mixture theory is a set of governing equations that includes a balance of mass and a
balance of momentum. A second continuity equation is needed to compute the volume fraction
of the secondary phase. These equations result from phase averaging, followed by ensemble
averaging. Equations similar to the Reynolds-Averaged Navier Stokes (RANS) equation for
single phase flows are obtained. The equivalent of the Reynolds stress found in the momentum
equation is closed using the realizable k-ε turbulent model; the coefficients used in the closure
are those found for single phase flows. The averaged continuity equation for the air-water
mixture is given as
  (mix umix )  0

(3.1)

with the mean mixture density

mix  f f  pp

(3.2)

and the mean mixture velocity

umix  (f f uf  ppup ) / mix

(3.3)

where ρ is the density; α is the volume fraction; u is the instant velocity; the subscript “f” and “p”
denote the carrier phase and the disperse phase. In this work, the carrier phase is air and the
disperse phase is water droplets.
The drift flux model allows different phases to be interpenetrable with a volume fraction
characterizing each phase. The summation of the volume fraction of all phases satisfies
n

 k  1

(3.4)

k 1

For a gas/liquid two-phase mist flow, the liquid is treated as the disperse phase. The distribution
of the liquid phase is determined by the value of the liquid volume fraction αp in the flow domain.
54

The media is liquid for αp=1 and gas for αp=0. It is a mixture for the liquid volume fraction
between 0 and 1. The liquid volume fraction is obtained from a transport equation derived from
mass balance and is given as

 (p u mix )   (p udr,p )

(3.5)

where udr,p is called the drift velocity for the liquid phase which needs to be closed through
models. It defined as the velocity difference between the droplet velocity and the mixture
velocity, given as
udr,p  up  umix

(3.6)

To model the above drift velocity, one relates it by convention to the phase slip velocity,
obtained
udr,p  uf ,p (1  G p )

(3.7)

The slip velocity uf,p is a difference between the gas phase and the liquid phase defined as
uf ,p  up  uf

(3.8)

Gp is the mass fraction of the water phase defined as

G p  pp / mix

(3.9)

The momentum equation as presented in Manninen et al. (1996) for the mixture is given as

  (mix u mix u mix )  P  mix g    ( τ m  τ Tm )
   [mix G p (1  G p )uf ,puf ,p ]

(3.10)

where τm and τTm are the viscous shear stress and the turbulent shear stress which are closed by
using the following constitutive equations according to Johansen et al. (1990):
τ m  mix (umix  uT )
mix

(3.11)
55

T
τ Tm  mix T (uf  uf )  p Dp (uf p  uf ,pp )

(3.12)

where the mixture dynamic viscosity is defined as

mix  f f  pp

(3.13)

According to Picart et al (1986), the diffusion coefficient for the liquid phase Dp is given as:
2

uf ,p 

Dp  T 1  0.85

2k / 3 



1/2

(3.14)

where νT is the turbulent viscosity for the gas phase, defined as

k2
T  C


(3.15)

An essential prerequisite of using the drift flux model is to close the slip velocity in Eq. (3.7),
(3.10), (3.12), and (3.14). Ishii (1975) developed the algebraic slip model based on an
assumption of a local equilibrium. The constitutive equation for the slip velocity is given as
uf ,p 

p (p  mix )
fdrag

p

(g  umix u mix )

(3.16)

The droplet relaxation time τp is given in Eq. (2.5). The disperse liquid phase is treated as
solid spherical droplets with a single diameter dp. Drop breakup, drop coalescence, and droplet
interaction are ignored. Under this assumption, the drag function fdrag is proportional to the
Reynolds number of the mixture Remix obtained from the empirical measurement of Schiller and
Naumann (1935) for a turbulent flow passing a spherical particle, given as

fdrag  0.0183Remix

(3.17)

56

A simplified version of the mixture model can be achieved for a homogenous flow which
assumes that the velocity difference between the phases u f ,p equals to zero. This kind of flow is
out of the scope in this work since the slip velocity is significant due to the gravity and the
centrifugal force induced by a bend.

3.2.2 Turbulence closure
To complete the formulation of the turbulent mixture momentum equations, transport
equations need to be solved each for the turbulent kinetic energy k and the dissipation rate ε,
respectively. The realizable k-ε two-equation model [Shih et al., 1995] is used in this work based
on the Boussinesq viscosity assumption [Boussinesq, 1877]. Any model based on an eddyviscosity approach has shortcomings for flows with strong streamline curvatures but this option
the only available choice that could converge easily with EWT. Using the approach of
Pourahmadi (1982), the transport equations are extended to a two-phase flow. The equations are
presented in Appendix A for the sake of brevity.

3.2.3 Near-wall treatment
From Chapter 2, the EWT was proven to be the most accurate near-wall treatment method in
dealing with particle deposition in curved pipes. Since studies will be focused on droplet
deposition and liquid film formation on pipe wall of a mist air-water flow, the EWT is used in
this work for the near-wall treatment.

57

3.3 Numerical methodologies
The governing equations are solved in a steady and adiabatic state in the commercial code
ANSYS FLUENT 12 using a SIMPLE algorithm (Patankar and Spalding 1972). The transport
equations for k and ε are spatially discretized using a second-order upwind scheme. A correction
equation derived from the continuity equation and the momentum equation for pressure is
discretized using the same PRESTO! scheme as we did in the previous chapter. The Quick
Scheme (quadratic upstream interpolation convective kinetics), a higher-order unwinding method
proposed by Perng and Street (1989), is employed for computing the volume fraction of the
droplet phase αp. The stability of the numerical computations is improved by reducing the underrelaxation factor to 0.2 for the pressure, 0.5 for the momentum and 0.001 for the slip velocity.
-4

The residual criterion for convergence (the “monitors”) is set to 10 for all variables. Meanwhile,
three monitors including the difference of the mass flow rates at the inlet and the outlet of a bend,
the area-average weighted pressure and the flow rate of the secondary phase at the bend exit are
used to ensure a converged solution.

3.4 Geometry and meshing
The geometric configuration and the mesh construction can be referred to Fig. 2.1 and Fig. 2.2
in Chapter 2. Compared to the one-way DPM, the drift flux model is more computationally
expensive for more equations need to be solved. To reduce the computational effort, a fully
developed velocity profile as well as the turbulent intensity profile based on the flow Reynolds
number is applied to the entrance of the pipe as the boundary conditions in all simulations. At the
pipe wall, the no-slip boundary condition is applied. The boundary condition for the outflow is

58

set as pressure outlet with the static pressure equal to zero. The applied boundary conditions are
summarized in Table 3.1. No heat transfer is considered in the process.

3.5 Results and discussion
The work consists of two parts. Part one is to assess the performance of the numerical
approach. In part two, a parametric study is conducted to investigate the impacts of the droplet
size and the volume fraction. The model is used for a bend geometry design based on the
numerical results of droplet deposition and pressure drop by changing the curvature ratios of a
bend.

Table 3.1: Boundary conditions used in all simulations

Outlet

velocity inlet [m/s]
turbulent intensity [%]
turbulent viscosity ratio
pressure outlet [Pa]

Wall

no-slip boundary [m/s]

Inlet

Fully developed profiles
Fully developed profiles
Fully developed profiles
Pstatic=0
uwall=0

3.5.1 Model assessment
Empirical studies show that air-water flows passing a bend not only create complex flow
patterns, such as secondary flow and film inversion but also cause a large pressure drop. The
mixture model was checked qualitatively to see whether or not the model is capable of capturing
the complex flow patterns and film formation. Moreover, to investigate its accuracy, the model
was tested quantitatively by comparing the numerical results of the pressure loss with the
empirical models from literature.

59

(a) Flow patterns and liquid film modeling
Fig. 3.1 shows the numerical results of an air-water flow passing through a vertical 180˚ bend
(U-bend) with the diameter equal to 0.0508m and the curvature ratio δ equal to 5. The flow is in
5

the turbulent region with Re=1.25× and with the volume fraction for the water phase α=0.06.
10
The droplet size for this case is 30μm. The liquid distribution (expressed by the volume fraction
for the water phase) and the streamline of the continuous phase are shown at the middle plane.
Furthermore, the liquid distribution and the streamlines for the gas and the liquid phase (“G” and
“L” for short) are shown along the bend at the cross sections of 0° 45° 90° 135° 180°
,
,
,
,
deflections and the location at downstream 2d from the bend exit.
The plots show that the streamlines for both phases exhibit very similar patterns at the cross
section of 0°deflection and the location at downstream 2d which is reasonable since at those
locations flows are only under the gravity field. Around the bend, the results show that the air
phase and the water phase have different behavior. Due to the centrifugal force, the air
streamline and the water streamline present a secondary flow pattern at the cross section of 45°
deflection. Meanwhile, the streamlines ended at the wall illustrate that the disperse phase tends to
deposit and start to accumulate to form a liquid film on the outer bend wall. Further downstream
in the bend, while simple streamlines of the water phase show that the droplets are traveling from
the inner bend to the outer bend due to the centrifugal force, the streamline for the air phase
exhibit more complex patterns due to the interaction between the phases. With respect to the
liquid film, it can be seen from the diagram that the structure presents a stratified feature along
the outer bend. This statement is consistent with the observation of Hoang and Davis (1984).

60

αp

Figure 3.1: The air and water streamlines show the development of the flow patterns and the
interaction between two phases. The contour of αp shows the liquid film formation along the
180˚ bend. “G” and “L” denotes gas and liquid, respective.

61

(b) Film inversion
In most situations of a horizontal bend, a liquid film forms on the pipe wall and travels along
the outer bends due to the centrifugal force. Under certain flow condition, however, liquid film
can travel along the inner bends regardless of existence of the centrifugal force. This
extraordinary phenomenon is referred as “film inversion” and was observed in a helical pipe by
Banerjee et al. in 1967. They stated that film inversion only occurs at low liquid but high gas
flow rates in a horizontal bend. Hart et al. (1988) derived a criterion for this phenomenon based
on the static equilibrium of moment, given as
mf / mp  195(f / p )1/2 (f / p )3/4

(3.19)

where mf , mp are the air and the liquid mass flow rate; f , p denote the gas and liquid density;

f , p indicate the air and the liquid dynamic viscosity. The ratio of the air mass flow rate and
the liquid mass flow rate is the function of their density and the volume fraction, obtained
 (1   p )
 f
mp
f  p
mf

(3.20)

Eq. (3.19) and Eq. (3.20) are solved to obtain the following criterion for an air-water flow

p  0.0037

(3.21)

To verify the models, the case is run for a mist flow with αp=0.003 passing through a horizontal
U-bend with d=0.0508m and δ=30. The gravity is in a negative z direction. The flow is in the
turbulent region with Re=16800. The droplet size uses the Sauter mean diameter d32=149μm.
This is solved by Eq. (3.22) (Azzopardi 1985)

G
d32
15.4

 3.5 LE
0.58

p U f
We

(3.22)

62

2
where    / pg is a length scale; We  p Uf  /  is a Weber number; GLE = ρpUp is the

entrained liquid mass flux ; the surface tension σ for air and water is equal to 0.072; the inlet air
velocity Uf and the inlet water velocity Up are 25m/s and 0.075m/s.

αp

Figure 3.2: Water film forms at the inner bend wall affected by the secondary flow pattern.

Fig. 3.2 shows the numerical result of the film location at the outlet cross-section of the Ubend. The gravity is in a minus z direction. It can be seen that the film locates at -23°from the
vertical centerline which indicates film inversion occurs and the model is able to capture the
phenomenon. The liquid film is thin due to the low water flow rate. The secondary flow patterns
are not symmetric to the horizontal axis due to the gravity. The film inversion is enhanced by this
secondary flow which transfers the liquid film from the outer wall to the inner wall. In addition,
a small vortex is generated due to the liquid film. From Fig. 3.3, more detail is covered for the
volume fraction distribution of the water phase at the pipe wall for the same cross section. The

63

volume fractions are almost the same around the outer bend except -270°where it is as low as
0.05. This may be caused by the secondary flow of the air phase due to the centrifugal force. In
addition, liquid film is observed at the inner bend between 0°and -30°where the volume fraction
is equal to unit which indicates that film inversion occurs.

Figure 3.3: Volume fraction distribution of the water phase shows the location of liquid film at
the pipe wall of the outlet cross section.

(c) Pressure drop prediction
Much work has been done on measuring the pressure loss of turbulent gas/liquid flow in
horizontal 90° bends since 1960s (Fitzsimmons 1964, Sookprasong 1980, Norstebo 1986,
Mandal and Das 2001, etc.). Based on the experimental data, different empirical models were

64

developed to predict the pressure drop of two-phase flows. Lockhart and Martinelli (1949)
computed the pressure loss of a gas/liquid flow by using a multiplier concept. This method has
been wildly used in developing empirical models by many authors (Chisholm 1971,
Sookprasoog 1980 and Paliwoda 1991). An alternative proposed by Mandal and Das (2001) is
the friction factor method based on the Fanning friction factor equation. A detail review on
empirical models of pressure drop was given by Azzi et al. (2000).

Figure 3.4: Comparison shows good agreement of the numerical results with the Paliwoda’s
empirical model

Simulations are run for a mist flow passing through a horizontal 90°bend with d=0.0508 and
δ=5. The air mass quality is from 0.056 to 0.16. Numerical results of the pressure drop between

65

the entrance and the exit of the bend are compared to different empirical models from literature
and shown in Fig. 3.4. The plot shows that the numerical results agree with Paliwoda’s empirical
model very well. It can be also seen that Chisholm’s B-type predicts a very high pressure drop
compared to other empirical models. This deviation is consistent with Norstebo’s observation.
He obtained a 110%+ discrepancy when he compared his experimental data with the Chisholm’s
method which indicates that Chisholm’s empirical model may overpredict the pressure drop.
Another observation of the plot is that the deviation of the pressure drop predicted by other
empirical models is getting greater when the air flow mass quality is increased. Since deviation
of the pressure drop predicted by different empirical models is large, it is uncertain whether the
mixture model is accurate. But it can be concluded that the model predicts results in a reasonable
range.

3.5.2 Applications
The test cases for parametric studies are summarized in Table 3.2. Case 1 and case 2 are to
investigate the effect of the droplet sizes and the volume fraction of the water phase on the flow
patterns, droplet deposition and pressure drop. A 90°bend with the pipe diameter d=0.0508m
and with the curvature ratio δ=5 is used. In addition, a single-phase air flow with the same Re

Table 3.2: Summary of the parameters used in the test cases
Test cases
Droplet size, dp [μm]
Volume fraction, αp[%]
5
Re (× )
10
Bend angle, θ [°
]
Incline angle, β [°]
Curvature ratio, δ

Case 1
1~50
6
1.25
90
90
5

Case 2
30
3~15
1.1~2.3
90
90
5
66

Case 3
30
6
1.25
180
0~90
5

Case 4
30
6
1.25
90,180
0
1.5~9.0

number is modeled in the same bend geometry. The results are compared those from air-water
flows. In Case 3, the effect of the incline angle on the pressure as well as the deposition
efficiency is investigated. In case 4, the mixture model is used for optimizing the bend geometry
for low pressure drop as well as high droplet deposition efficiency in a horizontal bend by
changing the curvature ratio.

(a) Parametric studies
Fig. 3.5 shows the comparison between the air flow and the air-water flows for the
streamlines of the continuous phase at the exit of a vertical 90° bend. The droplet sizes are
changed from 1μm to 50μm. The volume fraction equal to 6% is used for all droplet sizes. For
droplets as small as 1μm (Fig. 3.5c), the secondary flow streamline looks similar to the singlephase flow (Fig. 3.5b), which indicates that droplets with a small size like this have a marginally
effect on the continuous flow. However, with the increase of the droplet size, the effect of the
water phase becomes significant. Fig. 3.5e shows that the secondary flow of the continuous
phase is reversed by the disperse phase and rotates clockwise. With a larger droplet size, more
complex patterns like multi-pair secondary flows occur (Fig. 3.5f).

67

Figure 3.5: Comparison between the single-phase air flow (b) and the air-water flows (c-f) on the
secondary flow patterns at the cross section of 90°deflection (a) changing with different droplet
size

The different effects of the droplet size are attributed to the droplet behavior and can be
studied using the St number. The definition of this dimensionless number and how it affects the
particle behavior have been discussed in Chapter 2. Fig. 3.6 shows the volume fraction
distribution of the water phase at the bend exit. Here we assume that a liquid film is formed
when αp≥0.8. It can be seen from Fig. 3.6a that for a small St number the droplets follow the air
flow closely and no separation between the phases occurs. A small St number means the disperse
phase responds to the continuous phase quickly which results in a homogenous type of flow. For
St=0.027, droplets start moving from the inner bend to the outer bend which indicates separation
between two phases occurs (Fig. 3.6b). With a larger St number, droplets deposit on the wall and
accumulate to form a liquid film. This film becomes thicker with more droplets deposited on the
wall (Fig. 3.6c-f). Finally, droplets with large enough St number will all deposit on the wall to
68

thicken the liquid film (Fig. 3.6e). From the above discussion, it is found that the conclusion
made for a dilute flow on the effect of the St number on the particle behavior is still valid for
dilute flows with a large mass loading.
Plots from Fig. 3.7 show that the continuous secondary flow patterns vary with the increase of
αp at the bend exit. A single size equal to 30μm is used for the droplets. It can be seen that the
streamlines of the secondary flow patterns at the cross section are resemblant, which indicates
that the volume fraction of water phase has marginally effect on the flow patterns and the flow
intensity.
As regards the pressure drop, Fig. 3.8 shows the comparison of the pressure drop caused by
the bend component between the air flow and the air/liquid flow. The pressure drop in the air
flow is labeled with a circle. One can observe that the pressure drop in the air-water flow is much
greater than that in the air flow. This observation is consistent with Hoang and Davis’ statement
(1984). From Fig. 3.8a, it is shown that the pressure drop is increased with the increase of the
droplet size. It then decreases for droplets larger than 12μm because more droplets escape from
the core flow and deposit on the bend wall. For droplets greater than 40μm, the pressure drop
tends to be stable since those cases have the same situation where almost all the droplets deposit
on the pipe wall. In Fig. 3.8b, nevertheless, it can be seen that the pressure drop is simply
proportional to the volume fraction of the water phase.

69

Figure 3.6: Volume fraction distribution shows the liquid film location at the cross section of 90°
deflection changing with different droplet size

Figure 3.7: Secondary flow patterns and flow intensity at the cross section of 90°deflection for
flow with a different αp

70

Figure 3.8: Pressure drop of the bend section changing with (a) the droplet size, (b) volume
fraction of the water phase

71

Figure 3.9: Droplet deposition efficiency changing with (a) the droplet size, (b) volume fraction
of the water phase
72

The droplet deposition efficiency is defined as the ratio of the mass flow rate for the liquid
film to the total mass flow rate for the liquid at the bend exit. Fig. 3.9 shows the tendency of the
droplet deposition efficiency changes with the increase of the droplet size and the volume
fraction of the water phase, respectively. It can be seen from Fig. 3.9a that no deposition
efficiency is obtained for the droplet size less than 5μm due to a small St number. For droplets
larger than 5μm, the droplet deposition efficiency is proportional to the droplet size. The
deposition efficiency reaches 0.92 and becomes stable when droplets are larger than 40μm. This
can be treated as the maximal deposition efficiency that this bend geometry can obtain. With
regards to the influence of αp, Fig. 3.9b shows that αp only has a marginal effect on deposition
efficiency.

(b) Incline angles
The effect of the incline angle on the pressure drop and the deposition efficiency is
investigated. Fig. 3.10 demonstrates a U-bend changing its incline angle β from 0° (horizontal) to
90°(vertical) with an increment 15° Comparison of the pressure drop and droplet deposition
.
efficiency is shown in Table 3.3. Result shows that the bend incline angle does not affect the
pressure drop and the droplet deposition efficiency.

Table 3.3: Pressure drop and deposition efficiency in a U-bend with different incline angles
Angle β[°]
η [%]
∆P [Pa]

0
89.6
1196

15
89.8
1194

30
90.0
1194

73

45
90.1
1196

60
90.2
1200

75
90.1
1205

90
90.2
1200

Figure 3.10: Schematic of a U-bend changing from horizontal position to vertical position with
an increment of 15o. For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this dissertation.

(c) Bend optimization
Since the pressure drop is simply proportional to αp and the droplet deposition efficiency is
not affected by αp. A moderate droplet size with a single volume fraction can be chosen for a
bend geometry optimization based on the pressure drop and the deposition efficiency. The word
“moderate” here means the droplets have to be large enough to allow the occurrence of
separation between the phases but not too large so that not all the droplets deposit on the bend
wall regardless of the geometry of bends. In industrial transportation, a large range of droplet
sizes would exist in gas/liquid mist flow. The idea is that if we can obtain a geometric
configuration to improve the deposition efficiency by using a moderate size, the configuration

74

could be employed to improve the deposition efficiency for practical applications. In this work,
30μm is chosen as the droplet size with the volume fraction equal to 6%.
Fig. 3.11 shows the bend pressure drop and the droplet deposition efficiency changed due to
the varying curvature ratio of a 90°bend and a 180°bend with the same pipe diameter equal to
0.0508m at the horizontal plane. Results show that the tendency looks similar for both bends. It
can be seen from the plots that a tight bend causes extremely high pressure drop and low droplet
deposition efficiency. This extraordinarily high pressure drop is caused partially due to flow
detachment. Secondly, it is because the flow is “chocked” by a recirculation flow formed after
the tight bend due to the extremely sharp curvature streamline. On the other hand, the low
deposition efficiency is due to the short experience of centrifugal force and film entrainment. At
the same curvature ratio, pressure drop and deposition efficiency in a 180°bend is higher than
those in a 90°bend. The pressure drop is decreased with the increase of the curvature ratio at the
beginning as a sharp curvature streamline can cause flow detachment and lead to a high pressure
drop. With a higher curvature ratio (e.g. δ ≥ 7), however, the pressure drop smoothly increases
with the increase of the curvature ratio. The effect of the residence time at this curvature ratio
becomes more important than that of the curvature streamline as more wall friction can be
expected due to a longer residence time for the flow to stay in a bend. On the other hand, the
droplet deposition efficiency is proportional to the curvature ratio at the beginning and becomes
flat or increases slowly after the curvature ratio is larger than 5 and 7 for the 90°bend and 180°

75

Figure 3.10: Tendency of bend pressure drop and droplet deposition efficiency in a horizontal 90°
bend and a 180°bend shows the bend geometry is optimal at curvature ratio equal to 5 and 7,
respectively.

76

bend, respectively. Therefore, it can be concluded that a bend with the curvature ratio equal to 5
and 7 is the optimal geometry respectively for the 90° bend and 180° bend to minimize the
pressure drop as well as obtain relatively high droplet deposition efficiency.

3.6 Summary
Numerical simulations are performed using the commercial code ANSYS FLUENT to model
turbulent mist air-water flows passing through circular 90°and 180°bends. The flow patterns,
pressure drop and droplet deposition efficiency of the bend component are studied in this work.
The simulations are performed using the mixture model closed by the modified Realizable k-ε
model and using the enhanced wall treatment to solve the flow at the near-wall region.
It is observed that a two-way coupling mist flow is much more complicated than a one-way
coupling flow. The drift flux model studied in this Chapter is proven to be capable of capturing
the complex flow patterns observed experimentally. For example, secondary flows of the
continuous phase are induced by the centrifugal force and impacted by the disperse phase. The
liquid film formed by droplet deposition is usually stratified along the outer bend but the film can
be inversed to the inner bend under certain conditions. In addition, the comparison between the
numerical results and the empirical models shows the pressure drop estimated by the model
agree well with Paliwoda’s work. Results also show that the pressure drop and the droplet
deposition efficiency are significantly affected by a bend. Based on these computations, bend
geometry design was performed. It is found that 90°and 180°bends with the curvature ratios
equal to 5 and 7 respectively can be used to achieve a minimal pressure drop as well as relatively
high deposition efficiency regardless of the bend angles.

77

CHAPTER 4
SIMULATION OF FLUID-STRUCTURE INTERACTION USING AN ELASTIC
FORCING METHOD BASED ON AN IMMERSED BOUNDARY METHOD

In this Chapter, a numerical technique called the immersed boundary (IB) method is discussed
and modified for the simulation of a solid structure immersed in a fluid flow. The purpose of this
technique is to tackle the fluid-structure interaction and use a simple Cartesian rectangular mesh
to solve the fluid flow.

4.1 Introduction
Fluid-structure interactions are commonplace in biological systems and industrial applications.
Examples include blood flow in hearts, red blood cells in arteries, aquatic animal locomotion,
and vehicles or aircrafts in turbulent air flows. The numerical simulation of a fluid interacting
with a structure is a very challenging problem. Consider a structure immersed in a fluid flow
shown as in Fig. 4.1a. Conventional approaches of Direct Numerical Simulation (DNS) of fluidstructure interaction use a body-fitted mesh. The grid is first constructed on the object surface

S and extended to the fluid flow Ωf. There is no grid located inside the structure. This method
allows to easily imposing the no-slip boundary condition on the surface of the immersed
structures. However, if the structure is complex in geometry, deformable, and movable in the
flow, an adaptive mesh is inevitable. This requires an intense computational effort and frequent
human input as the flow is matching with time.
To avoid using an adaptive mesh, an alternative approach called immersed boundary (IB)
method can be used. The method considers the impact of the immersed structure in a force
density term of the fluid equation and solves the equation for the entire domain using a Cartesian

78

rectangular mesh. In general, the structure surface does not coincide with the Cartesian mesh
(Fig. 4.1b). Nevertheless, communication of variables, such as velocity and force density
between the Cartesian mesh and the structure surface could be achieved by using an interpolation
or an extrapolation process. The focus at this moment is on how the force density term is
obtained. Numerical approaches used for rigid structure simulation in literature can be classified
into two categories: virtual boundary methods (Peskin 1972, Goldstein et al. 1993) and fictitious
domain methods (Mohd-Yusof 1997, Fadlun et al. 2000, Kim et al. 2001, and Uhlmann 2005). A
virtual boundary method uses a set of discrete points to mark the structure surface and evaluate
the force at those points. The points are treated as virtual points existing in the flow domain. The
points representing the structure surface can be tracked and moved. The virtual boundary method
was first proposed by Peskin in 1972 to study a blood flow interacting with heart valves. By
assuming that the structure surface is an elastic boundary, the Lagrangian points are linked with
springs. The springs are dampened or stretched by the Lagragian points moving with the
surrounding flow. The force exerted on the points is then a function of the spring deformation
and its stiffness coefficient. This method is also referred to as the elastic forcing method.
Goldstein et al. (1993) used the similar idea for a stationary rigid structure by applying large
stiffness coefficients. However, Lee (2003) stated that a large stiffness coefficient may cause
numerical instability if the time step used is larger than the characteristic time scales of the
spring oscillation. Instead of evaluating the force on the Lagrangian points with an elastic
boundary assumption, Mohd-Yusof (1997) obtained the force density term from the momentum
equations by imposing a desired velocity at the points wherever the surface is intersected with
the Eulerian mesh via a B-Spine function. This approach was referred to as a fictitious domain
method. Fadlun et al. (2000) applied the method to flow interaction with complex geometries.

79

The desired velocity on the boundary was obtained via volume fraction weighing or a linear
interpolation. The method was proved to be able to overcome the drawback of an elastic forcing
method on time step restriction. Kim et al. (2001) improved the approximation scheme for
stability and accuracy by using a second-order-accurate interpolation scheme. However, force
oscillations occur at the boundary when the above methods were used to move the structure
boundaries. Uhlmann (2003) stated that force oscillations are not appropriate for particulate flow
simulation. Later, Uhlmann (2005) improved the stability of the direct forcing method by
combining the Dirac delta function of Peskin (1972) to smooth the force on the boundary of a
rigid body.

Figure 4.1: Schematic configuration shows a Eulerian mesh for the fluid domain Ωf and a
Lagrangian mesh around the fluid-object interface S .

80

By treating the impact of the immersed boundary in a force density term of the fluid equations,
the IB method possesses great flexibility and is compatible with most of the numerical schemes.
Recently, it becomes popular to combine the technique with the lattice Boltzmann method. The
lattice Boltzmann method is proved to be simple, algebraic, and intrinsically parallelized (Wu
and Shu 2010). Fogelson and Peskin (1988) implemented the IB method based on a spectral
method for rigid particle simulation in a Stokes flow. Glowinski et al. (1999) combined the IB
method with a Finite Element Method for particle sedimentation. The IB method based on a
finite volume method was used by Sharma and Patankar (2000) for rigid particulate flows. In this
work, a finite different method is used. It is very common that the incompressible N-S equation
is used with the IB method in dealing with fluid-structure interaction (Lai and Peskin 2000,
Hö
fler and Schwarzer 2000, Fadlun et al. 2000, Uhlmann 2005). However, solving the
incompressible N-S equation usually ends up with finding a solution to a system of equations or
the Poisson equation. The elliptic behavior of the Poisson equation needs the equations to be
solved for the whole domain simultaneously through an iteration process. This requires a large
computational effort each time step. Although the overall computational effort could be reduced
by using a large time step, a large time step is not suitable for research on fluid-structure
interaction. Structure moving with a large time step would lead to information missing and thus
cause inaccuracy in predicting its trajectory. In this work, in order to improve the computational
efficiency, compressible N-S equation is solved by using an explicit MacCormack scheme
(MacCormack 1969). The MacCormack scheme with a second-order accuracy in space and time
has been widely applied to compressible flows as well as incompressible flows due to the
advantages of easy implementation and friendly parallelizing computation. However, literature

81

on using the IB method based on the MacCormack scheme for fluid-structure interaction
simulation cannot be found.
In this chapter, an algorithm is developed for combining the MacCormack scheme with the IB
method and used for fluid-structure interaction. The equations are discretized in time and space.
The performance of the numerical approach on simulating fixed and movable rigid structure is
investigated. The numerical approach is tested in the following 2-D cases (1) pressure driven
flow, (2) uniform flow passing a stationary cylinder, and (3) settling rigid body impinging on a
wall.

4.2 The MacCormack/IB method
To implement the IB method based on the MacCormack scheme, we consider the impact of
the immersed structure to its surrounding flow in a force density term f of the compressible N-S
equation. For a 2-D case, the equations for the MacCormack/IB method are given as

U E F


f
t x y

(4.1)

f (x, t)   Ff ( X, t)h (x  X(s, t))ds

(4.2)

X(s, t)
 U( X(s, t), t)
t

(4.3)

=  u( x, t)h ( x  X(s, t))dx
where U, E, and F are in a vector form, given as

 
G  u 
 
v 
 

(4.4)

82

u

 2

E  u  p  xx 
uv  

xy



(4.5)

v




F  uv  xy


v 2  p   yy 



(4.6)

f x 
f  
f y 

(4.7)

The shear stress in the above vectors are given by
2  u v 
 xx    2  
3  x y 
2  v u 
 yy    2  
3  y x 

(4.8)

 u v 
 xy   yx     
 y x 

The pressure p is related to the flow density through the equation of state
2

p = ρc

(4.9)

where c is the speed of sound in the media. For isothermal flow with a low Mach number
(Ma=||u||/c: ratio of the flow speed and the sound speed), the flow can be approximated to be
incompressible (Kundu and Cohen 2004). The explicit MacCormack scheme uses a predictor
step and a corrector step to solve Eq. (4.1) (MacCormack 1969). The equation is discretized in
time as follows
Predictor:
 En F n 
n
G  G  t 

  tf
y 
 x
*

n

(4.10)

83

Corrector:
G n 1 


 E* F* 
1 n
*
G  G*  t 


  tf 
 x
2
y 






(4.11)

*

In the predictor step, the variables with are preliminary results computed using the variables
at the current time step n. In the corrector step, the variables at time step n+1 are obtained from
the preliminary results. The detail of equation discretization is presented in Section 4.3.
*

One should notice that the force density term (f and f ) in the above equations should be
computed first before Eq. (4.10)-(4.11) can be solved. Eq. (4.2) converts the force density from
the Lagrangian points Ff (X,t) to the Eulerian mesh f (x,t) through a Dirac delta distribution δh(xX(s,t)). When the velocity field is updated by solving the fluid equation (Eq. (4.1)), the second
row of Eq. (4.3) is used to interpolate the velocity from the Eulerian mesh u(x, t) to the
Lagrangian points U(X(s,t)) using the same Dirac delta function. The Lagrangian points
representing the structure are then moved to a new location by using the first row of Eq. (4.3).
The approach used to compute the force density on the Lagrangian points F f and the Dirac delta
function δh used in this work are covered in the following sections.

4.2.1 Elastic forcing method
The elastic forcing method introduced by Peskin (1972) assumes that the structure surface
marked by Lagrangian points is an elastic boundary. This can be accomplished by connecting the
Lagrangian points with springs. Each spring possesses a stiffness coefficient and a length
decided by the material property and the relative location of two Lagrangian points, respectively.
The surrounding flow interacting with the structure could move the Lagrangian points, which
84

could lead to spring stretching or damping. Therefore, the force generated on the Lagrangian
points is a function of the spring deformation and the material elasticity. The method has been
proven to be well-suited in deformable structure applications. Examples are blood flows in a
beating heart (Peskin 1972, Kovacs et al. 2001, McQueen and Peskin 2001, Vigmond et al 2003),
blood cells transporting in arteries (Pacull 2006, Kim et al. 2009, Adib et al. 2010), and
parachute dynamics (Kim and Peskin 2006, Kim and Peskin 2009, Miyoshi et al. 2009). In this
work, the technique is employed to construct a rigid body based on a feed-back mechanism.
Depending on whether a rigid structure is stationary or movable, the springs are constructed
differently.

(a) Stationary rigid body
The elastic forcing method is quite handy in stationary rigid body modeling. According to Lai
and Peskin (2000), each Lagrangian point on the structure surface is linked to its equilibrium
location by a spring with a sufficiently large stiffness coefficient κ. All the springs obey the
Hooks’ law and the forcing on the kth Lagrangian Fk(t) is given as

Fk (t)  (Xk (t)  Xe )
k

(4.12)

Where Xk (t) is the location of the kth Lagrangian point at time t; Xe denotes the equilibrium
k
location of the kth Lagrangian point. An equilibrium location is the initial location of a
Lagrangian point ( Xe  Xk (t)
k

t 0

). To enforce the boundary staying close to its equilibrium

location, a large value for the stiffness coefficient κ of the springs is required. Since a Lagrangian
point only links to its equilibrium point and the relative location to other Lagrangian point does
not matter, this method can be applied to a stationary structure with a complicated boundary.

85

(b)

Movable rigid body

Fogelson and Peskin (1988) applied the elastic forcing method to deal with rigid particle
settling in a Stokes flow. Under the assumption that the structure surface is an elastic boundary,
the rigid body is actually allowed to be slightly deformable. The springs are connected to resist
deformation in the following manner: each Lagrangian point was linked to its adjacent points.
An assistant point located at the center of the rigid body was linked to all Lagrangian points on
the surface. A large stiffness coefficient is used so that the sprints can resist the deformation
caused by the surrounding flow. In this work, instead of using the center point, we proposed a
method to constructe a rigid body using a triangle structure. Triangle is a geometry known to be
with a reliable support to a structure. Fig. 4.2 shows how the springs connect for a circle and an
ellipse with k Lagrangian points. The blue circles are Lagrangian points used to represent the
boundaries. Only a few Lagrangrian points are used in the diagrams for simplicity. To ensure
using the Largrangian points only once for the triangle connection, the number of Lagrangian
points is set to be divided exactly by 3. The red dash lines are springs for triangles. The blue dash
lines are springs connecting the neighbor points. The stiffness coefficient for all the springs is
sufficiently large to resist surface deformation.
Each spring linking to two Lagrangian points generates an opposite force with the same
magnitude. For example, the force generated at point k due to the spring k-3 is computed from

Fk3 (t)  (l k-3 (t)  l e )
k-3

(4.13)

where l k-3 (t) is the length of spring k-3 in two perpendicular directions equal to ||Xk-X3||; l e
k-3
is the reference length of spring k-3 in two perpendicular directions; κ is the spring stiffness
which is sufficient large so that Δlk-3 will keep a length close to its reference length. The force
86

generated at point k is a combination of the forces contributed by spring k-1, k-3, k-6, and k-(k1). This method can be extended to a 3-D case using triangular pyramids.

Figure 4.2: Schematic of spring construction for a movable rigid structure using a Lagrangian
forcing method.

According to Forgelson and Peskin (1988), if a rigid body is under an external force, such as
the gravitational, magnetic, and electrical force, the force acting on the rigid body can be forced
directly on the Lagrangian points by using Ff  Fint  Fext , where Fint is the internal force due to
the springs connecting to the Lagrangian points and Fext is the external forces acting on the rigid
body.

4.2.2 Construction of the Dirac delta function δh
The immersed boundary method involves two mesh systems, including a Eulerian mesh and a
Lagrangian mesh. Shown in Fig. 4.3, the Lagrangian mesh is a set of discrete points used to track

87

the structure surface. Since the Lagrangian mesh does not concise with the Eulerian mesh, a
distribution function is needed for information communication between two mesh systems.
Several distribution functions could be used in the past, such as a narrow Gaussian
distribution (Goldstein et al. 1993), a discrete hat function (Saiki & Biringen 1996), and various
versions of the Dirac direct function (Beyer & Leveque 1992, Lai & Peskin 2000, and Peskin
2002). In this work, a discrete Dirac delta function δh developed by Peskin (2002) is used. The
expression of δh for a 2-D case is assumed to be
h (r) 

1
(r1 )(r2 )
h2

where r1 

(4.14)

x1
x
; r2  2 ; x1 and x2 are two different space components; h denotes the Eulerian
h
h

grid size. The Eulerian grid size h is set to be at least two times larger than the Lagrangian grid
size H to avoid missing when information is converted between mesh systems. The scalar
function (r) has to obey the following criterions of a delta function in a discrete form:
1.

(r) is continuous for all real r,

2.

(r)  0 for | r | 2 ,

3.

 j(r  j)  r
j

for all real r,

Based on the above conditions, the expression for (r) is given as
(r)  0, | r | 2
1
= (3  2 | r |  1  4 | r | 4r 2 ), 0 | r |  1
8
1
= (5  2 | r |  7  12 | r | 4r 2 ), 1 | r |  2
8

88

(4.15)

(r) is a continuous and symmetric distribution plotted in Fig. 4.4. When |r| is larger than 2, (r)

is δh has no contribution to information communication between two mesh systems. This
generates a smooth but somewhat shape distribution at the vicinity of the immersed surface.
Inclusion of more points blurs the boundary and reduces the accuracy while involvement of too
few points might cause an extremely shape distribution which worsen the stability.
Fig. 4.3 shows how the delta distribution is used for information communication between two
mesh systems. Take the Lagrangian point highlighted by a blue circle for example, in the
spreading procedure, the force on the Lagrangian point is distributed through Eq. (4.14) to the
blue dots on the Eulerian mesh. In the interpolation procedure, the same delta distribution
function is used to interpolate the velocity on the blue dots of the Eulerian mesh to the velocity
of the Lagrangian point circled blue.

Figure 4.3: Schematic configuration shows two mesh systems and the Eulerian points involved in
a spreading procedure and an interpolation procedure for one Lagrangian point.
89

Figure 4.4: Distribution of (r) shows how the Dirac delta function works.

4.3 Algorithm
In this section, the algorithm of using the IB method based on the MacCormack scheme is
provided for 2-D cases as the follows.
*

*

1. The predictor process solves the preliminary variables u , v , and ρ
n

n

*

from the
n

compressible N-S equation using the variables at the current time step u , v , and ρ .
n
G* j  Gi, j 
i,

t n
t n
n
n
(Ei 1, j  Ei, j )  (Fi, j1  Fi, j )  tf n
x
y
*

(4.16)
*

2. The preliminary velocity at the Lagrangian points U , V is solved by Eq. (4.17) through
an interpolation procedure. The Lagrangian points are moved to a preliminary location
*

X(s) from the location at the current time X(s) using Eq. (4.18).

90

U* ( Xk ) 

 u* (x)h (x  Xn )h 2
k
xg

(4.17)

h

X*  Xn  tU(Xk )*
k
k

(4.18)

3. Depending on whether the structure is stationary or movable in the fluid flow, the
n

*

preliminary Lagrangian forcing Ff* is obtained from X , X , and κ using Eq (4.12) or Eq
*

(4.13). The preliminary force density term f is solved by Eq. (4.20) given as
Nb

f (x) =  Ff ( Xk , t)* h (x - Xk )H
*

(4.20)

k 1

4. Substitute the source term f * into Eq. (4.11) which is solved for the velocity and density
at the next time step u
n
Gi,1 
j

n+1

n+1

,v

and ρ

n+1

 t *
1 n
t *
t *
*
*
*
Gi, j  Gi, j  x (Ei, j  Ei 1, j )  y (Fi, j  Fi, j1 )  + 2 f
2


(4.21)

5. The location of the Lagrangian points at the time step n+1, Xn 1 is obtained by using the
same equations in Step 2

Un 1 ( Xk ) 

n
 un1(x)h (x  Xk )h 2
xg

(4.22)

h

Xn 1  Xn  tU(Xk )n 1
k
k

(4.23)

6. Again the Lagrangian forcing Ff
term for the next time step f

n+1

n+1

is computed and used to obtain the force density

through Eq. (4.24)

Nb

f (x)n 1 =  Ff ( Xk , t)n 1 h (x - Xk )H

(4.24)

k 1

91

The steps above are repeated as time is matching. To reduce the numerical viscosity, one can
use different combinations of forward/backward difference schemes for the space derivative
discretization. In this work, the combinations are circled in an order shown in Table 4.1.

Table 4.1: Differencing sequence for the space derivative discretization
step
1
2
3
4

predictor

 / x
F
F
B
B

 / y

corrector

 / x

 / y

F
B
B
B
F
B
B
F
F
F
B
F
F: forward B: backward

It is suggested by Tannehill et al. (1997) that circling those arrangements in the chronological
order can reduce the truncation error. For particulate flow simulation or particle problems,
Perrin and Hu (2006) found the arrangement helps avoid asymmetries due to the truncation error.
To maintain second-order accuracy, the derivative terms in E and F are discretized differently.
The terms with derivative in a different direction of E / x and F / x are always differenced
with a central difference while the terms with derivative in the same direction of E / x and
F / x are differenced in an opposite direction. For example, consider the term E3 in the

momentum equation in the y direction:

E3  uv  

u
v

y
x

(4.25)

The term should be differenced in the predictor procedure as follows
n
un  un
vn  vi1,j
n  (uv)n  n i,j1 i,j1   n i,j
E
i,j
i,j
i,j
i,j
2y
x

92

(4.26)

The MacCormack scheme is an explicit time-marching approach. To obtain a converged solution,
the time step has to satisfy the CFL condition. The following semi-empirical stability criterion is
given by Tannehill et al. (1997)
| u | | v |

1
1 
t 

c
 2

(1  2 / Re )  x y
x 2 x 

1

(4.27)

where σ is a safety constant equal to 0.9; Re∆ =min(ρ|u|∆x/μ, ρ|v|∆y/μ) is the minimal mesh
Reynolds number. According to this criterion, the MacCormack scheme will be inefficient for a
small Reynolds number (Re<10). For a moderate Reynolds number, Perrin and Hu (2006)
suggested using t  0.5h / c as the time step. In this work, in order to obtain simulation
convergence, we conduct a dependence study for both the grid size and the time step size.

4.4 Results and discussion
In order to test the performance of the numerical approach discussed above, the following 2-D
cases are studied: (1) Poiseuille flow, (2) flow passing a stationary cylinder, and (3) rigid body
impinging on a wall.

4.4.1 Poiseuille flow
Fig. 4.5 demonstrates the boundary conditions used for a simple Poiseuille flow and a
Poiseuille flow with a stationary cylinder immersed. In Case 2, a cylinder with a diameter equal
to L/6 (dp=0.17) stays at the location of (5/12L, L/2). The computational domain for both cases is
L (L=1) in width and 4L in length. The bounds at the top and at the bottom are wall using a noslip boundary condition (u=0, v=0). The pressure at the left bound (inlet) and at the right bound
(outlet) are P1=101 and P1=101, respectively. The flow viscosity μ is equal to 0.01. The grid
93

number in the channel length Nx and the channel width Ny is [60,240] which creates a uniform
Cartesian mesh for the fluid flow.
Analytical solution to the velocity profile is available for Poiseuille flows. For an arbitrary
cross-section in the stream-wise direction, the velocity profiles are the same and given as
u(y)  cy(1  y)

(4.28)

where c is a constant. Comparison of the velocity profile at the centerline of the channel length
between the analytical solution and the numerical result is shown in Fig. 4.6. The plot shows that
the MacCormack scheme is accurate in predicting the flow velocity.

Figure 4.5: Schematic of a flow passing two parallel planes (a) without an immersed cylinder and
(b) with an immersed cylinder

Fig. 4.7 demonstrates the comparison of the vector and vorticity distribution between the two
cases. The vorticity is a measure of the local ‘spin’ or ‘rotation’ of the fluid and defined as the
curl of the velocity field,    u . For a 2-D flow,
94



v u

x y

(4.29)

In Fig. 4.7a, the vector shows a parabolic profile cross the channel width. Since the flow is
unidirectional, no vortex is observed in the flow domain. The figure on the top of Fig. 4.7b
shows that the stationary cylinder located in the flow domain split the flow into two streams
which are recovered downstream. Meanwhile, the velocity at the vicinity of the stationary
structure is nearly zero. The figure on the bottom of Fig. 4.7b demonstrates that two vortices in
an opposite direction are generated due to the presence of the cylinder.

Figure 4.6: Comparison of the velocity profiles between the analytical solution and the numerical
results cross the channel shows a very good agreement.

95

(a)

(b)
Figure 4.7: Comparison of the vector and the vorticity between (a) a Poiseuille flow (b) a
Poiseuille flow interacted with a stationary structure.
96

Figure 4.8: Comparison of the velocity profiles along Line AB between two cases

Figure 4.9: Comparison of the pressure distribution along Line CD between two cases
97

To study the impact of the stationary cylinder on the fluid flow, the velocity and the pressure
along the line AB and CD in the above diagram are plotted in Fig. 4.8 and Fig. 4.9, respectively.
The solid line in Fig. 4.8 shows a maximum velocity at the center when no cylinder is
interrupting. The dash line with red dots shows that the maximum velocity is reduced due to the
presence of the cylinder. In reality, velocity is equal to zero at the boundary and inside the
cylinder. However, the result shows that the boundary is blurry for the Dirac delta function is
used to smooth the structure surface. Accuracy would be improved by using a finer mesh. A
mesh dependence study is conducted in the following section. In Fig. 4.9, the solid line shows a
linear decrease along the stream-wise direction. For the other case, the presence of the stationary
cylinder causes a high pressure in the front and a low pressure behind the cylinder. These
phenomena are consistent with the observation in experiments.

4.4.2 Flow passing a stationary cylinder
Uniform flow passing a two-dimension stationary circular cylinder is often used as a
benchmark case to evaluate the performance of a numerical method. Lai and Peskin (2000), and
Uhlmann (2005) tackled the stationary cylinder by using the IB boundary method based on a
projection method. Perrin and Hu (2006) enforced the no-slip boundary on the cylinder using a
Taylor expansion and solved the fluid equation by the MacCormack scheme. Experiment of flow
passing a stationary cylinder shows that the flow may present different flow patterns for different
Reynolds number. Flow with a low Reynolds number generates two stable and symmetrical
vortices behind the cylinder. Flow is unstable for Reynolds number in a certain range (typically
5

50<Re<10 ), which leads to a periodic alternate shedding of vortices. This well-known
phenomenon is referred to as “Kármán Vortex Street”.

98

To characterize the flow, dimensionless numbers, such as Reynolds number Re, Strouhal number
Sr, drag coefficient CD, and the lift coefficient CL are introduced. The Reynolds number is
defined as

Re 

f u  d
f

(4.30)

where u∞ is the far-field velocity, and d is the cylinder diameter. The Strouhal number is a
measure to the vortex shedding frequency, defined as

Sr  f c

d
u

(4.31)

where fc is the vortex shedding frequency. The drag coefficient CD and the lift coefficient CL are
to quantify the drag force or the lift force when an exterior flow is passing a blunt body. They are
obtained respectively from

CD 

Fx
2
1/ 2f u d

CL 

Fy

(4.32)

(4.33)

2
1/ 2f u  d

where the drag force Fx and the lift force Fy are the force components exerted on the cylinder in
the direction of the flow velocity and in the direction perpendicular the flow velocity,
respectively.
Fig. 4.10 demonstrates the boundary condition used for the computational domain used for [0,
0] × [40d, 40d]. A circular cylinder with d=0.2m in diameter is located at (6d, 20d). The
boundary conditions on the left, bottom and top bounds are imposed with velocity u = (u∞, 0)

99

where the far-field velocity u∞ is equal to 1. The Neumann boundary condition is applied to the
right bound using
(u ) (v)

0
x
x

(4.34)

The Reynolds number for the flow studied is equal to 100. The flow viscosity µf is 0.002. A
low Mach number (Ma= 0.1) is maintained by using 10 for the speed of sound. Fig. 4.11 plots
the vorticity patterns proceeding with time. At the beginning, two symmetric vortices rotating in
the opposite direction are generated (Fig. 4.11a). With the increase of the time, the vortices are
stretched further downstream (Fig. 4.11b). Fig. 4.11c shows that as the vortices grow with time,
they are distorted due to the instability of the status. The vortices begin detaching from the
cylinder at t=9.0s and form an alternative vortex shedding (Fig. 4.11d). This well-known
phenomenon is referred to as the “Kármán Vortex Street”. With the time growing, the Ká n
rmá
Vortex Street becomes periodic in time (Fig. 4.11e). The above phenomena have a good
agreement with the observation in experiments qualitatively.

100

Figure 4.10: Schematic of the computational domain used for flow passing a stationary cylinder

101

Figure 4.11: The developing vorticity patterns with a changing time shows the flow eventually
breaks into a periodic Ká n Vortex Street.
rmá
102

Table 4.2 compares the drag coefficient CD between the present scheme and the reference data
provided by Liu et al. (1998). The reference drag coefficient is 1.35 for flow with Re=100. A
mesh dependence study is conducted at first by using different grid sizes h equal to 1/32, 1/64,
and 1/128. The time step uses the one suggested by Perrin and Hu (2006), which is t  0.5g h / c .
The table shows that the discrepancy of CD predicted by the elastic forcing method is reduced
drastically from 94.8% to 11.8% when the grid size is decreased from 1/32 to 1/128. Using a
smaller time step 2/5∆t improves the accuracy for the cases with the grid size 1/64 by 14%. Since
no much change is obtained from 1/48 to 1/64 for time step 2/5∆t, it can be concluded that the
solution is mesh dependence. For time step 1/5∆t with the mesh size 1/48, the discrepancy is
reduced from 18.5% to 11.9%. Therefore, the time step size 1/5∆t is recommended to avoid
using a dense mesh. Although dependence study was performed in time and space for
computational convergence, the scheme is inaccurate in predicting CD with a discrepancy of
11.9%.

Table 4.2: Mesh dependence study in space and time for flow with Re=100
Time step
Mesh Size
1/32
1/48
1/64
1/128

t  0.5g h / c
CD
2.63
1.78
1.51

Discrepancy
94.8%
31.9%
11.9%

2/5∆t
CD
1.87
1.60
1.59
-

Discrepancy
38.5%
18.5%
17.8%
-

1/5∆t
CD
1.51
-

Discrepancy
11.9%
-

Reference CD
Lin et al.(1998)
1.35

4.4.3 Particle impinging on a wall
In this section, a rigid disc settling under gravity and impinging on a stationary wall is
simulated using the elastic forcing method. The computational domain is shown in Fig. 4.12 with

103

the dimension of [0 0] × [48dp 48dp]. A periodic boundary condition is imposed on all the
bounds. The domain has grid points of [480 480] to ensure 10 grid points are across the disc. At
t=0, the disc is located at (38.4dp, 24dp). A stationary wall is located below the disc at the cross
section of y=6dp. The disc is equal to 1 in diameter (dp=1). The ratio of the disc density and the
fluid density is p / f  1.1 and the flow viscosity μ is 0.1. Both the wall and the disc are treated
as immersed boundaries. The disc is released under the gravity g= (0,-9.8). The maximum
Reynolds number Remax for this case is Remax≈10.

Figure 4.12: Schematic of computational domain for a disc settling under gravity and a stationary
wall

104

The curves in Fig. 4.13 show the vertical velocity and the vertical location at the center of a
disc as well as the wall location changing with time under gravity. The vertical axis on the left is
velocity while on the right is location. The disc experiences four different stages as the plot
shown. Snapshots of the velocity contour during this process are also shown in Fig. 4.14.

Figure 4.13: The vertical velocity and center location of the disc changing with time shows the
settling disc is (1) accelerated, (2) with constant velocity, (3) decelerated, and (4) at rest.

At the beginning, the velocity of the disc speeds up very fast due to the gravitational force.
The acceleration is slowed down when the drag force is enhanced due to an increasing disc
Reynolds number. When the increasing drag force and the gravity become equilibrium at t=20s,

105

the disc reaches the maximum velocity called the “terminal velocity”. Shown in Fig. 4.14b, the
terminal velocity is 0.97m/s in the negative y direction. In stage 2, the disc is falling with the
terminal velocity until the equilibrium between the drag force and the gravity is broken at t≈30s
due to the presence of the wall. When a particle comes near a wall, the drag force is inversely
proportional to the distance from the disc to the wall (Brenner 1961). This force is also referred
to as “lubrication force” from the wall which is so large when the disc is getting close to the wall
that the disc velocity is reduced dramatically in a short time (stage 3). Finally, the disc is resting
on the top the wall. The curves in Fig. 4.13 show that the disc center and the wall location are
6.502m and 6m, respectively. Therefore, it is reasonable to say that the bottom of the disc is
located right on the stationary wall. This can also be seen in Fig. 4.14d. Although no collision
scheme is used when the disc is approaching to the wall, the elastic forcing method handles the
boundary well. The gravity is balanced through the fluid flow in the vicinity of the disc
contributed by the deformation of the stationary wall.

106

Figure 4.14: Schematic of velocity contour shows a disc is settling under gravity and impinging
on a stationary wall.

4.5 Summary
An algorithm is developed for combining the MacCormack scheme with the IB method. The
technique tracks the structure surface using a set of discrete Lagrangian points. The force
exerted on the Lagrangian points is computed via an elastic forcing method. Considering the
reacting force on the fluid flow in a source term of the momentum equation, the fluid equations
were solved in a rectangular Cartesian mesh. The method is flexible and easy to implement. In
this chapter, the method is tested in the cases of stationary disc in a Poiseuille flow, uniform flow
passing a stationary cylinder, and particle impinging on a stationary wall.

107

For stationary bodies, result shows that the no-slip boundary is handled properly by the
numerical approach. For example, using the elastic forcing method for a sealed stationary
structure can avoid flow penetrating into the body. It is observed that blurry boundary due to the
application of the Dirac delta function can cause numerical viscosity. This would be overcome
by using a finer mesh in the vicinity of the immersed boundary. When two rigid bodies are
approaching closely to each other, the hydrodynamic force generated is handled automatically by
the method without using a collision scheme. This is a great benefit of using a virtual boundary
method. The numerical method is promising and able to qualitatively capture important physical
phenomena observed in experiment, such as Ká n Vortex Street, particle impinging on a wall,
rmá
etc. However, quantitative study of flow passing a stationary cylinder indicates that the
numerical approach is not accurate compared to other models in estimating the drag coefficient
and the lift coefficient. A numerical viscosity is caused due to the intrinsic defect of treating a
rigid surface as an elastic boundary. It seems that to reduce the numerical viscosity, one can use
an extremely large stiffness coefficient for the springs. However using such a large stiffness
coefficient could cause stability issue. To ensure the stability, the time step size ∆t has to be
larger than the characteristic time of spring damping and stretching ts. A very large stiffness
coefficient could dramatically reduce ts and thus put a severe restriction on the time step used. In
order to increase the accuracy while maintain the stability, a fictitious domain method is
introduced to simulate particulate-laden flows in the following chapter.

108

CHAPTER 5
DIRECT NUMERICAL SIMULATION OF PARTICULATE-LADEN FLOW USING AN
FICFICIOUS DOMAIN METHOD BASED ON THE IMMERSED BOUNDARY
METHOD

In Chapter 4, it was found that the elastic forcing method can predict important physical
phenomena, such as Ká n Vortex Street behind a cylinder and particles impinging on a wall.
rmá
However, since the method treats a rigid body as a slightly deformable boundary using a large
stiffness, the technique may exhibit stability problems during the solution. In this chapter, instead
of using the elastic forcing method, the fictitious domain method developed by Uhlmann (2005)
is employed to compute the force density terms used to enforce the no slip boundary condition.
The technique is combined with the MacCormack scheme to simulate particulate-laden flows. A
collision scheme introduced by Glowinski et al. (1999) is used to account for rigid body collision.
Simulations of flow passing a stationary cylinder, particle sedimentation, and inertia-induced
particle sedimentation are performed to test the numerical method. Finally, the method is used to
validate the mixture model in the commercial software ANSYS FLUENT. Study focuses on the
behavior of neutrally buoyant particles in a Couette flow and a Poiseuille flow.

5.1 The fictitious domain method
The immersed boundary method considers the impact of an immersed structure on a fluid
flow in a force density term added to the momentum equation (Mohd-Yusof 1997). For a time
derivative term evaluated at an intermediate time level, n+1/2, then the time-discretized
incompressible N-S equation is given by

un 1  un
 rhsn 1/2  f n 1/2
t

(5.1)

109

where rhs

n+1/2

combines the convective, pressure and viscous terms. To satisfy the no-slip

boundary condition, the desired velocity (v) of a flow imposed on the immersed boundary
generates the force density term given by

f n 1/2 

v n 1  un
 rhsn 1/2
t

(5.2)

Uhlmann (2005) states that Eq. (5.2) can be evaluated at the Lagrangian points representing
the moving body through Dirac delta functions. The corresponding expression is given as

F( Xk )n 1/2 

V( Xk )n 1  U( Xk )

(5.3)

t

where U is the Lagrangian form of the preliminary velocity u
u  un  t  rhsn 1/2

(5.4)

A compressible N-S equation can also be solved for u by using the MacCormack scheme
introduced in Chapter 4. The desired velocity at the kth Lagrangian point V(Xk) can be obtained
from

V(Xk )  uc  ωc  (Xk  xc )

(5.5)

In Eq. (5.5), xc is the coordinate of the geometric center. uc is the translational velocity. For a
2-D case, the angular velocity ωc is a scalar. For situation of stationary rigid boundaries, both the
translational velocity and the angular velocity are set to be zero. When considering a moving
rigid boundary, uc and ωc are determined from the force balance equation and the torque balance
equation, given as

Vc (p  f )uc  f  Fdx  Vc (p  f )g

(5.6)

S

110


Ic (1  f )c  f  ( X  xc )  Fdx
p
S

(5.7)

The derivation for Eq. (5.6)-(5.7) can be found in Appendix B. The discrete form of the above
equations for a 2-D case is
n
n
uc 1  uc
Vc (p  f )
 f  k F( Xk )Vk  Vc (p  f )g
t

(5.8)

n
 n 1  c
Ic (1  f ) c
 f  k ( Xk  xc )  F( Xk )Vk
p
t

(5.9)

ΔVk is the control volume of the kth Lagrangian point equal to the square of the Lagrangian grid
2

size H .
The preliminary velocity on the Lagrangian points, U( Xk , t) , and the force density on the
Eulerian mesh, f(x,t) are obtained through intergrating the Dirac delta function from Eq. (4.15).
We have


xg

u(x, t)h (x - Xk )h 2

(5.10)

f (x, t)=  Fk ( Xk , t)h (x - Xk ) Vk
k 1

(5.11)

U( Xk , t)=

h

Nb

5.2 Collision scheme
Unlike the elastic forcing method which directly moves the Lagrangian points, the fictitious
domain method moves and rotates the geometric center of a rigid structure through the
translational velocity uc and the angular velocity ωc. Therefore, there is a possibility for surface
superposition when two particles approach to each other. This violates the famous Stokes
111

paradox saying that two rigid bodies can never make contact in a finite time in a viscous fluid
because of the infinite “lubrication force” when the distance approaches zero at the last moment
of contact (Zhang et al. 2005). Therefore a collision scheme is necessary to prevent collision
when two rigid structures approach closely to each other.
For particulate-laden flow, we use the collision scheme developed by Glowinski et al. (1999).
Fig. 5.1 demonstrates the method used for two spherical particles and for a particle moving
toward a stationary wall. In Fig. 5.1a, d12 denotes the distance between the center of two
particles P1 and P2. A safety zone with the thickness ls is imposed on each particle to avoid
collision. For particles in radius R1 and R2, the repulsive force Fr is generated when d12 is
smaller than the critical distance ds= R1+R2+ls. For d12> ds, the repulsive force Fr is zero. The
above method can also be used to avoid collision when a rigid particle is approaching to a
stationary wall by adding a ghost particle in the same radius in the wall (shown as in Fig. 5.1b).
The equation for the collision scheme for m rigid particles is given as
0

p
Frip   Fri,j   m 1
2
ji
   ( Xi  X j )(ds  dij )
 ji p

dij >ds

m

(5.12)

dij <ds

where εp is a small stiffness constant; the critical distance between particle Pi and Pj is
ds=Ri+Rj+ls; the distance between the center of Pi and Pj is computed by |Xi-Xj|. The force
acting on the particle Pi is the summation of force exerted by all other particle Pj (j≠i). For
example, Fig. 5.2 plots the repulsive force for two particles with the same radius 1. The thickness
-7

of the safety zone is 0.1. The stiffness constant εp uses a small value 1× . The curve shows
10

112

that when the distance of two particles is larger than 2.1, the repulsive force acting on the particle
is equal to zero. Within the safety zone, the force is activated and is increased exponentially with
the decrease of the particle distance. This force is so large than two particles coming to together
will be push back to their opposite directions. For elliptical particles, the scheme can be applied
by using an equivalent diameter. The equivalent diameter for a ellipsoid is the diameter for a
sphere possessing the same volume of this ellipsoid.

Figure 5.1: Schematic of the collision scheme used to prevent collision between (a) two particles
(b) a particle and a stationary wall

113

Figure 5.2: A curve shows the repulsive force acting on a rigid particle is decreased dramatically
as the increase of particle distance within the safety zone and is zero out of the safety zone.

5.3 Test cases
The performance of the fictitious domain method based on the explicit MacCormack scheme
is tested in the following cases: flow around a stationary cylinder, particle sedimentation, and
inertia-induced particle migration.

5.3.1 Flow passing a stationary cylinder
In Chapter 4, simulation of flow passing a stationary rigid cylinder is performed by using
the elastic forcing method. As a benchmark case, we use it to validate the fictitious domain
method based on the MacCormack scheme. The application of the fictitious domain method on a
114

stationary rigid structure can be simply achieved through setting the desired velocity V(X) (Eq.
(5.5)) equals to zero. The Reynolds number for this case is 100. One can refer to Section 4.4.2
for the detail of the case studied. Numerical results of the drag coefficient CD, the drag
coefficient fluctuation CD', the lift coefficient CL, and the Strouhal number Sr are compared to
the elastic forcing method and other sources in Table 5.1. The mesh and the time step size used
in this case are 1/48 and 1/5∆t (∆t=0.5h/c), respectively.
Compared to the reference data (Liu et al. 1998), the fictitious domain method performs the
best in estimating the drag coefficient with a discrepancy of 0.7%. Lai & Peskin (2000) and
Uhlmann (2005) used the elastic forcing method and the fictitious domain method based on a
project method, respectively. The discrepancies of CD from their schemes are both 7.4%. The
reason that our scheme (the fictitious domain method based on the MacCormack solver)
outmatches the other two models may be due to the small time step used.

Table 5.1: Comparison of results from the current schemes and other methods (Re=100)

Elastic forcing method
Fictitious domain method
Lai & Peskin (2000)
Uhlmann (2005)
Perrin and Hu (2006)
Liu et al. (1998)

CD
1.51
1.34
1.45
1.45
2.0
1.35

Discrepancy [%]
11.9
0.7
7.4
7.4
48.1
-

C D'
±
0.014
±
0.013
±
0.011
±
0.012

CL
±
0.342
±
0.337
±
0.329
±
0.339
±
0.339

Sr
0.165
0.164
0.165
0.169
0.233
0.165

5.3.2 Particle sedimentation
In this section, simulations of disc settling under gravity are conducted. The cases include
single disc sedimentation, two discs drafting-kissing-tumbling (DKT) and 240 particles settling
in a sealed container.
115

(a) Single disc setting under gravity
A single circular disk falling under gravity is simulated using the MacCormack solver with
the elastic forcing method and the fictitious domain method. The circular disc is 0.2m in
diameter (dp= 0.2m). It is located at the original point of the computational domain [-12dp, 12dp]
× [-65dp, 7dp]. The ratio of the disc density ρp and the flow density ρf is 1.5. The mesh size uses
240×
720. The boundary conditions on the top and the bottom are periodic. On the left and right
bounds, the outflow boundary condition from Eq. (4.34) is used. At t=0, the velocity for the disc
and the fluid flow are at rest. At t>0, the disc is released under the gravity g= (0, -9.8m/s2). The
computation is stopped when the disc reaches y=-7m.
Fig. 5.3a shows the comparison of the elastic forcing method and the fictitious domain
method in estimating the center position of the disc settling. The maximum Reynolds number
obtained for this case is 100. At the beginning of the settling process, it can be seen that no
horizontal displacement ∆y is generated for two methods. Due to the increase of the particle
velocity, the Re number is increased and the elastic forcing method predicts a ∆y for t>2.5s. Fig.
5.3b shows the comparison of the disc velocity from two methods. It can be seen that the
velocity of the disc increases with time until the gravity is balanced with the increasing drag
force. Thereafter, the disc velocity becomes a constant called the terminal velocity. Analytical
solution to the terminal velocity is available for comparison. For Re=100, the analytical solution
to the terminal velocity is 1.08m/s (shown as a dot line in Fig. 5.3b). Fig. 5.3b shows that the
terminal velocity from the fictitious domain method has a good agreement with the analytical
solution where a deviation of 3% is obtained. The elastic forcing method underestimates the
terminal velocity with a deviation of 20%. Moreover, the elastic forcing method could cause

116

velocity fluctuation in the vertical direction ux' and a horizontal velocity uy when t>2.5s. The
reason of causing the fluctuation may be due to the participation of large amount of springs to
mimic a rigid body. For the fictitious domain method, no ux' and uy are observed in the results.
This indicates that the method has a good stability.
Fig. 5.4 shows the comparison of the disc center position and the disc velocity between two
methods for Re≈10. For a small Reynolds number, both the elastic forcing method and the
fictitious domain method predict the disc falling on the centerline of the flow domain, which
indicates both the methods are stable for flow with a low Reynolds number. The terminal
velocity from the elastic force method is 14% lower compared to the result from the fictitious
domain method. We know from the last section that the elastic force method tends to
overestimate the drag coefficient CD and thus reduces the terminal velocity. Since the fictitious
domain method performs better in stability and accuracy on rigid particle simulation, we use this
method for particle flow simulation for the rest cases.

117

(a)

Figure 5.3: Comparison of (a) disc position and (b) disc velocity between the elastic forcing
method and the fictitious domain method shows that the fictitious domain method (short dash)
outperforms the elastic forcing method (solid line) in estimating the terminal velocity of a disc
falling in a flow with Re=100.

118

Figure 5.3 (cont’d)
(b)

119

(a)

Figure 5.4: Comparison of (a) disc position and (b) disc velocity between the elastic forcing
method and the fictitious domain method shows that both methods are numerically symmetric in
estimating the terminal velocity of a disc falling in a flow with small Reynolds number (Re=10
for the fictitious domain method).

120

Figure 5.4 (cont’d)
(b)

(b) Particles Drafting-Kissing-Tumbling
In this section, two discs settling under gravity are modeled via the fictitious domain method
based on the MacComack scheme. The discs are with an identical size dp=0.2m. The density
ratio of the discs and the flow p / f is 1.5. The flow viscosity is µ=0.01kg/m· At t=0s, Disc 1
s.
is located at the original point of the computational domain [-12dp, 12dp] × [-7dp, 65dp]. The
mesh size for this case is 240×
720. The center of Disc 2 is located 2dp below the center of Disc 1.
2

When t>0s, two discs are released simultaneously under the gravity g= (0, - 9.8m/s ).
Fig. 5.5 shows the velocity contour when two discs are falling under gravity. When a disc is
settling, a low pressure wake is generated behind the disc. If another disc is caught in this zone,
the disc could speed up due to the reduced drag force. This phenomenon is referred to as
121

“drafting”. In Fig. 5.5a, Disc 1 is drafting behind Disc 2 at t=2.25s. The velocity of Disc 1 (uc1=
0.78m/s) is larger than the velocity of Disc 2 (uc2= 0.70m/s). Two discs are drafting for a while
until the lagging disc catches up with the lead one, called “kissing” shown in Fig. 5.5b. In Fig.
5.5b, two kissing discs travel as a single long body in a position parallel to the flow at an
identical velocity (uc1= uc2= 0.97m/s). The long body is not stable traveling at this position
which would tumble to a position perpendicular to the flow stream (Fig. 5.5c). Two discs are
pushed apart until a stable separation distance is established (Fig. 5.5d). The above phenomena
obtained from the numerical results are consistent with the observation in the experiment
conducted by Joseph et al. (1987).

122

Figure 5.5: Schematic demonstration of two particles (a) drafting, (b) kissing, (c) tumbling and
(d) being apart after tumbling.

(c) Particles settling in a sealed container
240 circular particles settling in a sealed container under gravity is investigated using the
fictitious domain method based on the MacCormack scheme. The collision scheme is used to
avoid collision among the particles and between a particle and the container wall. Fig. 5.6 shows
the velocity vector and the particle position at time t= 0s, t=2.5s, t=5s, t=7s, t=10s, and t=20s,
respectively. The diameter dp of all the particles is 0.2. The density ratio of the particle and the

123

fluid flow p / f equals to 2. The flow viscosity μ is 0.0075. In Fig. 5.6a, the discs are initially
located on the top of a square container with [24.2d p 24.2dp] in dimension. For a sealed
container, the no-slip boundary condition is enforced on all four bounds. The Eulerian mesh size
h uses 1/20d so that there are 20 grids across each particle. The gap between two particles is 3h.
The mesh grids for the computational domain in the x direction and y direction is 484×
484. At
t=0s, the particles at the bottom are marked in red to highlight the interface between two phases.
When t>0s, the particles are settling under the gravity. Since the light phase is located at the
bottom of the computational domain, the heavy particles on the top are pushed by flow which
causes interface instability. This phenomenon was referred to as Rayleigh-Taylor instability

(a) t=0s
Figure 5.6: Schematic demonstration of 240 particles doing sedimentation under gravity in a
sealed container at (a) t=0s, (b) t=2.5s, (c) t=5s, (d) t=7s, and (e) t=20s.

124

Figure 5.6 (cont’d)

(b) t=2.5s

(c) t=5s

125

Figure 5.6 (cont’d)

(d) t=7s

(e) t=20s

126

(Sharp, D. H. 1984). Fig. 5.6b shows the interface is distorted due to Rayleigh-Taylor instability
at t=2.5s. The interface instability keeps growing with two vortices generated close to the wall
pulling down the particles on the sides. These two vortices develop with time and push those
particles around the container center upward (Fig. 5.6c - d). Finally, a stable state is achieved
when all the particles settle down at the bottom of the container (Fig. 5.6e).

5.3.3 Particle behavior
(a) Single elliptical fiber
Jeffery (1922) solved the translation and rotation of a neutrally buoyant ellipsoid in a shear
Stokes flow. The case is used here to validate the present scheme in a 2-D scenario. Consider an
elliptical disc located at the center of a square domain [L L] in Fig. 5.7. The elliptical disc has the
aspect ratio re=1.5. The angle between the major axis and the horizontal direction is ϕ=π/2. The
wall on the top is moving to the right with the velocity Uw=1 while the wall at the bottom is
moving to the left with the same velocity. The shear rate is  =0.4. The particle Reynolds number
e

Rep for an elliptical disc is defined through an equivalent diameter dp and is obtained by
Rep 

U w d pe

(5.13)



To compare with the Jeffery’s solution, simulation of particle Reynolds number equal to 0.4,
4, and 20 are performed. The angular velocity given by Jeffery’s solution is

 2
c ()  2 (re sin 2   cos 2 )
re  1

(5.14)

127

The period is obtained from
2
2  re  1 
T


  re 



(5.15)

Fig. 5.8 shows the comparison of the angular velocity and the corresponding period between
the DNS method and Jeffery’s analytical solution. The angular velocity ω c(ϕ) is plotted against
the angle ϕ. The negative sign indicates that the particle is rotating clockwise. When Rep=0.4,
the curve of the DNS data matches exactly that of Jeffery’s solution. The difference of the period
between two methods is marginal with 0.01s. This indicates that the fictitious domain method
based on the MacCormack scheme is accurate. It can also be concluded that the flow with
Rep=0.4 is small enough to be treated as a Stokes flow. The DNS method underestimates the
angular velocity for the cases Rep=4 and Rep=20, which leads to a longer period time. Jeffery’s
solution was derived based on Stokes flow ( Rep

1 ) and thus cannot handle flow with a large

Rep.

Figure 5.7: Schematic of computational domain for a neutrally buoyant elliptical disc moving in
a Stokes shear flow.
128

Figure 5.8: Comparison of results from the DNS method for different Rep and Jeffery’s solution
shows that the angular velocity and the period are well predicted by the present DNS scheme for
Rep=0.4.

(b) Particle migration
One important phenomenon that could cause non-uniform distribution of a uniform
suspension flow is particle moving across streamlines, referred to as particle migration. The
phenomenon could be resulted from the effects of fluid inertia, particle inertia, particle-particle
interaction, and particle Brownian motion. The effect of particle inertia was discussed in Chapter
2 and Chapter 3. The Brownian motion can be ignored as the size studied in this work is
sufficiently large. Research of interest in this section focuses on neutrally buoyant particle
129

migration induced by fluid inertia. A neutrally buoyant particle is a particle in a condition that its
density is equal to the density of the flow in which it is immersed.
Fig. 5.9 shows the computational domain and the initial location of a neutrally buoyant
particle in a Couette flow. The boundary condition for the top and the bottom bounds are walls
moving with Uw=7 in the opposition directions. The flow Reynolds number Re is 700. The
minor radius of the disc ae is 0.1. The aspect ratio re is equal to 1.5. The particle density and the
fluid density equal to 1. The effect of the particle size is investigated by using larger minor
e

radiuses equal to 1.5ae and 2ae. A circular disc with the equivalent diameter dp of the elliptical
e

disc is also used to study the effect of the particle shape on the particle location. d p is the
diameter when the area of the circular disc is equal to the area of the elliptical disc with ae. Fig.
5.10 plots the dimensionless vertical location (y/L2) with the dimensionless time (t· w/L2) for
U
comparison. The non-dimensional centerline is located at 0.5. The curves show that all the
particles would move to the centerline of the shear flow eventually regardless of the size and the
e

shape of the particles. The curve of the circular disc with dp is close to that of the elliptical disc
which may indicate that the shape effect could be ignored if an equivalent diameter of an
elliptical disc is used for a circular particle. Result also shows that the elliptical discs with a
larger minor radius are faster to reach the equilibrium location as a faster angular velocity is
obtained for a larger particle under a simple shear flow. Rubinow & Keller (1961) claimed that
the lateral force (lift force) is function of the angular velocity.
Fig. 5.11 shows the configuration of the computational domain for a single particle moving in
a Poiseuille flow. A periodic boundary condition is used for the left and the right bounds. A

130

pressure gradient is imposed on the domain which could generate a flow with Re=700. The noslip boundary condition is used to the top and the bottom bounds. Simulations are run for
particles with the same shape and size as the last case. It can be seen from the curves in Fig. 5.12
that the particles regardless of their shape and size reach the vertical location around 0.2 which is
60% of the distance from the centerline of the flow to the bottom wall. This result is consistent
with the observation in experiment of Segré & Silberberg (1962a,b) who stated that, for
undisturbed flow, lateral migration of a neutrally buoyant sphere induced by fluid inertia can
finally attain an equilibrium location about 60% from the centerline to the wall in a Poiseuille
flow. The establishment of the equilibrium position is likely due to the balance between the
lubrication force from the wall and the lift force induced by the fluid inertia. Autal et al. (1991)
stated that when a rigid body is moving close to the wall, the lubrication force is inversely
proportional to the distance between the body and the wall. Unlike the shear flow, the plot shows
that the vertical location the particle is fluctuating around the equilibrium location in the
pressure-driven flow, which indicates that the lift force induced by the fluid inertia and the
lubrication force is fighting back and forth to find a balance.

131

Figure 5.9: Schematic of the computational domain dimension [L1 L2] and the initial location of
an elliptical disc located at (1/4L1, 1/5L2) in a Couette flow with Re=700.

Figure 5.10: Comparison of the particle location changing with time in the vertical direction
shows that particles attain the centerline in a Couette flow regardless of the size and shape of a
particle

132

Figure 5.11: Schematic of the computational domain dimension [L1 L2] and the initial location
of an elliptical disc located at (1/4L1, 7/16L2) in a Poiseuille flow with Re=700.

Figure 5.12: Comparison of the particle location changing with time in the vertical direction
shows that particles attain approximately the location 60% from the centerline to the wall in a
Poiseuille flow regardless of the size and shape of a particle.
133

Fig. 5.13 shows the velocity contour and the concentration distribution of a uniform
particulate-laden flow changing with time in a Couette flow. At t=0s, elliptical particles are
uniformly distributed in the flow domain. The flow is dilute with the volume fraction αp=10%.
Monodisperse particles are used with the minor radius ae=0.1 and an aspect ratio equal to 1.5.
The computational domain is square with the mesh size 250×
250. The Reynolds number Re is
equal to 700. Fig. 5.13b shows that although the particles are moving with a different
translational velocity along an arbitrary vertical cross-section, all the particles are rotating with
an identical angular velocity caused by the single shear rate of a Couette flow. In Fig. 5.13c, the
plot shows that the particles are moving toward the centerline due to the phenomenon of particle
migration across streamlines. This leads to high concentration at the flow center. In Fig. 5.13d-e,
due to the high concentration at the flow center, interaction between the elliptical particles
becomes outstanding. It can be seen that both the velocity profile and the particle concentration
are affected by particle-particle interaction. Finally, Fig. 5.13f shows that the flow forms clusters.
It can be seen from those particles marked in the white circles that certain part of the particles
overlap, which indicates that the collision scheme designed for body in spherical shape is not
working properly for ellipsoids.
Fig. 5.14a shows the velocity contour and the concentration distribution of a uniform
particulate-laden flow changing with time in a Poiseuille flow. The particles in the pressuredriven flow act totally differently from those in the shear flow. Since a Poiseuille flow possesses
a parabolic flow pattern, along an arbitrary vertical cross-section the particles are rotating in a
different angular velocity due to the non-uniform shear rate. Moreover, the particles in the upper
half part of the domain are rotating in an opposite direction of those in the lower half shown as in
Fig. 5.14b. Due to particle migration in the Poiseuille flow, most particles are moving toward the
134

wall (Fig. 5.14c-e) except for those located at the centerline. The particles at the centerline
encounter shear in both directions. They stick near the centerline for a long time before they fall
into the either side of the centerline and move toward the wall. Eventually, all the particles
concentrate somewhere between the centerline and the walls. Fig 5.14f shows that the particles
are clustering which forms larger structures near the wall.

Figure 5.13: Schematic demonstration of the velocity contour and the concentration distribution
of a uniform particulate-laden flow with αp=10% in a Couette flow at t=0s, 5.5s, 100s, 125s,
135s and 160s.

135

Figure 5.14: Schematic demonstration of the velocity contour and the concentration distribution
of a uniform particulate-laden flow with αp=10% in a Poiseuille flow at t=0s, 1s, 30s, 60s, 100s
and 416.0 s

5.4 Validation of the mixture model in FLUENT
Based on the above study, the numerical approach, fictitious domain method with the
MacCormack scheme, is accurate in predicting particle behavior in a viscous flow. The result can
be used as a reliable reference to validate the mixture model in FLUENT. In this section,
particulate-laden flows with a non-uniform concentration and a uniform concentration are
investigated. Focus is put on the performance of the mixture model on predicting the
concentration distribution of neutrally buoyant particles under a simple shear flow.
136

5.4.1 Transport of a non-uniform shear flow
Fig. 5.15 shows the initial concentration distribution of the disperse phase in the
computational domain of [3 1]. The mesh size uses 750×
250 in two methods. Both the carrier
phase and the disperse phase are at rest at t=0s. The disperse phase concentrates in a square
domain with the volume fraction αp equal to 46%. The color in Fig. 5.15a shows the phase
distribution. Red means the region for the highest concentration and blue indicates the location
with zero concentration of the disperse phase. The disperse phase is neutrally buoyant and has an
identical density as the carrier phase (ρp= ρf=1). The wall on the top is moving to the right with
the velocity Uw=1 while on the bottom is moving to the left with the same velocity. A periodical
boundary condition is imposed on the left and the right bounds. The flow is within the laminar
regime with Re=100.
Fig. 5.16 shows the comparison of the disperse phase distribution at t=10s. The concentration
predicted by the mixture model is similar to the result from DNS. Under the shear, the square
shape distribution is changed. It can be seen that the disperse phase closer to the wall is moving
faster due to a higher flow velocity. In Fig. 5.16b, the edge on four corners becomes blurry with
a lighter color which indicates a lower concentration at those places. This distribution is
consistent with the result in Fig. 5.16a, showing as a longer distance among the particles on four
corners. In Fig. 5.17, the distribution is stretched by the flow and forms a stripe at t=25s. Due to
the periodic boundary condition, the disperse phase that passes the right bound will enter the left
bound and vice versa. At t=50s, Fig. 5.18b shows that the concentrated disperse phase is
stretched further with five stripes. The concentration is higher close to the centerline due to

137

slower velocity at the center. This distribution predicted by the mixture model is consistent
qualitatively with the result by the fictitious domain method shown in Fig.5.18a.

Figure 5.15: Initial condition shows the concentration of the dispersed phase in the
computational domain for (a) the present scheme (b) the mixture model

138

Figure 5.16: Comparison of the concentration of the dispersed phase at t=10s for (a) the present
scheme (b) the mixture model

139

Figure 5.17: Comparison of the concentration of the dispersed phase at t=25s for (a) the present
scheme (b) the mixture model

140

Figure 5.18: Comparison of the concentration of the dispersed phase at t=60s for (a) the present
scheme (b) the mixture model

5.4.2 Transport of a uniform unidirectional flow
To validate the performance of the mixture model in transporting a uniform shear flow,
simulations are run for a uniform particulate-laden flow in a Couette flow and a Poieuille flow.
Fig. 5.19 shows a square domain [L L] in which neutrally buoyant particles are uniformly
distributed. The volume fraction αp is 10%. The mesh size uses 250×
250. For the Couette flow,
two walls are moving at a velocity Uw=7, the Reynolds number is 700. Fig. 5.20 shows the
comparison of volume fraction distribution along an arbitrary vertical cross-section between the
fictitious method based on the MacCormack scheme and the mixture model at t=50s and t=100s.
The horizontal axis is a dimensionless location (y/L) and the vertical axis is the flow volume
141

fraction. Since the flow is symmetric to the centerline, only half of the line is plotted for
simplicity. Result from the present scheme shows that at t=50s, the maximum flow concentration
is 16.5% located at 0.2. At t=100s, a higher concentration is observed equal to 19%. The location
for maximum concentration is moved from 0.2 to 0.4 which is closer to the centerline (0.5). On
the other hand, the result shows that the distribution predicted by the mixture model stay at 10%
along the cross-section which indicates that the model is not able to predict the distribution
correctly. Fig. 5.21 shows the comparison of the flow volume fraction between two methods for
the Poiseuille flow with Re=700. Result from the present scheme at t=20s and t=40s shows that
the flow concentration becomes lower close to the centerline (0.5) and higher close to the wall
(0.0). Same conclusion can be draw from the curves that the mixture model is not able to predict
the phenomenon of disperse phase across streamline due to fluid inertia. The mixture model
considers the impact of the disperse phase in a momentum term contributed by the slip velocity.
We know from Chapter 3 that the slip velocity is the velocity between two phases which needs to
be closed by model. The only model in ANSYS FLUENT 12 is called the algebraic slip model
by Ishii (1975) (Eq. 3.16). It can be seen from the equation that the effect of the lateral force
generated due to fluid inertia is not included in the equation which cause the model fail to predict
the phenomenon studied in this work.

142

Figure 5.19: Schematic of the computational domain for a uniform flow with αp=10%.

Figure 5.20: Comparison of volume fraction distribution along a vertical cross-section between
two methods at t=50s and t=100s shows that the mixture model fails to predict a higher
concentration near the center for a uniform particulate-laden shear flow.

143

Figure 5.21: Comparison of volume fraction distribution along a vertical cross-section between
two methods at t=20s and t=40s shows that the mixture model fails to predict a higher
concentration close to the wall for a uniform particulate-laden pressure-driven flow (Re=700).

5.5 Summary
To improve the stability and accuracy of using the immersed boundary method for rigid
structures immersed in a fluid flow, the fictitious domain method was used. The technique was
combined with the MacCormack solver to simulate a flow passing a stationary rigid cylinder.
Result shows that the method is accurate in estimating the drag coefficient. Later, the method,
combined with a collision scheme, was applied to particulate-laden flow simulation. Result
shows that the terminal velocity of a single particle settling was well predicted by the numerical
method. In addition, compared to the elastic forcing method, the present method is more suitable
for a movable rigid body since it can tackle flow with a higher Reynolds number. Result also
144

indicates that the present method can capture two particles DKT (drafting-kissing-tumbling) and
the phenomenon of Rayleigh-Taylor instability at two phase interface. Subsequently, the method
was used to study the behavior of a neutrally buoyant particle in both Couette flow and Poiseuille
flow. Result shows that for a Couette flow, the particle would move toward the centerline while
for a Poiseuille flow, the particle would finally attain an equilibrium location that is 60% of the
distance from the centerline to the wall. Finally, the numerical approach was employed to
validate the mixture model in predicting the disperse phase concentration in a simple shear flow
and a pressure-driven flow. It can be concluded from the result that the mixture model is not able
to predict the concentration distribution correctly. To predict the phenomenon, the constitutive
equation used in the mixture model to close the slip velocity should include the effect of the
lateral force induced by fluid inertia.

145

CHAPTER 6
SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS

6.1 Summary and Conclusions
The work in this dissertation can be divided into two parts. The first part uses CFD modeling
to study turbulent particulate-laden flows passing through curved pipes. The one-way coupling
DPM in FLUENT is used for dilute flow with a low mass loading. The study focuses on the
performance of the turbulent closures based on the RANS equation with different near-wall
treatments on estimating the pressure drop and the particle deposition efficiency in bends. The
effects of the particle Stokes number as well as the bend configuration, including the bend angle,
the bend curvature ratio, and the bend diameter, are investigated. For dilute flows with a high
mass loading, the mixture model in FLUENT is employed to study the effect of the disperse
phase on the carrier phase. The flow patterns, pressure drop and liquid films are investigated.
Bend design is conducted for performance related to pressure drop and deposition efficiency. In
the second part, DNS of rigid structures immersed in a viscous fluid are performed at low
Reynolds number using the IB method. The elastic forcing method and the fictitious domain
method are used to compute the force density term contributed by the immersed structures. The
IB method is combined with the explicit MacCormack scheme to solve the fluid equations. The
performance is tested in the cases including flow passing a stationary cylinder, particle
impinging on a wall, particle sedimentation, and particle behavior in unidirectional flows. Finally,
DNS result is used to validate the mixture model in predicting concentration distribution of
neutrally buoyant particle flows in a Couette flow and a Poiseuille flow. Important observations
associated with these studies are presented below.

146

Turbulent air flow passing through a curved pipe possesses complicated flow patterns that
change with bend configuration and flow Reynolds number. A tight bend with a curvature ratio δ
less than 3 is usually associated with a recirculation zone located at the inner bend wall. Flow
passing a curved pipe generates secondary vortices and the intensity of these secondary vortices
is stronger at the inner wall and can become as high as 35% of the bulk velocity. For flows with
a high mass loading, the patterns are significantly affected by the dispersed phase and become
even more complex. For example, the disperse phase with a large size would reverse the
direction of the secondary vortices or cause multi-pair secondary vortices in a bend cross-section.
When examining the pressure drop of a turbulent air flow through curved pipes, computed
results are obtained from different near-wall treatments based on the k-ε model and the RSM; all
results are close to experimental measurements provided that the mesh meets the y+ requirement
for the wall treatment selected except the local pressure at the outer and inner pipe of the U-bend.
For a flow with a high mass loading, the dispersed phase affects the pressure drop significantly.
The pressure drop is proportional to the volume fraction of a dispersed phase. Comparison
between the numerical results and the empirical models shows the pressure drop estimated by the
model agree well with Paliwoda’s work. Since the deviation of the pressure drop predicted by
different empirical models is large, it is uncertain whether the mixture model is accurate. But it
can be concluded that the model predicts results in a reasonable range.
To increase particle deposition, using one-way coupling flow simulations through curved
pipes, it was found that one can either decrease the curvature ratio, increase the residence time,
or decrease the duct diameter, holding the other two parameters constant. The grade efficiency
for this flow is well predicted by the one-way coupling DPM provided that RSM, EWT, and the
stochastic model are used. The grade efficiency is found to be related to the particle St number,

147

the bend angle θ, and the bend curvature ratio δ and can be estimated through the empirical
model (Eq. (2.15)) developed in this work. For flows with a high mass loading, the liquid film
formed by droplet deposition could be predicted by the mixture model. The film is usually
stratified along the outer bend but it can be inversed to the inner bend of a horizontal pipe when
the flow is with a very small volume fraction αp<0.0037. Both the one-way DPM and the
mixture model are capable of predicting particles moving across streamlines due to high inertia
of a large particle.
As regards DNS of flow interacting with rigid structures, the IB method based on the
MacComack scheme appears to be promising. The elastic forcing method is able to capture
important physical phenomena qualitatively. In addition, the hydrodynamic force generated by
two approaching bodies can be automatically taken care by this method without using a collision
scheme. Nevertheless, approximation used in this method of treating a rigid body as a slightly
deformable structure may lead to numerical viscosity and instability. The accuracy and the
stability of the IB method are improved by using the fictitious domain method. The results
accurately match Jeffery’s solution in predicting the motion of an ellipsoid in Stokes flow and
agree well with the observation in experiment for particle migration in a unidirectional flow.
Finally, it can be concluded from the comparison between the DNS method and the mixture
model that the latter is not able to correctly predict the concentration distribution of neutrally
buoyant particles in a uniform unidirectional flow. To capture the phenomenon, the closure for
the slip velocity should be modified to include the effect of the lateral force caused by fluid
inertia.

148

6.2 Recommendations for future work
In the present study, the fictitious domain method with the MacCormack scheme was used to
simulate particulate-laden flow and used to validate the mixture model in FLUENT. Suggestions
are listed in the following for potential extensions of the present study.
1. To avoid collision of two approaching particles, Glowinski’s collision scheme is used to
provide a repulsive force. The method is designed for rigid spherical bodies. Using an
equivalent diameter, one can apply the method to an ellipsoid. This is a good
approximation for ellipsoids with a small aspect ratio. However, for ellipsoids with a
large aspect ratio or rigid bodies in an arbitrary shape, the method can cause problems,
such as particles superposition, and incorrect particle behavior. To deal with particles of
arbitrary shape, a more sophisticated collision strategy such as a lubrication model by
Maury (1997) is needed.
2. To predict particle migrating across streamlines in a unidirectional flow using the mixture
model in FLUENT, the constitutive equation for the slip velocity has to be modified to
include the lateral force caused by fluid inertia. Rubinow & Keller (1961) suggested to
compute this force from

FL  (dp / 2)3 ωc  uc

(6.1)

The parameters in the above equation are all known except the angular velocity ω c . This
unknown quantity can be modeled using the DNS data obtained from the numerical
method studied in this work.
3. Code implemented in this work is applied to 2-D cases. There is a possibility to extend
the code to 3-D application in the turbulent regime which would be more valuable to
validate the mixture model in FLUENT.
149

4. To increase the computational efficiency, parallel computing technique is necessary. One
can use the parallel computing strategy developed by Wang et al. (2008) to implement
immersed boundary method discussed in this work to simulate particulate-laden flows.

150

APPENDICES

151

APPENDICES

A. Realizable k-ε equation for two-phase flows

The modified transport equation for k is given as

   ua k     ( T k)  2 TSuf 
k



p
f  w

(2 p k

g T p
generation of k
due to the gravity

w
 D pu f ,p p )  
 w  l

(A-1)

dissipation due to the water droplets

The modified transport equation for ε is presented as

  uf     ( T ) 


 



p
f  w

( p 

w
)
 w  l

dissipation due to the water droplets

 C2
( 
g T p
)
k  T
generation of 
due to the gravity

(A-2)

where τl is the Lagrangian time-scale, given as
l  0.3

k


(A-3)

The turbulent viscosity for the air phase is defined as

k2
T  C


(A-4)

To ensure the realizability (that is positivity of normal stresses in the flow domain), the turbulent
coefficient Cμ is computed as follows:

152

C 

1
A 0  As

k
S


A0  4.04; As  6 cos 
1
  cos1( 6W)
3

W

SijS jkSki
S3

1  u j u i 

, S  SijSij , Sij  

2  x i x j 



The k and ε transport equations for two-phase flows are similar to those for single-phase flows
except two extra two terms due to existence of the water droplets and the external field such as
gravity.

153

B. Analysis of force acting on a particle

p Vcuc  f
Icωc  f

 τ  nd  (p  f )Vcg

(B-1)

S

 (X  xc )  (τ  n)d

(B-2)

S

where where ρf and ρp are the fluid density and the interior structure density, respectively; Ic is
the moment of inertia of the object; Vc is the volume of the rigid body; τ is the shear stress tensor
and n denotes the outward normal vector of the rigid structure. According to the divergence
theorem, the surface integral on the right hand size of Eq. (B-1) and Eq. (B-2) are

d

 τ  nd   Fdx  dt  udx

S

S

(B-3)

S

d

 (X  xc )  (τ  n)d   (X  xc )  Fdx  dt  (X  xc )  udx

S

S

(B-4)

S

The first term on the right hand size of Eq. (B-3) and Eq. (B-4) are simply obtained by summing
up the force F (Xk) and the torque (Xk-xc)× k for all the Lagrangian points. According to
F
Uhlmann (2003), the second term of Eq. (5.8) is the rate-of-change term which satisfies a rigidbody motion on the structure surface, obtained

d
 udx  Vcuc
dt S

(B-5)

The second term of Eq. (B-5) is the angular momentum changed due to non-rigid motion of fluid
inside the object domain, dΩ. Under the assumption of rigid-body motion inside the object
domain, this term is approximated to

154

I
d
S (X  xc )  udx  c ωc
dt
p

(B-6)

155

BIBLIOGRAPHYS

156

BIBLIOGRAPHYS

[1] Abuzeid, S., Busnaina, A.A., and Ahmadi, G. (1991) Wall deposition of aerosol particles in a
turbulent channel flow. Journal of Aerosol Science, 22 (1), 43-62.
[2] Adib, M.A.H.M., Osman, K. and Jong, R.P. (2010) Analysis of blood flow into the main
artery via mitral valve: fluid structure interaction model, 2010 International Conference on
Science and Social Research.
[3] ANSYS FLUENT 12.1 Documentation. (2009). Theory Guide.
[4] Autal, S.P., Lahey, R.T., and Flaherty, J.E., (1991) Analysis of phase distribution in fully
developed laminar bubbly two-phase flow, International Journal of Multiphase Flow, 7, 635652.
[5] Azzopardi, B.J. (1985) Drop sizes in annular two-phase flow, Experiments in Fluids, 3, 5359.
[6] Azzi, A., Friedel, L.and Belaadi, S. (2000) Two-phase gas/liquid flow pressure loss in bends,
Forschung im Ingenieurwesen, 65, 309-318.
[7] Banerjee, S., Rhodes, E. and Scott, D.S. (1967) Film inversion of concurrent two-phase flow
in helical coils, AIChE Journal, 13 (1), 189-191.
[8] Berger, S.A., Talbot, A., and Yao, L.-S. (1983) Flow in curved pipes. Annual Review of
Fluid Mechanics, 15, 461-512.
[9] Beyer R.P. (1992) A computational model of the cochlea using the immersed boundary
method, Journal of Computational Physics, 98, 145-162.
[10] Berrouk A.B., and Laurence D. (2008) Stochastic modeling of aerosol deposition for LES of
90°bend turbulent flow. International Journal of Heat and Fluid Flow, 29, 1010-1028.
[11] Boussinesq, J. (1877) “Essai sur la théorie des eaux courantes,” Mémoires présentés par
divers savants à l’Acad. Sci. Inst. Nat. France, 23, 46-50.
[12] Brenner, H. (1961) The slow motion of a sphere through a viscous fluid towards a plane
surface, Chem Eng Sci, 16, 242–251.
[13] Breuer, M., Baytekin, H.T., and Matida, E.A. (2006) Prediction of aerosol deposition in 90°
bends using LES and an efficient Llagrangian tracking method. Journal of Aerosol Science,
37, 1407-1428.
157

[14] Brockmann, J.E. (1993) Sampling and Transport of Aerosol, in Aerosol Measurement:
Principles, Techniques, and Applications. Van Nostrand Reinhold: New York.
[15] Cash, J.R., and Karp, A.H. (1990) A variable order Runge-Kutta method for initial value
problems with rapidly varying right-hand sides. ACM Transactions on Mathematical
Software, 16, 201-222.
[16] Chang, S.M., Humphrey, J.A.C., and Modavi, A. (1983) Turbulent flow in a strongly curved
u-bend and downstream tangent of square cross-sections. PCH Physico-chemical
Hydrodynamics, 4, 243-269.
[17] Chisholm, D. (1971) Prediction of pressure drop at pipe fittings during two-phase flow, 13th
Int. Cong. of Refrigeration, Washington, 2, 781-789.
[18] Clarke, D., Salas, M., and Hasan, H. (1996) Euler calculations for multi-element airfoils
using Cartesian grids. AIAA, 24, 1128-1135.
[19] Cheng, Y-S., and Wang C-S. (1975) Inertial deposition of particles in a bend. Journal of
Aerosol Science, 6, 139-145.
[20] Crane V. (1957) Flow of fluids through valves, fittings and pipe, Technical Paper, No. 410.
[21] Crane (1957) Flow of fluids through valves, fittings and pipe, Technical Paper No. 410.
[22] Crowe C. Sommerfeld, M., and Tsuji, Y. (1998) Multiphase flows with droplets and particles,
CRC Press.
[23] Cunningham, E. (1910) On the velocity of steady fall of spherical particles through fluid
medium. Proceedings of the Royal Society, 83, 357-65.
[24] Dean, W.R. (1928) Fluid motion in a curved channel. Proceedings of the Royal SocietyLondon Series A, 121, 402-420.
[25] Elghobashi, S. (1991) Particle-laden turbulent flows: direct simulation and closure models,
Applied Scientific Research, 48, 301-314.
[26] Eustice, J. (1910) Flow of water in curved pipes. Proceedings of the Royal Society-London
Series A, 84, 107-118.
[27] Fitzsimmons, D.E. (1964) Two-phase pressure drop in piping components, HW-80970, Rev.
1, General Electric Hanford Laboratories, Richland, OH.
[28] Fadlun, E.A., Verzicco, R., Orlandi, P., and Mohd-Yusof, J. (2000) Combined immersedboundary finite-difference methods for three-dimensional complex flow simulations, Journal
of Computational Physics, 161, 35-60.

158

[29] Fogelson, A. and Peskin, C. (1988) A fast numerical method for solving the threedimensional Stokes’ equations in the presence of suspended particles, Journal of
Computational Physics, 79, 50–69.
[30] Gibson, M.M., and Launder, B.E. (1978) Ground effects on pressure fluctuations in the
atmospheric boundary layer. Journal of Fluid Mechanics, 86, 491-511.
[31] Glowinski, R., Pan, T.-W., Hesla, T., and Joseph, D. (1999) A distributed Lagrange
multiplier/fictitious domain method for particulate flows, International Journal of Multiphase
Flow, 25, 755–794.
[32] Goldstein, D. Handler, R., Sirovich, L. (1993) Modeling a no-slip flow boundary with an
external force field, Journal of Computational Physics, 105, 354-366.
[33] Gore, R.A. and Crowe, C.T. (1989) Effect of particle size on modulating turbulent intensity,
International Journal of Multiphase Flow, 15, 279- 285.
[34] Hart, J., Ellenberger, J. and Hamersma, P.J. (1988) Single- and two-phase flow through
helically coiled tubes, Chemical Engineering Science, 43 (4), 775-783.
[35] Hetsroni, G. (1989) Particle-turbulence interaction, International Journal of Multiphase Flow,
15, 735- 746.
[36] Hewitt, G.F. and Robertson, D.N. (1969) Studies of Two-Phase Flow Patterns by
Simultaneous X-ray and Flash Photography, Rept AERE-M2159, UKAEA, Harwell.
[37] Ho, B.P., and Leal, L.G. (1973) Inertial migration of rigid spheres in two-dimensional
unidirectional flows, Journal of Fluid Mechanics, 65, 365-400.
[38] Hoang, K. and Davis, M.R. (1984) Flow structure and pressure loss for two-phase flow in
return bends, ASMEJ. Fluids Eng., 106, 30-37.
[39] Hossain, A. and Naser, J. (2004) CFD investigation of particle deposition around bends in a
turbulent flow, 15th Australasian Fluid Mechanics Conference, the University of Sydney, 1317.
[40] Hö
fler, K., and Schwarzer, S. (2000) Navier-Stokes simulation with constraint forces: Finitedifference method for particle-laden flows and complex geometries, Physical Review E, 61
(6), 7146-7160.
[41] Ishii, M. (1975) Thermo-fluid dynamic theory of two-phase flow (book), Paris: Eyrolles.
[42] Ito, H. (1959) Friction factors for turbulent flow in curved pipes. ASME Journal of Basic
Engineering, 81, 123-134.
[43] Ito, H. (1987) Flow in curved pipes. JSME International Journal, 30, 543-552.
159

[44] Jang, D.S., and Acharya, S. (1988) Calculation of particle dispersion due to turbulence in
elliptic flows. AICHE, 34, 514 - 518.
[45] Jeffery, G.B. (1922) The motion of ellipsoidal particles immersed in a viscous fluid,
Proceedings of the Royal Society A, 102, 161-179.
[46] Johansen, S.T., Anderson, N.M. and de Silva, S.R. (1990) A two-phase model for particle
local equilibrium applied to air classification of powers, Power Technology, 63, 121-132.
[47] Joseph, D.D., Fortes, A.F., Lundgreen, T.S., and Singh, P. (1987) Nonlinear mechanics of
fluidization of beds of spheres cylinders and disks in water, in Advances in Multiphase Flow
and Related Problems (Papanicolau, G. ed.), SIAM, 101–122.
[48] Kim, S.-E., and Choudhury, D. (1995) A near-wall treatment using wall functions sensitized
to pressure gradient, ASME, 217, Separated and Complex Flows.

[49] Kim, J., Kim, D. and Choi, H. (2001) An immersed-boundary finite-volume method for
simulations of flow in complex geometries, Journal of Computational Physics, 171 132–150.
[50] Kim, Y., Lim, S., Raman, S., Simonetti, O.P. and Friedman, A. (2009) Blood flow in a
compliant vessel by the immersed boundary method, Annals of Biomedical Engineering, 37
(5), 927-942.
[51] Kim, Y. and Peskin, C.S. (2006) 2-D parachute simulation by the immersed boundary
method. The SIAM Journal on Scientific Computing, 28 (6), 2294-2312.
[52] Kim, Y. and Peskin, C.S. (2009) 3-D parachute simuation by the immersed boundary method.
Computers and Fluids, 38, 1080-1090.
[53] Kovacs, S. J., McQueen, D. M., and Peskin, C. S. (2001) Modeling cardiac fluid dynamics
and diastolic function, Philos. Trans. R. Soc. Lond. Ser. A, 359, 1299–1314.
[54] Kundu, P.K., and Cohen, M.I. (2004) “Fluid Mechanics,” third ed., 686-687.
[55] Lai, M.-C., and Peskin, C.S. (2000) An immersed boundary method with formal secondorder accuracy and reduced numerical viscosity, Journal of Computational Physics, 160, 705719.
[56] Launder, B.E., and Sharma, B.I. (1972) Application of the energy dissipation model of
turbulent to the calculation of flow near a spinning disc. Letters in Heat and Mass Transfer, 1
(2), 131-138.
[57] Launder,B.E., and Spalding, D.B. (1974) The numerical computation of turbulent flows.
Computer Methods in Applied Mechanics and Engineering, 3, 269-289.
160

[58] Launder, B.E., Reece, G.J., and Rodi, W. (1975) Progress in the development of a Reynoldsstress turbulence closure. Journal of Fluid Mechanics, 68 (3), 537-566.
[59] Lee, C. (2003) Stability characteristics of the virtual boundary method in three-dimensional
applications, Journal of Computational Physics, 184, 559–591.
[60] Liu, C., Zheng, X. and Sung, C. (1998) Preconditioned multigrid methods for unsteady
incompressible flows, Journal of Computational Physics, 139, 35-57.
[61] Lockhart, R.W. and Martinelli, R.C. (1949) Proposed correlation of data for isothermal twophase flow, two-component flow in pipes, Chem. Eng. Prog., 45(1), 39-45.
[62] Longest, P.W., and Oldham, M.J. (2008) Numerical and experimental deposition of fine
respiratory aerosols: Development of a two-phase drift flux model with near-wall velocity
corrections. Journal of Aerosol Science, 39, 48-70.
[63] MacCormack, R.W. (1969) The effect of viscosity in hypervelocity impact cratering, AIAA
paper, 69-354.
[64] Majumdar, S., Iaccarino, G., and Durbin P.A. (2001) RANS solver with adaptive structured
boundary non-conforming grids, Annual Reviews Briefs, Center for Turbulence Research
353–364.
[65] Mandal S.N. and Das, S.K. (2001) Pressure losses in bends during two-phase gas-Newtonian
liquid flow, Ind. Eng. Chem. Res., 40, 2340-2351.
[66] Manninen, M., Taivassalo, V., Kallio, Sirpa, and Akademi, Å.(1996) On the mixture model
for multiphase flow, VTT Publications.
[67] Matida, E.A., Finlay, W.H., Lange, C.F., and Grgic, B. (2004) Improved numerical
simulation of aerosol deposition in an idealized mouth-throat. Journal of Aerosol Science, 35,
1-19.
[68] Marcus, R., Leung, L., Klinzing, G., and Rizk, F. (1990) Pneumatic Conveying of Solids,
Chapman and Hall, London.
[69] Maury, B. (1997) A many-body lubrication model, C.R. Acad. Sci., Paris, Serie IMathematics, 325, 1053-1058.
[70] Maxey, M.R., and Riley, J.J. (1983) Equation of motion for a small rigid sphere in a
nonuniform flow. Physics of Fluids, 26 (4), 883-889.

161

[71] McFarland, R.A., Gong, H., Muyshondt, A., Wente, W.B., and Anand, N.K. (1997) Aerosol
deposition in bends with turbulent flow. Environmental Science and Technology, 31, 33713377.
[72] McQueen, D. M. and Peskin, C. S. (2001) Heart simulation by an immersed boundary
method with formal second-order accuracy and reduced numerical viscosity, in Mechanics
for a New Millennium, Proceedings of the International Conference on Theoretical and
Applied Mechanics (ICTAM) 2000, H. Aref and J. W. Phillips, eds., Kluwer Academic
Publishers, Norwell, MA.
[73] Miyoshi, M., Mori, K. , and Nakamura, Y. (2009) Numerical Simulation of Parachute
Inflation Process by IB Method, Journal of the Japan Society for Aeronautical and Space
Sciences, 57 (670), 419-425.
[74] Mohd-Yosuf, J. (1997) Combined immersed boundary/B-spline methods for simulation of
flow in complex geometries. Annual Research Briefs, Center for Turbulence Research, 317328.
[75] Morsi, S.A. and Alexander, A.J. (1972). An investigation of particle trajectories in two-phase
flow systems. Journal of Fluid Mechanics, 55 (2), 193-208.
[76] Norstebo, (1986) A. Pressure drop in pipe components in two-phase gas-liquid flow,
Doctoral Thesis, The Norwegian Institute of Technology, NTH, Division of R.
[77] Pacull, F. (2006) “A numerical study of the immersed boundary method and application to
blood flow,” Ph.D. Thesis, University of Houston.
[78] Paliwoda A. (1992) Generalized method of pressure drop of refrigerants. Rev. Int. Froid, 15
(2), 120-126.
[79] Patankar, N., Singh, P., Joseph, D., Glowinski, R., and Pan, T.-W. (2000) A new formulation
of the distributed Lagrange multiplier/fictitious domain method for particulate flows,
International Journal of Multiphase Flow, 26, 1509–1524.
[80] Patankar, S.V., and Spalding, D.B. (1972) A calculation procedure for heat, mass and
momentum transfer in three-dimensional parabolic flows. International Journal of Heat Mass
Transfer, 15, 1787-1806.
[81] Peskin, C.S. (1972) Flow patterns around heart valves: a digital computer method for solving
the equations of motion. Ph.D. Thesis, Albert Einstein College of Medicine.
[82] Peskin, C.S. (2002) The immersed boundary method, Cambridge University Press: Acta
Numerica 2002, 11, 479-518.
[83] Peters, T.M., and Leith, D. (2004) Particle deposition in industrial duct bends. The Annals of
Occupational Hygiene, 48 (5), 483-490.
162

[84] Perng, C.Y., Street, R.L. (1989) Three-dimensional unsteady flow simulations: Alternative
strategies for a volume-averaged calculation. International Journal for Numerical Methods in
Fluids, 9 (3), 341-362.
[85] Perrin, A. and Hu, H.H. (2006) An explicit finite-difference scheme for simulation of moving
particles, Journal of Computational Physics, 212, 166-187.
[86] Picart, A., Berlemont, A. and Gouesbet, G. (1986) Modelling and predicting turbulence fields
and the dispersion of discretre particles transported by turbulent flows, Int. J. Multiphase
Flow, 12 ( 2), 237-261.
[87] Pourahmadi, F. (1982) Turbulence modeling of single-and two-phase curved-channel flows,
Ph. D. Thesis, University of California-Berkeley.
[88] Pui, D.Y.H., Romay-Novas, F., and Liu, B.Y.H. (1987) Experimental study of particle
deposition in bends of circular cross section. Aerosol Science and Technology, 7, 301-315.
[89] Saiki, E. and Biringen, S. (1996) Numerical simulation of a cylinder in uniform flow:
application of a virtual boundary method, Journal of Computational Physics, 123, 450–465.
[90] Schiller, L. and Naumann, Z. (1935) A drag coefficient correlation. Z. Ver. Deutsch. Ing. 77,
318-327.
[91] Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., and Zhu, J. (1995) A New Eddy-Viscosity
Model for High Reynolds Number Turbulent Flows-Model Development and Validation,
Computers and Fluids, 24, 227–238.
[92] Sharp, D.H. (1984) An overview of Rayleigh-Taylor instability. Physica D: Nonlinear
Phenomena, 12, 3-18.
[93] Segré G., Silberberg, A. (1962a.) Behaviour of macroscopic rigid spheres in Poiseuille flow.
,
Part 1, Journal of Fluid Mech. 14, 115-135.
[94] Segré G., Silberberg, A. (1962b.) Behaviour of macroscopic rigid spheres in Poiseuille flow.
,
Part 2, Journal of Fluid Mech. 14, 136-157.
[95] Spurk, J.H. (1997). Fluid Mechanics. Springer: New York.
[96] Starkey, T.V. (1956) Viscosity and action, British Journal of Applied Physics, 7, 448-450.
[97] Sookprasong, P. (1980) Two-phase flow in piping components, MSc Thesis, University of
Tulsa.

163

[98] Repetti, R.V. and Leonard, E.F. (1966) C~hem Engng Frog. Symp. 8cr. 62, 80.
[99] Rubinow, S.I. and Keller, J.B. (1961) The transverse force on a spinning sphere moving in a
viscous fluid, Journal of Fluid Mechanics, 11, 447-459.
[100] Wolfshtein, M. (1969) The velocity and temperature distribution of one-dimensional flow
with turbulence augmentation and pressure gradient, Int. J. Heat Mass Transfer, 12, 301-318.
[101] Sudo, K., Sumida, M., and Hibara, H. (1998) Experimental investigation on turbulent flow in
a circular-sectioned 90-degree bend. Experiments in Fluids, 25, 42-49.
[102] Sudo, K., Sumida, M., and Hibara, H. (2000) Experimental investigation on turbulent flow
through a circular-sectioned 180°bend. Experiments in Fluids, 28, 51-57.
[103] Theofanous, T.G. and Sullivan, J. (1982) Turbulence in Two-Phase Dispersed Flows, Journal
of Fluid Mechanics, 116, 343-362.
[104] Tannehill, J.C., Anderson, D.A., and Pletcher, R.H. (1997) “Computational Fluid Mechanics
and Heat Transfer,” second ed., Taylor & Francis, Washington, DC.
[105] Thomson, J. (1876) On the origin of windings of rivers in alluvial plains, with remarks on the
flow of water round bends in pipes. Proceedings of Royal Society-London Series A, 25, 5-8.
[106] Tunstall, M.J., and Harvey, J.K. (1968) On the effect of a sharp bend in a fully developed
turbulent pipe-flow. Journal of Fluid Mechanics, 34, 595-608.
[107] Uhlmann, M. (2003) First experiments with the simulation of particulate flows, Technical
Report No. 1020, CIEMAT, Madrid, Spain, ISSN 1135-9420.
[108] Uhlmann, M. (2005) An immersed boundary method with direct forcing for the simulation of
particulate flows, Journal of Computational Physics, 209, 448- 476.
[109] Vigmond, E. J., McQueen, D. M. and Peskin, C. S. (2003) The cable as a basis for a
mechanoelectric whole heart model, Internat. J. Bioelectromagnetism, 5, 183-184.
[110] Wang, Z., Fan J., and Luo, K. (2008) Parallel computing strategy for the simulation of
particulate flows with immersed boundary method, Science in China Series E: Technological
Sciences, 51, 1169-1176.
[111] Wilson, R.E., Mcadams, W.H., and Seltzer, M. (1922) The flow of fluids through
commercial pipelines. The Journal of Industrial and Engineering Chemistry, 14, 105-119.
[112] Wolfshtein, M. (1969) The velocity and temperature distribution of one-dimensional flow
with turbulence augmentation and pressure gradient. International Journal of Heat Mass
Transfer, 12, 301-318.

164

[113] Wu, J., and Shu, C. (2010) An improved immersed boundary-lattice Boltzmann method for
simulating three-dimensional incompressible flows, Journal of Computational Physics, 229,
5022- 5042.
[114] Zhang, P., Roberts, R.M. and Bé
nard, A. (2012) Computational Guidelines and an Empirical
Model for Particle Deposition in Curved Pipes Using an Eulerian-Lagrangian Approach,
Journal of Aerosol Science, 53,1-20.
[115] Zhang, W., Noda, R., Horio, M. (2005) Evaluation of lubrication force on colliding particles
for DEM simulation of fluidized bed, 158, 92-101.

165