Discrete Element Simulation for Fundamental Process Understanding

The pharmaceutical industry handles large volumes of granular materials such as powder blends for tablet production and filled capsules everyday (1, 2). Slight changes in ingredient properties or process operation conditions can have a major effect on a finished product's quality. Given the market and regulatory uncertainties that are commonly associated with drug product development, pharmaceutical companies typically have several drugs in various developmental stages at the same time. Because of this volume, the industry must have computer-based rapid-prototyping tools that can efficiently capture and resolve the technical aspects of drug product development so that companies can confidently make decisions about drug portfolio management and planning (3, 4).

Rapid prototyping offers several potential benefits for design exploration and manufacturing process validation. Methods such as computer-aided equipment design, computer-aided manufacturing systems, mechatronics, enterprise software, and computer-graphics–based imaging have become widespread in health-related fields (5–7) such as device design and high-end surgery.

Partly as a result of advances in the process analytical technology movement, computer-aided modeling also can be used to track the motion of granular entities that govern the flow and mixing of the solid ingredients in drug manufacturing. Although practitioners have attempted to handle powder-flow problems using empirical approaches in the past, it is clear that enhanced predictive power—rooted in more detailed studies of powder-flow physics—is required to develop strategies for the systematic design, optimization, and control of pharmaceutical manufacturing processes. Over past 10 years, considerable work has been done in this direction. For example, Barker (8) and Bertrand (9) have conducted comprehensive surveys about computer simulations of the dynamic behavior of granular material, structure, flow, and mixing.

This article describes current, state-of-the-art modeling techniques and discrete element methods that are increasingly used to model and simulate pharmaceutical granular flows. We describe obtainable simulation results for free-flowing and cohesive granular material using rotary blenders as an example and emphasize current limitations and challenges that must be explored in future work.

Numerical models

Granular-flow models can be divided into three broad categories: continuum, kinetic theory, and discrete. Continuum models neglect the discrete nature of grains and assume a continuous variation of matter that obeys mass and momentum conservation laws (10–11). Thus, the material's behavior is described by constitutive equations, frequently ad hoc, that relate kinematic, mechanical, and thermal field variables. Some continuum models incorporate microstructural parameters that describe characteristics and mechanisms that are distinctive of granular materials such as the solid fraction distribution and dilatancy phenomena (12).

Kinetic theory models exploit similarities between interacting grains and colliding molecules in a dense gas (13). They incorporate energy dissipation and take into account fluctuational velocities with a Maxwellian distribution. They solve continuum hydrodynamic forms of mass, momentum, and kinetic energy conservation laws (14).

Discrete models have numerous subclassifications, but they all assume that the constituent grains are distinct and move according to prescribed rules. Examples of discrete models include Monte Carlo methods (15), which apply probabilistic rules; cellular automata (16), which use deterministic (and often simplified) rules; and particle dynamics, for which rules are derived from momentum conservation. Two simulation methods typically used for granular materials are based on hard- or soft-particle approaches.

Hard-sphere models assume that particles interact by instantaneous binary collisions and collisional operators to dissipate energy. They are used for diluted systems and thus are more appropriate for rapid granular flow (17).

In soft-sphere models, collisions are not instantaneous and last for a finite amount of time (i.e., enduring collisions). These models are more appropriate for dense granular flows (i.e., volume fractions larger than, or equal to 50%), thereby allowing the handling of multiple enduring contacts (18–20). In this situation, particles are allowed to have deformations, which are used to compute restoring elastic, plastic, and frictional forces.

When two particles collide, they usually undergo a deformation that may fall somewhere between the extremes of perfectly elastic or perfectly inelastic states. In the soft-particle model, collisions are inelastic, but the deformation is an overlap distance— the larger the value, the larger the repulsive force, and as a consequence, the larger the kinetic energy loss. In this process, the particles may also slide relative to one another, resulting in frictional forces. These considerations have been taken into account in the soft-particle method used in the study described in this article. The method also can be used to simulate chute flow (21), heap formation (22), hopper discharge (23), and flows in rotating drums (24).

In the present study, the granular material was a collection of frictional and inelastic spherical particles. Each particle may interact with its neighbors or with the tumbler walls under normal and tangential forces. The normal forces that develop between particles-in-contact are calculated using the Walton and Braun model's "partially latching spring" (25). Tangential forces are calculated using Walton's three-dimensional "incrementally slipping" model (26).

Previous work investigating cohesive flow involved the addition of liquid to granular systems (27). Nonetheless, the dynamics of pharmaceutical materials with dry cohesion are substantially different from those of wet cohesive materials. More recently, the cohesive force between particles was simulated using a square-well potential. To compare simulations with various numbers and particle sizes, the force's magnitude was represented in terms of the dimensionless parameter,

K = F_cohes ÷ mg

