Modelling the solids inflow and solids conveying of single-screw extruders using the discrete element method

Powder Technology 153 (2005) 95 – 107 www.elsevier.com/locate/powtec

Modelling the solids inflow and solids conveying of single-screw extruders using the discrete element method P.A. Moysey, M.R. Thompson* Department of Chemical Engineering, McMaster University, Hamilton, Ontario, Canada L8S 4L7 Received 30 June 2004; received in revised form 9 November 2004; accepted 1 March 2005 Available online 20 April 2005

Abstract A new solids-conveying model for the single-screw extruder based on the Discrete Element Method (DEM) is proposed in this work. The polymer solids are treated as spherical particles moving in a 3-D environment which includes the feed hopper, the solids-inflow zone, and the solids-conveying region of an extruder, without inclusion of the plug flow assumption common to continuum models. Normal and tangential forces resulting from inelastic collisions with neighboring particles and surfaces dictate how each polymer pellet is conveyed through the model extruder. The DEM technique was implemented in this work to allow fundamental study of the local transport phenomena within the screw channel. Reported in this paper are results examining the cross- and down-channel velocity profile of solids in the screw; the residence time distribution; the cross-channel temperature profile; and the coordination number distribution. Two exit conditions were evaluated by the model: i) the open-discharge case where no compaction of the solids occurred; and ii) the restricted case where the axial pressure increased as the solids flowed towards the barrel exit. The predictions of the DEM simulations allowed for detailed observations of the solids movement in the screw, providing insight into the inherent flow fluctuations of extrusion systems. D 2005 Elsevier B.V. All rights reserved. Keywords: DEM; Simulation; Three-dimensional; Extrusion; Solids mixing

1. Introduction Bulk transport of solids within screw conveyors and feeders is based on the Archimedean principle of screw conveying, with these devices representing the most common process operations in solids handling plants due to their reliability and continuous flow characteristics [1]. Conveyors are usually operated starved at approximately 45% of capacity, whereas screw feeders which are often located at the bottom of storage bins and silos, are flood-fed yielding the maximum feed rate [2]. A more complex screw conveying device is the plasticating single-screw extruder. In the plastics industry, single-screw extruders are usually flood-fed solid polymer pellets which eventually melt due to viscous dissipation and conduction from the heated barrel. Unlike screw conveyors, the screw channel of an extruder is * Corresponding author. Tel.: +1 905 525 9140x23213; fax: +1 905 521 1350. E-mail address: [email protected] (M.R. Thompson). 0032-5910/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.powtec.2005.03.001

comparatively shallow in depth and includes a flight that is relatively thick in comparison. The restriction at the die end of an extruder generates high pressures (in the order of 10– 30 MPa), leading to compaction of the polymer pellets within the solids conveying zone [3]. It is this pressure generation within the solid assembly which principally differentiates single-screw extruders from screw conveyors or feeders. The majority of models found in the literature [3 –6] for solids conveying within single-screw extruders make the plug flow assumption in their derivations. This simplifies the analysis, allowing the output rate and pressure development to be determined from force and torque balances on an elastic plug. The result is a useful tool giving reasonably good agreement with output rate data provided the appropriate coefficient of friction is selected; however, the predicted pressure is greater than experimental data, often by several orders of magnitude [7]. Models based on this approach lack the ability to describe local solids motion within the screw channel and the dynamic nature of the flow which is observed in practice [8]. There are a few notable

96

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

exceptions to models using the plug flow assumption [7,9], with reported improvements in predicted output rate and pressure in comparison to experimental results. All of these models, analytical or numerical, are based on continuum mechanics, and therefore, are constrained to behave according a priori knowledge of solids motion in the system. The advantage of using a discrete approach, instead, is that the motion of each particle in a granular assembly is modeled separately, eliminating the need for the continuum, steady state, or plug flow assumptions. A new approach to modeling solids transport within an extruder is proposed in this paper, using the discrete element method to develop a non-isothermal, three-dimensional simulation of the solids conveying zone. The feed hopper and inflow zone are included in the simulation, eliminating the need to specify the pressure at the beginning of the conveying zone. The solids assembly of polymer pellets is modeled as spherical particles experiencing compaction without permanent deformation. The focus of this paper will be a description of the extruder model, and discussion of localized phenomena within the screw channel. Comparison of this DEM solids-conveying model with experimental data has been shown to give good agreement and is presented elsewhere [10].

2. The discrete element method The Discrete (or Distinct) Element Method (DEM) considers each particle in a granular assembly as a separate entity, which can interact with other particles or surfaces through collisions or lasting contacts. The method was first proposed by Cundall and Strack [11] to model particles of soil under shear, and has more recently been applied in the simulations of a wide array of applications, such as tumblers [12,13], centrifugal mills [13], and the solids-inflow zone of a single screw extruder [14]. Collision detection and determination of the contact forces is a major undertaking in any DEM model, as the simulation must keep track of the location and velocities of other particles and surfaces in reference to the particle of interest. As multiple contacts are possible with neighboring particles at any instance in time, the sum of all forces acting on a particle must be determined for each time step. The resultant acceleration due to the sum of all forces (contact and body) acting on a particle is calculated from Newton’s Second Law of Motion. Integration of the acceleration allows the velocity and position for the next time step to be calculated. By this approach it is possible to model the dynamic behavior of a granular flow in a complex 3-D geometry. Details of the DEM model used in this work are given in the following sections.

utilizing a variable latching spring model for the normal forces, and a non-linear incremental slipping model for the tangential (frictional) forces. For this discussion, Fi_ refers to the current particle, while Fj_ denotes a neighboring particle. 2.1.1. Normal force The force of a collision is represented by the degree of particle overlap (a) that occurs rather than considering actual deformation of the solids. This methodology is known as the ‘‘soft particle’’ approach. The magnitude of the normal force ( Fn) is given as: Fn ¼ K1 a

ða increasingÞ

ð1aÞ

Fn ¼ K2 ða  a0 Þ

ða decreasingÞ

ð1bÞ

