Particle Image Velocimetry Experiment and Discrete-Element

Jan 20, 2007 - Voidage Wave Instability in a Vibrated Liquid-Fluidized Bed ... waves that consist of alternating regions of high and low solid concent...
0 downloads 0 Views 2MB Size
Ind. Eng. Chem. Res. 2007, 46, 1375-1389

1375

Particle Image Velocimetry Experiment and Discrete-Element Simulation of Voidage Wave Instability in a Vibrated Liquid-Fluidized Bed Eldin Wee Chuan Lim,† Yee Sun Wong,‡ and Chi-Hwa Wang*,†,‡ Department of Chemical and Biomolecular Engineering, National UniVersity of Singapore, 4 Engineering DriVe 4, Singapore 117576, and Singapore-MIT Alliance, 4 Engineering DriVe 3, Singapore 117576

The nature of one-dimensional voidage waves in a liquid-fluidized bed that has been subjected to external perturbations and is exhibiting instabilities has been investigated both experimentally and numerically. Voidage waves that consist of alternating regions of high and low solid concentrations were observed to form and travel in a coherent manner along the fluidized bed. Solid particles were observed to move upward when a dense phase of the wave passed through their positions; otherwise, the particles settled downward. The voidage waves formed as a result of instability in such liquid-fluidized-bed systems are traveling waves with dense and dilute phases being convected along the bed. However, the motion of individual particles was observed to be highly restricted to a small region in space over short time scales. A diffusive type of behavior was observed where particles drifted gradually away from their initial positions within the bed. This type of motion was adequately described by a simple dispersion model used in the present study. Introduction Liquid-solid fluidization systems have evolved from its early industrial application as a way of separating solids of different sizes in the mining industry to its present day usage in hydrometallurgy, food technology, biochemical processing, and water treatment. In a liquid-fluidized bed, a suspension of solid particles is supported by a fluid flowing upward through a porous base plate.1 At low fluid flow rates, the bed is quiescent and the fluid flows through the interstitial spaces between the particles. When the fluid flow is increased, the drag force that is acting on the particles increases until it exactly balances their net gravitational force. The so-called minimum fluidization then occurs at the point where the particles are totally supported by the fluid. Further increase in fluid flow can lead to bed expansion and subsequently voidage instabilities over the entire bed. In practice, fluidized beds show different regimes of hydrodynamic behavior above the minimum fluidization velocity. The stable, uniform, and homogeneous state of fluidization is seldom found in practical applications. For example, gas-fluidized beds are usually very unstable with excess gas rising in the form of bubbles.2 On the other hand, it is well-established that liquidfluidized beds exhibit more-stable fluidization behaviors but may develop a type of instability known as voidage-wave instability under certain conditions. Such instabilities are of considerable interest to many researchers and have been studied using both experimental3-5 and theoretical approaches.6-11 Voidage waves were observed to be one-dimensional in narrow beds.4,12-14 In such cases, concentration plane waves propagated along the bed in the form of alternating dense and dilute layers. However, in wider beds, transverse structures were observed to develop as a result of secondary instabilities.15 Anderson et al.16 and Batchelor17 conjectured that bubbles originated from such secondary instabilities. Bubbles are important features in fluidized beds for inducing solid motion and promoting mixing. However, to date, the origin * To whom correspondence should be addressed. Tel.: 65-65165079. Fax: 65-6779-1936. E-mail: [email protected]. † Department of Chemical and Biomolecular Engineering, National University of Singapore. ‡ Singapore-MIT Alliance.

of bubbles is still not well-understood. As mentioned previously, they are conjectured to evolve from secondary instabilities, which, in turn, are believed to appear due to primary voidagewave instabilities.6,7,17 As such, investigations on such voidage waves are essential to providing us with a fundamental understanding of bubble formation. Therefore, the focus of the present study is on the first type of voidage instability observed in liquid-fluidized beds. A combination of experimental study using particle image velocimetry (PIV) and numerical investigation using computational fluid dynamics (CFD) and the discrete element method (DEM) is applied. The objective is to verify the usefulness of such a numerical approach in reproducing computationally the important features in this type of phenomenon. Experimental Section The fluidized-bed apparatus in this study consists of four major components: a vertical cylindrical column, rotameters, a centrifugal pump, and a liquid tank. The experimental setup is shown schematically in Figure 1a. The cylindrical column is a glass tube with a diameter of 2 cm and a height of 1 m. The bed is fluidized by water under ambient conditions. During the experiments, the bed is clamped in position to maintain its vertical and straight orientation. A centrifugal pump is used to circulate the water from the liquid tank through rotameters to the bed in a closed loop. The bed suspension is supported by a specially designed piston, which consists of a 2-mm-thick sintered stainless steel plate with 50-µm orifice holes. These small orifice holes ensure a large pressure drop and uniform distribution of flow.12 The piston can be set to oscillate vertically at a given frequency and amplitude, driven by a stepper motor. Such a pistondistributor setup would allow us to examine the response of the bed toward external perturbations. The amplitude of the piston was carefully chosen as 1.5 times greater than the particle diameter, to avoid close packing in the near-distributor region.4 The ratio of bed to particle diameters should be greater than 10, but should not exceed 25, to ensure smooth propagation of one-dimensional voidage-wave structures along the bed.4,12 In the present study, the granular material used was glass beads

10.1021/ie060864e CCC: $37.00 © 2007 American Chemical Society Published on Web 01/20/2007

1376

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

piston. To overcome this problem, the PIV system was operated under external trigger mode with a sampling frequency of 5 Hz. The motions of particles in one full oscillation cycle of the base could then be extracted from the captured images. Numerical Simulation

Figure 1. (a) Schematic diagram of the liquid-fluidized bed setup. Legend is as follows: 1, vertical cylindrical bed; 2, piston-like distributor; 3, rotameters; 4, centrifugal pump; and 5, liquid tank. (b) Schematic diagram of the velocity data acquisition system (PIV system). Legend is as follows: 1, test section; 2, PIV camera; 3, New Wave Nd:YAG laser; 4, TSI synchronizer; and 5, computer for data post-processing.

