Particle pressure and marginal stability limits for a homogeneous monodisperse gas-fluidized bed: kinetic theory and numerical simulations

A linear stability analysis is performed for the homogeneous state of a monodisperse gas-fluidized bed of spherical particles undergoing hydrodynamic interactions and solid-body collisions at small particle Reynolds number and finite Stokes number. A prerequisite for the stability analysis is the determination of the particle velocity variance which controls the particle-phase pressure. In the absence of an imposed shear, this velocity variance arises solely due to the hydrodynamic interactions among the particles. Since the uniform state of these suspensions is unstable over a wide range of values of particle volume fraction φ and Stokes number St, full dynamic simulations cannot be used in general to characterize the properties of the homogeneous state. Instead, we use an asymptotic analysis for large Stokes numbers together with numerical simulations of the hydrodynamic interactions among particles with specified velocities to determine the hydrodynamic sources and sinks of particle-phase energy. In this limit, the velocity distribution to leading order is Maxwellian and therefore standard kinetic theories for granular/hard-sphere molecular systems can be used to predict the particle-phase pressure and rheology of the bed once the velocity variance of the particles is determined. The analysis is then extended to moderately large Stokes numbers for which the anisotropy of the velocity distribution is considerable by using a kinetic theory which combines the theoretical analysis of Koch (1990) for dilute suspensions (φ [Lt ] 1) with numerical simulation results for non-dilute suspensions at large Stokes numbers. A linear stability analysis of the resulting equations of motion provides the first a priori predictions of the marginal stability limits for the homogeneous state of a gas-fluidized bed. Dynamical simulations following the detailed motions of the particles in small periodic unit cells confirm the theoretical predictions for the particle velocity variance. Simulations using larger unit cells exhibit an inhomogeneous structure consistent with the predicted instability of the homogeneous gas–solid suspension.

A linear stability analysis is performed for the homogeneous state of a monodisperse gas-fluidized bed of spherical particles undergoing hydrodynamic interactions and solid-body collisions at small particle Reynolds number and finite Stokes number. A prerequisite for the stability analysis is the determination of the particle velocity variance which controls the particle-phase pressure. In the absence of an imposed shear, this velocity variance arises solely due to the hydrodynamic interactions among the particles. Since the uniform state of these suspensions is unstable over a wide range of values of particle volume fraction φ and Stokes number St, full dynamic simulations cannot be used in general to characterize the properties of the homogeneous state. Instead, we use an asymptotic analysis for large Stokes numbers together with numerical simulations of the hydrodynamic interactions among particles with specified velocities to determine the hydrodynamic sources and sinks of particle-phase energy. In this limit, the velocity distribution to leading order is Maxwellian and therefore standard kinetic theories for granular/hard-sphere molecular systems can be used to predict the particle-phase pressure and rheology of the bed once the velocity variance of the particles is determined. The analysis is then extended to moderately large Stokes numbers for which the anisotropy of the velocity distribution is considerable by using a kinetic theory which combines the theoretical analysis of Koch (1990) for dilute suspensions (φ 1) with numerical simulation results for non-dilute suspensions at large Stokes numbers. A linear stability analysis of the resulting equations of motion provides the first a priori predictions of the marginal stability limits for the homogeneous state of a gas-fluidized bed. Dynamical simulations following the detailed motions of the particles in small periodic unit cells confirm the theoretical predictions for the particle velocity variance. Simulations using larger unit cells exhibit an inhomogeneous structure consistent with the predicted instability of the homogeneous gas-solid suspension.

Introduction
One of the simplest and yet most stringent challenges in the modelling of gasparticulate flows is to predict the conditions under which a homogeneous fluidized bed will be unstable to volume fraction variations. Jackson (1963) showed that equations of motion which include particle-phase inertia and a drag coefficient that is a function of volume fraction lead to a prediction that the bed will always be unstable. This instability arises because the volume-fraction-dependent drag yields a kinematic wave speed that differs from the mean velocity of the particles. Therefore, volume-fraction waves propagate relative to a test particle and lead to a time-dependent drag force acting on the particle. The particle's inertia causes it to overshoot this forcing and move toward regions of higher volume fractions. Later studies (Anderson & Jackson 1968) introduced the concept of a particulate-phase pressure that tends to resist variations in volume fraction and thereby has a stabilizing influence. By a suitable choice of the particle pressure as a function of solids volume fraction it was possible to explain the stability of fluidized beds. However, the resulting theory is not predictive in the absence of a method of calculating the particle pressure from first principles. Batchelor (1988) suggested that the particle pressure be related to the hydrodynamic particle self-diffusivity in a low Reynolds number liquid-solid suspension and used measurements and plausible reasoning concerning this quantity to suggest a model for the particle pressure. While this suggestion has some physical basis, it ignores inertial effects and is most applicable to low-to-moderate Reynolds number liquid-fluidized beds. More recent studies (Didwania & Homsy 1982;Anderson, Sundaresan & Jackson 1995;Glasser, Kevrekidis & Sundaresan 1997;Goz 1995;Harris & Crighton 1994;Batchelor 1993) have explored the complex dynamics that occur beyond the onset of the initial instability without seriously re-examining the assumed form of equations of motion.
Recent progress in the kinetic theory of granular materials and theories and numerical simulations of gas-solid suspensions raises the prospect that the particle pressure may be derived from first principles allowing an a priori determination of the marginal stability limits for a monodisperse fluidized bed. The kinetic theory of granular materials or dense gases (Lun et al. 1984;Jenkins & Richman 1985) shows that a particle pressure arises when momentum is transported by the random translational motion of the particles and when interparticle collisions instantaneously transport momentum from the centre of one particle to another. A crucial step in determining the particle pressure is to predict the magnitude of the particle velocity fluctuations. In the theory of rapid granular flows, the particle temperature (or onethird the velocity variance) is determined by a balance of the work done to shear the fluid with the energy dissipated by inelastic interparticle collisions.  incorporated the effects of the interstitial gas in a study of rapid shearing motion of a particulate suspension. Exploiting the large O(2000) ratio of the particle to gas density, the inertia of the particles was retained while the inertia of the gas was neglected. The viscous dissipation due to the interstitial flow supplemented the dissipation due to inelastic collisions lowering the particle temperature relative to that for particles in vacuum. The rate of dissipation was determined from a simulation for a suspension of particles with a Maxwellian velocity distribution. A kinetic theory which incorporated this viscous dissipation rate was able to predict most of the rheological properties observed in full dynamic simulations of sheared suspensions with hydrodynamic interactions among the particles.
In a homogeneous fluidized bed, there is no shearing motion to provide the source in the fluctuation energy budget for the particulate phase. Koch (1990) noted that an alternative source of energy could arise from the random forces acting on a particle due to the hydrodynamic disturbances produced by its neighbours as they settle relative to the gas. Koch determined the magnitude of this source in dilute suspensions, φ 1, of very (St φ −3/2 ) massive particles for which the meanfree time between collisions is smaller than the particle viscous relaxation time m/(6πµa) and moderately (φ −3/2 St φ −3/4 ) massive particles for which the viscous relaxation time is shorter than the collision time. Here, φ is the particle volume fraction, St = mU t /(6πµa 2 ) is the Stokes number, m and a are the mass and radius of the particles, U t = mg/(6πµa) is the terminal velocity of an isolated particle, g is the acceleration due to gravity, and µ is the gas viscosity. Using this source to determine the particle temperature and pressure, Koch considered the limits of stability for a dilute sedimenting gas-solid suspension. A region of stability was predicted but only for very dilute suspensions φ 6 O(St −2/3 ). In this paper, Koch's (1990) theoretical predictions for the particle pressure and the marginal stability limits will be extended to higher particle volume fractions. At high Stokes numbers the source of energy is related to a time integral of the autocorrelation of the force felt by a test particle. Since the homogeneous fluidized state is generally unstable, this correlation cannot be determined by straightforward dynamic simulation of the particles' motion in a gas. We exploit the fact that in this limit the velocities of all the particles are nearly equal with the distribution of the small velocity fluctuation being nearly Maxwellian. The autocorrelation of the force can therefore be evaluated by computing the force on the particles having the same speed but moving according to a hard-sphere Maxwellian distribution. In § 2 we use a method of multipole expansion described in  and  for computing the force on the particles given their velocities to determine the force autocorrelation and hence the source of velocity fluctuations for φ ranging from 0.03 to about 0.56.
In § 3, we develop a kinetic theory for concentrated suspensions of moderate-Stokes-number particles whose velocity variance is significantly anisotropic. Koch (1990) determined the anisotropic source in this regime for dilute suspensions where the collisional distribution of fluctuation energy into vertical and horizontal velocity moments is negligible. We modify his analysis to account for the collisional redistribution of energy and the effect of particle interactions at finite φ on the source of energy.
Knowing the velocity variance, granular flow theory may be used to determine the particle-phase pressure and rheology. The resulting set of equations for the particle phase is analysed in § 4 to determine the criterion for the stability of homogeneous fluidized suspensions to small perturbations. It is found that these suspensions are unstable over a wide range of values of φ and St. While this is generally well known, it must be appreciated that the present study is the first one based on equations of motion in which all the important terms are derived from first principles.
Section 5 presents the results of dynamic simulations following the motion of individual particles. We first consider relatively dilute suspensions with a small system size for which we expect the suspension to be stable. The simulations for this case allow us to validate the theory for velocity variance over a wide range of Stokes numbers. Next, we consider more concentrated suspensions where the suspension is expected to be unstable even for small system sizes. Our simulations show clear evidence for the instabilities arising in these dense suspensions.