where K 1 and K 2 are stiffness constants [15], and a 0, is the value of the overlap at zero normal force. The secondary stiffness constant, K 2, is taken to be constant in the model (although it can be a function of the impact velocity [16]), which gives a constant coefficient of restitution (e): rffiffiffiffiffiffi K1 e¼ : ð2Þ K2 Due to the difference in the calculated normal force for loading and unloading events, this method is known as the hysteresis model. 2.1.2. Tangential forces While the normal force only incurs translational acceleration, the tangential component of a contact force results in both translational and rotational acceleration. The magnitude of the tangential (or shear) force is based on the amount of tangential slip over the last time step and the history of the contact, and is limited by the coefficient of friction, l, to lF n. In the Walton –Braun model [16], the shear forces build up non-linearly, until the Coulomb friction limit is reached, at which point the incremental frictional force becomes zero, and full sliding occurs. The tangential slip, DS, is decomposed into components that are parallel and normal to the previous shear force, DS ‹ and DS –, respectively. The incremental slip perpendicular to the shear force from the previous time step (projected onto the current tangential plane), is assumed to have no previous surface strain [15], so the current shear force, F ts, is given by: Fst ¼ Fst1 þ ks DS‹ þ kso DS8

ð3Þ

where k s is the effective tangential stiffness and k so is the initial tangential stiffness. 2.2. Particle – wall interactions

2.1. Particle –particle interactions The particle – particle interactions in our model were first proposed by Walton [15] and Walton and Braun [16],

Within the DEM model, the walls making up the three dimensional geometry of our extruder are represented as two-dimensional surfaces. In order to determine whether a

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

contact occurs, the minimum distance between the surface and the pellet center must be calculated. If this distance is less than the pellet radius, then the contact force of the collision may be determined. The difficulty in determining this minimum distance is based on the complexity of the system being simulated. In this paper, at least twenty surfaces and ten edges are modeled to describe the geometry of the extruder. The surfaces make up the barrel, hopper, feed throat, screw core, and helical flight. Fig. 1 shows a graphical representation of the extruder geometry. Edges that exist when two surfaces intersect, present the greatest difficulty in modeling without generating erroneous contact detection, and will be discussed at greater length in the section pertaining to the screw flight. Particle –wall contact forces are modeled in a similar fashion to the particle – particle forces described in the previous section, being based on the approach of Walton [15,16]. Contact detection at the helical flight requires the most computational operations in the model, due to the nature of the helix. For this reason, detection with the helix face, the helix edge, and the top of the flight is discussed in detail in this paper. 2.2.1. Contact detection and force resolution with the helix face In order to describe the flight structure clearly, certain nomenclature must be adopted. The screw flight can be described as a helical ribbon built up from the screw root which makes up two retaining walls for the screw channel; the barrel and screw root make up the other two walls retaining surfaces for the rectangular channel cross-section. The equation of a helix, in parametric form, is given by: rðhÞ ¼ RcosðhÞi þ RsinðhÞj þ Chk

ð4Þ

where is R the radius of the helix, C is a parameter dependent on the pitch (C = Ls / 2p), and h is a parametric variable. The conveying force on the solids is provided by the advancing wall of the flight which is referred to in this work as the advancing helix face. The flight wall on the other side of the screw channel is referred to as the retreating helix face. Fig. 2 shows the relevant distances,

Torpedo (for restricted case)

Hopper Feed Ø Opening

e

Hf y x z

Barrel Surface

Screw Root

Db Ls

Fig. 1. Schematic diagram of the feed hopper and screw geometry.

97

Fig. 2. Schematic diagram showing some relevant distances to the determination of possible contacts with a helical flight.

denoted by d 1, d 2, and d 3, between a pellet (exaggerated size) and the advancing helix face of a flight. The distance d 1 is therefore: d1 ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  2ffi xi  x0  Rxy cosðh1  hN Þ þ yi  y0  Rxy sinðh1  hN Þ

ð5Þ where h N = 2pNt sim represents the cumulative rotation of the screw, N is the screw speed, tsim is the current simulation time, and h 1 = z i/C is determined from the current z-coordinate of the pellet. The coordinates x 0, y 0, and z 0 indicate the origin of the system which is located on the screw center-line, in the center of the feed opening. The distance d 2 is determined by the magnitude of d 1 and R xy:   d1 1 d2 ¼ 2C I sin : ð6Þ 2Rxy A contact with the helix occurs if d 3 is less than the pellet radius at the current R xy. The helix angle, /, is given by:   Ls / ¼ tan1 : ð7Þ 2pRxy Finally, the distance d 3 is dependent on d 2 and the helix angle, and from the viewpoint of the advancing flight may be given by: d3 ¼ d2 cosð/Þ:

ð8Þ

The distance d 3 must be less than the pellet radius, R p, for a contact with the advancing flight to occur. Eq. (8) is similar to Eq. (7) in [17], although the distance d 2 is calculated in a different manner in our work. If this condition is satisfied, the unit vector, nh, may be determined from the center of the pellet to the surface of the flight, which is equivalent to the negative of the binormal vector, B, of a helix. The parametric variable, h 3, at the contact point is:    d1  2 h3 ¼ h1 þ 2sin1 ð9Þ sin /  1 : 2Rxy The binormal vector is defined as T  N, where T and N are the tangential and normal vectors respectively, of the

98

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

contact point on the helix. The resulting equation for the binormal vector is: B ¼ ssinðh3 Þi  scosðh3 Þj þ jk ¼  nx i  ny j  nz k ¼  n h

ð10Þ

R C j ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; s ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 2 2 2 Rxy þ C Rxy þ C 2 The normal force imparted to the pellet by the flight ( F h) is calculated by the approach developed by Walton [15], and is based on the amount of overlap, a = R p  d 3, between the pellet and helix face during a collision (Eqs. (1a), (1b)). The force magnitude, F n, is projected onto the binormal vector: ð11Þ

The equations for contact detection are similar for the retreating flight, except that a parametric variable, h e, which represents the thickness of the flight in the radial direction, has to be taken into account. The variable h e is given by:   e 1 he ¼ 2sin ð12Þ sinð/b ÞDb where /b ¼ tan1



Ls pDb

 ð13Þ

The distance from the center of the sphere to the retreating flight at the current z coordinate is therefore: d1 ¼