with a density of 2500 kg/m3 and a diameter of 1 mm (umf ≈ 0.01 m/s). The corresponding ratio of bed size to particle diameter then was 20. The column was filled with glass beads up to a height of 12 cm and fluidized by water at flow velocities in the range of 0.018-0.030 m/s. External harmonic perturbations were introduced into the bed by oscillating the piston distributor with an amplitude of 1.5 mm and frequencies in the range of 1-2 Hz. Typically, four regimes could be observed during experiments: packed, worming, planar wave, and turbulent regimes. The scope of the present study is limited to the planar wave regime. Solid concentrations were determined using a light scattering method at various bed levels. The basic principle of this technique is that under column backlighting conditions, the light intensity transmitted through the bed suspension varies strongly as a function of solid concentration. Here, a stabilized He-Ne laser (25 mW) was used as the light source and a low-energy photodiode (Newport, USA) was used to detect the transmitted intensity signals. The velocity field data of the bed was obtained using a PowerView PIV system from TSI Company (see Figure 1b). The column was placed inside a square glass tube with dimensions of 2.7 cm × 2.7 cm, filled with water to above the height of the bed within the column, to facilitate more-accurate measurements using the PIV system. New-Wave Nd:YAG lasers operating at ∼15 mJ/pulse were introduced from the side wall of the square tube, to illuminate particles near the central region of the column. Images of the illuminated particles were acquired by a TSI PowerView camera (resolution of 2K × 2K) at an angle perpendicular to the laser sheet. The PIV sampling frequency does not match the vibrating frequency of the moving

In the present study, the nature of voidage-wave instability was also investigated computationally using the DEM methodology coupled with CFD. This combined CFD-DEM approach has been used by many research workers recently for numerical studies of various types of solid-fluid multiphase systems. Feng et al.18 applied such an approach for the study of segregation and mixing of binary mixtures of particles in a gas-fluidized bed. Pandit et al.19 investigated the fluidization behavior of Geldart Group A and B powders in the presence of an externally imposed cohesive interparticle force. Following the experimental setup, the geometry of the simulated fluidization system consisted of a pseudo-threedimensional narrow channel 2 cm wide and 1.0 mm thick, containing 2500 glass beads as the solid phase and water as the interstitial fluid. It was deemed justifiable to perform a pseudothree-dimensional simulation of the fluidization system as the physical structures of interest in this case are one-dimensional traveling voidage waves. Feng et al.18 applied periodic boundary conditions to the front and rear walls of a rectangular fluidized bed with two side walls in their numerical study of particle mixing and segregation behaviors. This was to reduce wall effects and allow three-dimensional motion of particles without the use of a large number of particles. This aspect of the simulation was considered to be important to the particular phenomenon being investigated, because it has been previously reported that the packing structure of a bed of particles is strongly dependent on the size of the particles that constitute the bed.20 However, the fluidized-bed system of interest in the present study is basically a geometrically constrained type of system, which may be viewed as a pseudo-one-dimensional bed. Such a type of bed is essential for the observation of traveling voidage waves that propagate along the “single” dimension of the bed and do not exhibit any transverse structures along a second dimension. The good agreement between various experimentally observed and numerically simulated characteristics of such voidage waves, as will be shown in the remainder of this paper, also indicates that the pseudo-three-dimensional CFD-DEM model applied in the present study is sufficiently adequate for capturing most of the physics associated with such phenomena. Each glass bead in the CFD-DEM simulation had a diameter of 1.0 mm and a density of 2500 kg/m3. The coefficients of friction and restitution used were 0.3 and 0.9, respectively. The base of the fluidization system was allowed to undergo simple harmonic motion when desired to facilitate the study of the effects of external perturbations on the stability of the bed. To the knowledge of the present authors, this is the first computational study of vibrated liquid-fluidized beds whereby the external forcing (in the form of an oscillating base) imposed experimentally in previous work4 and the present work is also incorporated as part of the simulation. This reproduces the physical conditions present in such experiments in a realistic manner and thus facilitates corresponding studies of the underlying physics through this approach. The amplitude and frequency applied when a vibrating base was simulated were 1.5 times greater than the diameter of a glass bead (i.e., 1.5 mm) and 2 Hz, respectively, to have a basis for comparison with actual experimental observations. The superficial velocities

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1377

of the liquid used were 0.018 and 0.030 m/s. In simulations that involve an oscillating base, the moving piston was assumed not to impart a significant amount of momentum to the fluid4 and the superficial velocities used were maintained at a constant value. The method of implementation of the CFD-DEM model in the present study follows that used by Lim and co-workers21,22 in their numerical study of pneumatic conveying of granular materials through vertical, horizontal, and inclined pipes. The translational and rotational motions of individual solid particles are governed by Newton’s laws of motion:

dVi mi

dt

N

)

(fc,ij + fd,ij) + mig + ff,i + fB ∑ j)1 dωi Ii

dt

(1)

N

)

Tij ∑ j)1

(2)

where mi and Vi are the mass and velocity of the ith particle, respectively; N is the number of particles in contact with the ith particle; fc,ij and fd,ij are the contact and viscous contact damping forces, respectively; ff,i and fB are the fluid drag and buoyancy force due to the interstitial fluid, respectively; Ii is the moment of inertia of the ith particle; ωi is its angular velocity; and Tij is the torque that results from contact forces, which causes the particle to rotate. The effects of added mass force and Basset force are deemed to be of secondary importance only toward the simulation of voidage wave instabilities and, therefore, have been neglected in this calculation. Contact and damping forces must be calculated using forcedisplacement models that relate such forces to the relative positions, velocities, and angular velocities of the colliding particles. Following the research of Lim and co-workers,21,22 a linear spring-and-dashpot model was implemented for the calculation of these collision forces. With such a closure, interparticle collisions are modeled as compressions of a perfectly elastic spring while the inelasticities associated with such collisions are modeled by the damping of energy in the dashpot component of the model. Collisions between particles and a wall may be handled in a similar manner, but with the latter not incurring any change in its momentum. In other words, a wall at the point of contact with a particle may be treated as another particle but with an infinite amount of inertia. The normal (fcn,ij, fdn,ij) and tangential (fct,ij, fdt,ij) components of the contact and damping forces are calculated according to the following equations:

fcn,ij ) -(κn,iδn,ij)ni

(3)

fct,ij ) -(κt,iδt,ij)ti

(4)

fdn,ij ) -ηn,i(Vrni)ni

(5)

fdt,ij ) -ηt,i[(Vrti)ti + (ωi × Ri - ωj × Rj)]

(6)

where κn,i, δn,ij, ni, ηn,i and κt,i, δt,ij, ti, ηt,i are the spring constants, displacements between particles, unit vectors, and viscous contact damping coefficients in the normal and tangential directions, respectively; Vr is the relative velocity between particles and Ri and Rj are the radii of particles i and j, respectively. If |fct,ij| > |fcn,ij| tan φ, then “slippage” between two contacting surfaces is simulated, based on Coulomb-type friction law, i.e. |fct,ij| ) |fcn,ij| tan φ, where tan φ is analogous to the coefficient of friction.

