- Email: [email protected]
Longitudinal and transverse mixing in rotary kilns: A discrete element method approach
Chemical Engineering Science 60 (2005) 4083 – 4091 www.elsevier.com/locate/ces
Longitudinal and transverse mixing in rotary kilns: A discrete element method approach G.J. Finniea , N.P. Kruytb,∗ , M. Yec , C. Zeilstrac , J.A.M. Kuipersc a Department of Chemical Engineering, Soshanguve Campus, Tshwane University of Technology, Private Bag X07, Pretoria North, Republic of South Africa b Department of Mechanical Engineering, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands c Department of Chemical Engineering, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands
Received 13 September 2004; received in revised form 20 October 2004; accepted 4 December 2004 Available online 7 April 2005
Abstract Longitudinal and transverse mixing in horizontal rotary kilns has been studied, using a three-dimensional Discrete Element Method approach. The focus is on the effect of the main operating conditions, i.e., the filling degree and the rotational speed of the rotary kiln, on quantitative measures of longitudinal and transverse mixing. Discrete Element Method simulations have been performed for various values of the main operating conditions. Flow regimes have been determined from visualisations of the simulations. It is shown that the mixing in the longitudinal direction can be described by a diffusion equation. The dependence of the diffusion coefficient on the operating conditions has been determined from the results of the simulations. It is found that the diffusion coefficient increases linearly with rotational speed, while the influence of the filling degree is relatively small. The mixing in the transverse plane is quantified by an entropy-like mixing index. The results of the simulations show that this mixing index varies exponentially with the number of revolutions of the rotary kiln. The dependence on the operating conditions of the transverse mixing speed, as determined from the exponential behaviour of the mixing index, has been determined. This mixing speed decreases with increasing rotational speed and with increasing filling degree. The transverse mixing speed that has been determined from the results of additional two-dimensional simulations is generally larger than that observed in the corresponding three-dimensional simulations. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Rotary kiln; Granular material; Mixing; Discrete Element Method; Diffusion; Simulation
1. Introduction Rotary kilns are widely used in chemical and mechanical engineering for processing of granular materials for various purposes, such as for mixing, drying and coating of materials, for calcining of limestone, for reduction of iron ore and for waste incineration. Rotary kilns are able to handle varying feedstocks, in particular with
∗ Corresponding author. Tel.: +31 53 489 2528; fax: +31 53 489 3695.
E-mail addresses: fi[email protected] (G.J. Finnie), [email protected] (N.P. Kruyt), [email protected] (M. Ye), [email protected] (C. Zeilstra), [email protected] (J.A.M. Kuipers). 0009-2509/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.12.048
wide particle-size distributions and significant differences in physical properties. The operating conditions of a rotary kiln (as sketched in Fig. 1) involve rotational speed around its axis, radius R, length L in longitudinal direction, inclination with respect to the horizontal direction and filling degree J (i.e., the fraction of the volume of the rotary kiln that is filled with granular material, including voids between particles), as well as parameters related to the feedstock. In the sequel the focus is on horizontal rotary kilns, i.e., = 0◦ . Inclined rotary kilns were studied by Sudah et al. (2002). Generally, the main operating conditions of rotary kilns are expressed in terms of non-dimensional numbers by the filling degree J and the Froude number Fr, which is
4084
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
nificant segregation is expected. A recent overview of mixing and segregation of granular materials is given by Ottino and Khakhar (2000). A computational DEM study of mixing of binary particle mixtures (including density differences between the two sizes) is reported by Dury and Ristow (1999). At present, gas flows in the rotary kilns are not considered here. Separate descriptions are formulated for the longitudinal and transverse mixing. The results of the DEM simulations are analysed to determine quantitatively the influence of nondimensional operating conditions of the rotary kiln on the parameters describing the longitudinal and transverse mixing. The outline of this study is as follows. Firstly, the DEM simulations are described in Section 2. This is followed in Section 3 by the formulation of measures to quantify longitudinal and transverse mixing. The results of the simulations are analysed in Section 4. Finally, findings from this study are discussed.
Fig. 1. Sketch of the geometry of the horizontal rotary kiln.
defined by
R , g
2. Discrete element method simulations
2
Fr =
(1)
where g is the acceleration due to gravity. Two different types of mixing can be distinguished in rotary kilns. The first type occurs in the plane transverse to the axis of rotation, while the second type occurs in longitudinal direction. The first type of mixing is vigorous, while the second is much slower. Depending on the non-dimensional operating conditions, several flow regimes may be observed in the plane transverse to the axis of rotation, namely slipping, slumping, rolling, cascading, cataracting and centrifuging (Henein et al., 1983a,b; Mellmann, 2001). In industrial operations the rolling regime is mostly encountered, since it promotes good mixing of particles. This flow regime is characterised by two distinct regions: (i) a relatively thin active layer is formed at the top in which the material flows along the inclined slope and (ii) a thick passive region in which the material rotates as a rigid body with the angular velocity of the rotary kiln (see for example, Boateng and Barr, 1997). Scaling relationships for rotary kilns in the rolling and slumping regime have been studied by Ding et al. (2001). The present study is a basic computational study of longitudinal and transverse mixing in horizontal rotary kilns, focusing on the effect of the main operating conditions on quantitative measures of longitudinal and transverse mixing. A recent overview of computational methods for studying mixing of granular materials is given by McCarthy et al. (2000). The method used here to study the longitudinal and transverse mixing processes is that of three-dimensional Discrete Element Method (DEM) simulations, as proposed by Cundall and Strack (1979a,b), for a number of nondimensional operating conditions. The particle-size distribution employed in the simulations is fairly narrow, so no sig-
The computational method employed here is the threedimensional DEM (also called Discrete Particle Method), as proposed by Cundall and Strack (1979a,b). In this method the motion of a large number of particles is determined numerically by integrating the equations of motion, using a suitable model, the so-called contact constitutive relation, for the particle–particle interactions and the particle–wall interactions. Here, equations of motions are solved for particle translations as well as for particle rotations, contrary to what was done by Schutyser et al. (2001). Particle shapes are restricted to spheres (with varying radii) for reasons of computational efficiency. The contact constitutive relation employed here is basically that of Cundall and Strack (1979a,b). It involves (linear) elastic behaviour, Coulomb frictional behaviour and viscous contact damping to model dissipative collisions between particles. A linear spring is acting in the direction normal to the contact and a tangential spring is acting in the direction of the tangential component of the relative displacement at the contact. The respective spring constants are denoted by kn and kt . Coulomb friction is modelled by an interparticle friction coefficient , such that at each contact c we have ftc |fnc | where fnc is the (scalar) normal component of the contact force vector f c , i.e., fnc = f c · n c with n c the unit normal vector at the contact c, ftc is the magnitude of the tangential component of the force, i.e., ftc = f c − fnc n c and |fnc | is the absolute value of the normal component of the contact force. Viscous contact damping is implemented in terms of a coefficient of restitution cr (see for example, Hoomans, 2000). A similar contact constitutive relation is used for the interactions between particles and cylinder walls and end walls of the rotary kiln, but with a different value for the friction coefficient . Unlike in the simulations by
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
Yang et al. (2003), no rolling resistance was included here in the contact constitutive relation, considering the problem of establishing a manner for determining the material parameters involved in the rolling resistance model. The radii r of the spherical particles follow a (truncated) lognormal particle-size distribution with rmax /rmin = 1.4. The ratio of the standard deviation over the mean of the particle-size distribution equals 0.04. Hence, the particlesize distribution is fairly narrow, so no significant size segregation is expected. Furthermore, all particles have the same density, so no density segregation will occur. The elastic parameter kn has been selected such that the “overlap” between particles is small, yet the time step in the (explicit) numerical integration of the equations of motion does not become exceedingly small. The stiffness ratio kt /kn = 0.29. A fairly high interparticle friction coefficient is 0.5 is selected, since real particles are not perfectly spherical (as is assumed here). Non-spherical particles have a higher effective (continuum) friction, in comparison to spherical particles with the same interparticle friction coefficient, due to suppression of rolling (Rothenburg and Bathurst, 1991). A high friction coefficient of 1.5 is used for the interactions between particles and cylinder wall, corresponding to a roughened wall. The friction coefficient of the end walls of the rotary kiln is 0.01. The selected restitution coefficient cr is 0.71. Note that the value of the (microscopic) interparticle friction coefficient and the so-called (macroscopic) internal friction coefficient of the granular material as a continuum (related to the static angle of repose of the granular material) are generally not equal. In fact, there is no direct relation between these two friction coefficients (Skinner, 1969). Similarly, the value of the particle–wall friction coefficient is generally not equal to the wall friction coefficient of the granular material as a continuum. The geometry of the rotary kiln (see also Fig. 1) is such that its aspect ratio, L/R, equals 2.0. The ratio of the kiln radius R over the average particle radius r¯ is R/¯r = 40.0. Several assemblies were simulated, with the number of particles in the simulations, Np , depending on the filling degree J of the kiln: for J = 20%, Np = 10, 667, for J = 30%, Np = 16, 000 and for J = 40%, Np = 21, 333. These numbers of particles were determined from considerations of acceptable computer simulation times. The time step in the explicit numerical integration of the equations of motion is such that it is a factor 30 smaller than the maximum time step allowed for stability of the numerical method. Five values for the rotational speed of the kiln, , have been considered. These correspond to Froude numbers, see (1), of Fr = 0.25, 0.56, 1.00, 1.56, 2.25. In total, 15 simulations have been performed for the five values of Froude number Fr, and for each Froude number with the three filling degrees J = 20%, 30%, 40%. For each simulation the simulated time corresponds to 40 revolutions of the rotary kiln.
4085
3. Description of longitudinal and transverse mixing Two different types of mixing occur in rotary kilns. The first type is mixing in the transverse plane. This type of mixing is vigorous and quickly leads to mixed-out states. The second, much slower type is mixing in the longitudinal direction. Separate descriptions are formulated for these two types of mixing. 3.1. Longitudinal mixing In horizontal rotary kilns there is no net transport in the longitudinal direction. Motion of the particles in longitudinal direction is the result of interparticle collisions induced by the rotation of the rotary kiln. This motion is therefore expected to be completely random. Consequently, the mixing process is considered as analogous to diffusion, and hence Fick’s law should apply. This leads to the one-dimensional diffusion equation (Hogg et al., 1966) j j jc(x, t) c(x, t) = D∗ jt jx jx
(2)
for practically identical particles, where c(x, t) is the relative concentration of the mixture after time t at position x in longitudinal direction and D ∗ (with units m2 s−1 ) is the diffusion coefficient that is assumed to be spatially uniform in the sequel. By performing the transformation n = (t)/(2), this equation can be expressed in terms of the (non-dimensional) number of revolutions, n, of the rotary kiln instead of time t. The resulting equation becomes j j jc(x, n) c(x, n) = D , jn jx jx
(3)
where the units of D are m2 . The diffusion coefficients D ∗ and D are related by D=
2 ∗ D .
(4)
With initial conditions c(x, 0) = 0 for −L/2 < x < 0 and c(x, 0)=1 for 0 x < L/2 (corresponding to “blue” particles on the left and “red” particles to the right of the interface at x = 0; see also the sketch of the initial configuration for the longitudinal mixing on the left-hand side of Fig. 2) and boundary conditions jc(x, n)/jx|x=±L/2 =0 (corresponding to zero flux at the boundaries), the solution of the diffusion equation (3) is (Hogg et al., 1966) ∞ 1 2 1 −(2j − 1)2 2 Dn c(x, n) = + exp 2 2j − 1 L2 j =1
(2j − 1)x × sin . L
(5)
4086
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 Longitudinal
“blue”
3.2. Transverse mixing
Transverse
“red”
“gray”
“black”
Fig. 2. Sketch of initial configurations for longitudinal mixing (left) and transverse mixing (right).
The (time-dependent) non-dimensional mean position of the “red” particles is defined by L/2 1 xc(x, n) dx. (6)
c(n) = 2 L −L/2 Using (5), it follows (Kohring, 1995) that c(n) is given by ∞ 8 (−1)j +1 2 2 Dn
c(n) = 3 exp −(2j − 1) 2 . (7) L (2j − 1)3 j =1
After some algebra it follows that for small values of the Fourier number (Dn)/L2 , this expression can be approximated by 1 Dn
c(n) − 2 2 4 L
for
Dn 1. L2
(8)
Thus, for small Fourier number, the relationship between mean position of the “red” particles and (non-dimensional) time is linear. It is expected that the diffusion coefficient D depends on the process parameters and J , the gravitational acceleration g, as well as on the geometry of the rotary kiln, i.e., its radius R and length L. Furthermore, it assumed that of the material properties only the average particle size r¯ , density and frictional material parameters (collectively denoted by here) are important. This means that the influence of elastic and damping parameters of the granular material is excluded. Thus we obtain the following functional relationship: D = F (, J, g, R, L, r¯ , , ).
(9)
From dimensional analysis we find the functional relationships for D and D ∗ D R L , = G Fr, J, , , R2 r¯ R D∗ R L ∗ , . (10) Fr, J, = G , R 2 r¯ R An analogous dimensional analysis for the motion in the transverse plane was performed by Ding et al. (2001).
In analogy to the case of longitudinal mixing, as described in Section 3.1, Heydenrych (2001) solved a simplified convection-diffusion equation for the mixing in the transverse plane. This complex approach was not pursued here, considering the spatial variations in velocity that have to be taken into account in the convective part of the convectiondiffusion equation for the concentration c(x, n) and the possible spatial variations of the diffusion coefficient D in the transverse plane. Instead, transverse mixing is quantified here using an “entropy”-like quantity, as proposed by Schutyser et al. (2001, 2002). In this approach a grid with dimensions (I, J ) is defined that overlaps the region of interest. For each grid cell (i, j ), an “entropy” Si,j is calculated as Si,j = −kV i,j [pi,j ln pi,j + qi,j ln qi,j ],
(11)
where k is a constant, Vi,j is the volume occupied by the particles in cell (i, j ), pi,j is the concentration (fraction) of “black” particles, while qi,j = 1 − pi,j is the concentration (fraction) of “grey” particles. By definition, pi,j is the ratio of the total volume in cell (i, j ) occupied by “black” particles over the total volume in cell (i, j ) occupied by “black” and “grey” particles. Hence, the voids between particles do not affect the concentrations pi,j and qi,j . Note that many particles, with differing longitudinal position, may be mapped into the same cell in the transverse plane. Initially, the colours of the particles are set to “grey” when they are located in the left-hand side (when viewed along the kiln axis) of the rotary kiln, and to “black” when they are located in the right-hand side; see also the sketch of the initial configuration for the transverse mixing on the right-hand side of Fig. 2. Note that this “entropy” Si,j is not related to the thermodynamical entropy of the granular material. Such information-theoretical entropies are used in many applications (Katz, 1967; Kapur and Kesavan, 1992). A “mixing index” M(n) is now defined by M(n) =
Si,j (n) i,j
S∞
,
(12)
where the “entropy” of the perfectly mixed state, S∞ , is given by S∞ = −kV [p0 ln p0 + q0 ln q0 ], V is the total volume occupied by the particles, p0 is the concentration (fraction) of “black” particles in the system and q0 =1−p0 is the concentration (fraction) of “grey” particles in the system. In the initial state the mixing index M(0) = 0, while in the perfectly mixed state that is attained as n → ∞, the mixing index equals M(n → ∞)=1. Note that the constant k in Eq. (11) cancels in expression (12) for the mixing index M(n). In the determination of the concentrations pi,j , it was taken into account that particles may occupy multiple cells. Thus, the fractional volumes of particles inside the cells were used to determine pi,j and qi,j , instead of basing their determination solely on the coordinates of particle centres.
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 Kiln speed, Ω⋅R
Fig. 3. Employed fractional volumes in the determination of the cell concentrations pi,j and qi,j .
This procedure is shown schematically in Fig. 3 for disks. The (approximate) computation of the fractional volumes is straightforward (Hoomans, 2000), provided that the length of the cells is larger than the maximum particle diameter. This procedure reduces the dependence of the computed “entropy” on the grid sizes (especially when the number of particles in each cell is fairly small, such as in twodimensional simulations) that was reported by Schutyser et al. (2001).
4. Results The results of the three-dimensional DEM simulations have been analysed with respect to occurring flow regimes, longitudinal mixing and transverse mixing. Additional twodimensional DEM simulations have been performed in order to compare the transverse mixing in two-dimensional and three-dimensional simulations. 4.1. Flow regimes Depending on the non-dimensional process parameters, such as Froude number Fr and filling degree J, various flow regimes may be distinguished: sliding, surging, slumping, rolling, cascading, cataracting and centrifuging. An overview of the flow regimes, together with transition criteria between the flow regimes, is given by Mellmann (2001).
4087 Kiln speed, Ω⋅R
Fig. 4. Velocity of particles in a transverse cross-section. Left: rolling regime for J = 40% and Fr = 0.25; right: cataracting regime for J = 40% and Fr = 2.25. Note that the magnitude of the velocity vectors has been scaled for visual clarity. For comparison, the kiln speed, R, has also been plotted at the top of the figures.
Based on analyses of the particle-velocity profiles and on visualisations of the results of the simulations, the flow regimes were identified, see Table 1. Examples of the observed velocity profiles are shown in Fig. 4 for the rolling regime (left) and for the cataracting regime (right). A direct quantitative comparison of these results with the transition criteria given by Mellmann (2001) is not possible, since these criteria involve the continuum internal friction and wall friction coefficients. As noted in Section 2, these continuum friction coefficients are not (necessarily) equal to the particle–particle and the particle–wall friction coefficients, respectively. Note that the transition from the cataracting regime to the centrifuging regime does not have to correspond to Fr=1, as is discussed by Mellmann (2001). The centrifuging regime has not been clearly observed in the simulations. This is attributed to the low effective friction between the granular material and the wall of the rotary kiln that is caused by the rolling tendency of spheres and disks (see for example, Rothenburg and Bathurst, 1991). Shinbrot et al. (1999) experimentally observed chaotic transverse mixing for very fine particles. Since the particles considered here are fairly large, this chaotic mixing was not observed here. 4.2. Longitudinal mixing In the initial configuration at n=0, all particles on the left were “blue” and all particles to the right of the interface at
Table 1 Flow regimes in three-dimensional simulations
J = 20% J = 30% J = 40%
Fr = 0.25
Fr = 0.56
Fr = 1.00
Fr = 1.56
Fr = 2.25
Rolling Rolling Rolling
Rolling Rolling Rolling
Rolling/cataracting Rolling/cataracting Rolling/cataracting
Cataracting Cataracting Cataracting
Cataracting Cataracting Cataracting
4088
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
0.25
1 0.9
0.245
0.8 0.7
0.24
c(x,n)
(n)
0.6 0.235
0.5 0.4
0.23
0.3 0.2
0.225
n=10 revs n=20 revs n=30 revs n=40 revs
0.1 0.22 0
5
10
15
20
25
30
35
0 -0.5
40
n
x =0 were “red” (see also the left-hand side of Fig. 2). Since the evolution of all particle positions is determined in the DEM simulations, concentration profiles can be determined by dividing the longitudinal extent of the rotary kiln into slices and computing, for each time instant, the volume fraction of “red” particles in each slice. Thus, time-dependent (with time expressed in terms of the number of revolutions, n) concentrations profiles c(x, n) were obtained, as well as time-dependent mean positions c(n) of the “red” particles. The diffusion coefficients D were determined from the simulations by considering the evolution of the mean position c(n) with number of revolutions n. According to Eq. (8), c(n) depends linearly on Fourier number (Dn)/L2 for small (Dn)/L2 . The results of the simulations confirm this linear relationship between c(n) and n, see Fig. 5. According to the theory presented in Section 3.1 for the longitudinal mixing, the concentration profiles are given by Eq. (5). A typical comparison of the actual and the theoretical concentration profiles is given in Fig. 6. This Figure shows that the actual concentration profiles, as determined from the simulations, closely correspond to the theoretical concentration profiles that follow from the solution of the diffusion equation (3). This confirms the correctness of describing the longitudinal mixing process by the diffusion equation (3). Note that all theoretical concentration profiles shown in Fig. 6 (for a single operating condition) employ the same value for the diffusion coefficient. The dependence of the diffusion coefficient D, as determined from the results of the simulations using Eq. (8), on the non-dimensional operating parameters Fr and J is shown in Fig. 7. Note that the diffusion coefficient is presented in a non-dimensional form (compare Eq. (10)). Fig. 7 shows that the influence of Fr (i.e., the rotational speed) on the diffusion coefficient D is small. This means that the diffusion coefficient D ∗ will be proportional to , as
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
x/L
Fig. 6. Concentration profiles for various time instants (expressed in terms of the number of revolutions); results for Fr = 1.00 and J = 30%. The markers give the actual concentration as determined from the simulations, while the solid lines indicate the theoretical concentration profiles, see Eq. (5), with a (single) diffusion coefficient D.
-3
2
x 10
1.8 1.6 1.4 1.2
D /R 2
Fig. 5. Evolution of mean position of the “red” particles with time (expressed in terms of the number of revolutions, n). Results for Fr=1.00 and J = 30%.
-0.4
1 0.8 0.6 0.4 J=2 0% J=3 0% J=4 0%
0.2 0
0
0.5
1
1.5
2
2.5
Fr
Fig. 7. Dependence of non-dimensional diffusion coefficient on Froude number Fr and filling degree J.
follows from Eq. (4). Hence there is no diffusion when the rotary kiln does not rotate, = 0, as there is no equivalent of thermal movement of the particles in this case. For filling degree J = 20% and Fr 1.5, the diffusion coefficients D are higher than for J = 30% and 40%, while for Fr 1.0 they are practically the same. The diffusion coefficients for J = 30% and 40% are practically identical. Results for the dependence of the diffusion coefficient D ∗ on operating conditions were given by Dury and Ristow (1999) (using DEM simulations with J 55% and a binary mixture of grain sizes), Rutgers (1965) and Rao et al. (1991) (based on experiments) and Sheritt et al. (2003) (based on experiments and an analysis of data from the literature).
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 1
4089
3D simulations
0.9
J=20% J=30% J=40%
0.6
0.8
0.5
0.7
0.4
0.5
S
M(n)
0.6
0.3
0.4 0.3
0.2 0.2
0.1
0.1 0
0
5
10
15
20
25
30
35
40
n
0
0.5
1
1.5
2
2.5
Fr
Fig. 8. Evolution of mixing index M(n). The markers denote the mixing index as determined from the simulations, while the solid line corresponds to the exponential form (13) with a fitted parameter S. Results for Fr=1.00 and J = 30%.
Caution must be exercised in quantitative comparison of these results, since they have been obtained for different materials. The reported dependence of diffusion coefficient D ∗ on rotational speed is contradictory: Dury and Ristow (1999) give D ∗ ∝ 2.0 , Rao et al. (1991) find D ∗ ∝ 1.0 , while Rutgers (1965) and Sheritt et al. (2003) obtain D ∗ ∝ 0.5 and D ∗ ∝ 0.4 , respectively. The dependence observed here, D ∗ ∝ 1.0 , is identical to that of Rao et al. (1991). Rutgers (1965), Rao et al. (1991) and Sheritt et al. (2003) find that D ∗ generally decreases with increasing filling degree J , as is also noted here. In particular, Rutgers (1965) and Sheritt et al. (2003) give D ∗ ∝ J −0.5 in the rolling and cascading flow regime. This dependence is stronger than what is observed here. 4.3. Transverse mixing The evolution with number of revolutions n of the index M(n) for mixing in the transverse plane, as defined in Eq. (12), is illustrated in Fig. 8. This evolution can be accurately described by an exponential form M(n) 1 − exp(−Sn),
0
(13)
where the parameter S is a measure of the speed of the mixing process in the transverse plane (large S corresponds to quick mixing). Note that this transverse mixing is expressed in terms of the number of revolutions n of the rotary kiln, not in time t. The parameter S was determined by fitting the data from the simulations with the exponential form (13). The dependence of the parameter S on the nondimensional operating parameters Froude number Fr and filling degree J is shown in Fig. 9. This transverse mixing speed S (in terms of number of revolutions!) decreases
Fig. 9. Dependence of transverse mixing speed S, see Eq. (13), on operating parameters Fr and J in the three-dimensional DEM simulations.
with increasing rotational speed. This means that slow rotation, in comparison to faster rotation, enhances mixing in the transverse plane at identical number of revolutions (not identical time). Similarly, low filling degrees J also enhance transverse mixing. 4.4. Two-dimensional versus three-dimensional simulations: transverse mixing Additional two-dimensional DEM simulations, with far fewer particles but identical particle-size distribution and particle properties as described in Section 2, have been performed in which particle motions were restricted to a twodimensional transverse plane. In these two-dimensional simulations the filling degree is defined as the fraction of the area of the transverse plane that is filled with granular material (including voids between particles). The resulting flow regimes and transverse mixing speeds in these two-dimensional simulations have been studied, analogous to those for the three-dimensional case (see Sections 3.2, 4.1 and 4.3). This additional investigation has been performed to determine the effect of the dimensionality of the simulations (two-dimensional versus three-dimensional) on the flow regimes and the mixing speed. Note that the accuracy of the results in the two-dimensional case is smaller than in the three-dimensional case, since the number of particles in the simulation is fairly small (Np 500 for J =40%). The observed flow regimes in the two-dimensional simulations, see Table 2, differ from those observed in the threedimensional simulations, see Table 1. In the two-dimensional simulations the transitions, from rolling to cataracting and from cataracting to centrifuging, seem to occur at lower values of Froude number Fr than in the corresponding threedimensional simulations.
4090
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
Table 2 Flow regimes in two-dimensional simulations
J = 20% J = 30% J = 40%
Fr = 0.25
Fr = 0.56
Fr = 1.00
Fr = 1.56
Fr = 2.25
Rolling Rolling Rolling
Rolling Rolling Rolling
Rolling/cataracting Rolling/cataracting Cataracting
Cataracting Cataracting Cataracting/centrifuging
Cataracting Cataracting Centrifuging
2D simulations J=20% J=30% J=40%
1
0.8
S
0.6
0.4
0.2
0
0
0.5
1
1.5
2
2.5
Fr
Fig. 10. Dependence of transverse mixing speed S, see Eq. (13), on operating parameters Fr and J in the two-dimensional DEM simulations.
In comparison to the transverse mixing speed S in threedimensional DEM simulations (see Fig. 9), the mixing speed S in two-dimensional simulations is higher (see Fig. 10), although the trends are the same: decrease of mixing speed with increasing rotational speed and increasing filling degree. In the two-dimensional case particles have fewer neighbours (in comparison to the three-dimensional case) and hence are less constrained by other particles. This increased mobility is seen as the cause of the changes in the transitions between the flow regimes and of the increase in the transverse mixing speed S.
5. Discussion A computational study, using three-dimensional DEM simulations, has been performed of mixing in horizontal rotary kilns. Separate descriptions have been formulated for mixing in longitudinal and transverse directions. The transverse mixing is vigorous, while the longitudinal mixing is slow. DEM simulations have been performed for various values of the primary operating conditions, i.e., the rotational speed and the filling degree. The results of the simulations show that longitudinal mixing can be described by a one-dimensional diffusion equa-
tion. The dependence of the corresponding diffusion coefficient on the operating conditions has been determined. The diffusion coefficient D ∗ is proportional to the rotational speed , while its dependence on filling degree J is weak. Mixing in the transverse plane has been characterised by an “entropy”-like mixing index M(n). The results of the simulations show that the evolution with number of revolutions of this mixing index is exponential. From this exponential evolution, a transverse mixing speed S has been defined. The dependence of this mixing speed on the operating conditions has been determined from the simulations. It has been shown that, at identical number of revolutions, the speed of mixing in the transverse plane decreases with increasing rotational speed and filling degree J. A comparison of the transverse mixing speed S in two-dimensional and three-dimensional DEM simulations showed that this mixing speed is higher in the twodimensional case. Furthermore, transitions between flow regimes occur at different values of Froude number Fr. Thus, caution should be exercised when using two-dimensional simulations to quantity (three-dimensional) transverse mixing. The current investigation focused on the effect of the main operating conditions on the mixing behaviour, while keeping the particle properties constant. The result of the dimensional analysis for the diffusion coefficient, see Eq. (10), suggests that the non-dimensional parameters R/¯r and may also be of importance. Evidence for the influence of the interparticle friction coefficient is presented by Dury and Ristow (1999), while evidence for the influence of the parameter R/¯r is given by Rao et al. (1991). Further research on the influence of such parameters is recommended. A study of the effect of R/¯r will make it necessary to perform simulations with a larger number of particles. If this is done, a detailed study of the velocity profile in the active layer in the rolling regime can also be performed (see Boateng and Barr, 1997). The present number of particles in the simulations is not sufficient for this purpose. It is recommended to extend the DEM simulations by also considering non-spherical particles, since rotary kilns generally operate with feedstock consisting of angular particles. Furthermore, with non-spherical particles the rolling mechanism exhibited by spherical particles is suppressed. DEM simulations with non-spherical particles would lead to an increased effective friction coefficient (at an identical interparticle friction). A major disadvantage of DEM simulations with non-spherical particles is the increased computational
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
complexity. A promising approach to this problem has recently been reported by Kuhn (2003), who uses approximate axisymmetric ellipsoids, the so-called “ovoids”. The computational complexity of a DEM simulation employing these ovoidal particles is about twice that of a similar DEM simulation using spherical particles. Future work will also deal with inclined rotary kilns, with models for gas flow through rotary kilns and with experimental studies of longitudinal and transverse mixing in rotary kilns. Acknowledgements The authors acknowledge financial support of the first author, during his stay at the University of Twente, by IMPACT, the research institute of the University of Twente on Mechanics, Processes and Control. References Boateng, A.A., Barr, P.V., 1997. Granular flow behaviour in the transverse plane of a partially filled rotating cylinder. Journal of Fluid Mechanics 330, 233–249. Cundall, P.A., Strack, O.D.L., 1979a. A discrete numerical model for granular assemblies. Géotechnique 9, 47–65. Cundall, P.A., Strack, O.D.L., 1979b. The distinct element method as a tool for research in granular media, Part II. Report to the National Science Foundation, ENG76-20711, Department of Civil and Mineral Engineering, University of Minnesota, MN, USA. Ding, Y.L., Forster, R.N., Seville, J.P.K., Parker, D.J., 2001. Scaling relationships for rotating drums. Chemical Engineering Science 56, 3737–3750. Dury, C.M., Ristow, G.H., 1999. Axial particle diffusion in rotating cylinders. Granular Matter 1, 151–161. Henein, H., Brimacombe, J.K., Watkinson, A.P., 1983a. Experimental study of transverse bed motion in rotary kilns. Metallurgical Transactions B 14B, 191–205. Henein, H., Brimacombe, J.K., Watkinson, A.P., 1983b. The modelling of transverse solids motion in rotary kilns. Metallurgical Transactions B 14B, 207–220. Heydenrych, M.D., 2001. Modelling of rotary kilns. Ph.D. Thesis, Department of Chemical Engineering, University of Twente, The Netherlands. Hogg, R., Cahn, D.S., Healy, T.W., Fuerstenau, D.W., 1966. Diffusional mixing in an ideal system. Chemical Engineering Science 21, 1025–1038.
4091
Hoomans, B.P.B., 2000. Granular dynamics of gas-solid two-phase flows. Ph.D. Thesis, Department of Chemical Engineering, University of Twente, The Netherlands. Kapur, J.N., Kesavan, H.K., 1992. Entropy Optimization Principles with Applications. Academic Press, San Diego, CA, USA. Katz, A., 1967. Principles of Statistical Mechanics: The Information Theory Approach. W.H. Freeman and Company, San Francisco, CA, USA. Kohring, G.A., 1995. Studies of diffusional mixing in rotating drums via computer simulations. Journal de Physique 15, 1551–1571. Kuhn, M.R., 2003. Smooth convex three-dimensional particle for the Discrete Element Method. Journal of Engineering Mechanics (Transactions of the ASCE) 129, 539–547. McCarthy, J.J., Khakhar, D.V., Ottino, J.M., 2000. Computational studies of granular mixing. Powder Technology 109, 72–82. Mellmann, J., 2001. The transverse motion of solids in rotating cylinders—forms of motion and transition behaviour. Powder Technology 118, 251–270. Ottino, J.M., Khakhar, D.V., 2000. Mixing and segregation of granular materials. Annual Review of Fluid Mechanics 32, 55–91. Rao, S.J., Bhatia, S.K., Khakhar, D.V., 1991. Axial transport of granular solids in rotating cylinders. Part. 2: Experiments in a non-flow system. Powder Technology 67, 153–162. Rothenburg, L., Bathurst, R.J., 1991. Numerical simulation of idealised granular assemblies with plane elliptical particles. Computers and Geotechnics 11, 315–329. Rutgers, R., 1965. Longitudinal mixing of granular materials flowing though a rotating cylinder. Chemical Engineering Science 20, 1097–1100. Schutyser, M.A.I., Padding, J.T., Weber, F.J., Briels, W.J., Rinzema, A., Boom, R., 2001. Discrete particle simulations predicting mixing behaviour of solid substrate in a rotating drum fermenter. Biotechnology and Bioengineering 75, 666–675. Schutyser, M.A.I., Weber, F.J., Briels, W.J., Boom, R.M., Rinzema, A., 2002. Three-dimensional simulation of grain mixing in three different rotating drum designs for solid-state fermentation. Biotechnology and Bioengineering 79, 284–294. Sheritt, R.G., Chaouki, J., Mehrotra, A.K., Behie, L.A., 2003. Axial dispersion in the three-dimensional mixing of particles in a rotating drum reactor. Chemical Engineering Science 58, 401–415. Shinbrot, T., Alexander, A., Muzzio, F.J., 1999. Spontaneous chaotic granular mixing. Nature 397, 675–678. Skinner, A.E., 1969. A note on the influence of interparticle friction on the shearing strength of a random assembly of spherical particles. Géotechnique 1, 150–157. Sudah, O.S., Chester, A.W., Kowalski, J.A., Beeckman, J.W., Muzzio, F.J., 2002. Quantitative characterization of mixing processes in rotary calciners. Powder Technology 126, 166–173. Yang, R.Y., Zou, R.P., Yu, A.B., 2003. Microdynamic analysis of particle flow in a horizontal rotating drum. Powder Technology 130, 138–146.
Longitudinal and transverse mixing in rotary kilns: A discrete element method approach G.J. Finniea , N.P. Kruytb,∗ , M. Yec , C. Zeilstrac , J.A.M. Kuipersc a Department of Chemical Engineering, Soshanguve Campus, Tshwane University of Technology, Private Bag X07, Pretoria North, Republic of South Africa b Department of Mechanical Engineering, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands c Department of Chemical Engineering, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands
Received 13 September 2004; received in revised form 20 October 2004; accepted 4 December 2004 Available online 7 April 2005
Abstract Longitudinal and transverse mixing in horizontal rotary kilns has been studied, using a three-dimensional Discrete Element Method approach. The focus is on the effect of the main operating conditions, i.e., the filling degree and the rotational speed of the rotary kiln, on quantitative measures of longitudinal and transverse mixing. Discrete Element Method simulations have been performed for various values of the main operating conditions. Flow regimes have been determined from visualisations of the simulations. It is shown that the mixing in the longitudinal direction can be described by a diffusion equation. The dependence of the diffusion coefficient on the operating conditions has been determined from the results of the simulations. It is found that the diffusion coefficient increases linearly with rotational speed, while the influence of the filling degree is relatively small. The mixing in the transverse plane is quantified by an entropy-like mixing index. The results of the simulations show that this mixing index varies exponentially with the number of revolutions of the rotary kiln. The dependence on the operating conditions of the transverse mixing speed, as determined from the exponential behaviour of the mixing index, has been determined. This mixing speed decreases with increasing rotational speed and with increasing filling degree. The transverse mixing speed that has been determined from the results of additional two-dimensional simulations is generally larger than that observed in the corresponding three-dimensional simulations. 䉷 2005 Elsevier Ltd. All rights reserved. Keywords: Rotary kiln; Granular material; Mixing; Discrete Element Method; Diffusion; Simulation
1. Introduction Rotary kilns are widely used in chemical and mechanical engineering for processing of granular materials for various purposes, such as for mixing, drying and coating of materials, for calcining of limestone, for reduction of iron ore and for waste incineration. Rotary kilns are able to handle varying feedstocks, in particular with
∗ Corresponding author. Tel.: +31 53 489 2528; fax: +31 53 489 3695.
E-mail addresses: fi[email protected] (G.J. Finnie), [email protected] (N.P. Kruyt), [email protected] (M. Ye), [email protected] (C. Zeilstra), [email protected] (J.A.M. Kuipers). 0009-2509/$ - see front matter 䉷 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.12.048
wide particle-size distributions and significant differences in physical properties. The operating conditions of a rotary kiln (as sketched in Fig. 1) involve rotational speed around its axis, radius R, length L in longitudinal direction, inclination with respect to the horizontal direction and filling degree J (i.e., the fraction of the volume of the rotary kiln that is filled with granular material, including voids between particles), as well as parameters related to the feedstock. In the sequel the focus is on horizontal rotary kilns, i.e., = 0◦ . Inclined rotary kilns were studied by Sudah et al. (2002). Generally, the main operating conditions of rotary kilns are expressed in terms of non-dimensional numbers by the filling degree J and the Froude number Fr, which is
4084
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
nificant segregation is expected. A recent overview of mixing and segregation of granular materials is given by Ottino and Khakhar (2000). A computational DEM study of mixing of binary particle mixtures (including density differences between the two sizes) is reported by Dury and Ristow (1999). At present, gas flows in the rotary kilns are not considered here. Separate descriptions are formulated for the longitudinal and transverse mixing. The results of the DEM simulations are analysed to determine quantitatively the influence of nondimensional operating conditions of the rotary kiln on the parameters describing the longitudinal and transverse mixing. The outline of this study is as follows. Firstly, the DEM simulations are described in Section 2. This is followed in Section 3 by the formulation of measures to quantify longitudinal and transverse mixing. The results of the simulations are analysed in Section 4. Finally, findings from this study are discussed.
Fig. 1. Sketch of the geometry of the horizontal rotary kiln.
defined by
R , g
2. Discrete element method simulations
2
Fr =
(1)
where g is the acceleration due to gravity. Two different types of mixing can be distinguished in rotary kilns. The first type occurs in the plane transverse to the axis of rotation, while the second type occurs in longitudinal direction. The first type of mixing is vigorous, while the second is much slower. Depending on the non-dimensional operating conditions, several flow regimes may be observed in the plane transverse to the axis of rotation, namely slipping, slumping, rolling, cascading, cataracting and centrifuging (Henein et al., 1983a,b; Mellmann, 2001). In industrial operations the rolling regime is mostly encountered, since it promotes good mixing of particles. This flow regime is characterised by two distinct regions: (i) a relatively thin active layer is formed at the top in which the material flows along the inclined slope and (ii) a thick passive region in which the material rotates as a rigid body with the angular velocity of the rotary kiln (see for example, Boateng and Barr, 1997). Scaling relationships for rotary kilns in the rolling and slumping regime have been studied by Ding et al. (2001). The present study is a basic computational study of longitudinal and transverse mixing in horizontal rotary kilns, focusing on the effect of the main operating conditions on quantitative measures of longitudinal and transverse mixing. A recent overview of computational methods for studying mixing of granular materials is given by McCarthy et al. (2000). The method used here to study the longitudinal and transverse mixing processes is that of three-dimensional Discrete Element Method (DEM) simulations, as proposed by Cundall and Strack (1979a,b), for a number of nondimensional operating conditions. The particle-size distribution employed in the simulations is fairly narrow, so no sig-
The computational method employed here is the threedimensional DEM (also called Discrete Particle Method), as proposed by Cundall and Strack (1979a,b). In this method the motion of a large number of particles is determined numerically by integrating the equations of motion, using a suitable model, the so-called contact constitutive relation, for the particle–particle interactions and the particle–wall interactions. Here, equations of motions are solved for particle translations as well as for particle rotations, contrary to what was done by Schutyser et al. (2001). Particle shapes are restricted to spheres (with varying radii) for reasons of computational efficiency. The contact constitutive relation employed here is basically that of Cundall and Strack (1979a,b). It involves (linear) elastic behaviour, Coulomb frictional behaviour and viscous contact damping to model dissipative collisions between particles. A linear spring is acting in the direction normal to the contact and a tangential spring is acting in the direction of the tangential component of the relative displacement at the contact. The respective spring constants are denoted by kn and kt . Coulomb friction is modelled by an interparticle friction coefficient , such that at each contact c we have ftc |fnc | where fnc is the (scalar) normal component of the contact force vector f c , i.e., fnc = f c · n c with n c the unit normal vector at the contact c, ftc is the magnitude of the tangential component of the force, i.e., ftc = f c − fnc n c and |fnc | is the absolute value of the normal component of the contact force. Viscous contact damping is implemented in terms of a coefficient of restitution cr (see for example, Hoomans, 2000). A similar contact constitutive relation is used for the interactions between particles and cylinder walls and end walls of the rotary kiln, but with a different value for the friction coefficient . Unlike in the simulations by
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
Yang et al. (2003), no rolling resistance was included here in the contact constitutive relation, considering the problem of establishing a manner for determining the material parameters involved in the rolling resistance model. The radii r of the spherical particles follow a (truncated) lognormal particle-size distribution with rmax /rmin = 1.4. The ratio of the standard deviation over the mean of the particle-size distribution equals 0.04. Hence, the particlesize distribution is fairly narrow, so no significant size segregation is expected. Furthermore, all particles have the same density, so no density segregation will occur. The elastic parameter kn has been selected such that the “overlap” between particles is small, yet the time step in the (explicit) numerical integration of the equations of motion does not become exceedingly small. The stiffness ratio kt /kn = 0.29. A fairly high interparticle friction coefficient is 0.5 is selected, since real particles are not perfectly spherical (as is assumed here). Non-spherical particles have a higher effective (continuum) friction, in comparison to spherical particles with the same interparticle friction coefficient, due to suppression of rolling (Rothenburg and Bathurst, 1991). A high friction coefficient of 1.5 is used for the interactions between particles and cylinder wall, corresponding to a roughened wall. The friction coefficient of the end walls of the rotary kiln is 0.01. The selected restitution coefficient cr is 0.71. Note that the value of the (microscopic) interparticle friction coefficient and the so-called (macroscopic) internal friction coefficient of the granular material as a continuum (related to the static angle of repose of the granular material) are generally not equal. In fact, there is no direct relation between these two friction coefficients (Skinner, 1969). Similarly, the value of the particle–wall friction coefficient is generally not equal to the wall friction coefficient of the granular material as a continuum. The geometry of the rotary kiln (see also Fig. 1) is such that its aspect ratio, L/R, equals 2.0. The ratio of the kiln radius R over the average particle radius r¯ is R/¯r = 40.0. Several assemblies were simulated, with the number of particles in the simulations, Np , depending on the filling degree J of the kiln: for J = 20%, Np = 10, 667, for J = 30%, Np = 16, 000 and for J = 40%, Np = 21, 333. These numbers of particles were determined from considerations of acceptable computer simulation times. The time step in the explicit numerical integration of the equations of motion is such that it is a factor 30 smaller than the maximum time step allowed for stability of the numerical method. Five values for the rotational speed of the kiln, , have been considered. These correspond to Froude numbers, see (1), of Fr = 0.25, 0.56, 1.00, 1.56, 2.25. In total, 15 simulations have been performed for the five values of Froude number Fr, and for each Froude number with the three filling degrees J = 20%, 30%, 40%. For each simulation the simulated time corresponds to 40 revolutions of the rotary kiln.
4085
3. Description of longitudinal and transverse mixing Two different types of mixing occur in rotary kilns. The first type is mixing in the transverse plane. This type of mixing is vigorous and quickly leads to mixed-out states. The second, much slower type is mixing in the longitudinal direction. Separate descriptions are formulated for these two types of mixing. 3.1. Longitudinal mixing In horizontal rotary kilns there is no net transport in the longitudinal direction. Motion of the particles in longitudinal direction is the result of interparticle collisions induced by the rotation of the rotary kiln. This motion is therefore expected to be completely random. Consequently, the mixing process is considered as analogous to diffusion, and hence Fick’s law should apply. This leads to the one-dimensional diffusion equation (Hogg et al., 1966) j j jc(x, t) c(x, t) = D∗ jt jx jx
(2)
for practically identical particles, where c(x, t) is the relative concentration of the mixture after time t at position x in longitudinal direction and D ∗ (with units m2 s−1 ) is the diffusion coefficient that is assumed to be spatially uniform in the sequel. By performing the transformation n = (t)/(2), this equation can be expressed in terms of the (non-dimensional) number of revolutions, n, of the rotary kiln instead of time t. The resulting equation becomes j j jc(x, n) c(x, n) = D , jn jx jx
(3)
where the units of D are m2 . The diffusion coefficients D ∗ and D are related by D=
2 ∗ D .
(4)
With initial conditions c(x, 0) = 0 for −L/2 < x < 0 and c(x, 0)=1 for 0 x < L/2 (corresponding to “blue” particles on the left and “red” particles to the right of the interface at x = 0; see also the sketch of the initial configuration for the longitudinal mixing on the left-hand side of Fig. 2) and boundary conditions jc(x, n)/jx|x=±L/2 =0 (corresponding to zero flux at the boundaries), the solution of the diffusion equation (3) is (Hogg et al., 1966) ∞ 1 2 1 −(2j − 1)2 2 Dn c(x, n) = + exp 2 2j − 1 L2 j =1
(2j − 1)x × sin . L
(5)
4086
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 Longitudinal
“blue”
3.2. Transverse mixing
Transverse
“red”
“gray”
“black”
Fig. 2. Sketch of initial configurations for longitudinal mixing (left) and transverse mixing (right).
The (time-dependent) non-dimensional mean position of the “red” particles is defined by L/2 1 xc(x, n) dx. (6)
c(n) = 2 L −L/2 Using (5), it follows (Kohring, 1995) that c(n) is given by ∞ 8 (−1)j +1 2 2 Dn
c(n) = 3 exp −(2j − 1) 2 . (7) L (2j − 1)3 j =1
After some algebra it follows that for small values of the Fourier number (Dn)/L2 , this expression can be approximated by 1 Dn
c(n) − 2 2 4 L
for
Dn 1. L2
(8)
Thus, for small Fourier number, the relationship between mean position of the “red” particles and (non-dimensional) time is linear. It is expected that the diffusion coefficient D depends on the process parameters and J , the gravitational acceleration g, as well as on the geometry of the rotary kiln, i.e., its radius R and length L. Furthermore, it assumed that of the material properties only the average particle size r¯ , density and frictional material parameters (collectively denoted by here) are important. This means that the influence of elastic and damping parameters of the granular material is excluded. Thus we obtain the following functional relationship: D = F (, J, g, R, L, r¯ , , ).
(9)
From dimensional analysis we find the functional relationships for D and D ∗ D R L , = G Fr, J, , , R2 r¯ R D∗ R L ∗ , . (10) Fr, J, = G , R 2 r¯ R An analogous dimensional analysis for the motion in the transverse plane was performed by Ding et al. (2001).
In analogy to the case of longitudinal mixing, as described in Section 3.1, Heydenrych (2001) solved a simplified convection-diffusion equation for the mixing in the transverse plane. This complex approach was not pursued here, considering the spatial variations in velocity that have to be taken into account in the convective part of the convectiondiffusion equation for the concentration c(x, n) and the possible spatial variations of the diffusion coefficient D in the transverse plane. Instead, transverse mixing is quantified here using an “entropy”-like quantity, as proposed by Schutyser et al. (2001, 2002). In this approach a grid with dimensions (I, J ) is defined that overlaps the region of interest. For each grid cell (i, j ), an “entropy” Si,j is calculated as Si,j = −kV i,j [pi,j ln pi,j + qi,j ln qi,j ],
(11)
where k is a constant, Vi,j is the volume occupied by the particles in cell (i, j ), pi,j is the concentration (fraction) of “black” particles, while qi,j = 1 − pi,j is the concentration (fraction) of “grey” particles. By definition, pi,j is the ratio of the total volume in cell (i, j ) occupied by “black” particles over the total volume in cell (i, j ) occupied by “black” and “grey” particles. Hence, the voids between particles do not affect the concentrations pi,j and qi,j . Note that many particles, with differing longitudinal position, may be mapped into the same cell in the transverse plane. Initially, the colours of the particles are set to “grey” when they are located in the left-hand side (when viewed along the kiln axis) of the rotary kiln, and to “black” when they are located in the right-hand side; see also the sketch of the initial configuration for the transverse mixing on the right-hand side of Fig. 2. Note that this “entropy” Si,j is not related to the thermodynamical entropy of the granular material. Such information-theoretical entropies are used in many applications (Katz, 1967; Kapur and Kesavan, 1992). A “mixing index” M(n) is now defined by M(n) =
Si,j (n) i,j
S∞
,
(12)
where the “entropy” of the perfectly mixed state, S∞ , is given by S∞ = −kV [p0 ln p0 + q0 ln q0 ], V is the total volume occupied by the particles, p0 is the concentration (fraction) of “black” particles in the system and q0 =1−p0 is the concentration (fraction) of “grey” particles in the system. In the initial state the mixing index M(0) = 0, while in the perfectly mixed state that is attained as n → ∞, the mixing index equals M(n → ∞)=1. Note that the constant k in Eq. (11) cancels in expression (12) for the mixing index M(n). In the determination of the concentrations pi,j , it was taken into account that particles may occupy multiple cells. Thus, the fractional volumes of particles inside the cells were used to determine pi,j and qi,j , instead of basing their determination solely on the coordinates of particle centres.
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 Kiln speed, Ω⋅R
Fig. 3. Employed fractional volumes in the determination of the cell concentrations pi,j and qi,j .
This procedure is shown schematically in Fig. 3 for disks. The (approximate) computation of the fractional volumes is straightforward (Hoomans, 2000), provided that the length of the cells is larger than the maximum particle diameter. This procedure reduces the dependence of the computed “entropy” on the grid sizes (especially when the number of particles in each cell is fairly small, such as in twodimensional simulations) that was reported by Schutyser et al. (2001).
4. Results The results of the three-dimensional DEM simulations have been analysed with respect to occurring flow regimes, longitudinal mixing and transverse mixing. Additional twodimensional DEM simulations have been performed in order to compare the transverse mixing in two-dimensional and three-dimensional simulations. 4.1. Flow regimes Depending on the non-dimensional process parameters, such as Froude number Fr and filling degree J, various flow regimes may be distinguished: sliding, surging, slumping, rolling, cascading, cataracting and centrifuging. An overview of the flow regimes, together with transition criteria between the flow regimes, is given by Mellmann (2001).
4087 Kiln speed, Ω⋅R
Fig. 4. Velocity of particles in a transverse cross-section. Left: rolling regime for J = 40% and Fr = 0.25; right: cataracting regime for J = 40% and Fr = 2.25. Note that the magnitude of the velocity vectors has been scaled for visual clarity. For comparison, the kiln speed, R, has also been plotted at the top of the figures.
Based on analyses of the particle-velocity profiles and on visualisations of the results of the simulations, the flow regimes were identified, see Table 1. Examples of the observed velocity profiles are shown in Fig. 4 for the rolling regime (left) and for the cataracting regime (right). A direct quantitative comparison of these results with the transition criteria given by Mellmann (2001) is not possible, since these criteria involve the continuum internal friction and wall friction coefficients. As noted in Section 2, these continuum friction coefficients are not (necessarily) equal to the particle–particle and the particle–wall friction coefficients, respectively. Note that the transition from the cataracting regime to the centrifuging regime does not have to correspond to Fr=1, as is discussed by Mellmann (2001). The centrifuging regime has not been clearly observed in the simulations. This is attributed to the low effective friction between the granular material and the wall of the rotary kiln that is caused by the rolling tendency of spheres and disks (see for example, Rothenburg and Bathurst, 1991). Shinbrot et al. (1999) experimentally observed chaotic transverse mixing for very fine particles. Since the particles considered here are fairly large, this chaotic mixing was not observed here. 4.2. Longitudinal mixing In the initial configuration at n=0, all particles on the left were “blue” and all particles to the right of the interface at
Table 1 Flow regimes in three-dimensional simulations
J = 20% J = 30% J = 40%
Fr = 0.25
Fr = 0.56
Fr = 1.00
Fr = 1.56
Fr = 2.25
Rolling Rolling Rolling
Rolling Rolling Rolling
Rolling/cataracting Rolling/cataracting Rolling/cataracting
Cataracting Cataracting Cataracting
Cataracting Cataracting Cataracting
4088
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
0.25
1 0.9
0.245
0.8 0.7
0.24
c(x,n)
0.6 0.235
0.5 0.4
0.23
0.3 0.2
0.225
n=10 revs n=20 revs n=30 revs n=40 revs
0.1 0.22 0
5
10
15
20
25
30
35
0 -0.5
40
n
x =0 were “red” (see also the left-hand side of Fig. 2). Since the evolution of all particle positions is determined in the DEM simulations, concentration profiles can be determined by dividing the longitudinal extent of the rotary kiln into slices and computing, for each time instant, the volume fraction of “red” particles in each slice. Thus, time-dependent (with time expressed in terms of the number of revolutions, n) concentrations profiles c(x, n) were obtained, as well as time-dependent mean positions c(n) of the “red” particles. The diffusion coefficients D were determined from the simulations by considering the evolution of the mean position c(n) with number of revolutions n. According to Eq. (8), c(n) depends linearly on Fourier number (Dn)/L2 for small (Dn)/L2 . The results of the simulations confirm this linear relationship between c(n) and n, see Fig. 5. According to the theory presented in Section 3.1 for the longitudinal mixing, the concentration profiles are given by Eq. (5). A typical comparison of the actual and the theoretical concentration profiles is given in Fig. 6. This Figure shows that the actual concentration profiles, as determined from the simulations, closely correspond to the theoretical concentration profiles that follow from the solution of the diffusion equation (3). This confirms the correctness of describing the longitudinal mixing process by the diffusion equation (3). Note that all theoretical concentration profiles shown in Fig. 6 (for a single operating condition) employ the same value for the diffusion coefficient. The dependence of the diffusion coefficient D, as determined from the results of the simulations using Eq. (8), on the non-dimensional operating parameters Fr and J is shown in Fig. 7. Note that the diffusion coefficient is presented in a non-dimensional form (compare Eq. (10)). Fig. 7 shows that the influence of Fr (i.e., the rotational speed) on the diffusion coefficient D is small. This means that the diffusion coefficient D ∗ will be proportional to , as
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
x/L
Fig. 6. Concentration profiles for various time instants (expressed in terms of the number of revolutions); results for Fr = 1.00 and J = 30%. The markers give the actual concentration as determined from the simulations, while the solid lines indicate the theoretical concentration profiles, see Eq. (5), with a (single) diffusion coefficient D.
-3
2
x 10
1.8 1.6 1.4 1.2
D /R 2
Fig. 5. Evolution of mean position of the “red” particles with time (expressed in terms of the number of revolutions, n). Results for Fr=1.00 and J = 30%.
-0.4
1 0.8 0.6 0.4 J=2 0% J=3 0% J=4 0%
0.2 0
0
0.5
1
1.5
2
2.5
Fr
Fig. 7. Dependence of non-dimensional diffusion coefficient on Froude number Fr and filling degree J.
follows from Eq. (4). Hence there is no diffusion when the rotary kiln does not rotate, = 0, as there is no equivalent of thermal movement of the particles in this case. For filling degree J = 20% and Fr 1.5, the diffusion coefficients D are higher than for J = 30% and 40%, while for Fr 1.0 they are practically the same. The diffusion coefficients for J = 30% and 40% are practically identical. Results for the dependence of the diffusion coefficient D ∗ on operating conditions were given by Dury and Ristow (1999) (using DEM simulations with J 55% and a binary mixture of grain sizes), Rutgers (1965) and Rao et al. (1991) (based on experiments) and Sheritt et al. (2003) (based on experiments and an analysis of data from the literature).
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091 1
4089
3D simulations
0.9
J=20% J=30% J=40%
0.6
0.8
0.5
0.7
0.4
0.5
S
M(n)
0.6
0.3
0.4 0.3
0.2 0.2
0.1
0.1 0
0
5
10
15
20
25
30
35
40
n
0
0.5
1
1.5
2
2.5
Fr
Fig. 8. Evolution of mixing index M(n). The markers denote the mixing index as determined from the simulations, while the solid line corresponds to the exponential form (13) with a fitted parameter S. Results for Fr=1.00 and J = 30%.
Caution must be exercised in quantitative comparison of these results, since they have been obtained for different materials. The reported dependence of diffusion coefficient D ∗ on rotational speed is contradictory: Dury and Ristow (1999) give D ∗ ∝ 2.0 , Rao et al. (1991) find D ∗ ∝ 1.0 , while Rutgers (1965) and Sheritt et al. (2003) obtain D ∗ ∝ 0.5 and D ∗ ∝ 0.4 , respectively. The dependence observed here, D ∗ ∝ 1.0 , is identical to that of Rao et al. (1991). Rutgers (1965), Rao et al. (1991) and Sheritt et al. (2003) find that D ∗ generally decreases with increasing filling degree J , as is also noted here. In particular, Rutgers (1965) and Sheritt et al. (2003) give D ∗ ∝ J −0.5 in the rolling and cascading flow regime. This dependence is stronger than what is observed here. 4.3. Transverse mixing The evolution with number of revolutions n of the index M(n) for mixing in the transverse plane, as defined in Eq. (12), is illustrated in Fig. 8. This evolution can be accurately described by an exponential form M(n) 1 − exp(−Sn),
0
(13)
where the parameter S is a measure of the speed of the mixing process in the transverse plane (large S corresponds to quick mixing). Note that this transverse mixing is expressed in terms of the number of revolutions n of the rotary kiln, not in time t. The parameter S was determined by fitting the data from the simulations with the exponential form (13). The dependence of the parameter S on the nondimensional operating parameters Froude number Fr and filling degree J is shown in Fig. 9. This transverse mixing speed S (in terms of number of revolutions!) decreases
Fig. 9. Dependence of transverse mixing speed S, see Eq. (13), on operating parameters Fr and J in the three-dimensional DEM simulations.
with increasing rotational speed. This means that slow rotation, in comparison to faster rotation, enhances mixing in the transverse plane at identical number of revolutions (not identical time). Similarly, low filling degrees J also enhance transverse mixing. 4.4. Two-dimensional versus three-dimensional simulations: transverse mixing Additional two-dimensional DEM simulations, with far fewer particles but identical particle-size distribution and particle properties as described in Section 2, have been performed in which particle motions were restricted to a twodimensional transverse plane. In these two-dimensional simulations the filling degree is defined as the fraction of the area of the transverse plane that is filled with granular material (including voids between particles). The resulting flow regimes and transverse mixing speeds in these two-dimensional simulations have been studied, analogous to those for the three-dimensional case (see Sections 3.2, 4.1 and 4.3). This additional investigation has been performed to determine the effect of the dimensionality of the simulations (two-dimensional versus three-dimensional) on the flow regimes and the mixing speed. Note that the accuracy of the results in the two-dimensional case is smaller than in the three-dimensional case, since the number of particles in the simulation is fairly small (Np 500 for J =40%). The observed flow regimes in the two-dimensional simulations, see Table 2, differ from those observed in the threedimensional simulations, see Table 1. In the two-dimensional simulations the transitions, from rolling to cataracting and from cataracting to centrifuging, seem to occur at lower values of Froude number Fr than in the corresponding threedimensional simulations.
4090
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
Table 2 Flow regimes in two-dimensional simulations
J = 20% J = 30% J = 40%
Fr = 0.25
Fr = 0.56
Fr = 1.00
Fr = 1.56
Fr = 2.25
Rolling Rolling Rolling
Rolling Rolling Rolling
Rolling/cataracting Rolling/cataracting Cataracting
Cataracting Cataracting Cataracting/centrifuging
Cataracting Cataracting Centrifuging
2D simulations J=20% J=30% J=40%
1
0.8
S
0.6
0.4
0.2
0
0
0.5
1
1.5
2
2.5
Fr
Fig. 10. Dependence of transverse mixing speed S, see Eq. (13), on operating parameters Fr and J in the two-dimensional DEM simulations.
In comparison to the transverse mixing speed S in threedimensional DEM simulations (see Fig. 9), the mixing speed S in two-dimensional simulations is higher (see Fig. 10), although the trends are the same: decrease of mixing speed with increasing rotational speed and increasing filling degree. In the two-dimensional case particles have fewer neighbours (in comparison to the three-dimensional case) and hence are less constrained by other particles. This increased mobility is seen as the cause of the changes in the transitions between the flow regimes and of the increase in the transverse mixing speed S.
5. Discussion A computational study, using three-dimensional DEM simulations, has been performed of mixing in horizontal rotary kilns. Separate descriptions have been formulated for mixing in longitudinal and transverse directions. The transverse mixing is vigorous, while the longitudinal mixing is slow. DEM simulations have been performed for various values of the primary operating conditions, i.e., the rotational speed and the filling degree. The results of the simulations show that longitudinal mixing can be described by a one-dimensional diffusion equa-
tion. The dependence of the corresponding diffusion coefficient on the operating conditions has been determined. The diffusion coefficient D ∗ is proportional to the rotational speed , while its dependence on filling degree J is weak. Mixing in the transverse plane has been characterised by an “entropy”-like mixing index M(n). The results of the simulations show that the evolution with number of revolutions of this mixing index is exponential. From this exponential evolution, a transverse mixing speed S has been defined. The dependence of this mixing speed on the operating conditions has been determined from the simulations. It has been shown that, at identical number of revolutions, the speed of mixing in the transverse plane decreases with increasing rotational speed and filling degree J. A comparison of the transverse mixing speed S in two-dimensional and three-dimensional DEM simulations showed that this mixing speed is higher in the twodimensional case. Furthermore, transitions between flow regimes occur at different values of Froude number Fr. Thus, caution should be exercised when using two-dimensional simulations to quantity (three-dimensional) transverse mixing. The current investigation focused on the effect of the main operating conditions on the mixing behaviour, while keeping the particle properties constant. The result of the dimensional analysis for the diffusion coefficient, see Eq. (10), suggests that the non-dimensional parameters R/¯r and may also be of importance. Evidence for the influence of the interparticle friction coefficient is presented by Dury and Ristow (1999), while evidence for the influence of the parameter R/¯r is given by Rao et al. (1991). Further research on the influence of such parameters is recommended. A study of the effect of R/¯r will make it necessary to perform simulations with a larger number of particles. If this is done, a detailed study of the velocity profile in the active layer in the rolling regime can also be performed (see Boateng and Barr, 1997). The present number of particles in the simulations is not sufficient for this purpose. It is recommended to extend the DEM simulations by also considering non-spherical particles, since rotary kilns generally operate with feedstock consisting of angular particles. Furthermore, with non-spherical particles the rolling mechanism exhibited by spherical particles is suppressed. DEM simulations with non-spherical particles would lead to an increased effective friction coefficient (at an identical interparticle friction). A major disadvantage of DEM simulations with non-spherical particles is the increased computational
G.J. Finnie et al. / Chemical Engineering Science 60 (2005) 4083 – 4091
complexity. A promising approach to this problem has recently been reported by Kuhn (2003), who uses approximate axisymmetric ellipsoids, the so-called “ovoids”. The computational complexity of a DEM simulation employing these ovoidal particles is about twice that of a similar DEM simulation using spherical particles. Future work will also deal with inclined rotary kilns, with models for gas flow through rotary kilns and with experimental studies of longitudinal and transverse mixing in rotary kilns. Acknowledgements The authors acknowledge financial support of the first author, during his stay at the University of Twente, by IMPACT, the research institute of the University of Twente on Mechanics, Processes and Control. References Boateng, A.A., Barr, P.V., 1997. Granular flow behaviour in the transverse plane of a partially filled rotating cylinder. Journal of Fluid Mechanics 330, 233–249. Cundall, P.A., Strack, O.D.L., 1979a. A discrete numerical model for granular assemblies. Géotechnique 9, 47–65. Cundall, P.A., Strack, O.D.L., 1979b. The distinct element method as a tool for research in granular media, Part II. Report to the National Science Foundation, ENG76-20711, Department of Civil and Mineral Engineering, University of Minnesota, MN, USA. Ding, Y.L., Forster, R.N., Seville, J.P.K., Parker, D.J., 2001. Scaling relationships for rotating drums. Chemical Engineering Science 56, 3737–3750. Dury, C.M., Ristow, G.H., 1999. Axial particle diffusion in rotating cylinders. Granular Matter 1, 151–161. Henein, H., Brimacombe, J.K., Watkinson, A.P., 1983a. Experimental study of transverse bed motion in rotary kilns. Metallurgical Transactions B 14B, 191–205. Henein, H., Brimacombe, J.K., Watkinson, A.P., 1983b. The modelling of transverse solids motion in rotary kilns. Metallurgical Transactions B 14B, 207–220. Heydenrych, M.D., 2001. Modelling of rotary kilns. Ph.D. Thesis, Department of Chemical Engineering, University of Twente, The Netherlands. Hogg, R., Cahn, D.S., Healy, T.W., Fuerstenau, D.W., 1966. Diffusional mixing in an ideal system. Chemical Engineering Science 21, 1025–1038.
4091
Hoomans, B.P.B., 2000. Granular dynamics of gas-solid two-phase flows. Ph.D. Thesis, Department of Chemical Engineering, University of Twente, The Netherlands. Kapur, J.N., Kesavan, H.K., 1992. Entropy Optimization Principles with Applications. Academic Press, San Diego, CA, USA. Katz, A., 1967. Principles of Statistical Mechanics: The Information Theory Approach. W.H. Freeman and Company, San Francisco, CA, USA. Kohring, G.A., 1995. Studies of diffusional mixing in rotating drums via computer simulations. Journal de Physique 15, 1551–1571. Kuhn, M.R., 2003. Smooth convex three-dimensional particle for the Discrete Element Method. Journal of Engineering Mechanics (Transactions of the ASCE) 129, 539–547. McCarthy, J.J., Khakhar, D.V., Ottino, J.M., 2000. Computational studies of granular mixing. Powder Technology 109, 72–82. Mellmann, J., 2001. The transverse motion of solids in rotating cylinders—forms of motion and transition behaviour. Powder Technology 118, 251–270. Ottino, J.M., Khakhar, D.V., 2000. Mixing and segregation of granular materials. Annual Review of Fluid Mechanics 32, 55–91. Rao, S.J., Bhatia, S.K., Khakhar, D.V., 1991. Axial transport of granular solids in rotating cylinders. Part. 2: Experiments in a non-flow system. Powder Technology 67, 153–162. Rothenburg, L., Bathurst, R.J., 1991. Numerical simulation of idealised granular assemblies with plane elliptical particles. Computers and Geotechnics 11, 315–329. Rutgers, R., 1965. Longitudinal mixing of granular materials flowing though a rotating cylinder. Chemical Engineering Science 20, 1097–1100. Schutyser, M.A.I., Padding, J.T., Weber, F.J., Briels, W.J., Rinzema, A., Boom, R., 2001. Discrete particle simulations predicting mixing behaviour of solid substrate in a rotating drum fermenter. Biotechnology and Bioengineering 75, 666–675. Schutyser, M.A.I., Weber, F.J., Briels, W.J., Boom, R.M., Rinzema, A., 2002. Three-dimensional simulation of grain mixing in three different rotating drum designs for solid-state fermentation. Biotechnology and Bioengineering 79, 284–294. Sheritt, R.G., Chaouki, J., Mehrotra, A.K., Behie, L.A., 2003. Axial dispersion in the three-dimensional mixing of particles in a rotating drum reactor. Chemical Engineering Science 58, 401–415. Shinbrot, T., Alexander, A., Muzzio, F.J., 1999. Spontaneous chaotic granular mixing. Nature 397, 675–678. Skinner, A.E., 1969. A note on the influence of interparticle friction on the shearing strength of a random assembly of spherical particles. Géotechnique 1, 150–157. Sudah, O.S., Chester, A.W., Kowalski, J.A., Beeckman, J.W., Muzzio, F.J., 2002. Quantitative characterization of mixing processes in rotary calciners. Powder Technology 126, 166–173. Yang, R.Y., Zou, R.P., Yu, A.B., 2003. Microdynamic analysis of particle flow in a horizontal rotating drum. Powder Technology 130, 138–146.