- Email: [email protected]
Numerical analysis of particle mixing characteristics in a single helical ribbon agitator using DEM simulation
Powder Technology 108 Ž2000. 55–64 www.elsevier.comrlocaterpowtec
Numerical analysis of particle mixing characteristics in a single helical ribbon agitator using DEM simulation Yasunobu Kaneko a
a,)
, Takeo Shiojima a , Masayuki Horio
b,1
Production Technology Center, Idemitsu Petrochemical Co. Ltd., Ichihara, 1-1 Anesaki-kaigan, Ichihara, Chiba 299-0193, Japan b Department of Chemical Engineering, Tokyo UniÕersity of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan Received 24 May 1999; received in revised form 11 October 1999; accepted 11 October 1999
Abstract Numerical analysis of three-dimensional motion of particles in a single helical ribbon agitator was carried out by means of the Discrete Element Method ŽDEM.. To validate the computed results experiments were carried out with a cold scale model of 0.3 m inside diameter. Circulation time of particles in the agitator and the horizontal particle velocity distribution in the core region predicted by the present simulation agreed well with those obtained by experiments. Based on DEM simulation, the particle circulation and mixing characteristics in the agitator vessel were investigated in detail. Vertical mixing of particles was found rather poor during upward and downward flows through the blade and core regions, respectively. However, in the core region, particles flow in the manner of funnel flow. Additionally, in the blade region the vertical velocity decreases toward the wall and reduces to zero near the boundary between the blade and core regions. Particles were carried up to the bed surface by the helical ribbon flow down into the center of the core, and those not carried completely to the surface flow into the outer periphery of the core. The effect of bed height on circulation and mixing was analyzed numerically in terms of tracer particle concentration response as well as particle velocity and angular velocity distributions, which are quite difficult to obtain by experiment. The bed height was found to be the most important factor for mixing and circulation. When the bed height is appropriate, particles circulate rather regularly. If the bed height is lower than the blade top, particles mix rapidly because of the high fall of particles from the blade top. On the other hand, if the bed height is higher than the blade top, the velocity of particles in the bed surface region is much lower than that of the other region and accordingly a semi-stagnant eddy is formed. q 2000 Elsevier Science S.A. All rights reserved. Keywords: Discrete element method; Particle mixing; Helical ribbon agitator
1. Introduction Helical ribbon agitators have been widely used for mixing liquids of high viscosity. There are many reports providing basic information of helical ribbon agitators concerning power consumption of agitators, flow pattern andror mixing characteristics w1–4x. This type of agitators has also been applied to blend powders. However, there have been only few published reports on helical ribbon agitation of powders. Moreover, their literature focus has been only on the overall mixing characteristics and power consumption for agitation. ) Corresponding author. Tel.: q81-436-60-1850; fax: q81-436-601121; e-mail: [email protected] 1 E-mail: [email protected].
Fig. 1 illustrates a single helical ribbon agitator and a typical flow pattern of particles in the vessel. Particles are scraped by the inclined blade, and lifted up in the blade or upflow region by the ribbon. This is because particles cannot rotate at the same speed as the ribbon’s due to the friction between particles and the wall. At the top of the bed the particles fall down into the core or downflow region of the bed. Particles are thus circulated in the axial direction and mixed by the particle exchange between the core and the blade regions. Dead zone is scarcely formed and the power consumption is rather low. Accordingly, in the field of gas phase propylene polymerization reactors helical ribbon agitators have been applied successfully w5x. However, since polymerization is exothermic with a rather high heat of reaction and since the difference between the reaction temperature and melting temperature of
0032-5910r00r$ - see front matter q 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 3 2 - 5 9 1 0 Ž 9 9 . 0 0 2 5 1 - X
56
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
a rotary vessel mixer w18x and in a centrifugal tumbling granulator w19,20x. Wightman et al. w21x and Sakamoto w22x investigated the motion of particles in a rotating drum, and Matsusaka et al. w23x and Shimosaka et al. w24x applied DEM simulation to estimate particle behavior under vibrating conditions. However, the DEM simulation of particle motion in the vessel with an agitator as well as with the helical ribbon agitator has not been carried out yet. In the present study, an attempt was made to develop a three-dimensional DEM code to simulate particle motion in a helical ribbon agitator vessel. To evaluate the accuracy of the calculated results, experiments were carried out in a mixer cold test rig with the same geometry and size. Fig. 1. Schematic of a helical ribbon blade and flow pattern in a single ribbon agitator.
2. Theoretical model polypropylene is small, the temperature of polymer particles tends to rise and exceeds the melting temperature. If the heat of reaction generated in each particle is not effectively removed, the polymer particles tend to melt and form lumps, which may hinder the stable operation of the reactor. Recent progress in catalyst activity has been quite remarkable. On the other hand, the demand for random copolymerization with a rather high ethylene content has been increasing to lower the melting temperature for easier film production. Since these reduce the difference between reaction temperature and the melting temperature mentioned above, the reactor control for safe operation is becoming more important. Accordingly, more attention should be focused on the study of detailed particle mixing behavior. In the literature, Reichert et al. w6x studied the mechanical aspects, such as power consumption of the agitator, shape of the vortex and circulation time for a large-scale agitator design presenting a dimensionless correlation. Cooker and Nedderman w7x and Cooker et al. w8x studied the probability distribution of circulation times by means of a radio pill tracking system. They further proposed a theory to predict the power consumption for agitation and the circulation rate of particles based on stress analysis with appropriate modifications of Janssen’s method w9,10x. However, no reports can be found concerning the detailed particle mixing behavior in the helical ribbon agitator vessel. This is probably because of the difficulties in both experiment and mathematical analysis to pursue the motion of particles. The recent rapid growth of numerical power in the field of particle technology is now becoming quite promising to obtain valuable information on a variety of powder processes. Discrete element method ŽDEM. with soft sphere model initially proposed by Cundall and Strack w11x for quasi-static deformation of a particle bed has been successfully applied to analyze particle motions in densely packed beds such as fluidized beds w12–15x, andror hoppers w16,17x. In the field of mechanical particle agitation, Muguruma et al. applied DEM to simulate particle flows in
2.1. Assumption A mathematical model was developed based on the following assumptions for simplification: 1. Particles are spherical in shape and uniform in diameter and move in a three dimensional space. 2. The gas-to-particle interaction and the particle-to-particle adhesion are neglected. Forces considered are those for gravity and collisions between particles, between particle and wall and between particle and blade. 3. Particles are supposed to have a soft sphere interaction expressed by a Hookean linear spring, a dash pot and a friction slider with Coulomb’s law of friction. 4. Based on justifications in previous work w13,25x, the spring constant is adjusted only from the viewpoint of numerical economy. 2.2. Equation of particle motion Both translational and rotational motions are considered in equations of motion. The equations for individual particles can be written as follows: x¨ s Ž Fn H Ft . rm q g
Ž 1.
v ˙ s TcrI
Ž 2.
The above equations are the same as those of Tsuji et al. w25x, where x is the position vector of the center of the particle, Fn and Ft are the sum of the collision forces respectively for the normal and tangential directions acting on colliding particles, m is the mass of particles, g is the gravity acceleration, v is the angular velocity vector of a particle, Tc is the summation of torque caused by particle collision and I is the moment of inertia. The Euler scheme is used to integrate Eqs. Ž1. and Ž2. to obtain the new velocity, position and angular velocity of a particle after a time step. Fig. 2 illustrates the definitions of variables under the colliding condition of particles. When particle i is in
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
57
zr is the velocity vector of particle i relative to particle j. Subscripts ‘n’ and ‘t’ denote the normal and tangential directions and subscripts i and j denote the particle i and j, respectively. The direction of Fn is the one from point c to the center of particle j, and the direction of Ft is the same as the one of the grade crossing between the surface including the vector zr and the vector zn and the surface of solid crossing of both particles.
2.3. Description of the blade and the judgment of particleto-blade collisions The following equations with a parameter u are used to express the spiral curve of a blade shown in Fig. 3 with radius R and pitch Hp :
Fig. 2. Schematic illustration of particle collision.
contact with particle j at point x c , the normal and tangential components of the collision force acting on particle i are given respectively by: Fn s yk dn y h zn
Ž 3. Ž 4a .
Ft s yk dt y h zt < Ft < F m < Fn < zt < Ft < ) m < Fn < Ft s ym < Fn < < zt <
Ž 4b .
where zn s Ž zr P n . n zt s zr = n zr s zj y zi q Ž v i q v j . = Ž x i y x c .
Ž 5. Ž 6. Ž 7.
where k is the spring constant, d is the particle displacement, h is the damping coefficient, z is the particle velocity vector, m is the friction coefficient, n is the unit vector from the center of particle i to that of particle j and
x s Rcos u y s Rsin u z s bu Ž b ' Hpr2p .
Ž 8.
A helical ribbon blade can be represented by changing R from the inner radius of the blade to its outer radius. Point Ž x b , y b , z b . on the blade, where there is a possibility of contact with particle i having its center at Ž x i , yi , z i ., is expressed by the following equations: x b s ri cos u y b s ri sin u zb s bŽ u y v b t . ri ' Ž x i2 q yi2 .
Ž 9.
0.5
where v b is the angular velocity of the agitator rotation.
Fig. 3. Mathematical expression of the single helical agitator.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
58
The distance l between particle i Ž x i , yi , z i . and a point on the blade Ž x b , y b , z b . is obtained as: 2
2
l 2 s Ž x i y x b . q Ž yi y y b . q Ž z i y z b . 2
s Ž x i y ri cos u . q Ž yi y ri sin u . qŽ zi y b Ž u y v b t . .
2
2
2
Ž 10 .
Differentiating Eq. Ž10. by u , we obtain: dŽ l 2 . du
s 2 ri Ž x i sin u y yi cos u . y 2 b Ž z i y b Ž u y v b t . .
Ž 11 . Setting dŽ l 2 .rd u s 0 and solving Eq. Ž11. for u , we obtain l min , the minimum distance between particle i and the blade, from Eq. Ž9.. In the present computation the Newton–Raphson method was used for iteration. If l min is smaller than particle radius, it is judged that the particle i and the blade are in collision.
Table 1 Geometry of the agitator vessel and property of particles Vessel Height Inside diameter Blade Total height Inside diameter Outside diameter Pitch Height of plough Particle Mean diameter Apparent density
Hv Dv
0.7 m 0.31 m
Hb Di Do Hp hp
0.495 m 0.18 m 0.28 m 0.11 m 0.08 m
dp ra
1.0 mm 550 kgrm3
distribution. In all experiments the level of the powder bed was adjusted to the top edge of the agitator blade for each mixing condition. 3.2. Circulation time
3. Experimental 3.1. Apparatus The single helical ribbon agitator vessel used for the present experiment is illustrated in Fig. 4. Table 1 shows the geometry of the agitator vessel and the properties of particles. The motor and converter assembly was able to successively vary the rotation speed of the agitator up to 150 rpm. The powder was polypropylene homopolymer having approximately 1 mm in diameter with a narrow size
One of the most important factors in estimating the performance of an agitator is the average circulation time. To obtain the circulation time the following two methods were adopted: Ž1. Particles are agitated with a colored rubber ball having 20 mm in diameter. If the ball reached near the bed surface, the ball was able to be distinguished from the polypropylene particles because they were white. The time required for the ball to emerge again near the bed surface was measured more than fifty times and the particle circulation time was determined. Ž2. One kilogram of polypropylene pellets of about 3 mm in diameter, i.e., the weight ratio of the pellets to the whole bed weight was about 0.05, were poured into the center of the bed surface as tracer particles. A small amount of the bed material was collected from the center of the bed surface every 2 s for about 1 min. The weight concentration of the pellets in the sampled particles was obtained by sieving. From the oscillatory plot of pellet concentration against elapsed time, the time interval between the neighboring peaks was determined as the particle circulation time. Simulation was done for conditions similar to method Ž2. of the above two experiments. 3.3. Horizontal distribution of Õertical Õelocity of particles in the core region
Fig. 4. Experimental apparatus.
Fig. 5 shows the set-up for measurement of vertical velocity of particles in the core region. A thin glass tube of 10 mm inside diameter was immersed in the bed. A funnel was installed at the top of the tube to easily pour particles into the glass tube. The tube was perpendicularly fixed by the tube guide which was fastened on the base plate with bolts so that it can be horizontally moved. Polypropylene particles, the same as the bed material, were filled up in the tube before the start of agitation.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
59
8 mm. For this computation condition the CPU time required for computation of real 1 s was approximately 8 h. The time step, an important factor for the stable calculation, was determined by the following equation proposed by Tsuji et al. w25x: DtF
2p 5
(
m
Ž 12 .
k
In this study the spring constant, k, was assumed to be 800 Nrm, so that the time step for integration D t was as large as 1.0 = 10y4 s.
5. Results and discussion 5.1. Comparison between the results of experiment and calculation
Fig. 5. Set-up to determine the vertical velocity of particles in the core region.
Right after the start, particles in the glass tube also start to move down at a constant rate, which is supposed to correspond to a vertical flow rate of particles at a horizontal position of the bed Žat the tube exit.. From the descending rate of the upper surface of particles in the straight part of the tube, the vertical rate of particles was determined. The same procedure was repeated over different horizontal positions, but at a fixed height Ž0.4 m from the bottom. in the bed. Another tube of 20 mm I.D. was also applied but no appreciable difference was observed.
4. Computation conditions Computation conditions and parameters used in the present study are shown in Table 2. The column size and the agitator dimensions were identical except for geometry of the vessel bottom, i.e., in the computation a flat bottom was assumed for simplification. Computation was performed with a personal computer ŽFujitsu, FMV-6266T6.. The particle number in the present experiment is roughly 4.0 = 10 7, for which the numerical simulation with the present DEM model cannot be practically carried out. To reduce the computational load the number of particles was reduced to 1.0 = 10 5 by increasing particle diameter up to
Fig. 6 shows the typical snapshot of the mixing behavior under the condition of 3.0p radrs of the agitation speed. The black and gray particles indicate having upward and downward velocity vectors, respectively. As can be clearly seen in the numerical results of Fig. 6, particles in a helical ribbon agitator circulate by being lifted up along the wall by the rotating blade and descending in the core. The revolution of powder sinking at the center causes a concave bed surface. The vortex depth, i.e., the difference of height between the bed surface at the center and that at the wall, generally increased with the rotating speed. Fig. 7 shows the snapshots of mixing behavior at vertical cross section for three rotational speeds. From this figure, it is shown that the maximum vortex depth gradually increases with the rotating speed. This result is in agreement with the general observation of the vortex profile of helical ribbon agitators. To validate the computed particle circulation times they were compared with the experimental values. Computed circulation time was obtained by the following procedure; Ži. After reaching a stable mixing, particles located within 0.1 m from the top surface of the core region shown in Fig. 8 were defined as tracer particles. The above region was also defined as the inspection volume. Žii. The concentration of tracer particles in the inspection volume was
Table 2 Computation conditions Number of particles Particle diameter Spring constant Damping coefficient Friction coefficient Particle density Time step Rotational speed
np dp k h m rp Dt vb
7.0, 9.0, 10.0=10 4 8 mm 800 Nrm 0.6 0.3 800 kgrm3 1.0=10y4 s 2.0p , 3.0p , 4.0p radrs
60
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
Fig. 6. Typical snapshot of computed particle behavior at t s 9.9 s.
calculated as a function of time. Žiii. There are several peaks in a concentration vs. time curve. The time interval between the adjacent peaks is defined as the particle circulation time as shown in the lower half of Fig. 8. Fig. 9 shows that the particle circulation time as a function of rotation speed. Experimental values were obtained by the one ball tracer test and particle tracer test. As can be clearly seen in Fig. 9, the computed circulation time agreed well with that of the experiments.
The responses of tracer particle concentration for both experiment and calculation are compared in Fig. 10. The period of the initial low concentration corresponds to the time of first circulation. The time at which the tracer particles enter again into the inspection volume is approximately same for both cases. In Fig. 11, the computed horizontal distribution of particle velocity Žvertical component. is compared with that from the experiment. As defined in Fig. 11, the
Fig. 7. Effect of rotational speed on particle behavior.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
61
Fig. 8. Method of determining the circulation time from calculation.
inspection volume, i.e., the region between 0.4 and 0.45 m from the bottom, was divided coaxially into 15 sections Žnine for the core and six for the blade region.. The computed values were obtained by averaging the vertical velocities of particles in each section. In the previous investigations w6,9,26x, the flow pattern of the core region was treated as the plug flow. However, the present result both from experiment and computation indicates that the flow in the core is a sort of funnel flow, having a nonuniform velocity distribution. In addition, in the blade region the vertical velocity decreases toward the wall because a finite clearance existed between the blade tip and the wall. On the other hand, getting closer to the boundary between the blade and the core regions, the vertical velocity component reduces to zero due to the momentum exchange between the two regions.
5.2. Detail particle flow and the effect of bed height on particle motion
Fig. 9. Model validation in terms of circulation time.
Fig. 10. Tracer concentration vs. time.
Fig. 12 shows the snapshots of particle locations, where particles are painted in different gray levels depending on their initial locations. From this figure it can be seen that vertical mixing of particles is rather poor during upward and downward motion through the blade and core regions, respectively. Particle mixing takes place near the bed
62
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
Fig. 11. Vertical velocity distribution.
surface by the particle flow into the cavity of the core top. Moreover, from Fig. 12Žd. the particles located at the middle height on the blade at t s 0 s flow through the outer side of the core region. Particles that were lifted not completely up to the bed surface tend to descend through the outside section of the core. Fig. 13 shows the snapshots of particle motion for three different bed heights. The shape of vortex is quite different for different bed heights. The effect of bed height on particle mixing seems considerable even though the rotational speed is the same. To make this point clearer, the mixing process is examined by tracking the tracer particles at three different regions, i.e., the top region of the bed, the
middle region of the core section and the bed bottom region. Fig. 14 shows the concentration of tracer particles in three regions Žor inspection volumes. for three cases of bed heights. Regions 2 and 3 have the same heights for all three cases, but the definition of region 1 is different for each case so that the tracer concentration in the surface region can be examined. The initial position of tracer particles, inspection volumes and particle velocity vectors in x–z plane are also depicted in the upper half of Fig. 14. In case Žb., a rather sharp peak can be found in all three regions and intervals between adjacent peaks are the same for all regions indicating the solid circulation is quite through from the bottom to the top. In case Ža. with a low bed height, solid circulation only took place between regions 2 and 3. In the top region the tracer particle concentration did not change after some of them were carried up from region 2. Furthermore, tracer concentration in all
Fig. 12. Mixing process of particles.
Fig. 13. Effect of bed height on particle mixing.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
63
Fig. 14. Effect of bed height on particle mixing process.
regions approached to an equilibrium value rather quickly. The higher vertical distance between the blade top and the core surface should have helped the tracer mixing. When the bed surface is much higher than the blade top Žcase Žc.., the particle velocity near the bed surface region is rather low as shown in the upper half of Fig. 14. This indicates that the near surface region, where no forced circulation exists, acts as a semi-stagnant eddy region. The formation of this semi-stagnant eddy region helps particles
mixing by exchange but disturbs the clear circulation. Fig. 15 shows the angular velocity of particles both in the core and blade regions for cases Žb. and Žc. of Fig. 14. By experiment the angular velocity in the blade region usually cannot be easily obtained. This indicates a clear advantage of DEM simulation. Obviously, the near surface region of case Žc. rotates much more slowly than the lower section. Furthermore, the rotation up to the height of blade top is almost the same for different bed heights. Bed height is clearly one of the most important operating parameters for the helical ribbon agitators.
6. Conclusions
Fig. 15. Vertical profiles of angular velocity Ž v b s 4p wradrsx..
Three-dimensional motion of particles mixed in a single helical ribbon agitator vessel was numerically analyzed by the DEM simulation. To reduce the computation load, a spring constant smaller than the reality and larger particles were adopted for computation. Experiments were conducted and computed results were successfully compared with experimental data. By tracking particle motion in detail, the circulation and mixing characteristics were analyzed in terms of tracer particle concentration response, particle trajectories and
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
64
angular velocity distribution, some of which are quite difficult to obtain by experiment. It was demonstrated that the distance between the blade top and the bed surface is the most important parameter to attain good mixing or good circulation. When bed surface is lower than the blade top, particles lifted to the blade top fall into the center of the core cavity, enhancing rapid mixing. When bed surface is higher than the blade top, a semi-stagnant eddy was formed on the bed surface. If the bed height is appropriate, a proper circulation throughout the bed height is created.
7. List of symbols dp Dt Do Dv Fn Ft g hp Hb Hp Hv I k l m n np r t Dt Tc z zr x d h m
Particle diameter wmx Inside diameter of blade wmx Outside diameter of blade wmx Vessel diameter wmx Normal soft sphere interaction force at contact wNx Tangential soft sphere interaction force at contact wNx Gravity acceleration wmrs 2 x Height of plough wmx Total blade height wmx Blade pitch wmx Vessel height wmx Moment of inertia wkg m2 x Spring constant wNrmx Distance between particle and blade wmx Mass of particle wkgx Unit vector from the center of particle i to that of particle j Number of particles Radial coordinate wmx Time wsx Time step for calculation wsx Torque wN mx Particle velocity vector wmrsx Relative velocity wmrsx Position vector of the center of particle wmx Overlap distance wmx Damping coefficient Friction coefficient
ra rp v vb vp
Apparent density of particle wkgrm3 x Particle density wkgrm3 x Angular velocity vector of particle wradrsx Rotational speed of agitator wradrsx Angular velocity of particle in the core and blade regions wradrsx
References w1x E. Oshima, K. Yuge, Kagaku Kogaku 34 Ž1970. 779–785. w2x V.V. Chavan, J. Ulbrecht, Ind. Eng. Chem. Proc. Des. Dev. 12 Ž1973. 472–476. w3x R. Chowdhury, K.K. Tiwari, Ind. Eng. Chem. Proc. Des. Dev. 18 Ž1979. 227–231. w4x W.I. Patterson, P.J. Carreau, C.Y. Yap, AIChE J. 25 Ž1979. 508–516. w5x G. Schweier, Polym. Prepr. 38 Ž1997. 450–451. w6x H. Reichert, F. Vock, E. Kolk, R. Sinn, Ger. Chem. Eng. 1 Ž1978. 82–88. w7x B. Cooker, R.M. Nedderman, Powder Technol. 46 Ž1986. 263–269. w8x B. Cooker, F.R.G. Mitchell, R.M. Nedderman, Chem. Eng., May, 1983, 81–83. w9x B. Cooker, R.M. Nedderman, Powder Technol. 50 Ž1987. 1–13. w10x B. Cooker, R.M. Nedderman, Powder Technol. 52 Ž1987. 117–129. w11x P.A. Cundall, O.D.L. Strack, Geotechnique 29 Ž1979. 47–65. w12x T. Mikami, H. Kamiya, M. Horio, Chem. Eng. Sci. 53 Ž1998. 1927–1940. w13x Y. Kaneko, T. Shiojima, M. Horio, Chem. Eng. Sci. 54 Ž1999. 5809–5821. w14x T. Kawaguchi, T. Tanaka, Y. Tsuji, Powder Technol. 96 Ž1998. 129–138. w15x B.H. Xu, A.B. Yu, Chem. Eng. Sci. 52 Ž1997. 2785–2809. w16x P.A. Langston, U. Tuzun, D.M. Heyes, Chem. Eng. Sci. 51 Ž1996. 873–891. w17x S. Yuu, T. Abe, T. Saitoh, T. Umekage, Adv. Powder Technol. 6 Ž1995. 259–269. w18x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, Powder Technol. 93 Ž1997. 261–266. w19x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, J. Mech. Soc. Jpn., B 63r616 Ž1997. 3876–3882. w20x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, J. Mech. Soc. Jpn., B 64r619 Ž1998. 662–669. w21x C. Wightman, M. Moakher, F.J. Muzzio, O. Walton, AIChE J. 44 Ž1998. 1266–1276. w22x S. Sakamoto, J. Soc. Powder Technol. Jpn. 35 Ž1998. 508–513. w23x S. Matsusaka, M. Furutate, H. Masuda, J. Soc. Powder Technol. Jpn. 35 Ž1998. 168–173. w24x A. Shimosaka, S. Higashihara, J. Hidaka, J. Soc. Powder Technol. Jpn. 35 Ž1998. 242–249. w25x Y. Tsuji, T. Kawaguchi, T. Tanaka, Powder Technol. 77 Ž1993. 79–87. w26x H. Derrenbecher, F.R. Faulhaber, H.P. Kurz, K.H. Sartor, Ger. Chem. Eng. 7 Ž1984. 55–61.
Numerical analysis of particle mixing characteristics in a single helical ribbon agitator using DEM simulation Yasunobu Kaneko a
a,)
, Takeo Shiojima a , Masayuki Horio
b,1
Production Technology Center, Idemitsu Petrochemical Co. Ltd., Ichihara, 1-1 Anesaki-kaigan, Ichihara, Chiba 299-0193, Japan b Department of Chemical Engineering, Tokyo UniÕersity of Agriculture and Technology, Koganei, Tokyo 184-8588, Japan Received 24 May 1999; received in revised form 11 October 1999; accepted 11 October 1999
Abstract Numerical analysis of three-dimensional motion of particles in a single helical ribbon agitator was carried out by means of the Discrete Element Method ŽDEM.. To validate the computed results experiments were carried out with a cold scale model of 0.3 m inside diameter. Circulation time of particles in the agitator and the horizontal particle velocity distribution in the core region predicted by the present simulation agreed well with those obtained by experiments. Based on DEM simulation, the particle circulation and mixing characteristics in the agitator vessel were investigated in detail. Vertical mixing of particles was found rather poor during upward and downward flows through the blade and core regions, respectively. However, in the core region, particles flow in the manner of funnel flow. Additionally, in the blade region the vertical velocity decreases toward the wall and reduces to zero near the boundary between the blade and core regions. Particles were carried up to the bed surface by the helical ribbon flow down into the center of the core, and those not carried completely to the surface flow into the outer periphery of the core. The effect of bed height on circulation and mixing was analyzed numerically in terms of tracer particle concentration response as well as particle velocity and angular velocity distributions, which are quite difficult to obtain by experiment. The bed height was found to be the most important factor for mixing and circulation. When the bed height is appropriate, particles circulate rather regularly. If the bed height is lower than the blade top, particles mix rapidly because of the high fall of particles from the blade top. On the other hand, if the bed height is higher than the blade top, the velocity of particles in the bed surface region is much lower than that of the other region and accordingly a semi-stagnant eddy is formed. q 2000 Elsevier Science S.A. All rights reserved. Keywords: Discrete element method; Particle mixing; Helical ribbon agitator
1. Introduction Helical ribbon agitators have been widely used for mixing liquids of high viscosity. There are many reports providing basic information of helical ribbon agitators concerning power consumption of agitators, flow pattern andror mixing characteristics w1–4x. This type of agitators has also been applied to blend powders. However, there have been only few published reports on helical ribbon agitation of powders. Moreover, their literature focus has been only on the overall mixing characteristics and power consumption for agitation. ) Corresponding author. Tel.: q81-436-60-1850; fax: q81-436-601121; e-mail: [email protected] 1 E-mail: [email protected].
Fig. 1 illustrates a single helical ribbon agitator and a typical flow pattern of particles in the vessel. Particles are scraped by the inclined blade, and lifted up in the blade or upflow region by the ribbon. This is because particles cannot rotate at the same speed as the ribbon’s due to the friction between particles and the wall. At the top of the bed the particles fall down into the core or downflow region of the bed. Particles are thus circulated in the axial direction and mixed by the particle exchange between the core and the blade regions. Dead zone is scarcely formed and the power consumption is rather low. Accordingly, in the field of gas phase propylene polymerization reactors helical ribbon agitators have been applied successfully w5x. However, since polymerization is exothermic with a rather high heat of reaction and since the difference between the reaction temperature and melting temperature of
0032-5910r00r$ - see front matter q 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 3 2 - 5 9 1 0 Ž 9 9 . 0 0 2 5 1 - X
56
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
a rotary vessel mixer w18x and in a centrifugal tumbling granulator w19,20x. Wightman et al. w21x and Sakamoto w22x investigated the motion of particles in a rotating drum, and Matsusaka et al. w23x and Shimosaka et al. w24x applied DEM simulation to estimate particle behavior under vibrating conditions. However, the DEM simulation of particle motion in the vessel with an agitator as well as with the helical ribbon agitator has not been carried out yet. In the present study, an attempt was made to develop a three-dimensional DEM code to simulate particle motion in a helical ribbon agitator vessel. To evaluate the accuracy of the calculated results, experiments were carried out in a mixer cold test rig with the same geometry and size. Fig. 1. Schematic of a helical ribbon blade and flow pattern in a single ribbon agitator.
2. Theoretical model polypropylene is small, the temperature of polymer particles tends to rise and exceeds the melting temperature. If the heat of reaction generated in each particle is not effectively removed, the polymer particles tend to melt and form lumps, which may hinder the stable operation of the reactor. Recent progress in catalyst activity has been quite remarkable. On the other hand, the demand for random copolymerization with a rather high ethylene content has been increasing to lower the melting temperature for easier film production. Since these reduce the difference between reaction temperature and the melting temperature mentioned above, the reactor control for safe operation is becoming more important. Accordingly, more attention should be focused on the study of detailed particle mixing behavior. In the literature, Reichert et al. w6x studied the mechanical aspects, such as power consumption of the agitator, shape of the vortex and circulation time for a large-scale agitator design presenting a dimensionless correlation. Cooker and Nedderman w7x and Cooker et al. w8x studied the probability distribution of circulation times by means of a radio pill tracking system. They further proposed a theory to predict the power consumption for agitation and the circulation rate of particles based on stress analysis with appropriate modifications of Janssen’s method w9,10x. However, no reports can be found concerning the detailed particle mixing behavior in the helical ribbon agitator vessel. This is probably because of the difficulties in both experiment and mathematical analysis to pursue the motion of particles. The recent rapid growth of numerical power in the field of particle technology is now becoming quite promising to obtain valuable information on a variety of powder processes. Discrete element method ŽDEM. with soft sphere model initially proposed by Cundall and Strack w11x for quasi-static deformation of a particle bed has been successfully applied to analyze particle motions in densely packed beds such as fluidized beds w12–15x, andror hoppers w16,17x. In the field of mechanical particle agitation, Muguruma et al. applied DEM to simulate particle flows in
2.1. Assumption A mathematical model was developed based on the following assumptions for simplification: 1. Particles are spherical in shape and uniform in diameter and move in a three dimensional space. 2. The gas-to-particle interaction and the particle-to-particle adhesion are neglected. Forces considered are those for gravity and collisions between particles, between particle and wall and between particle and blade. 3. Particles are supposed to have a soft sphere interaction expressed by a Hookean linear spring, a dash pot and a friction slider with Coulomb’s law of friction. 4. Based on justifications in previous work w13,25x, the spring constant is adjusted only from the viewpoint of numerical economy. 2.2. Equation of particle motion Both translational and rotational motions are considered in equations of motion. The equations for individual particles can be written as follows: x¨ s Ž Fn H Ft . rm q g
Ž 1.
v ˙ s TcrI
Ž 2.
The above equations are the same as those of Tsuji et al. w25x, where x is the position vector of the center of the particle, Fn and Ft are the sum of the collision forces respectively for the normal and tangential directions acting on colliding particles, m is the mass of particles, g is the gravity acceleration, v is the angular velocity vector of a particle, Tc is the summation of torque caused by particle collision and I is the moment of inertia. The Euler scheme is used to integrate Eqs. Ž1. and Ž2. to obtain the new velocity, position and angular velocity of a particle after a time step. Fig. 2 illustrates the definitions of variables under the colliding condition of particles. When particle i is in
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
57
zr is the velocity vector of particle i relative to particle j. Subscripts ‘n’ and ‘t’ denote the normal and tangential directions and subscripts i and j denote the particle i and j, respectively. The direction of Fn is the one from point c to the center of particle j, and the direction of Ft is the same as the one of the grade crossing between the surface including the vector zr and the vector zn and the surface of solid crossing of both particles.
2.3. Description of the blade and the judgment of particleto-blade collisions The following equations with a parameter u are used to express the spiral curve of a blade shown in Fig. 3 with radius R and pitch Hp :
Fig. 2. Schematic illustration of particle collision.
contact with particle j at point x c , the normal and tangential components of the collision force acting on particle i are given respectively by: Fn s yk dn y h zn
Ž 3. Ž 4a .
Ft s yk dt y h zt < Ft < F m < Fn < zt < Ft < ) m < Fn < Ft s ym < Fn < < zt <
Ž 4b .
where zn s Ž zr P n . n zt s zr = n zr s zj y zi q Ž v i q v j . = Ž x i y x c .
Ž 5. Ž 6. Ž 7.
where k is the spring constant, d is the particle displacement, h is the damping coefficient, z is the particle velocity vector, m is the friction coefficient, n is the unit vector from the center of particle i to that of particle j and
x s Rcos u y s Rsin u z s bu Ž b ' Hpr2p .
Ž 8.
A helical ribbon blade can be represented by changing R from the inner radius of the blade to its outer radius. Point Ž x b , y b , z b . on the blade, where there is a possibility of contact with particle i having its center at Ž x i , yi , z i ., is expressed by the following equations: x b s ri cos u y b s ri sin u zb s bŽ u y v b t . ri ' Ž x i2 q yi2 .
Ž 9.
0.5
where v b is the angular velocity of the agitator rotation.
Fig. 3. Mathematical expression of the single helical agitator.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
58
The distance l between particle i Ž x i , yi , z i . and a point on the blade Ž x b , y b , z b . is obtained as: 2
2
l 2 s Ž x i y x b . q Ž yi y y b . q Ž z i y z b . 2
s Ž x i y ri cos u . q Ž yi y ri sin u . qŽ zi y b Ž u y v b t . .
2
2
2
Ž 10 .
Differentiating Eq. Ž10. by u , we obtain: dŽ l 2 . du
s 2 ri Ž x i sin u y yi cos u . y 2 b Ž z i y b Ž u y v b t . .
Ž 11 . Setting dŽ l 2 .rd u s 0 and solving Eq. Ž11. for u , we obtain l min , the minimum distance between particle i and the blade, from Eq. Ž9.. In the present computation the Newton–Raphson method was used for iteration. If l min is smaller than particle radius, it is judged that the particle i and the blade are in collision.
Table 1 Geometry of the agitator vessel and property of particles Vessel Height Inside diameter Blade Total height Inside diameter Outside diameter Pitch Height of plough Particle Mean diameter Apparent density
Hv Dv
0.7 m 0.31 m
Hb Di Do Hp hp
0.495 m 0.18 m 0.28 m 0.11 m 0.08 m
dp ra
1.0 mm 550 kgrm3
distribution. In all experiments the level of the powder bed was adjusted to the top edge of the agitator blade for each mixing condition. 3.2. Circulation time
3. Experimental 3.1. Apparatus The single helical ribbon agitator vessel used for the present experiment is illustrated in Fig. 4. Table 1 shows the geometry of the agitator vessel and the properties of particles. The motor and converter assembly was able to successively vary the rotation speed of the agitator up to 150 rpm. The powder was polypropylene homopolymer having approximately 1 mm in diameter with a narrow size
One of the most important factors in estimating the performance of an agitator is the average circulation time. To obtain the circulation time the following two methods were adopted: Ž1. Particles are agitated with a colored rubber ball having 20 mm in diameter. If the ball reached near the bed surface, the ball was able to be distinguished from the polypropylene particles because they were white. The time required for the ball to emerge again near the bed surface was measured more than fifty times and the particle circulation time was determined. Ž2. One kilogram of polypropylene pellets of about 3 mm in diameter, i.e., the weight ratio of the pellets to the whole bed weight was about 0.05, were poured into the center of the bed surface as tracer particles. A small amount of the bed material was collected from the center of the bed surface every 2 s for about 1 min. The weight concentration of the pellets in the sampled particles was obtained by sieving. From the oscillatory plot of pellet concentration against elapsed time, the time interval between the neighboring peaks was determined as the particle circulation time. Simulation was done for conditions similar to method Ž2. of the above two experiments. 3.3. Horizontal distribution of Õertical Õelocity of particles in the core region
Fig. 4. Experimental apparatus.
Fig. 5 shows the set-up for measurement of vertical velocity of particles in the core region. A thin glass tube of 10 mm inside diameter was immersed in the bed. A funnel was installed at the top of the tube to easily pour particles into the glass tube. The tube was perpendicularly fixed by the tube guide which was fastened on the base plate with bolts so that it can be horizontally moved. Polypropylene particles, the same as the bed material, were filled up in the tube before the start of agitation.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
59
8 mm. For this computation condition the CPU time required for computation of real 1 s was approximately 8 h. The time step, an important factor for the stable calculation, was determined by the following equation proposed by Tsuji et al. w25x: DtF
2p 5
(
m
Ž 12 .
k
In this study the spring constant, k, was assumed to be 800 Nrm, so that the time step for integration D t was as large as 1.0 = 10y4 s.
5. Results and discussion 5.1. Comparison between the results of experiment and calculation
Fig. 5. Set-up to determine the vertical velocity of particles in the core region.
Right after the start, particles in the glass tube also start to move down at a constant rate, which is supposed to correspond to a vertical flow rate of particles at a horizontal position of the bed Žat the tube exit.. From the descending rate of the upper surface of particles in the straight part of the tube, the vertical rate of particles was determined. The same procedure was repeated over different horizontal positions, but at a fixed height Ž0.4 m from the bottom. in the bed. Another tube of 20 mm I.D. was also applied but no appreciable difference was observed.
4. Computation conditions Computation conditions and parameters used in the present study are shown in Table 2. The column size and the agitator dimensions were identical except for geometry of the vessel bottom, i.e., in the computation a flat bottom was assumed for simplification. Computation was performed with a personal computer ŽFujitsu, FMV-6266T6.. The particle number in the present experiment is roughly 4.0 = 10 7, for which the numerical simulation with the present DEM model cannot be practically carried out. To reduce the computational load the number of particles was reduced to 1.0 = 10 5 by increasing particle diameter up to
Fig. 6 shows the typical snapshot of the mixing behavior under the condition of 3.0p radrs of the agitation speed. The black and gray particles indicate having upward and downward velocity vectors, respectively. As can be clearly seen in the numerical results of Fig. 6, particles in a helical ribbon agitator circulate by being lifted up along the wall by the rotating blade and descending in the core. The revolution of powder sinking at the center causes a concave bed surface. The vortex depth, i.e., the difference of height between the bed surface at the center and that at the wall, generally increased with the rotating speed. Fig. 7 shows the snapshots of mixing behavior at vertical cross section for three rotational speeds. From this figure, it is shown that the maximum vortex depth gradually increases with the rotating speed. This result is in agreement with the general observation of the vortex profile of helical ribbon agitators. To validate the computed particle circulation times they were compared with the experimental values. Computed circulation time was obtained by the following procedure; Ži. After reaching a stable mixing, particles located within 0.1 m from the top surface of the core region shown in Fig. 8 were defined as tracer particles. The above region was also defined as the inspection volume. Žii. The concentration of tracer particles in the inspection volume was
Table 2 Computation conditions Number of particles Particle diameter Spring constant Damping coefficient Friction coefficient Particle density Time step Rotational speed
np dp k h m rp Dt vb
7.0, 9.0, 10.0=10 4 8 mm 800 Nrm 0.6 0.3 800 kgrm3 1.0=10y4 s 2.0p , 3.0p , 4.0p radrs
60
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
Fig. 6. Typical snapshot of computed particle behavior at t s 9.9 s.
calculated as a function of time. Žiii. There are several peaks in a concentration vs. time curve. The time interval between the adjacent peaks is defined as the particle circulation time as shown in the lower half of Fig. 8. Fig. 9 shows that the particle circulation time as a function of rotation speed. Experimental values were obtained by the one ball tracer test and particle tracer test. As can be clearly seen in Fig. 9, the computed circulation time agreed well with that of the experiments.
The responses of tracer particle concentration for both experiment and calculation are compared in Fig. 10. The period of the initial low concentration corresponds to the time of first circulation. The time at which the tracer particles enter again into the inspection volume is approximately same for both cases. In Fig. 11, the computed horizontal distribution of particle velocity Žvertical component. is compared with that from the experiment. As defined in Fig. 11, the
Fig. 7. Effect of rotational speed on particle behavior.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
61
Fig. 8. Method of determining the circulation time from calculation.
inspection volume, i.e., the region between 0.4 and 0.45 m from the bottom, was divided coaxially into 15 sections Žnine for the core and six for the blade region.. The computed values were obtained by averaging the vertical velocities of particles in each section. In the previous investigations w6,9,26x, the flow pattern of the core region was treated as the plug flow. However, the present result both from experiment and computation indicates that the flow in the core is a sort of funnel flow, having a nonuniform velocity distribution. In addition, in the blade region the vertical velocity decreases toward the wall because a finite clearance existed between the blade tip and the wall. On the other hand, getting closer to the boundary between the blade and the core regions, the vertical velocity component reduces to zero due to the momentum exchange between the two regions.
5.2. Detail particle flow and the effect of bed height on particle motion
Fig. 9. Model validation in terms of circulation time.
Fig. 10. Tracer concentration vs. time.
Fig. 12 shows the snapshots of particle locations, where particles are painted in different gray levels depending on their initial locations. From this figure it can be seen that vertical mixing of particles is rather poor during upward and downward motion through the blade and core regions, respectively. Particle mixing takes place near the bed
62
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
Fig. 11. Vertical velocity distribution.
surface by the particle flow into the cavity of the core top. Moreover, from Fig. 12Žd. the particles located at the middle height on the blade at t s 0 s flow through the outer side of the core region. Particles that were lifted not completely up to the bed surface tend to descend through the outside section of the core. Fig. 13 shows the snapshots of particle motion for three different bed heights. The shape of vortex is quite different for different bed heights. The effect of bed height on particle mixing seems considerable even though the rotational speed is the same. To make this point clearer, the mixing process is examined by tracking the tracer particles at three different regions, i.e., the top region of the bed, the
middle region of the core section and the bed bottom region. Fig. 14 shows the concentration of tracer particles in three regions Žor inspection volumes. for three cases of bed heights. Regions 2 and 3 have the same heights for all three cases, but the definition of region 1 is different for each case so that the tracer concentration in the surface region can be examined. The initial position of tracer particles, inspection volumes and particle velocity vectors in x–z plane are also depicted in the upper half of Fig. 14. In case Žb., a rather sharp peak can be found in all three regions and intervals between adjacent peaks are the same for all regions indicating the solid circulation is quite through from the bottom to the top. In case Ža. with a low bed height, solid circulation only took place between regions 2 and 3. In the top region the tracer particle concentration did not change after some of them were carried up from region 2. Furthermore, tracer concentration in all
Fig. 12. Mixing process of particles.
Fig. 13. Effect of bed height on particle mixing.
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
63
Fig. 14. Effect of bed height on particle mixing process.
regions approached to an equilibrium value rather quickly. The higher vertical distance between the blade top and the core surface should have helped the tracer mixing. When the bed surface is much higher than the blade top Žcase Žc.., the particle velocity near the bed surface region is rather low as shown in the upper half of Fig. 14. This indicates that the near surface region, where no forced circulation exists, acts as a semi-stagnant eddy region. The formation of this semi-stagnant eddy region helps particles
mixing by exchange but disturbs the clear circulation. Fig. 15 shows the angular velocity of particles both in the core and blade regions for cases Žb. and Žc. of Fig. 14. By experiment the angular velocity in the blade region usually cannot be easily obtained. This indicates a clear advantage of DEM simulation. Obviously, the near surface region of case Žc. rotates much more slowly than the lower section. Furthermore, the rotation up to the height of blade top is almost the same for different bed heights. Bed height is clearly one of the most important operating parameters for the helical ribbon agitators.
6. Conclusions
Fig. 15. Vertical profiles of angular velocity Ž v b s 4p wradrsx..
Three-dimensional motion of particles mixed in a single helical ribbon agitator vessel was numerically analyzed by the DEM simulation. To reduce the computation load, a spring constant smaller than the reality and larger particles were adopted for computation. Experiments were conducted and computed results were successfully compared with experimental data. By tracking particle motion in detail, the circulation and mixing characteristics were analyzed in terms of tracer particle concentration response, particle trajectories and
Y. Kaneko et al.r Powder Technology 108 (2000) 55–64
64
angular velocity distribution, some of which are quite difficult to obtain by experiment. It was demonstrated that the distance between the blade top and the bed surface is the most important parameter to attain good mixing or good circulation. When bed surface is lower than the blade top, particles lifted to the blade top fall into the center of the core cavity, enhancing rapid mixing. When bed surface is higher than the blade top, a semi-stagnant eddy was formed on the bed surface. If the bed height is appropriate, a proper circulation throughout the bed height is created.
7. List of symbols dp Dt Do Dv Fn Ft g hp Hb Hp Hv I k l m n np r t Dt Tc z zr x d h m
Particle diameter wmx Inside diameter of blade wmx Outside diameter of blade wmx Vessel diameter wmx Normal soft sphere interaction force at contact wNx Tangential soft sphere interaction force at contact wNx Gravity acceleration wmrs 2 x Height of plough wmx Total blade height wmx Blade pitch wmx Vessel height wmx Moment of inertia wkg m2 x Spring constant wNrmx Distance between particle and blade wmx Mass of particle wkgx Unit vector from the center of particle i to that of particle j Number of particles Radial coordinate wmx Time wsx Time step for calculation wsx Torque wN mx Particle velocity vector wmrsx Relative velocity wmrsx Position vector of the center of particle wmx Overlap distance wmx Damping coefficient Friction coefficient
ra rp v vb vp
Apparent density of particle wkgrm3 x Particle density wkgrm3 x Angular velocity vector of particle wradrsx Rotational speed of agitator wradrsx Angular velocity of particle in the core and blade regions wradrsx
References w1x E. Oshima, K. Yuge, Kagaku Kogaku 34 Ž1970. 779–785. w2x V.V. Chavan, J. Ulbrecht, Ind. Eng. Chem. Proc. Des. Dev. 12 Ž1973. 472–476. w3x R. Chowdhury, K.K. Tiwari, Ind. Eng. Chem. Proc. Des. Dev. 18 Ž1979. 227–231. w4x W.I. Patterson, P.J. Carreau, C.Y. Yap, AIChE J. 25 Ž1979. 508–516. w5x G. Schweier, Polym. Prepr. 38 Ž1997. 450–451. w6x H. Reichert, F. Vock, E. Kolk, R. Sinn, Ger. Chem. Eng. 1 Ž1978. 82–88. w7x B. Cooker, R.M. Nedderman, Powder Technol. 46 Ž1986. 263–269. w8x B. Cooker, F.R.G. Mitchell, R.M. Nedderman, Chem. Eng., May, 1983, 81–83. w9x B. Cooker, R.M. Nedderman, Powder Technol. 50 Ž1987. 1–13. w10x B. Cooker, R.M. Nedderman, Powder Technol. 52 Ž1987. 117–129. w11x P.A. Cundall, O.D.L. Strack, Geotechnique 29 Ž1979. 47–65. w12x T. Mikami, H. Kamiya, M. Horio, Chem. Eng. Sci. 53 Ž1998. 1927–1940. w13x Y. Kaneko, T. Shiojima, M. Horio, Chem. Eng. Sci. 54 Ž1999. 5809–5821. w14x T. Kawaguchi, T. Tanaka, Y. Tsuji, Powder Technol. 96 Ž1998. 129–138. w15x B.H. Xu, A.B. Yu, Chem. Eng. Sci. 52 Ž1997. 2785–2809. w16x P.A. Langston, U. Tuzun, D.M. Heyes, Chem. Eng. Sci. 51 Ž1996. 873–891. w17x S. Yuu, T. Abe, T. Saitoh, T. Umekage, Adv. Powder Technol. 6 Ž1995. 259–269. w18x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, Powder Technol. 93 Ž1997. 261–266. w19x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, J. Mech. Soc. Jpn., B 63r616 Ž1997. 3876–3882. w20x Y. Muguruma, T. Tanaka, S. Kawatake, Y. Tsuji, J. Mech. Soc. Jpn., B 64r619 Ž1998. 662–669. w21x C. Wightman, M. Moakher, F.J. Muzzio, O. Walton, AIChE J. 44 Ž1998. 1266–1276. w22x S. Sakamoto, J. Soc. Powder Technol. Jpn. 35 Ž1998. 508–513. w23x S. Matsusaka, M. Furutate, H. Masuda, J. Soc. Powder Technol. Jpn. 35 Ž1998. 168–173. w24x A. Shimosaka, S. Higashihara, J. Hidaka, J. Soc. Powder Technol. Jpn. 35 Ž1998. 242–249. w25x Y. Tsuji, T. Kawaguchi, T. Tanaka, Powder Technol. 77 Ž1993. 79–87. w26x H. Derrenbecher, F.R. Faulhaber, H.P. Kurz, K.H. Sartor, Ger. Chem. Eng. 7 Ž1984. 55–61.