The influence of DEM simulation parameters on the particle behaviour in a V-mixer

Chemical Engineering Science 57 (2002) 3621 – 3638

www.elsevier.com/locate/ces

The in!uence of DEM simulation parameters on the particle behaviour in a V-mixer H. P. Kuoa , P. C. Knighta , D. J. Parkerb , Y. Tsujic , M. J. Adamsd , J. P. K. Sevillea; b; ∗ a Department

of Chemical Engineering, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK Imaging Centre, School of Physics and Astronomy, The University of Birmingham, Birmingham B15 2TT, UK c Department of Mechanical Engineering, Osaka University, Yamada-oka 2-1, Suita, Osaka 565-0871, Japan d Unilever Research, Port Sunlight Laboratory, Bebington, Wirral CH63 3JW, UK

b Positron

Received 22 August 2001; received in revised form 19 December 2001; accepted 30 January 2002

Abstract Results are described of simulations based on the discrete element method (DEM) using a code developed by Tsuji, Kawaguchi, and Tanaka (Discrete particles simulation of 2-dimensional !uidized bed. Powder Technology 77 (1993) 79–87). The mechanical interactions between particles and also between particles and the walls in granular !ows are modelled by linear springs, dash-pots and friction sliders. The simulation parameters are the restitution coe:cient, normal sti;ness, friction coe:cient between particles and between particles and the walls, and two ratios which relate the normal and tangential sti;ness and damping coe:cients. Their in!uence on particle motion in a V-mixer has been evaluated and compared with radioactive tracer measurements of particle motion. A number of quantitative methods for comparing DEM and experimental data were developed. Given the simpli>ed nature of the modelled interactions, the agreement between the predicted and measured data is remarkably close for restitution coe:cient values of 0.7 and 0.9, internal friction coe:cient values of 0.3 and 0.6 and wall friction coe:cient values of 0 and 0.3. The internal and wall friction coe:cients are important in determining the initiation of particle movement, while the value of the restitution coe:cient has a larger in!uence on particles in a dynamic state. The simulation of the fully elastic case (coe:cient of restitution = 1:0) with zero internal and wall friction, gives results that are very di;erent from the experiment data. ? 2002 Elsevier Science Ltd. All rights reserved. Keywords: Granular material; Mixing; Powder technology; Simulation; Solid mechanics; V-mixer

1. Introduction Simulations based on the discrete element method (DEM) have been used to study granular !ows in di;erent systems, including pneumatic conveying (Tsuji, Tanaka, & Ishida, 1992), !uidisation (Tsuji, Kawaguchi, & Tanaka, 1993), rotating drums (Yamane, Nakagawa, Altobelli, Tanaka, & Tsuji, 1998), and spouted beds (Kawaguchi, Sakamoto, Tanaka, & Tsuji, 2000). Soft-sphere simulations are commonly used where a >nite value of the contact overlap is allowed. This approach works well in slow-shearing sliding !ows of assemblies at high values of packing density. Hard-sphere simulations, in which there is instantaneous detachment of point contacts during particle collisions, ∗ Corresponding author. Tel.: +44-121-414-53-22; fax: +44-121-414-53-77. E-mail address: [email protected] (J. P. K. Seville).

are more appropriate for a collisional !ow regime where the assembly packing fraction is low. It has been shown by the current authors that the !ow patterns and other important characteristics associated with a widely used mixer known as a V-mixer, which were computed using the soft-sphere simulation code developed by Tsuji et al. (1993), are consistent with experimental data using glass ballotini (Kuo, Knight, Burbidge, Parker, Tsuji, & Seville, 2001). A V-mixer was chosen for this study because it is one of the simplest mixing devices and yet displays the division and recombination steps that are characteristic of solids mixing equipment. The current paper describes an extension of this work. It was aimed at elucidating the sensitivity of the simulation results to the selection of the numerical values of the parameters that are used in the code to represent the mechanical interactions of particles with other particles and at equipment walls. It is important to establish the sensitivity to the numerical values of the parameters used since simple spring and

0009-2509/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 0 9 - 2 5 0 9 ( 0 2 ) 0 0 0 8 6 - 6

3622

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