Particle velocity variance in a homogeneous suspension
We begin our analysis with the derivation of a general expression for the source of fluctuation energy in a homogeneous suspension. The velocity of a representative particle α sedimenting through a gas satisfies where m is the mass of the particle, g is the gravitational acceleration, F α v is the force resulting from the viscous stresses in the gas around the particle, and F α c is the collision force experienced by the particle when it comes in contact with another particle in the suspension. Note that the buoyancy force is neglected in the above expression since the density of the gas is typically much smaller than the density of the particle. The viscous force will in general depend on the velocity and position of all the particles in the suspension. We shall denote the average over all the particles in the suspension by an angular bracket. Since the sum of the collisional forces on the colliding particles vanishes, averaging (1) yields Note that the above applies strictly to homogeneous suspensions for which the spatial gradients in the average velocity and volume fraction vanish. The spatial derivative terms to account for slow variations in these quantities will be added to the above equation later in § 4 where we consider the stability of homogeneous suspensions to small spatial gradients of particle volume fraction. An equation for the fluctuation energy in homogeneous suspensions can be obtained now by multiplying (1) by v α and subtracting from its average an equation obtained by multiplying (2) by v . This yields where T = ( 1 3 )( v α2 − v 2 ) is the particle-phase temperature. At steady state the terms on both sides of the above equation must vanish for homogenous suspensions. It is convenient, to decompose the right-hand side into a number of terms representing the sources and sinks of fluctuation energy by various mechanisms. For example, in rapidly sheared gas-solid suspensions examined in , and in the rapid granular flow literature (e.g. Lun et al. 1984), the source of energy is expressed in terms of a product of the particle-phase stress and the strain rate while the sink is expressed in terms of energy dissipation by viscous forces acting on the particles and due to inelastic collisions. Expressing the right-hand side of (3) by S − Γ with Γ representing the sink of energy that occurs in the absence of shearing motion and relative motion between the phases, we obtain an expression for the source of energy fluctuations as (4) The dissipation rate is Γ = Γ vis + Γ inelas with the viscous and inelastic sinks given, respectively, by where χ is the value of the radial distribution function at contact. (Note that Γ in the present paper is equivalent to Γ /n in the notation of .) The above two expressions correspond to the energy dissipation per particle in a suspension with an isotropic Maxwellian velocity distribution with a zero mean velocity and 3T variance when the Reynolds number based on particle radius a and characteristic fluctuation velocity T 1/2 is small compared with unity and when e, the coefficient of restitution for particle collisions, is close to unity (nearly elastic collisions). R diss (φ) accounts for the particle interactions in Stokes flow. Expressions for this quantity and χ will be given later in this section.
With Γ specified, one can, in principle, compute S from dynamic simulations of sedimenting particles by evaluating the right-hand side of (4) using the force and velocity of each particle computed at every time step in simulations. This procedure is suitable provided that the suspension remains homogeneous. Since we are primarily interested in investigating the stability of the homogeneous state itself, this procedure cannot be adopted. It also does not allow one to develop a simple analytical expression for the source of energy as a function of Stokes number and other relevant variables. Therefore, we shall use a perturbation technique to determine S in the limit of large Stokes number for nearly elastic particles. An approximate theory will then be developed for the case of moderate Stokes numbers.

