Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Comparison between DNS and DEMCFD coupling mesoscopic simulation
for 2D spouted fluidized bed
Takuya Tsuji*, Hirotaka Yada, Kaoru Yoshikawa and Toshitsugu Tanaka
Department of Mechanical Engineering, Osaka University
21 Yamadaoka, Suita, Osaka 5650871, Japan
tak @mech.eng.osakau.ac.jp
Keywords: GasSolid Flow, Fluidized Bed, Direct Numerical Simulation, Immersed Boundary Method, Discrete Element
Method, DEMCFD coupling mesoscopic model
Abstract
The direct simulation of flows including dense solid particles such as that in gasfluidized beds is a challenging problem. A
coupling method between discrete element method (DEM) and immersed boundary ( IB ) method is developed in this paper.
In gasfluidized beds, particles take complex arrangements by the interactions inbetween particles, particleswall in addition
to particlegas flows. Gas flows go through narrow gaps inbetween particles and it is still difficult to capture these
microscopic flows accurately while it has not been discussed well up to the present. In this paper, resolutiondependency
studies are performed for several test problems in that hydrodynamic interactions between particles are important. Besides, a
direct simulation of a twodimensional gasfluidized bed under a spouting condition is performed using IBDEM method. To
enable direct comparison, a DEMCFD mesoscopic model simulation is also performed under the same condition and results
are compared.
Introduction
Gasfluidized beds are widely used in industrial processes.
Flows in the beds have a multiscale structure. The
concentration of solid particle is high and interactions
inbetween dense particles, particlewall and particlegas
flow induce a spontaneous formation of internal
characteristic flow structures such as bubbles that are far
larger than the particle. These internal structures make the
flows unstable and complex and are influential for the entire
behavior of the flow. To enable advanced engineering
designs, it is important to know the behavior of internal
flow structures in the bed. The enhancement of our
knowledge on phenomena occurring in each scale under an
overall multiscale structure shall help the essential
understanding of this type of flows. Due to the existence
of dense particles, however, it is still difficult to observe the
internal flows at each scale directly and a reliable numerical
model has been desired. The coupling method between
discrete element method and computational fluid dynamics
( DEMCFD ) originally proposed by Tsuji et al. (1993) has
been applied to this type of flows successfully. From its
concept, this model can be regarded as a mesoscopic one.
The size of fluid calculation cell used in the model is larger
than particles and enough smaller than mesoscopic
characteristic structures such as bubbles. Hence, only flow
phenomena existing in meso and macroscales are directly
resolved and microscopic flows occurring around each
particle are not. In the model, the effects due to
microscopic flows are taken into account indirectly by using
empirical equations.
The momentum, heat and mass transfers between gassolid
phases are occurring through each particle's surface mainly.
In addition to meso and macroscopic behaviors, it is also
important to advance our understandings on microscopic
flow phenomena. A numerical technique which enables a
direct observation of phenomena occurring in a microscopic
scale is required.
Direct simulations of the flow including solid particles
become feasible (e.g., Kajishima et al., 2001; Pan et al.,
2002 and Uhlmann, 2005). The simulations of a
gasfluidized bed including dense solid particles are tried by
several researchers recently (van der Hoef et al., 2008 and
Kuwagi et al., 2009), however, it is still a challenging
problem.
The purpose of our study is to establish a reliable direct
simulation technique which enables detailed observations of
the flow including microscopic phenomena. In the present
study, a coupling method between discrete element method
and immersed boundary method ( DEMIB ) is developed.
In IB approaches, it is usual to use a fixed Cartesian grid
system and the accuracy of calculation heavily depends on
the resolution we adopted. In dense flows, particles are
interacting with each other and complex arrangements are
formed during a fluidization. It is usual that an
interparticle distance is smaller than a particle diameter and
it is needed to resolve the flows existing in narrow gaps
inbetween particles. The resolution requirement might be
more crucial comparing to dilute particleladen flows,
however, it has not been discussed well up to the present.
In the present study, resolutiondependency studies are
performed for DEMIB coupling method at first. In
addition to fluid drag forces working on a single particle,
resolutiondependency studies are performed for adjacent
paired particles and particles in randomlypacked bed.
Besides, by using DEMIB method, a calculation of a
twodimensional gasfluidized bed under a spouting
condition is performed. For comparisons, a calculation
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
using DEMCFD coupling mesoscopic model is also
performed.
Governing equations
The governing equations of particleladen incompressible
Newtonian flow are the continuity and NavierStokes
equations.
V.uf = 0, (1)
Du,
P D = Vr+pfg (2)
where pf denotes the fluid density, ufthe fluid velocity, rthe
stress tensor:
T= plI+p Vu+ (Vu) (3)
where p the static pressure, p, the fluid viscosity. pf and yf
are assumed to be constant. The last term of the righthand
side of Eq. (2) represents external forces such as the gravity.
The motion of a spherical solid particle is given by an
equation of momentum and angular momentum in a
Lagrangian frame of reference,
d(mpUp) Fc + Ff +G (4)
dt
dl( IM w, )
d(t = M +M +NP, (5)
where up is the translational velocity of the particle, mp is
the angular velocity of particle rotation, mp the mass, Ip the
inertia tensor given by Ip=(2/5)a2mj for the particle of
radius a. F, is a contact force with other particle and wall,
Ff a fluid force and Gp an external force. Me, Mf and Np
are the corresponding torques, respectively. The
summations in Eqs. (4), (5) are performed for all particles
and walls in contact.
Contact force by discrete element method
In the flows including dense solid particles, contact forces
inbetween particles and particlewall become important.
Cundall & Strack (1979) proposed a model in which a
contact force is modeled by the combination of simple
mechanical models such as a spring, a dashpot and a friction
slider as shown in Fig. 1.
The contact force is divided into the normal force Fc, and
tangential force Fc that are given by
Fen = kdn rlv, (6)
v,n =(v, n)n, (7)
and tangential directions, respectively. k, and k, are the
stiffness of the spring for each direction, respectively. n is
the unit vector in the normal and outward direction at the
surface. v, is the relative velocity vector and 7 is the
coefficient of viscous dissipation. When the following
relation is satisfied,
IF, > f IFc,l (10)
where /tis the coefficient of friction, sliding takes place and
the tangential force is replaced by
Fc, = f IFct. (11)
t is the unit vector defined by t = v/vl.
Bodyforce type IB method by Kajishima et al.
(2001)
Kajishima et al. (2001) proposed a coupling method for a
moving particlefluid system. In the method, cell size of
the calculation is smaller than particles. A system
consisted from particles and fluid flow is calculated
assuming the fluid occupies the entire flow field and the
effect of particles is expressed by introducing an additional
force which constrains the fluid velocity field to meet the
boundary condition at the particles surface. The following
particle volumeweighted velocity u is introduced.
u av +(1 a)uf
(12)
where a is the volumetric fraction of particle at a targeted
cell. It takes a =0 for a fluidonly cell, a =1 for a cell
contained in a particle and 0< a <1 for a cell across the
interface. v, is the velocity inside of the particle:
Vp = u, +rx %,p
(13)
where r represents the relative position from the center of
the particle to the particle surface. Noslip and
nopermeable conditions are imposed on its interface,
because the particle is rigid. Thus, the continuity
restriction is also assured on u in the whole domain,
Vu = 0.
Fc = kd, rv,
Vt = Vr Vn
(14)
Equation for u is introduced as
where d. and d, are the particle displacements in the normal
Paper No
Paper No
au 1
S Vp + H +fp,
at pf
H = uVu +v f V.[Vu+(Vu)T]+g.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
of continuity and momentum derived by Anderson &
Jackson (1967).
0 +V.(U,)=0 (23)
f is the kinematic viscosity of fluid. Eq. (15) is similar to
the NavierStokes equation (2) excepting for the last term.
fp is the force to constrain the predicted flow field to satisfy
the boundary condition on the particles' surface.
Timemarching calculation of Eq. (15) consists of two steps.
In the first step, we predict the velocity by using
= u" +At(Vp /p+H)
regardless of a. The superscript represents the time and At
the time increment. The predicted velocity should be
modified byfp to meet the definition of u"+'. For the cells
inside the particles (a =1), fp =(up u)/At gives u =vp.
For the cells occupied with fluid (a =0), on the other hand,
Eq. (15) is identical to Eq. (2) because fp vanishes. The
added termf, is modeled with a linear interpolation of a
Hence as the second step, the predicted velocity i is
constrained to the particle volumeweighted velocity u by
using Eq. (18).
Fluid force and torque working on a particle can be obtained
by integrating the fluid stress and torque on the particle
surface,
F,= frTndS+G,
Sp
Mf = frx(Tn)dS.
fp in Eq. (18) can be interpreted as the momentum exchange
term through the interface. In this method, fluid force
acting on the particle surface is replaced by the volume
integral of mutual interaction force, such as
F1 = p, ffdV,
VP
V
C +
p +f
Pf
(24)
Quantities such as pressure p and velocity uf are averaged
in the cell using a weight function. In the meso and
macro scales directly treated in DEMCFD approach, the
influence of viscosity is negligible and the inviscid
behavior is assumed in Eq. (24) (Tsuji et al., 1993). The
void fraction of each cell e can be defined by the number
of particles existing in the cell at each time step. The last
term in Eq. (24) denotes the effect of particles on fluid
motion which is given by
f=P Uf) (25)
where up is the particle velocity vector averaged in a cell.
When the void fraction is less than 0.8, the coefficient P is
deduced from the wellknown Ergun equation (Ergun,
1952). When c 20.8, that from Wen & Yu equation is
used (Wen & Yu, 1966).
_U' }
1 (1e) 1 p
I150 +1.75p UP
dp8 dp
(e <0.8)
3 Upu p(10) 2
4 D d
(26)
(e > 0.8)
Drag coefficient CD is obtained by using Shiller &
Naumann and Newton equations depending on the particle
Reynolds number.
=24(l+0.15Re687)/IRe (Re < 1000) (27)
S 0.43 (Re > 1000)
where
Re= u uf psEdp.p.
(28)
In DEMIB coupling method, fluid and particle motions
are obtained by solving Eqs. (4), (5), (14), (15).
DEMCFD coupling mesoscopic model by Tsuji et
al. (1993)
In DEMCFD coupling mesoscopic model (Tsuji et al.,
1993), the cell size used in calculations is larger than
particles and enough smaller than characteristics
mesoscale characteristic structures. Fluid motion is
obtained by solving the locallyphase averaged equations
Fluid force working on a particle Ff can be obtained as a
combination of fluid drag and pressure gradient forces,
(29)
Ff,= (uuf) VP}m/pp
In DEMCFD mesoscopic model, fluid and particle
motions are obtained by solving Eqs. (4), (5), (23), (24).
where
st u)+ V.(su/U)
fp =c (v fi)At.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
inflow
Figure 2: Setup for drag force calculation of
an isolated particle
Fluid force working on a single particle
In the IB method adopted in this paper, a solid particle is
represented by its volume fraction distributions. The
reproducibility of particle's surface geometry and flow field
around the particle in calculation heavily depends on the
resolution. Resolutiondependency studies for a single
particle in flows have been extensively conducted such as
fluid drag force, Saffman and Magnus lift forces and viscous
torque (Kajishima et al., 2001; Tsuji et al., 2003; Kuwagi et
al., 2009 and Yada, 2010).
In this section, the results of steady drag force working on a
single spherical particle are shown. Hereafter, the size of
computational cell used in DEMIB calculations is uniform
through the study. The setup of calculations is shown in
Fig. 2. A particle is fixed in the domain. Uniform inflow
and convective outflow conditions are imposed on the
upstream (x =0) and downstream (x =30d,) boundaries,
respectively. For side boundaries, a tractionfree condition
is used. The resolution is varied as d/Ax = 4, 8, 16, 32 and
results are compared with following empirical equations;
Figure 3: Drag force of a single particle
White (1974):
24 6
CD = + +0.4,
Rep 1+ Rep
Shiller and Naumann (1933):
CD 24 (1+0.15Re.687)
Re
(Rep <800),
and
Morsi and Alexander (1972):
CD = C + C + C2
Re, Re2
(30)
(31)
(32)
where Co, cl and c2 are the model parameters.
Fig. 3 shows the relation between the particle Reynolds
number and drag coefficient. We can confirm that a
general tendency of all resolution cases agrees with the
empirical equations. Especially for Re, =10, 50, 100, a
resolutiondependency is not apparent and the results agree
well with empirical equations quantitatively. In the regime
100 < Rep < 1000, on the contrary, a dependency on the
resolution becomes apparent. From Fig. 3, it is confirmed
that a fluid drag force is estimated larger as the resolution
becomes lower. It is known that wake structures formed
behind the particle give large influences on the fluid drag
force working on a particle. Wake structures becomes
smaller as the particle Reynolds number becomes larger.
Hence, a higher resolution is required in high particle
Reynolds number cases. The empirical equation by White
(1974) shows larger values comparing to Shiller and
Naumann (1933) and Morsi and Alexandar (1972) in this
region. Excepting for d/Ax = 4 case, the results are
located inbetween these empirical equations. When Rep
=1000, the result of d/Ax = 4 case is 76 and 33 % larger
comparing to Morsi and Alexandar (1972) and White (1974)
equations, respectively.
Fluid force working on paired two particles
In dense flows, many particles come close each other and
have contacts. In such situation, hydrodynamic
interactions due to the flows through the narrow gaps
existing inbetween particles become important. In this
section, fluid drag forces working on paired two particles
are investigated. Particles are arranged in tandem and side
by side normal to the main stream. In addition to the
resolution of calculation, an interparticle distance is varied
and its influence is investigated. In tandem arrangements,
the same domain size with the single particle case is used.
In sidebyside arrangements, a crosssectional size of
calculation domain is enlarged to keep the distance between
a particle surface and side boundary to 5.75dp through the
cases.
Fig. 4 shows the relation between normalized interparticle
distance and drag coefficient of the trailing particle in
tandem arrangements. The particle Reynolds number in
gasfluidized beds can be O(103). Investigations in paired
Paper No
White (1991)
Schiller and Nauman (1933)
orsi and Alexander (1972)
x d/Ax 4
p
A d/Ax 8
p
d/Ax 16
p
d* /Ax 32
`VX
Paper No
0.7
0.6
0.5
0.4
0.3,
0.2
0.1 k
X d/Ax
P
p
A d/Ax
p
V d/Ax
P
* d/Ax
p
0 1 2 3
l/d
Figure 4: Drag force of the trailing particle for two
particles aligned streamwise
1.25
X d/Ax 4
p
A dAx 8
P
v d/Ax 16
* d/Ax 32
P
0 1 2 3
l/d
p
Figure 5: Drag force of particles aligned side by side
arrangements also should be conducted up to 0 (103),
however, the number of experimental results is restricted in
high Reynolds number regions and the results of Rep= 106
case is discussed in this section. Interparticle distance 1 is
defined as a distance between particle surfaces. Drag
coefficient is normalized by that of a single particle under
the same Reynolds number condition. Experimental
results by Zhu et al. (1994) are also shown. A general
tendency of all resolution cases agrees with the experimental
result and the drag forces decrease when l/d, becomes
smaller. When two particles are touching (l/d,=O), the drag
coefficient becomes 23% of single particle in the case
d/Ax=32. In the case of /d, = 3, the hydrodynamic
interaction between two particle is still clearly observed.
When the interparticle distance becomes larger than 2, a
resolutiondependency is not observed clearly. Excepting
for the case d/Ax=4, a resolutiondependency is also not
clear in the region /d, < 1.
Fig. 5 shows the relation between interparticle distance and
drag coefficient working on the particles placed side by side
normal to the main stream. The particle Reynolds number
is set to 102. An experimental result by Chen & Lu (1999)
is also shown. A general tendency of all resolution cases
also agrees with the experimental result in sidebyside
arrangements. Comparing to the tandem arrangements, the
effect of adjacent particle is not so intense even if the two
particles are almost touching (l/d,=0.3). In is confirmed
that the results of d,/Ax=16, 32 cases are almost converged
in the range l/d, >1. Similar to the single particle and two
particles in tandem arrangement, the result of d/Ax=4 case
is largely deviated from others.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table 1: Conditions for a packed bed
Case 1 2
Domain size
(L i siL 60 x 300 x 60
(Lx, L,, L,) [mm]
Number of mesh [] 50 x 750x 150 300 x 1500 x 300
Particle diameter:
[mm] 3.20
d, [mm]
Number [] 455
dAx [] 8 16
Void fraction: e [] 0.4261
Viscosity: u [Pa s] 1.81x105
Density: pf[kg/m3] 1.205
Reynolds number:
1, 5, 10, 50, 100, 500
Re, (=uodJv)
Fluid force working on particles in a
randomlypacked bed
The fluid drag force working on particles in
randomlypacked bed is verified. Conditions of
calculation are shown in Table 1. Due to the
computational time limitation, only results of d/Ax=8 and
16 cases are compared. A packed bed including
randomlydistributed particles is obtained as follows.
(1) An arbitrary random velocity is given to all particles
arranged regularly at first.
(2) By imposing the gravity force, particles are settled
down from an appropriate height. During the settling,
routines for the fluid calculation are skipped and only
DEM calculation is performed.
(3) Due to the artificial random velocity, particles have
contacts with other particles. Potential and kinetic
energy particles obtained are dissipated gradually due
to the dashpot. After the steady state is reached,
positions of particles are fixed.
For inflow and outflow boundaries, an uniform inflow and
convective outflow condition are imposed, respectively.
Periodic boundary condition is used for horizontal (x, z)
directions. To avoid the effect of artificial forced inflow,
all particles are kept in 60 mm higher position from the
inflow boundary. The fluid drag forces working on all
particles in the domain are averaged and results are
compared with the drag force equation obtained from Ergun
equation.
Ergun (1952) shows that the pressure drop in a packed bed
consisted from spherical particles having the same diameter
dp is expressed as
Ap 1 150(1p +1.75pu0
L d, dp
(33)
where L, uo are the dimension and superficial velocity of the
packed bed, respectively. The fluid drag force working on
a particle, FD, is expressed as
S(34)
F, = CDA 2p'fof (34)
Paper No
p
104 V v rd/Ax16
102
101
10
100 10' 102 103
Figure 6: Drag force working on a particle in
randomlypacked bed (e =0.427).
where CD is the drag coefficient, A is the crosssectional
area of the particle. When the pressure drop is balanced
with the fluid drag force working on N particles existing in
the domain V,
FD xN= xV.
L
(35)
From Eqs. (33), (34), (35), an empirical equation for the
drag force working on a particle in a packed bed is obtained.
CD = 150,u 1 +1.75 1
3 d,u ' C3
Fig. 6 shows the relation between the particle Reynolds
number and the drag coefficient defined by Eq. (36). The
averaged void fraction of the packed bed we obtained by
the settling method is 0.4261. In the region 1_ Rep, 50,
a resolutiondependency is not apparent and the behavior
of both cases agrees with Eq. (35) in the qualitative sense.
Quantitatively, the values obtained by calculations are 23
times larger from Ergun equation. In the region Rep >
100, the results of d/Ax=8 become relatively large
comparing to the Eq. (35).
2D Gasfluidized bed under spouting condition
In this section, a calculation of flow inside of a fluidized
bed is performed by using DEMIB method. For
comparisons, a calculation by using DEMCFD
mesoscopic model is also performed. The conditions of
calculations are shown in Table 2. 3D calculations are
performed for socalled 2D fluidized bed in which the
depth of the bed is restricted comparing to other
dimensions. Initially, particles are packed randomly in
the bed using the same technique shown in the previous
section. The same initial condition is used for both
calculations. Calculations are performed under a
spouting condition where a gas jet is injected only from the
bottom center of the bed.
Noslip condition is used for side and bottom walls in the
DEMIB approach. On the contrary, slip condition is
used in the DEMCFD mesoscopic model approach
following Tsuji et al. (1993). For the upper boundary, the
convective outflow condition is used in both approaches.
In DEMIB calculation, a numerical code is parallelized by
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Table 2: Conditions for a fluidized bed under
a spouting condition
Container
Size of container (Lx, L, L,) [mm] 96 x 160 x 16
Number of mesh (DEMIB) [] 240 x 400 x 40
Number of mesh (DEMCFD) [] 12 x 16 x 2
Inflow crosssection size [mm] 16 x 16
Particle
Particle diameter: d, [mm] 3.20
Number [] 3104
Density: p, [kg/m3] 910
Initial bed height [mm] 64
Gas
Superficial velocity: uo [m/s] 1.0
Minimum fluidization velocity: 1.037
Umf [m/s]
Viscosity: p [Pa s] 1.81x105
Kinetic viscosity: v [m2/s]
Density: p [kg/m3]
Gravitational constant: g [m/s2]
DEM
Normal spring constant: k, [N/s]
Tangential spring constant: kt [N/s]
Coefficient of restitution: e []
Coefficient of friction
(particle particle): p, []
Coefficient of friction
(particlewall): uli, []
Time
Time increment (IBDEM): At [s]
Time increment (IBCFD): At [s]
Total time: T [sl
1.45 105
1.205
9.8
800
200
0.43
0.88
Ixl06
2.0 x 105
1.8
using a standard 1D domain decomposition technique.
For the data transfers between different domains, MPI
library is used. All calculations are performed using a
cluster computer consisted from 128 Intel Xeon Woodcrest
processors (3.0 GHz) in Cybermedia center, Osaka
University.
As we observed already, an accuracy of IB calculation
heavily depends on the resolution. In dense case, the
resolution requirement is depending on the interparticle
distance (void fraction) in addition to the particle Reynolds
number. Void fraction and particle Reynolds number are
expected to vary depending on flow conditions and
positions in a gasfluidized bed under a spouting condition.
In the DEMIB calculation demonstrated in this section, a
resolution is set to d/Ax= 8. Under the conditions shown
in Table 2, the particle Reynolds number is expected to
exceed 103 in the vicinity of gas inlet section and the
averaged void fraction at the initial condition is less than
0.45. We admit that d/Ax 8 is not sufficient and higher
resolution is needed to enable quantitative observations
and predictions. Even in case of d,/Ax=8, however, it
requires more than 40 days to obtain 1.8 s results by using
16 cores. We restrict our discussions in this section to
qualitative observations using d,/Ax=8 only.
Fig. 7 shows the temporal development of the bed for 0.01
to 0.50 s. The results are obtained by using DEMIB
Paper No
method. The position and velocity of the particles,
isosurface of gas velocity and streamwise gas velocity
component at the bed center (z = 8 mm) and a plane near
the front wall (z = 0.2 mm) are shown.
From Fig. 7 (a), it is confirmed that a bubble is formed just
above the gas inlet and particles are start to fluidize. The
investigations are performed under a spouting condition
and we can confirm that fluidized and unfluidized regions
are clearly separated. As observed in previous
experimental studies, particles in the region just above the
gas inlet is wellfluidized and it is not in the region close to
the bottom covers of the bed. As Eq. (36) shows, the
fluid drag force working on a particle in a packed bed
changes drastically depending on the void fraction. In
actual gasfluidized beds, the void fraction is not uniformly
distributed and can change drastically depending on the
time. Fig. 7 (b) shows the isosurface of gas velocity
vector. The same value with inflow gas velocity is used
for the visualization. Fig. 8 show the fluid velocity
vectors and particle's volume fraction distributions. Only
the results at t = 0.01 and 0.50 s are shown. At the initial
stage (t = 0.01 s), void fraction distributions of the bed is
relatively uniform excepting for the regions near the side
walls. Gas flows are trying to find a route in which a
resistance due to particles is minimum. As can be seen
from Figs. 7 (b), (c) and 8 (a), gas flow paths spread over
the whole bed and result in a very complex structure. In a
plane near the front wall (Figs. 7 (d), 8 (b)), gas velocity
becomes larger comparing to the center because void
fraction is relatively high. Once the bubble is formed
(0.05 to 0.5 s), gas flows do not spread over the bed and
follow the bubble because the void fraction is large inside
of the bubble and fluid resistance is small.
As shown in the previous section, when the resolution is
not enough, the fluid drag force working on the particles in
packed bed is overestimated especially in a high particle
Reynolds number region. It means that particles start to
fluidize under a smaller superficial velocity. In the
calculations performed here, the superficial velocity of the
calculations is set to 1.0 m/s which is smaller than the
minimum fluidization velocity 1.037 m/s. The result of
DEMCFD mesoscopic model under the same conditions is
shown in Fig. 9. Particles in the region just above the gas
inlet are raised and high void fraction region is formed,
however, fluid drag force working on particles is estimated
from Ergun and Wen & Yu equations (Eq. (26)) in
DEMCFD mesoscopic model and it does not start to
fluidize and reach a steady state as we expected. In case
of the DEMC mesoscopic model, it requires only 3 hours
to obtain the results by use a desktop PC.
Conclusions
Toward an establishment of reliable direct simulation
technique for the flows including dense solid particles such
as that in a gasfluidized bed, DEMIB coupling method
was presented in this paper. In addition to a standard
steady drag force problem, resolutiondependency studies
were conducted for the drag force problems in that
hydrodynamic interactions due to gas flows existing
inbetween particles are important. Besides, by using
DEMIB coupling method, a numerical simulation of a 2D
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
gasfluidized bed was performed. For comparison, a
simulation by DEMCFD mesoscopic model under the same
condition was also performed.
If the resolution is not sufficient, fluid drag force is
overestimated in IB calculations. The trend is common to
the problems investigated in this paper: drag forces working
on a single particle, paired particles in tandem and
sidebyside arrangements and particles in randomlypacked
bed. This becomes apparent when the particle Reynolds
number is high and interparticle distance is small.
Due to the computational time restriction, the resolution
used in the calculation of a gasfluidized bed under a
spouting condition is still insufficient for quantitative
observations and predictions of flows. Actually, it starts to
fluidize with the superficial velocity condition which is
smaller than the minimum fluidization velocity. By using
the direct simulation as we demonstrated in this paper, it is
possible to observe the microscopic flows in a particlelevel.
It helps essential understandings of the behavior of flows
including dense solid particles. We expect that the results
obtained by DEMIB method will approach to experimental
results when the resolution of calculation becomes higher.
In gasfluidized beds, void fraction changes drastically and
has distributions in the bed. Resolutionrequirement in
most dense region will be a criterion to enable quantitative
observations and predictions.
Acknowledgements
We would like to show our acknowledgements to New
Energy and Industrial. Technology Development
Organization (NEDO) and GrandinAid for Young
Scientists (B), Japan Society for the Promotion of
Science. We also would like to show our
acknowledgement to Cybermedia center, Osaka University.
References
Anderson, T.B. & Jackson, R. A fluid mehcanical
description of fluidized bed, I & EC Fundamentals, Vol. 6,
527539 (1967)
Chen, R.C. & Lu, Y.N. The flow characteristics of an
interactive particle at low Reynolds numbers. Int. J.
Multiphase Flow, Vol. 25, 16451655 (1999)
Cundall, PA. & Strack, O.D.L. Discrete numerical model
for granular assemblies, Geotechnique, Vol. 29 (1), 4765
(1979)
Ergun, S. Fluid flow through packed columns, Chem. Eng.
Progress, Vol. 48 (2), 8994 (1952)
Kajishima, T. Takiguchi, S. Hamasaki, H. & Miyake, Y.
Turbulence structure of particleladen flow in a vertical
plane channel due to vortex shedding. JSME Int. B Vol. 244,
526535 (2001)
Kuwagi, K. Shimoyama, Y Hirano H. & Takami T.
Estimation of lift force and Viscous torque with Immersed
boundary (IB) method. Proc. 15th SCEJ Sympo.
Paper No
Fluidization & Particle processing, 8386 (2009) (In
Japanese)
Morsi, S.A. & Alexander, A.J. An investigation of particle
trajectories in twophase flow systems, J. Fluid Mech., Vol.
55, 193208 (1972)
Pan, T.W. Joseph D.D. BAI, R. Glowinski, R. and Sarin, V
Fluidization of 1204 spheres: simulation and experiment. J.
Fluid Mech., Vol. 451, 169191 (2002)
Schiller, L. & Naumann, A. Uber die grundlegenden
berechungen bei der schwerkraftaufbereitung, Vereines
Deutscher Ingenieure, Vol. 77, 318320 (1933)
Tsuji, T. Narutomi, R. Yokomine, T. Ebara, S. & Shimizu A.
Unsteady threedimensional simulation of interactions
between tow and two particles, Int. J. Multiphase flow, Vol.
29 (9), 14311450 (2003)
Tsuji, Y Kawaguchi, T. & Tanaka, T. Discrete Particle
Simulation of TwoDimensional Fluidized Bed. Powder
Technol., Vol. 77, 7987 (1993)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Uhlmann, M. An immersed boundary method with direct
forcing for the simulation of particulate flows. J. Comp.
Phys., Vol. 209, 448476 (2005)
van der Hoef, M.A. van Sint Annaland, M. Deen, N.G. &
Kuipers, J.A.M. Numerical simulation of dense gassolid
fluidized beds: a multiscale modeling strategy. Annu. Rev.
Fluid Mech.,Vol. 40, 4770 (2008)
Wen, C.Y. & Yu, Y.H. Mechanics of fluidization. Chem. Eng.
Prog. Sympo. Ser., Vol.62 (62), 100111 (1966)
White, FM. Viscous Fluid Flow, McGrawHill, New York
(1974)
Yada, H. Direct numerical simulation of flows containing
dense solid particles. Master thesis, Osaka University
(2010) (In Japanese)
Zhu, C. Liang, S.C. & Fan, L.S. Particle wake effects on
the drag force of an interactive particle. Int. J. Multiphase
Flow, Vol. 20, 117129 (1994)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.01 s
0.05 s
0.15 s
0.35 s
0.50 s
(a) Particles' position
and velocity
Figure 7:
Fluid Velocity
Io n 5 1n I52n
Fluid Velocity
[r/si
105 0 5 101520
Ims
Fluid Velocity
Fm/si
Fluid Velocity
LlLLJJU521L
(b) Isosurface of gas (c) Streamwise gas
velocity (lu=6.0m/s) velocity at z = 8 mm
Development of flow field by DEMIB method
Fluid Velocity
Fluid Velocity I
Fluid Velocity
Fluid Velocity
Fluid Velocity l
(d) Streamwise gas
velocity at z = 0.2 mm
Paper No
Paper No
0.01 s
0.50 s
(a) Center (z = 8 mm) (b) In a plane close to the front wall (z = 0.2 mm)
Figure 8: Instantaneous fluid velocity distributions obtained by DEMIB at z = 8 mm and t = 0.50 s.
Figure 9: Particles' position and velocity
obtained by DEMCFD mesoscopic model
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
ii:
t
t
... 1 1' 1.; "'
; ,, '
' *'. ..',
i
; /.
[**, *' ,'