np  nr ¼ npx sinðh1  hN þ he Þ  npy cosðh1  hN þ he Þ ð17Þ

where j is the curvature, and s is the torsion of the helix:

Fh ¼  Fn nh :

while the cross product with the retreating flight is:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 xi  x0  Rxy cosðh1  hN þ he Þ þ yi  y0  Rxy sinðh1  hN þ he Þ

ð14Þ

Eqs. (6) –(8) apply to the retreating flight, although the parametric variable at the contact point is now defined as:    d1  ð15Þ h3 ¼ h1 þ 2sin1 1  sin2 / 2Rxy

where np is the unit vector pointing from the center of the screw to the pellet center, and na and nr are the unit vectors from the screw center to the advancing and retreating helix faces at the current z-coordinate, respectively. A pellet can come into contact with the surface at the top of the flight only if np  na < 0 and np  nr > 0. Furthermore, a contact can only occur if R xy  R b < R p, where R b is the outer radius of the screw. Besides possible contacts with the helix faces and the top of the flight, the edges where the two helix faces and the top of the flight intersect must also be considered, which represents a complicated problem in the model. These edges are referred to as helix edges according to Fig. 2. Only the equations for determining a contact with the advancing helix edge will be discussed here, as they are both similar. The distance from the helix edge at the current z-position to the center of the sphere, d 4, is: d4

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi  x0  Rb cosðh1  hN ÞÞ2 þ ðyi  y0  Rb sinðh1  hN ÞÞ2

ð18Þ The parametric variable, h 5, represents the difference at the current z-position (h 1) and the point on the helix edge closest to the pellet: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sinð2/b Þ h 5 ¼ h1  ð19Þ d42  DR2 2C where, DR, is the difference between the current distance from the pellet center to the center of the screw, R xy, and the screw radius (which is equal to the barrel radius, R b, when a zero gap clearance is assumed between the flight and barrel). The coordinates of the point on the edge closest to the pellet (x e, y e, z e) are then: xe ¼ Rb cosðh5  hN Þ þ x0 ye ¼ Rb sinðh5  hN Þ þ y0 ze ¼ Ch5

2.2.2. Contact detection with the helix edge and top of the flight Due to the finite thickness of the flight used in this model, particle contact with the top of the flight, as well as the helix faces (both advancing and retreating) must be considered. For a particle contact to occur with the top of the flight, the pellet center must be located axially between the advancing and retreating helix faces. This condition can be determined by evaluating the cross product vector with respect to a pellet and the two helix faces. The cross product with the advancing flight is: np  na ¼ npx sinðh1  hN Þ  npy cosðh1  hN Þ

ð16Þ

ð20Þ

Therefore, the distance from the helix edge to the center of the pellet, d e, is: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð21Þ de ¼ ðxi  xe Þ2 þ ðyi  ye Þ2 þ ðzi  ze Þ2 A contact occurs if this distance is less than the particle radius. Modeling edges in DEM is critical for ensuring the stability of the simulation, in order to prevent excessive overlap and unrealistic solids behavior. For example, the particle-to-screw distance, R xy, must be less than the barrel radius, R b, for a contact to occur with the helical flight. If a particle is within one radius of this distance, a contact with the helix edge is possible. If the edge was not included, a

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

99

particle could migrate into the flight, producing a very large overlap. The relatively large contact force produced would result in an unrealistic acceleration, potentially leading to numerical instability.

nation of the area of contact through which the heat may flow. The heat conducted between two collided spherical particles i and j may be approximated by [22]:   q ¼ 2ka tj  Ti ð23Þ

2.3. Particle motion

where k is the thermal conductivity of the solid, and a is the contact radius. In the case of collisions with curved surfaces, such as the barrel wall, the contact area may not be circular. However, in our model the contact area is assumed to be circular since the degree of overlap is small, allowing determination of the contact radius, a, according to: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð24Þ a ¼ 2aRp  a2

The net normal and tangential forces resulting from particle –particle and particle – wall interactions are summed for the current spherical particle. The translational and rotational acceleration of the particle are then determined from these forces according to Newton’s Second Law of Motion, and integrated to determine the linear and angular velocity in the next time step. The linear velocity is integrated further to provide the future position of the sphere. While second-order and higher predictor – corrector integration schemes can be applied to granular dynamics [18,19], a straight-forward Fleap frog_ algorithm [20] is adopted for this simulation so that a greater number of particles can be simulated. The size of the time step has been shown by previous authors [21] to have critical importance to the stability of the discrete element method. The time step, Dt, is a function of the smallest particle mass, m, and the stiffness constant, K 1: rffiffiffiffiffiffi m Dt ð22Þ K1

The temperature of each particle is updated at the end of a time step according to the sum of heat flow in and out of that particle due to immediate contacts. For accuracy of this method, the time step must be chosen with consideration of the rate of heat transfer so that changes in temperature do not influence solids beyond the immediate particle neighbors at any instance in time. Fortunately, the appropriate size for the integration time step was found to be several orders of magnitudes greater than the time step required to ensure stability of the solids flow in the DEM model. Therefore, inclusion of heat transfer in the simulations was done at minimal computational expense. 2.5. Simulation parameters and details

-5

In these simulations, a time step of 4.72  10 s was used for screw speeds below 200 RPM. However, the time step was lowered for higher speeds, in order to account for the increasing flight velocity. 2.4. Heat transfer within the solids assembly Incorporation of heat transfer into the DEM extruder model allows for a reasonable estimation of the temperature profile within the solids assembly based on the degree of compaction of the solids assembly. It is assumed in our model that heat transfer within the solids takes place solely by conduction and since the thermal conductivity of the solid polymer pellet is much greater than that of the stagnant interstitial medium, then only heat transfer between particle contacts needs to be considered. It is also assumed that the temperature throughout each pellet is spatially uniform. This latter assumption is made for computational simplicity and is reasonable for the small diameter of the particles studied. Indeed, Vargas and McCarthy [22] found this approach to be effective in quantitatively determining the temperature field and effective thermal conductivity of packed beds of spherical particles under varying loads. The amount of heat transferred across the particle – particle or particle – wall interface is dependent on the contact force imparted during a collision. The contact force dictates the degree of particle overlap, a, in the DEM model (as already discussed) and hence, allows for the determi-

