1. Introduction
Possibly the biggest hindrance to an understanding of slow solids mixing is that there is no accepted set of governing equations -- as in the companion case of liquid mixing (Ottino 1990). Moreover, due to the fact that stress in granular materials is principally carried in 'stress chains' (Liu et al., 1995), it is not even assured that a continuum approach is generally applicable. This lack of a universal mathematical description can be attributed to the intrinsic physical complexity of the problems, and also, in part, to the difficulty in experimentally measuring the bulk properties -- stress, strain, voidage, etc. -- which would be necessary in a continuum description of a granular flow. While, there have been notable advances in non-invasive experimental methods (positron tomography, Broadbelt et al., 1993; nuclear magnetic resonance, Nakagawa et al., 1993; and gamma ray tomography, Nikitidis et al., 1994), which circumvent problems associated with bed opacity and allow for the measurement of valuable information -- particle positions, velocity profiles, and void fraction -- it is apparent that a continuum approach has achieved only limited success in mixing problems, the primary field of use being restricted to 'fast flow' situations (Campbell, 1989). Similarly, other approaches to modeling granular flow -- statistical mechanics (Mehta and Edwards, 1989), and cellular automata (Baxter and Behringer, 1991) -- while useful in certain regimes, are not general for all types of granular flow, and have not been successfully used as a predictive tool in granular mixing.
In contrast, the interactions of individual grains pose a much more tractable problem, both experimentally and theoretically, and much has been learned in this area (Hertz, 1881; Mindlin, 1949). The ever increasing power and efficiency of modern computers, coupled with the fact that particle-particle interactions are believed to be relatively well understood, has recently caused investigators to turn to molecular dynamics-like methods of investigation (Walton, 1993; Tsuji, 1993; and Lee and Herrmann, 1993). This type of model seems tailor-made for mixing applications, where particle properties are allowed to vary on a particle-by-particle basis and detailed, mixed structure is easily captured and visualized. The application of these methods to mixing, however, is still in its infancy and much remains to be done.
In this paper, we develop a model of granular mixing, verify it with experiments, and show how to extend the model in order that it may: (1) handle complicated geometries, (2) be applicable for 3-D mixers, (3) rapidly test mixing enhancement strategies, (4) incorporate different particle properties. Our philosophy throughout this paper will be to take the simplest description consistent with experiment and progressively expand its scope. Section 2 describes tumbling mixers and the geometrical model. Section 3 applies the model to a variety of simple 2-D geometries, comparing against experiment for mixing patterns, mixing rates, and mixing efficiencies. Section 4 explores 3-D models. Section 5 uses this method to test strategies for mixing enhancement in 2- and 3-D containers, stretching the ideas to include baffles and general non-convex domains, and time-varying avalanches. Section 6 points out some physical effects not captured by geometric modeling alone and discusses some limitations of the purely geometric viewpoint. After a brief review of particle dynamics computational methods in section 7, section 8 proposes a novel way to combine particle dynamics with geometric modeling to create a hybrid computational method that is highly efficient and, moreover, gets around many of the limitations highlighted in section 6. We conclude in section 9 with an outlook on promising areas of inquiry in granular mixing.
2. Mixing By Tumbling; A Geometrical View
In a slowly rotated container, the bulk of the granular material moves as a solid body, while a wedge-shaped region periodically avalanches. When the granular material exceeds its maximum angle of stability, often termed its dynamic angle of repose, the elevated material falls down the free surface until the material reaches its static angle of repose. The net effect of this avalanche is that a quantity of material -- a 'wedge' -- is moved from the upper portion of the surface to the lower portion of the surface (see figure 1). Although the exact mechanism of the motion during an avalanche and the detailed dynamics involved are not well understood (this will be discussed more fully in section 8), this simple geometrical fact yields valuable insights. The key observation is that the net transport in a slowly rotated container can be separated into two parts: transport of wedges (geometrical); and transport within wedges (dynamical). By decomposing the problem in this way, it is possible to study the implications of the geometry and dynamics separately and to add complexities in a controlled fashion. This allows even seemingly difficult problems (i.e., mixing in non-convex geometries, and mixing of dissimilar particles) to be probed in a simple yet methodical way.
3. Simple Geometries
Figure 2 illustrates the behavior predicted by geometry alone. As a mixer slowly rotates, the material inside moves with the mixer as a solid body until it reaches the dynamic angle of repose. At this point, a wedge of material avalanches so that the material surface returns to a stable configuration. As previously mentioned, for two similar, non-cohesive powders a random map can be used to determine where the individual particles go within the restabilized wedge. This simple model captures the essential behavior of the system, which is characterized by f , the fill level, i.e., the height to which the mixer has been filled with material. A comparison of mixing simulations with experiments in circular and triangular mixers is shown in figure 3.
The key to understanding the mixing process lies in the quadrilateral intersections between wedges. For each fill level, one can draw a series of wedge shapes to follow the motion of a wedge of material as it travels clockwise from its original position to its avalanche position. Doing this demonstrates the relevance of the fill level. At low fill levels, there are many wedge intersections, so a portion of any wedge -- a quadrilateral region -- eventually finds itself being part of another wedge, thus enhancing inter-wedge mixing. As the fill level increases, a point is reached where the wedge intersections vanish. The mixing is slowed since no wedge of material can mix with any other wedge. In fact, this fill level, whose value depends on the geometry of the mixer, constitutes a minimum in the mixing rate (Wolf, 1995). Once past this fill level, another phenomenon occurs: core formation. Here the wedges do not encompass the entire breadth of the material, and, as a result, a region of material -- the core -- does not get included in the avalanching process. This region only rotates with the mixer itself and does not get mixed at all. The core size increases with the fill level until finally the container is full, at which point no avalanching can occur and mixing stops.
From these observations of the model behavior, one can deduce that the efficiency of the mixer depends upon two important geometric entities: the fill level and the center of mass of the mixer. Core onset, the point at which there are no wedge intersections, occurs when the surface of the material remains above the center of mass of the mixer throughout the entire rotation. Obviously, material confined to the core does not participate in mixing, however, the material around the core will mix faster than at core onset as the fill level increases and wedge intersections re-emerge.
The degree of mixing can be measured in terms of interfacial length per unit area (Mohanty et al., 1982). As mixing proceeds, the number of contacts between dissimilar particles increases, hence the interfacial length per unit area -- or in 3-D, the interfacial area per unit volume -- increases. In the simulations, this length refers to the contact between particles on the imposed computational grid while in the experiments, this length is measured as the contact between pixels of a digitized image. The maximum value of the interfacial length per unit area for random mixing is one (Wolf, 1995).
For each fill level, the growth of the interface can be fit to an exponential curve. If the rate constant for the growth is multiplied by the amount of material being mixed, the product can be used as a measure of the overall mixing efficiency. Figure 4 shows computer simulation data of how the efficiency varies with fill level in the case of a circular mixer -- these data are in agreement with experimental results (Metcalfe et al., 1995). For a circular cross-section one can recognize core onset by the minimum, where mixing efficiency is nearly zero, which occurs at f = 0.5. It is interesting to note that the efficiency will only reach zero in the idealized case where the number of avalanches per revolution is an integer and the interface between the species lies on a wedge boundary. In real systems, neither condition will be met exactly and mixing will proceed slowly outward from the original wedge which contains the material interface. The data illustrate the presence of two maxima, one before and one after core onset. Clearly, the maximum which precedes core onset is more efficient and in general, except at fill levels very close to empty or to core onset, mixing is more efficient before core formation than afterwards. This trend holds true for all the other researched geometries as well. A complete account, including asymmetric geometries, appears in Wolf (1995). Let us consider here a few of the most important points.
As the symmetry of the geometry decreases, the amount of material mixed below core onset increases, increasing the overall efficiency. The results for a computer simulation of an equilateral triangle appear in Figure 5. The efficiency of the maximum before core onset is almost ten percent higher than that for the circle at an approximately equal value of the mass fraction. Thus, for mixers of equal area filled to the same mass fraction, the triangle mixes more efficiently than the standard circle found in industry. A general survey of results for the shapes researched appears in table 1.
4. Three Dimensional Geometries
Several studies have dealt with axial mixing of mechanically identical particles (Hirosue, 1980; Rutgers, 1965; Carley-Macauly and Donald, 1962; Hogg et al., 1966; Hwang and Hogg, 1980; and Kaye and Sparrow, 1964). These investigations have demonstrated that this process is well described by a cross-sectionally averaged diffusion equation. Thus for the case of mechanically identical particles, a 3-D map of the mixing within a wedge should have two parts, each described by random distributions. The radial mixing will be a uniform random distribution, while the axial mixing will be a Gaussian random distribution (figure 6).
By making a straightforward modification of the 'random map' of Section 3, diffusion parallel to the axis of rotation can be included and realistic simulation of 3-D wedges can be achieved. While the particles within the wedge still move randomly in the radial direction, the possible axial positions of the particles are weighted so that the probability of a particle moving a distance dz along the mixer axis obeys a Gaussian distribution. This model can be verified as follows: consider a cylindrical mixer whose initial condition has a marker particle concentration of 1.0 for -L < z < 0 and 0.0 for 0 < z < L. The centroid of the marker material begins at z = -0.5L and eventually decays to z = 0. By calculation of the centroid of either component, a quantitative comparison between theory, experiment, and simulation can be made. Figure 7 shows the position of the centroid of the marker particles versus non-dimensional time as calculated from a mixing simulation -- using this model, the analytical solution of the diffusion equation (Hogg et al., 1966), and experimental data (Hogg et al.,1966) -- time is made dimensionless with the mixer length and the diffusion coefficient from the analytical solution.
5. Strategies for Mixing Enhancement
Very little is known from a theoretical viewpoint about the effects of baffles on mixing of solids, even in cases restricted to 2-D situations. Our experiments indicate that baffles provide negligible function in most 2-D solids mixing applications -- with a subtle exception. Judiciously deployed baffles can erode the unmixed core through a curious mechanism. Figure 8 shows a pictorial representation of the observed effect. An even number of uniformly spaced baffles have no effect on the core; similarly baffles intruding into the mixing vessel do not affect the core behavior. However, an odd number of uniformly spaced, outwardly protruding baffles -- or any number of non-uniformly spaced baffles -- generate the following sequence of events. First, as shown in Fig's 8(a)-(b), the protruding baffle partially fills with grains following an avalanche. However a void (arrow) can be created when part of the baffle obstructs the flow. This void is typically sustained by arching of granular material above the baffle as the mixer continues to rotate (Fig's 8(b)-(c)). Finally, the arch collapses, and at that point the bulk of material in the container - including the core - shifts toward the baffle, as indicated in Fig. 8(d).
In controlled experiments in our laboratory, we have confirmed that although the cores of mixing containers with an odd number of uniformly spaced and protruding baffles steadily erode, no such effect is seen in vessels with simple intruding baffles, and their cores are static over time. Interestingly, vessels with uniformly spaced, even numbers of protruding baffles, do not erode their cores either. The reason for this is that in these vessels, each shifting event in one direction is paired with a shifting event in the opposite direction. So although the core may be slightly smaller than it would otherwise be, there is no steady motion of the core, and it remains intact for essentially all time.
A final test of the mechanism illustrated in Fig. 8 is depicted in Fig. 9. If one constructs a chiral mixing container -- such as the mixer in figure 9, it is possible to obtain different mixing behaviors for different directions of rotation of the container. If we rotate the container shown in Fig. 9 clockwise (top), the baffles should fill completely following an avalanche, while if we rotate it counterclockwise (bottom) a void can form within the baffles, causing the core to shift as shown. Experiments confirm that this effect does occur as predicted (figure 10).
Tumbling mixers have been fashioned in a variety of shapes and configurations over the years in an attempt to produce more efficient mixing (Wang and Fan, 1974; Carley-Macauly and Donald, 1962 and 1964). In most instances designs are ad hoc. However, with a full 3-D description of the geometrical effects in a tumbling mixer, new insights into the behavior of real slowly rotated mixers are available, and design improvements can be readily verified.
As was shown in section 4, axial transport in a tumbler limits the rate at which material can be mixed in a 3-D mixer. There have been, over the past decade, ample demonstrations of the beneficial effects of time periodic operation in improving mixing efficiency in viscous fluids (Ottino et al., 1992). Time periodic operation often leaves unmixed islands which may be destroyed by other modes of operation (Franjione et al., 1989) or combinations of two frequencies. Immediate parallels can be drawn with the present case of powder mixing (Wightman et al., 1995). Thus one method of overcoming the obstacle of slow axial mixing is to 'wobble' the tumbler -- for example, for every revolution of the tumbler, allow the axis of rotation of the mixer to sinusoidally move vertically. Figure 11 shows a schematic of this motion. In addition to avalanches perpendicular to the axis of rotation, which are caused by rotation alone, wobbling causes avalanches to occur parallel to the axis of rotation. These extra avalanches move material from one end of the tumbler to the other, whereby the rapid radial mixing will disperse the material more quickly than axial diffusion alone could, effectively enhancing axial transport.
Figure 12 shows the relative enhancement of axial transport as a function of the ratio of the number of avalanches -- wobbling to rotational. The cylinder has an aspect ratio of 2 to 1, and the mixer is filled to a fill level of f = 0.4. The wobbling motion is discrete -- i.e., in the one-to-one case, the cylinder is rotated until an avalanche occurs and then it is wobbled until an avalanche occurs and so on. It was found that this mixing protocol produces maximum mixedness at one wobbling avalanche per two rotational avalanches. This wobbling scheme is but one example of a mixer improvement whose effectiveness can be investigated quickly and easily with a 3-D model of a tumbling mixer.
6. Limitations of the Geometrical Viewpoint: Precession
Nevertheless, as we shall show, the core displays physical effects not contained in the simplest version of the model. In this section we show how the simplest geometric ideas about the core may be modified by inclusion of a surface boundary layer and a precessional rotation.
Further geometrical modeling can easily include the effect of the boundary layer. From geometrical considerations we can calculate how the core grows with f. By assuming a boundary layer of width d, we can calculate the reduction in size and a delay in onset due to the boundary layer. Then, by measurements of the core area vs. f we can back calculate d. Interestingly, for the half dozen or so materials we have examined in this way, the boundary layer is consistently about 6 particle diameters. Agreement between the calculated d and the visually observed boundary depth is a satisfying self-consistency check. Figure 14 shows a comparison of two square-shaped mixers, a modified simulation -- which includes a boundary layer -- and an experiment. In this simulation, the boundary layer is included as follows: A typical boundary layer depth is chosen (e.g., six particle diameters), the material down-slope of the avalanching wedge and within this depth is allowed to mix within itself, but not with the newly formed wedge, as it is sheared by the passing avalanche. There is one mechanism, however -- core precession -- which cannot be captured by geometry alone.
Core precession is illustrated in figure 15. After every full revolution of the container, the line demarking the initial interface between the colored grains should return to its original orientation. For a small number of revolutions, it does appear to do so. However, for a large number of revolutions -- the exact number of revolutions differs for each material, yet it is generally well beyond a practical number of revolutions for mixing -- the phenomenon is clearly evident. The line rotates past its original position in the direction of the container rotation. While precession does not affect mixing in a circular mixer, it is nevertheless important in non-circular containers whose cores are not rotationally symmetric. Precession rotates these cores into the mixing zones where the corners are sheared off, ultimately decreasing the size of the core.
Preliminary experiments have lead to three observations regarding core precession. (1) Precession depends on the container shape. Precession happens in containers of about equal width to height ratios -- the core does not precess in long, skinny containers. (2) The precession rate, i.e. how many degrees, Q, the core precesses per container revolution, depends on the fill level, f -- the rate is highest just above half full and falls to zero for a completely full container. (3) The precession rate depends on the material being mixed.
Figure 16 shows precession measurements for several materials. Precession is linear. (At least at the sampling rate used. We cannot rule out the possibility that finer measurements would reveal qp moving in small jumps.) Therefore, the precession rates can be obtained from the slopes of the lines in figure 16. There is an order of magnitude difference in the rate of change of the angle, dQ/dt, between small, cubic salt cubes and larger, spherical sugar balls. Since the precession rate depends on the material, core precession is an effect that does not seem explainable from geometry alone.
The physical mechanism driving core precession is at present obscure, but capturing the behavior may offer a validating test for future dynamical models of granular mixing, including the particle dynamics models we describe below.
7. Particle Dynamics
Particle dynamics simulations capture the macroscopic flow of the granular material by calculation of the individual particle trajectories. The time evolution of these trajectories then determines the flow of the granular material. Depending on the density and character of the flow to be modeled, different methods of calculating the trajectories are employed: a rigid particle model for low density, fast flow (the grain-inertia regime) (Campbell and Brennen, 1985); or a soft particle model for high density, slow flow (the quasi-static regime) (Walton, 1984).
In the quasi-static regime, such as that encountered in a slowly rotated container, lasting contacts are important and particle density is high, thus a 'soft particle' approach is used (Cundall and Strack, 1979). In this model, particle trajectories are determined by solving Newton's equation of motion for each particle, where the forces on the particles are described by a contact-dependent force law. The collisions/contacts in this model occur simultaneously and have finite duration. The particles are allowed to deform, or computationally overlap, and the restoring force of the deformed particle is generally chosen to agree with Hertz's (1881) contact theory, with a damping term to provide inelasticity. Similarly, the force on a particle subject to a tangential load is modeled after the work by Mindlin (1949). These force models have been essentially unchanged since Cundall and Strack's (1979) pioneering work. A notable exception is the model employed by Walton and Braun (1993). In their work, a partially latched spring is used in the normal direction, and an incrementally slipping friction model, fashioned after the work of Mindlin (1949), is used in the tangential direction.
Numerous investigators have utilized a 'soft particle' method in modeling a variety of applications, such as chute flow (Walton, 1993), heap formation (Buchholtz and Poschel, 1994; Lee and Herrmann, 1993), avalanching (Lee 1993), hopper flow (Tsuji 1993), and rotating containers (Walton and Braun, 1993; Tsuji, 1993). While more computationally intensive than the 'rigid particle' model, the 'soft particle' method is quite versatile, capturing both rapid and quasi-static flows with accuracy (Tsuji et al., 1992).
In what follows we employ force laws similar to those of Cundall and Strack (1979). Forces in the normal and tangential directions are modeled as non-linear and linear springs, respectively. Hertz (1881) found that the restoring force of two deformed elastic spheres in contact is given by:
where a is the difference between the particle separation and the sum of
their radii, and kn is the normal spring constant -- which is related to the
Young's modulus, Y, and the Poisson ratio, s, by:
where a is the radius of the particle. Similarly, a simplified version of
Mindlin's (1948) relation between tangential force and displacement is given
by:
where d is the relative tangential displacement of the particles since they
first came in contact with each other and kt is the tangential spring
constant -- which is related to the normal spring constant by:
where K is a proportionality constant which depends on the Poisson ratio.
Dissipation terms are then added to these forces to provide realistic
particle responses. A dashpot (viscous dissipation term) in the normal
direction provides collision inelasticity -- the form is chosen such that
normal collisions agree with experimental observations, while tangential
dissipation is provided by limiting the tangential force to the Coulomb
friction law. The complete force laws for the normal, fn, and tangential,
ft, directions then become:
where the negative sign in the tangential force expression denotes the direction opposite of d.
A schematic of these force laws is shown in figure 17. These forces depend on particle collisions/contacts and are integrated for each particle to determine the particle's trajectory. The collision/contact search is performed using an indexed nearest-neighbor algorithm, such that the computational time increases with the number of particles, N, approximately as N*log(N).
8. Merging Geometrical and Dynamical Approaches
In a slowly rotated container, where particle motion consists of discrete avalanches, non-trivial dynamics are principally confined to the avalanching wedge, and the large scale aspects of the problem can be easily handled by a geometrical approach. Confining the dynamics only where it matters -- the avalanching wedge of material -- results in substantial savings of computing time with respect to a full particle dynamics simulation.
This hybrid simulation technique utilizes the geometrical model of section 2 with the particle dynamics techniques of section 8. First, the gross motion of the particles -- solid body rotation -- and identification of the particles within the avalanching wedge are accomplished using the geometrical model alone. At that point, a particle dynamics simulation, which includes the particles within the wedge as well as a boundary layer of particles -- chosen such that increasing the size of the boundary layer caused no change in the result -- is used to determine the small scale dynamics of the avalanche itself. Once the avalanche is completed, the wedge particles are re-assimilated into the bulk, and the geometrical model is used once again. By including only those particles predetermined by geometry to be within the wedge (including a shear boundary layer of approximately 5-10 particles), the number of particles involved in a particle dynamics simulation is typically 1/15 that of the full mixer simulation -- this number depends on the fill level, f.
This increase in computational speed allows simulation of as many as 75,000 particles to be accomplished in a short period of time, while maintaining the versatility of a particle dynamics simulation. In addition, different particle properties, such as size/shape/morphology, are easily incorporated. Figure 18 shows a comparison of the mixing in an experiment and simulation. The experiment consists of small, cubic, dense particles and larger, spherical, less dense particles, the simulation, on the other hand, involves small, dense particles and large, less dense particles. Due to the fact that the parameters necessary to match the simulation to the experiment are not known -- i.e., Young's modulus, Poisson ratio, coefficient of restitution, angle of internal friction, etc. -- no attempt was made to achieve this match. Instead, a comparison is made here in order to show that it is indeed possible to model segregation effects similar to those found experimentally by using this method. It is apparent that in both the simulation and the experiment the combination of materials does not mix well -- radial segregation of the particles is clearly evident. This simulation technique easily captures the large scale, geometrical features, such as the core, as well as the smaller scale, dynamical characteristics like radial segregation.
9. Outlook
This model can be also be used to examine and screen methods of mixer optimization. Some of the findings are new: even numbers of uniformly spaced, outwardly protruding baffles -- or any number of inwardly protruding baffles -- do not affect mixing in 2-D systems; however, non-uniformly spaced, outwardly protruding baffles -- or odd numbers of uniformly spaced, outwardly protruding baffles -- can be used to erode the unmixed core. Even more unexpected is the finding that some chiral geometries -- such as "bent stars" -- may mix very differently according to the sense of rotation.
It has been demonstrated as well that axial mixing in a drum mixer can be greatly enhanced by altering the protocol of motion of the drum. By wobbling the drum -- sinusoidally moving the drum axis vertically -- axial mixing rates can be increased by as much as a factor of 25.
While the geometrical approach alone yields valuable insight into the mixing of granular materials, there are a number of effects which this technique fails to capture. Precession is but one of the material-dependent mechanisms which cannot be incorporated in a purely geometrical description of granular mixing. Similarly, particle segregation is not captured with this method either. However, by combining particle dynamics techniques with the geometrical approach -- in essence, by focusing the particle dynamics simulation only where it is needed -- a new hybrid method of mixer simulation, which is both more accurate than the geometrical method and much faster than the particle dynamics method, can be achieved. This novel numerical technique can easily and quickly examine a wide range of particle properties, so that the relative importance of each can be easily assessed. In doing so it may be possible to propose a criterion which will aid in determining, a priori, which combinations of material will be successfully mixed in tumbling containers.
Acknowledgment: This work was supported by the Department of Energy, Division of Basic Energy Sciences and by ALCOA. We wish to thank the support of Dr. Phil Hseih.