Suspensions of high-Stokes-number, nearly elastic particles
In a high-Stokes-number suspension, the particles respond only weakly to hydrodynamic velocity disturbances in the gas. Since these hydrodynamic interactions are the only source of velocity variance in a monodisperse, sedimenting suspension, the leading-order approximation to the particle velocity is simply a common mean velocity U = v obtained by solving (2). Koch (1990) showed that the particle velocities in a high-Stokes-number sedimenting suspension of perfectly elastic particles can be expanded in inverse powers of St −1/3 . Therefore, we expand the velocity of a particle in a series . . . It will be shown that, for large St and e = 1, v II /U and v III /v II are O(St −1/3 ). Similarly, we expand the viscous and collisional forces on particle α in a series Here, F α v,I represents the viscous force on particle α when all the particles are moving with the velocity U , F α v,II the force on the particle given the velocity of particle γ to be v γ II , γ = 1, 2, . . . , N, N being the number of particles in the suspension, and so on. Once again, F v,II /F v,I and F v,III /F v,II are expected to be of O(St −1/3 ) for elastic particles.
The mean velocity is controlled by a balance of the gravitational and drag forces acting on the particles as well as the particle inertia in cases of unsteady flows. Since the particles have a common velocity, this leading-order velocity does not contribute to the collisional force. The mean velocity is obtained from the average momentum balance for the particles (2).
The characteristic value of the collisional force F α c,II averaged over a small time interval is mv 2 II /a since the change in momentum upon a collision is O(mv II ) and the collision frequency is O(v II /a). Thus, the velocity fluctuations v γ II induce collisional contributions to the momentum balance that are O(µaUSt 1/3 ) and the leading approximation to the momentum balance is The viscous and gravity forces do not contribute to the above equation of motion indicating that to leading order the particle fluctuation velocity is similar to that of a hard-sphere molecular system. Together with the condition v α II = 0, (9) implies that the velocity distribution f(v II ) is an isotropic Maxwellian. If we further require that the variance in v α be the same as the variance in v α II , we have where T is the granular temperature.
Although the leading-order momentum balance (9) establishes that the particle velocity distribution is nearly Maxwellian, it is not sufficient for determining the granular temperature T . This is because (9) includes only conservative forces. The O(µaU) contributions to the momentum balance (1) are , which represents the fluctuation in the force driven by the mean velocity U of the particles. This force fluctuates because of variations in the spatial configuration of the particles. In a dynamic suspension, F α v,I − F v,I will vary with time as a particle translates relative to its neighbours. To summarize, the particle velocity v α consists to leading-order of a common mean velocity U , which may be determined by the mean momentum balance for the particle phase (2). The O(USt −1/3 ) correction to this velocity has a Maxwellian distribution, because the leading-order momentum balance (9) for the individual particles is dominated by elastic interparticle collisions. Finally, an O(USt −2/3 ) correction is driven by the fluctuating viscous force associated with the gas flow induced by the mean motion of the particle assembly motion. This fluctuating force changes with time as the particles move relative to one another with a motion that can be approximated using hard-sphere dynamics.
To see how v III determines the granular temperature of the suspension, we must examine the energy balance of the particle phase (3). The right-hand side of (3) is given to leading order by where we have used the fact that v α II is uncorrelated with the particles' spatial configuration and, hence, with F α v,I . The first term in the above expression simply represents the energy dissipation per particle in a suspension with the particle velocities given by an isotropic Maxwellian, and therefore represents −Γ (cf. (5) and (6)). Comparing the above with (4), we see that to leading order the source of energy fluctuation is given by whereF α v,I = F α v,I − F v,I and v α III can be obtained by integrating (11) over time. Denoting the right-hand side of (11) by F α , the above expression for the source can be expressed as Thus, we see that the source is related to the autocorrelation of the force felt by a representative particle. We expect the correlation time to be comparable with the O(a/T 1/2 ) time over which the spatial configuration of particles changes significantly.
We therefore write the source as the product of the square of the O(6πµaU) force and the O(a/T 1/2 ) correlation time scale, i.e.
In deriving (15), we have assumed that the average velocity of the gas-solid mixture is zero. However, this expression may be generalized to non-zero mixture velocities by replacing U by the relative velocity of the two phases. S * is the non-dimensional source which will be determined using numerical simulations in § 2.2. In summary, the average momentum and fluctuation energy equations for a homogenous fluidized bed in the limit of high Stokes number and nearly elastic particles are given by The steady-state particle-phase temperature for elastic particles is given by when e = 1 in agreement with our assumption that T /U 2 is small. In writing the above expression we have made use of the relation obtained from (16) at steady state and mg = 6πµaU t . Substituting this estimate of T into the expression (final term on the right-hand side of (17)) for the energy dissipated by inelastic collisions, we find that the particles may be treated as elastic provided 1 − e St −2/3 . If the dissipation by inelastic collisions is much larger than the viscous dissipation, then we obtain which also vanishes as St → ∞. In this case the expansion in the velocity field (7) proceeds in inverse powers of St −1/2 but the resulting expression for the source remains the same. For fixed e we see that the non-dimensional particle temperature will undergo a transition from a St −2/3 scaling at moderately large St to a St −1 scaling at very large St.
The drag coefficient R drag in (18) is determined to leading order from the average of F v,I . Since the average drag on a particle in a suspension where all the particles are moving with the same velocity U is the same as in the case of a fixed bed with an average velocity of the fluid −U , we note that R drag to leading order is the same as the drag coefficient R fb for fixed beds. Thus, for the case of elastic particles, The error of O(St −2/3 ) in the above expression arises from the contribution to the total drag from v III .
2.2. Results for energy source S * and mean drag R drag We have seen that to leading order the relative velocity v II of the particles is the same as that of a hard-sphere system. Thus, to determine S * (φ) we carried out simulations in which the evolution of particle positions was determined by a molecular dynamics code. The leading-order viscous force F α v,I is that associated with the gas motion produced by the mean velocity U of all the particles. However, this force varies with time (non-dimensionalized by a/T 1/2 ) due to the changes in particle configuration obtained from a hard-sphere dynamics code with a non-dimensional temperature of unity. At time t = 0, v α III is taken to be zero for all the particles. At subsequent times v III is determined by integrating (9) using Simpson's rule. At collision both v α III and v α II instantaneously change their values according to elastic collision rules. We found that the straightforward integration of (9) as mentioned above led to large fluctuations in v III ·F v,I . To reduce these statistical fluctuations, a small damping term −6πµa R drag v α III was therefore added to the right-hand side of (9). The simulations were carried out typically with = 0.03. Simulations for different values of in the range 0.01-0.05 were also carried out for selected values of φ, and the variations in were found to have no statistically significant effect on S * .
The forces on the particles were evaluated using the method described in  in which N particles are placed in a unit cell periodically repeated in space and the velocity disturbance due to each particle is expressed in terms of a series of multipoles. Table 1 shows results of simulations for three different values of volume fraction φ. N s indicates the largest order of the multipoles retained in the expansion. The total number of multipoles is then (3N 2 s − 1)N. We see that even though both S * and the variance of the viscous force, F v,I , change considerably with N s , the ratio of the two is approximately constant for a given φ. If we interpret S * as a product of the force variance and a correlation time, then the results of table 1 suggest that the correlation time is essentially independent of N s . Calculation of S * requires averaging over several thousand time steps in order to determine the correlation time with reasonable accuracy and the simulations with high N s are, as a consequence, computationally very intensive. We can, however, take advantage of the fact that the correlation time is relatively insensitive to N s and determine its value from simulations with lower values of N s . The force variance must be computed using higher N s , but its value can be determined from a relatively small number of independent hard-sphere configurations. Table 2 gives results for the force variance as a function of φ and N s , which were obtained by averaging over thirty independent hard-sphere configurations. Both the average drag force F v,I and the variance are sensitive to the volume fraction but the ratio of variance to average drag force squared is nearly constant at about 0.085 for φ greater than about 0.15.
The analysis of Koch (1990) for dilute suspensions gave S * = 1/(2π 1/2 ). Therefore we present the results for S * in figure 1 in terms of R s defined by R s represents the energy source due to a specified mean force acting on the particles whereas S * is the source for a specified mean relative velocity of the particles and gas. Note that R s → 1 as φ → 0. For φ > 0.45 we determined the fixed-bed drag coefficient R drag using N = 16 and N s = 6 while the correlation time was estimated using N s = 3 or 4 and N = 64. We found the correlation time at such high volume fractions to be approximately proportional to the average collision time for the hard-sphere systems. For example, the ratio of collision times at φ = 0.53 for simulations with N = 16 and N = 64 was 1.28 and the ratio of correlation times for the same values of N was about 1.25. The spatial configurations of hard-sphere systems change significantly at such high volume fractions with N, and this makes it necessary to use larger values of N at high volume fractions. On the other hand, Koch's analysis suggests that at very low volume fractions the correlation time is proportional to the time it takes for a pair of particles to separate from each other by the Brinkman length of O(a/φ 1/2 ). This requires that the unit cell be large compared with the Brinkman length, or that N be large enough. The results for φ = 0.03 and 0.05 were therefore obtained with, respectively, N = 256 and N = 128. The number of multipoles ((3N 2 s − 1)N) required in determining the forces on the particles is large for these situations, and we used the O(N) algorithm described in  to enhance the efficiency of the computations for φ 6 0.05 and φ > 0.45.
As shown in figure 1 the numerical results for R s can be represented well by the expression where χ is the value at contact, r = 2a, of the radial distribution function for a hard-sphere system which is given by Ma & Ahmadi (1988) as The rationale for fitting the numerical results according to (23) is as follows. At high volume fraction, the correlation time appears to be proportional to the average collision time, which, in turn, is expected to be inversely proportional to χ. The above expression for χ by Ma & Ahmadi suggests that this correlation time will vanish at φ = 0.64356 which is the maximum packing density for spheres in a random array. The first variation of R s from 1 at low volume fractions is expected to be O(φ 1/2 ) because of the dependence of the correlation time on Brinkman screening length. Results for the fixed-bed drag coefficient R drag were obtained previously using a multipole method by Ladd (1990) for φ 6 0.45. Our results are in very good agreement with his as shown in table 3. For higher volume fractions, the empirical correlation by Carman (1937) Table 3. R drag as a function of φ. The result for φ = 0.61 is taken from . Shown also are the corresponding results obtained by Ladd (1990), equation (26), and the Carman correlation, (23). is known to provide reasonably accurate estimates. Our results for φ > 0.4 are indeed in good agreement with Carman's expression.
For the stability analysis it will be necessary to have an analytical fit for the numerical results for R drag as well. We shall use The above expression agrees to O(φ) with the theoretical expression for the drag in dilute fixed beds (Howells 1974;Hinch 1977;Kim & Russel 1985) and with the results of numerical simulations as can be seen from figure 2. For φ > 0.4 we shall use the much simpler expression given by Carman (25) with a small addition of 0.7 to its right-hand side which renders the values of R drag given by the two expressions equal at φ = 0.4. Figure 2 shows a comparison between the numerical results and the above expressions for R drag . Finally, we shall also need R diss (φ) to compute the particle velocity variance. The viscous energy dissipation in a suspension with a Maxwellian particle velocity distribution was evaluated by  who obtained The above expression agrees that for R drag up to terms of O(φ ln φ). At O(φ) the two expressions differ since the O(φ) coefficient depends on the details of pair interactions and the velocities of the particles are different in the two problems. R diss is determined for Maxwell-distributed particle velocities with zero mean and R drag for particles with a common velocity relative to the gas. The energy dissipated in the gap between two particles having a non-zero relative velocity and separated by a gap width of increases as 1/ as → 0 if one assumes that the standard incompressible, continuum lubrication flow between smooth, rigid, spherical particles is applicable at all particle separations. The resulting R diss would be infinite.  allowed for the breakdown of the usual lubrication analysis for the particles with the gap width less than 2 m a. The resulting R diss therefore diverges logarithmically with m . For gas-solid suspensions the most likely cause of the lubrication breakdown is non-continuum effects in the gas since the mean free path λ of the gas molecules is typically about 1000Å. Sundararajakumar & Koch (1996) have determined the energy dissipation in the gap between colliding particles accounting for the finite mean free path of the molecules and their analysis indicates that using in expression (27) for R diss will yield the correct rate of dissipation. Gopinath, Chen & Koch (1997) have shown that m may also depend on the compressibility of the gas for larger particles and larger velocity differences between the particles. However, the qualitative behaviour of the suspension will not be very sensitive to the exact nature of lubrication breakdown since it only contributes a weak logarithmic factor.