dash-pot interaction laws are used in the code. This has the advantage of computational e:ciency, but it is at the expense of a complete physical representation of the interactions compared with models based on contact mechanics theory. Consequently, one of the aims of the work was to investigate critically the implications of using such laws on the results of the simulations. The calculated results were compared with radioactive tracer experiments and a number of test procedures have been devised. 2. Interaction laws Cundall and Strack (1979) >rst proposed DEM for application to granular systems. It has been developed by a number of workers who have implemented a range of interaction laws. Some workers (e.g. Thornton & Yin, 1991; Ning & Thornton, 1993; Thornton, 1997) have employed interaction laws based on theoretical contact mechanics solutions, for example Hertzian contact theory (Johnson, 1985) for the elastic collision of two solid spheres and the Mindlin and Deresiewicz (1953) model for the tangential forces (e.g. Langston, TJuzJun, & Heyes, 1994, 1995). Other workers have employed computationally more e:cient simpler contact models as discussed below. When an inelastic normal impact occurs between two particles, the kinetic energy may be dissipated through internal friction, viscoelastic=plastic deformation and elastic waves. In the model employed in the current work, such an interaction is represented by a linear spring and dash-pot in parallel. This method of representing the dynamic response of a viscoelastic material is known as the Kelvin=Voigt model in mechanics. A friction slider, which is correlated with the frictional coe:cient of the particles and the wall, is used to model the magnitude of the frictional forces. Amontons’s law is applied such that the frictional force is proportional to the normal loading when gross slip occurs. The DEM simulation model developed by Tsuji et al. (1993), which was employed here, utilises two additional parameters: the normal=tangential sti;ness ratio and the normal=tangential damping coe:cient ratio. This model has been widely applied in the simulation of granular !ows (Tsuji et al., 1993; Yamane et al., 1998; Horio, Mikami, & Kamiya, 1998; Kawaguchi et al., 2000). Other authors have employed di;erent forms of contact model. Walton and Braun (1986) used two linear springs with di;erent values of sti;ness to model the normal impact and rebound between particles. Kuwabara and Kono (1987) applied a non-linear spring–damping model in order to simulate the interaction between two spheres. In the model, for a restitution coe:cient, e, close to unity, the value of (1−e) was proportional to (incident velocity)1=5 , which was consistent with experimental near-elastic impact results. Tsuji et al. (1992) employed a non-linear Hertzian spring and dash-pot contact model for which the damping term was a function of both the displacement and displacement rate. Realistic !ow

patterns were obtained for the plug !ow of cohesionless particles in a horizontal pipe. The damping term in this model was di;erent from that used by Kuwabara and Kono (1987). The linear spring and damping model may be written in the following form: mxJ + x˙ + Kx = 0;

(1)

where m is the mass of the particle, x is the overlap distance,  is a damping coe:cient, K is the sti;ness of the linear spring and the dot denotes the time di;erential. Such models and more complex systems of springs and dash-pots may provide an adequate response to an imposed strain for a viscoelastic material in the sense that the elastic modulus and the loss tangent representing the damping will show the expected strain rate sensitivity. Applying the initial conditions: x = 0; x˙ = V0 at t = 0, where V0 is the impact velocity, the solution of Eq. (1) is then
V0 x=
exp(− !t) sin( 1 − 2 !t); (2) ! 1 − 2
dx = V0 exp(− !t) cos( 1 − 2 !t) dt
V0 −
exp(− !t) sin( 1 − 2 !t); (3) 1 − 2 where  !=

K m

(4a)

and

 : (4b) 2m! The coe:cient of restitution is given by the following expression:      −   e = exp
(5) :  4mK − 2  =

It may be seen that the coe:cient of restitution is independent of the impact velocity, which arises because the vibrational frequency of the system is a constant controlled by the sti;ness of the spring. This is contrary to the expected behaviour of either viscoelastic or plastically deforming bodies where, in general, the coe:cient of restitution is a function of impact velocity (Thornton, 1997). Although it has been shown that the relative approach velocity (Labous, Rosato, & Dave, 1997) and impact angle have some e;ect (Gorham & Kharaz, 2000), under near-elastic conditions the coe:cient of restitution of a particle is approximately constant. Zhang and Whiten (1996) compared the linear spring– damping model developed by Tsuji et al. (1993) (i.e. Eq. (1)) and the non-linear spring–damping model by Tsuji et al. (1992) which may be represented by the following expression: mxJ + x1=4 x˙ + Kx3=2 = 0:

(6)

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

3623

Table 1 Relationship between coe:cient of restitution e and 

e 

1 0

0.9 0.03

0.8 0.07

0.7 0.11

0.6 0.16

0.5 0.22

In the non-linear model, because the damping force is a function of both displacement and the rate of displacement, the magnitude of the initial contact force is zero, unlike the linear model. The non-linear model gives a qualitatively acceptable form for the inter-particle force at initial contact. However, the force becomes attractive before the particles separate (x = 0), represented by a change in its direction; this is not physically acceptable behaviour since an adhesion function is not embedded in the model. Therefore, they suggested that the particles should separate at the time at which the force is zero, although this corresponds to a non-zero displacement. Zhang and Whiten (1996) normalised Eq. (1) as d 2 xˆ d xˆ + xˆ = 0; + 2 2 d tˆ d tˆ where xˆ x= ;   m ˆ t=t K and

(7)

(8a) (8b)

 = √ : (8c) 2 Km They plotted the normalised position, velocity and acceleration (force) as a function of the normalised time using di;erent values of  in Eq. (8c), ranging from 0 to ∼ 1:5. In the notation of the paper by Tsuji et al. (1993), the parameter  can be correlated with the coe:cient of restitution by  ; (9) = √ 1 + 2 where −ln e = : (10)

Thus, the parameter  can be calculated using the values of e given in Table 1. For the type of particles commonly simulated using DEM, values of e are ¿ 0:7, thus appropriate values for  are smaller than 0.1. The range of  used by Zhang and Whiten (1996) was, therefore, too large and included the over-damped condition. Fig. 1 shows the displacement, velocity and force (acceleration) as a function of elapsed time for di;erent values of coe:cient of restitution e. For e¡1:0, the inter-particle force has a non-zero value at initial contact and it becomes attractive before detachment occurs. For near-elastic collisions (e.g. e ¿ 0:8), these e;ects are small in magnitude and the linear model shows physically reasonable results. However, for smaller values of e (e.g. e¡0:8), these e;ects