The fluid drag force model proposed by Di Felice,23 which is applicable over a wide range of particle Reynolds numbers, was used to evaluate the drag force. Because this empirical model was developed by fitting with previously reported experimental data for both fixed and fluidized beds, it was assumed that the various interaction effects, such as those due to added mass and history forces, which must necessarily be present in the physical experiments used for correlation, have been incorporated empirically in the model. As such, these effects are not further computed explicitly through a separate set of equations in the present simulations. It may also be noted that Derksen and Sundaresan24 reported in their recent direct numerical simulation of liquid-fluidization systems that the contribution from both added mass force and history force is negligible, with the former being no more than 0.4% of the total hydrodynamic force present. The equations in the Di Felice fluid drag force model include the following:

ff,i ) ff0,i i-χ

(7)

ff0,i ) 0.5cd0,i Ff πRi2|ui - Vi|(ui - Vi)

(8)

[

]

(1.5 - log10 Rep,i)2 χ ) 3.7 - 0.65 exp 2

(

cd0,i ) 0.63 + Rep,i )

4.8 Rep,i0.5

)

2FfRi|ui - Vi| µf

2

(9) (10)

(11)

where ff0,i is the fluid drag force on particle i in the absence of other particles, χ an empirical parameter, i the local average porosity in the vicinity of particle i, cd0,i the drag coefficient, Rep,i the Reynolds number based on particle diameter, Ff the fluid density, µf the fluid viscosity, and ui the fluid velocity. The motion of the continuum liquid phase is governed by the Navier-Stokes equations, with interphase interactions taken into account as an additional source term in the momentum equation. The computational domain was divided into uniform 2 mm × 2 mm computational cells, and all quantities such as velocities and pressure were assumed to be constant over each cell. According to the local averaging technique, the continuity and momentum equations are expressed in terms of local mean variables over each computational cell:

∂ + ∇‚(u) ) 0 ∂t

(12)

∂(Ffu) + ∇‚(Ffuu) ) -∇P + ∇‚(µf∇u) + Ffg - F ∂t (13) where u is the velocity vector,  the local average porosity, P the fluid pressure, and F the source term due to fluid-particle interaction. The source term F for a particular computational cell was calculated by summing the fluid drag forces on all particles present within that cell: n

F)

ff,i ∑ i)1 ∆V

(14)

where ∆V is the volume of a computational cell and n is the number of particles present in the cell.

1378

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

Figure 2. Voidage waves in a liquid-fluidized bed operating under the following conditions: liquid superficial velocity at inlet, 0.030 m/s; vibrating amplitude of base, 1.5 mm; and frequency of base, 2 Hz. Time interval between each frame shown is 0.05 s. Dimensions of the system are 16 cm (height) × 2 cm (width). The center of each major dense region of the voidage wave is enclosed in a dashed circle to aid in visualizing the propagation of the wave along the bed. Particles are color-coded according to the vertical velocity, increasing from blue (-0.030 m/s) to green to red (0.030 m/s).

Results and Discussion CFD-DEM Simulation of Voidage Waves. Figure 2 shows four consecutive snapshots of the fluidized-bed system obtained from the CFD-DEM simulation. The time interval between each snapshot is 0.05 s, and the dimensions are 16 cm (height) by 2 cm (width). Alternating bands of dense and dilute solid concentrations may be clearly discerned from the figure. There are approximately two major dense and dilute phases present in the system, under the present set of operating conditions. The coherent motion of these phases of the voidage wave up along the fluidized bed may also be observed. As such, the present simulation has been successful in reproducing the main qualitative feature of the phenomenon associated with voidage wave instability in a liquid fluidization system. It may also be stated that, in the absence of a vibrating base, the fluidized bed was observed to expand slightly upon introduction of liquid with a superficial velocity and remain homogeneously fluidized with minimal tendency to develop such voidage wave instability (data not shown for brevity). This implied that the system was

intrinsically stable in the absence of external perturbations while any internal noises were not sufficiently significant to cause instability. This was true for both liquid superficial velocities investigated in the initial phase of the present study (0.018 and 0.030 m/s). At the lower liquid superficial velocity of 0.018 m/s and in the presence of a vibrating base, a small amount of voidage waves could be observed in the system. These formed at the vibrating base but were propagated only a short distance up the bed. As the fluidized bed at this low superficial velocity was only expanded slightly and close to a packed condition, the likely reason for attenuation of the voidage waves could be the high effective solid viscosity. In contrast, when the liquid superficial velocity applied was 0.030 m/s, such that the bed was expanded to a larger extent, voidage instabilities in the form of waves of high and low solid concentrations could be observed traveling up the expanded bed, as shown in Figure 2. These clearly show the unstable nature of the system toward external perturbations and the convective characteristic of the resulting instability.

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1379

Figure 3. Voidage waves in a liquid-fluidized bed containing inelastic solid particles with coefficient of restitution of 0.1. Other parameters and operating conditions are as described in the caption of Figure 2.

It may be noted that Duru et al.4 demonstrated in their experimental investigation of such voidage wave instabilities that the shape of these waves evolves along the bed and may reach a saturated state after a sufficient distance of wave propagation. Nicolas et al.14 applied a cnoidal function to describe quantitatively the shapes of such waves after attaining the state of saturation. In the present study, a shorter bed height of ∼12 cm was used in the CFD-DEM simulations in consideration of the computational requirements involved. A similar bed height was also used in the experiments, to have a good basis for comparison between the two approaches. As seen in Figure 2 previously, the structure of the voidage wave formed is such that approximately two wavelengths of the wave are developed over the entire bed. In other words, the voidage wave propagates for a distance of approximately two wavelengths along the bed and therefore may not reach the state of saturation. This may also be the reason, in the absence of an imposed external perturbation, where the natural or intrinsic perturbations

due to flow of liquid at the base would then be only infinitesimally small in such cases, no voidage wave was observed to develop over the entire fluidized bed that was simulated. However, it may be noted that the state of saturation of the voidage waves (or lack thereof) does not affect the interpretation of the results obtained throughout this entire study. Duru et al.4 observed empirically the bed height required for saturation of voidage waves to occur by measuring experimentally the shapes of the waves at various vertical positions along the bed. Nevertheless, this aspect of the voidage waves observed in the present fluidization system used will be discussed more quantitatively in the following section of this paper. The nature of particle-particle and particle-wall collisions may be an important factor that affects the structure of the voidage waves formed, especially in the relatively dense fluidization system used in the present study, where every particle is almost certainly always in contact with one or more other particles at any instant of time. In particular, the elasticity