in which K is the bond number and is the measure of cohesiveness that is independent of particle size, F_cohes is the cohesive force between particles, and mg is each particles' weight. This constant force may represent short-range van der Waals or electrostatic effects.

Granular materials flow and mixing simulations

The soft-particle dynamics algorithm has been used to report the first three-dimensional simulation of a tumbling (double-cone) blender for free-flowing particles (28). Previous studies of mixing and segregation in the double-cone blender provided the experimental data needed to validate the models (29). According to these studies, blenders with reflectional planes of symmetry orthogonal to the rotation axis suffer impeded transport across the symmetry plane. A synopsis of early results appears in Figure 1.

Figure 1: (Top) Front view of a time sequence from a simulation of monodispersed particles differing in color in a double-cone blender. (Bottom) Corresponding experimental snapshots for a transparent blender with identical operational conditions. N denotes the number of revolutions.

In this article, we followed 15,000 identical 2-mm diameter spheres (44% of total fill volume) through 6 blender revolutions at 15 rpm. Particles were loaded first from above in random, nonoverlapping positions and were allowed to fall to rest before the simulation began. The simulations (Figure 1, top row) showed the essential convective transfer (Figure 1, bottom row). After 6 revolutions, the yellow particles appeared uniformly dispersed with the right half of the tumbler, but there was very little convective transport between the two halves of the tumbler. In both experiments and the simulation, particle motion consisted of an azimuthal convective flow accompanied by distinct axial diffusive transport.

Figure 2(a) shows the prediction of segregation patterns and displays a top view of a simulation using a bidisperse mixture of approximately equal masses of large and small beads. First, 11,500 small (6-mm diameter, blue) spheres were mixed randomly with 1500 large (12-mm diameter, red) spheres. The blend was rotated at 30 rpm for 12 revolutions. For comparison, a snapshot was taken of an experiment in which an equal mass blend of 1.6-mm (blue) and 4-mm (red) diameter glass beads were rotated at 16 rpm in a 1-qt double-cone blender (see Figure 2[b]). In both the experiment and the simulation, the tumbler was filled initially to 50% of its total volume.

Figure 2: (a) Top view of the segregation of 12-mm (red) and 6-mm (blue) spheres in a double-cone tumbler after 12 revolutions at 30 rpm. (b) Top view for comparison, using 4-mm (red) and 1.6-mm (blue) glass spheres. (c) Slice through the vertical plane passing through the blender´s axis of rotation, showing interior structure in the simulation. (d) Slice through the corresponding vertical plane in the experiment.

The same segregation pattern occurred in the tumbler's first few rotations: Differently sized particles predominated in a shell and core pattern along the blender's axis of rotation. In the experiment, the pattern changed little after the first few revolutions. Figure 2(c–d) shows the segregation patterns in the tumbler's interior for the simulation and experiment respectively. These figures show a sideview of the blender intersecting a vertical plane through the blender's rotation axis (indicated by green arrows in Figure 2[a, b]). The simulations and experiments showed that the small grains formed a contiguous core parallel to the rotation axis.

To produce the experimental view (see Figure 2[d]), a thin, metal plate was inserted along the vertical plane of interest to divide the blender into equal halves. Melted paraffin wax was poured on the top half of the blender. Then, the grains were vacuumed away from the other half (see Figure 2[d]). Once the paraffin settled, it held the grains in place so that the tumbler could be tipped with the metal plate facing up. Then, the plate was removed, and the topmost exposed layer was vacuumed to reveal the interior structure.

The previously used discrete element model was expanded later to include cohesive and frictional forces to characterize the behavior of cohesive systems. The model system used in the remaining examples shown here consisted of 20,000 2-mm diameter particles in a cylindrical vessel with a 9-cm diameter, 1-cm length, and frictionless end caps. Although the number of particles was much smaller than that used in experiments, the intensity of cohesive forces per unit volume was matched by direct observation of flow features. The model captured the main experimental observations.

In other experiments, the cylinder was loaded with the powder to 50% by volume and rotated at 7 rpm. The model displayed discrete, well-defined avalanches (see Figure 3) for the moderate cohesion scenario, with a K value of 60. The flow of particles in the simulation compares well with the flow of Avicel 102. The largest avalanches were at the onset of flow, settling to a smaller, nearly constant avalanche magnitude as the bed dilated. The volume-averaged cohesive forces decreased because of a smaller density of interparticle contacts.

Figure 3: Typical avalanching behavior, captured in a time sequence, displayed by the DEM model (20 rpm) for K = 60. The bottom row shows experimental snapshots (7 rpm) for A102 showing essentially identical behavior for comparison.