The geometry of the feed hopper and solids conveying zone is illustrated in Fig. 1. The solids conveying zone extends five barrel diameters (5 L/D) in front of the feed opening, which is 2 L/D in length. The rectangular feed opening is 0.98 L/D wide in this model. The geometric and material parameters used in the simulation runs are summarized in Table 1. The material studied was high density polyethylene (HDPE) in pellet form, modeled as 3.0 mm diameter spherical particles. This size was chosen to reflect the common diameter of commercially available polymer pellets. The density of HDPE was 945 kg/m3, and the coefficients of friction (l s and l b) were set to 0.25 and Table 1 Geometric and material parameters used in the DEM simulations

HDPE

Extruder

Parameter

Value

Diameter of pellet (D p) Density of a pellet (q s) External coefficient of friction (barrel) External coefficient of friction (screw) Stiffness coefficient (K 1) Poisson ratio (r) Coefficient of restitution (e) Barrel inner diameter (D b) Screw channel depth (H f) Flight lead (L s) Flight thickness (e) Zone length (includes length under feed opening)

3.0 mm 945 kg/m3 0.28 0.25 375 N/m 0.46 0.5 50.8 mm 10.16 mm 50.8 mm 5.08 mm 7 L/D

100

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

0.28 for the screw and barrel surfaces, respectively. These values for the external coefficient of friction fall within the range of 0.15 –0.30 for HDPE found in [23]. The internal (inter-particle) coefficient of friction was set to 0.29, which corresponds to an angle of repose of ¨16-. The initial stiffness coefficient, K 1, was 375 N/m for inter-particle collisions. In comparison to stiffness values used for metals and soil, this value is rather small; however, it provides realistic results in comparison to experimental data [10] and is similar to the value used for a 2-D simulation of a shear cell of couette type with photoelastic material [24]. More than 20,000 particles are necessary in order to fill the entire length of the extruder screw during flood-fed operation. The temperature of the screw and barrel are held at a constant 80 -C during the simulation runs, while the solids are fed in at a temperature of 25 -C. Under these conditions, 20 s of simulation time required about three days to run on a 2.0 GHz Pentium 4 processor with 1024 MB of DDR SDRAM. Due to the long computation time required for meaningful data, the DEM technique is not suited to commercial usage at this time; however, for research purposes, it is invaluable. This method can provide detailed information of the solids inflow and solids conveying zones that cannot be obtained with traditional plug-flow models.

restriction created by the torpedo caused the particles to completely fill the screw channel all the way back to the feed opening and was found to produce an axial pressure which exponentially increased towards the exit. Since this DEM model does not consider permanent deformation of the particles, the gap around the torpedo was set so that the pressure developed by the restriction was relatively small (15 kPa) compared to typical extrusion conditions. The restricted case was tested to show the influence of increased particle contacts (i.e. minor compaction) on the flow behavior. 3.1. Behavior of solids in the extruder Fig. 3 presents images of the solids-inflow zone of the model extruder demonstrating the time evolution of solids flow at 75 RPM for the open discharge case. The lower section of the hopper and the rectangular chute directly below, which is known as the feed opening, make up the

2.6. Measurement of variables In this paper, several variables are reported which includes pressure, temperature, cross- and down-channel velocity, residence time, and the coordination number. Only the coordination number and residence time are examined on an individual particle basis. The other variables are determined by a Fmixing-cup_ method. By this method, the reported variables are determined by time-averaging the individual property for each particle entering a control volume found within the discretized cross-section of the channel. Sampling of the properties of the particles is done at 2.5 L/D forward of the feed opening for the velocities, and at 4 L/D for temperature and residence time. For pressure measurements, the value is calculated by determining the total force exerted by the particles against the barrel at its top most point; the barrel location having a surface diameter similar to a commercial pressure transducer.

3. Results and discussion Two exit conditions have been examined in this work, namely the open-discharge and restricted case. For the exit condition of open discharge, the pellets were allowed to freely leave the end of the extruder, and axial pressure was observed to decrease from the minor pressure present at the feed opening to zero at the extruder exit. For the restricted case, a torpedo geometry introduced at the end of the screw (as shown in Fig. 1) produced a narrow gap between the screw and barrel through which the particles had to exit. The

Fig. 3. Simulation images at 75 RPM showing a side view of the feed hopper and the first section of the screw. The time (in seconds) for each image is: (a) 2 (b) 2.5; (c) 3; (d) 3.5; (e) 4; (f) 4.5; (g) 5; (h) 5.5.

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

‘‘solids-inflow’’ zone of the machine. To visualize the flow behavior of the solids within the model extruder, a tracer was introduced within the feed opening as a layer of black pellets (shown in Fig. 3(a)). Each subsequent image shown in Fig. 3 was taken at a 0.5 s interval. As the screw flight passed across the feed opening, an empty space was formed at the back of the feed zone of the screw behind the flight as illustrated in Fig. 3(c) and (f). As a consequence of the empty region in the screw being formed, the pellets were observed to be preferentially drawn into the screw channel from the back of the feed hopper; an observation similarly found from experimental trials with extruders [10], and with constant-pitch screw feeders [2]. The periodic filling of the screw channel by solids caused the bulk flow rate to pulsate. The frequency of these fluctuations in flow corresponds to the speed of the screw. The fluctuations are referred to in this work as ‘‘solids pulsing’’, and may have a significant contribution to a industrial phenomenon is known as ‘‘screw beat’’ which is an important characteristic of extrusion with respect to process stability. The time evolution of the solids flow shown in Fig. 3 also revealed a recirculation pattern in the feed opening. By tracking the path of the tracer particles, it was observed that some of the solids located at the back of the hopper (left-most side to each illustration in Fig. 3) did not fall into the screw as might be expected, but rather stayed in the feed opening and were pulled to the front. At the front of the feed opening, these particles were pushed upwards towards the hopper. This ‘‘recirculation’’ phenomenon has been attributed to the motion of the screw flight as it passes across the bottom of the feed opening [10,25]. The influence of this recirculation on the output stability of an extruder has been shown to be negligible at screw speeds below 200 RPM; however, at higher speeds, Potente and coworkers [14,25,26] have noted that it may lead to starvation of the screw. Fig. 4 shows images for the same simulation run presented in Fig. 3, focusing instead on the solids conveying zone forward of the feed opening. The black pellets present are those of the tracer introduced in the feed opening, shown in Fig. 3. Observation of the tracer pellets revealed that the solids in the center of the channel were insensitive to the flow fluctuations of the screw (i.e. solids pulsing), and no mixing was evident. The pellets in this region could be described as a plug or wedge. The pellets close to the advancing and retreating helix faces appeared to be freeflowing, not participating in the plug-like behavior of the solids in the center of the channel. Tracking the black pellets in this region of the screw channel revealed that the flow near the flights was moving in surges, at a frequency which corresponded to the screw speed (i.e. solids pulsing) [10]. The presence of the tracer particles also revealed that the flow near the flight generated axial dispersion within the transported solids. The impact of this Fflight effect_ on the mixedness of the solids in comparison to the recirculation found in the solids-inflow zone, can be quantified by determining the cumulative residence time distribution

