470
Ind. Eng. Chem. Res. 2008, 47, 470-480
Computational Investigation of Horizontal Slug Flow in Pneumatic Conveying S. B. Kuang,† K. W. Chu,† A. B. Yu,*,† Z. S. Zou,‡ and Y. Q. Feng§ Centre for Simulation and Modelling of Particulate Systems, School of Materials Science and Engineering, The UniVersity of New South Wales, Sydney NSW 2052, Australia, College of Materials and Metallurgy, Northeastern UniVersity, Shenyang, Liaoning 11004, People’s Republic of China, and CSIRO Minerals, P.O. Box 312, Clayton South, Vic 3169, Australia
Dense-phase pneumatic transportation of bulk materials in the form of slug flow has become a very important technology in industry. In order to understand the underlying mechanisms of slug flow, this paper presents a numerical study of slug flow in horizontal pneumatic conveying by means of discrete particle simulation. To be computationally efficient, the motion of discrete particles is three-dimensional and the flow of continuum gas is two-dimensional, and periodic boundary conditions are applied to both gas and solid phases horizontally. The proposed numerical model is qualitatively verified by comparing the calculated and measured results in terms of particle flow pattern and gas pressure drop. Then the influence of operational conditions such as gas and solid flowrates on slug behavior is numerically investigated. It is shown that slug velocity linearly increases with gas flowrate but is not sensitive to solid flowrate, and slug length increases with both gas and solid flowrates. The results qualitatively agree with the experimental observations. Finally, forces governing the gas-solid flow are analyzed on different length scales. It is shown that the movement of a slug is macroscopically controlled by the axial particle-fluid and particle-wall interactions, whereas the particleparticle interaction microscopically causes a slug to sweep up particles in a settled layer. The magnitudes of these interaction forces increase with gas and solid flowrates. 1. Introduction Low-velocity, dense gas-solid flow in pipelines in the form of solid slug can be found in many industries such as chemical, mineral processing, and agricultural industries to transport bulk materials from one location to another, because of its low product degradation and plant wear, efficient power utilization, and small gas-solid separating unit. Therefore, this flow system has attracted an increasing interest in recent years. This involves the use of various techniques such as physical1-10 and mathematical11-22 modeling. Mathematical modeling is mainly based on the so-called two-fluid model (TFM),20,22 combined computational fluid dynamics (for gas phase) and discrete element model (for particles) (CFD-DEM),11-19 and direct numerical simulation (DNS).21 Slug flow is the main feature in dense gas-solid pneumatic transport and is controlled by particle-particle, particle-wall, and particle-fluid interactions. At this stage of development, experimentation can only access particle-wall interaction in slug flow by means of special stress transducers,7,8,10 and information about particle-particle and particle-fluid interactions is not available. The CFD-DEM approach can overcome this difficulty.11,15,23 The approach has been used in the study of dense pneumatic conveying in recent years. According to the dimensions (i.e., one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D)) used for gas and solid phases, the CFD-DEM models thus-far proposed can be categorized into three groups: 1D CFD + 3D DEM model,11,12,16 2D CFD + 2D DEM model,13-15,17 and 3D CFD + 3D DEM.18,19 While useful, there are various deficiencies associated with the previous CFD-DEM studies. For example, slug flow in a
horizontal pipe is largely 2D because of the gravity effect,24 so a 1D CFD model is not good enough to describe gas flow. In fact, the radial pressure characteristic is missing by this model. On the other hand, the existing 2D DEM models restrict the motions of particles in a 2D duct with arbitrary treatment of porosity,13,17 and in the so-called 3D simulation,18,19 the gas flow is simply described by the Darcy law to avoid CFD effort. Moreover, in order to reduce computational effort, periodic boundary conditions for particles in the flow direction were usually adopted in the study of slug flow, whereas the corresponding gas velocity at the inlet was defined with a fixed profile,11,15,17 or a periodic boundary condition was used in a 1D CFD supposing that gas pressure drop linearly varied along a pipeline.12 However, the profile of inlet velocity would change with the movement of slug once periodic boundary conditions are applied to particles, and pressure drop along conveying pipelines in slug flow is not linear. Another deficiency is that the previous CFD-DEM studies mainly focus on reproducing solid behaviors in slug flow, with limited effort made to analyze the particle-particle and particle-fluid interaction forces that are directly related to the underlying mechanisms. In this paper, a 2D CFD + 3D DEM model, facilitated with suitable periodic boundary conditions for both gas and solid phases, is developed to investigate horizontal slug flow. The model overcomes the problems associated with the previous studies as discussed above. It is applied to study the behavior of particles in slug flow under the experimental conditions similar to those in ref 8. The influence of gas and solid flowrates is also examined. Finally, forces acting on particles are analyzed for better understanding the nature of slug flow. 2. Model Description
* Corresponding author. Tel.: +61 2 93854429. Fax: +61 2 93855856. E-mail:
[email protected]. † The University of New South Wales. ‡ Northeastern University, Shenyang. § CSIRO Minerals.
The current model is an extension of our previous CFD-DEM models in the study of gas-solid flow in fluidization.23,25 For brevity, below we only describe the key feature of this model. However, new developments will be emphasized.
10.1021/ie070991q CCC: $40.75 © 2008 American Chemical Society Published on Web 12/20/2007
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 471 Table 1. Components of Forces and Torques Acting on Particle i forces and torques normal forces
tangential forces
symbol contact
fcn,ij
4 - E*xR*δn3/2nˆ 3
damping
fdn,ij
contact damping
fct,ij fdt,ij
-γn(6mijE*xR*δn)1/2vn,ij
friction forces fluid drag forces gravity torque by tangential forces rolling friction torque
where
equations
-µs|fcn,ij|(1 - (1 - δt,ij/δt,ij,max)3/2)δˆ t (δt,ij < δt,ij,max )
-γt(6µsmij|fcn,ij|x1 - δt,ij/δt,ij,max /δt,ij,max)1/2vt,ij(δt,ij < δt,ij,max ) -µs|fcn,ij|δˆ t (δt,ij > δt,ij,max) 0.5cd0,iFfπRi2|ui - vi|(ui - vi)i1-χ mig Rij × (fct,ij + fdt,ij) µr,ijdi|fn,ij|ω ˆ ˆi
ft,ij ff,i fg,i Tt,ij Tr,ij
ωt,ij δt Ri 2-V 1 1 1 E + , E* ) ,ω ˆ ) , δˆ ) ,δ δ , v ) vj - vi + ωj × Rj - ωi × Ri, vn,ij ) ) ) µs , nˆ ) |Ri| t,ij |ωt,ij| t |δt| t,ij,max 2(1 - V) n ij R* |Ri| |Rj| 2(1 - V2)
(
(vij‚nˆ )nˆ , vt,ij ) vij - vn,ij, cd0,i ) 0.63 +
2.1. Governing Equations for Particle Flow. The particle flow is treated as a discrete phase that is described by DEM. Hence, the translational and rotational motions of particles are determined by Newton’s law of motion, which can be written as
dvi mi
dt
ki
) fp-f,i +
(fc,ij + fd,ij) + mig ∑ j)1
(1)
and
dωi Ii
dt
ki
)
Tij ∑ j)1
(2)
where mi, Ii, vi, and ωi are, respectively, the mass, momentum of rotational inertia, and translational and rotational velocities of particle i; ff-p,i and mig are the force between the particle and fluid and the gravitational force, respectively; and fc,ij, fd,ij, and Ti,j are the contact force, the viscous contact damping force, and the torque between particles i and j, respectively. These interparticle forces and torques are summed over the ki particles in contact with particle i. The particle-particle or particle-wall contact force is calculated according to the equations commonly used in DEM, as recently reviewed by Zhu et al.26 The fluid drag force is achieved according to Di Felice’s correlation,27 which covers a wide range of flow regimes and porosities. The equations used to calculate the forces and torques involved in eqs 1 and 2 are listed in Table 1. 2.2. Governing Equations for Gas Flow. The gas flow is treated as a continuous phase and modeled in a way similar to the one in the conventional two-fluid modeling.28 Thus, its governing equations are the conservation of mass and momentum in terms of local mean variables over a computational cell, given by
∂(Ff) + ∇‚(Ffu) ) 0 ∂t
(3)
and
∂(Ffu) + ∇‚(Ffuu) ) -∇P - Fp-f + ∇‚(τ) + Ffg ∂t
)
[
2FfRii|ui - vi| (1.5 - log10Rep,i) 4.8 2 , Rep,i ) , χ ) 3.7 - 0.65 exp 0.5 ηf 2 Rep,i
(4)
]
2
where Ff, u, and P are, respectively, the fluid density, velocity, and pressure; τ, , and Fp-f are the fluid viscous stress tensor, the porosity, and the volumetric forces between particles and fluid, respectively. 2.3. Coupling between DEM and CFD. The modeling of particle flow by DEM is at an individual particle scale, whereas the gas flow by CFD is at a computational cell scale, which is supposed to be much lager than a particle. Their coupling is numerically achieved as follows. At each time step, DEM will give information, such as the positions and velocities of individual particles, for the evaluation of porosity and volumetric particle-fluid forces in a computational cell. CFD will then use these data to determine the gas flow field, yielding the particle-fluid forces acting on individual particles. Incorporation of the resulting forces into DEM will produce information about the motion of individual particles for the next time step. The particle-fluid forces acting on individual particles will react on the fluid phase so that Newton’s third law of motion is satisfied. This coupling principle has been well-established in the literature.29 Some techniques useful to its implementation in simulation are detailed below. 2.3.1. Property Mapping from Computational Cell Center to Individual Particle. In a CFD cell, there are many particles that may reside at any point of the cell, whereas properties of the cell, such as gas velocity and pressure drop, are usually stored in the center of the cell. To correctly calculate the particle-fluid force acting on a particle, a mapping technique that transforms a cell-based property into a particle-based property is required. Different from our previous bilinear interpolation employed by Feng et al.,25 a generalized interpolation approach, i.e., least-square interpolation,30 is adopted in this work. Hence, a gas property at a particle position can be expressed as
φp ) φcell + ∇φcell‚∆r
(5)
where φp and φcell are the gas properties at the particle position and the cell center, respectively; and ∆r is the distance vector directed from the cell center to the particle position. The gradient of the gas property ∇φcell at the cell center is evaluated by a weighted least-square method,31 and the weighted factor in this work is set to 1/|∆r|2. 2.3.2. Calculation of Porosity and Source Terms. Fluid drag force and the governing equations for gas phase are closely associated with porosity, which is usually calculated at a
472
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
Figure 1. Spheres used in the calculation of porosity over computational cells.
computational cell level. Exact evaluation of porosity in such a cell is sometimes time-consuming, particularly when cells of irregular shape are involved. Numerical instability may occur because of the rapid temporal and spatial change in porosity and associated properties. Some methods have been proposed to overcome these problems.15,32,33 Following the approach used in the local averaging of granular materials,34 we here propose a more general method for this purpose. In this method, a property at any point is calculated using a spherical cell, as schematically shown in Figure 1. The cell size may vary, depending on the need in a study. Therefore, the local porosity corresponding to a particle, usually defined as the porosity value at the center of the computational cell where the particle is residing, can be expressed as Np,sphere
)1-
∑ i)1
asphere,iVp,i
∆Vsphere
(6)
where ∆Vsphere is the volume of a spherical cell; Np,sphere and asphere,i are, respectively, the number of particles and the fractional volume of particle i residing in the spherical cell. Correspondingly, the volumetric particle-fluid interaction force in the considered cell is calculated by Np,sphere
Fp-f )
∑ i)1
asphere,ifp-f,i ∆Vsphere
(7)
Equations 6 and 7 also apply to spherical cells close to the wall boundary, where the region out of the computational domain is not considered.34 At this stage of development, the size of such a spherical cell has to be determined case by case.26,34 In this work, after some trials, the size of a spherical cell is set to 3 particle diameters. 2.4. Boundary Conditions and Numerical Solution. The interaction between a particle and a wall is calculated by the force models between particles, assuming that the wall is a rigid sphere with an infinite diameter, having no displacement and movement attributed by particle-wall interaction. A no-slip condition is applied to the gas phase as used in a CFD study. In the flow direction, periodic boundary conditions for both particle and gas are implemented. Therefore, a particle reenters a pipe with its velocity and radial position being the same as
Figure 2. Schematic of (a) the geometry used in this work and (b) the initial setup of particles in the pipe.
those at the outlet of the pipe. A similar consideration applies to the gas phase, so that
u(x) ) u(x + La)
(8a)
P(x) - P(x + La) ) P(x + La) - P(x + 2La)
(8b)
and
where u and P are the velocity and pressure of gas; x stands for a point in the considered domain; and La is the length of the pipe chosen as the computational domain. Equation 8 can be directly applied to the inlet (x ) 0) and the outlet (x ) La). There are two methods available to numerically implement the boundary condition given by eq 8. One was suggested by Patankar et al.35 In this algorithm, the pressure terms in the gas momentum equations are required to be decomposed into the average pressure gradient and the reduced pressure on condition that average pressure varies linearly along the fluid flow direction; then the average pressure is adjusted in the solution procedure to accommodate algebraic equations. The other method directly implements periodic boundary by mutual replacement of filed variables at the inlet and outlet regions.36-38 Its implementation can be different because of the various periodic lengths under consideration. Considering nonlinear pressure drop in the flow direction of slug flow, the method suggested by Wang and Tao36 is adopted in this work. The method takes one cycle as the solution domain. Its interpolation of fluid velocities at the domain inlet and outlet in the iterations of one time step is modified to take into account the presence of particles in fluid. This gives the following equation, * * M-1 )/2 u11 ) uMM ) (u*2 *2 + uM-1
(9)
where * represents the last iteration; subscripts 1 and 2 stand for the arrangements of computational nodes in the flow direction at the inlet as shown in Figure 1. The same treatment as denoted by M and M - 1 is applied to the node arrangement at the outlet. In addition, the inlet velocity needs to be adjusted according to gas mass conservation right after the velocities at the inlet and outlet are replaced according to eq 9. In the direction of thickness (y-axis), as shown in Figure 2a, one control volume with symmetric boundary condition at the rear and front planes is taken for the gas phase, which allows us to resolve a 2D gas flow field by an existing 3D CFD code. At the same time, periodic boundary conditions
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 473 Table 2. Parameters Used in Present Simulations solid phase
gas phase
materials
polyethylene
shape particle diameter, dp (m) particle density, Fp (kg‚m-3) Young’s modulus, E (Pa) Poisson ratio, ν (Nm-2) sliding friction coefficient, µs rolling friction coefficient, µr damping coefficient, γ particle number, N solid flowrate (kg‚s-1) time step, ∆t (s)
spherical 4 × 10-3 922 1 × 107 0.33 0.3 0.01 0.1 11000a (8000-12000) 0.15a (0.071-0.467) 2.5 × 10-6
a
gas used
air
viscosity, η (kg‚m-1‚s-1) density, Ff (kg‚m-3) pipe geometry CFD cell superficial air velocity, u (m‚s-1) gas flowrate (kg‚s-1) time step, ∆t (s)
width, Lra (m) thickness, Lth (m) length, La (m) width (m) thickness (m) length (m)
1 × 10-5 1.205 0.05 0.013 1.8 0.007 0.013 0.01 2.09a (2.09-3.56) 0.0016a (0.0016-0.0028) 2.5 × 10-6
Base value and, if applicable, varying range in bracket.
are applied to the front and rear walls for the solid phase. This treatment can help produce reasonable porosity and account for 3D particle-particle interaction with minimum computational effort. Facilitated with the above boundary conditions, eqs 1 and 2 for particle flow can be solved by an explicit time-integration method, while eqs 3 and 4 for gas flow can be solved by using the finite volume method with pressure-velocity coupling handled by the well-established SIMPLE algorithm on a nonstaggered grid arrangements.39,40 The discretizations of the convective transport terms and the diffusive transport terms in the gas momentum equations are, respectively, obtained by the QUICK scheme41 and the central difference scheme, while the partial derivative in time is replaced by an implicit Euler approximation. Strongly implicit procedure (SIP)42 then serves as a solver for the discretized governing equations of the gas phase. 2.5. Simulation Conditions. Table 2 lists the parameters used in the present study. Two variables are considered to generate results compared to experimental ones: gas and solid flowrates. The effect of each variable is quantified, while others are fixed as listed in Table 2. The geometry of the conveying pipeline as illustrated in Figure 2a and the type of particles are based on the experimental work of Vasquez et al.8 so that a meaningful comparison between the simulated and experimental results can be made. The damping coefficient used corresponds to the restitution coefficient of 0.8. Before introducing gas at 0.4 s, the pipe is initially filled with a settled layer of particles of height 0.016 m and a slug of length 0.8 m shown in Figure 2b. Note that a side view like Figure 2b will be applied to all figures in the following sections to illustrate the spatial distribution of a considered property. Each simulation lasts for 15 s of physical time, and the simulation results are recorded every 0.015 s. Because macroscopically stable slug flow is considered in this work, the numerical results during the start-up period are not included in the present analysis. For example, the flow after 2.0 s in the case corresponding to Figure 3 is regarded to reach its stable motion because the overall pressure drop of gas in the pipe keeps fluctuating around a constant value. The data before 2.0 s are not included in the analysis. Gas-solid flow in the start-up period will be studied in the future work. 3. Results and Discussion 3.1. Basic Flow Characteristics. For convenience, particles in the pipe at t ) 12.236 s are divided into two groups according to the solid concentration at a given pipe cross section, i.e.,
slug and settled layers, respectively, marked with red and blue color, as illustrated in Figure 4a. The process of particle exchange between the two layers is then traced. It can be seen from parts a-c of Figure 4 that, as a slug moves, it picks up particles in the front settled layer and leaves behind a new settled layer. As a result, particles in the slug keep updating with time. Interestingly, a trapezoidalshape particle configuration colored in blue on the front of the slug, as shown in Figure 4c, is observed. The flow pattern is in line with the video records by Klinzing,44 who conducted experiments under conditions close to the present simulation conditions. It indicates that particles in the center of a settled layer move into the upper part of a slug while particles in the lower part of a settled layer move into the lower area of the slug, as recorded by Takei at al.43 using a high-speed camera. This phenomenon is interesting but has not been reported in the previous CFD-DEM simulations. It has also been observed that some particles are discharged to the gas phase while moving into the upper part of a slug from a settled layer.14,43,44 This phenomenon can also be reproduced in the present CFD-DEM simulations, as seen in these figures. Parts d-i of Figure 4 are the spatial distributions of porosity, particle velocity, gas velocity, gas streamlines, and gas pressure drop along the pipeline, respectively, corresponding to the particle configuration given in Figure 4b. The concept of solid concentration is here the same as packing density ()1-porosity). As seen from Figure 4d, in the axial direction, lower solid concentration emerges on the head and tail of a slug, and in the radial direction, the solid concentration in the center of a slug is denser than that near the wall. In addition, a region with high solid concentration at the bottom of a settled layer in front of a slug is observed and is believed to be caused by the compression from the slug. Although the distribution of solid concentration in the pipe is not uniform, the average solid concentration of a slug over the pipe cross-sectional area fluctuates around a constant, as shown in Figure 5. Furthermore, the mean value of the solid concentration is ∼0.555, very close to 0.54 for a stable slug as experimentally measured by Vasquez et al.8 The result manifests that the current 3D DEM model is capable of reproducing the porosity in the domain of interest although it only employs a relatively small number of particles. Figure 5 also shows that the average solid concentrations are quite different in slug and settled layers. To be quantitative, in the present analysis, solid concentration is averaged over a series of pipe cross sections along the pipeline, and a section is considered to belong to a slug if its average solid concentration is >0.4. This treatment has already been used in Figure 4a to identify the slug and settled layers.
474
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
Figure 3. Variation of overall pressure drop with time.
Parts e and f of Figure 4 show the contours of the radial and axial particle velocities, which are calculated by the following equation,
1
Np,cell
∑ vi
Np,cell i)1
(10)
where Np,cell is the number of particles in a computational cell. It can be seen that the axial particle velocity is relatively uniform within a slug but close to zero within a settled layer. A region with a high particle velocity at the front of the slug is observed. Overall, the radial particle velocity in a slug is much smaller than the axial one. The variation of particle velocity in the whole pipe suggests that particles experience three different stages when passing through a slug, i.e., acceleration, uniform motion, and deceleration. Figure 4g shows that gas mainly flows in the empty part of the pipe over a settled layer before encountering a slug, redistributes itself to cover the entire cross-sectional area at the tail of a slug, and then becomes a partial flow again after passing through the slug. When gas flowrate is low and particles in a settled layer are stationary, a backflow of gas can be seen inside settled layers before and after a slug, as shown in Figure 4h. However, such backflow disappears when particles in a settled layer start to move forward with a small velocity when gas flowrate is higher than a certain value. These results confirm that a fixed velocity profile for gas phase at the inlet of a pipe is incorrect once a periodic boundary condition is applied to particles of slug flow. Figure 4i shows the plot of the pressure drop averaged over the pipe cross section along the pipeline. It can be seen that the overall pressure drop in the pipe is mainly contributed by pressure difference across the slug, as expected. In other words, the overall pressure drop is proportional to slug length and the pressure drop in the space between slugs can be ignored, as suggested by Tomita et al.1 Evidently, this nonlinear pressure drop along the pipeline can justify our choice of method for the implementation of a periodic boundary condition of gas. To examine temporal features of pressure drop, the radial pressure drop at point x ) 1.0 m from the pipe inlet and the axial pressure drop between two points at x ) 1.0 m and x ) 1.04 m are monitored. Figure 5a shows that the axial pressure difference is almost zero when the two monitoring points are located outside the slug but can be ∼600 Pa when they are both inside the slug. This trapezoidal-shape signal of pressure drop qualitatively agrees with the experimental observation of Li et al.6 Similar to the axial pressure drop, the radial pressure drop is significant only when a slug passes the monitored points, as shown in Figure 5b. Gas flow is very sensitive to local porosity or solid concentration, which varies because of the chaotic
motion of individual particles, although, macroscopically, the motion of a slug is stable. This is believed to be the major factor responsible for the gas pressure fluctuations observed in Figure 5a. Note that the above two pressure characteristics have also been reported in the CFD-DEM studies of Li et al.14 and McGlinchey et al.45 These investigators have not used periodic boundary conditions, either for gas or solid phase, in their simulations. This is understandable because periodic boundary conditions can only improve computational efficiency. Once the pipeline as the computational domain is long enough, the inlet conditions should not affect the final outcomes much. 3.2. Influence of Flow Conditions on Slug Behavior. Slug length and slug velocity are two key parameters in the description of slug behavior. As part of this work, the influence of flow conditions, i.e., gas and solid flowrates, on slug length and slug velocity are studied. In this work, the number of particles needs to be given before a simulation starts as a result of the periodic boundary conditions employed. The solid flowrate can then be obtained by counting particles passing through a pipe cross-sectional area. Figure 6a shows the so obtained relationship between the number of particles and the solid flowrate for different gas velocities under the present simulation conditions. Theoretically, the number of particles should be related to gas and solid flowrates. However, such a relationship seems difficult to produce. This may explain why the previous CFDDEM studies11,14,15,17 simply used the concept of particle number rather than solid flowrate. Here, we attempt to produce an equation to estimate this relationship. Since the volume of solids flowing into the system is NpVp and the total volume of the gas and solids is Vpipe, by definition solid concentration p can be calculated by
p )
NpVp Vpipe
(11)
where Np is the particle number. On the other side, solid concentration should be related to gas and solid flowrates wf and wp. By definition, in a time interval ∆t, the volumes of gas and solids flowing into a pipe should be wf∆t/Ff and wp∆t/Fp, respectively. However, the volume occupied by the solids should be larger than wp∆t/Fp, and we assume it can be calculated by multiplying a modifying factor kp. A similar consideration is also applied to the gas phase. In this case, the volumes occupied by gas and solid phases are, respectively, given by kpwp∆t/Fp and kfwf∆t/Ff. Thus, p can be expressed as
p )
kpwp/Fp kfwf/Ff + kpwp/Fp
(12)
Combining eqs 11 and 12 and rearranging the resultant equation gives
Np )
wpFf Vpipe ‚ kwfFp + wpFf Vp
(13a)
where
k ) kf/kp
(13b)
The modifying factors kf and kp are associated with the flow structures of gas and solids in a pneumatic conveying pipe. At this stage of development, their theoretical determination is difficult. By fitting eq 13 to the data in Figure 6a, we found
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 475
Figure 4. Snapshots showing the spatial distributions of flow properties when the number of particles ) 11 000 and the gas velocity ) 2.09 m/s: (a-c) particle configurations at 14.236 s, 14.498 s, and 14.760 s, respectively; (d) porosity; (e) axial particle velocity; (f) radial particle velocity; (g) gas velocity profile; (h) gas streamlines; and (i) axial pressure drop; (d-i) are corresponding to (b).
that the ratio k at a given gas velocity can be approximately treated as a constant and increases with gas velocity described by the following equation:
k ) -0.848 + 0.763u - 0.118u2
(14)
As shown in Figure 6, the two equations can be used to provide a first approximation to the relationship between solid flowrate and the number of particles in CFD-DEM simulation of slug flow. On the other hand, it is noted that the above equations actually relate to a more general question: what is the relationship between solid concentration and operational or flowing conditions? Consequently, they can be linked to the previous studies of pneumatic conveying. For example, by definition, average particle velocity Vjp can be calculated by wp/(pFpA), hence combining with eqs 12 and 14 gives
Vjp ) (-0.848 + 0.763u - 0.118u2)u +
wp FpA
(15)
where A is the cross-sectional area. In the past, many efforts have been made to empirically determine Vjp.46-48 The previous studies are all concerned with dilute flow, where the solid concentration is up to 0.1. The resulting correlations indicate that Vjp should be a function of u, Fp, A (or pipe diameter D),
and other parameters. The present work confirms this is indeed the case. Moreover, it points out that solid flowrate wp is also an important variable for dense flow. At this stage, comparison between the previous and present studies is inadequate because they are dealing with different flow regimes and conditions. However, work has been planned to address this problem, which will be reported in the future publication. Figure 7a presents the variations of slug velocity and slug length under different gas velocities when the number of particles is fixed. In this work, slug velocity is calculated according to the distance a slug has traveled in the whole simulation time for a given case, and slug length is a timeaveraged value determined from the information of solid concentration. It can be seen that, as gas velocity is increased, slug velocity linearly increases, but the slug length rapidly increases first and then slows down. The trends are qualitatively comparable to the experimental observation.8,11 Figure 7b shows the influence of solid flowrate on slug velocity and length when gas velocity is fixed. The results demonstrate that slug length linearly increases with solid flowrate. However, slug velocity is not sensitive to solid flowrate, with only a slight decrease observed. The numerical results again agree qualitatively with the experimental observation.8 The velocity of a slug is determined by the velocities of particles in the slug, and the original driving force of particle
476
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
Figure 6. (a) Dependence of solid flowrate on the number of particles for different gas velocities and (b) relationship between ratio k in eq 13 and gas velocity. Figure 5. Variation of (a) axial pressure drop between x ) 1.0 m and x ) 1.04 m from inlet and (b) radial pressure drop at x ) 1.0 with time.
motion is provided by the gas phase through the particle-fluid interaction. Hence, there should be some relationships between slug velocity, particle velocity, and superficial gas velocity. Following the treatment of Vasquez et al.,8 the ratio of the particle velocity in a slug to the corresponding slug velocity is plotted as a function of superficial gas velocity. Here, the mean velocity of particles in a slug is calculated according to the following expression,
|vjp,slug| )
1
∑
(
1
∆t ∆t Np,slug
)
∑ |vi| i)1
Np,slug
(16)
where |vp,slug| is the magnitude of particle velocity, Np,slug is the number of particles in a slug, and ∆t is the time interval under consideration. Figure 8 shows that the predicted relationship between superficial gas velocity, slug velocity, and particle velocity is in good agreement with the experimental one.8 It can be seen that the ratio of particle velocity to slug velocity infinitely approaches its limit 0.9 under the present experimental conditions as gas velocity is increased. The result confirms that particle velocity is always lower than slug velocity. 3.3. Force Analysis. Analysis of particle-particle, particlewall, and particle-fluid interaction forces acting on particles can generate some insights into particles behavior, hence providing a better understanding of slug flow. Such analysis is difficult to achieve experimentally but can be readily performed based on the CFD-DEM simulations. Parts a and b of Figure 9 show the spatial distributions of particle-fluid interaction forces, corresponding to the particle configuration given in Figure 4b. The magnitudes of these forces
Figure 7. Slug velocity and slug length as a function of (a) gas velocity and (b) solid flowrate.
acting on a particle are expressed relative to its gravity force. This is also applied to other forces presented in this work. It can be seen from Figure 9a that the axial particle-fluid interaction force is larger than particle gravity in a slug but smaller than particle gravity in a settled layer. This is because,
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 477
Figure 8. Variation of particle velocity relative to slug velocity with gas velocity.
when gas passes through a slug with a relatively uniform velocity distribution, it experiences intensive energy exchange with particles through particle-fluid interaction, whereas gas outside a slug mainly flows over a settled layer with little energy exchange between gas and particles. On a slug head, however, the axial particle-fluid interaction force is relatively small because of the low solid concentration, as discussed in Section 3.1. Moreover, the axial particle-fluid interaction force on a slug tail is much smaller than those on other parts of a slug. This is attributed by two factors: the low solid concentration and the change of gas flow direction. Figure 9b shows the spatial distribution of the particle-fluid interaction force in the radial direction. It can be seen that the force is much smaller than the particle gravity. Unlike the axial particle-fluid interaction force whose direction is forward for most of particles, the direction of the radial particle-fluid interaction force is uncertain in the middle of the slug. In general, it is positive on a slug head and negative on a slug tail as gas flow changes its direction in these two regions. The radial particle-fluid interaction force on a slug head is not big enough to pick up particles in the settled layer in front of the slug because its value is much smaller than the gravity of a particle, as given in Figure 9b. In fact, as discussed below, it is particleparticle interaction that drives a slug to sweep up particles in a settled layer. In the present simulation, particle-particle interaction includes the normal and tangential contact forces and the sliding and rolling friction forces. It has been illustrated that, among these forces, the normal contact force is the most important because other contact forces are all related to it.49,50 Therefore, the present analysis is focused on the normal contact forces.
Figure 10. Different average forces on particles in a slug vs time: (a) particle-wall force; (b) particle-particle force; (c) particle-fluid force; (d) overall forces, corresponding to Figure 4.
Figure 9c shows the structure of the contact forces, where each stick represents one connection of two particle centers and its thickness represents the magnitude of the normal contact force. It can be observed that particles experience large normal forces in a slug. The direction of the normal contact forces is almost vertical in the slug except for the front part of a slug where the force direction is inclined to the horizontal direction and directed downward. Evidently, the inclined normal contact forces on a slug head enable the slug to pick up particles in the settled layer in front. On the other hand, it is also found in Figure 9c that the normal contact forces decrease rapidly on a slug tail. In order to investigate forces on particles under different operational conditions, the concept of the average force on particles in a slug is introduced. Because the movement of a slug is dominated by the axial forces acting on particles, only the axial forces are considered in the following. Here, the average particle-fluid force at a given time is calculated by
1 Np,slug
Np,slug
fi
∑ i)1 m
(17)
p,ig
The average particle-particle and particle-wall interaction forces are both given by
1 Np,slug
Np,slug
1
ki
∑ ∑fcn,ij + fdn,ij + fct,ij + fdt,ij i)1 m g j)1 p,i
(
)
(18)
The overall force on a particle is the sum of the above forces. Figure 10 shows the variation of the different average forces
Figure 9. Snapshots showing spatial distributions of (a) axial particle-fluid force; (b) radial particle-fluid force; and (c) network of normal contact forces, corresponding to Figure 4b.
478
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
Figure 11. Average particle-fluid and particle-wall interaction forces on a particle in a slug as a function of (a) gas velocity and (b) solid flowrate.
acting on a particle in a slug with time. It can be seen that the positive particle-fluid force as the driving force of slug flow and the negative particle-wall force as the resistance force balance each other, as expected. As a result, the mean overall force is zero for a stable motion of a slug, if the fluctuation is not considered. Parts a and b of Figure 11 show that the timeaveraged particle-wall and particle-fluid forces are both increased with gas and solid flowrates and the balance of two forces is always maintained. The result agrees with the theory of force balance for a stable slug, which is widely used to predict pressure drop of slug flow in industrial application.2,10,51 Therefore, the discrete particle simulation and continuum modeling are consistent at a slug scale. Note that, according to eq 18, the forces between particles, as internal interactions, interact with each other and can be cancelled out. Figures 10 and 11 are, therefore, the results of slug-scale analysis. To be at a particle scale, the average particle-particle interaction force can be determined as
1
∑
(
1
∆t ∆t Np,slug
Np,slug
1
∑ i)1 m
p,ig
ki
(∑ j)1
|fcn,ij + fdn,ij + fct,ij + fdt,ij|
))
(19)
This force represents the mean total force acting on a particle in a slug. Figure 12 shows that this force increases with both gas velocity and solid flowrate. This is in line with the average particle-fluid interaction force, shown in Figure 11, which is also at a particle scale according to eq 17. The results indicate that increasing gas velocity does not change the solid concentration in a slug much but increases the relative velocity between gas and particles. This will increase the particle-fluid interaction force and correspondingly the particle-wall interaction force. On the other hand, increasing solid flowrate can increase the solid concentration in a slug or the computational domain as a
Figure 12. Average particle-particle interaction force and solid concentration in a slug as a function of (a) gas velocity and (b) solid flowrate.
whole. This can enhance the particle-particle interaction and, at the same time, increase the particle-fluid interaction force, finally leading to a stronger particle-wall interaction. Therefore, the force results at different length scales are interrelated. In particular, the results on a particle scale well-explain the experimental observations on a macroscopic, slug scale. 4. Conclusions A 2D CFD and 3D DEM model has been developed to simulate the slug flow in horizontal pneumatic conveying. Some new modeling techniques have successfully been implemented in the developed model to improve numerical efficiency and stability. These include: (1) suitable periodic boundary conditions for both gas and solid phases in the flow direction, (2) the least-square interpolation of gas properties, and (3) the concept of spherical cell in the calculation of porosity and particle-fluid interaction force. The numerical results indicate that the current model can capture the flow behaviors observed in experiments, such as the exchange of particles between a slug and a settled layer, the gas pressure characteristic, and the interrelationship between gas velocity, slug velocity, and particle velocity. They also show that slug velocity linearly increases with gas flowrate but is not sensitive to solid flowrate and slug length increases with both gas and solid flowrates. An empirical correlation is proposed to approximate the relationship between the number of particles or solid concentration and solid flowrate for different gas flowrates, which is useful to relate CFD-DEM simulations of slug flow to engineering application. Forces governing the slug flow in pneumatic conveying have been analyzed at different length scales. It is shown that the axial particle-fluid force in slug flow is much bigger than the radial one. The movement of a slug is macroscopically controlled by the axial particle-fluid and particle-wall interac-
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008 479
tions, whereas the particle-particle interaction microscopically causes a slug to sweep up particles in a settled layer. The magnitudes of these interaction forces increase with gas and solid flowrates. Acknowledgment The authors are grateful to the Australia Research Council (ARC) for the financial support of this work.
n ) normal component p ) particle p-f ) particle-fluid r ) rolling friction ra ) radial direction s ) sliding sphere ) spherical cell used in the calculation of porosity t ) tangential component th ) thickness of pipe
List of Symbols a ) fractional volume of a particle A ) cross-sectional area of pipe, m2 cd0 ) fluid drag coefficient on an isolated particle d ) particle diameter, m D ) pipe diameter, m E ) Young’s modulus, Pa f ) particle scale forces, N F ) volumetric force, N‚m-3 g ) gravitational acceleration, m‚s-2 I ) moment of inertia of particle, kg‚m2 k ) empirical factor ki ) number of particles in contact with particle i L ) pipe geometry, m m ) mass of particle, kg N ) number of particles P ) pressure, Pa ∆r ) vector from cell center to particle center position, m R ) vector from the mass center of the particle to the contact point, m Rep,i ) Reynolds number of particle i t ) time, s ∆t ) time step, s T ) torque, N‚m u ) gas velocity, m‚s-1 v ) particle translational velocity, m‚s-1 V ) volume, m3 w ) mass flowrate, kg‚s-1 x ) position vector in the computational domain, m Greek δt ) vector of the accumulated tangential displacement between particles i and j φ ) general variable ) porosity γ ) damping coefficient η ) gas viscosity, kg‚m-1‚s-1 F ) density, kg‚m-3 ω ) particle angular velocity, s-1 ω ˆ ) unit vector defined by ω ˆ ) ω/|ω| µ ) friction coefficient ν ) Poisson ratio τ ) fluid viscous stress tensor, kg‚m-1‚s-2 Subscripts 1, 2, 3, M ) index of computation cell a ) axial direction c ) contact cell ) computational cell d ) damping f ) fluid i ) particle i ij ) between particles i and j j ) particle j
Literature Cited (1) Tomita, Y.; Jotaki, T.; Hayashi, H. Wave-Like Motion of Particulate Slugs in a Horizontal Pneumatic Pipeline. Int. J. Multiphase Flow 1981, 7, 151. (2) Konrad, K. Dense-Phase Pneumatic ConveyingsA Review. Powder Technol. 1986, 49, 1. (3) Mi, B.; Wypch, P. W. Pressure-Drop Prediction in Low-Velocity Pneumatic Conveying. Powder Technol. 1994, 81, 125. (4) Molerus, O. Overview: Pneumatic transport of solids. Powder Technol. 1996, 88, 309. (5) Jaworski, A. J.; Dyakowski, T. Investigations of flow instabilities within the dense pneumatic conveying system. Powder Technol. 2002, 125, 279. (6) Li, J.; Pandiella, S. S.; Webb, C.; McGlinchey, D.; Cowell, A.; Xiang, J.; Knight, L.; Pugh, J. An experimental technique for the analysis of slug flows in pneumatic pipelines using pressure measurements. Particulate Sci. Technol. 2002, 20, 283. (7) Wypych, P. W.; Yi, J. Prediction of pressure drop for low-velocity pneumatic conveying. In Proceedings of the 4th World Congress of Particle Technology, Sydney, Australia, 2002. (8) Vasquez, N.; Sanchez, L.; Klinzing, G. E.; Dhodapkar, S. Friction measurement in dense phase plug flow analysis. Powder Technol. 2003, 137, 167. (9) Li, J.; Pandiella, S. S.; Webb, C.; Dyakowski, T.; Jones, M. G. Analysis of gas-solids feeding and slug formation in low-velocity pneumatic conveying. Particulate Sci. Technol. 2003, 21, 57. (10) Niederreiter, G.; Sommer, K. Modeling and experimental validation of pressure drop for pneumatic plug conveying. Granular Matter 2004, 6, 179. (11) Tsuji, Y.; Tanaka, T.; Ishida, T. Lagrangian Numerical-Simulation of Plug Flow of Cohesionless Particles in a Horizontal Pipe. Powder Technol. 1992, 71, 239. (12) Kawaguchi, T.; Tanaka, T.; Tsuji, Y. Numerical analysis of density wave in dense gas-solid flows in a vertical pipe. Prog. Theor. Phys. Suppl. 2000, 696. (13) Xiang, J. S.; McGlinchey, D. Numerical simulation of particle motion in dense phase pneumatic conveying. Granular Matter 2004, 6, 167. (14) Li, J.; Webb, C.; Pandiella, S. S.; Campbell, G. M.; Dyakowski, T.; Cowell, A.; McGlinchey, D. Solids deposition in low-velocity slug flow pneumatic conveying. Chem. Eng. Process. 2005, 44, 167. (15) Lim, E. W. C.; Wang, C. H.; Yu, A. B. Discrete element simulation for pneumatic conveying of granular material. AIChE J. 2006, 52, 496. (16) Fraige, F. Y.; Langston, P. A. Horizontal pneumatic conveying: A 3d distinct element model. Granular Matter 2006, 8, 67. (17) Ouyang, J.; Yu, A. B. Simulation of gas-solid flow in vertical pipe by hard-sphere model. Particulate Sci. Technol. 2005, 23, 47. (18) Strauss, M.; McNamara, S.; Herrmann, H. J.; Niederreiter, G.; Sommer, K. Plug conveying in a vertical tube. Powder Technol. 2006, 162, 16. (19) Strauss, M.; McNamara, S.; Herrmann, H. J. Plug conveying in a horizontal tube. Granular Matter 2007, 9, 35. (20) Lee, L. Y.; Quek, T. Y.; Deng, R. S.; Ray, M. B.; Wang, C. H. Pneumatic transport of granular materials through a 90 degrees bend. Chem. Eng. Sci. 2004, 59, 4637. (21) Singh, P.; Wang, A.; Rosato, A. D. Direct Simulation of Pressure Driven Flows of Dense Powders. In 14th U.S. Congress on Theoretical & Applied Mechanics, Blacksburg, VA, 2002. (22) Levy, A. Two-fluid approach for plug flow simulations in horizontal pneumatic conveying. Powder Technol. 2000, 112, 263. (23) Feng, Y. Q.; Yu, A. B. Microdynamic modelling and analysis of the mixing and segregation of binary mixtures of particles in gas fluidization. Chem. Eng. Sci. 2007, 62, 256. (24) Fan, L. S.; Chao, Z. Principles of Gas-Solid Flows; Cambridge University Press: New York, 1998.
480
Ind. Eng. Chem. Res., Vol. 47, No. 2, 2008
(25) Feng, Y. Q.; Xu, B. H.; Zhang, S. J.; Yu, A. B.; Zulli, P. Discrete particle simulation of gas fluidization of particle mixtures. AIChE J. 2004, 50, 1713. (26) Zhu, H. P.; Zhou, Z. Y.; Yang, R. Y.; Yu, A. B. Discrete particle simulation of particulate systems: Theoretical developments. Chem. Eng. Sci. 2007, 62, 3378. (27) Di Felice, R. The Voidage Function for Fluid Particle Interaction Systems. Int. J. Multiphase Flow 1994, 20, 153. (28) Gidaspow, D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions; Academic Press: Boston, MA, 1994. (29) Xu, B. H.; Yu, A. B. Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chem. Eng. Sci. 1997, 52, 2785. (30) Imaichi, K.; Ohmi, K. Numerical Processing of Flow-Visualization Pictures Measurement of Two-Dimensional Vortex Flow. J. Fluid Mech. 1983, 129, 283. (31) Jayantha, P. A.; Turner, I. W. A comparison of gradient approximations for use in finite-volume computational models for two-dimensional diffusion equations. Numer. Heat Transfer, Part B 2001, 40, 367. (32) Wu, C. L.; Zhan, J. M.; Li, Y. S.; Lam, K. S. Dense particulate flow model on unstructured mesh. Chem. Eng. Sci. 2006, 61, 5726. (33) Link, J. M.; Cuypers, L. A.; Deen, N. G.; Kuipers, J. A. M. Flow regimes in a spout-fluid bed: A combined experimental and simulation study. Chem. Eng. Sci. 2005, 60, 3425. (34) Zhu, H. P.; Yu, A. B. Averaging method of granular materials. Phys. ReV. E 2002, 66. (35) Patankar, S. V.; Liu, C. H.; Sparrow, E. M. Fully Developed Flow and Heat-Transfer in Ducts Having Streamwise-Periodic Variations of CrossSectional Area. J. Heat Transfer Trans. ASME 1977, 99, 180. (36) Wang, L. B.; Tao, W. Q. Heat-Transfer and Fluid-Flow Characteristics of Plate-Array Aligned at Angles to the Flow Direction. Int. J. Heat Mass Transfer 1995, 38, 3053. (37) Bai, L.; Mitra, N. K.; Fiebig, M.; Kost, A. A Multigrid Method for Predicting Periodically Fully-Developed Flow. Int. J. Numer. Methods Fluids 1994, 18, 843. (38) Amano, R. S. A Numerical Study of Laminar and Turbulent HeatTransfer in a Periodically Corrugated Wall Channel. J. Heat Transfer Trans. ASME 1985, 107, 564. (39) Ferziger, J.; Peric, M. Computational Methods for Fluid Dynamics, 3rd ed.; Springer: New York, 2002.
(40) Partankar, S. V. Numerical Heat Transfer and Fluid Flow; Hemisphere Publishing Co.: Washington, DC, 1980. (41) Hayase, T.; Humphrey, J. A. C.; Greif, R. A Consistently Formulated Quick Scheme for Fast and Stable Convergence Using FiniteVolume Iterative Calculation Procedures. J. Comput. Phys. 1992, 98, 108. (42) Stone, H. L. Iterative Solution of Implicit Approximations of Multidimensional Partial Differential Equations. Siam J. Numer. Anal. 1968, 5, 530. (43) Takei, M.; Ochi, M.; Saito, Y. Image extraction of particle concentration at the plug front using 3D wavelets and comparison with LDV. Powder Technol. 2004, 142, 70. (44) Klinzing, G. E. Dense Phase (Plug) ConveyingsObservations and Projections. In Handbook of ConVeying and Handling of Particlulate Solids; Elsevier Science: New York 2001. (45) McGlinchey, D.; Cowell, A.; Pugh, J. R.; Knight, E. A.; Xiang, J.; Li, J. Axial and radial pressure drops of dense phase plugs. Particulate Sci. Technol. 2005, 23, 215. (46) Plasynski, S. I.; Klinzing, G. E.; Mathur, M. P. High-Pressure Vertical Pneumatic Transport Investigation. Powder Technol. 1994, 79, 95. (47) Klinzing, G. E.; Myler, C. A.; Zaltash, A.; Dhodapkar, S. A Simplified Correlation for Solids Friction Factor in Horizontal Conveying Systems Based on Yangs Unified Theory. Powder Technol. 1989, 58, 187. (48) Matsumoto, S.; Harakawa, H.; Suzuki, M.; Ohtani, S. Solid ParticleVelocity in Vertical Gaseous Suspension Flows. Int. J. Multiphase Flow 1986, 12, 445. (49) Zhou, Y. C.; Yu, A. B.; Stewart, R. L.; Bridgwater, J. Microdynamic analysis of the particle flow in a cylindrical bladed mixer. Chem. Eng. Sci. 2004, 59, 1343. (50) Jayasundara, C. T.; Yang, R. Y.; Yu, A. B.; Curry, D. Discrete particle simulation of particle flow in the IsaMill process. Ind. Eng. Chem. Res. 2006, 45, 6349. (51) Pan, R.; Wypych, P. W. Pressure drop and slug velocity in lowvelocity pneumatic conveying of bulk solids. Powder Technol. 1997, 94, 123.
ReceiVed for reView July 20, 2007 ReVised manuscript receiVed October 6, 2007 Accepted October 11, 2007 IE070991Q
