5304
Ind. Eng. Chem. Res. 2010, 49, 5304–5309
Extraction of Transport Coefficients from Molecular Dynamics Simulations of Granular Flows: A Perspective Christine M. Hrenya* Department of Chemical and Biological Engineering, UniVersity of Colorado, Boulder, Colorado 80309-0424
Molecular dynamics (MD) simulations of solids flows provide a useful tool for extracting transport coefficients needed for the closure of continuum models. This “empirical” approach is particularly powerful for systems with complex interactionssnonspherical particles, cohesive particles, etc.sthat are difficult to handle with first-principles (e.g., kinetic-theory-based) approaches. Nonetheless, the extraction of such quantities is characterized by several subtleties which may lead to inaccuracies if ignored. In this work, four such complexities are illustrated: (i) the presence of clustering instabilities, (ii) the a priori need for the general form of the flux law, (iii) the presence of more than one transport coefficient in a given flux law, and (iv) the presence of Knudsen (higher-order) effects. 1. Introduction 1
Ever since the pioneering work of Cundall and Strack, molecular dynamics (MD) simulations have been a ubiquitous and valuable tool for bettering the understanding of a wide class of flows involving solids.2 Generally speaking, the types of studies pursued fall into two categories. The first category contains those studies in which MD simulations are used to simulate a specific unit operation or natural system. Typically the number of particles used in these simulations is far less than that of the actual system due to the high computational requirements. For example, a moderately sized, lab-scale bubbling bed (6 in. diameter and 10 in. bed height) operated with 100 µm particles contains ∼5 billion particles, which is well beyond the capabilities of current desktop computing. In the second category are those studies in which MD simulations of relatively simple geometries are used to probe certain phenomena on a smaller scale and/or generate benchmark data for purposes of validating continuum theories. These simulations involve a relatively modest number of particles, given that the primary goal is to determine averaged (constitutive) quantities. A contemporary variation on this latter theme is to use MD simulations of relatively simple systems to extract closures for the constitutive quantities, rather than just using this data as a testbed for existing continuum models. In other words, closures for continuum models can be derived empirically using MD data rather than from fundamental approaches like the kinetic theory analogy.3,4 Such an approach is particularly useful for systems which are difficult to handle from first-principles. Examples include rapid solids flows with cohesive forces between particles, nonspherical particles, etc., in which the collisional integrals arising from a kinetic-theory-based treatment are intractable without severe assumptions. Examples of this alternative (empirical) approach include simple shear flow simulations of cohesive5 and nonspherical6 particles to obtain closures for stress. The purpose of this perspective is to describe current challenges associated with the extraction of constitutive relations from MD simulations. The system considered here is the rapid granular flow (characterized by instantaneous and binary collisions) of inelastic, frictionless, monodisperse spheres. Albeit fairly ideal, this system serves as an illustrative example of the * To whom correspondence should be addressed. E-mail: hrenya@ colorado.edu.
challenges which are common to more complex systems as well such as cohesive interactions, polydisperse systems, nonspherical particles, etc. (Although cohesive flows involve agglomerates of particles that engage in sustained, many-particle contacts, such contacts have been approximated by a series of instantaneous, binary collisions. For example, a square-well approximation has been used to incorporate cohesive effects in MD simulations of granular5 and gas-solid7,8 flows.) In particular, four main challenges that arise in the MD-based extraction of constitutive quantities include (i) clustering instabilities, (ii) a priori knowledge of the general form of the flux law, (iii) determining multiple transport coefficients from a single flux law (more unknowns than equations), and (iv) Knudsen (higherorder) effects. 2. Governing (Continuum) Equations For the rapid granular flow of monodisperse spheres, the mass, momentum, and granular energy balances take the general form: Dn + n∇ · U ) 0 Dt
(1)
DU + ∇ · σ ) nF Dt
(2)
3 DT 3 ) -∇ · q + σ:∇U - nTζ n 2 Dt 2
(3)
F
where n is the number density; U is the solids velocity; T ) 1 /3m〈U′2〉 is the granular temperature, where m is the particle mass and U′ is the fluctuating component of the particle velocity (relative to the mean local velocity); F ) mn; σ is the solidphase stress tensor; F is a (conservative) external force like gravity; q is the heat flux of granular energy; and ζ is the cooling rate of granular energy due to dissipative (inelastic) collisions. Here, the stress σ, heat flux q, and cooling rate ζ are the constitutive quantities for which closures are required. Numerous kinetic-theory-based closures have been developed for this system over the past several decades,3,7-12 in which different assumptions used during the derivation process have led to different closures. The closures themselves take the form of equations of state and/or flux laws, as exemplified below. (The terms associated with divergence of the velocity field are
10.1021/ie901859v 2010 American Chemical Society Published on Web 02/16/2010
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
5305
omitted for purposes of simplicity and illustration in the analysis to follow. For further details on these terms, the interested reader is referred to Garzo´ and Dufty,9 for example.) σ ) pI - η[∇U + (∇U)T]
(4)
q ) -k∇T - µ∇n
(5)
ζ ) ζ(0)
(6)
where I is the identity matrix, p is the pressure, η is the shear viscosity, k is the conductivity, µ is the Dufour coefficient, and ζ(0) is the zeroth-order cooling rate. Consequently, closures for constitutive quantities σ, q, and ζ require expressions for p, η, k, µ, and ζ(0) in terms of the hydrodynamic variables (n, T) and physical parameters (particle diameter, mass, restitution coefficient, etc.), which are derived using the rigors of kinetic theory analysis. As an alternative approach, the direct extraction of the constitutive quantities σ, q, and ζ from MD simulations will be considered in turn below, illustrating each of the said complexities.
Figure 1. Sketch of a 2D slice of a (a) homogeneous cooling system (HCS) and (b) the same system at a later time, in which clustering instabilities have developed.
3. Extraction of Constitutive Quantities from MD-Based Simulations Hard-sphere MD simulations10 are used in this work since collisions are treated as binary and instantaneous in nature, which is coincident with the assumptions of rapid granular flows. The algorithms used to collect the kinetic and collisional contributions to σ and q, as well as ζ (collisional contributions only), have been described elsewhere11-13 and will not be repeated here. However, it is worthwhile to emphasize that only the constitutive quantities themselves (σ, q, and ζ) can be extracted directly from the MD simulations, whereas the corresponding transport coefficients (shear viscosity, conductivity, etc.) have to be “back-calculated”. This important distinction, which leads to some of the aforementioned complexities, is expounded upon below. For purposes of simplicity in the ensuing discussion, only dimensional quantities are considered. However, an appropriate nondimensionalization will reduce the size of the parameter space. This practice is standard when extracting constitutive quantities from MD simulations, and the interested reader is referred to the references included in the subsections below. Not surprisingly, certain geometries and flow types are better suited to the extraction of a given constitutive quantity than others. Accordingly, the following subsections are divided according to the constitutive quantity (σ, q, or ζ) being sought, along with the MD system being used for extraction and a discussion of the related complexities. 3.1. Cooling Rate. An ideal system for the MD extraction of cooling rate (ζ) is the homogeneous cooling system (HCS). As sketched in Figure 1a, this system is essentially a threedimensional box with periodic boundaries in all directions and zero mean flow. Particles are initially placed randomly in space and given a distribution of random velocities, while also ensuring zero mean flow in any direction. Within a few collisions per particle, the system achieves its characteristic cooling state.14 As expected, the granular temperature decreases with time due to the dissipative nature of the inelastic collisions. Accordingly, the cooling rate ζ also decreases with time, which precludes the collection of this quantity by tracking energy loss per collision and then dividing by the time-averaging period.12 Nonetheless, on the basis of a combination of Haff’s law (ln(T) decays linearly with collision number in HCS) and dimensional
Figure 2. Temperature decay obtained from a hard-sphere simulation of a HCS with an initial temperature of To, a solids volume fraction of 0.3, a restitution coefficient of 0.85, and 1000 particles. Also shown are the predictions from Haff’s law, characterized by a constant slope.
arguments, ζ can be determined accurately at any instant in time based on the value of T at that same instant.13 A clear advantage of the HCS system is that no spatial gradients exist for the hydrodynamic variables (n and T), and thus cooling rate is constant throughout system, making its extraction as described above straightforward. Namely, T is determined by appropriate averaging over the entire domain. Nonetheless, such HCS systems are known to display clustering instabilities,15,16 as sketched in Figure 1b. Clearly, different local values of n and T exist within and outside of a cluster (T is lower within the clustered region due to the increased number of dissipative collisions). It follows that ζ, which is a function of n and T, will also vary throughout the system. Hence, the extraction of ζ as described above (averaging over entire domain) would lead to inaccurate measurements when clusters are present in the system. Fortunately, two related options are available to ensure that measurements are extracted only from systems in which no clusters are present. The first is based on Haff’s law,2,4 which indicates that the ln(T) versus number of collisions is constant for HCS. As portrayed in Figure 2, a deviation from this law occurs as the clustering instability forms,17 and thus the extraction of ζ should be restricted to times prior to the onset of clustering. As a second option, the size of the system can be made small enough to suppress the formation of clusters.13 Either way, the resulting extraction of ζ from the HCS will result in “clean” data and an accurate constitutive expression. 3.2. Stress. For the extraction of stress (σ), an oft used system is that of simple shear flow (SSF), which can be achieved in MD simulations by means of Lees-Edwards boundary conditions.18 As displayed in Figure 3a, simple shear flow, unlike its Couette counterparts, has the advantage of a linear velocity
5306
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Figure 3. (a) Schematic of a two-dimensional slice of a simple shear flow system (SSF) and (b) snapshot obtained from a hard-sphere MD simulation of two-dimensional SSF, in which clustering instabilities are evident.
Figure 4. Variation of nondimensional shear stress obtained from a hardsphere MD simulation of SSF with number of particles (i.e., system size). The shear stress is nondimensionalized using the material density, particle diameter, and shear rate. Other system parameters are solids volume fraction of 0.3 and restitution coefficient of 0.85.
gradient perpendicular to the direction of mean flow. For such an ideal flow, the production of granular temperature via shear is constant throughout the system, as is its dissipation rate, resulting in a constant n and T across the domain. This lack of gradients in n and T allow for the collection of the kinetic and collisional contributions to the stress without the need for local averaging. In practice, however, this ideal flow is difficult to achieve again due to the presence of clustering instabilities, as seen in Figure 3b and originally documented by Hopkins and Louge.19 Similar to the HCS, the presence of clusters leads to local values of n and T which differ depending on whether the measurements are taken within or outside of a clustersnamely, the interior of a cluster is characterized by larger n and smaller T than the surrounding, diluter region.20 As was first shown by Liss and Glasser,21 an “effective” stress obtained by averaging across the domain (i.e., ignoring the local spatial gradients which develop) increases in magnitude with system size due to the presence of clusters, and reaches an asymptotic value for continued increases in the system size (or, equivalently, the number of particles). This behavior is demonstrated in Figure 4. Compared to the HCS system, techniques for avoiding this clustered state in SSF have not yet been reported, though the effects of clusters can be minimized by restricting the data collection to relatively small systems in which clusters are suppressed, and/or focusing on portions of the parameter space where the instabilities play a smaller role (e.g., nearly elastic particles). 3.3. Heat Flux. Unlike the HCS and SSF systems considered above, the system chosen for the extraction of heat flux (q) must have gradients in granular temperature, and perhaps the simplest
Figure 5. Diagram of a two-dimensional slice of a bounded conduction system in which the left and right walls are characterized by a cold (TC) and hot (TH) constant temperature, respectively, and all other boundaries are periodic. At steady state, the zero-mean-flow system displays concentration (n) and granular temperature (T) profiles with a local maximum and minimum, respectively, just inside of the cool wall.
system that fits this bill is the bounded conduction system depicted in Figure 5. For this system, two opposite walls are kept at constant granular temperatures, while the remaining four boundaries are periodic in nature. Zero mean flow is also maintained by setting the average system velocity to zero in the initial state and ensuring that no convective patterns develop through the course of the simulation. At steady state, a temperature profile develops between the walls, with its maximum located at the hotter wall (TH) and its minimum located near to, but not at, the colder wall (TC).22 The location of this minimum away from the cold wall is due to the dissipative nature of particle-particle collisions, whereas each wall, regardless of its temperature, serves as an energy source for the system. As expected, the concentration (n) profile is opposite in shape to that of the T profile, with the highest concentration of particles located in the low-T region, and vice versa. Accordingly, local averaging of MD-extracted quantities in the direction perpendicular to the constant-T boundaries is required. This averaging can be practically achieved by dividing the domain into strips parallel to said boundaries, with strip widths small enough to ensure that the resulting n, T, and q profiles are insensitive to further decreases in width size.23 Although a welcome advantage of the bounded conduction system compared to HCS and SSF is the absence of clustering instabilities (i.e., HCS and SSF have clusters which are transient in nature, whereas the bounded conduction has a steady, though nonuniform, distribution of particles), three additional complications arise. First, the form of the flux law must be known a priori. Although the form of the flux law for stress appears well agreed upon, the same is not true for the heat flux. In particular, two forms are common in the literature (see refs 9 and 24, for example), namely q ) -k∇T
(7)
q ) -k∇T - µ∇n
(8)
and
Recall that MD simulations can be used to obtain the constitutive quantity itself (q), as well as the hydrodynamic variables (n and T). With this in mind, it is clear that the value of the conductivity k calculated from the MD values for q, T, and n would be different depending on which of the two equations above were used. The results of MD simulations25 and experiments26 may provide an indication of the correct form. For example, a similar MD system was examined25 in which a nonzero flux is present for a system with ∇T ) 0, indicating that the µ∇n term is nonzero. Nonetheless, care must be taken to ensure that other effects, such as the Knudsen-wall effects
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
Figure 6. Diagram of a two-dimensional slice of a system similar to that of the bounded conduction show in Figure 5, but with a constant external force (Fext) added, which is acting right to left. As illustrated, the transverse location of ∇T ) 0 and ∇n ) 0 may no longer coincide for certain values of Fext.
described below, are not also playing a role, thereby muddying the interpretation. A second complication arises if a given flux law contains more than one transport coefficient, since the number of unknowns is then greater than the number of equations. For example, for the heat flux law given by eq 8, expressions for both k and µ are needed. Recall, however, that the MD simulations can be used to find q (or any constitutive quantity) directly, while values for k and µ need to be determined from this value. If the transport coefficients k and µ were constants, this task would be straightforward since two sets of conditions could be simulated (two different sets of wall temperatures, for example) resulting in two heat flux equations in two unknowns. Because k and µ are functions of the hydrodynamics variables (n and T), however, their determination is not straightforward. Several approaches to dealing with this complication of extracting q have been reported: (i) assume a form of the flux law which only contains one transport coefficient,27 (ii) assume k and µ are constant over a domain where n and T vary, and use a best-fit approach to determine k and µ,12 and (iii) extract µ from locations in the domain where ∇T ) 0 and ∇n * 0 such that the flux law reduces to q ) -µ∇n and thus only one unknown (µ) remains; once µ is determined, less restrictive systems can be examined to find k using the full form of the flux law (eq 8) and the known q and µ values. This latter approach, which is illustrated in Figure 6, is the most rigorous. Nonetheless, in-house simulations for this system indicate that this approach is not robust since the location at which ∇T ) 0 often corresponds to that at which ∇n ) 0, thereby precluding the isolated extraction of µ. The third complication associated with the extraction of q from the bounded conduction system involves the surprising role played by Knudsen effects. The Knudsen number (Kn) is defined as λ/L, where λ is the mean free path of the particles and L is a characteristic length scale. Two Knudsen numbers are relevant here: (i) Knwall ) λ/Lwall, where Lwall is the perpendicular distance from the wall, and (ii) Knbulk ) λ/Lgrad, where Lgrad is the a length scale characterizing the spatial gradient of hydrodynamic variablessa large Lgrad implies slowly changing variables (i.e., it takes a long distance for these variables to change their value by a given amount) and vice versa. The regions in which each of these Knudsen numbers are relevant is depicted in Figure 7. The Knudsen layer is that region of the flow where the role of particle-wall collisions cannot be ignored relative to particle-particle collisions, whereas particle-particle collisions dominate in the bulk interior. Extraction of constitutive quantities like q for purposes of formulating MD-generated flux laws is not appropriate in the Knudsen layer, since the near presence of the wall will alter the fluxes. To avoid the introduction of wall effects into constitutive relations not intended to include them, an estimate
5307
Figure 7. Knudsen layer and bulk interior for a wall-bounded flow of solid particles. Arrows indicate particle trajectory in near-wall region (Knudsen layer) is dictated by the wall and its material properties, e.g., roughness.
of the width of the Knudsen layer is needed. Such an estimate is possible by means of a heat-flux error analysis or via a ruleof-thumb indicating that wall effects occur for 1/Knwall < 2.5 (i.e., at distances less than 2.5 (local) mean free paths from the wall).23 Although the avoidance of Knudsen-layer effects in the MD extraction of q is thus fairly straightforward, the Knudsen effects present in the bulk interior present a bigger challenge. In this region, Knbulk is the relevant quantity since it dictates the form of the flux law. As a means of illustration, consider first that kinetic-theory-based models are typically based on a perturbative expansion about low Knudsen numbers (or equivalently “small gradients”): the Navier-Stokes-order constitutive relations contain only first-order gradients (as in eqs 4-6), the Burnettorder constitutive relations contain up to second-order gradients, super-Burnett-order constitutive relations contain up to thirdorder gradients, etc. Accordingly, the form of the flux law, and namely whether or not it contains “higher-order” effects (>1st order), is dictated by Knbulk. For the bounded conduction system under consideration here, the role of higher-order effects is found to be surprisingly high, leading to an error in q up to ∼40% if these effects are ignored.13 This point has relevance on MDextracted constitutive quantities since the form of the flux law must be known a priori, in order to convert measurements of q into the appropriate transport coefficients. It may be tempting to consider a less restrictive form of the flux law initially, that is, one which includes higher-order effects by adding higherorder gradients onto the expressions given in eqs 7-8. However, with such additional gradients also comes additional (unknown) transport coefficients, leading to an even greater disparity between the number of unknowns and equations. As discussed above, the handling of two transport coefficients in a flux law is nontrivial at best; handling more may become intractable on a practical level. 4. Concluding Remarks In this work, likely systems for the MD extraction of continuum (constitutive) quantities are examined in turn, in order to illustrate the subtleties which go hand-in-hand with such an extraction. In particular, a monodisperse, rapid flow of inelastic, frictionless spheres is considered, in which closures for the cooling rate (ζ), stress (σ), and conduction (q) are targeted. The systems considered for the extraction of these quantities are the homogeneous cooling system (HCS), simple shear flow (SSF), and bounded conduction, respectively. Four complexities associated with the MD measurements are shown to arise: (i) clustering instabilities, (ii) a priori knowledge of the flux law, (iii) flux laws with more than one transport coefficient, and (iv) Knudsen (higher-order) effects. The first complexity described is that associated with the clustering instability, which appears in both the HCS and SSF systems. The clusters create inhomogeneities in n and T, making the extraction of desired constitutive quantities (ζ and σ, respectively) nontrivial. In the HCS system, MD measurements
5308
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010
can be restricted to nonclustering systems by tracking the onset of clusters or performing simulations in a system small enough to suppress their onset. As a result, MD measurements of ζ are extremely “clean”. Similar techniques for completely eliminating clustering effects in the SSF system have not yet been reported, though the impact of clusters on the MD data extraction can be minimized by focusing on smaller systems and/or portions of the parameter space in which clusters play a smaller role (i.e., at higher coefficients of restitution). Nonetheless, if a higher level of accuracy is desired, more work is needed to best determine how to obtain MD-based σ data which is as clean as its ζ counterpart. Another complication which arises is an a priori need for the form of flux law, since the MD simulations provide direct measurements of the constitutive quantity only (ζ, σ, and q), rather than their corresponding transport coefficients. Although the form of flux laws for ζ and σ are generally agreed upon, the same is not true for the heat flux, in which two forms have been proposed for the monodisperse system considered here: q ) -k∇T and q ) -k∇T - µ∇n. Clearly, the value of the conductivity k determined from the MD value obtained for q will depend on which of these two forms are assumed. Both kinetic-theory and MD simulations can provide guidance as to which form is correct in a given situation, though caution should be exercised that other complexities (e.g., Knudsen effects) are not also playing a non-negligible role. The third difficulty associated with the MD extraction of transport coefficients occurs when the flux law contains more than one transport coefficient. Again, taking the heat flux law q ) -k∇T - µ∇n as an example, since only q is directly extracted from the MD simulations, the back-calculation of k and µ is straightforward only if a location where ∇T ) 0 and ∇n * 0 (or vice versa) can be identified over a range of desired of input parameters. Otherwise, the issue of more unknowns (k and µ) than equations (q) persists. Finally, Knudsen effects may also be present, which cannot be described by Navier-Stokes-order flux laws. Although Knudsen-layer (near-wall) effects can be effectively eliminated by avoiding collection of MD-extracted quantities in this region, Knudsen effects in the bulk (“large-gradient” effects) present a bigger challenge. As discussed for the bounded conduction system, the errors introduced if a Navier-Stokes-order heat flux law for q is incorrectly assumed can be up to 40%.13 Although the inclusion of higher-order spatial gradients in the flux law is a possible remedy, such an approach introduces additional transport coefficients (unknowns) into the flux law. Another possible remedy is the identification of a system(s) for the extraction of q in which such higher-order effects play a lesser role. This approach is promising since granular flows, like their molecular fluid counterparts, have been shown to be wellpredicted by Navier-Stokes-order theories even for systems in which higher-order effects are known to be present.28-30 In other words, the range of applicability of Navier-Stokes-order models is often wider than that expected based on the assumptions made during the derivation process. More work is needed in this area. Collectively, the systems examined here indicate that MD extraction of ζ is the most straightforward and accurate, the extraction of σ is less accurate though reasonable, while the extraction of q is the most tricky and subject to the greatest inaccuracies. As an aside, it is worthwhile to comment on the notion of geometrical independence of the constitutive relations. Note that the systems discussed here for the extraction of each constitutive quantity were originally chosen for the simplicity of the
corresponding extraction. For example, simple shear flow is chosen to extract the stress components since these components remain constant throughout the domain. The same is not true for the bounded conduction system, though stress can also be extracted from this system via the vertical-strip-averaging discussed above (for example, see Hrenya, Galvin, and Wildman,13 who extracted ζ, σ, and q from the bounded conduction system). Ideally, the corresponding transport coefficients obtained from each system would be equal for the same values of the hydrodynamic quantities (n and T). Practically speaking, however, since clustering instabilities are present in simple shear flow but not in bounded conduction, the transport coefficients also differ. A similar issue arises when one system contains Knudsen effects and the other does not. Hence, the complexities associated with the extraction of transport coefficients discussed above goes hand-in-hand with the (undesired) geometrical dependence of the extracted quantities. Finally, it is worth noting that the system considered heresthe rapid flow of identical, inelastic, frictionless spheressis among the simplest. A more complex, yet more realistic, system will have (i) additional constitutive quantities for which closures are needed and/or (ii) more complex forms of existing quantities. For example, in frictional31 or nonspherical particles, an additional balance for the “spin” granular energy is needed, in which closures for the flux of spin energy and couple stress are required. Similarly, for polydisperse flows,32 a closure is also needed for species mass flux, and the existing heat flux will depend on gradients in each of the species concentration, which gives rise to additional transport coefficients.33,34 In cohesive systems,35 an agglomeration rate is required to close a population balance. Despite such added complexities, the MD-based extraction of additional quantities like the agglomeration rate have been reported.36 Further progress along these lines is required for the continued progress of modeling such realistic systems, particularly when kinetic-theory-based approaches become essentially intractable. Acknowledgment The author is grateful to Daniel Cromer and Mike Pacella for generating the MD data for the HCS system. Special thanks are also extended to the funding agencies who provided support for portions of this work: the Department of Energy National Energy Technology Laboratory (DE-FC26-07NT43098) and the National Science Foundation (CBET-0318999). Literature Cited (1) Cundall, P. A.; Strack, O. D. L. A discrete numerical model for granular assemblies. Geotechnique 1979, 29, 47–65. (2) Po¨schel, T.; Schwager, T. Computational Granular Dynamics; Springer-Verlag: New York, 2005. (3) Gidaspow, D. Multiphase Flow and Fluidization; Academic Press: San Diego, CA, 1994. (4) Brilliantov, N. V.; Po¨schel, T. Kinetic Theory of Granular Gases; Oxford University Press: New York, 2004. (5) Weber, M. W.; Hoffman, D. K.; Hrenya, C. M. Discrete-particle simulations of cohesive granular flow using a square-well potential. Granular Matter 2004, 6, 239–254. (6) James, B.; Curtis, J. S. In The Effect of Particle Shape on ParticlePhase Stress, Metal Engineering International, Discrete Element Methods 07, Brisbane, Australia, August 27-29, 2007. (7) Weber, M. W.; Hrenya, C. M. Square-well model for cohesion in fluidized beds. Chem. Eng. Sci. 2006, 61, 4511–4527. (8) Weber, M. W.; Hrenya, C. M. Computational study of pressuredrop hysteresis in fluidized beds. Powder Technol. 2007, 177, 170–184. (9) Garzo´, V.; Dufty, J. Dense fluid transport for inelastic hard spheres. Phys. ReV. E 1999, 59, 5895–5911.
Ind. Eng. Chem. Res., Vol. 49, No. 11, 2010 (10) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (11) Campbell, C. S.; Gong, A. The stress tensor in a two-dimensional granular shear flow. J. Fluid Mech. 1986, 164, 107–125. (12) Herbst, O.; Mu¨ller, P.; Zippelius, A. Local heat flux and energy loss in a two-dimensional vibrated granular gas. Phys. ReV. E 2005, 72, 141303. (13) Hrenya, C. M.; Galvin, J. E.; Wildman, R. D. Evidence of higherorder effects in thermally-driven, rapid granular flows. J. Fluid Mech. 2008, 598, 429–450. (14) Dahl, S. R.; Hrenya, C. M.; Garzo´, V.; Dufty, J. W. Kinetic temperatures for a granular mixture. Phys. ReV. E 2002, 66, 041301. (15) Goldhirsch, I.; Zanetti, G. Clustering instability in dissipative gases. Phys. ReV. Lett. 1993, 70, 1619–1622. (16) Goldhirsch, I.; Tan, M.-L.; Zanetti, G. A molecular dynamical study of granular fluids I: The unforced granular gas in two dimensions. J. Sci. Comput. 1993, 8, 1–40. (17) Brito, R.; Ernst, M. H. Extension of Haff’s cooling law in granular flows. Europhys. Lett. 1998, 43, 497–502. (18) Lees, A. W.; Edwards, S. F. The computer study of transport processes under extreme conditions. J. Phys. C: Solid State Phys. 1972, 5, 1921–1929. (19) Hopkins, M.; Louge, M. Inelastic microstructure in rapid granular flows of smooth disks. Phys. Fluids A 1991, 3, 47–57. (20) Rice, R. B.; Hrenya, C. M. Characterization of clusters in rapid granular flows. Phys. ReV. E 2009, 79, 021304. (21) Liss, E. D.; Glasser, B. J. The influence of clusters on the stress in a sheared granular material. Powder Technol. 2001, 116, 116–132. (22) Galvin, J. E.; Dahl, S. R.; Hrenya, C. M. On the role of non-equipartition in the dynamics of rapidly-flowing, granular mixtures. J. Fluid Mech. 2005, 528, 207–232. (23) Galvin, J. E.; Hrenya, C. M.; Wildman, R. D. On the role of the Knudsen layer in rapid granular flows. J. Fluid Mech. 2007, 585, 73–92. (24) Jenkins, J. T., Kinetic theory for nearly elastic spheres. In Physics of Dry Granular Media; Hermann, H. J., Hovi, J. P., Luding, S., Eds.; Kluwer: Netherlands, 1998; pp 353-369. (25) Soto, R.; Mareschal, M.; Risso, D. Departure from Fourier’s law for fluidized granular media. Phys. ReV. Lett. 1999, 83, 5003–5006.
5309
(26) Martin, T. W.; Huntley, J. M.; Wildman, R. D. Hydrodynamic model for a vibrofluidized granular bed. J. Fluid Mech. 2006, 535, 325– 345. (27) Shattuck, M. D.; Bizon, C.; Swift, J. B.; Swinney, H. L. Computational test of kinetic theory of granular media. Phys. A 1999, 274, 158– 170. (28) Rericha, E. C.; Bizon, C.; Shattuck, M. D.; Swinney, H. L. Shocks in supersonic sand. Phys. ReV. Lett. 2002, 88, 014302. (29) Xu, H.; Louge, M.; Reeves, A. Solutions of the kinetic theory for bounded collisional granular flows. Continuum Mech. Thermodyn. 2003, 15, 321–349. (30) Wildman, R. D.; Martin, T. W.; Huntley, J. M.; Jenkins, J. T.; Viswanathan, H.; Fen, X.; Parker, D. J. Experimental investigation and kinetic-theory-based model of a rapid granular shear flow. J. Fluid Mech. 2008, 602, 63–79. (31) Jenkins, J. T.; Richman, M. W. Kinetic-theory for plane flows of a dense gas of identical, rough, inelastic, circular disks. Phys. Fluids 1985, 28, 3485–3494. (32) Hrenya, C. M., Kinetic theory for granular materials: Polydispersity. In Computational Gas-Solids Flows and Reacting Systems: Theory, Methods and Practice; Pannala, S., Syamlal, M., O’Brien, T., Eds.; IGI Global: Hershey, PA, in press. (33) Garzo´, V.; Dufty, J. W.; Hrenya, C. M. Enskog theory for polydisperse granular mixtures. I. Navier-Stokes order transport. Phys. ReV. E 2007, 76, 031303. (34) Garzo´, V.; Hrenya, C. M.; Dufty, J. W. Enskog theory for polydisperse granular mixtures. II. Sonine polynomial approximation. Phys. ReV. E 2007, 76, 031304. (35) Kim, H.; Arastoopour, H. Extension of kinetic theory to cohesive particle flow. Powder Technol. 2002, 122, 83–94. (36) Kantak, A. A.; Hrenya, C. M.; Davis, R. H. Initial rates of aggregation for dilute, granular flows of wet particles. Phys. Fluids 2009, 21, 023301.
ReceiVed for reView November 24, 2009 ReVised manuscript receiVed January 29, 2010 Accepted February 3, 2010 IE901859V