101

Fig. 4. Simulation images at 75 RPM showing a side view of the solids conveying zone after the feed opening. The time (in seconds) for each image is: (a) 2 (b) 2.5; (c) 3; (d) 3.5; (e) 4; (f) 4.5; (g) 5; (h) 5.5.

(RTD). A discussion of the cumulative RTD will be given in a later section. The pressure was monitored at locations along the top of the barrel at a rate of 10 measurements per screw rotation. Unlike other predicted parameters discussed in this paper where the cross-channel direction is taken perpendicular to the helix face, the pressure data presented is taken perpendicular to the screw center-line axis. This means of presenting the pressure data was done so that the results better represented the recorded measurements by a pressure transducer. Fig. 5 shows the cross-channel pressure profiles for the open discharge and restricted cases at three axial locations along the barrel. The pressures for the open discharge case were very close to zero, while the pressures for the restricted case showed significant development. Axially, the pressure developed exponentially for the restricted case. Similar for both exit conditions, though different in magnitude, it can be observed from the figure that the cross-channel pressure near the flights was always close to zero, and it increased towards the center of the channel. This cross-channel pressure profile matches with

102

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107 60 3 L/D, restricted case 4 L/D 5 L/D 3 L/D, open discharge 4 L/D 5 L/D

Pressure (kPa)

50 40 30 20 10 0 0.0

0.2

0.4

0.6

0.8

1.0

Normalized Cross-Channel Location Fig. 5. Cross-channel pressure distribution for the restricted flow (solid symbols) and open discharge (open symbols) cases at 100 RPM.

the visual observations made where the solids were freeflowing at the flights, and behaved as a plug in the center of the screw channel. The error bars shown in the figure depict the magnitude of the pressure fluctuations monitored over time (mostly attributed to ‘‘solids pulsing’’), demonstrating the ability of DEM to capture the dynamic nature of the solids within an extruder; dynamic phenomena not revealed by previous solids conveying models. 3.2. Coordination number The coordination number is defined as the total number of simultaneous contacts for a given particle, and the distribution shows the normalized frequency (or number fraction) of particles with a particular coordination number. The resulting normalized distributions for both open discharge and restricted cases at 100 RPM are shown in Fig. 6. The distributions shown reveal marked differences in the contact behavior of the solids based on the exit condition. Under open discharge conditions, an average of four simultaneous contacts was observed, with virtually no

particles experiencing more than eight simultaneous contacts. For the restricted case, the average coordination number increased to five simultaneous contacts, though a significant portion of the pellets experienced up to ten simultaneous contacts. The coordination number distribution can be viewed as an indicator of the level of solids compaction within the monitored section, with higher coordination numbers indicating higher compaction. Maximum compaction around a sphere without particle deformation allows for twelve simultaneous contacts based on a closest packing geometry. As would be expected for the restricted case compared to the open discharge case, the higher pressure within the solids assembly led to closer packing of the particles, and hence to a higher number of simultaneous contacts. The level of compaction in the assembly of solids may also be adversely affected by screw speed due to channel starvation—a phenomenon that has been observed above 200 RPM for extruders. With increasing screw speed above 200 RPM, it has been shown experimentally that the degree of channel fill decreases [14,25,26]. The research found a relationship between the recirculation flow in the feed opening and the degree of channel fill observed for any screw speed. Therefore, as the screw speed increases and the degree of channel fill decreases, it is reasonable to assume that the number of contacts any particle will experience will also decrease. To study whether the effect of screw speed has an equivalent impact on particle contacts compared to the exit condition, a series of simulations were conducted using the open discharge condition at screw speeds from 501120 RPM. The resulting coordination number distributions for screw speeds up to 1120 RPM are shown in Fig. 7. At 50 RPM, the distribution resembles the Gaussian shape typical of random packings [27], and as the screw speed increases, the distributions shift toward lower coordination numbers. By 1120 RPM, the greatest fraction of pellets were taking part in binary collisions, illustrating the collisional nature of solids transport at very high screw speeds. Clearly, screw 0.5

Restricted Case Open Discharge

0.25

Normalized Frequency

Normalized Frequency

0.30

0.20 0.15 0.10 0.05

50 RPM 300 RPM 700 RPM 1120 RPM

0.4

0.3

0.2

0.1

0.0

0.00

0

2

4

6

8

10

12

Coordination Number Fig. 6. Coordination number distribution near the exit for the restricted flow (solid symbols) and open discharge (open symbols) cases at 100 RPM.

0

2

4

6

8

10

12

Coordination Number Fig. 7. Coordination number distribution vs. screw speed, up to very high speeds.

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