1380

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

Figure 4. Voidage wave formed at a vibrating frequency of 1 Hz. Other operating parameters are as described for Figure 2. Particles are color-coded according to the vertical velocity, increasing from blue (-0.030 m/s) to green to red (0.030 m/s).

or inelasticity of such collisions may have an important role in such voidage wave phenomena. Derksen and Sundaresan24 performed direct numerical simulations of voidage waves in liquid fluidization systems similar to those considered in the present study by applying the lattice-Boltzmann method to compute the fluid flow field and a hard-sphere collision model to simulate solid particles. They reported that the coefficients of restitution and friction had only weak influences on the structure of the voidage waves formed when these parameter values were varied in the range of 0.8-1.0 and 0.0-0.2, respectively. In this study, to characterize the importance of elasticity of collisions on the voidage wave structure, the previous simulation where a liquid superficial velocity of 0.030 m/s was applied was repeated for the case where the coefficient of restitution of the solid particles was 0.1, representing inelastic particles. All other operating conditions were maintained identical to those applied in the previous case. Figure 3 shows the resulting structure of the voidage waves formed in the fluidization system. The dense and dilute phases and the coherent motions of these phases are similar to those observed

in Figure 2 previously for elastic particles with a coefficient of restitution of 0.9. As such, the elasticity of collisions among particles does not seem to have a significant effect, at least in a qualitative sense, on voidage wave formation in the present liquid fluidization system. However, there may be subtle quantitative differences such as a thicker layer of particles descending along the side walls of the column and the formation of streaks of particles that have relatively high velocities near the top of the bed, as may be discerned from the snapshots shown in Figure 3. Similarly, the frequency of vibration of the base also may have an important role in determining the structure of the voidage waves. It may be observed from Figure 4 that the macroscopic structure of the voidage wave has been changed with a reduction in the vibrating frequency to 1 Hz, while keeping all other parameter values identical to those used to produce the simulation result shown in Figure 2. Both the solid distribution in the fluidized bed and the velocity vectors of the particles shown in the figure suggest that there is only one major dense phase and one dilute phase present. Nevertheless, the

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1381

Figure 5. Calibration of light intensity transmitted against mean solid volume fraction, using a light scattering method. Signal intensities are normalized with respect to those obtained at a zero flow rate (packed-bed condition). Liquid superficial velocities required to achieve the respective mean volume fractions are indicated.

characteristic ascending and descending behaviors of particles in these dense and dilute phases, respectively, have remained unchanged. This aspect of voidage wave instabilities in liquidfluidized beds will not be within the scope of the present study. In the following sections, detailed quantitative analyses of the voidage wave phenomenon corresponding to the case where the liquid superficial velocity is 0.030 m/s, the coefficient of restitution of solid particles is 0.9, and the vibrating frequency of the base is 2 Hz would be provided. Solid Volume Fraction. Figure 5 shows the calibration of light intensity transmitted against the solid volume fraction obtained using the light scattering method. Each light intensity signal was averaged over 200 s and measurements were made at distances of 1, 5, and 10 cm above the base of the column. It may be observed that the transmitted-light intensities varied strongly with both solid volume fraction and flow rate. This calibration curve was used for all subsequent solid volume fraction measurements. Spatial averaging of solid fraction was performed at different instants of time by calculating the ratio of total volume of

particles within a 2 mm (height) × 2 cm (width) domain 5 cm above the vibrating base to the volume of that domain. An ensemble average was then calculated by superimposing together several such profiles, ensuring that the phases were matched to each other in this process. Figures 6a and 6b show the ensemble averaged solid fraction profile at a distance of 5 cm above the vibrating base, with respect to time, and the corresponding power spectrum obtained by fast Fourier transform (FFT) of the instantaneous solid fraction signal, respectively, obtained from CFD-DEM simulations. Because of high computational requirements, the ensemble averaging was performed using 20 sets of solid fraction data from the CFD-DEM simulations, each over a time interval of 3 s. Figure 6c shows the corresponding data obtained from the experiments conducted and ensemble averaged over 120 waveforms for the solid volume fraction. Similarly, Figure 6d shows quantitatively that the characteristic frequency of the voidage wave is also equal to that of the vibrating base. In addition, there may be higher or lower harmonics in the vicinity of this characteristic frequency (in the range of 1-5 Hz) (see Figure 6b). Figures 7 a and 7b show the corresponding power spectral densities of the solid fraction obtained at positions 1 and 10 cm above the vibrating base, respectively, from the CFD-DEM simulations. Upon comparison, it may be observed that the driving frequency of 2 Hz of the base has been manifested as the characteristic frequency of the voidage wave up to a height of 10 cm along the fluidized bed. However, the various harmonics mentioned are propagated a smaller distance. Other than the frequency, which is exactly twice (4 Hz) that of the characteristic frequency, these are no longer present in the solid fraction signal at a distance of 10 cm from the vibrating base. The corresponding power spectral densities obtained from experiments show a similar attenuation of the voidage wave along the height of the bed. A qualitatively comparable disappearance of harmonics may also be observed (see insets to Figures 7a and 7b). Duru et al.4 showed, in their stability diagram, that there exists a maximum limit of the frequency of an external perturbation for voidage wave instabilities to develop in such liquid-fluidized-bed systems. The higher harmonics of voidage waves previously discussed may correspond to higherfrequency modes of such instability waves, which have been shown by the previous workers to be more stable. These have

Figure 6. Ensemble-averaged variation of spatially averaged solid fraction at a distance of 5 cm above the vibrating base, with respect to time, obtained from (a) CFD-DEM simulations and (c) experiments. Also shown is the corresponding power spectral density of the time-varying solid fraction obtained from (b) CFD-DEM simulations and (d) experiments.

1382

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

Figure 8. Voidage structure obtained from (a) linear stability analysis of the continuum model8 and (b) CFD-DEM simulations. Solid fraction profile over a voidage wave obtained from (c) Glasser et al.,8 where ky (dimensionless) ) 0.204, and (d) CFD-DEM simulations using 3250 particles, where the corresponding dimensionless ky value is estimated to be ∼0.060.

Figure 7. Power spectral density of the solid fraction profile obtained at distances of (a) 1 cm and (b) 10 cm above the vibrating base, from CFDDEM simulations; insets show the corresponding power spectral densities obtained from experiments. (c) Evolution of voidage wave shape with vertical position along the bed.