0.4 0.28

0.3 0.36

0.2 0.46

0.1 0.59

0.001 0.91

become more signi>cant. This non-physical behaviour occurs because of the over-simpli>ed nature of the spring and dash-pot approximation. If alternative boundary conditions are considered as a zero initial force and a speci>ed incident velocity, the solutions for Eq. (7) are di;erent from those obtained by Zhang and Whiten (1996) and are given below:
xˆ = −2 exp(−tˆ) cos( 1 − 2 tˆ)
1 − 22 +
exp(−tˆ) sin( 1 − 2 tˆ); 1 − 2
d xˆ = exp(−tˆ) cos( 1 − 2 tˆ) d tˆ
 exp(−tˆ) sin( 1 − 2 tˆ); +
1 − 2 2
d xˆ 1 = −
exp(−tˆ) sin( 1 − 2 tˆ): 2 2 1− d tˆ The coe:cient of restitution is then     −  e = exp √ : 2 Km 

(11)

(12) (13)

(14)

Fig. 2 shows plots of the normalised displacement, velocity and force (acceleration) as a function of the normalised time for di;erent values of the coe:cient of restitution. At time zero, although the initial contact force decreases from zero, the displacement starts from a non-zero value, which is physically impossible. However, for near-elastic collisions, this spurious displacement is relatively small. Thus there is an analogous set of problems to those exhibited for the previous boundary conditions. In summary, the linear spring and damping model may be able to provide reasonable behaviour for nearly elastic materials, such as glass ballotini, but not necessarily for complex materials of practical interest. The time step and sti;ness used in the simulation are important parameters. Tsuji et al. (1993) suggested that the discrete time step, Rt, should be less than or equal to 1=10 of the oscillation time, which for the linear spring and damping system is given by 

m Rt 6 ; (15) 5 K where m is the mass of each sphere and K is the spring sti;ness. The magnitude of K cannot be unambiguously related to the elastic properties of a spherical particle; the sti;ness derived from the classic Hertzian contact theory is a non-linear function of the deformation (Johnson, 1985). Thus there is

3624

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

e=1.0

e=0.7

1

1

0.5

0.5

0

0 0

1

2

3

-0.5

0

1

2

3

-0.5

-1

-1 Normalised time

Normalised time

e=0.9

e=0.6

1

1

0.5

0.5

0

0 0

1

2

3

-0.5

0

1

2

3

-0.5

-1

-1 Normalised time

Normalised time

e=0.8

e=0.001

1

12

0.5

8.5

0

5 0

1

2

3

-0.5

1.5

-1

-2 Normalised time

0

1

2

3

Normalised time

Fig. 1. Normalised displacement (- - -), velocity (—) and force ( ) as a function of the normalised time for a range of restitution coe:cients calculated using Eq. (7) with boundary conditions: t = 0, displacement = 0, velocity = 1.

no physical basis for its selection. It has become normal practice to use values of K that are considered unrealistically small because the particle overlaps/deformations and contact times are greater than in real contacts. This is done for computational e:ciency. Yuu, Abe, Saitoh, and Umekage (1995) studied the particle behaviour in rectangular hoppers over several orders of magnitude of the time steps using a linear-spring and damping model. The time steps were 10−6 ; 10−5 and 10−4 s, which correspond to the values of sti;ness of 7×107 ; 7×105 and 7 × 103 N=m. It was found that the sti;ness values had relatively little e;ect on the particle behaviour even though

a range of several orders of magnitude was examined. Other examples of this approach include the work of Kawaguchi et al. (2000); Kuwagi, Mikami, and Horio (2000), and Kuo et al. (2001). Typical sti;ness values used in these calculations are of the order 103 N=m which gives a minimum time step to ca. 10 s. Horio et al. (1998) simulated particle motion inside a !uidised bed using di;erent values of the spring sti;ness and the restitution coe:cient for the linear springs, dash-pots and sliders; obvious di;erences between the results using di;erent parameters were not observed in the range examined. Cleary, Metcalfe, and Li;man (1998) studied solids

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

e=1.0

3625

e=0.7

1

1

0.5

0.5

0

0 0

1

2

0

3

1

2

3

-0.5

-0.5

-1

-1

Normalised time

Normalised time

e=0.9

e=0.6

1

1

0.5

0.5

0

0 0

1

2

3

-0.5

0

1

2

3

-0.5

-1

-1 Normalised time

Normalised time

e=0.8

1

e=0.001

1

0.5 0 0

1

2

3

-0.5 0

1

2

3

-0.5

-1

-2 Normalised time

Normalised time

Fig. 2. Normalised displacement (- - -), velocity (—) and force ( ) as a function of the normalised time for a range of restitution coe:cients calculated using Eq. (7) with boundary conditions: t = 0, velocity = 1, force = 0.

mixing and segregation in two-dimensional rotating drums using the linear springs, dash-pots and sliders model. In their simulations, di;erent particle–particle friction coe:cients, coe:cients of restitution, and size distributions were examined to assess the sensitivity of the quality of mixing to these parameters. The mixing rate was found to be sensitive to the coe:cient of restitution for values smaller than 0.5, to the magnitude of the frictional coe:cients and to size variations ¿ 20%. The rate of mixing increased when the restitution and friction coe:cients were decreased and the size variations increased. However, a detailed interpretation was not given. Hoomans, Kuipers, Briels, and van Swaaij (1996) used both measured and fully elastic frictionless collision parameters (e = 1 and p = w = 0) to study particle motion in a !uidised bed in two-dimensional hard-sphere simulations.