3.3. Down- and cross-channel velocity profiles Initially, it was thought that the solids in the channel may be circulating in the cross-channel direction. For this reason, the flow of the particles in the screw channel were originally discretized into two layers, one closer to the screw root and one closer to the barrel surface [10]. However, variation in the velocities with respect to depth was found to be insignificant, and so the two volume layers were combined into one with discretization of the screw channel only across the channel width. The normalized cross-channel velocity (V sx / V b, where V b = pD bN, the velocity at the barrel) distribution was determined under both exit conditions at 2.5 L/D forward of the feed opening for various screw speeds. A constant value of 0.28 was predicted across the cross-channel width for both exit conditions, with the largest variability (61%) in the normalized velocity seen closest to the advancing flight, and for the rest of the channel, the variability was less at approximately 42%. The variability in velocity was attributed to solids pulsing since the frequency of the velocity variations matched the rotation of the screw. The positive velocity signifies that, on average, the pellets were exclusively flowing towards the retreating flight in the cross-channel direction. The fact that this normalized velocity component was constant with respect to screw speed indicates that the extruder operated at constant specific throughput for the exit conditions examined. The normalized down-channel velocity (V sz / V b) was plotted in Fig. 8. Positive values indicate flow towards the exit, while negative values indicate back flow. It can be seen that there was a distinct influence of the flights on the downchannel velocity distribution. Near the advancing flight, the down-channel velocity was positive, while near the retreating flight, it was negative and of smaller magnitude. The down-channel velocity distribution away from the flights was relatively flat, which when considered with the crosschannel velocity distribution, indicates that the majority of the solids in the channel exhibit plug-like flow. The presence of plug-like solids transport gives some validation to the classical models for solids conveying. However, the noted problems with these models [3 – 6,28] in predicting output rate, pressure development, and screw torque suggests that DEM modeling may be a useful tool in improvements of our understanding of solids-conveying. It is possible that the effect the flight has on solids flow that

4

Normalized Down-Channel Velocity

speed has significant influence on the nature of particle collisions similar to that of the employed exit condition. Since both conditions of higher exit pressure and lower screw speed are typically attributed with improved process stability by an extruder, it may be inferred that the coordination number distribution is an indicator of process stability. A larger mean coordination number would suggest improved process stability, at least within the solids conveying region of the extruder.

103

50 RPM, restricted case 100 RPM 50 RPM, open discharge 100 RPM

3 2 1 0 -1 -2 -3 0.0

0.2

0.4

0.6

0.8

1.0

Normalized Cross-Channel Location Fig. 8. Normalized down-channel velocity (V sz / V b) vs. normalized crosschannel location.

was seen in this work, exists at higher levels of particle compaction-development of a contact force model capable of plastic deformation is therefore seen as the next step towards improvement of our fundamental knowledge of solids flow in an extruder. The variation in down-channel velocity over time, as shown by the large error bars, was more pronounced near the flights, which follows with visual observations that the solids in these regions appear to transmit the influence of ‘‘solids pulsing’’ to the latter zones of the extruder (i.e. melting and pumping). Screw speed and exit condition were examined for their influence on the down-channel velocity profile. The data shown in Fig. 8 suggests that the velocity distribution was not significantly affected by either factor, though, from the residence time studies in the following section, both screw speed and exit condition do impact the flow characteristics of the solids. 3.4. Residence time distribution The level of mixing in an extruder is a critical parameter for most industrial operations [29], even from the rather restricted viewpoint of only the solids, since it impacts product quality and the time required to switch between different resins. Previous discussions in this paper on the velocity distributions and visual observations have indicated that the solids do not move solely as a plug and that some mixing is occurring as they progress down the length of the extruder. A useful method for characterizing the level of mixing is the residence time distribution (RTD). Therefore, in this paper the cumulative RTD was evaluated to improve upon our fundamental understanding of solids motion within the screw channel, and to possibly determine the influence of the noted flow features (i.e. recirculation in the solids-inflow zone, and flight effects) on the level of mixing. The cumulative RTD, F(t), is defined as the fraction of flow that has exited the system by time t [30]. One of the benefits of the DEM approach is that the RTD can be determined directly from the simulation results; evaluating the time each tracer particle spends in the extruder model

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

F ðti Þ ¼ F ðti1 Þ þ

Ni NT

ð25Þ

where N T is the total number of tracer pellets, and F(t i-1), is the cumulative RTD up to t i-1. The function F(t) is plotted against the normalized residence time, k, which allows for comparisons between systems with different mean residence times. The residence time (t) is normalized by the mean residence time (t¯ ). The change in the shape of the cumulative RTD down the axial length of the screw was monitored by determining F(t) at five axial locations: immediately in front of the feed opening (z = 0), and at subsequent axial locations spaced 1 L/D apart. Fig. 9 shows the evolution of F(t) with respect to axial location for the restricted flow case in which the tracer pellets were introduced to the flow at the feed throat (similar to that shown in Fig. 3(a)). It can be seen that F(t) changes significantly in the axial direction, showing the increasing influence of the screw flight. Close to the feed opening, the cumulative RTD demonstrates plug flow-like characteristics with some distributive mixing evident as a result of recirculation in the solids-inflow zone. Progressively down the screw, each subsequent RTD indicated increased mixedness caused by the motion of solids at the flight. Similar results were obtained for the open discharge condition. Changing the location where a tracer is introduced into the system can help identify regions that display unique flow features [30]. In an attempt to resolve the individual contributions of the solids-inflow and solids conveying zones on F(t), the tracer was introduced into the flow at two locations: i) within the screw immediately beyond the feed opening, and ii) as a layer of tracer particles in the feed throat. Simulations at 50 RPM and 100 RPM under both open discharge and restricted flow conditions were carried out for comparison. Fig. 10(a) and (b) compare the

Cumulative Distribution Function, F(λ)

1.0

0.8

0.6

0 L/D 1 L/D 2 L/D 3 L/D 4 L/D

0.4

0.2

0.0 0.0

0.5

1.0

1.5

2.0

Normalized Time (λ) Fig. 9. Evolution of F(k) down the screw length for the restricted flow case, in which the tracer pellets were introduced at the feed throat.

(a) Cumulative Distribution Function, F(λ)

does not require any post-simulation calculations. To determine F(t), a time interval, Dt RTD, is designated, and the number of particles, N i, between t i and t i + Dt RTD are counted. The cumulative RTD, F(t i), is then:

1.2 Feed Throat, Restricted Case Screw Channel Feed Throat, Open Discharge Screw Channel Perfect Mixing Plug Flow