a tendency to be attenuated with vertical position along the bed. As such, the wave that persists at higher vertical positions, as observed in the present study, may be the dominantly unstable mode whose frequency matches that of the external perturbation imposed, while other more-stable modes are eventually damped out. Figure 7c shows the evolution of the shape of voidage waves obtained from the CFD-DEM simulations. It may be observed that the wave shape at a distance of 5 cm above the vibrating base is sinusoidal with a sharp peak, whereas near the top of the expanded bed at 14 cm above the base, it becomes flatter, showing signs of forming a plateau. This indicates evolution of the voidage waves toward the state of saturation as they propagate along the bed and is consistent with the observations reported by Duru et al.4 It may be noted that a complete state of saturation of the voidage waves may not be attained in the present fluidization system, as noted earlier. Nevertheless, the results shown here verify the adequacy of the

current CFD-DEM model in capturing the physics involved and associated with voidage wave evolution and saturation. Glasser et al.8 conducted linear stability analyses of a uniform base state of a solid-liquid suspension of 1.0 mm glass beads in water, using the continuum continuity and momentum equations to determine the most unstable one-dimensional eigenvalue and the corresponding eigenfunction for the solid fraction. The grayscale plot of this eigenfunction, corresponding to the volume fraction of solids obtained by these authors, is compared with a typical snapshot of the physical structure of a voidage wave obtained from the present CFD-DEM simulations (Figures 8a and 8b). It may be observed, in Figure 8b, that two dense phases (regions of high solid concentration) are present in the section of the fluidized bed shown with a relatively more dilute region in between them. On closer inspection, it may be easily discerned that this dilute region is further comprised of alternating layers of high and low solid concentrations. It is highly unlikely that such structures are voidage waves with higher wavenumbers corresponding to higher-frequency perturbation modes, because these are theoretically more stable (according to Duru et al.4). Furthermore, no characteristic peak at the characteristic frequency of such structures was observed in the power spectral density profiles shown earlier, even though they were also present in those simulation results. As such, these layered structures may represent intermediate stages of decompactification of an originally dense phase whereby particles start to move downward, causing the dense region to disintegrate into smaller layers of particles. These structures will not be considered further in the present study. In Figure 8b, the characteristic wavelength of the voidage waves (λy), which is given by the distance between the two dense phases indicated, is estimated to be ∼6 cm. Following Glasser et al.,8 the

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1383

Figure 9. Instantaneous particle velocity vector field in a 1.5 cm (height) × 1.0 cm (width) section at a distance of 5 cm above the vibrating base, obtained from (a) CFD-DEM simulations and (b) experiments. Snapshots are shown at intervals of 0.2 s.

characteristic dimensionless wavelength (λ′y) was obtained by nondimensionalizing, with respect to a characteristic length scale (given by L ) xAVt/(Fpg)) and the corresponding dimensionless wavenumber was calculated by ky ) 2π/λ′y. Here, A is a constant equal to 50 cP and Vt is the particle terminal velocity. The resulting value for ky was estimated to be ∼0.060, which is less than the value of 0.204 that was applied by Glasser et al.8 in their linear stability analysis. The corresponding solid fraction profiles are compared quantitatively in Figures 8c and 8d. There is a qualitative agreement in the voidage wave structure obtained from the two different calculation approaches. The present solid fraction profile shows a lower overall solid fraction along one wavelength of the voidage wave, in comparison with that obtained by Glasser et al.8 This may be due to the lower overall solid fraction used in the present simulation (∼0.5), compared with that used by the previous researchers (∼0.57). It may also be noted that the voidage waves computed by the previous researchers, following the principle of linear stability analysis, resulted from perturbations applied over a homogeneously fluidized system, whereas, in the present approach, an oscillating base was simulated to introduce external perturbations into the system to produce the type of voidage waves presented thus far. This difference in the type of boundary conditions imposed may have a strong effect on the eigenfunctions calculated via linear stability analysis, although the method to impose an oscillating solid boundary in such analyses remains elusive to date. Nevertheless, the qualitative agreement in the

structure of the voidage waves and their corresponding characteristic wavenumbers calculated using the two different approaches provides a strong indication that the combined CFD-DEM approach is also suitable for the numerical study of such voidage wave instabilities in liquid-fluidized-bed systems. Solid Velocity. Figure 9a shows the instantaneous velocity vectors of particles in a section of the system ∼5 cm above the vibrating base obtained from the CFD-DEM simulation. Correspondingly, Figure 9b shows the velocity vectors obtained from the experimental setup using the PIV system. Each frame of the figure shows approximately half a channel and has dimensions of 1.5 cm (height) × 1.0 cm (width) with one wall of the channel on the left and the centerline on the right. The frames were captured at intervals of 0.2 s. Each velocity vector was obtained by locally averaging over an interrogation window containing ∼3-4 solid particles. In principle, this method can be applied anywhere within the column, including close to the wall, as long as images of the solid particles can be captured by the PIV system. However, because of the presence of the wall, fewer solid particles are available within an interrogation window in this region for local averaging to calculate velocity vectors. Furthermore, the power of the laser had to be set sufficiently high so that it could penetrate the dense solid packing to illuminate those particles close to the center of the column. This resulted in over-reflection by some particles close to the wall which were nearer to the laser source. As such, the

1384

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

Figure 10. Ensemble-averaged variation of the spatially averaged vertical component of solid velocities at a distance of 5 cm above the vibrating base, with respect to time, obtained from (a) CFD-DEM simulations and (c) experiments. Also shown is the corresponding power spectral density of the time-varying solid velocities obtained from (b) CFD-DEM simulations and (d) experiments.