In the fully elastic case, the resulting simulations were different from physical reality: bubbling was completely absent and hence particle circulation motion was suppressed. Simulations with measured values corresponding to e = 0:9 and p = w = 0:3 showed realistic !ow behaviour for a powder which may be classi>ed as type D under the Geldart scheme (Geldart, 1973). Explosive bubble growth was observed in the simulations followed by a transition to slug !ow, which is typical for type-D !uidisation. The bubble growth, gas-phase hydrodynamics and slug formation in !uidised beds were sensitive to the values of the restitution and friction coe:cients. However, only two cases were included in their simulations. It should be emphasised that, because of the long computing times required, only limited studies of the in!uence of the interaction parameters have been carried out previously.

3626

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638 180°

0°/ 360°

Division 0° - 180°

Combination 180° - 360°

Fig. 4. Schematic representation of the blender position.

Fig. 3. Schematic drawing of the V-blender (coordinate Y is into the paper).

3. Simulations and experiments Particle motion in a V-mixer has been studied by only a few researchers (Sawahata, 1967; Inoue & Yamaguchi, 1969; Brone, Alexander, & Muzzio, 1998; Moakhar, Shinbrot, & Muzzio, 2000; Kuo et al., 2001) due to the di:culties in obtaining detailed information about the complex !ow behaviour. DEM simulations o;er a powerful technique for investigating such motion at a particulate level. The previous work done by the current authors (Kuo et al., 2001) compared the simulation results with tracer radioactive experiments using the positron emission particle tracking (PEPT) technique (Parker, Broadbent, Fowles, Hawkesworth, & McNeil, 1993). However, only one set of simulation parameters was used. In the current study, several sets of parameters were investigated. The geometry of the V-mixer used in the simulations and the experiments is given in Fig. 3. It was constructed from Perspex acrylic tubing and the particles utilised in the experiments were 3 mm glass spheres. In the simulations, uniform spherical particles with a diameter of 3 mm and a density of 2503 kg=m3 were used. The linear spring and dash-pot model with the boundary conditions shown in Fig. 1 was used to simulate particle contact. A >xed normal/tangential sti;ness ratio of 0.25 and a >xed normal/tangential damping coe:cient ratio of 0.5 were speci>ed. A >ll fraction

of 20% v/v was investigated, which corresponds to 11,168 particles, at four di;erent rotational speeds of 15, 30, 45 and 60 rpm. For these simulations, gravity was imposed as the only external force. The values of the interaction parameters investigated are summarised in Table 2. The maximum time step, Rt, can be obtained from Eq. (15). In the current study, it should correspond to ¡120 s. A value of 10 s was actually selected, i.e. an order of magnitude smaller than the maximum recommended time, and the total simulation time was ca. 8 s. The particle position and velocity information were stored at 0:1 s intervals. The computing facility used had two EV6 CPUs running at 550 MHz. 4. Results and discussion 4.1. The in;uence on the ;ow pattern in the arms Particle motion was computed relative to rotating Lagrangian coordinates which were >xed with respect to the body of the mixer. In the V-mixer, two distinct parts of the rotation are identi>ed: division and combination (Fig. 4). Figs. 5 and 6 show the 5 sub-steps of the combination and division steps, respectively, for simulations A, C and G. These cases are selected in order to visualise the in!uence of the coe:cient of restitution (A (e = 0:9) and C (e = 0:7)) and the friction coe:cient (A (p = w = 0:3) and G (p = w = 0)). The transition from the static

Table 2 Simulation parameters used in the V-mixer modelling ID

Restitution coe:cient

Friction coe:cient

Normal sti;ness (N/m)

A B C D E F G H

0.9 0.9 0.7 0.9 1.0 1.0 0.9 0.9

p p p p p p p p

1000 500 1000 1000 1000 1000 1000 1000

p : particle–particle friction coe:cient; w : particle–wall friction coe:cient.

= w = 0:3 = w = 0:3 = w = 0:3 = w = 0:6 = w = 0:0 = w = 0:3 = w = 0:0 = 0:3; w = 0:0

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

3627

Fig. 5. Particle motion during a combination step at 60 rpm (cases A, C, and G).

to the dynamic states for both steps is not a;ected by the coe:cient of restitution while there is a considerable sensitivity to the friction coe:cient. When particles are travelling freely down the tubes, particles with greater values of the restitution coe:cient (A) transport faster than those with a smaller value (C). For cases A and B with di;erent sti;ness values, no di;erence could be distinguished. The volume inside the mixer was discretised using cubic cells with a side of length 7:35 mm. The time average velocity for each cell was calculated by averaging the velocities within the cell for all particles visiting it. The calculated velocity components in the X –Z plane for cases A–H together with the values measured using the PEPT technique are shown in Fig. 7. The magnitude and the directions of the velocity vectors in the !ow pattern charts are similar when using di;erent parameters, except for the fully elastic and frictionless case E. In the last case, the cir-

