- Email: [email protected]
Effect of a wall on flow with dense particles
Advanced Powder Technology 24 (2013) 565–574
Contents lists available at SciVerse ScienceDirect
Advanced Powder Technology journal homepage: www.elsevier.com/locate/apt
Translated paper; SPTJ Best Paper Award 2011
Effect of a wall on flow with dense particles q Takuya Tsuji ⇑, Eiji Narita, Toshitsugu Tanaka Department of Mechanical Engineering, Graduate School of Engineering, 2-1 Yamada-Oka, Suita, Osaka 565-0871, Japan
a r t i c l e
i n f o
Article history: Received 6 November 2012 Accepted 12 November 2012 Available online 30 November 2012 Keywords: Wall effect Gas–solid flow Fluidized bed Direct numerical simulation Immersed boundary method
a b s t r a c t The behavior of dense gas–solid flows in engineering applications such as fluidized beds and pneumatic conveyers is highly complex and a reliable numerical model is required. Such flows are usually within solid walls that considerably affect the flow fields, and it is important to correctly include this effect in numerical models to improve their prediction capability. The observation of microscopic flows near walls can enhance our understanding of the flow behavior and assist in improving models. In this study, direct simulations are performed to investigate the effect of a wall on flow fields at a microscopic level. The effects of the bulk void fraction, particle Reynolds number, and particle diameter are investigated. The prediction performances of existing correlation equations usually used in mesoscopic model calculations are also investigated. It is found that the Ergun and Beetstra equations produce large discrepancies in the region within a distance equal to the particle diameter from the wall. Ó 2012 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights reserved.
1. Introduction Dense particles coexist with a gas in industrial equipment such as fluidized beds and pneumatic conveying systems. The flow behavior of the gas–particle mixture is significantly influenced by the interactions between the particles, between the particles and the wall and between the particles and the gas. The spontaneous formation of mesoscopic heterogeneous structures inside the equipment, such as bubbles that are far larger than the particles, makes the flow more complex. A better understanding of the flow behavior and design optimization is therefore desired. This has, however, not been successfully achieved owing to observation difficulties in experiments. Against this background, the development of numerical models has been extensively pursued until the present time. The majority of these models are mesoscopic in the sense that the size of the computational cells of the fluid motion is small enough to resolve characteristic mesoscopic structures but too large to directly capture the microscopic phenomena on the particle level [1–3]. The models have been applied to several flow problems, and the reproduction of size and frequency with which bubbles appear in fluidized beds have been successfully demonstrated. In practical engineering situations, solid walls that significantly influence the behavior of both phases are very common. It is important to correctly express the influence of the solid walls whenever a numerical simulation is performed. In mesoscopic model calcula-
q Japanese version published in JSPTJ, Vol. 48 No. 12 (2011); English version for APT received on November 6, 2012. ⇑ Corresponding author. Tel./fax: +81 6 6879 7317. E-mail address: [email protected] (T. Tsuji).
tions, the momentum exchange between solid particles and the gas in a computational cell is generally determined by using empirical equations that are a function of the void fraction and the particle Reynolds number. In the majority of existing models, empirical equations such as those of Ergun [4], Di Felice [5], and Beetstra et al. [6] are used without any modifications, including in computational cells near a solid wall. The existence of a wall directly exerts viscous friction on nearby fluids and affects the arrangement of nearby particles. As a result, the anisotropy of the flow in the vicinity of solid walls is increased. In the aforementioned empirical equations, the effect of the wall is not considered and it is questionable to use these equations for computational cells near solid walls. Moreover, the problem is one that has been recognized by researchers since early times and was mentioned by Beetstra et al. [6]. When microscopic flows are directly considered, fluid motions should obey the Navier–Stokes and continuity equations, and it is appropriate to use the no-slip boundary condition on the surface of solid walls. However, in mesoscopic model calculations based on locally phase averaged equations, microscopic flows occurring on a subgrid scale are averaged and it is questionable to use the no-slip boundary condition on the surface of solid walls on the mesoscopic cell scale. Currently, different boundary conditions are used depending on a study, and a general consensus is yet to be achieved. For further improvement of mesoscopic models, it is important to adopt a proper drag correlation in computational cells near a solid wall and a proper boundary condition on the wall. This cannot be accomplished without a thorough understanding of the behavior of microscopic flows near such walls, which never be captured by the mesoscopic model calculations. Although the effects of a solid wall on a flow that contains dense particles have been investigated [7–9], it is still difficult to directly observe microscopic flows
0921-8831/$ - see front matter Ó 2012 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights reserved. http://dx.doi.org/10.1016/j.apt.2012.11.006
566
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
Nomenclature a dp FD FL fp Gp g I Ip L Lx, y, z Mp mp Np n nx, y, z p Rep r Sp t U Up u
particle radius (m) particle diameter (m) fluid drag force (N) fluid lift force (N) interaction force (m/s2) external force acting on a particle (N) acceleration due to gravity (m/s2) unit tensor (–) moment of inertia (kg m2) streamwise domain size (m) domain size in x, y, and z directions (m) external moment acting on a particle (N m) particle mass (kg) number of particles in the computational domain (–) vector normal to the surface of a particle (–) number of grid points in x, y, and z directions (–) fluid pressure (Pa) particle Reynolds number (–) relative vector from the center of a particle (–) surface area of a particle (m2) time (s) inflow velocity (m/s) translating velocity of a particle (m/s) coupling velocity weighted by particle volume fraction (m/s)
on a particle scale by experimental methods. In the present study, the body-force-type immersed boundary method (IBM) is used to conduct microscopic numerical simulations of flows that contain dense solid particles near a solid wall, and the effects of the wall on the flow behavior are investigated in detail. In addition, the data obtained from the microscopic simulations are used to discuss the validity of the empirical drag equations popularly used in mesoscopic model calculations. In this paper, which is the first report, we restrict our discussions to the case of a fixed particle. 2. Calculation methods 2.1. Governing equations It is assumed that the fluid is incompressible and Newtonian and the particles are rigid spheres. The governing equations of the fluid flow are the continuity equation and Navier–Stokes equation, respectively, expressed below:
r uf ¼ 0
qf
Duf ¼ r s þ qf g Dt
s ¼ pI þ lf ½ruf þ ðruf Þ
dðmp U p Þ ¼ dt dðI p Xp Þ ¼ dt
s n dS þ Gp
Sp
r ðs nÞ dS þ M p
e qf lf mf s Xp
w hi h iyz
where Up and Xp respectively represent the translational and rotational velocities of the particle, mp is the mass of the particle, Ip is the moment of inertia given by Ip = (2/5)a2mpI for a particle with radius a, Gp and Mp are respectively the external force and moment, and r is the relative position vector from the center of a particle. 2.2. Body-force-type IBM [10] Kajishima et al. developed the body-force-type IBM for particle– fluid systems. Here, the assumption of a uniform Cartesian grid is used to simplify the explanation. The presence of particles is represented by the solid volume fraction of particles in each computational cell a. Kajishima et al. [10] introduced the following coupling velocity weighted by the solid volume fraction of particles:
u ¼ aup þ ð1 aÞuf
ð6Þ
ð2Þ
ru¼0
ð3Þ
ð4Þ
Sp
Z
Dp DpBeetstra DpDF DpErgun Dt Dx
ð1Þ
The particles are tracked individually in a Lagrangian manner by using the equations of translational and rotational motions:
Z
a
superficial velocity (m/s) fluid velocity (m/s) particle velocity (m/s) streamwise fluid velocity (m/s) velocity predicted as a fluid (m/s) volume of computational domain (m3) volume of a particle (m3) volume fraction of particle in each computational cell (–) pressure drop (Pa) pressure drop predicted by Beetstra equation (Pa) pressure difference imposed as a driving force (Pa) pressure drop predicted by Ergun equation (Pa) temporal increment (s) computational cell size (m) void fraction in whole region (–) density of fluid (kg/m3) viscosity of fluid (Pa s) kinetic viscosity of fluid (m2/s) stress tensor (Pa) rotational velocity of a particle (rad/s) nondimensional pressure drop (–) spatial average in whole region spatial average in y–z plane temporal average
where up is the particle velocity given by up = Up + r Xp. The particles are solid and no-slip and no-permeability conditions are imposed on their interfaces. Hence, the continuity restriction also applies to u:
where uf and qf are respectively the fluid velocity and density, g is the acceleration due to gravity, p is the pressure, lf is the fluid viscosity, and s is the stress tensor given by T
u0 uf up uy ~ u V Vp
ð5Þ
ð7Þ
The following equation for u is introduced:
@u 1 ¼ rp þ H þ f p @t qf
ð8Þ
where
H ¼ uru þ mf r ½ru þ ðruÞT þ g
ð9Þ
and mf is the kinematic viscosity of the fluid. In Eq. (8), fp is the interaction force term that enforces all the predicted velocity fields (including inside the particles) to the particle–fluid coupling velocity field. This is explained in terms of the time-marching procedure. If Eq. (8) marches in time with the explicit Euler method, it would be
567
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
unþ1 ¼ un þ Dt
1
qf
!n
rp þ H þ f p
ð10Þ
Rep (=Udp/vf)
where the superscript represents time and Dt is the time increment. We predict the velocity as the fluid phase once regardless of the cell is occupied by fluid (a = 0), particle (a = 1) and both (0 < a < 1).
~ ¼ un þ Dt u
1
qf
Table 1 Comparison of predicted drag coefficient with Zeng et al. [13].
10 200
CD
Difference (%)
Present
Zeng et al. [13]
4.988 0.8719
4.721 0.8156
5.655 6.902
!n
rp þ H
ð11Þ
The predicted velocity should be modified by fp to meet the definition of u in Eq. (6). For the cell inside the particle (a = 1), ~ Þ=Dt gives un+1 = up (inside the particle). However, for f p ¼ ðup u the cell occupied by fluid (a = 0), Eq. (8) is identical to Eq. (2) because fp vanishes. Hence, the added term fp is modeled with the linear interpolation of a as
~ Þ=Dt f p ¼ aðup u
ð12Þ
Table 2 Calculation conditions for packed bed without wall. Computational domain: Lx Ly Lz (mm) Number of grid points: nx ny nz (–) Particle diameter: dp (mm) Number of particles: Np (–) Voidage: e = h1 ai (–) Fluid density: qf (kg/m3) Kinetic viscosity: mf (Pa s)
6dp 10dp 6dp 96 160 96 1.0 360 0.4764 1.205 1.448 105
This force is used to modify the solution of Eq. (11) to u (Eq. (6)). Additionally, the term is used for determining the momentum exchange between two phases through the interface. The fluid force acting on a particle surface (Eqs. (4) and (5)) is obtained by integrating the mutual interaction force fp with respect to the volume:
dðmp U p Þ ¼ qf dt dðI p Xp Þ ¼ qf dt
Z
f p dV þ Gp
ð13Þ
Vp
Z
r f p dV þ M p
ð14Þ
Vp
The present method is applicable to moving particles. In all calculations in this study, the particles are fixed. The numerical calculations are performed using the finite difference method based on the fractional step method. A staggered uniform Cartesian grid is used. Fourth-order central difference and second-order Adams– Bashforth schemes are respectively used for the discretization of the spatial and temporal derivatives. The method proposed by Tsuji et al. [11] is used to obtain the solid volume fraction in each computational cell. 3. Verification of calculation accuracy In IBM calculations, solid particles are represented on a Cartesian grid and the accuracy of the calculations varies, depending on the particle resolution, which is the number of computational cells per particle diameter. Although improved accuracy can be expected with higher particle resolutions, excessive resolutions significantly increase computational load and should be avoided. In advance to the main calculations, the relationship between the particle resolution and the calculation accuracy should be investigated. A number of investigations have been carried out, including those on the fluid drag force acting on single and paired particles in a uniform flow, the Saffman lift force, the Magnus lift force acting on a particle, and the viscous torque acting on a rotating particle [11,12]. Our points of focus in the present study are the flow containing dense particles near a solid wall, verifications study are performed for the fluid drag force acting on a single particle moving steadily near a solid wall and the pressure drop of a fluid in a packed bed without solid walls. 3.1. Fluid drag force acting on a single particle near a solid wall The accuracy of present calculations is verified in this section. A single particle traveling through an acquiescent fluid near a solid wall with a constant velocity U is considered. The distance between the particle surface and the solid wall is fixed at 1/2dp. Based on the
Fig. 1. Comparison with Ergun and Beetstra equations.
results of previous resolution-dependency studies, the number of computational cells per particle diameter is set to dp/Dx = 16. The results are compared with those that Zeng et al. obtained using the spectral element method [13]. The particle Reynolds number, defined as Rep = Udp/mf, is set to 10 and 200, where dp is the particle diameter and mf is the dynamic viscosity of the fluid. The results are shown in Table 1. For both values of Rep, the agreement between the results of our calculation and those of Zeng et al.’s are fairly good. The values of our results are slightly higher but the differences are within 7%. The same particle resolution dp/Dx = 16 is hereafter used for our study. 3.2. Pressure drop in packed beds Calculations of the pressure drop in packed beds without solid walls are performed and the results are compared with those obtained with the empirical equations by Ergun [4] and Beetstra et al. [6]. Ergun [4] expressed the pressure drop in packed beds as
ð1 eÞlf u0 Dp 1 e 2 150 þ 1:75 q u ¼ f 0 dp e3 dp L
ð15Þ
where Dp is the pressure difference, L is the bed height, e is the void fraction, dp is the particle diameter, lf is the viscosity of the fluid, u0 is the superficial velocity, and qf is the density of the fluid. Meanwhile, Beetstra et al. [6] determined the fluid drag force FD acting on a single particle in packed beds from the results of the lattice Boltzmann simulation:
568
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574 Table 3 Calculation conditions. Computational domain: Lx Ly Lz (–) Number of grid points: nx ny nz (–) Particle diameter: dp (mm) Number of particles: Np (–) Reynolds number: Rep ¼ e< uy >dp =v f (–) dp/Dx (–) Voidage: e = h1 ai (–) Fluid density: qf (kg/m3) Kinetic viscosity: mf (Pa s)
6dp 10dp 6dp 96 160 96
8.25dp 13.5dp 8.25dp 132 216 132 1.0 360 0.1–54.6 16 0.802 1.205 1.448 105
0.503
y
10d p [mm]
13.5dp [mm]
x
8.25dp [mm]
6dp [mm] (a) ε = 0.503
(b) ε = 0.802 Fig. 2. Calculation domain.
F D ¼ 3plf dp u0 Fðe; Rep Þ Fðe; Rep Þ ¼
10ð1 eÞ
e
3
ð16Þ
pffiffiffiffiffiffiffiffiffiffiffi þ eð1 þ 1:5 1 eÞ
0:413Rep þ 24e3
! e1 þ 3eð1 eÞ þ 8:4Re0:343 p 1 þ 103ð1eÞ Rep2e2:5
w¼ ð17Þ
Let us assume that Np particles are present in a control volume V with height L. Considering that the pressure drop in the control volume is due only to the particles in the box, the total fluid drag force acting on the Np particles should be balanced by the pressure drop. Hence,
F D Np ¼
Dp V L
ð18Þ
The void fraction e in the box can be related to the volume of each particle, Vp, by
1e¼
Np V p V
ð19Þ
Hence, the pressure drop of the Beetstra equation is obtained with
DP Vp ¼ FD L 1e
Hereafter in this paper, the pressure drop is discussed in terms of the following non-dimensional form:
ð20Þ
DP dp L qf u20
ð21Þ
The calculation conditions are listed in Table 2. In the table, hi shows the spatial averaging over the whole calculation domain, including the inside of solid particles, and the bulk void fraction e can be calculated from the solid volume fraction at each computational cell using
e ¼ h1 ai
ð22Þ
In this section, the periodic boundary condition is used in all directions. To obtain a random arrangement of the particles, the position of each particle is fixed after tracking the motion of all the particles for a sufficient time imposing an initial random velocity and an artificial potential force between particles. During the packing, fluid calculations are not performed. The constant averaged pressure gradients in the y direction are the driving forces and the calculations are performed by varying the particle Reynolds number. The results are not included here; however, the sufficiency of the domain size was confirmed by preliminary calculations. The working fluid is air at room temperature. Fig. 1 shows the results, where Rep, the
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
569
4. Effects of a solid wall 4.1. Calculation conditions In this section, the influence of a wall on a flow field is investigated. The calculation conditions are listed in Table 3. As in the part discussed in the previous section, random packing of the particles near the solid wall at x = 0 is obtained by tracking all the particles without fluid calculations. An artificial potential force is also applied between the particles and the wall, in addition to between the particles. During particle packing, the periodic boundary condition is used in the y and z directions and solid walls are placed at x = 0 and x = Lx + dp, respectively. Flows near a single solid wall at x = 0 are of interest in the present study; however, the effect of the wall on the particle arrangement at x = Lx + dp would remain with above procedure. To minimize the effect, a sub-domain within Lx < x 6 Lx + dp is eliminated from the computations and the slip boundary condition is imposed instead at x = Lx. The flows are maintained by imposing a constant averaged pressure gradient in the y direction. After a steady state is attained, the flow field is observed. In this study, different bulk void fractions are achieved by varying the calculation domain size while maintaining the number of particles included in the domain. Snapshots of the particle arrangements for e = 0.503 and e = 0.802 are shown in Fig. 2a and b, respectively. Based on the results of the verification process discussed in the previous section, the investigations are restricted to the case of Rep < 60. Roughly speaking, this condition corresponds to that of free bubbling in beds with 500 lm particles. 4.2. Observation of flow field near a solid wall Fig. 3a shows the void fraction distribution h1 aiyz at e = 0.503. The result is averaged in the y–z planes. It is observed that the void fraction has an oscillatory profile and progressively approaches the
Fig. 3. Voidage distribution near the wall.
particle Reynolds number based on a superficial velocity, is defined by
Rep ¼
ehuy idp mf
ð23Þ
where — indicates temporal averaging. A fairly good agreement between the results of our calculations and those obtained with empirical equations is observed; the errors are less than 20% in the region of Rep < 60. The difference increases as Rep increases, which is a result of the resolution limitations of our calculations; we cannot capture the velocity gradient of the fluid on the surface of the particles. Improved accuracy is expected upon increasing the resolution of the simulations, although this would also lead to an increase in computational cost. In our study, we restrict our investigations to the region of Rep < 60, where the accuracy of the calculations is guaranteed. While the verification process of this section is conducted by comparing our calculation results with those obtained with empirical equations, it must be noted that there is considerable difference between the results by the Ergun and Beetstra equations, especially in the region of Rep > 100. Further investigations into these empirical equations may also be required.
Fig. 4. Instantaneous velocity vector distribution in an x–y plane (e = 0.503, Rep = 54.6).
570
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
x–z planes, respectively. It can be observed that the streamwise velocity component becomes relatively large in the regions with higher void fractions. This is especially apparent in the near-wall region where x/dp = 0 in Fig. 5b. Fig. 6 shows the distributions of the averaged streamwise velocity huy iyz for different particle Reynolds numbers. The void fraction distribution is also shown. Different particle Reynolds numbers are obtained by changing the value of the imposed averaged pressure gradient. The averaged streamwise velocity is normalized by the bulk streamwise velocity in the entire domain < uy >. Except in the vicinity of the wall where the influence of the no-slip boundary condition on the wall surface exists, the streamwise velocity is large whenever the void fraction is large. The results correspond to Fig. 4 and show that the flows preferentially flow through large gaps between the particles and between the particles and the wall. A dependency on Rep is only observed in the vicinity of the wall, where the streamwise velocity becomes relatively large under higher Rep conditions. This is not observed in the far-wall regions, i.e., when x/dp > 0.7 in Fig. 6a and x/dp > 1.6 in Fig. 6b. In the body-force type IBM, the fluid force acting on each particle can be easily obtained with Eq. (13). Fig. 7 shows the relationship between the distance of the particle’s center from the wall and the time-averaged drag and lift forces acting on each particle when
Fig. 5. Streamwise velocity component distributions (e = 0.503, Rep = 54.6).
bulk volume fraction with increasing distance from the wall. It shows that the arrangement of the particles in the near-wall region is directly influenced by the presence of the wall. The period of the void fraction distribution is almost the same as dp. Fig. 3b shows the results of an MRI measurement in a cylindrical container, taken by Sederman et al. [14,15]. Unfortunately, the bulk void fractions of Fig. 3a and b are different and a quantitative comparison is difficult. However, a very similar distribution is observed between the results of the calculations of our study and the experimental results, which validates our calculations. Fig. 4 shows a fluid velocity vector distribution in a x–y plane when e = 0.503 and Rep = 54.6. The result is obtained after the flow has fully developed. Fig. 5a and b shows the distributions of the corresponding streamwise fluid velocity components in x–y and
Fig. 6. Time-averaged streamwise velocity and voidage distributions near the wall.
571
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
2
Averaged streamwise velocity / [-]
dp=2mm 1.6
dp=4mm
1.4 1.2
yz
10
5
0.5
1
1.5
2
2.5
3
3.5
0.8 0.6
0
0.5
1
1.5
2
2.5
3
Distance from wall x/dp [-] Fig. 8. Effect of particle diameter on averaged streamwise velocity distribution (e = 0.503, Rep = 25.3).
4
tigations using the same bulk void fraction and relative particle arrangements, the arrangement for dp = 1 mm shown in Fig. 2 is adopted as the standard, and the whole domain, including the particles, is dilated to achieve a specific particle size in the other cases. The results shown here are for e = 0.503 and Rep = 25.3. It is observed that the averaged streamwise fluid velocity in the near-wall region (x/dp < 0.2) increases slightly with increasing particle size, while no significant differences as a result of particle size are observed in the region of x/dp > 0.7.
3 2 1 0 -1 -2
4.4. Influence of wall friction on flow field
-3 -4
1
0.4
0
Distance from wall x/dp [-] (a) Drag force
Lift force FL [μN]
dp=1mm
1.8
y
15
y
Drag force FD [μN]
20
0
0.5
1
1.5
2
2.5
3
3.5
Distance from wall x/dp [-] (b) Lift force Fig. 7. Drag and lift force acting on particles (e = 0.503, Rep = 54.6).
e = 0.503 and Rep = 54.6. It can be observed that there are wide variations in both the drag and the lift forces acting on particles located at almost the same distance from the wall. The differences are as much as three-fold. Meanwhile, no clear dependency on the distance from the wall is observed for both the drag and the lift forces. In the present study, the flows are driven by maintaining the average pressure gradient in the y direction and all investigations are performed in a steady state. This means that the force due to the averaged pressure gradient should be balanced by the resistance of the particles and the wall. Fluid drag force acting on particles does not become large in places where the streamwise averaged fluid velocity is relatively large, as in Fig. 6a. It is evident that the local arrangement of the particles differs with the location, which might be the reason for the variations in the drag and lift forces at the same distance in Fig. 7. Further detailed investigations on this issue are needed. 4.3. Effects of particle size on flow field To investigate the effects of the particle size, dp is varied to 1, 2, and 4 mm. Fig. 8 shows the averaged streamwise fluid velocity distributions for the different particle sizes. To carry out the inves-
In this part of our study, a number of numerical experiments are performed to observe the effects of wall friction on the flow field. The arrangement of the particles is the same as in Section 4.2 and the boundary condition at x/dp = 0 is changed to slip. Fig. 9 shows the results of the averaged streamwise fluid velocity. Here, the same averaged pressure gradients for the cases of Rep = 25.3 and Rep = 28.7 in Fig. 6 are imposed on the flows. In the slip boundary case, the viscous effect on the solid wall eventually vanishes and the fluid velocity increases in the vicinity of the boundary. The regions where large differences are observed with regard to the boundary condition are close to the boundary (0 < x/dp < 0.7 when e = 0.802, and 0 < x/dp < 1.4 when e = 0.802), while the results in the outer regions are almost the same. From these observations, it can be concluded that the effects of wall friction are limited to the near-wall region. The discussions in this section are nothing but for the effects of the boundary layer that is expected to form on the surface of the solid wall. The regions where a dependency on Rep is observed in Fig. 6 are almost the same as those where a dependency on the boundary conditions is observed in Fig. 9. It is expected that the Rep dependency in Fig. 6 is due to the boundary layer that formed on the surface of the solid wall. Furthermore, the regions where particle size dependency is observed in Fig. 8 also coincided with the above; we will further study the relationship between the thickness of the boundary layer on the wall and the particle size. 4.5. Prediction performance of pressure drop in mesoscopic model calculations As noted in Section 1, empirical equations such as those of Ergun [4], Di Felice [5], and Beetstra et al. [6] are used to calculate
572
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
1.2
1.5
0.8
yz
y
0.5
0
0.6
y
1
no-slip slip
1
Averaged streamwise velocity [m/s]
no-slip slip
2
yz
Averaged streamwise velocity [m/s]
2.5
0
0.5
1
1.5
0.4 0.2 0
2
0
0.5
Distance from wall x/dp [-]
1
1.5
2
2.5
3
Distance from wall x/d p [-]
(a) ε = 0.503
(b) ε = 0.802
Fig. 9. Effect of wall friction on averaged streamwise velocity distribution.
4
3.5
Rep=0.128 Rep=25.3 Rep=54.6
3.5 3
3
Rep=0.128 Rep=25.3 Rep=54.6
2.5
2.5 2 2
1.5
1.5
1
1
0.5
0.5 0
0
0.5
1
1.5
2
2.5
3
3.5
0
4
0
0.5
1
1.5
2
2.5
3
3.5
4
Distance from wall x/d p [-]
Distance from wall x/d p [-]
(a) Comparison with Ergun equation ( ε = 0.503)
(b) Comparison with Beetstra equation ( ε = 0.503) 11
30 27.5 25
10
Rep=0.143 Rep=3.98 Rep=28.7
22.5 20
Rep=0.143 Rep=3.98 Rep=28.7
9 8 7
17.5
6
15 12.5
5
10 7.5
4
5 2.5
2
3
1
0
0
0.5
1
1.5
2
2.5
3
3.5
Distance from wall x/d p [-]
(c) Comparison with Ergun equation ( ε = 0.802)
4
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Distance from wall x/d p [-]
(d) Comparison with Beetstra equation ( ε = 0.802)
Fig. 10. Pressure drop predictions at mesoscopic cell level.
573
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
Rep=0.128 Rep=25.3 Rep=54.6
Rep=0.128 Rep=25.3 Rep=54.6
Rep=0.143 Rep=3.98 Rep=28.7
Rep=0.143 Rep=3.98 Rep=28.7
Fig. 11. Effect of averaging cell size in the wall-normal direction.
momentum exchange between the particles and the gas in mesoscopic model calculations. Care is required in using these equations in computational cells where the effect of a wall is not negligible. In this part of our study, the applicability of these empirical equations is investigated. In the IBM calculations, information on the microscopic flows that cannot be obtained by mesoscopic model calculations is directly obtained. In this part of our study, the prediction performance of the aforementioned empirical equations is investigated. In addition to the computational cells used for direct simulation, a mesoscopic scale cell is introduced. The latter cell is hereafter called mesoscopic cell. The positions of all the particles are explicitly known and it is easy to obtain the void fraction in the mesoscopic cells. A locally phase averaged fluid velocity in the mesoscopic cells is calculated by directly averaging the velocity field obtained by the IBM simulations. An estimate of the pressure drop in the mesoscopic cell is obtained by substituting these two quantities into the empirical equations. The prediction performance of the empirical equations is verified by comparing the results obtained by the above procedure with the averaged pressure gradient imposed as a driving force in the microscopic simulations. The dimension of the mesoscopic cell in the wall-normal (x direc-
tion) is set to 0.5dp and the other dimensions are the same as those of the domain of the direct simulation. Fig. 10 shows the results of the Ergun and Beetstra equations, where DpDF is the averaged pressure difference imposed as a driving force in the direct simulations and DpErgun and DpBeetstra are the values respectively predicted by the Ergun and Beetstra equations. The prediction performance of the empirical equations improves as the ratios of DpErgun and DpBeetstra to DpDF approach unity. Regardless of the bulk void fraction, fairly good predictions are obtained in the region of x/dp P 1, i.e., the predicted values and the averaged pressure difference imposed as a driving force mostly agrees. In contrast, it is observed that the Ergun and Beetstra equations significantly underpredict the pressure drop in the near-wall region (0 < x/dp < 1). This is apparent for large void fractions (see Fig. 10c and d) and is more pronounced for low values of Rep. 4.6. Effect of cell size used in mesoscopic model calculations In the previous section, the dimension of the mesoscopic cell in the wall-normal direction was fixed at 0.5dp. In practical mesoscopic model calculations, the cells used are larger than the parti-
574
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
cle structures formed in computational cells. Empirically, it is at least 3dp [16]. In the previous section, the prediction performance of the Ergun and Beetstra equations were observed in the cells near the solid wall. The effect of the solid wall is reduced when the size of the mesoscopic cell is increased. In this section, a mesoscopic cell facing the solid wall is assumed. By changing the size of the cell in the wall-normal direction (Lw), we observe its influence on the accuracy of the pressure drop prediction. As in the previous section, the other dimensions of the cell are the same as those of the computational domain. The results are shown in Fig. 11. As expected, the errors are large when Lw/dp = 1, but gradually decreases as Lw increases. The investigations are conducted up to Lw/dp = 6 and Lw/dp = 8 for e = 0.503 and e = 0.802, respectively, and faster convergence is observed for higher values of Rep. It is observed that Rep does not converge to unity for larger Lw/dp within the range of our investigations. In addition to the effects of the wall, note should also be taken of the prediction accuracy of the empirical equations as discussed in Section 3.2. From the comparison with the case of Lw/dp = 3, which is the cell size typically used in mesoscopic model calculations, we observe that the Ergun equation predicts values that are twice those for a large void fraction e = 0.802 and small Reynolds numbers Rep = 0.143 and Rep = 3.98 (Fig. 11c). However, the difference is less than 1.5 times in the other cases.
with dimensions at least three times the particle diameter, which is the size of cells empirically adopted in mesoscopic calculations. All the particles were fixed in this study and the particle structure should be different from that of moving particles. Careful note should be taken of this fact. Investigations of the case of moving particles will be conducted and the results will be compared with those of this study. As shown in Fig. 6, microscopic flows, which cannot be captured by mesoscopic model calculations, exist in the vicinity of a solid wall. In a large number of applications such as fluidized beds, heat and mass transfers within the equipment are also as important as the flow dynamics. Concerning heat transfer for instance, it is expected that a fluid velocity profile near a solid wall will greatly affect the amount of convective heat transfer through the wall. A study of the sub-grid scale (SGS) model, which represents the flows that occur on the unresolved microscopic scale, is also required.
5. Conclusions
References
In the present study, the body-force type immersed boundary method was used to perform microscopic numerical simulations for the detailed investigation of the effects of a solid wall on flows containing dense particles. As a first report, investigations were performed under fixed particle conditions and the following conclusions are made. Investigations were conducted in relatively dense conditions (e = 0.503 and e = 0.803). In such flows, the wall affects the arrangement of the particles within a distance of three to four times the particle’s diameter from the wall, and the fluid preferentially flows through the large gaps between the particles and between the particles and the wall. The streamwise fluid velocity varies with the bulk velocity of the entire domain, the particle Reynolds number, the void fraction, and the particle size. However, a scaling based on the bulk velocity of the entire domain and the particle diameter was possible in the outer regions, at a distance of at least one particle diameter from the wall. In the region within one diameter from the wall, the effect of the wall friction was dominant and a dependency on the particle Reynolds number and particle diameter was observed. Further investigations into the relationship between the thickness of the viscous sub-layer and the particle diameter are required. Using the results of direct simulations, we also examined the prediction performance of empirical equations popularly used in mesoscopic model calculations. The Ergun and Beetstra equations significantly underestimated the pressure drop in the region within one particle diameter from the wall. However, the predictions were good enough in outer regions. The influence of a wall in mesoscopic model calculations was investigated by changing the size of the mesoscopic cell facing a solid wall. Regarding the pressure drop estimation, the wall does not create serious problems for cells
Acknowledgements This study was supported by the Japan Society for the Promotion of Science (JSPS) under a Grant-in-Aid for Young Scientists (B), Grant No. 22760128. The authors are very grateful for this support.
[1] Y. Tsuji, T. Kawaguchi, T. Tanaka, Discrete particle simulation of 2-dimensional fluidized-bed, Powder Technology 77 (1993) 79–87. [2] B.H. Xu, A.B. Yu, Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics, Chemical Engineering Science 52 (1997) 2785–2809. [3] B.P.B. Hoomans, J.A.M. Kuipers, W.J. Briels, W.P.M. vanSwaaij, Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach, Chemical Engineering Science 51 (1996) 99–118. [4] S. Ergun, Fluid flow through packed columns, Chemical Engineering Progress Symposium Series 48 (1952) 89–94. [5] R. Di Felice, The voidage function for fluid particle interaction systems, International Journal of Multiphase Flow 20 (1994) 153–159. [6] R. Beetstra, M.A. van der Hoef, J.A.M. Kuipers, Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE Journal 53 (2007) 489–501. [7] D. Mehta, M.C. Hawley, Wall effect in packed columns, Industrial & Engineering Chemistry Process Design and Development 8 (1969) 280–282. [8] B. Eisfeld, K. Schnitzlein, The influence of confining walls on the pressure drop in packed beds, Chemical Engineering Science 56 (2001) 4321–4329. [9] X. Ren, S. Stapf, B. Blümich, Magnetic resonance visualisation of flow and pore structure in packed beds with low aspect ratio, Chemical Engineering Technology 28 (2005) 219–225. [10] T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, Turbulence structure of particle-laden flow in a vertical plane channel due to vortex shedding, JSME International, B 44 (2001) 526–535. [11] T. Tsuji, R. Narutomi, T. Yokomine, S. Ebara, A. Shimizu, Unsteady threedimensional simulation of interactions between flow and two particles, International Journal of Multiphase Flow 29 (2003) 1431–1450. [12] H. Yada, Direct Numerical Simulation of Flow with Dense Particle, Master thesis, Osaka University, 2010. [13] L. Zeng, S. Balachandar, P. Fischer, Wall-induced forces on a rigid sphere at finite Reynolds number, Journal of Fluid Mechanics 536 (2005) 1–25. [14] A.J. Sederman, P. Alexander, L.F. Gladden, Structure of packed beds probed by magnetic resonance imaging, Powder Technology 117 (2001) 255–269. [15] L. Shuiqing, D. Yulong, W. Dongsheng, H. Yurong, Modelling of the behaviour of gas–solid two phase mixtures flowing through packed beds, Chemical Engineering Progress Symposium Series 61 (2006) 1922–1931. [16] T. Tsuji, K. Yabumoto, T. Tanaka, Spontaneous structures in three-dimensional bubbling gas-fluidized bed by parallel DEM–CFD coupling simulation, Powder Technology 184 (2008) 132–140.
Contents lists available at SciVerse ScienceDirect
Advanced Powder Technology journal homepage: www.elsevier.com/locate/apt
Translated paper; SPTJ Best Paper Award 2011
Effect of a wall on flow with dense particles q Takuya Tsuji ⇑, Eiji Narita, Toshitsugu Tanaka Department of Mechanical Engineering, Graduate School of Engineering, 2-1 Yamada-Oka, Suita, Osaka 565-0871, Japan
a r t i c l e
i n f o
Article history: Received 6 November 2012 Accepted 12 November 2012 Available online 30 November 2012 Keywords: Wall effect Gas–solid flow Fluidized bed Direct numerical simulation Immersed boundary method
a b s t r a c t The behavior of dense gas–solid flows in engineering applications such as fluidized beds and pneumatic conveyers is highly complex and a reliable numerical model is required. Such flows are usually within solid walls that considerably affect the flow fields, and it is important to correctly include this effect in numerical models to improve their prediction capability. The observation of microscopic flows near walls can enhance our understanding of the flow behavior and assist in improving models. In this study, direct simulations are performed to investigate the effect of a wall on flow fields at a microscopic level. The effects of the bulk void fraction, particle Reynolds number, and particle diameter are investigated. The prediction performances of existing correlation equations usually used in mesoscopic model calculations are also investigated. It is found that the Ergun and Beetstra equations produce large discrepancies in the region within a distance equal to the particle diameter from the wall. Ó 2012 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights reserved.
1. Introduction Dense particles coexist with a gas in industrial equipment such as fluidized beds and pneumatic conveying systems. The flow behavior of the gas–particle mixture is significantly influenced by the interactions between the particles, between the particles and the wall and between the particles and the gas. The spontaneous formation of mesoscopic heterogeneous structures inside the equipment, such as bubbles that are far larger than the particles, makes the flow more complex. A better understanding of the flow behavior and design optimization is therefore desired. This has, however, not been successfully achieved owing to observation difficulties in experiments. Against this background, the development of numerical models has been extensively pursued until the present time. The majority of these models are mesoscopic in the sense that the size of the computational cells of the fluid motion is small enough to resolve characteristic mesoscopic structures but too large to directly capture the microscopic phenomena on the particle level [1–3]. The models have been applied to several flow problems, and the reproduction of size and frequency with which bubbles appear in fluidized beds have been successfully demonstrated. In practical engineering situations, solid walls that significantly influence the behavior of both phases are very common. It is important to correctly express the influence of the solid walls whenever a numerical simulation is performed. In mesoscopic model calcula-
q Japanese version published in JSPTJ, Vol. 48 No. 12 (2011); English version for APT received on November 6, 2012. ⇑ Corresponding author. Tel./fax: +81 6 6879 7317. E-mail address: [email protected] (T. Tsuji).
tions, the momentum exchange between solid particles and the gas in a computational cell is generally determined by using empirical equations that are a function of the void fraction and the particle Reynolds number. In the majority of existing models, empirical equations such as those of Ergun [4], Di Felice [5], and Beetstra et al. [6] are used without any modifications, including in computational cells near a solid wall. The existence of a wall directly exerts viscous friction on nearby fluids and affects the arrangement of nearby particles. As a result, the anisotropy of the flow in the vicinity of solid walls is increased. In the aforementioned empirical equations, the effect of the wall is not considered and it is questionable to use these equations for computational cells near solid walls. Moreover, the problem is one that has been recognized by researchers since early times and was mentioned by Beetstra et al. [6]. When microscopic flows are directly considered, fluid motions should obey the Navier–Stokes and continuity equations, and it is appropriate to use the no-slip boundary condition on the surface of solid walls. However, in mesoscopic model calculations based on locally phase averaged equations, microscopic flows occurring on a subgrid scale are averaged and it is questionable to use the no-slip boundary condition on the surface of solid walls on the mesoscopic cell scale. Currently, different boundary conditions are used depending on a study, and a general consensus is yet to be achieved. For further improvement of mesoscopic models, it is important to adopt a proper drag correlation in computational cells near a solid wall and a proper boundary condition on the wall. This cannot be accomplished without a thorough understanding of the behavior of microscopic flows near such walls, which never be captured by the mesoscopic model calculations. Although the effects of a solid wall on a flow that contains dense particles have been investigated [7–9], it is still difficult to directly observe microscopic flows
0921-8831/$ - see front matter Ó 2012 The Society of Powder Technology Japan. Published by Elsevier B.V. and The Society of Powder Technology Japan. All rights reserved. http://dx.doi.org/10.1016/j.apt.2012.11.006
566
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
Nomenclature a dp FD FL fp Gp g I Ip L Lx, y, z Mp mp Np n nx, y, z p Rep r Sp t U Up u
particle radius (m) particle diameter (m) fluid drag force (N) fluid lift force (N) interaction force (m/s2) external force acting on a particle (N) acceleration due to gravity (m/s2) unit tensor (–) moment of inertia (kg m2) streamwise domain size (m) domain size in x, y, and z directions (m) external moment acting on a particle (N m) particle mass (kg) number of particles in the computational domain (–) vector normal to the surface of a particle (–) number of grid points in x, y, and z directions (–) fluid pressure (Pa) particle Reynolds number (–) relative vector from the center of a particle (–) surface area of a particle (m2) time (s) inflow velocity (m/s) translating velocity of a particle (m/s) coupling velocity weighted by particle volume fraction (m/s)
on a particle scale by experimental methods. In the present study, the body-force-type immersed boundary method (IBM) is used to conduct microscopic numerical simulations of flows that contain dense solid particles near a solid wall, and the effects of the wall on the flow behavior are investigated in detail. In addition, the data obtained from the microscopic simulations are used to discuss the validity of the empirical drag equations popularly used in mesoscopic model calculations. In this paper, which is the first report, we restrict our discussions to the case of a fixed particle. 2. Calculation methods 2.1. Governing equations It is assumed that the fluid is incompressible and Newtonian and the particles are rigid spheres. The governing equations of the fluid flow are the continuity equation and Navier–Stokes equation, respectively, expressed below:
r uf ¼ 0
qf
Duf ¼ r s þ qf g Dt
s ¼ pI þ lf ½ruf þ ðruf Þ
dðmp U p Þ ¼ dt dðI p Xp Þ ¼ dt
s n dS þ Gp
Sp
r ðs nÞ dS þ M p
e qf lf mf s Xp
w hi h iyz
where Up and Xp respectively represent the translational and rotational velocities of the particle, mp is the mass of the particle, Ip is the moment of inertia given by Ip = (2/5)a2mpI for a particle with radius a, Gp and Mp are respectively the external force and moment, and r is the relative position vector from the center of a particle. 2.2. Body-force-type IBM [10] Kajishima et al. developed the body-force-type IBM for particle– fluid systems. Here, the assumption of a uniform Cartesian grid is used to simplify the explanation. The presence of particles is represented by the solid volume fraction of particles in each computational cell a. Kajishima et al. [10] introduced the following coupling velocity weighted by the solid volume fraction of particles:
u ¼ aup þ ð1 aÞuf
ð6Þ
ð2Þ
ru¼0
ð3Þ
ð4Þ
Sp
Z
Dp DpBeetstra DpDF DpErgun Dt Dx
ð1Þ
The particles are tracked individually in a Lagrangian manner by using the equations of translational and rotational motions:
Z
a
superficial velocity (m/s) fluid velocity (m/s) particle velocity (m/s) streamwise fluid velocity (m/s) velocity predicted as a fluid (m/s) volume of computational domain (m3) volume of a particle (m3) volume fraction of particle in each computational cell (–) pressure drop (Pa) pressure drop predicted by Beetstra equation (Pa) pressure difference imposed as a driving force (Pa) pressure drop predicted by Ergun equation (Pa) temporal increment (s) computational cell size (m) void fraction in whole region (–) density of fluid (kg/m3) viscosity of fluid (Pa s) kinetic viscosity of fluid (m2/s) stress tensor (Pa) rotational velocity of a particle (rad/s) nondimensional pressure drop (–) spatial average in whole region spatial average in y–z plane temporal average
where up is the particle velocity given by up = Up + r Xp. The particles are solid and no-slip and no-permeability conditions are imposed on their interfaces. Hence, the continuity restriction also applies to u:
where uf and qf are respectively the fluid velocity and density, g is the acceleration due to gravity, p is the pressure, lf is the fluid viscosity, and s is the stress tensor given by T
u0 uf up uy ~ u V Vp
ð5Þ
ð7Þ
The following equation for u is introduced:
@u 1 ¼ rp þ H þ f p @t qf
ð8Þ
where
H ¼ uru þ mf r ½ru þ ðruÞT þ g
ð9Þ
and mf is the kinematic viscosity of the fluid. In Eq. (8), fp is the interaction force term that enforces all the predicted velocity fields (including inside the particles) to the particle–fluid coupling velocity field. This is explained in terms of the time-marching procedure. If Eq. (8) marches in time with the explicit Euler method, it would be
567
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
unþ1 ¼ un þ Dt
1
qf
!n
rp þ H þ f p
ð10Þ
Rep (=Udp/vf)
where the superscript represents time and Dt is the time increment. We predict the velocity as the fluid phase once regardless of the cell is occupied by fluid (a = 0), particle (a = 1) and both (0 < a < 1).
~ ¼ un þ Dt u
1
qf
Table 1 Comparison of predicted drag coefficient with Zeng et al. [13].
10 200
CD
Difference (%)
Present
Zeng et al. [13]
4.988 0.8719
4.721 0.8156
5.655 6.902
!n
rp þ H
ð11Þ
The predicted velocity should be modified by fp to meet the definition of u in Eq. (6). For the cell inside the particle (a = 1), ~ Þ=Dt gives un+1 = up (inside the particle). However, for f p ¼ ðup u the cell occupied by fluid (a = 0), Eq. (8) is identical to Eq. (2) because fp vanishes. Hence, the added term fp is modeled with the linear interpolation of a as
~ Þ=Dt f p ¼ aðup u
ð12Þ
Table 2 Calculation conditions for packed bed without wall. Computational domain: Lx Ly Lz (mm) Number of grid points: nx ny nz (–) Particle diameter: dp (mm) Number of particles: Np (–) Voidage: e = h1 ai (–) Fluid density: qf (kg/m3) Kinetic viscosity: mf (Pa s)
6dp 10dp 6dp 96 160 96 1.0 360 0.4764 1.205 1.448 105
This force is used to modify the solution of Eq. (11) to u (Eq. (6)). Additionally, the term is used for determining the momentum exchange between two phases through the interface. The fluid force acting on a particle surface (Eqs. (4) and (5)) is obtained by integrating the mutual interaction force fp with respect to the volume:
dðmp U p Þ ¼ qf dt dðI p Xp Þ ¼ qf dt
Z
f p dV þ Gp
ð13Þ
Vp
Z
r f p dV þ M p
ð14Þ
Vp
The present method is applicable to moving particles. In all calculations in this study, the particles are fixed. The numerical calculations are performed using the finite difference method based on the fractional step method. A staggered uniform Cartesian grid is used. Fourth-order central difference and second-order Adams– Bashforth schemes are respectively used for the discretization of the spatial and temporal derivatives. The method proposed by Tsuji et al. [11] is used to obtain the solid volume fraction in each computational cell. 3. Verification of calculation accuracy In IBM calculations, solid particles are represented on a Cartesian grid and the accuracy of the calculations varies, depending on the particle resolution, which is the number of computational cells per particle diameter. Although improved accuracy can be expected with higher particle resolutions, excessive resolutions significantly increase computational load and should be avoided. In advance to the main calculations, the relationship between the particle resolution and the calculation accuracy should be investigated. A number of investigations have been carried out, including those on the fluid drag force acting on single and paired particles in a uniform flow, the Saffman lift force, the Magnus lift force acting on a particle, and the viscous torque acting on a rotating particle [11,12]. Our points of focus in the present study are the flow containing dense particles near a solid wall, verifications study are performed for the fluid drag force acting on a single particle moving steadily near a solid wall and the pressure drop of a fluid in a packed bed without solid walls. 3.1. Fluid drag force acting on a single particle near a solid wall The accuracy of present calculations is verified in this section. A single particle traveling through an acquiescent fluid near a solid wall with a constant velocity U is considered. The distance between the particle surface and the solid wall is fixed at 1/2dp. Based on the
Fig. 1. Comparison with Ergun and Beetstra equations.
results of previous resolution-dependency studies, the number of computational cells per particle diameter is set to dp/Dx = 16. The results are compared with those that Zeng et al. obtained using the spectral element method [13]. The particle Reynolds number, defined as Rep = Udp/mf, is set to 10 and 200, where dp is the particle diameter and mf is the dynamic viscosity of the fluid. The results are shown in Table 1. For both values of Rep, the agreement between the results of our calculation and those of Zeng et al.’s are fairly good. The values of our results are slightly higher but the differences are within 7%. The same particle resolution dp/Dx = 16 is hereafter used for our study. 3.2. Pressure drop in packed beds Calculations of the pressure drop in packed beds without solid walls are performed and the results are compared with those obtained with the empirical equations by Ergun [4] and Beetstra et al. [6]. Ergun [4] expressed the pressure drop in packed beds as
ð1 eÞlf u0 Dp 1 e 2 150 þ 1:75 q u ¼ f 0 dp e3 dp L
ð15Þ
where Dp is the pressure difference, L is the bed height, e is the void fraction, dp is the particle diameter, lf is the viscosity of the fluid, u0 is the superficial velocity, and qf is the density of the fluid. Meanwhile, Beetstra et al. [6] determined the fluid drag force FD acting on a single particle in packed beds from the results of the lattice Boltzmann simulation:
568
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574 Table 3 Calculation conditions. Computational domain: Lx Ly Lz (–) Number of grid points: nx ny nz (–) Particle diameter: dp (mm) Number of particles: Np (–) Reynolds number: Rep ¼ e< uy >dp =v f (–) dp/Dx (–) Voidage: e = h1 ai (–) Fluid density: qf (kg/m3) Kinetic viscosity: mf (Pa s)
6dp 10dp 6dp 96 160 96
8.25dp 13.5dp 8.25dp 132 216 132 1.0 360 0.1–54.6 16 0.802 1.205 1.448 105
0.503
y
10d p [mm]
13.5dp [mm]
x
8.25dp [mm]
6dp [mm] (a) ε = 0.503
(b) ε = 0.802 Fig. 2. Calculation domain.
F D ¼ 3plf dp u0 Fðe; Rep Þ Fðe; Rep Þ ¼
10ð1 eÞ
e
3
ð16Þ
pffiffiffiffiffiffiffiffiffiffiffi þ eð1 þ 1:5 1 eÞ
0:413Rep þ 24e3
! e1 þ 3eð1 eÞ þ 8:4Re0:343 p 1 þ 103ð1eÞ Rep2e2:5
w¼ ð17Þ
Let us assume that Np particles are present in a control volume V with height L. Considering that the pressure drop in the control volume is due only to the particles in the box, the total fluid drag force acting on the Np particles should be balanced by the pressure drop. Hence,
F D Np ¼
Dp V L
ð18Þ
The void fraction e in the box can be related to the volume of each particle, Vp, by
1e¼
Np V p V
ð19Þ
Hence, the pressure drop of the Beetstra equation is obtained with
DP Vp ¼ FD L 1e
Hereafter in this paper, the pressure drop is discussed in terms of the following non-dimensional form:
ð20Þ
DP dp L qf u20
ð21Þ
The calculation conditions are listed in Table 2. In the table, hi shows the spatial averaging over the whole calculation domain, including the inside of solid particles, and the bulk void fraction e can be calculated from the solid volume fraction at each computational cell using
e ¼ h1 ai
ð22Þ
In this section, the periodic boundary condition is used in all directions. To obtain a random arrangement of the particles, the position of each particle is fixed after tracking the motion of all the particles for a sufficient time imposing an initial random velocity and an artificial potential force between particles. During the packing, fluid calculations are not performed. The constant averaged pressure gradients in the y direction are the driving forces and the calculations are performed by varying the particle Reynolds number. The results are not included here; however, the sufficiency of the domain size was confirmed by preliminary calculations. The working fluid is air at room temperature. Fig. 1 shows the results, where Rep, the
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
569
4. Effects of a solid wall 4.1. Calculation conditions In this section, the influence of a wall on a flow field is investigated. The calculation conditions are listed in Table 3. As in the part discussed in the previous section, random packing of the particles near the solid wall at x = 0 is obtained by tracking all the particles without fluid calculations. An artificial potential force is also applied between the particles and the wall, in addition to between the particles. During particle packing, the periodic boundary condition is used in the y and z directions and solid walls are placed at x = 0 and x = Lx + dp, respectively. Flows near a single solid wall at x = 0 are of interest in the present study; however, the effect of the wall on the particle arrangement at x = Lx + dp would remain with above procedure. To minimize the effect, a sub-domain within Lx < x 6 Lx + dp is eliminated from the computations and the slip boundary condition is imposed instead at x = Lx. The flows are maintained by imposing a constant averaged pressure gradient in the y direction. After a steady state is attained, the flow field is observed. In this study, different bulk void fractions are achieved by varying the calculation domain size while maintaining the number of particles included in the domain. Snapshots of the particle arrangements for e = 0.503 and e = 0.802 are shown in Fig. 2a and b, respectively. Based on the results of the verification process discussed in the previous section, the investigations are restricted to the case of Rep < 60. Roughly speaking, this condition corresponds to that of free bubbling in beds with 500 lm particles. 4.2. Observation of flow field near a solid wall Fig. 3a shows the void fraction distribution h1 aiyz at e = 0.503. The result is averaged in the y–z planes. It is observed that the void fraction has an oscillatory profile and progressively approaches the
Fig. 3. Voidage distribution near the wall.
particle Reynolds number based on a superficial velocity, is defined by
Rep ¼
ehuy idp mf
ð23Þ
where — indicates temporal averaging. A fairly good agreement between the results of our calculations and those obtained with empirical equations is observed; the errors are less than 20% in the region of Rep < 60. The difference increases as Rep increases, which is a result of the resolution limitations of our calculations; we cannot capture the velocity gradient of the fluid on the surface of the particles. Improved accuracy is expected upon increasing the resolution of the simulations, although this would also lead to an increase in computational cost. In our study, we restrict our investigations to the region of Rep < 60, where the accuracy of the calculations is guaranteed. While the verification process of this section is conducted by comparing our calculation results with those obtained with empirical equations, it must be noted that there is considerable difference between the results by the Ergun and Beetstra equations, especially in the region of Rep > 100. Further investigations into these empirical equations may also be required.
Fig. 4. Instantaneous velocity vector distribution in an x–y plane (e = 0.503, Rep = 54.6).
570
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
x–z planes, respectively. It can be observed that the streamwise velocity component becomes relatively large in the regions with higher void fractions. This is especially apparent in the near-wall region where x/dp = 0 in Fig. 5b. Fig. 6 shows the distributions of the averaged streamwise velocity huy iyz for different particle Reynolds numbers. The void fraction distribution is also shown. Different particle Reynolds numbers are obtained by changing the value of the imposed averaged pressure gradient. The averaged streamwise velocity is normalized by the bulk streamwise velocity in the entire domain < uy >. Except in the vicinity of the wall where the influence of the no-slip boundary condition on the wall surface exists, the streamwise velocity is large whenever the void fraction is large. The results correspond to Fig. 4 and show that the flows preferentially flow through large gaps between the particles and between the particles and the wall. A dependency on Rep is only observed in the vicinity of the wall, where the streamwise velocity becomes relatively large under higher Rep conditions. This is not observed in the far-wall regions, i.e., when x/dp > 0.7 in Fig. 6a and x/dp > 1.6 in Fig. 6b. In the body-force type IBM, the fluid force acting on each particle can be easily obtained with Eq. (13). Fig. 7 shows the relationship between the distance of the particle’s center from the wall and the time-averaged drag and lift forces acting on each particle when
Fig. 5. Streamwise velocity component distributions (e = 0.503, Rep = 54.6).
bulk volume fraction with increasing distance from the wall. It shows that the arrangement of the particles in the near-wall region is directly influenced by the presence of the wall. The period of the void fraction distribution is almost the same as dp. Fig. 3b shows the results of an MRI measurement in a cylindrical container, taken by Sederman et al. [14,15]. Unfortunately, the bulk void fractions of Fig. 3a and b are different and a quantitative comparison is difficult. However, a very similar distribution is observed between the results of the calculations of our study and the experimental results, which validates our calculations. Fig. 4 shows a fluid velocity vector distribution in a x–y plane when e = 0.503 and Rep = 54.6. The result is obtained after the flow has fully developed. Fig. 5a and b shows the distributions of the corresponding streamwise fluid velocity components in x–y and
Fig. 6. Time-averaged streamwise velocity and voidage distributions near the wall.
571
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
2
Averaged streamwise velocity / [-]
dp=2mm 1.6
dp=4mm
1.4 1.2
yz
10
5
0.5
1
1.5
2
2.5
3
3.5
0.8 0.6
0
0.5
1
1.5
2
2.5
3
Distance from wall x/dp [-] Fig. 8. Effect of particle diameter on averaged streamwise velocity distribution (e = 0.503, Rep = 25.3).
4
tigations using the same bulk void fraction and relative particle arrangements, the arrangement for dp = 1 mm shown in Fig. 2 is adopted as the standard, and the whole domain, including the particles, is dilated to achieve a specific particle size in the other cases. The results shown here are for e = 0.503 and Rep = 25.3. It is observed that the averaged streamwise fluid velocity in the near-wall region (x/dp < 0.2) increases slightly with increasing particle size, while no significant differences as a result of particle size are observed in the region of x/dp > 0.7.
3 2 1 0 -1 -2
4.4. Influence of wall friction on flow field
-3 -4
1
0.4
0
Distance from wall x/dp [-] (a) Drag force
Lift force FL [μN]
dp=1mm
1.8
y
15
y
Drag force FD [μN]
20
0
0.5
1
1.5
2
2.5
3
3.5
Distance from wall x/dp [-] (b) Lift force Fig. 7. Drag and lift force acting on particles (e = 0.503, Rep = 54.6).
e = 0.503 and Rep = 54.6. It can be observed that there are wide variations in both the drag and the lift forces acting on particles located at almost the same distance from the wall. The differences are as much as three-fold. Meanwhile, no clear dependency on the distance from the wall is observed for both the drag and the lift forces. In the present study, the flows are driven by maintaining the average pressure gradient in the y direction and all investigations are performed in a steady state. This means that the force due to the averaged pressure gradient should be balanced by the resistance of the particles and the wall. Fluid drag force acting on particles does not become large in places where the streamwise averaged fluid velocity is relatively large, as in Fig. 6a. It is evident that the local arrangement of the particles differs with the location, which might be the reason for the variations in the drag and lift forces at the same distance in Fig. 7. Further detailed investigations on this issue are needed. 4.3. Effects of particle size on flow field To investigate the effects of the particle size, dp is varied to 1, 2, and 4 mm. Fig. 8 shows the averaged streamwise fluid velocity distributions for the different particle sizes. To carry out the inves-
In this part of our study, a number of numerical experiments are performed to observe the effects of wall friction on the flow field. The arrangement of the particles is the same as in Section 4.2 and the boundary condition at x/dp = 0 is changed to slip. Fig. 9 shows the results of the averaged streamwise fluid velocity. Here, the same averaged pressure gradients for the cases of Rep = 25.3 and Rep = 28.7 in Fig. 6 are imposed on the flows. In the slip boundary case, the viscous effect on the solid wall eventually vanishes and the fluid velocity increases in the vicinity of the boundary. The regions where large differences are observed with regard to the boundary condition are close to the boundary (0 < x/dp < 0.7 when e = 0.802, and 0 < x/dp < 1.4 when e = 0.802), while the results in the outer regions are almost the same. From these observations, it can be concluded that the effects of wall friction are limited to the near-wall region. The discussions in this section are nothing but for the effects of the boundary layer that is expected to form on the surface of the solid wall. The regions where a dependency on Rep is observed in Fig. 6 are almost the same as those where a dependency on the boundary conditions is observed in Fig. 9. It is expected that the Rep dependency in Fig. 6 is due to the boundary layer that formed on the surface of the solid wall. Furthermore, the regions where particle size dependency is observed in Fig. 8 also coincided with the above; we will further study the relationship between the thickness of the boundary layer on the wall and the particle size. 4.5. Prediction performance of pressure drop in mesoscopic model calculations As noted in Section 1, empirical equations such as those of Ergun [4], Di Felice [5], and Beetstra et al. [6] are used to calculate
572
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
1.2
1.5
0.8
yz
y
0.5
0
0.6
y
1
no-slip slip
1
Averaged streamwise velocity [m/s]
no-slip slip
2
yz
Averaged streamwise velocity [m/s]
2.5
0
0.5
1
1.5
0.4 0.2 0
2
0
0.5
Distance from wall x/dp [-]
1
1.5
2
2.5
3
Distance from wall x/d p [-]
(a) ε = 0.503
(b) ε = 0.802
Fig. 9. Effect of wall friction on averaged streamwise velocity distribution.
4
3.5
Rep=0.128 Rep=25.3 Rep=54.6
3.5 3
3
Rep=0.128 Rep=25.3 Rep=54.6
2.5
2.5 2 2
1.5
1.5
1
1
0.5
0.5 0
0
0.5
1
1.5
2
2.5
3
3.5
0
4
0
0.5
1
1.5
2
2.5
3
3.5
4
Distance from wall x/d p [-]
Distance from wall x/d p [-]
(a) Comparison with Ergun equation ( ε = 0.503)
(b) Comparison with Beetstra equation ( ε = 0.503) 11
30 27.5 25
10
Rep=0.143 Rep=3.98 Rep=28.7
22.5 20
Rep=0.143 Rep=3.98 Rep=28.7
9 8 7
17.5
6
15 12.5
5
10 7.5
4
5 2.5
2
3
1
0
0
0.5
1
1.5
2
2.5
3
3.5
Distance from wall x/d p [-]
(c) Comparison with Ergun equation ( ε = 0.802)
4
0
0
0.5
1
1.5
2
2.5
3
3.5
4
Distance from wall x/d p [-]
(d) Comparison with Beetstra equation ( ε = 0.802)
Fig. 10. Pressure drop predictions at mesoscopic cell level.
573
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
Rep=0.128 Rep=25.3 Rep=54.6
Rep=0.128 Rep=25.3 Rep=54.6
Rep=0.143 Rep=3.98 Rep=28.7
Rep=0.143 Rep=3.98 Rep=28.7
Fig. 11. Effect of averaging cell size in the wall-normal direction.
momentum exchange between the particles and the gas in mesoscopic model calculations. Care is required in using these equations in computational cells where the effect of a wall is not negligible. In this part of our study, the applicability of these empirical equations is investigated. In the IBM calculations, information on the microscopic flows that cannot be obtained by mesoscopic model calculations is directly obtained. In this part of our study, the prediction performance of the aforementioned empirical equations is investigated. In addition to the computational cells used for direct simulation, a mesoscopic scale cell is introduced. The latter cell is hereafter called mesoscopic cell. The positions of all the particles are explicitly known and it is easy to obtain the void fraction in the mesoscopic cells. A locally phase averaged fluid velocity in the mesoscopic cells is calculated by directly averaging the velocity field obtained by the IBM simulations. An estimate of the pressure drop in the mesoscopic cell is obtained by substituting these two quantities into the empirical equations. The prediction performance of the empirical equations is verified by comparing the results obtained by the above procedure with the averaged pressure gradient imposed as a driving force in the microscopic simulations. The dimension of the mesoscopic cell in the wall-normal (x direc-
tion) is set to 0.5dp and the other dimensions are the same as those of the domain of the direct simulation. Fig. 10 shows the results of the Ergun and Beetstra equations, where DpDF is the averaged pressure difference imposed as a driving force in the direct simulations and DpErgun and DpBeetstra are the values respectively predicted by the Ergun and Beetstra equations. The prediction performance of the empirical equations improves as the ratios of DpErgun and DpBeetstra to DpDF approach unity. Regardless of the bulk void fraction, fairly good predictions are obtained in the region of x/dp P 1, i.e., the predicted values and the averaged pressure difference imposed as a driving force mostly agrees. In contrast, it is observed that the Ergun and Beetstra equations significantly underpredict the pressure drop in the near-wall region (0 < x/dp < 1). This is apparent for large void fractions (see Fig. 10c and d) and is more pronounced for low values of Rep. 4.6. Effect of cell size used in mesoscopic model calculations In the previous section, the dimension of the mesoscopic cell in the wall-normal direction was fixed at 0.5dp. In practical mesoscopic model calculations, the cells used are larger than the parti-
574
T. Tsuji et al. / Advanced Powder Technology 24 (2013) 565–574
cle structures formed in computational cells. Empirically, it is at least 3dp [16]. In the previous section, the prediction performance of the Ergun and Beetstra equations were observed in the cells near the solid wall. The effect of the solid wall is reduced when the size of the mesoscopic cell is increased. In this section, a mesoscopic cell facing the solid wall is assumed. By changing the size of the cell in the wall-normal direction (Lw), we observe its influence on the accuracy of the pressure drop prediction. As in the previous section, the other dimensions of the cell are the same as those of the computational domain. The results are shown in Fig. 11. As expected, the errors are large when Lw/dp = 1, but gradually decreases as Lw increases. The investigations are conducted up to Lw/dp = 6 and Lw/dp = 8 for e = 0.503 and e = 0.802, respectively, and faster convergence is observed for higher values of Rep. It is observed that Rep does not converge to unity for larger Lw/dp within the range of our investigations. In addition to the effects of the wall, note should also be taken of the prediction accuracy of the empirical equations as discussed in Section 3.2. From the comparison with the case of Lw/dp = 3, which is the cell size typically used in mesoscopic model calculations, we observe that the Ergun equation predicts values that are twice those for a large void fraction e = 0.802 and small Reynolds numbers Rep = 0.143 and Rep = 3.98 (Fig. 11c). However, the difference is less than 1.5 times in the other cases.
with dimensions at least three times the particle diameter, which is the size of cells empirically adopted in mesoscopic calculations. All the particles were fixed in this study and the particle structure should be different from that of moving particles. Careful note should be taken of this fact. Investigations of the case of moving particles will be conducted and the results will be compared with those of this study. As shown in Fig. 6, microscopic flows, which cannot be captured by mesoscopic model calculations, exist in the vicinity of a solid wall. In a large number of applications such as fluidized beds, heat and mass transfers within the equipment are also as important as the flow dynamics. Concerning heat transfer for instance, it is expected that a fluid velocity profile near a solid wall will greatly affect the amount of convective heat transfer through the wall. A study of the sub-grid scale (SGS) model, which represents the flows that occur on the unresolved microscopic scale, is also required.
5. Conclusions
References
In the present study, the body-force type immersed boundary method was used to perform microscopic numerical simulations for the detailed investigation of the effects of a solid wall on flows containing dense particles. As a first report, investigations were performed under fixed particle conditions and the following conclusions are made. Investigations were conducted in relatively dense conditions (e = 0.503 and e = 0.803). In such flows, the wall affects the arrangement of the particles within a distance of three to four times the particle’s diameter from the wall, and the fluid preferentially flows through the large gaps between the particles and between the particles and the wall. The streamwise fluid velocity varies with the bulk velocity of the entire domain, the particle Reynolds number, the void fraction, and the particle size. However, a scaling based on the bulk velocity of the entire domain and the particle diameter was possible in the outer regions, at a distance of at least one particle diameter from the wall. In the region within one diameter from the wall, the effect of the wall friction was dominant and a dependency on the particle Reynolds number and particle diameter was observed. Further investigations into the relationship between the thickness of the viscous sub-layer and the particle diameter are required. Using the results of direct simulations, we also examined the prediction performance of empirical equations popularly used in mesoscopic model calculations. The Ergun and Beetstra equations significantly underestimated the pressure drop in the region within one particle diameter from the wall. However, the predictions were good enough in outer regions. The influence of a wall in mesoscopic model calculations was investigated by changing the size of the mesoscopic cell facing a solid wall. Regarding the pressure drop estimation, the wall does not create serious problems for cells
Acknowledgements This study was supported by the Japan Society for the Promotion of Science (JSPS) under a Grant-in-Aid for Young Scientists (B), Grant No. 22760128. The authors are very grateful for this support.
[1] Y. Tsuji, T. Kawaguchi, T. Tanaka, Discrete particle simulation of 2-dimensional fluidized-bed, Powder Technology 77 (1993) 79–87. [2] B.H. Xu, A.B. Yu, Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics, Chemical Engineering Science 52 (1997) 2785–2809. [3] B.P.B. Hoomans, J.A.M. Kuipers, W.J. Briels, W.P.M. vanSwaaij, Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach, Chemical Engineering Science 51 (1996) 99–118. [4] S. Ergun, Fluid flow through packed columns, Chemical Engineering Progress Symposium Series 48 (1952) 89–94. [5] R. Di Felice, The voidage function for fluid particle interaction systems, International Journal of Multiphase Flow 20 (1994) 153–159. [6] R. Beetstra, M.A. van der Hoef, J.A.M. Kuipers, Drag force of intermediate Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE Journal 53 (2007) 489–501. [7] D. Mehta, M.C. Hawley, Wall effect in packed columns, Industrial & Engineering Chemistry Process Design and Development 8 (1969) 280–282. [8] B. Eisfeld, K. Schnitzlein, The influence of confining walls on the pressure drop in packed beds, Chemical Engineering Science 56 (2001) 4321–4329. [9] X. Ren, S. Stapf, B. Blümich, Magnetic resonance visualisation of flow and pore structure in packed beds with low aspect ratio, Chemical Engineering Technology 28 (2005) 219–225. [10] T. Kajishima, S. Takiguchi, H. Hamasaki, Y. Miyake, Turbulence structure of particle-laden flow in a vertical plane channel due to vortex shedding, JSME International, B 44 (2001) 526–535. [11] T. Tsuji, R. Narutomi, T. Yokomine, S. Ebara, A. Shimizu, Unsteady threedimensional simulation of interactions between flow and two particles, International Journal of Multiphase Flow 29 (2003) 1431–1450. [12] H. Yada, Direct Numerical Simulation of Flow with Dense Particle, Master thesis, Osaka University, 2010. [13] L. Zeng, S. Balachandar, P. Fischer, Wall-induced forces on a rigid sphere at finite Reynolds number, Journal of Fluid Mechanics 536 (2005) 1–25. [14] A.J. Sederman, P. Alexander, L.F. Gladden, Structure of packed beds probed by magnetic resonance imaging, Powder Technology 117 (2001) 255–269. [15] L. Shuiqing, D. Yulong, W. Dongsheng, H. Yurong, Modelling of the behaviour of gas–solid two phase mixtures flowing through packed beds, Chemical Engineering Progress Symposium Series 61 (2006) 1922–1931. [16] T. Tsuji, K. Yabumoto, T. Tanaka, Spontaneous structures in three-dimensional bubbling gas-fluidized bed by parallel DEM–CFD coupling simulation, Powder Technology 184 (2008) 132–140.