region close to the column wall was not included in the calculation of velocity vectors, to avoid possible inaccuracies. Nevertheless, the figure illustrates the unique behavior of solid particles in a liquid-fluidization system in the presence of voidage waves. Here, adopting an Eulerian point of view, it may be observed that particles switch periodically between generally upward and downward motions. These respectively correspond to the passage of dense and dilute phases of the voidage wave through the particles. In other words, when a dense phase of the wave propagates through a section of the bed, particles in that section were observed to be moving in the upward direction and vice versa. For the present case studied, the frequency and amplitude of the vibrating base were 2 Hz and 1.5 mm, respectively. Both Figures 9a and 9b show that the characteristic frequency of the periodic motion of the solid particles is also ∼2 Hz. Figure 10a shows the ensemble-averaged profile of the vertical component of solid velocities at a distance of 5 cm above the vibrating base, with respect to time. As with the case for solid fraction data described previously, the ensemble averaging was performed using 20 sets of solid velocity data from the CFD-DEM simulations, each over a time interval of 3 s. This profile quantitatively shows the oscillatory behavior of particle motion over six cycles. Figure 10b shows the corresponding power spectrum obtained from the original instantaneous velocity signals. Instantaneous velocity data were collected at intervals of 0.01 s from the CFD-DEM simulation for the calculation of the power spectrum. It may be observed that the characteristic frequency of the oscillatory solid velocity matches that of the vibrating base, as mentioned earlier. Figure 10c shows the waveform for solid velocity obtained by ensemble averaging over 120 sets of PIV experimental data. Each set of instantaneous solid velocity data was taken over a duration of 3 s, and the periodicity of each waveform was observed to be fairly uniform and reproducible (data not shown for brevity). These individual sets of instantaneous solid velocity waveforms were then superimposed on each other, making sure that the phases were matched (peak to peak and trough to trough) before ensemble averaging was performed. The number of experimental data sets used in the ensemble averaging has been verified to be sufficiently large to ensure statistical invariance of the resulting waveform. To the knowledge of the authors, this is the first report of such a statistical representation of the solid

Figure 11. Granular temperature profiles of solids at various positions above the vibrating base obtained from (a) CFD-DEM simulations and (b) experiments. The inset to panel (a) shows the granular temperature plot over the entire fluidized bed, from the CFD-DEM simulation. The range of the color scale used is 2.1 × 10-5 m2/s2 (blue) to 3.2 × 10-4 m2/s2 (red). Granular temperatures are higher at the bottom of the bed and decrease with vertical position along the bed. The origin of the horizontal position is the left lateral wall of the bed.

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1385

Figure 12. Positions of four representative particles at intervals of 1 s; the vibrating frequency of the base applied is (a) 2 Hz and (b) 1 Hz. Panels (c) and (d) show the corresponding data obtained from experiments.

velocity in a liquid-fluidized bed exhibiting voidage wave instability. It may be noted that the instantaneous solid vertical velocity waveforms show both positive and negative values, in correspondence with the direction of the PIV velocity vectors observed in Figure 9b. However, such instantaneous velocity waveforms are not symmetrical about the zero velocity value. Instead, the magnitudes of the positive velocities are usually larger than those of their negative counterparts. As such, the solid velocity waveform after ensemble averaging only exhibits negative solid velocity values at the trough positions and positive values elsewhere. Nevertheless, the power spectral densities for both simulation and experimental data (Figures 10b and 10d) show that the dominant oscillating frequency of solid velocity is equal to that of the driving frequency of the external perturbation imposed. Granular Temperature. Figures 11a and 11b show the granular temperature profiles of the solids in the liquid-fluidized bed at distances of 1, 5, and 10 cm above the base, obtained from the CFD-DEM simulations and experiments, respectively. The local granular temperature was calculated based on the average of the sum of the squared fluctuating velocity components over a domain with a size of 2 mm (height) × 1 mm (width):

T)

1

N

∑ [u(i) - uj]2 + [V(i) - Vj]2

2 i)1

(15)

where N is the number of particles in the local domain; u(i) and V(i) are the ith horizontal and vertical velocity components, respectively; and uj and Vj are mean horizontal and vertical velocity components, respectively. It may be observed that granular temperatures are higher at the bottom of the bed and decrease as the bed heights increase. The granular temperatures

calculated from both simulation and experimental data are of similar orders of magnitude (∼10-5 m2/s2). There is closer agreement at larger bed heights and some deviations between the two sets of profiles near the bottom of the bed. Specifically, the granular temperature at 10 cm above the base is approximately in the range of (3.0-5.0) × 10-5 m2/s2 for both simulation and experimental data. At a distance of 5 cm above the base, the corresponding values are ∼10.0 × 10-5 and ∼5.0 × 10-5 m2/s2, respectively, whereas at a distance of 1 cm above the base, the average granular temperatures obtained from simulation and experiments are ∼12.0 × 10-5 and ∼6.0 × 10-5 m2/s2, respectively. The granular temperature profile over the entire bed can also be observed from a plot constructed from the simulation data (see inset of Figure 11a). The plot shows that the bed has higher energy in the near-base region and the energy is dissipated as one progresses upward along the bed. This characteristic of the system, whereby granular temperature is a function of bed height over the entire fluidized bed, also supports the point made earlier that voidage waves formed in the fluidized-bed system used in the present study do not reach a state of saturation. On the other hand, granular temperature values would be expected to reach saturated values beyond a certain bed height, where the voidage waves do reach saturation. The granular temperature profile shown in the inset of Figure 11a also indicates the presence of high and low Reynolds stresses near the bottom and top of the bed, respectively, and seems to suggest possible correlations between such stresses and the onset of instability. Instabilities in the form of voidage waves originate at the vibrating base where the granular temperature is highest and are gradually attenuated as they travel upward along the bed toward regions of lower granular temperatures. This observation is consistent with the conclusion made by Ham et al.,12 whose flow visualization experiments

1386

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

Figure 13. Positions of particles exhibiting localized motion over a period of 60 s at intervals of 1 s, observed in CFD-DEM simulations.

and detailed quantitative measurements showed Reynolds stresses due to particle velocity fluctuations to be the main contributor to the formation of voidage instabilities. Localization of Solids. Despite the convective nature of the voidage waves, which originate at the vibrating base and travel upward along the length of the fluidized bed, one interesting feature of such instabilities observed in the present CFD-DEM simulations is the highly localized motion of individual particles. As discussed earlier, solid particles move upward when a dense phase of the voidage wave passes through and settle downward in the dilute phase of the wave. However, the overall motion of some particles was observed to be highly restricted to a small region within the system. Taking a Lagrangian point of view, the positions of four representative particles were tracked for a period of 6 s, corresponding to 12 cycles of the vibrating base. Figure 12a shows the position of each particle at intervals of 1 s. Regardless of the precise starting location of the particle within the system, it may be seen that there is very restricted motion in both the axial (vertical) and lateral (horizontal) directions. Each particle seems to be “trapped” within a small region in space. When compared to the structure of the voidage wave presented in Figure 2, the dimensions of such regions are