culation !ow pattern is weak. The major di;erence between the other cases is that the velocities near the walls are more pronounced for cases G and H, for which the particle–wall friction coe:cient is zero. This result is as might be reasonably expected. The data shown in Fig. 7 show circulating !ow patterns roughly corresponding to the axial plane of symmetry, although some velocity vectors are located outside the physical boundary of the mixer. Apparent locations outside the physical boundaries result from a combination of a small !uctuation of the rotational speed of the vessel and the limitations of the spatial resolution of the experimental technique. However, these errors are not signi>cant. In the following comparisons, the locations of the particle recorded outside the physical boundaries of the mixer are ignored. In order to make a more quantitative comparison between simulation results and experiments, values of axial and

3628

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

Fig. 6. Particle motion during a division step at 60 rpm (cases A, C, and G).

radial velocity were computed at the cross-sectional plane whose position is shown in Fig. 7(A). This plane is located at half of the length of the arm, subtracting the join part of the mixer, and is sliced into 9 sections as indicated in Fig. 7(A). The radial and axial directions are de>ned as the X and Z directions shown in Fig. 3. Fig. 8 shows the relationship computed from the simulations and compares them with those obtained from PEPT. Considering the axial velocities >rst, agreement between the computed and measured velocities is satisfactory, except for case E (fully elastic and frictionless) and cases G and H (zero wall friction coe:cient). The radial velocity data show greater divergence between simulation and experiment. The computed radial velocity is sensitive to the magnitude of the wall friction coe:cient and the best >t is obtained with case G (zero wall friction coe:cient).

4.2. The in;uence on the speed distribution and the velocity ;uctuation in the arms The particle average speeds calculated from the three velocity components in each cell are shown in Fig. 9. Cases A–D, and F all show qualitative agreement with the PEPT result. The results for cases E, G and H, which all have a zero wall friction coe:cient, are di;erent from the other cases and the experimental observations. In these cases, as might be expected from the absence of a frictional resistance at the walls, the particles move faster in this region. In case D, the size of the high-speed region near the centre is slightly smaller than in the other cases. This is attributable to the local high collision frequency and high friction coe:cient. The most marked e;ect occurs in the fully elastic case E, in which the high-speed

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

3629

Fig. 7. Average velocity distributions for cases A–H together with the experimental data at 60 rpm. The length of the arrows is proportional to the magnitude of the velocity.

regions correspond to the wall boundaries of the vessel; this is quite di;erent from the measured speed distribution (note: a di;erent grey scale is used in case E to display the data).

Fig. 10 shows the values of the speed at the cross-sectional plane shown in Fig. 7(A) with the radial positions numbered from 1 to 9. Cases A and B give a better prediction of the speed at this plane compared to other cases. The speeds

3630

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

be seen that the speed is reduced signi>cantly at low radial positions. The speed distribution allows a more quantitative comparison with PEPT than is possible from the vectors in Fig. 7 since only the magnitude is involved. The di;erences between the calculated and measured speed values for each cell were calculated and summarised. The rank of the agreement in the speed comparison between DEM and PEPT is in the following order: B ≈ A ≈ F ¿ C ≈ D  E  H ≈ G. The velocity !uctuation inside each cell was compared for the di;erent cases. The comparison was based on the average of a sum of the squares of the three !uctuating velocity components:

Radial velocity, mm/s

10

-10

-30

-50

n

T= -70 0

2

4

6

8

10

Radial position

(a)

A

B

C

D G

E H

F PEPT

Axial velocity, mm/s

200

0

1 * (vi − v) V 2; 3

(16)

i=1

where T is the velocity !uctuation inside each cell * (m2 =s2 ); vi is the instantaneous particle velocity (m/s) for the ith data point, vV is the average velocity inside each cell (m/s), and n is the number of data points and each one a di;erent “pass” by the tracer in the experiments or the particle located inside each cell. Fig. 11 shows the calculated velocity !uctuation distributions for the di;erent cases. Except for case E, the computed distributions and the measured distributions are quite similar, with the greatest values near the centres of the arms. Cases G (p = 0) and H (p = 0:3) have di;erent particle–particle friction coef>cients. Case H, with the physically more realistic particle friction coe:cient value of 0.3, gives a better agreement with the experimental result. 4.3. The in;uence on the exchange between the two arms

-200

-400 0

2

4

6

8

10

Radial position

(b)

A D G

B E H

C F PEPT

Fig. 8. The velocities as a function of the radial position at the cross-sectional plane shown in Fig. 7(A): (a) radial velocity comparison and (b) axial velocity comparison.

predicted by cases G and H, which all have a zero wall friction coe:cient, are higher than the PEPT measurement close to the wall (positions 1, 8 and 9). Case D with a particle– wall friction coe:cient value of 0.6 is di;erent from the PEPT; however, the disagreement is not as large as in cases G and H with zero wall friction coe:cients. Case C has the lowest value of the coe:cient of restitution and it can