1.0 0.8

50 RPM Screw Speed

0.6 0.4 0.2 0.0 0.0

(b)

Cumulative Distribution Function, F(λ)

104

0.5

1.0

Normalized Time, λ

1.5

2.0

1.5

2.0

1.2 Feed Throat, Restricted Case Screw Channel Feed Throat, Open Discharge Screw Channel Perfect Mixing Plug Flow

1.0 0.8

100 RPM Screw Speed

0.6 0.4 0.2 0.0 0.0

0.5

1.0

Normalized Time, λ Fig. 10. Comparison of the cumulative residence time distribution, F(k), for the restricted flow (solid symbols) and open discharge (open symbols) cases with the tracers introduced at different locations. All data collected at distance 4 L/D in front of feed opening for (a) 50 RPM and (b) 100 RPM.

cumulative RTDs produced from the two tracer injection locations at a position z = 4 L/D forward of the feed opening for screw speeds of 50 and 100 RPM, respectively. The cumulative RTD for the ideal cases of plug flow and perfect mixing are included for reference to represent the extremes of mixedness. It can be seen from the figures that the location in which the tracer pellets are introduced has a marked affect on F(t). The distributions shown in the figures reveal that the majority of mixing occurs prior to the solids entering the solids-conveying zone. Distributive mixing experienced by the solids exiting the extruder appears to come largely from the solids in-flow zone. The evolution of mixedness down the screw channel length that was shown in Fig. 9 appears to indicate that distributive mixing was also present in the screw channel itself, though to a lesser extent than shown within the solids in-flow zone. Evidently, the flight effect and recirculation phenomenon noted earlier in this paper have a significant impact on the mixing within the solids as they are transported towards the melting zone of the extruder. With respect to the influences of screw speed and exit condition on the cumulative residence time, examination of both plots in Fig. 10 revealed that these factors had greater

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107 1.0

3.5. Temperature distribution in the screw channel 30

30

30

30

29 30

29

27

29

29

28

31

28 28

27 29 28

30

28

27

0.4

30

29

28

28

28 29

27 28 29

29

28

0.2 29

0.0 0.0

30

29

0.2

0.4

30

30

0.6

0.8

1.0

Normalized Cross-Channel Width Fig. 11. Temperature contour plot for the open discharge case at 100 RPM, and z = 4L/D.

significance on the solids flow than revealed by the downchannel velocity distribution. The distributions in the figure showed a minor improvement in mixedness with decreased screw speed (demonstrated by the RTDs produced from tracer injection in the screw channel). It can be seen for the 50 RPM case that even by k = 2.0 (twice the mean residence time) not all of the tracer pellets had progressed forward to the monitored location (4 L/D) yet in the case of 100 RPM all pellets had passed. For those distributions produced using tracer injected at the feed throat, the influence was either screw speed or restriction were not observable, likely due to the dispersive behavior of particles in the solidsinflow zone. The influence of screw speed can be attributed to the degree of channel fill which has previously been shown to decrease with increasing screw speed [10], thus lowering contact forces which would be experienced by those particles close to the flight faces and resulting in lower velocities. The influence of exit condition on mixedness of the system was less apparent in the Fig. 10. Overall, the results show that increased velocities of the solids next to the flight helix faces produce greater axial dispersion within the system and therefore increased mixedness that may be shown by the cumulative RTD. Increased mixedness of the solids within the screw channel is beneficial for uniformity of the temperature throughout the solids bed, which in turn leads to higher melting rates and a reduction of unmelted polymer in the extruded product—a fact clearly demonstrated in studies on starve-fed extrusion [31]. While distributive mixing throughout the solids bed would be necessary to bring about a uniform temperature distribution, the unique mixing observed in the RTD studies near the screw flights is likely to have an impact of the cross-channel temperature profile.

The temperature distribution of the solids in the solidsconveying zone is an important parameter for determining the onset of melting and the rate of melting once it begins. Unfortunately, there is no known non-intrusive instrument which allows researchers to determine the cross-channel temperature distribution throughout the solids as it evolves down the length of the solids conveying zone. The discrete nature of the present model offers a unique opportunity to explore some of the aspects of the temperature profile in the solid bed. The influences of screw speed (50 –125 RPM) and the exit condition (i.e. open discharge versus restricted flow) were studied in this work. Contour plots of the cross-channel temperature distribution for the solids predicted by the DEM model are shown in Figs. 11 and 12. Fig. 11 shows the contour plot for the crosschannel temperature distribution in the screw at an axial location 4 L/D forward of the feed opening for the open discharge case at 100 RPM, which is compared with Fig. 12 for the restricted case with all other conditions being the same. The origin for the two axes in each plot corresponds to the location where the advancing face of the flight intersects with the screw root. As seen in both figures, the temperature at the center of the solids bed (i.e. its core) existed at several degrees below the outer regions (i.e. near the flight, screw root, and barrel surfaces). The observed temperature gradient is typical of a slab-like body experiencing transient heat conduction, supporting the earlier observations that the majority of the solids in the screw channel exhibit plug-like flow. With increased screw speed, as would be expected, the solids temperature across the channel width was lower for any axial location due to the reduced residence time of the solids in the extruder. 1.0

32

Normalized Screw Channel Depth

Normalized Screw Channel Depth

31

0.8

0.6

105

0.8

32

31

31 30

31 30

29

29

29

29

30

30

28

28

28

27

0.6

29 27

27 29

28 27

0.4

27 28

0.2

29

29

0.2

0.4

30

30

29 29

29

28

28

30

0.0 0.0

28

31

0.6

31 32

0.8

1.0

Normalized Cross-Channel Width Fig. 12. Temperature contour plot for the restricted flow case at 100 RPM, and z = 4L/D.

106

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107