also comparable to the length scale of each dense phase of the wave. However, when the base is vibrated at a frequency of 1 Hz, with all other operating parameters unchanged, the size of the region in which each particle is trapped remains substantially similar (see Figure 12b). The extent of motion of individual particles does not seem to be significantly affected by the frequency of the vibrating base. Similar observations were made from the experiments conducted where a few arbitrarily selected particles were dyed and tracked visually for a period of time (see Figures 12c and 12d). For characterization of such localized particle motion on a relatively short time scale, the time interval of 6 s as used here was deemed to be sufficiently long, compared to the characteristic periodic times of 0.5 and 1.0 s of the system when the base was vibrated at 2 and 1 Hz, respectively. This phenomenon of localized solid motion in a liquid-fluidized bed exhibiting voidage wave instabilities does not seem to have been reported previously. Such behaviors of solids in a liquidfluidized bed system may have important implications in actual applications, such as with regard to the type of mixing behavior that can be achieved during the operations. This study of localized solid motion was extended to include its characterization over longer time scales. It was observed in

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1387

Similarly, those captured clearly using a high-speed video camera and observed to exhibit localized motion in Figure 14 are in close proximity to the near wall of the tube. The six dyed particles, labeled “1”-“6” in each of the panels in the figure, are observed to remain close to their original positions within the bed throughout 60 s of observation. Such observations seem to imply that two types of recirculation patterns of solid particles may be present with the liquid-fluidized bed system, exhibiting voidage wave instabilities. A global recirculation pattern results in solid particles traversing throughout the entire height of the bed while a local recirculation pattern results in the localization behavior of solid particles reported in the present study. As mentioned previously, such behaviors may lead to uneven mixing, poor heat- or mass-transfer rates, and even unexpected dead zones within a liquid-fluidized bed and therefore should be of practical relevance to industrial applications. Dispersion of Solids. Following the previous qualitative discussion on localization of solid motion, a more quantitative approach was used to characterize the nature of such motions. Solids drifted away from their initial positions gradually, much like the diffusion of material from a region of high concentration to another of lower concentration. To characterize this phenomenon quantitatively, the concept of dispersion was used to examine this drifting behavior of particles within the liquidfluidized bed. The following dispersion model developed in analogy to diffusion in unbounded systems was applied:

|r(t) - r(0)|2 ) 2Dt

Figure 14. Snapshots of a section of the fluidized bed at a distance of 5.0 cm above the vibrating base, containing some dyed particles captured using a high-speed video camera at intervals of 10 s. In both cases, the vibrating frequency of the base applied is 2 Hz.

both the CFD-DEM simulations and experiments that some particles within the fluidized bed continued to exhibit such localized motion over a period of 60 s, corresponding to 120 cycles of the vibrating base when its frequency was set at 2 Hz (see Figures 13 and 14, respectively). Figure 13 shows examples of the periodic type of movements that such particles undergo, which results in their localization. It may be observed that particles have a tendency to ascend near the center of the column and descend along the wall, which results in a local recirculation pattern. The spatial extension over which such restricted motion occurs is now determined to be larger than those observed previously over shorter time scales. In addition, most of the particles illustrated here seem to lie fairly close to the tube walls.

(16)

where r is the position vector of the tagged particle, D the dispersion coefficient, and t the time over which the tagged particle has been tracked. Here, data from the CFD-DEM simulations were used to calculate the dispersion coefficient. Figure 15a shows the mean squared vertical displacement of an arbitrarily selected particle over an interval of 3 s. The excellent fit between the simulation data and the dispersion model used confirms that particles in the bed indeed exhibited a diffusion type of motion. To the knowledge of the present authors, this has not been reported in previous studies of voidage wave instabilities in liquid-fluidized beds. The corresponding dispersion coefficient can be calculated from the gradient of the straight line fitted to the simulation data. Figure 15b shows that the same type of diffusion behavior of individual particles could be observed at various other positions of the fluidized bed exhibiting voidage wave instabilities. Furthermore, dispersion coefficients were determined to be larger at higher bed levels, in contrast to the trend observed previously for granular temperatures (see Figure 15c). However, as with the study of granular temperature profiles discussed earlier, this feature also suggests that the voidage waves formed in the present fluidizedbed system have not reached the state of saturation after propagating over the entire bed height. The values of dispersion coefficients may also be expected to become saturated beyond a sufficient bed height after the voidage waves become saturated. Thus, a consistent conclusion regarding the state of saturation of the voidage waves developed within the present system may be inferred from qualitative inspections of the wave structure obtained from CFD-DEM simulations, as well as quantitative analyses of the granular temperature and dispersion coefficient profiles. The inter-relationships between granular temperature, dispersion coefficient and their influences on the formation, propagation, and attenuation of voidage wave instabilities in a liquid-fluidized bed may be the basis of an interesting study beyond the present work.

1388

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007

of such oscillatory motions of solid particles was obtained by a fast Fourier transform of the vertical component of solid velocities and observed to match the vibrating frequency of the base. Similarly, the characteristic frequencies of the solid fraction, which also exhibited periodic variations, with respect to time, were determined to be equal to the vibrating frequency of the base. The voidage waves formed as a result of instability in such liquid-fluidized-bed systems are traveling waves with dense and dilute phases being convected along the bed. The effective wavelength of such waves, in terms of the distance between two adjacent major dense or dilute phases, was determined to be dependent substantially on the frequency of vibration of the base. However, the motion of individual particles was observed to be highly restricted to a small region in space whose dimensions did not seem to be dependent significantly on the vibrating frequency. Despite the convective nature of such waves, the phenomenon is a highly localized one at the individual particle level. This may have important implications for the operations of such liquid-fluidization systems, with regard to the effectiveness with which solids mixing can be accomplished in such systems and others. Quantitatively, such motions of individual particles were determined to be diffusive in nature. This type of behavior was adequately described by a simple dispersion model used in the present study. Acknowledgment This study has been supported by the SERC (A*STAR) and the Singapore-MIT Alliance (under Grant Nos. R-279-000208-305 and MEBCS-003-Project003, respectively). The authors thank Prof. Sankaran Sundaresan for many helpful discussions on this project. Nomenclature

Figure 15. Variation of mean squared vertical displacement with time of (a) an arbitrarily selected particle and (b) particles with different initial positions (1.8, 3.0, 4.0, 5.0, 6.5, 7.5, 8.5, and 9.5 cm above the base) within the bed. Panel (c) shows the variation of the particle dispersion coefficient, relative to different initial positions above the base.

Conclusions The characteristics of one-dimensional voidage waves in a liquid-fluidized bed have been investigated both experimentally and computationally in the present study. The computer simulations that combine the discrete element method with computational fluid dynamics calculations have successfully reproduced the phenomenon of voidage wave formation in a liquid-fluidized bed system subjected to external perturbations. A uniform inlet velocity was applied and the base of the system was allowed to undergo simple harmonic motions, to simulate the introduction of small external disturbances. Voidage waves that consist of alternating regions of high and low solid concentrations were observed to form and travel in a coherent manner along the fluidized bed. Solid particles were observed to move upward when a dense phase of the wave passed through their positions and settle downward otherwise. The characteristic frequency