The data were analysed to determine the exchange rate, de>ned as the probability of each single particle crossing to the other arm of the V-mixer per rotation. Data for rotational speeds of 15, 30, 45 and 60 rpm were analysed in this way and the results for each set of interaction parameters are shown in Fig. 12. The exchange rate provides a quantitative test of the ability of DEM to predict mixing behaviour. The results reveal that the mixing is sensitive to the parameters used in the simulation. The exchange rate for the fully elastic simulation E is about one order of magnitude greater than the other cases. The calculated results are similar to the experimental values except for G (no inter-particle or wall friction) which is signi>cantly greater and H (no wall friction) which is much smaller. The large di;erence between the results for A and H, which di;er only in wall friction, suggests that the exchange rate is particularly sensitive to this parameter rather than the restitution coe:cient. With a zero wall friction value, the simulation results show the smallest exchange rates due to the particles falling more rapidly to the two ends of the arms during the division step and reaching the join between the arms faster during the combination step. Intuitively, the exchange rate should decrease with the time the particles spend in the vicinity of the plane

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

3631

Fig. 9. Speed distributions for cases A–H together with the experimental data at 60 rpm. The grey scale of case E is di;erent from the other cases. Units are mm/s.

3632

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

increase in the exchange rate during the combination step. These calculated results suggest that increasing the particle– wall friction may increase the exchange rate at the division step and a small wall friction promotes an increase in the exchange rate at the combination step.

400

Speed, mm/s

300

4.4. The in;uence on the circulation time in the arms

200

100

0 0

2

4

6

8

10

Radial position

A D

B E

G

H

C F PEPT

Fig. 10. The speed as a function of the radial position at the cross-sectional plane shown in Fig. 7(a).

of symmetry, in which they may be dispersed to the other arm. Thus, with a zero particle–wall friction coe:cient, the exchange should decrease at the division step, because particles will spend less time in the exchange region, and increase at the combination step. It has been shown elsewhere by the current authors that the exchange between the two arms of a V-mixer is more likely to occur in the division step rather than the combination step (Kuo et al., 2001). The number of crossings of the plane of symmetry during each half revolution was distinguished in the DEM calculations and is shown in Fig. 13. The data shown in the >gure represent the number of crossings in four revolutions (four times each in the division and combination steps). Fig. 13(a) shows the in!uence of the particle–wall friction (A: w = 0:3; D: w = 0:6, and H: w = 0). Case H shows a signi>cant reduction in the exchange rate compared to cases A and D at the division step. The magnitude of the reduction is more signi>cant than the increase in the rate at the combination step, and this accounts for the smallest exchange rate for all the cases observed in Fig. 12. Also in Fig. 13(a) for case D, with the greatest particle–wall friction coe:cient, the exchange rate increases during the division step which results in a higher overall exchange rate. Fig. 13(b) compares the in!uence of the particle–particle friction and the particle–wall friction (A: p = w = 0:3 and G: p =w =0). In both cases, the division step is responsible for the larger contribution to the overall exchange rate and the removal of internal friction and external friction has a negligible e;ect. However, it does result in a substantial

The circulation times in the two arms may be compared by recording the times when particles cross a pair of imaging planes, shown as the two levels in Fig. 3. The circulation time is de>ned as the time required for particles to travel from the >rst plane to the second plane, and return to the >rst plane. Although there are inevitable experimental errors in the location obtained from PEPT, the time corresponding to that location is extremely accurate. Consequently, comparisons of the circulation time can be made with more con>dence than direct comparisons of velocities derived from locations. The circulation times are shown as a function of the vessel rotation period in Fig. 14 for the eight cases. The major di;erence is observed for case E, which does not exhibit a monotonic increase in the time with increasing vessel. All the other combinations give results that are of the same order of magnitude as the experimental values; however, none of the DEM simulation cases can give quantitative agreement in the circulation time. The particles in the experiment take longer to circulate in the mixer than the DEM predictions, especially at short vessel rotation periods. 4.5. The in;uence on the dispersion in the arms Martin (1999) proposed a method for assessing dispersion within a system. It is based on analysing all particle trajectories that have an origin in a common cell. It involves a calculation of the current position as a function of the time since the particle visited the cell of interest. A dispersion variance provides a measure of the extent of local mixing. A comparison of the computed dispersion of particles with the experimental values has been made previously (Kuo et al., 2001). For the PEPT experiments, the dispersion analysis must be based on trajectories that are independent of time. However, for a tumbling mixer, the motion is repeated every cycle, but there must be a separate analysis of the half-cycles because of the !ow reversal. The particle motion is independent of time approximately every half revolution and the dispersion in the mixer was therefore considered in the combination and division steps separately. Fig. 15 shows the results of the dispersion analysis during the combination step at 60 rpm. The 3D cubic grid was that de>ned above and the time duration for the particles to disperse from the starting cell was set at 200 ms. During this period of time, the particles can travel ca. 20 mm, which corresponds to a third of the diameter of the tube. Using this limit ensures that the local dispersion value is a;ected only by the

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

3633

Fig. 11. Velocity !uctuation distributions for cases A–H together with the experimental data at 60 rpm. The grey scales are di;erent to characterise each plot. Units are m2 =s2 .

3634

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638 0.70

0.14

B D G PEPT

0.60

0.10

0.50

0.06

0.40

600

Number of crossings

A C F H E

Exchange rate for case E, rev −1

Exchange rate, rev −1

0.18

350

0.30

0.02 0

15

30

45

60

100

75

0

Rotational speed, rpm

Fig. 12. The in!uence of simulation parameters on the exchange rate for cases A–H together with the experimental data at 60 rpm.