Moderate-Stokes-number suspensions of slightly inelastic particles
In the preceding development, we have assumed that the velocity distribution function is controlled by the elastic collisions among the particles and is therefore an isotropic Maxwellian. This requires that the mean-free time between collisions be much smaller than the viscous relaxation time. In a dilute suspension, this occurs for very massive particles with St φ −3/2 . In the regime φ −3/4 St φ −3/2 , corresponding to moderately massive particles, the particle velocities do not relax to the local fluid velocity but the velocity distribution function is controlled by hydrodynamic interactions rather than collisions. Hydrodynamic interactions feed more energy into the vertical particle velocity fluctuations than the horizontal. As a consequence, the anisotropy in velocity fluctuations is quite significant in moderateto-low-Stokes-number suspensions. As in the case of  and Kang et al. (1997) we shall develop an approximate kinetic theory for the moderate-Stokes-number case using a Grad approximation for evaluating the collisional terms in the balance equation. The hydrodynamic source of particle fluctuating energy will be approximated as the anisotropic tensor derived by Koch (1990) modulated by a function of particle volume fraction computed from the high-Stokes-number simulations presented in the previous section.
To determine the anisotropy in the velocity variance we need balance equations for the second moments of the fluctuating velocity C α i = v α i − v i . An equation for the second moment T ij ≡ C α i C α j can be readily obtained from (1) and (2) and is given by Now we decompose the right-hand side into the source and sink tensors resulting from various mechanisms. For example, the sink due to viscous dissipation is expressed as  used the above sink term in their kinetic theory for sheared gassolid suspensions and found that the even though R diss was determined for high-Stokesnumber suspensions with isotropic Maxwellian distributions, the resulting theory predicted the velocity moments reasonably accurately even when Stokes number was not large. The source due to hydrodynamic interactions will be expressed as Koch (1990) derived the vertical (S * ) and horizontal (S * ⊥ ) components of the source tensor in the above expression for dilute sedimenting suspensions. We shall assume that the same expressions hold for non-dilute suspensions when his expressions are multiplied by the factor R s (φ)R 2 drag (φ) that accounts for higher-order particle interactions in high-Stokes-number suspensions, i.e. we assume where β 2 = T /(T − T ⊥ ) and T and T ⊥ are, respectively, the variance in the velocity components parallel and perpendicular to U . The collisional change of the second moment of velocity in (29) can be evaluated using a Grad approximation as was done by Jenkins & Richman (1985) and . Setting the velocity gradient terms in (4.21) of  to zero we obtain The trace of the above expression for nearly elastic particles (0 < 1 − e 1) is the same as the expression (6) for −Γ inelas given earlier. The terms inside the square brackets in (34) reduce to T δ ij − T ij when the coefficient of restitution is unity indicating the tendency of the collisions to drive the velocity distribution towards an isotropic Maxwellian.
Collecting all the terms, the balance equation for the second moments becomes In the limit of large St and small 1−e the collision term forces the velocity distribution to become isotropic (β → ∞), and the trace of the above equation reduces to (17) 10 -3  Koch (1990). His analysis gave Note that the anisotropy of the variance is very large in this limit. The variance can be estimated for arbitrary φ and St by keeping all the terms in (35) and solving the resulting equations numerically. The theoretical results for φ = 0.02 and e = 1 are shown in figures 3 and 4 where they are compared with the results obtained by the direct numerical simulations to be described in § 5. At low Stokes number, the simulated velocity variance reaches a constant value characteristic of an inertialess suspension, cf. figure 3. As the Stokes number is increased the Brinkman screening of the hydrodynamic interactions among the particles leads to a velocity variance that decreases in proportion to St −2/3 . In view of this high-Stokes-number scaling, we plot the products T St 2/3 and T ⊥ St 2/3 as a function of the Stokes number in figure 4. At low Stokes numbers, the collisions have relatively little effect on the velocity variance and the anisotropy of the hydrodynamicinteraction source leads to a large ratio of T /T ⊥ . As St increases, collisions become increasingly rapid compared with viscous relaxation and the velocity distribution becomes increasingly isotropic. Figure 5 illustrates the effects of inelastic interparticle collisions on the velocity variance. Results are presented for φ = 0.05 with coefficients of restitution e = 0.8, 0.9 and 1. The dashed lines indicate the high-Stokes-number asymptotes (18) and (20). It can be seen that the variance for perfectly elastic particles scales like St −2/3 over the full range of Stokes numbers, St > 10 shown. On the other hand, the inelastic particles exhibit a transition from a St −2/3 scaling at moderate St to a St −1 scaling at very high Stokes numbers. The circle indicates the velocity variance obtained in a dynamic simulation with 16 particles in a unit cell. The variance is approximately 40% higher than the theory. The additional variance observed in the dynamic simulation can be attributed to the inhomogeneous structure that develops in this unstable suspension. This will be discussed further in the forthcoming sections.