cd0,i ) drag coefficient D ) dispersion coefficient fc,ij ) contact force fcn,ij ) normal component of contact force fct,ij ) tangential component of contact force fd,ij ) viscous contact damping force fdn,ij ) normal component of viscous contact damping force fdt,ij ) tangential component of viscous contact damping force ff,i ) fluid drag force ff0,i ) fluid drag force in the absence of other particles F ) source term due to fluid-particle interaction g ) gravitational acceleration Ii ) moment of inertia ky ) dimensionless wavenumber L ) characteristic length mi ) particle mass n ) number of particles in a computational cell ni ) unit vector in normal direction N ) number of particles P ) fluid pressure r ) position vector of particle Ri ) radius of particle Rep,i ) Reynolds number based on particle diameter t ) time ti ) unit vector in the tangential direction T ) granular temperature Tij ) torque ui ) fluid velocity Vi ) particle velocity Vt ) terminal velocity

Ind. Eng. Chem. Res., Vol. 46, No. 4, 2007 1389

u ) mean horizontal velocity component V ) mean vertical velocity component ∆V ) volume of a computational cell Greek Symbols δn,ij ) displacements between particles in the normal direction δt,ij ) displacements between particles in the tangential direction i ) local average porosity κn,i ) spring constant for normal collisions κt,i ) spring constant for tangential collisions λy ) wavelength λ′y ) dimensionless wavelength µf ) fluid viscosity ηn,i ) viscous contact damping coefficients in the normal direction ηt,i ) viscous contact damping coefficients in the tangential direction Ff ) fluid density Fp ) particle density χ ) empirical parameter ωi ) angular velocity Literature Cited (1) Davidson, J. F., Clift, R., Harrison, D., Eds. Fluidization; Academic Press: London, 1985. (2) Clift, R.; Grace, J. R. Continuous bubbling and slugging. In Fluidization; Davidson, J. F., Clift, R., Harrison, D., Eds.; Academic Press: London, 1985. (3) Didwania, A. K.; Homsy, G. M. Flow regime and flow transitions in liquid fluidized beds. Int. J. Multiphase Flow 1981, 7, 563. (4) Duru, P.; Nicolas, M.; Hinch, J.; Guazzelli, EÄ . Constitutive laws in liquid-fluidized beds. J. Fluid Mech. 2002, 452, 371. (5) El-Kaissy, M. M.; Homsy, G. M. Instability waves and the origin of bubbles in fluidized beds. Part I: Experiments. Int. J. Multiphase Flow 1976, 2, 379. (6) Batchelor, G. K. A new theory of the instability of a uniform fluidized bed. J. Fluid Mech. 1988, 193, 75. (7) Batchelor, G. K.; Nitsche, J. M. Instability of stationary unbounded stratified fluid. J. Fluid Mech. 1991, 227, 357. (8) Glasser, B. J.; Kevrekidis, I. G.; Sundaresan, S. Fully developed travelling wave solutions and bubble formation in fluidized beds. J. Fluid Mech. 1997, 334, 157.

(9) Go¨z, M. F.; Sundaresan, S. The growth, saturation and scaling behavior of one- and two-dimensional disturbances in fluidized beds. J. Fluid Mech. 1998, 362, 83. (10) Homsy, G. M.; El-Kaissy, M. M.; Didwania, A. Instability waves and the origin of bubbles in fluidized bedssII. Comparison with theory. Int. J. Multiphase Flow 1980, 6, 305. (11) Nicolas, M.; Chomaz, J.-M.; Guazzelli, EÄ . Absolute and convective instabilities of fluidized beds. Phys. Fluids 1994, 6, 3936. (12) Ham, J. M.; Thomas, S.; Guazzelli, EÄ .; Homsy, G. M.; Anselmet, M.-C. An experimental study of the instability of liquid-fluidized beds. Int. J. Multiphase Flow 1990, 16, 171. (13) Nicolas, M.; Chomaz, J.-M.; Vallet, D.; Guazzelli, EÄ . Experimental investigations on the nature of the first wavy instability in liquid-fluidized beds. Phys. Fluids 1996, 8, 1987. (14) Nicolas, M.; Hinch, J.; Guazzelli, EÄ . Wavy instability in liquidfluidized beds. Ind. Eng. Chem. Res. 1999, 38, 799. (15) Duru, P.; Guazzelli, EÄ . Experimental investigation on the secondary instability of liquid-fluidized beds and the formation of bubbles. J. Fluid Mech. 2002, 470, 359. (16) Anderson, K.; Sundaresan, S.; Jackson, R. Instabilities and the formation of bubbles in fluidized beds. J. Fluid Mech. 1995, 303, 327. (17) Batchelor, G. K. Secondary instability of a gas-fluidized bed. J. Fluid Mech. 1993, 257, 359. (18) Feng, Y. Q.; Xu, B. H.; Zhang, S. J.; Yu, A. B.; Zulli, P. Discrete particle simulation of gas fluidization of particle mixtures. AIChE J. 2004, 50, 1713. (19) Pandit, J. K.; Wang, X. S.; Rhodes, M. J. On Geldart Group A behaviour in fluidized beds with and without cohesive interparticle forces: A DEM study. Powder Technol. 2006, 164, 130. (20) Yang, R. Y.; Zou, R. P.; Yu, A. B. Computer simulation of the packing of fine particles. Phys. ReV. E 2000, 62, 3900. (21) Lim, E. W. C.; Wang, C. H.; Yu, A. B. Discrete element simulation for pneumatic conveying of granular material. AIChE J. 2006, 52, 496. (22) Lim, E. W. C.; Zhang, Y.; Wang, C. H. Effects of an electrostatic field in pneumatic conveying of granular materials through inclined and vertical pipes. Chem. Eng. Sci. 2006, 61, 7889. (23) Di Felice, R. The voidage function for fluid-particle interaction systems. Int. J. Multiphase Flow 1994, 20, 153. (24) Derksen, J. J.; Sundaresan, S. DNS of dense suspensions: wave instabilities in liquid-fluidized beds. Submitted to J. Fluid Mech., 2006.

ReceiVed for reView July 5, 2006 ReVised manuscript receiVed December 3, 2006 Accepted December 12, 2006 IE060864E