5

4

2 3 Times

A-combination

D-combination

A-division

D-division

H-combination

H-division

Number of crossings

600

350

100 0

1

2

(b)

3 Times

4

5

A-combination

G-combination

A-division

G-division

Fig. 13. The in!uence of the particle–particle friction and the particle–wall friction on the exchange rate (60 rpm): (a) the in!uence of the particle–wall friction (A: w = 0:3; D: w = 0:6; H: w = 0:0); (b) the in!uence of the particle–particle friction and the particle–wall friction (A: p = w = 0:3; G: p = w = 0:0).

4.0 Circulation time, sec

local conditions. During the combination step, particles travel away from the two ends of the arms and move towards the centre plane. In cases A and B, the change in the normal sti;ness did not signi>cantly in!uence the results. For case C, there is an enhancement in the dispersion when particles travel along the tube. This implies that with a low restitution coe:cient, particles exhibit greater dispersion within the arms. Case D corresponds to particles with high particle–particle and particle–wall friction coe:cients and the dispersion is greatest when they move from a static state. This implies that friction may increase the degree of freedom of motion locally thus improving the local quality of mixing. Case E, with the fully elastic collision parameters, shows unrealistic results, in that the particles are dispersed to an almost uniform extent. Case F (e = 1:0) shows a similar distribution to cases A (e = 0:9) and C (e = 0:7), which indicates the relative unimportance of the coe:cient of restitution when the particles are close packed. However, case G (p = w = 0) is di;erent from cases A (p = w = 0:3) and D (p = w = 0:6). With zero friction coe:cients, the particles travel faster from the two ends of the arms to the central region and the better dispersed areas shift nearer to the symmetrical plane. Fig. 16 shows the dispersion analysis at the division step for a dispersion time of 200 ms in all the cases. Particles are moving from the join toward the ends of the two arms. The changes in sti;ness in cases A and B again do not show an e;ect. In case C, the well-dispersed region near the central plane breaks into two small sections which is similar to the experimental results. In case D, with the greater particle–particle and particle–wall friction coe:cients, a larger well-dispersed region is found near the join, where the particles are accelerating from a static to a dynamic state. In case E, with the fully elastic collision parameters, there is again the unrealistic result that the particles are well dispersed inside the vessel. The zero particle–wall friction associated with cases G and H results in the particles

1

(a)

3.0 2.0 A C E G PEPT

1.0 0.0 0.0

1.0

2.0

3.0

4.0

B D F H

5.0

Vessel rotation period, sec

Fig. 14. The in!uence of simulation parameters on circulation time.

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

Fig. 15. Dispersion variance for cases A–H together with the experimental data at 60 rpm during the combination step. Units are mm2 .

3635

3636

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

Fig. 16. Dispersion variance for cases A–H together with the experimental data at 60 rpm during the division step. Units are mm2 .

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

being able to travel further with a greater dispersion at the ends of the arms. These comparisons in both the division and combination steps show that the restitution coe:cient is the most important factor controlling the dynamics of the particles, while the friction coe:cient is a critical factor in controlling the transition from the static to the dynamic states. 5. Conclusions The in!uence of the DEM simulation parameters on particle behaviour in a V-mixer has been studied. Linear springs, damping dash-pots, and frictional sliders were used to model the normal and tangential contact between particles and particles and walls. A number of tests have been made to evaluate the in!uence of the parameters associated with the interaction laws. A comparison of the e;ects of using di;erent spring sti;nesses, restitution and friction coe:cients was made based on the velocity, speed and velocity !uctuation distributions, the exchange rate between the two arms, the circulation time in the mixer, and dispersion at division and combination steps. Quantitative comparisons between DEM and PEPT in the velocity and speed distributions using physically reasonable values of the particle–particle friction coe:cient 0.3, the particle–wall friction coe:cient 0:3 and the coe:cient of restitution 0:9 were made. We conclude that physically sensible values of the coe:cient of restitution, the particle–particle friction coe:cient and the particle–wall friction coe:cient give semi-quantitative agreement with experiments in the operating conditions studied. Quantitative di;erences were generally observed between DEM and PEPT when using a particle–wall friction coef>cient value of zero and a restitution coe:cient of unity. Changes in the friction coe:cient (0:3–0:6) and restitution coe:cient (0:7–0:9) simulation parameters give di;erent predictions. The friction coe:cient a;ects the particle transition from the static to the dynamic state more signi>cantly and the restitution coe:cient values have a greater in!uence on the behaviour of the moving particles. Future work should be made to investigate the e;ect of the calculation time step and the sti;ness values. Acknowledgement The authors are grateful for >nancial support from Unilever Research, the British Council (for supporting travel between Birmingham and Osaka) and the University of Birmingham (for a Bridon Postgraduate Award to H. P. Kuo). References Brone, D., Alexander, A., & Muzzio, F. J. (1998). Quantitative characterization of mixing of dry powders in V-blenders. A.I.Ch.E. Journal, 44, 271–278.

3637