Stability analysis
In this section, we will introduce equations of motion applicable to an inhomogeneous fluidized bed and perform a linear stability analysis to determine the conditions under which a homogeneous bed is unstable to waves of particle concentration. Previous studies (Didwania & Homsy 1982;Anderson et al. 1995;Glasser et al. 1997;Goz 1995;Harris & Crighton 1994;Batchelor 1988Batchelor , 1993 of fluidized bed stability have indicated that waves with the wave vector parallel to gravity (vertical waves) have the largest growth rate. Thus, the initial instability involves vertical waves, although horizontal structure may set in subsequently. We will restrict ourselves to determining the form of the initial instability and will therefore consider only vertical waves in which the average particle velocity U = e 3 U is in the vertical direction and the dependent variables U, φ, T , and T are functions of the vertical coordinate z = e 3 · x and time t only. The average mixture velocity u will be taken to be zero.
The traditional particle-phase equations of motion for fluidized beds have included mass and momentum conservation equations. The momentum conservation equation included a particle pressure which was assumed to be a function of the volume fraction only. In our treatment based on the kinetic theory of granular materials, one finds that the particle pressure is related to a particle-phase temperature as well as the volume fraction. Thus, we require additional conservation equations for the kinetic energy associated with the particle velocity fluctuations. We found in § 3 that the anisotropy of the velocity variance is often quite significant and so we will consider separate conservation equations for the vertical velocity variance T and the overall variance 3T = T + 2T ⊥ .
Non-dimensionalizing U and T 1/2 by the terminal velocity U t , z with the particle radius a, and t by a/U t , the equations of mass and momentum conservation reduce to ∂φ ∂t The above equations are the same as those given in our treatment of a sheared gas-solid suspension (cf. (4.13)-(4.16) and (4.22) from  except that the momentum equation contains a body force, φ/St, due to gravity. In (38) P m = (φ + 8B/5)T + (12/5)BT , The left-hand side of (38) contains the particle-phase inertia terms. The second term on the right-hand side is the mean drag per unit volume exerted by the gas on the particles; this force is a function of volume fraction owing to the hydrodynamic interactions among the particles and the dependence of R drag on φ plays a crucial role in the instability. The third term is the divergence of a stress P m given by (39). This stress contains terms proportional to the particle velocity variances T and T . These variances are different from those obtained for a homogeneous suspension in the previous section and contain contributions proportional to ∂U/∂z (see (42) and (43) below). Thus, the third term on the right-hand size of (38) includes both pressure gradient and viscous stress terms. The fourth term on the right-hand side of (38) is an additional viscous stress contribution due to the collisional stress.
The solution of (39) and (35) for the vertical particle pressure in a homogeneous fluidized bed indicates that P m scales like St −2/3 in the limit of large Stokes number. The pressure is non-dimensionalized with ρ p U 2 t where ρ p is the mass density of a particle. Since U t is proportional to a 2 and St ∼ a 3 , the dimensional pressure is predicted to be proportional to the particle radius squared. This scaling of the particle pressure with particle radius (at a constant ratio of the superficial to the minimum fluidization velocity) has been observed in experimental studies by Cody et al. (1996) in non-bubbling gas fluidized beds. On the other hand, Campbell & Wang (1991) showed that the particle pressure in a bubbling gas-fluidized bed was well correlated with the bubble-diameter suggesting that the pressure is produced primarily by the bubble-induced flow in this inhomogeneous situation. Figure 6 illustrates the variation of P m St 2/3 with volume fraction for St = 10, 100, 1000, and ∞. In the dilute limit, the source and dissipation of energy approach constant values. Thus, the granular temperature is finite as φ → 0 and the pressure is proportional to φ. As the particle concentration increases, two competing factors influence the particle pressure. For a given particle temperature, the particles in a dense bed undergo more frequent collisions. However, the source of energy decreases and the viscous dissipation of energy increases with increasing volume fraction. The net result is that the particle pressure passes through a maximum at φ ≈ 0.5 and then diminishes as one approaches the close-packed limit. Zenit, Hunt & Brennen (1997) have observed a similar qualitative variation of the particle pressure with volume fraction in liquid-fluidized beds at large Stokes and Reynolds numbers. However, the absolute magnitude of the particle pressure was considerably larger than that predicted here for a homogeneous, low-Reynolds-number, gas-fluidized bed. The scaled vertical pressure P m St 2/3 is high at moderate Stokes numbers owing to the anisotropy of the source of particle fluctuation energy and eventually reaches a finite high-Stokesnumber asymptote corresponding to a state where interparticle collisions efficiently exchange energy between the vertical and horizontal particle motions. Figure 7 illustrates the effects of inelastic interparticle collisions on the granular pressure. The vertical pressure P m is plotted as a function of particle volume fraction for St = 100 and three coefficients of restitution e = 0.8, 0.9, and 1. The inelastic collisions lead to a modest decrease in the granular pressure that is most significant at intermediate volume fractions.
In the dilute limit, particle collisions are infrequent compared with interparticle hydrodynamic interactions (which occur at interparticle separations comparable with the Brinkman screening length), so the dissipation due to inelastic collisions is unimportant. At high volume fractions, the granular temperature becomes quite small. Since the viscous dissipation is proportional to T while the loss due to inelastic collision scales like T 3/2 , the small granular temperature makes viscous dissipation more important than the energy lost due to solid-body collisions.
Equations for the second moments of the velocity distribution in an inhomogeneous suspension can be obtained using (4.8), (4.10), (4.21), and (4.22) of Sangani et al.  1996). Thus, the conservation equations for the vertical velocity variance T and the overall variance 3T are Sangani et al. only presented results for a spatially homogeneous suspension. When we generalize their equation (4.20) to inhomogeneous time-dependent suspensions, we encounter terms involving third moments C 3 3 and C 3 C 2 associated with the fluxes of vertical and overall temperatures. To obtain closed equations of motion, we have assumed that the fluxes of vertical and total energy can be written as −( 1 3 )κ(∂T /∂z) and −κ(∂T /∂z), where κ = 8 π 1/2 φ 2 χT 1/2 1 + 25π 512φ 2 χ 2 (44) is the thermal conductivity for an isotropic gas with temperature T . We use the vertical temperature in (44) since we are interested in the vertical flux of energy. The left-hand sides of (42) and (43) represent the convection of vertical and overall fluctuation energy, respectively. The first terms on the right-hand sides are the divergences of the thermal fluxes. The second and third terms are the work done by the stress, i.e. pressure work and viscous dissipation, where The last two terms on the right-hand sides of (42) and (43) are the sources and dissipations of vertical and overall energy: with S * = S * + 2S * ⊥ . The effect of hydrodynamic interaction forces was purely dissipative in Sangani et al.'s analysis of sheared gas-solid suspensions in the absence of gravity and it led to the second terms on the right-hand sides of (46) and (47). However, the fluctuating force associated with the gas-velocity disturbance induced by the particles' gravitational settling creates sources of energy represented by the first term. The third term on the right-hand side of (47) is the dissipation of total energy due to inelastic collisions, whereas the corresponding term in (46) involves a collisional exchange of energy between vertical and horizontal fluctuations as well as an inelastic dissipation.
In the absence of spatial and temporal gradients of average properties, the equations of motion have a unique solution for velocity, temperature and T /T given by (18), (19), and the implicit result obtained by setting the right-hand side of (35) to zero. To assess the stability of the homogeneous state, we substitute a solution of the form into the equations of motion (37), (38), (42), and (43). Here, Ψ represents the dependent variables φ, U, T , and T , σ is the growth rate, k is the wavenumber, Ψ 0 is the base state, and is the amplitude of the perturbation. Linearizing the resulting equations for 1, we can solve for the relative amplitudes U/φ, T /φ, and T /φ and obtain a fourth-order polynomial dispersion equation for the growth rate σ(k). This dispersion equation is rather long and will not be presented here. We find that three of the solutions for σ remain negative for all values of k. The fourth root, which controls the stability starts from the origin, i.e. σ = 0 at k = 0.
The real part of this root is plotted as a function of k for perfectly elastic particles with a volume fraction of 0.5 and several values of the Stokes number in figure 8. Below a critical Stokes number of St c = 4.36, the growth rate is always negative indicating that perturbations of all wavelengths decay and the suspension is stable. At higher St, the growth rate first increases taking on positive values, then passes through a maximum, and eventually becomes negative. At the larger wavenumbers, the perturbations are damped primarily by the particle-phase viscosity and conductivity. Thus, above the critical Stokes number St c , an unbounded suspension is unstable and the instability involves perturbations with wavenumbers between 0 and some maximum unstable wavenumber k m .
At small volume fractions a slightly different behaviour is observed as illustrated in figure 9 for perfectly elastic particles with a volume fraction of 0.02. In this case, the growth rate is negative for all wavenumbers for St < 16. For a small range of intermediate Stokes numbers 17 < St < 21.2, the growth rate starts from zero at k = 0, becomes negative, passes through a minimum, becomes positive and finally passes through a maximum and becomes negative again. At St > 21.2, the behaviour of σ(k) is similar to that observed above the critical Stokes number in the dense suspension φ = 0.5 discussed above. A simple criterion for the marginal stability of an unbounded suspension can be obtained by considering the small-wavenumber behaviour of σ(k). Near k = 0, the growth rate can be expanded in a Taylor series of the form: where σ 1 and σ 2 are purely real. The first term ikσ 1 represents a propagation of the wave with wave speed σ 1 and the stability of the suspension in the long-wavelength limit is determined by the sign of σ 2 . A positive value of σ 2 provides a sufficient but not a necessary condition for instability of the unbounded suspension. For dense suspensions, the Stokes number at which σ 2 becomes positive corresponds exactly to the marginal stability limit. In dilute suspensions, there is a narrow range of Stokes numbers for which σ 2 < 0 but the growth rate becomes positive at higher wavenumbers. Nonetheless, we will adopt the condition σ 2 = 0 as an approximate criterion for the marginal stability of the homogeneous suspension. In the case φ = 0.02, this leads to an estimate St c = 21.2 whereas the actual critical Stokes number is between 16 and 17. Thus, the critical Stokes number will be estimated by substituting (49) into the full dispersion relation and setting terms of O(k) and O(k 2 ) to zero. This procedure indicates that the thermal convection and conduction terms have no effect on σ 1 and σ 2 and the energy equations reduce to the simple algebraic equations: In other words, the time for thermal convection or conduction over large wavelengths is sufficiently long that a local balance of sources and sinks of energy controls the particle-phase temperatures. The linearized momentum and mass conservation equations take the form where α = σ + ikU 0 and the partial derivatives of pressure in (52) can be evaluated using and the partial derivatives of T = T (φ, U) and T = T (φ, U) can be evaluated using the algebraic energy balances (50). The dispersion relationship takes the form Substituting (49) into (55), we obtain a kinematic wave speed: and the growth rate for large wavelengths: To confirm that this local energy balance formulation captured the correct long-wavelength behaviour, we performed a small-k asymptotic analysis on the dispersion equation resulting from the full problem and obtained the same results for σ 1 and σ 2 . Since all suspensions with σ 2 > 0 are unstable and most cases with σ 2 < 0 are stable at all wavenumbers, (57) provides a condition that can be used to estimate the critical Stokes number above which the suspension is unstable. The first term on the right-hand side of (57) results from the particle inertia and the dependence of the drag on volume fraction. This is the same destabilizing term found in previous stability analyses (Jackson 1963). The second term indicates that the increase of pressure with increasing volume fraction tends to stabilize the suspension by driving a flow from dense regions to dilute regions of the bed. This is similar to the stabilizing mechanism identified previously (Anderson & Jackson 1968). In previous analyses, the pressure was considered to be a function of the volume fraction alone and no energy balance was considered. In our analysis, the partial derivative of pressure with respect to volume fraction involves both a direct dependence and a contribution from the change in the particle temperature as determined by the local energy balance. An additional destabilizing term results from the dependence of the pressure on velocity. An increase in the mean velocity of the particles increases the source of energy due to hydrodynamic interactions and thereby increases the particle temperature and pressure. This tends to drive particles from dilute regions with fast moving particles into dense regions with slow moving particles. This effect is not incorporated in  Figure 10. The leading real contribution to the growth rate in the limit of small wavenumber, σ 2 , plotted as a function of the particle volume fraction is indicated by the solid curve. The dashed, dash-dot and dash-dot-dot curves indicate the first, second and third contributions defined in (57). The Stokes number is 50 and the coefficient of restitution is 1.
previous studies which assumed no dependence of the particle pressure on relative velocity. The contributions of these three terms to σ 2 are illustrated in figure 10 for perfectly elastic particles with St = 50 and a range of volume fractions. The destabilizing term resulting from the dependence of R drag on volume fraction (term 1, dashed line) is zero at φ = 0, grows rapidly with φ, passes through a maximum at φ = .046, and decreases at higher volume fractions.
The destabilizing term due to the dependence of particle pressure on U (term 3, dash-dot-dot line) shows a similar behaviour but with a smaller amplitude. Thus, the traditional destabilizing term is still the most significant one. The resistance of the suspension to compression, i.e. −∂P m /∂φ (term 2, dash-dot line) is negative and provides a small window of stability at low volume fractions φ < 0.0125. However, this stabilizing effect is diminished at higher volume fractions. In fact, −∂P m /∂φ actually becomes positive for φ > 0.575. As illustrated previously in figures 6 and 7, the particle pressure decreases with increasing φ in very dense beds. This occurs because the increased viscous dissipation and decreased source of fluctuation energy lead to a diminution of the particle temperature as one approaches the close-packed limit. The destabilizing effect of this negative compressibility leads to a 6.5-fold increase in σ 2 between φ = 0.575 and 0.64 (see insert of figure 10).
The critical Stokes number above which the homogeneous unbounded suspension is unstable is plotted as a function of volume fraction for several values of the coefficient of restitution in figure 11. In the dilute limit, the critical Stokes number is 0.149φ −3/2 . However, this dilute asymptote is only quantitatively accurate (relative errors of less than 5%) for volume fractions less than 0.001. For values of φ greater than about 0.05, the critical Stokes number for an unbounded suspension is sufficiently small, St c < 5, that the high-Stokes-number theory is not applicable. Thus, except in the dilute limit, an unbounded, high-Stokes-number suspension undergoing hydrodynamic interactions and short-duration solid-body collisions is unstable. This observation is consistent with the experimental observation of Tsinotides & Jackson (1993) that frictional forces are associated with the interval of stable fluidization in the gas-fluidized beds examined in their experimental studies. It is interesting to note that the critical Stokes number for elastic particles passes through a minimum of 1.41 at φ = 0.237 and a maximum of 4.98 at φ = 0.479 before diminishing once again. Over a small range of volume fractions, 0.518 < φ < 0.52, the theory predicts an unstable suspension at low Stokes numbers followed by a narrow window of stability and an unstable suspension at high Stokes numbers. Above φ = 0.52, the suspension is unstable at all Stokes numbers. However, it should be noted that the details of the behaviour of St c (φ) in dense suspensions may not have physical significance because the high-Stokes-number theory is not applicable at the low-to-moderate Stokes numbers at which the marginal stability limit occurs. Therefore, the main conclusion that can be drawn is that high-Stokesnumber suspensions are unstable.
Most analyses of gas-fluidized beds, other gas-solid suspensions, and granular flows assume that the particle velocity fluctuations are nearly isotropic. In the present theory, we have adopted separate energy balances for the vertical and horizontal velocity fluctuations with an energy exchange between these modes resulting from solid-body collisions. We found that the velocity variance exhibits substantial anisotropy owing to the larger source of fluctuating energy in the vertical direction, cf. figure 4. It is interesting then to examine the effect that the anisotropy of the velocity variance has on the stability predictions by comparing the results for the full analysis with those for a simplified analysis based on a single overall energy balance. This comparison is given in figure 12 where the critical Stokes number is plotted as a function of φ. The theory based on an isotropic velocity variance gives St c ∼ 0.0788/φ 3/2 in the limit φ → 0 in agreement with the asymptotic analysis of Koch (1990). This dilute suspension behaviour is similar to that obtained from the full anisotropic variance calculations but St c is decreased by a factor of about 2. The anisotropic theory takes into account the fact that most of the energy is contained in vertical velocity  Figure 12. The approximate critical Stokes number for perfectly elastic particles obtained from the full theory (upper solid curve) is compared with that obtained from a model with an isotropic granular temperature (lower solid curve). The dashed line is the dilute asymptote 0.0788 St/φ 3/2 for an isotropic suspension derived by Koch (1990).
fluctuations and it is these fluctuations that tend to stabilize the suspension to vertical waves. The critical Stokes number obtained for an isotropic granular temperature goes to zero at φ = 0.0794, whereas the anisotropic-variance analysis gives a small but finite St c in more concentrated suspensions. When numerical simulations are performed in a unit cell with periodic boundary conditions, this will restrict the wavelength of the disturbances that can affect the system. The stability analysis predicts a minimum wavelength (maximum wavenumber) for instability at any St > St c as illustrated in figures 8 and 9. In practice, we have found that a periodic unit cell that is large enough to contain enough particles to approximate a continuum will be large enough to admit unstable modes except at Stokes numbers that are very close to the critical value. For example, for φ = 0.02 and St = 22 (close to the critical Stokes number of 16-17), the maximum unstable wavenumber k = 0.42 corresponds to a cubic unit cell containing 15 particles.
The plots of growth rate as a function of wavenumber in figures 8 and 9 show that the growth rate passes through a maximum at an intermediate value of the wavenumber. In figure 13, we present the maximum growth rate σ m as a function of volume fraction for perfectly elastic particles at St = 5, 10, 50, and 100. In very dilute suspensions, the maximum growth rate is zero and the suspension is stable. As the volume fraction increases, σ m grows, passes through a maximum at φ ≈ 0.1 and decreases again. The maximum growth rate becomes quite small (σ m = O(10 −3 )) for φ ≈ 0.45-0.55. This suggests the possibility that the instability of a fluidized bed could be missed in numerical integrations of the equations of motion if insufficient time is allowed for the perturbations to grow. However, it should be noted that the dimensional growth rate remains sufficiently high that the instability would be readily observed in experiments. For example, the dimensional growth rate for the most dangerous mode for 60 µm diameter particles with density of 1 g cm −3 at atmospheric conditions would always be greater than 1.4 s −1 . The dimensional wave speed under these conditions is about 0.8 cm s −1 , so that the disturbance would only propagate about 1 cm while increasing in amplitude by a factor of e. Comparing figures 10 and 13, it may be seen that the qualitative dependence of the maximum growth rate on volume fraction mirrors the behaviour of the growth rate in the long-wavelength limit. The wavelength of the fastest growing wave non-dimensionalized by the particle radius is plotted as a function of φ for St = 100 in figure 14. The wavelength is very large in dilute suspensions and approaches infinity at the critical volume fraction for instability. At moderate volume fractions, the wavelength is about 10-15 times the particle radius. This value is somewhat smaller than the wavelength (about 40a) of the void fraction waves observed in liquid-fluidized beds by Didwania & Homsy (1981). The maximum growth rate increases rapidly with volume fraction for very dense suspensions, 0.57 < φ < 0.64, cf. figure 13. This sharp increase coincides with the regime in which the particle pressure decreases with increasing volume fraction so that the suspension is mechanically unstable even in the absence of a volume-fractiondependent mean drag coefficient. It should be noted that the high growth rates predicted for φ > 0.57 may not be obtained in practice because of the inapplicability of the assumed form of the equations of motion at these high concentrations. Dense suspensions may exhibit a solid-like behaviour when there is insufficient free volume to allow the particles to translate relative to one another. Deviations of the behaviour of dynamical simulations of hard spheres from the predictions of the kinetic theory for dense gases are typically observed for φ > 0.55 (Woodcock 1981).

Dynamic simulations
In the previous sections, we presented equations of motion for a gas-fluidized bed and examined their linear stability. The source of particle fluctuation energy arising due to the particles' motion relative to the gas was computed using simulations of interparticle hydrodynamic interactions that approximated the particle motions as those in a homogeneous hard-sphere system. In the present section, we will discuss results of dynamic simulations obtained by solving the full equations of motion for the particles and the viscous gas. These simulations will be used to test certain aspects of the theoretical predictions of particle temperature and the instability of a homogeneous gas-fluidized bed. The simulations for systems with fewer than N = 64 particles were conducted using the method described in  and  while those for larger systems were based on the O(N) algorithm described in  with a few modifications to be specified later in this section. All of the simulations used a lubrication cut-off of m = 0.01 to account for the finite mean-free path of the interstitial gas.
Results for the velocity variance as a function of St in a dilute suspension with φ = 0.02 and N = 16 were presented earlier in figure 3. A more detailed account of these simulations is given in table 4. The initial condition for each simulation was specified as a hard-sphere spatial distribution with each particle's velocity equal to the steady-state mean sedimentation velocity. The equations of motion for the hydrodynamically interacting particles were then solved for several thousand Stokes settling times a/U t to reach a statistical steady state. Thereafter, the simulations were run for the times specified in table 4 and time-average properties of the suspension were computed.
Two measures of the spatial configuration of the particles were monitored. The value of the radial distribution function g(r) at r = 2.1a was used to quantify short-range clustering tendencies. As noted previously (Ladd 1996;Koch & Shaqfeh 1991), the probability of close pairs in a dilute sedimenting suspension without particle inertia is higher than that in a hard-sphere fluid. As the Stokes number is increased to St = O(10), particle inertia allows particles to come closer to one another and the tendency for small-scale clustering increases. At higher Stokes numbers, however, particles have sufficient inertia to bounce off one another and separate after a collision and g(2.1a) decreases and approaches the hard-sphere value. The small-scale clustering at moderate St results in an increase in the mean sedimentation velocity because two spheres in close proximity fall faster than a single sphere. As the Stokes number is increased further, the mean velocity passes through a maximum and decreases as a result of the decreased probability of close pairs. The stronger (fixed-bed-like) hydrodynamic interactions between high-Stokes-number particles also decrease the mean velocity. For St > 500, the simulated mean velocities are comparable with the value 0.683U t obtained for a fixed bed, cf. (19).
The second measure of spatial configurations is the structure factor defined by This is used to quantify microstructural changes occurring on the lengthscale of the periodic box used in the simulations. In particular, since the instability is predicted to occur first for the largest wavelength possible, the structure factors corresponding to the smallest vertical and horizontal wavenumbers (largest wavelengths) that fit in the box are reported in table 4. For hard-sphere systems with φ = 0.02 and N = 16 these structure factors are approximately equal to 0.85. We see that the structure factors for finite-Stokes-number suspensions are comparable in magnitude: variations seen in table 4 are statistically insignificant. The stability analysis predicts that the homogeneous suspension is unstable for the cases with St > 17 listed in table 4. However, the box length is smaller than the mean-free path for the small number of particles (16) used in the simulations and the continuum description of volume fraction variations is inapplicable. These small boxes provide a means of approximating the behaviour of a homogeneous suspension and it was shown in figures 3 and 4 that the velocity variance obtained in the dynamic simulations was in agreement with the theory. To find evidence of the inhomogeneity induced by the predicted instability, we performed a series of simulations (summarized in table 5a) with St = 50, φ = 0.05 and a range of simulation cell sizes N = 16-512. The theory predicts that the homogeneous suspension is unstable for these conditions. The structure factor for a homogeneous hard-sphere distribution is about 0.7 and is nearly independent of N. The structure factors obtained from the dynamic simulations grow with increasing box size and are consistently higher than those for a hard-sphere distribution. This is illustrated in figure 15(a) where the vertical and horizontal structure factors are plotted as circles and squares, respectively, and the structure factor for a hard-sphere distribution is indicated by the triangles. In the largest unit cells N = 512, the structure factor is 2-4. The structure factor may be expected to grow with N if the structure consists of a macroscopic inhomogeneity rather than local particle-scale clustering. The inhomogeneous structure leads to a spatially varying gravitational  force that drives velocity fluctuations in the suspension. These are seen in the form of a velocity variance in the vertical direction that grows from 0.027 for N = 16 to 0.084 for N = 512. For comparison, the velocity variance predicted for a homogeneous suspension is 0.0129. It is likely that the enhanced particle pressure created by the shear work associated with this macroscopic shearing motion limits the extent of the suspension structure formed. The mean velocity of the particles also grows with box size. This may be attributed to the fact that more particles will be found in dense downward flowing portions of the suspension than in less dense regions with upward mean mixture velocities.
The simulations discussed above were performed with the multipole order for individual particles corresponding to N s = 2. Large-N simulations using the O(N) algorithm described in  also require a specification of N sp , the multipole order used in representing the far-field velocity induced by groups of particles. The results presented in table 5 for N > 64 employed N sp = 3. Computations were also made with N sp = 4 for somewhat smaller time duration to assess the effect of varying N sp , and it was found that N sp = 3 and N sp = 4 gave essentially the same results.
Simulations for a more dense suspension (φ = 0.35 and St = 100) that is predicted to exhibit inhomogeneous structure are summarized in table 5(b). Since the source of velocity fluctuations at φ = 0.35 for N s = 2 is much smaller than for N s = 3, these simulations were made with N s = 3. At higher volume fractions the error introduced by grouping the particles in the O(N) algorithm is significant unless sufficiently large N sp is used. For simulations with N = 1024 this required that N sp > 5 be chosen to minimize the error. The computing time increases roughly as N 4 sp , and therefore N sp = 5 requires about 35 times more computing time than the value N sp = 2 used in . In that study, the same N sp was used for groups of different sizes. We found that the computing time can be reduced significantly without sacrificing accuracy if a smaller N sp was used for smaller groups of particles and a larger N sp for larger groups. Thus, the results for N = 1024 were obtained by using N sp of 4 for level 2 particle groups and N sp of 5 for level 1 groups.
The results for φ = 0.35 are similar to those for φ = 0.05. The box-length structure factor reaches a value of 0.4-1.3 at N = 1024 that is much higher than the value (0.1) for a hard-sphere distribution. This is illustrated in figure 15(b) where the parallel and perpendicular structure factors for the suspension are plotted as circles and squares respectively and the hard-sphere structure factor is given by the triangles. An alternative way of visualizing the suspension structure is to divide the simulation cell for N = 1024 into 64 boxes so that on average each box contains 16 particles. The variation of the mean volume fraction among these boxes then gives an indication of the inhomogeneity of the suspension and the tendency to form regions devoid of particles which may represent the first stages of bubble formation. The minimum and maximum volume fractions φ min and φ max among the boxes is plotted as a function of time in figure 16 for the gas-fluidized bed (lines) and a hard-sphere system with the same T (plusses and crosses). The smallest φ min attained in the gas-solid suspension correspond to 6 particles in a box whereas the most dilute regions for the hard-sphere system contained 10 particles. As in the more dilute suspensions, the velocity variance at φ = 0.35 increases with increasing system size. The variance grows from a value T = 6.3 × 10 −4 at N = 16 that is 40% higher than that for a homogeneous suspension to a value T = 1.5 × 10 −3 for N = 1024 that is more than three times larger than the homogeneous value.

Conclusions
In this paper, we have derived equations of motion for gas-solid suspensions with a mean relative motion between the gas and particulate phases. The equations presented in § 4 apply when all the spatial gradients are parallel to gravity. A more complete set of equations of motion allowing for arbitrary spatial gradients is summarized in the Appendix. The theory of rapid granular flows as derived, for example, by Lun et al. (1984) describes the flow of particles in a vacuum by incorporating inelastic particle collisions in the theory of dense gases.  extended this work to sheared gas-solid suspensions by determining the viscous dissipation of energy. Here, we have generalized the work of Sangani et al. to allow for a mean relative motion of the two phases. This relative motion gives rise to a mean drag force which, for a high-Stokes-number suspension, is equivalent to that acting on a fixed bed of particles, cf. figure 2 and (26). In addition, the fluctuating gas velocity field produced by the mean motion of the randomly positioned particles exerts fluctuating forces on the particles that result in a source of particulate-phase fluctuating energy given by (22) and (23) and figure 1. The present theory like those of Lun et al. and Sangani et al. assumes that the particles undergo instantaneous solid-body collisions but do not experience enduring solid-body contacts.
As applications of these equations of motion, we calculated the particle velocity variance (or granular temperature) and particle pressure in a homogeneous sedimenting gas-solid suspension or gas-fluidized bed and determined the marginal stability limits of the homogeneous state. The theoretical results were compared with dynamical simulations that solved the detailed equations of motion for an array of N = 16-1024 particles in a viscous gas. At moderate-to-high Stokes numbers, the velocity variance non-dimensionalized by the square of the particles' terminal velocity decays like St −2/3 with increasing Stokes number as illustrated in figure 3. The sedimenting particles produce stronger velocity fluctuations in the vertical than the horizontal direction and this makes the source of fluctuating energy highly anisotropic. As a result, the vertical velocity variance is significantly higher than the horizontal except at very high Stokes numbers where the solid-body collisions can efficiently transfer energy from vertical to horizontal motions, cf. figure 4.
The particle pressure (presented in figures 6 and 7) grows at first with increasing volume fraction due to the increasing particle concentration and the higher frequency of interparticle collisions. However, this growth is slower than would be obtained for a constant-temperature molecular gas because the temperature in the gas-solid suspension decreases with increasing volume fraction. This decreased temperature arises because of the stronger viscous dissipation occurring in a dense suspension and because the hydrodynamic source of energy decreases with increasing volume fraction. In very dense suspensions, φ > 0.57, the particle pressure actually decreases with increasing volume fraction.
Prediction of the marginal stability limits for the homogeneous state of a gasfluidized bed has long been considered an important test for proposed equations of motion for particulate flows. The standard theory identifies the coupled effects of particle inertia and the dependence of the drag coefficient on volume fraction as the destabilizing mechanism, while the particle-phase pressure was supposed to stabilize the homogeneous suspension under certain conditions. However, up to the present time, there has been no rigorous derivation of the magnitude of the particlephase pressure and so a predictive theory for the marginal stability limits could not be obtained. We have carried this line of theoretical investigation to its logical conclusion by using kinetic theory and Stokes flow simulation methods to derive the particlephase pressure and rheology in a gas-solid suspension over a wide range of particle volume fractions and Stokes numbers. An important part of this derivation is the determination of the particle-phase velocity variance tensor including hydrodynamic sources and sinks of particle-phase fluctuating energy. Thus, the equations of motion in the present investigation differ from those in much of the previous fluidizedbed literature in that they incorporate explicit energy equations for the vertical and horizontal particle velocity fluctuations. The linear stability analysis of the resulting equations of motion indicates that the homogeneous state is unstable to vertical waves for Stokes numbers greater than about 5 except in very dilute suspensions, where the critical Stokes number is proportional to φ −3/2 , cf. figure 11. This prediction is corroborated by our dynamic simulations which exhibit substantial deviations from the hard-sphere structure factor whenever the periodic unit cell is sufficiently large to admit waves that can be described by the continuum equations of motion for the suspension. The simulations also exhibit a velocity variance that grows with the simulation cell size; this variance may be attributed to the macroscopic shearing motions induced by the inhomogeneous particle distribution. Therefore, we conclude that the homogeneous state of a high-Stokes-number, dense, gas-fluidized bed is always unstable under circumstances in which the suspension may be described as a 'fluid' obeying equations of motion that can be derived on the basis of a microscale description of particles interacting via instantaneous hard-sphere collisions and hydrodynamic interactions.
This finding is consistent with the conclusion drawn by Tsinotides & Jackson (1993) based on their experimental observations that the homogeneous state of many (and possibly all) dense gas-fluidized beds is associated with a solid-like state of the suspension. To gain a better understanding of this phenomenon, it would be desirable to formulate a theoretical description of the transition from the solid-like behaviour that must occur in the dense-packed limit and the fluid state described in the present work.
It may be difficult to obtain a direct experimental test of our theoretical predictions beyond the observation that high-Stokes-number, dense, fluid-like gas-fluidized beds are unstable. A typical system of 70 µm diameter particles with a density of 2.5 g cm −3 has a Stokes number of 400. Thus, a homogeneous suspension of these particles would only be expected to be stable if it were very dilute φ < 0.004. To observe this marginal stability limit, great care would need to be taken to ensure that the suspension is initially well mixed so that any observed inhomogeneities can be attributed exclusively to inherent instabilities of the suspension. The marginal stability limits for high volume fractions are not of practical relevance to terrestrial gas-solid suspensions, because any particles small enough to exhibit the small Stokes numbers required would also be affected by electrostatic, van der Waals, and/or Brownian forces. It may be possible to produce a moderate-Stokes-number gas-solid suspension interacting via collisional and hydrodynamic forces by using large particles under microgravity conditions. A less exotic system that can exhibit Stokes numbers in the range 1-10 is a suspension of particles fluidized by a viscous liquid. The present theory is not directly applicable to this situation because it neglects the inertia of the fluid (in a liquid-solid system the Stokes and Reynolds numbers are comparable in value). In addition, the quantitative accuracy of the high-Stokes-number theory is questionable for St = 1-10. Nonetheless, the present work may elucidate some of the relevant physical considerations governing the stability of such a system.
The authors gratefully acknowledge the support of the National Aeronautics and Space Administration through grant NAG3-1853 and the National Science Foundation through grants CTS-9607723 (to A. S. S.) and CTS-9526149 (to D. L. K.). The computations were performed using the facilities of the supercomputing centres at Cornell University and University of Illinois at Urbana-Champaign.

Appendix. Equations of motion for a sheared gas-solid suspension with relative motion between the phases
In § 4, we presented equations of motion for a gas-solid suspension obtained by generalizing the analysis of  to allow for a relative motion of the gas and solid phases. Because we were interested in the behaviour of particle volume fraction waves with spatial variations parallel to the direction of the mean flow, we simplified the equations for this special case. In this Appendix, the full equations of motion for a gas-solid suspension with relative velocity between the phases and spatial gradients of velocity and volume fraction in arbitrary directions will be summarized. These consist of mass and momentum conservation equations for the gas and solid phases and an equation for the second moment of the particle velocity.
The mass conservation equations for the particle and gas phases are ∂φ ∂t + ∇ · (U φ) = 0 (A 1) where V is the mean velocity of the gas phase. Alternatively, one can consider the mixture average velocity W = V (1 − φ) + φU as a dependent variable in place of V . The incompressibility of the mixture implies that ∇ · W = 0. (A 3) divergence of the third velocity moment arising in the equation for the second moment must be approximated to obtained a closed set of equations. We have adopted the assumption that the conductive transport of each second moment may be expressed using Fourier's law for a hard-sphere gas. In the application in § 4 to void fraction waves we used T in the expression for κ based on the insight that the gradients of granular temperature and the relevant flux of granular energy were directed parallel to gravity. In the more general situation, it may be more appropriate to approximate κ by using the isotropic temperature T (one third of the velocity variance) in (44).