Comparing the two exit conditions, it appears that the level of solids compaction had no significant effect on the core temperature—both conditions exhibited similar core temperatures for the same screw speed. For clarity in this discussion, the term Fcentral region_ will refer to those solids within the core and near the screw root and barrel surfaces (i.e. particles located between 0.3 – 0.7 of the normalized cross-channel width). The more notable effect of the exit condition (and level of compaction) on the temperature distribution was observed at the walls. The highest temperatures in the screw channel were observed near the barrel and screw root surfaces for the restricted case (seen in Fig. 12). Interestingly, the solid particles next to the advancing and retreating flights (normalized cross-channel locations of 0 and 1, respectively) did not heat up as quickly as those next to the barrel and screw root, despite that fact at all surfaces were maintained at the same temperature. The lowest particle temperatures for solids located next to a heated surface were noted for the restricted flow case near the advancing helix face. Explanation of these observations can be attributed to two factors: level of solids compaction, and residence (contact) time. Higher compaction of the solids led to larger particle overlaps and increased contact area, particularly with the central region at the barrel and screw root surfaces (as indicted by the pressure distribution in Fig. 5), and therefore, allowed for larger heat transfer. Considering the relatively large pressures generated near the exit for the restricted case compared to the open-discharge case, the heat transfer was correspondingly larger for that exit condition due to increased contact area within the solids. The greater heat transfer observed at the walls of the central region was attributed to not only the larger pressures at that location but also to the near-zero down-channel velocities present. Slower down-channel velocities equate to longer contact times for heat transfer within the solids assembly. The direct influence of down-channel velocity on the cross-channel temperature distribution can be shown comparing the temperature of solids next to the two flight walls where the pressure was negligible. The higher temperature noted for the retreating face was attributed to the back flow of the particles, which bring solids of longer residence time and greater temperature back from further down the length of the screw. Based on these findings, it would be difficult to achieve a goal of improving uniformity of the cross-channel temperature distribution within the solids, based on the observed flow behavior in this paper. Compaction of the solids leads to greater contact area for heat transfer but also yields pluglike flow which yields a colder core of solids within the channel due to the insulative properties of polymers. The motion of the solids next to the flight helix faces showed the potential for improving uniformity of the temperature distribution; however, the influence of these regions was too small in comparison to the larger plug of solids near the center of the screw channel.

4. Conclusions A new 3-D model for solids conveying in a single screw extruder has been developed using DEM. The model has been shown to be an excellent tool for studying the local phenomenon of solids flow within the screw channel. Visual observations and evaluation of the cross- and down-channel velocity profiles has revealed two flow phenomena which have a marked affect on the transport behavior of solids within a single-screw extruder. The solids were observed to exhibit recirculation in the solids-inflow zone of the extruder which influenced the mixedness of the solids entering the screw and the residence time of the system. The other strong influence on the solids transport was the motion of particles close to the advancing and retreating faces of the screw flight. Those solids moving next to the flight contributed to the mixedness of the system. Despite the mixing observed in the solids flow within the solidsconveying region of the extruder, heat transfer into the particle bed was dominated by the plug-like nature of the system.

Acknowledgements The authors wish to thank the Natural Sciences and Engineering Research Council of Canada for their funding support of this work.

References [1] C.R. Woodcock, J.S. Mason, Bulk Solids Handling: An Introduction to the Practice and Technology, Chapman and Hall, 1987. [2] H. Colijn, Mechanical Conveyors for Bulk Solids, Elsevier, 1985. [3] Z. Tadmor, I. Klein, Engineering Principles of Plasticating Extrusion, Krieger Publishing Co., 1978. [4] W.H. Darnell, E.A.J. Mol, SPE J. (1956 April) 20. [5] G.A. Campbell, N. Dontula, Int. Polym. Process. 10 (1995) 30. [6] K.S. Hyun, M.A. Spalding, SPE ANTEC 43 (1997) 211. [7] S. Fang, L. Chen, F. Zhu, Polym. Eng. Sci. 31 (1991) 1117. [8] F. Zhu, L. Chen, Polym. Eng. Sci. 31 (1991) 1113. [9] S. Zhang, V. Sernas, SPE ANTEC 46 (2000) 94. [10] P.A. Moysey, M.R. Thompson, Polym. Eng. Sci. 44 (2004) 2203. [11] P.A. Cundall, O.D.L. Strack, Geotechnique 29 (1979) 47. [12] M. Moakher, T. Shinbrot, F.J. Muzzio, Powder Technol. 109 (2000) 58. [13] P.W. Cleary, Powder Technol. 109 (2000) 83. [14] H. Potente, T.C. Pohl, SPE ANTEC 47 (2001) 185. [15] O.R. Walton, Mech. Mater. 16 (1993) 239. [16] O.R. Walton, R.L. Braun, J. Rheol. 30 (1986) 949. [17] Y. Shimizu, P.A. Cundall, J. Engr. Mech. (2001 September) 864. [18] P.W. Cleary, Powder Technol. 109 (2000) 83. [19] X.M. Zheng, J.M. Hill, Powder Technol. 86 (1996) 219. [20] M. Moakher, T. Shinbrot, F.J. Muzzio, Powder Technol. 109 (2000) 58. [21] P.A. Cundall, O.D.L. Strack, The Distinct Element Method as a Tool for Research in Granular Media: Part 2, Technical Report, Department of Chemical Engineering, University of Minnesota, 1979. [22] W.L. Vargas, J.J. McCarthy, AIChE J. 47 (2001) 1052. [23] E. Gamache, P.G. Lafleur, C. Peiti, B. Vergnes, Polym. Eng. Sci. 39 (1999) 1604.

P.A. Moysey, M.R. Thompson / Powder Technology 153 (2005) 95 – 107 [24] [25] [26] [27] [28]

S. Schollmann, Phys. Rev. E. 59 (1999) 889. H. Potente, T.C. Pohl, Int. Polym. Process. XVII (2002) 11. H. Potente, T.C. Pohl, T. Vinke, Kunst. Plast. Eur. 91 (6) (2001) 22. Y.F. Cheng, S.J. Guo, H.Y. Lai, Powder Technol. 107 (2000) 123. I. Klein, SPE ANTEC (1977) 444.

107

[29] C. Rauwendaal, P. Gramann, SPE ANTEC 46 (2000) 111. [30] O. Levenspiel, Chemical Reaction Engineering, John Wiley & Sons, 1972. [31] M.R. Thompson, G. Donoian, J.P. Christiano, Polym. Eng. Sci. 40 (2002) 2014.