Cleary, P. W., Metcalfe, G., & Li;man, K. (1998). How well do discrete element granular !ow models capture the essentials of mixing processes? Applied Mathematical Modelling, 22, 995–1008. Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. G>eotechnique, 1, 47–65. Geldart, D. (1973). Types of gas !uidisation. Powder Technology, 7, 285–292. Gorham, D. A., & Kharaz, A. H. (2000). The measurement of particle rebound characteristics. Powder Technology, 112, 193–202. Hoomans, B. P. B., Kuipers, J. A. M., Briels, W. J., & van Swaaij, W. P. M. (1996). Discrete particle simulation of bubble and slug formation in a 2-dimensional gas-!uidised bed: A hard-sphere approach. Chemical Engineering Science, 51, 99–118. Horio, M., Mikami, T., & Kamiya, H. (1998). Discrete element simulation of cohesive particle !uidization. In: World Congress of Particle Technology, 3, Brighton, UK. Inoue, I., & Yamaguchi, K. (1969). Particle motion in mixer and mixing process—mixing in a two-dimensional V-type mixer. Kagaku Kogaku, 33, 286–292 (in Japanese). Johnson, K. L. (1985). Contact mechanics. Cambridge: Cambridge University Press. Kawaguchi, T., Sakamoto, M., Tanaka, T., & Tsuji, Y. (2000). Quasi-three-dimensional numerical simulation of spouted beds in cylinder. Powder Technology, 109, 3–12. Kuo, H. P., Knight, P. C., Burbidge, A. S., Parker, D. J., Tsuji, Y., & Seville, J. P. K. (2001). Particle motion in a V-mixer: Experimental validation of the DEM simulations. A.I.Ch.E. Journal, submitted for publication. Kuwabara, G., & Kono, K. (1987). Restitution coe:cient in a collision between two spheres. Japanese Journal of Applied Physics, 26(8), 1230–1233. Kuwagi, K., Mikami, T., & Horio, M. (2000). Numerical simulation of metallic solid bridging particles in a !uidized bed at high temperature. Powder Technology, 109, 27–40. Labous, L., Rosato, A. D., & Dave, R. N. (1997). Measurements of collisional properties of spheres using high-speed video analysis. Physical Review E, 56(6), 5717–5725. Langston, P. A., TJuzJun, U., & Heyes, D. M. (1994). Continuous potential discrete particle simulations of stress and velocity >elds in hoppers: Transition from !uid to granular !ow. Chemical Engineering Science, 49, 1259–1275. Langston, P. A., TJuzJun, U., & Heyes, D. M. (1995). Discrete element simulation of granular !ow in 2D and 3D hoppers: Dependence of discharge rate and wall stress on particle interactions. Chemical Engineering Science, 50, 967–987. Martin, T. W. (1999). Studies of particle motion in mixers. Ph.D. Thesis, The University of Birmingham, UK. Mindlin, R. D., & Deresiewicz, H. (1953). Elastic spheres in contact under varying oblique forces. Transactions of the ASME Journal of Applied Mechanics, 20, 327–344. Moakhar, M., Shinbrot, T., & Muzzio, F. J. (2000). Experimentally validated computations of !ow, mixing and segregation of non-cohesive grains in 3D tumbling blenders. Powder Technology, 109, 58–71. Ning, Z., & Thornton, C. (1993). Elastic-plastic impact of >ne particles with a surface. In C. Thornton (Ed.), Powders and grains, 93, Rotterdam (pp. 33–38). Parker, D. J., Broadbent, C. J., Fowles, P., Hawkesworth, M. R., & McNeil, P. (1993). Positron emission particle tracking—a technique for studying !ow within engineering equipment. Nuclear Instruments and Methods, A 326, 592–607. Sawahata, Y. (1967). On the process of mixing and circulation of particles in a V-type mixer. Kagaku Kogaku, 31, 1212–1218 (in Japanese). Thornton, C. (1997). Coe:cient of restitution for collinear collisions of elastic-perfectly plastic spheres. Transactions of the ASME Journal of Applied Mechanics, 64, 383–386. Thornton, C., & Yin, K. K. (1991). Impact of elastic spheres with and without adhesion. Powder Technology, 65, 153–166.

3638

H. P. Kuo et al. / Chemical Engineering Science 57 (2002) 3621–3638

Tsuji, Y., Kawaguchi, T., & Tanaka, T. (1993). Discrete particles simulation of 2-dimensional !uidized bed. Powder Technology, 77, 79–87. Tsuji, Y., Tanaka, T., & Ishida, T. (1992). Lagrangian numerical simulation of plug !ow of cohesionless particles in a horizontal pipe. Powder Technology, 71, 239–250. Walton, O. R., & Braun, R. L. (1986). Viscosity, granular-temperature, and stress calculations for shearing assemblies of inelastic, frictional disks. Journal of Rheology, 30(6), 949–980. Yamane, K., Nakagawa, M., Altobelli, S. A., Tanaka, T., & Tsuji, Y. (1998). Steady particulate !ows in a horizontal rotating cylinder. The Physics of Fluids, 10, 1419–1427.

Yuu, S., Abe, T., Saitoh, T., & Umekage, T. (1995). Three-dimensional numerical simulation of the motion of particles discharging from a rectangular hopper using discrete element method and comparison with experimental data (e;ect of time steps and material properties). Advanced Powder Technology, 6, 259–269. Zhang, D., & Whiten, W. J. (1996). The calculation of contact forces between particles using spring and damping models. Powder Technology, 88, 59–64.