To test the model further, the material's dilation was measured in the simulation and compared with that of the experimental results. A custom-designed pixel-counting program was used to determine the relative volumes of powder to void. The bed dilation was uniform along the cylinder axis; visual observations support this assumption. Pixel data were analyzed for snapshots taken every 0.3 s for simulations and every 1 s for experiments. The bed's initial volume at time (t) = 0 s is denoted as Volume_initial and as Volume_new at later times.

% increase in bed volume = (Volume_new – Volumeinitial) ÷ Volume _initial

Dilation was quantified as the percent increase in bed volume relative to the initial volume. The last frames were taken at t = 0 before the onset of powder movement.

The dilation of a few other powders of differing cohesion were measured and compared with simulations of various K values. In Figure 4, the three powders of increasing cohesion yielded dilation values of approximately 10, 15, and 24; with the simulations at K = 45, 60, and 90, dilation results of 10, 14, and 20% were obtained, respectively. Both the simulations and the experiments show the same trend of increasing dilation with increasing cohesion. The slight differences in the magnitude of dilation may originate from the very different vessel–particle size ratios used in the experiments and simulations.

Figure 4: Effect of cohesion on dilation. Experiment (top) performed using powdered lactose, Avicel-102, and Fast-Flo Lactose. Simulation (bottom) corresponds to K = 45, 60, and 90.

The effect of cohesion on the mixing and segregation of the uniform binary mixture also can be captured by simulations. The computational model system for this case study consisted of 10,000 green and 10,000 red identical particles of 2-mm diameter, loaded initially side-by-side along the axis. The mixing drum in the simulation had a 9-cm diameter and a 1-cm length; It was rotated at 20 rpm. Four cohesion values were simulated with our model system. Snapshots for each value of K were compared (see Figure 5). Snapshots were taken at t = 0, 1, 3, and 5 revolutions.

Figure 5: This figure shows the effect of cohesion on mixing. The mixing of highly cohesive materials (K = 90) is slower than for cohesion-free materials (K = 0). For mildly cohesive materials (K = 0.1), the mixing is faster than for the free-flowing materials (K = 0).

The free-flowing mixture (K = 0) blended faster than the moderately cohesive (K = 60) and highly cohesive (K = 90) cases, for which red particles moved in unison, forming a stretchable band. The very mild cohesive case (K = 0.1) showed faster mixing than the free-flowing case, which is in agreement with previous experiments (30).

Conclusions

Particle-dynamics simulations were applied to characterize the flow and mixing in partially filled industrial tumblers with uniform free-flowing and cohesive binary mixtures. Our results show that dry cohesive powders dilate substantially and exhibit avalanches. Their size is governed by the intensity of the cohesive forces. The discrete element model captures the main phenomena observed in experiments. Mixing slows down for highly filled mixers. Axial mixing is much slower than radial and angular mixing. Particles of different sizes segregate in a shell-and-core pattern. Low-intensity cohesion enhances mixing. Despite these results, high cohesion values (typical of commonly used pharmaceutical powders) show slower mixing, which is attributed to the formation of plug flow in the cascading layer, causing sluggish thinning of striations of similar particles.

Although the examples used in this article focus on blenders, discrete element models also can be used to simulate flow, mixing, and segregation for other applications. By appropriately changing the boundaries and tuning of friction-cohesion forces, the model can simulate die filling, compaction and decompaction of the cohesive granular systems in tablet-manufacturing processes, hopper flow, wet and dry milling, agitated dryers, motion through transitions and feeders, and numerous other applications, leading to an enhanced understanding of the effects of particle properties, system geometry, and operating conditions on the outcome of manufacturing processes. The field is advancing rapidly and substantial progress in anticipated in the next few years.

M.S. Tomassone* is an assistant professor, B. Chaudhuri is a postdoctoral associate, A. Faqih and A. Mehrotra are graduate students, and F.J. Muzzio is a professor in the department of chemical and biochemical engineering at Rutgers University, Piscataway, NJ 08855, tel. 732.445.2972, fax 732.445.2581, silvina@soemail.rutgers.edu Dr. Muzzio also is a member of Pharmaceutical Technology's editorial advisory board.

*To whom all correspondence should be addressed.

References

1. M.E. Aulton, Pharmaceutics the Science of Dosage Form Design (Churchill Livingstone, Spain, 2002, 2d ed.).

2. G. Alderborn and C. Nyström,Pharmaceutical Powder Compaction Technology, Drugs and Pharmaceutical Sciences (New York, NY, Marcel Dekker Inc., 1996).

3. S. Karri et al., "Biopharmaceutical Process Development, Part III: A Framework to Assist Decision Making," Biopharm. Eur. (Sept.) 76–82, 2001.

4. M.A. Mustafa et al., "A Software Took to Assist Business Process Decision Making in the Biopharmaceutical Industry," Biotechnol. Prog. 20 (4), 1096–1102 (2004).

5. G.G.K. Jacob, C.C. Kai, and T. Mei, "Development of a New Rapid Prototyping Interface," Comput. Ind. 39 (1), 61–70, 1999.

6. M. Deppe et al., "Rapid Prototyping of Real Time Control Laws for Complex Mechatronic System: A Case Study," J. Syst. & Softwares,70 (3), 263 (2004).

7. S. Girod et al., "Computer Aided 3-D Simulation and Prediction of Craniofacial Surgery: A New Approach," J. Cranio-Maxillofacial Surgery 29 (2), 156–158 (2001).

8. G.C. Barker, "Computer Simulations of Granular Materials" in Granular Matter: An Interdisciplinary Approach, A. Mehta, Ed.(Springer-Verlag, New York, NY, 1994), p. 35.

9. F. Bretrand, L.A. Leclaire, and G. Levecque, "DEM-Based Models for the Mixing of Granular Materials," Chem. Eng. Sci. 60, 2517–2531, (2005).

10. M.A. Goodman and S.C. Cowin, "A Continuum Theory for Granular Materials," Arch. Rational Mech. Anal. 44, 249–266 (1972).

11. A.J.M. Spencer, "Deformation of an Ideal Granular Material," in Mechanics of Solids, H.G. Hopkins and M.J. Sewell, Eds., (Pergamon Press, 1982), pp. 610–652.

12. A. Mehta and S.F. Edwards, "A Phenomenological Approach to Relaxation in Powders," Physica A. 168 (2), 714, (1990).

13. S. Chapman and T.G. Cowling, The Mathematical Theory of Nonuniform Gases," (New York, NY, Cambridge University Press, 1970).

14. J.T. Jenkis and M.W. Richman, "Kinetic Theory for Plane Flows of a Dense Gas of Identical, Rough, Inelastic, Circular Disk," Phys. Fluids, 28, 3485–3493, (1985).

15. A. Rosato et al., "Monte Carlo Simulation of Particulate Matter Segregation," Powder Technol. 49 (1), 59–69 (1986).

16. T. Shinbrot, D.V. Khakkar, and J.M. Ottino, "A Simple Model for Granular Convection," Phys. Rev. Lett. 79 (5), 829–832, (1997).

17. C.S. Campbell and C.E. Brennen, "Computer Simulation of Granular Shear Flows," J. Fluid Mechan. 151, 167–188. (1985).

18. P.A. Cundall and O.D.L. Strack, "A Discrete Numerical Model for Granular Assemblies," Geotechnique 29 (1), 47–65 (1979).

19. C. Thornton and S. Sun, "Uniaxial Compression of Granular Media: Numerical Simulation and Physical Experiment," Powders and Grains, 93, 129 (1993).

20. G. Lian, C. Thornton, and M.J. Adams, "Discrete Particle Simulation of Agglomerate Impact Coalescence," Chem. Eng. Sci. 53 (19), 3381–3391(1998).

21. S. Dippel, G.G. Batrouni, and D.E. Wolf, "Collision Induced Friction in the Motion of a Single Particle on a Bumpy Inclined Plane," Phys. Rev. E 54 (6), 6845 (1996).

22. S. Luding, "Stress Distribution in Static Two-Dimensional Granular Model Medial in the Absence of Friction," Phys. Rev. E 55 (4), 4720 (1997).

23. G. Ristow, "Density Patterns in Two-Dimensional Hoppers," Phys. Rev. E. 50 (1), R5–R8 (1994).

24. C. Wightman et al., "Use of the Discrete Element Methods to Characterize Particles Are Flowing and Mixing in a Rotating and Rocking Cylinder," J. AICHE 44 (6), 1266–1276 (1998).

25. O.R.Walton and R.L. Braun, "Viscosity, Granular Temperature and Stress Calculation for Shearing Assembles of Inelastic, Frictional Disks," J. Rheol. 30 (6), 949–980 (1986).

26. O.R. Walton, "Numerical Simulation of Inclined Chute Flows of Mono-disperse, Inelastic Frictional Spheres," Mechan. Materials 16 (5), 239–247 (1993).

27. J.J. McCarthy, "Micro-Modeling of Cohesive Mixing Processes," Powder Technol. 138 (1), 63 (2003).

28. M. Moakher, T. Shinbrot, and F.J. Muzzio, "Experimentally Validated Computations of Flow, Mixing and Segregation of Non-cohesive Grains in 3D Tumbling Blenders," Powder Technol. 109 (1), 58–71 (2000).

29. D. Brone and F.J. Muzzio, "Enhanced Mixing in Double-Cone Blenders," Powder Technol. 110 (3), 179–189 (2000).

30. T. Shinbrot, A. Alexander and F.J. Muzzio, "Spontaneous Chaotic Granular Mixing," Lett. Nature, 397 (6721), 675–678 (1999).