- Email: [email protected]
Square-well model for cohesion in fluidized beds
Chemical Engineering Science 61 (2006) 4511 – 4527 www.elsevier.com/locate/ces
Square-well model for cohesion in fluidized beds Michael W. Weber, Christine M. Hrenya ∗ Department of Chemical and Biological Engineering, University of Colorado, Boulder, CO 80309, USA Received 23 September 2005; received in revised form 10 December 2005; accepted 5 February 2006 Available online 8 May 2006
Abstract Cohesive forces are implemented into a discrete-particle, fluidized-bed simulation using a square-well potential. The square-well description treats cohesive interactions as instantaneous, binary events, thereby making it a viable option for the incorporation of cohesion into a kinetictheory-based continuum model. Cohesive forces are also incorporated into the simulation using the more elaborate Hamaker description of van der Waals forces in order to provide a basis for assessing the square-well model. Both cohesion models are implemented in the discrete-particle framework of the MFIX software package. A mapping method is also developed to convert material-specific Hamaker constants into equivalent square-well parameters. The corresponding results from the two models are compared both qualitatively and quantitatively. The predictions of the square-well model are on par with the Hamaker model with respect to mixing level, particle mobility and minimum fluidization velocity. Subtle differences are observed between the two models in cases that involved such high levels of cohesion that the particle bed could not fully fluidize. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Fluidization; Mathematical modeling; Multiphase flow; Particle; Cohesion
1. Introduction An important consideration in fluidized-bed operation is the potential impact of cohesive forces between particles. The presence of cohesion can lead to inefficient operation due to reduced gas–solid contacting and can also result in undesired behaviors such as plugging and channeling. Cohesive forces may arise from a variety of sources including liquid bridges, van der Waals forces, electrostatic forces and magnetic forces. Mathematical models (Mikami et al., 1998; Rhodes et al., 2001a,b; Tatemoto et al., 2005) for cohesive systems provide an important complement to experimental investigations since experimentally controlling the level of cohesion, independent of changing other system conditions, is non-trivial. For example, particle cohesion has been controlled experimentally by adding a liquid to the system such as water, oil or surfactant to control the volume of the liquid bridges formed between particles (Forsyth et al., 2002; McLaughlin and Rhodes, 2001; Nase et al., 2001). However, such changes may have secondary ∗ Corresponding author. Tel.: +1 303 492 7689; fax: +1 303 492 4341.
E-mail address: [email protected] (C.M. Hrenya). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.02.008
effects on other properties. For example, humidification of the air changes the density and viscosity of the fluidizing gas, as well as the mass of the particles due to adsorbed water (Tsinontides and Jackson, 1993). Computational models, on the other hand, allow for independent control of cohesive forces. In models for fluidized-bed systems, the fluid is treated as a continuum and thus a single mass and momentum balance is used for this phase. With regard to the solid phase, both discrete and continuum treatments have been utilized. The discrete approach (Lagrangian) involves the solution of a force balance for every particle present in the system, whereas the continuum approach (Eulerian) involves the solution of a single mass and momentum balance. For systems of industrial size, the Eulerian approach is more computationally efficient than the Lagrangian approach, which is currently limited to approximately 100,000 particles for desktop machines. However, the incorporation of complex particle physics (e.g., cohesion) is a more difficult task with Eulerian models. In particular, the averaging technique used to derive the continuum balances gives rise to terms, like the solid-phase stress tensor, which require a constitutive relation. The incorporation of cohesion into such constitutive quantities is considerably less straightforward than its incorporation into a discrete-particle model.
4512
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
To date, Eulerian models that include the effects of cohesion are extremely limited. A kinetic-theory-based description requires that particle interactions are binary and instantaneous in nature, whereas cohesive forces are known to vary as a function of separation distance and often result in enduring contacts. Kim and Arastoopour (2002) used the kinetic theory framework to account for cohesion via a contact bonding energy. Their model accounts for the formation of spherical particle agglomerates but does not allow for deagglomeration. A less rigorous approach was taken by McKeen and Pugsley (2003) who developed an Eulerian model for the fluidization of cohesive, fluid-catalytic-cracking particles in which cohesion is accounted for by increasing the effective particle diameter used in the drag-force relation. Efforts to include cohesion in Lagrangian models are more numerous since, as mentioned earlier, the incorporation of such physics is more straightforward. Mikami et al. (1998) incorporated a liquid-bridge model into a fluidized-bed simulation. This effort produced pressure data illustrating an increase in the minimum fluidization velocity for wet particles relative to a dry system. Kuwagi et al. (2000) applied a sintering model and studied its effects on the formation of dead zones in the fluidized bed. Rhodes et al. (2001a,b) used the same fluidized-bed simulation developed by Mikami et al., but incorporated a simpler cohesion model. Instead of the liquid-bridge model, they devised a cohesion model in which cohesive forces are only implemented when the particles are in physical contact. The magnitude of the cohesive force is set to a constant and equal to an integer multiple of the particle weight. More recently, Wang and Rhodes (2004) developed new parameters to quantify the flowability of the fluidized bed. Using parameters such as the particle mobility and the average particle speed, they defined quantitative criteria for the transition from Geldart Type A to Type C particles. Xu et al. (2002) implemented a model for van der Waals forces into a fluidized-bed simulation. The results of this effort show the formation of static force chains within the fluidized bed. Ye et al. (2004) incorporated van der Waals forces into a soft-sphere discrete-particle fluidized-bed model using the Hamaker theory. They observed the formation of channels reminiscent of Geldart Type C behavior for high levels of cohesion. Tatemoto et al. (2005) studied the effects of cohesion in a discrete-element simulation of vibrated fluidized beds. They examined the effect of vibration frequency and strength on the pressure losses in the fluidized bed. To summarize, numerous Lagrangian models for cohesive particles in fluidized beds have been shown to compare well with experimental observations. Eulerian models for cohesion are relatively scarce, and their development has been hampered by the challenge of incorporating cohesion into a framework that requires binary and instantaneous interactions between particles, whereas cohesive interactions are often enduring in nature. For example, the Hamaker description of van der Waals is characterized by enduring contacts. With these ideas in mind, a primary goal of this effort is to investigate an alternative description of cohesive forces, namely the squarewell potential. The square-well description treats cohesion as arising from impulsive events between two particles, and thus
represents a possible means of incorporating cohesion into continuum (kinetic-theory) models. This model has previously been applied to granular flows (no interstitial fluid) where it was shown to produce qualitatively correct predictions (Weber et al., 2004). The performance of the square-well approximation for cohesive forces in fluidized beds is investigated here via discrete-particle simulations. For purposes of comparison, cohesive forces are also implemented in discrete-particle simulations using the more complex Hamaker description of van der Waals forces (recall that such enduring contacts cannot be incorporated directly into kinetic theory models). A mapping method is developed to determine the value of square-well inputs given the physical parameters required by the Hamaker description. The resulting discrete-particle simulations indicate the effectiveness of the simpler, square-well model relative to the more elaborate Hamaker model. 2. Computational method The simulation domain consists of a two-dimensional rectangular system with the origin of the x- and y- axis in the bottom left corner. Spheres are restricted to motion in a twodimensional plane (i.e., thickness of fluidized bed is a single particle diameter). The bed typically used is 0.08 m wide by 0.2 m tall. For this length and width, 27 computational cells in the x-direction and 30 computational cells in the y-direction are used. In the following text, any changes from this bed size are noted. The governing equations for the complete fluidizedbed simulation are composed of fluid-phase and particle-phase balances. The fluid-phase continuum equations are handled in a standard manner, the details of which are given in Appendix A. The boundary conditions for the gas phase consist of no-slip, impermeable walls on the vertical sides of the bed. For the outflow boundary condition at the top of the bed, zero gauge pressure is specified across the entire width. Thirteen separate inlet jets are used to specify the inlet boundary conditions. The boundaries of each inlet jet are required to coincide with the boundary of a computational cell, limiting the range of sizes that can be used for the inlet jets. In most cases, each jet consists of one computational cell (2.96 × 10−4 m in width). At each inlet jet, the void fraction is set to unity and the gas velocity is set to achieve the desired superficial velocity in the bed. This inlet jet configuration provides the most uniform fluidization with the least amount of dead space (unfluidized particles) in the bed. Between each jet, the distributor plate is represented by impermeable, no-slip walls. The particle-phase equations are described below. For further detail on the computational method, the reader is referred to Weber (2004). Table 1 includes a listing of all parameters used in the simulation, unless otherwise noted. 2.1. Particle advancement The discrete-particle simulation operates using a timestepped algorithm (Allen and Tildesley, 1987). The particles are modeled as soft spheres confined to movement in the x- and
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Table 1 Input parameters used in fluidized bed simulation Parameter
Value
Number of particles (n) Fluidized-bed width (xmax ) (m) Fluidized-bed height (ymax ) (m) Computational cells in x-direction (mmax ) Computational cells in y-direction (nmax ) Solid-phase time step (ts ) (s) Gas-phase time step (tf ) (s) Particle diameter (DP ) (m) Particle density (s ) (kg/m3 ) Normal spring constant (Kn ) (kg/s2 ) Tangential spring constant (Kt ) (kg/s2 ) Normal damping coefficient (n ) (kg/s) Tangential damping coefficient (n ) (kg/s) Friction coefficient () Coefficient of restitution (e) Collision duration (tc ) (s) Normal wall spring constant (Kn-wall ) (kg/s2 ) Tangential wall spring constant (Kt -wall ) (kg/s2 ) Normal wall damping coefficient (n-wall ) (kg/s) Tangential wall damping coefficient(n-wall ) (kg/s) Wall friction coefficient (wall ) Gas-phase viscosity (g ) (kg m s) Gas-phase density (g ) (kg/m3 )
4000 0.08 0.2 27 30 2.25 × 10−5 4.50 × 10−4 0.001 2650 800 800 0.79 × 10−3 0.79 × 10−3 0.3 0.94 9.25 × 10−5 800 800 0.79 × 10−3 0.79 × 10−3 0.3 1.75 × 10−5 1.205
y- directions only (i.e., two-dimensional motion). The velocity change of the ith particle ( vi ) within a given solid-phase time step (ts ) is determined from the acceleration ai as follows: vi = ai ts .
(1)
The acceleration is determined using the Newtonian equation of motion ai =
Fi , m
(2)
where m is the particle mass and the net force of the ith particle (Fi ) is a sum of the forces induced by the fluid (Ff,i ); gravitational forces (Fg,i ); (repulsive) collisional contacts with other particles and/or the walls (Fcoll,i ); and (attractive) cohesive contacts with other particles and/or the walls (Fcoh,i ): Fi = Ff,i + Fg,i + Fcoll,i + Fcoh,i .
(3)
The first three of the forces appearing on the right-hand side of Eq. (3) are treated in a typical manner, as described in Appendix B. The following section describes the implementation of cohesive forces. 2.2. Cohesive forces (Fcoh,i ) Two different cohesion models are applied to the fluidizedbed simulation in this work, namely a “simple” square-well model and a more physical, though more complex, Hamaker model for van der Waals forces. The following sections outline the mathematical implementation of these models in this fluidized-bed simulation.
4513
2.2.1. Hamaker model for van der Waals forces The Hamaker theory is used to estimate the van der Waals force between two equal-sized, spherical particles according to the following equation (Israelachvili, 1991): FvdW ,ij =
Ar inner , 12Hij2
(4)
where rinner is the particle radius, A is the Hamaker constant, which is specific to a given material and has typical values on the order of 10−20 J, and Hij is surface-to-surface separation distance between particles i and j: Hij = (xi − xj )2 + (yi − yj )2 − 2rinner , (5) where xi and yi are the x- and y- positions of the ith particle. The corresponding cohesive force between a spherical particle and a flat wall is (Israelachvili, 1991) FvdW ,iw =
Ar inner , 2 6Hiw
(6)
where Hiw is the shortest surface-to-surface distance between the particle and the wall. For both of these expressions, the cohesive force approaches infinity as the separation distance approaches zero. This unrealistic singularity incurred at particle contact is avoided by introducing a “cutoff” distance, Hcut . For this work, a minimum cutoff distance of 4 × 10−10 m is used, which is based on the average interatomic distance for many common solids (Seville, 1997). For separation distances below this cutoff distance, the interparticle cohesive force is given by a surface adhesion force (Fad ) model (Seville, 1997): Fad = 2rinner ,
(7)
where is a constant surface energy per unit area, the value of which is calculated such that the van der Waals force is continuous at the cutoff distance. Thus, the cohesive force is maintained at a constant value for any separation distances (based on Eq. (5)) below the minimum cutoff distance, which includes any “negative” separation distances that occur during actual particle contact. The direction of the Hamaker force is along the normal unit vector connecting the particle centers. The total cohesive force is the sum of all attractive contacts experienced by the given particle in the current time step. This leads to the following expression for the total cohesive force on the ith particle if the Hamaker model is active: Fcoh,i =
n vdW
FvdW ,ij kij ,
(8)
j =1
where nvdw is the number of Hamaker interactions experienced by the ith particle and kij is the normal unit vector pointing from the center of the ith particle to the center of the j th particle. This treatment is depicted graphically in Fig. 1 along with the square-well model, which is described below. 2.2.2. Square-well model Similar to the Hamaker model for van der Waals forces, the square-well model is characterized by two input parameters,
4514
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
Hamaker model
Hamaker model
Force Forc
Square-well model
H 2router− 2rinner
Hcut (a)
H
Potential Energy
(to ∞)
(b)
D
Square-well model Hcut
2router− 2rinner
Fig. 1. (a) Cohesive force and (b) potential energy for Hamaker and square-well models as a function of surface-to-surface separation distance (H).
>2router
No Interaction
Agglomerate
2router
Formation 2rinner
2router
Capture-Cohesive Interaction
2router ApproachingCohesive Interaction
Inelastic Collision
Escaping-Cohesive Interactions
Fig. 2. Square-well interaction sequence.
namely the length scale router and the well depth D. Based on this treatment, cohesion is incorporated into the simulation as instantaneous forces whenever the separation distance is equal to 2router (where router > rinner ), as displayed in Fig. 1. A typical interaction sequence is shown in Fig. 2. The approachingcohesive interaction occurs when distant particles become separated by 2router . At this point, particles experience an instantaneous, attractive force (which corresponds to an increase in the normal component of the relative velocity between particles). Particles then move in a straight-line trajectory which leads to a collisional (repulsive) interaction if the particle separation distance reduces to 2rinner . Following the collisional interaction, particles travel in a straight path until the departingcohesive interaction, which also occurs at a separation distance of 2router . If the magnitude of the relative, normal velocity between particles is not large enough, an internal reflection of the particles about the well occurs, which is referred to as a capture-cohesive interaction. In this case, the trajectory is symmetric about the outer wall of the square well, which may lead to another dissipative collision. Presuming no collisions with a third particle occur, this process continues to repeat.
Alternatively, if a departing interaction between particles is characterized by a normal, relative velocity which is great enough to overcome the cohesive well, then an escapingcohesive interaction occurs at a separation distance of 2router . At this instant, the particles experience a reduction in their normal, relative velocity and continue to separate. Although the physical picture of the interaction sequence shown in Fig. 2 is limited to two particles (for purposes of simplicity), the square-well description also allows for agglomeration and de-agglomeration of more than two particles. Furthermore, the agglomerates formed via the square-well approach are dynamic in nature. Animations of the simulated system indicate that agglomerates continuously form, rearrange and break apart. It is worthwhile to note that the dynamics of agglomerates formed using the square-well potential (i.e., a “bouncing” between physical contact and the outer, cohesive well) are different than those of more elaborate models (like a Hamaker model for van der Waals forces), which are characterized by cohesive interactions in which particles “stick” for extended periods of time. However, the maximum separation distance
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
between particles in a square-well agglomerate is kept small (∼ 0.1% of particle radius) via the chosen router value. Furthermore, the normal component of relative velocity is reduced after each inelastic collision, and hence the motion between the agglomerated particles becomes more tangential in nature. As a result, the dynamic appearance of the square-well agglomerates closely mimics that of more elaborate models. Namely, the particles remain extremely close and appear to rotate around a common axis. A mathematical description of these cohesive interactions is given below. Changes in particle velocity due to square-well interactions are determined via solution of the momentum and kinetic energy balance equations for the given two particles. For cohesive interactions, kinetic energy can be gained, lost or stay the same depending on the type of cohesive interaction experienced. In particular, (i) potential energy is converted to kinetic energy (resulting in kinetic energy gain) for approaching-cohesive interactions, (ii) kinetic energy is converted to potential energy (resulting in kinetic energy loss) for escaping-cohesive interactions and (iii) kinetic energy remains the same for capturecohesive interactions. The corresponding momentum and energy balances are detailed elsewhere (Weber et al., 2004); only a summary of the final results is given below for the sake of brevity. For each type of interaction, solution of the relevant equations leads to an expression for the momentum impulse imparted to each particle during the interaction (Jcoh ). This momentum impulse is converted to a force (Fsw ) by dividing by the time step used for particle advancement (ts ): Fsw,ij = Jcoh,∗,ij /ts .
(9)
This attractive force acts along the direction of the unit vector (kij ) pointing from the center of the ith particle to the center of the j th particle leading to the following expression for the total cohesive force on the ith particle: Fcoh,i = Fsw,ij kij .
(10)
Because the kinetic energy change is different for each type of interaction, different expressions are obtained for the momentum impulse for each type of interaction. For an approachingcohesive interaction, this treatment leads to the relation: 4D m 2 Jcoh,app,ij = kij ,1 · vij ,a − + (kij ,1 · vij ,a ) , (11) 2 m
4515
where Jcoh,esc,ij is the magnitude of the momentum impulse transferred to each particle during the escaping cohesive interaction, kij ,3 is the unit vector pointing from the center of the ith particle to the center of the j th particle when the departing particles are separated by 2rinner and vij ,c is the relative velocity after the inelastic collision and prior to the cohesive interaction. However, if the value inside the square root in Eq. (12) is negative, particle escape from the attractive well is not possible. Instead, internal reflection from the wall of the square well occurs and the particles form a dynamic agglomerate. In this case, the normal, relative velocity between the two particles does not change in magnitude, but simply changes its direction. The magnitude of the momentum impulse imposed (Jcoh,cap ) on both particles during a capture cohesive interaction is given by Jcoh,cap,ij = m(kij ,3 · vij ,c ).
(13)
Particle–wall contacts are modeled like particle–particle contacts except that the wall is assigned an infinite mass (which comes into effect in Eqs. (11)–(13)). Furthermore, to account for the additional attractive force resulting from the interaction of a plane with a sphere rather than between two spheres, the well depth is doubled in all particle–wall interactions. The doubled cohesive-well depth at the wall is consistent with the Hamaker model of van der Waals forces in which the cohesive force between a sphere and a plane (Eq. (6)) is twice that of the force between two spheres (Eq. (4)). To identify square-well (instantaneous) interactions in the context of the soft-sphere algorithm (time-stepped), particles are checked for square-well width (router ) overlap at each time step. Once a square-well overlap is detected, an approaching cohesive interaction is implemented at that time step. For each particle, a list of overlapped particles is kept so that once a particle leaves the square well, an escaping-cohesive or capturecohesive interaction can be implemented in the form of a cohesive force at that time step. The exact separation distance at which interactions are implemented is determined by the end of the time step. However, the time step is set small enough that the extra distance traveled by the particles is typically an order of magnitude smaller than the distance between the particle surface and the outer edge of the square well (i.e., router − rinner ). 2.3. Model verification and validation
where Jcoh,app,ij is the magnitude of the momentum impulse transferred to each particle, m is the particle mass, vij ,a is the relative velocity of the particles ( vi − vj ) before the cohesive interaction and kij ,1 is the unit vector pointing from the center of the ith particle to the center of the j th particle when the approaching particles are separated by 2router . For the departing interaction, two possibilities exist. First, the particles can continue to separate. In this case, the following relation is obtained: m 4D 2 Jcoh,esc,ij = kij ,3 · vij ,c + (kij ,3 · vij ,c ) − , (12) 2 m
Various test cases have been examined for the purposes of model verification and validation. Verification refers to the process of ensuring numerical accuracy of the model solution, where as validation refers to the ability of the model to accurately predict system behavior (i.e, verification deals with numerics while validation deals with physics) (Grace and Taghipour, 2004). In this work, verification and validation efforts focused on the discrete-element portion of the simulation, since the fluid-phase portion of MFIX has been thoroughly tested and documented (McKeen and Pugsley, 2003; Syamlal and O’Brien, 2003, Xie et al., 2004).
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
2.3.1. Verification methods To ensure that the fluidized-bed results are relatively insensitive to numerical parameters, the time step associated with the dynamic solution of both phases and the size of the spatial grid were systematically varied. The parameters varied include the gas-phase time step (tf ), the solid-phase time step (ts ), the number of computational cells in the x-direction (mmax ), and the number of computational cells in the y-direction (nmax ). For all parameters investigated, the simulation results do not vary significantly as the settings were changed from the default values (Table 1); see Weber (2004) for details. The discrete-particle portion of the fluidized-bed model was verified by simulating two- particle collisions and comparing the post-collisional velocities with calculations performed in a spreadsheet. This procedure was carried out for non-cohesive particles, as well as for simulations using the square-well and Hamaker model. For the cohesive cases, three-particle simulations were also compared to spreadsheet calculations to ensure that agglomeration and breakup of two-particle agglomerates by a third particle are mathematically correct. 2.3.2. Model validation Both cohesive and non-cohesive simulations exhibited the pressure drop hysteresis that is frequently observed in experimental fluidized beds at the point of fluidization (see Fig. 3, which is detailed later, for a sketch of hysteretic behavior). Starting with the bed in a bubbling state (Us = 1.20 m/s), the gas velocity was slowly decreased in 37 equal steps to a superficial velocity of Us = 0.56 m/s over the course of 3.7 s (i.e., each gas velocity was maintained constant for 0.1 s). The superficial gas velocity (Us ) is calculated by summing the overall volumetric flow rate through the inlet jets and dividing this quantity by the cross-sectional area of the bed (Ab ) as follows: njets Us =
i=1 Ajet Vjet
Ab
,
(14)
where the volumetric flow rate through each jet is the specified gas velocity at that jet (Vjet ) multiplied by the cross-sectional area (Ajet ) of that jet. Three additional steps were used to take the superficial velocity to zero in 0.3 s. The inlet gas velocity was then gradually increased in stepped increments, using the reverse of the defluidization procedure. At each step, the gas velocity was maintained long enough to allow the bed to come to a statistical steady state (0.1 s). The steady-state pressure drop was then measured at each inlet velocity to compile a plot of pressure drop versus superficial velocity. The pressure drop across the bed was estimated by the pressure calculated at the lowest computational node of the fluidized-bed (i.e., the first computational node above the distributor plate) and averaged over time. The resulting pressure drop versus superficial gas velocity curve displayed a hysteresis over a wide range of conditions (see Weber, 2004 for details), thereby mimicking the experimentally observed behavior (Rhodes, 1998) and providing a qualitative benchmark for the simulation. An important quantitative benchmark is the minimum fluidization velocity. The minimum fluidization generally refers to
B
A
C
W Ab
∆P
4516
Us Fig. 3. Points identified as minimum fluidization velocity in previous studies.
the velocity above which further increases in the superficial velocity do not cause a further increase in the bed pressure drop. In previous efforts, however, some ambiguity exists in the exact procedure used to determine the minimum fluidization velocity. Fig. 3 depicts the typical defluidization and fluidization branches of the pressure drop (P ) versus superficial velocity (Us ) plot. A dotted line in the plot represents the weight of the particle bed (W ) divided by the bed cross-sectional area (Ab ). As the superficial velocity is increased through the packed bed regime, a point will be reached at which the pressure drop is equal to the weight of the particles (point A in Fig. 3). Some efforts have identified this velocity as the minimum fluidization velocity (Rhodes et al., 2001b; Xu et al., 2002) even though the particles are not fluidized at this point. As the superficial velocity is increased further, a point is reached at which the pressure overshoot is sufficient to free the particles from their packed configuration (point B in Fig. 3). This point has been identified as the minimum fluidization velocity by some authors (Rhodes, 1998), although the pressure drop will decrease slightly beyond this point as the particle bed expands. As the superficial velocity is decreased from the bubbling state, a velocity is reached where the pressure drop starts to decrease (point C in Fig. 3). At this point, the fluidized particles are gradually transitioning to the packed state. Some references have identified this point as the minimum fluidization velocity (Gera et al., 1998; Kawaguchi et al., 1998; Mikami et al., 1998; Wright and Raper, 1998ba). While all three of these velocities loosely satisfy the vague definition of the minimum fluidization velocity given earlier, the associated numerical values are different. For the purposes of this project, point A is referred to as the theoretical minimum fluidization velocity (Umf ,th ), point B is referred to as the incipient fluidization velocity (Umf ,in ), and point C is referred to as the defluidization velocity (Umf ,df ). The procedure for calculating these critical velocities involves graphical analysis of the pressure drop versus superficial velocity curve for increasing superficial velocity (Umf ,im and Umf ,th )
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Table 2 Non-cohesive critical velocities (m/s) Umf ,th
Umf ,in
Umf ,df
0.719
0.743
0.73
or decreasing (Umf ,df ) superficial velocity. The details of these procedures are described elsewhere (Weber, 2004). Correlations for the minimum fluidization velocity are welldocumented (Geankopolis, 1983; Rhodes, 1998). These correlations are developed by equating the weight of the particles in the bed (divided by the bed cross-sectional area) with empirical relations for the pressure drop in a packed bed such as the Ergun equation (Rhodes, 1998). This intersection corresponds to the theoretical minimum fluidization velocity (Umf ,th ) described above. As a means of validating the simulation used in this effort, predictions based on the Ergun equation were compared to values of the theoretical minimum fluidization velocity. These results, which are not presented here (Weber, 2004) show good agreement between the simulation results and the correlation predictions. From this point on, pressure drop results will be presented as a normalized pressure drop (P ∗ ), which is the pressure drop divided by the weight of particles above the lowest computational node normalized by the bed area (W /Ab ). In order to be consistent with the location of the pressure measurement only the weight of particles above the lowest computational node (W ) are considered. Furthermore, critical fluidization veloci∗ ) which is the ties will be presented in a normalized form (Umf given critical velocity divided by the corresponding velocity for the non-cohesive case. The values of the non-cohesive critical fluidization velocities used for normalization are presented in Table 2. 3. Parameter mapping from Hamaker model to square-well model As stated earlier, the impulsive nature of the square-well model makes it a viable option to incorporate cohesive forces into the kinetic theory, thereby affording it an advantage over more elaborate cohesion models (with non-instantaneous interactions). However, the square-well model has a major disadvantage in that no experimental data exists to characterize cohesive forces in terms of square-well parameters. In order to make the square-well model a viable, predictive tool, a method is needed to relate square-well parameters to tabulated cohesion properties (e.g., Hamaker constant A). A basic gauge of the level of particle cohesion is the magnitude of the cohesive force experienced by both particles. Unfortunately, using this force as a common ground to determine equivalent parameters for the square-well and Hamaker models is not plausible. As displayed in Fig. 1a, the square-well model exhibits an infinite force, whereas the Hamaker model is finite. Another possible mapping method is to set the parameters in each model such that the overall potential energy change in a given particle–particle interaction is the same for each model.
4517
When using the square-well model, the overall potential energy change is equal to the depth of the square well. However, the overall potential energy change induced by the Hamaker model is more complicated (see Fig. 1b). An analytical examination (Weber, 2004) of the collision dynamics associated with the Hamaker model shows that the overall change in potential energy is a function of the relative, normal velocity, whereas no such dependency exists for the square-well model. Hence, a mapping between the two sets of parameters would require the choice of a single, representative velocity, whereas the fluidized bed will exhibit a range of relative velocities between particles. A mapping technique which overcomes the aforementioned limitations is one based on the critical velocity needed for agglomerate formation. With both models, for given particle trajectories, a unique (incoming) relative velocity exists that demarcates the onset of escaping-cohesive interactions (no agglomeration) from capture-cohesive interactions (agglomeration). More specifically, for head-on (normal) collisions, the parameter that controls the magnitude of the cohesive force (A in the Hamaker model and D in the square-well model) dictates this critical “escape” velocity. As expected, higher levels of cohesion (higher A or D) result in higher escape velocities. The influence of length scale (rinner and Hcut ) on oblique collisions will be detailed in the results section. For the square-well model, the escape velocity can be calculated analytically. During an interaction between the ith and j th particles, the criterion for agglomerate formation is defined by the definition of the square-well combined with the kinetic energy and momentum balance for the two-particle system, resulting in the following inequality: 4D (kij ,3 · vij ,c )2 − 0, (15) m where kij ,3 is the unit vector between particle centers at the instant of the departing cohesive interaction, vij ,c is the relative velocity after the inelastic collision and m is the particle mass. This expression defines the minimum departing velocity necessary to avoid agglomeration. By combining Eq. (15) with the energy balance employed during an approaching cohesive interaction (see Eq. (11)), the minimum approach velocity for escape (vij ,esc ) in a head-on collision can be obtained as follows (Weber, 2004): 4D 1 vij ,esc = − 1 , (16) m e2 where e is the coefficient of restitution:
−0 e = exp q
(17)
and the intermediate variables are defined in terms of the normal spring constant Kn and damping coefficient n as follows: 0 = 2Kn /m, (18) n , 2mK n q = 0 1 − 2 .
= √
(19) (20)
4518
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
Table 3 Escape velocities obtained using normal collisions for both the square-well and Hamaker models D (J)
vij ,esc (m/s)
A (J)
vij ,esc (m/s)
1.00 × 10−11 5.00 × 10−11 1.00 × 10−10
1.79 × 10−3 4.00 × 10−3 5.66 × 10−3
3.34 × 10−19 7.49 × 10−19 1.06 × 10−18
1.79 × 10−3 4.00 × 10−3 5.66 × 10−3
This set of Eqs. (16)–(20) can be used to analytically calculate the escape velocity for a given square-well depth. For the Hamaker model, an analytical determination of the escape velocity is not possible. Nonetheless, the escape velocity can be obtained numerically via a series of two-particle simulations. In these simulations, two particles are separated by a finite distance in the x-direction with the same y-coordinates. Both particles are given a zero velocity in the y-direction with equal and opposite velocities in the x-direction. The simulation is allowed to proceed in time until a determination is made as to whether the particles escape or agglomerate. To determine the escape velocity, a series of these simulations are run and the x-component of velocity of each particle is gradually increased (or decreased) in each simulation until the transition from agglomeration to escaping (or vice versa) is achieved. For the purposes of this work, Eq. (16) is used to develop a table of escape velocities for various cohesive well depths. Twoparticle simulations are used to develop a table of escape velocities associated with different levels of cohesion with Hamaker model. A comparison between both tables is then used to match equivalent cohesion parameters from each model (D and A). As detailed in the next section, many-particle simulations are then run with both models. The mapped parameter values used in this assessment process are presented in Table 3. 4. Results and discussion The “escape-velocity” mapping is derived from the behavior of two-particle interactions. The mapping method is now tested in many-particle systems to determine if the square-well model and the Hamaker model produce the similar results on a bulk scale when using equivalent parameters. These comparisons serve not only to assess the effectiveness of the escape-velocity mapping, but also to determine if the square-well model is capable of predicting similar hydrodynamics as the more elaborate Hamaker model. The range of cohesive forces investigated is chosen to obtain the maximum cohesive forces possible, without exhibiting plugging and/or channeling. The useful range of parameters for each model that produce these conditions provides an initial means of testing the escape-velocity mapping method. For these simulations, router /rinner is kept constant at 1.01 and Hcut is kept constant at 4 × 10−10 m. Again, the latter choice is based on the known interatomic distance for common solids, whereas the impact of the former choice is detailed below.
4.1. Qualitative comparison: defluidization From a practical standpoint, the most undesirable consequence of cohesion in fluidized beds is the inability to fluidize in a conventional manner (homogenous or bubbling) due to the formation of plugs or channels. In terms of the Geldart classification (Geldart, 1973), this behavior is exhibited by Type C (“cohesive”, or relatively small) particles. The simulation developed here is capable of predicting this type of behavior with both cohesion models. With the square-well model, plugging and channeling occur at roughly D > 10−10 J, while simulations using the Hamaker model exhibit plugging and channeling at approximately A > 10−18 J. These parameters map to escape velocities that are roughly same (see Table 3). It is worthwhile to note, though, that the transition from Type A (fluidizable) to Type C (not fluidizable) behavior is not sufficiently abrupt to be used as a quantitative benchmark between the two cohesion models. While defluidization occurs at relatively similar levels of cohesion (i.e., mapped to the same escape velocity) for both the models investigated here, the nature of the defluidization is different. As revealed in Fig. 4a, the square-well model forms plugs that extend across the bed diameter as a single mass. Plugs formed with the Hamaker model are not as strong and tend to break apart into stable channels as illustrated in Fig. 4b. For higher levels of cohesion, simulations using the Hamaker model exhibit the formation of a completely static bed structure. This phenomenon is not observed with simulations using the square-well potential. Video images of fluidized-bed experiments (Rhodes, 2001) with Type C particles indicate that the solid-plugs predicted by the square-well model more closely approximate the experimental system. However, some experiments have shown the formation of stable, expanded bed structures due to cohesive effects (Mutsers and Rietema, 1977). So, although it is clear that the square-well model is capable of predicting defluidization, an assessment of its ability to predict the more subtle aspects of defluidization is inconclusive due to the discrepancies in previous experimental observations. As an aside, the level of cohesion considered here (A = 10−18 J) in conjunction with the particle diameter (Dp = 0.001 m) corresponds to a cohesive-force-to-particleweight ratio of approximately 20. Such high levels of van der Waals force are not typical with such a large particle diameter. Data provided by Molerus (1982) indicate that most particles which exhibit plugging (Geldart Type C) have diameters between 4.0 × 10−6 and 4.5 × 10−5 m with a cohesive force to particle weight ratio of approximately 1900. To use such a small diameter would require more particles than is computationally tractable or the use of a reduced system size. The use of unnaturally large levels of cohesion with large particles in fluidized-bed simulations is a common means of examining the effects of cohesion in a computationally efficient manner. Rhodes et al. (2001b) observed Type C behavior at cohesiveforce-to-particle-weight ratios greater than 46 (Dp = 0.001 m) with a simplified surface cohesion model. Ye et al. (2004) observed this transition at a ratio of 56 (Dp = 0.0001 m) with a van der Waals model. Therefore, while the cohesive force to
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
(a)
(b)
4519
t = 0.4 sec
t = 0.5 sec
t = 0.6 sec
t = 0.7 sec
U*s =
U*s =
U*s =
1.10
U*s = 1.18
0.94
1.02
t = 0.4 sec
t = 0.5 sec
t = 0.6 sec
t = 0.7 sec
U*s = 0.94
U*s = 1.02
U*s = 1.10
U*s = 1.18
Fig. 4. Simulation snapshots displaying plugging using the (a) square-well model (D = 1 × 10−10 J) and the (b) Hamaker model (A = 1 × 10−18 J).
4.2. Quantitative comparison: mixing
4
vij,esc= 0.00566 m /s (A =1.06x10-18 J)
3.5 3 Mixing Index
particle weight ratios used in this study are made greater than that expected for particles of the size considered here, the system mimics the behavior of fluidized beds composed of (smaller) particles in which cohesive forces play a nonnegligible role.
2.5
vij,esc= 0.00179 m/s (A =3.34x10-19 J)
2 1.5 1
I=
h1 , h2
0
0
10
20
30
Time (s)
(a) 4 3.5
vij,esc= 0.00566 m/s (D =1x10-10 J)
3 2.5
vij,esc= 0.00179 m/s (D =1x10-11 J)
2 1.5 1
No cohesion
0.5 0
(21)
0 (b)
n1 1 h1 = yi , n1
(22)
n2 1 yi , n2
(23)
i=1
h2 =
No cohesion
0.5
Mixing Index
In systems that do not exhibit plugging and channeling, mild levels of cohesion may have a significant effect on the mixing behavior. In order to further investigate the nature of mixing in both models, a mixing index is used to quantify the extent of mixing (Wu and Baeyens, 1998). The mixing index used here is defined by dividing all the particles in the bed into two “families”. Family 1 includes all the particles that have a height equal to or greater than the median particle height for the bed in its initial, packed state (i.e., top half of bed). Family 2 includes all the particles with a height less than the median particle height (i.e., bottom half of bed). The mixing index, I , is defined as the ratio of average height of all the particles in family 1, h1 , to the average height of all the particles in family 2, h2 , as defined below:
i=1
where n1 is the number particles in family 1, n2 is the number of particles in family 2, and yi is the height of the ith
10
20
30
Time (s)
Fig. 5. Evolution of mixing index in time for simulations using the (a) square-well cohesion model and (b) Hamaker model.
particle. For a system that remains segregated, the mixing index has a value of 3, while a perfectly mixed system exhibits a mixing index of 1. In Fig. 5a, the evolution of the mixing index over time is shown for several different levels of
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
cohesion using the square-well model. Fig. 5b displays similar results for simulations using the Hamaker model, with input parameters mapped according to the escape-velocity criterion. This figure also shows the escape velocity (vij ,esc ) associated with each level of cohesion. A comparison between the two figures indicates that the mapping leads to similar mixing characteristics in the fluidized bed.
1.12 1.10 square well Hamaker
1.08 U*mf,df
4520
1.06 1.04
4.3. Quantitative comparison: minimum fluidization velocity
1.02
Experimental results have shown that cohesive forces can lead to an increase in the minimum fluidization velocity (Passos and Mujumdar, 2000; Seville and Clift, 1984, Wright and Raper, 1998a,b). In the current section, the effects of cohesion are examined in terms of its influence on the three velocities associated with minimum fluidization defined previously, namely the theoretical minimum fluidization velocity (Umf ,th ), the incipient fluidization velocity (Umf ,in ) and the defluidization velocity (Umf ,df ). Both the square-well model and the Hamaker model are used to simulate defluidization and fluidization curves, with parameters chosen such that both models map to the same escape velocity. In these simulations, the fluidized-bed starts in the bubbling state and the defluidization cycle is implemented first, before the fluidization cycle. The effects of cohesion on the defluidization velocity, theoretical minimum fluidization velocity, and incipient fluidization velocity are displayed in Fig. 6 for simulations using both the Hamaker model and the square-well model. The simulation parameters outlined in Table 1 were used in all of these simulations. For each model, three levels of cohesion were implemented using the A and D values given in Table 3. In order to gauge the effectiveness of the mapping, results from both models are presented on the same plots as a function of the escape velocity. Three replicates were run for each parameter set and the standard deviation of the replicates is represented by error bars at each data point. Additional replicates did not lead to a significant decrease in this standard deviation. The parameter mapping works well for the defluidization velocity, as is displayed in Fig. 6a, with results from both cohesion models showing reasonable agreement at each escape velocity. Furthermore, experimental efforts by Wright and Raper (1998a,b) and simulation efforts by Mikami et al. (1998) and Kobayashi et al. (2003) all examined the behavior of the defluidization velocity (as defined previously). In each of these efforts an increase in the defluidization velocity was observed as the magnitude of cohesion was increased, which is consistent with the simulation results presented here (since an increase in escape velocity corresponds to an increase in the level of cohesion, as prescribed by A and D). Results from the increasing-velocity branch of the fluidization curve are presented in Figs. 6b and c. These plots illustrate the increase in both the theoretical and incipient fluidization velocities as the level of cohesion is increased in both models. However, the mapping exhibits some limitations at the higher levels of cohesion. The small discrepancies between the cohesion models for fluidization curve (Figs. 6b and c) are not surprising since the initial state for this data involves particle
1.00 0
0.002
0.004
0.006
0.004
0.006
vij,esc (m/s)
(a) 1.08 square well Hamaker
U*mf,th
1.06
1.04
1.02
1.00 0
0.002 vij,esc (m/s)
(b) 1.12 square well Hamaker
U*mf,in
1.08
1.04
1.00
0.96 0 (c)
0.002
0.004
0.006
vij,esc (m /s)
Fig. 6. (a) Defluidization velocity, (b) theoretical fluidization velocity, and (c) incipient fluidization velocity as a function of cohesive forces for both the Hamaker and square-well models.
contacts that are multi-particle and enduring in nature. Such enduring contacts are inherently accounted for by the Hamaker model but not by the square-well model (in which the contacts are treated as instantaneous and binary). On the other hand, for mixing, defluidization velocity, and average particle displacement (see next section), the system is fluidized at its initial state and the square-well model provides results very close to the Hamaker model at similar values of escape velocity. A comparison of these results to previous experimental results, unfortunately, is not straightforward. A review of both experimental and simulation studies reveals some ambiguity in
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
0.0005
0.0003
square well
0.0002
Hamaker no cohesion
avg
(m)
0.0004
∆l
the observed behavior of the critical fluidization/defluidization velocities. The general increase in the minimum fluidization velocity with an increase in the magnitude of cohesion is consistent with many experimental (Passos and Mujumdar, 2000; Seville and Clift, 1984, Wright and Raper, 1998a,b) and simulation (Kobayashi et al., 2003; Mikami et al., 1998) studies. These results agree with the results shown here. Unfortunately, very few efforts have provided enough detail to determine which critical velocity has been measured (theoretical or incipient). Therefore, it is difficult to comment on the validity of the results presented here for the theoretical minimum fluidization (Fig. 6b) and incipient fluidization (Fig. 6c).
0.0001 0
0
0.002
0.004 v
4.4. Quantitative comparison: average particle displacement To further evaluate the effectiveness of the escape-velocity mapping, particle displacement over a set time interval is examined, similar to that used by Wang and Rhodes (2004). An extreme example of small displacement is the plugging and channeling exhibited by Type C particles. On the other end of the spectrum is a vigorous, bubbling fluidization, in which net particle displacement is relatively high. The level of particle displacement is calculated by recording the distance moved by each particle (ti ) over the time increment t1 − t2 : li = (xt1 − xt2 )2 + (yt1 − yt2 )2 , (24) where xt1 and yt1 are the x- and y-positions of a given particle at time t1 , and xt2 and yt2 are the x- and y-positions of a given particle at time t2 . This quantity is then averaged over every particle in the simulation to obtain the average particle displacement (lavg ) in the given time increment: 1 li , lavg = n n
(25)
i=1
where n is the number of particles in the simulation. The values reported in Fig. 7 represent a cumulative average of 30 temporal data collection intervals that each lasts 0.3 s. These results were taken using the simulation parameters presented in Table 1 and the cohesion parameters presented in Table 3 with a constant superficial velocity of 0.77 m/s. Furthermore, the system was allowed to come to a steady state for 0.65 s before data was collected and three replicates were run for each case. Fig. 7 illustrates that the square-well model predicts average particle displacement similar to that of the Hamaker at similar escape velocities. 4.5. Effect of length scale of cohesive force One drawback of the escape-velocity mapping is that it does not account for the length scale of the square-well (router ) and Hamaker (Hcut ) models. For the square-well model, a change in the well-width-to-particle-radius ratio (router /rinner ) does not affect the escape velocity for head-on (normal) collisions (see Eq. (15)), but does change the escape velocity for oblique impacts. For the Hamaker model, changes in the cutoff distance
4521
ij,esc
0.006
(m /s)
Fig. 7. Average particle displacement with both cohesion models compared to escape velocity.
2
Approach angle
1
1
2
Fig. 8. Layout of two-particle simulations used to determine effect of square-well width on escape velocity.
(Hcut ) will affect the head-on escape velocity. Two approaches are taken to assess the effect of the characteristic length scale in each model. First, two-particle simulations are used to directly predict the escape velocity for various angles of approach (i.e., oblique collisions). Second, many-particle (4000) simulations are run using the square-well model with various values for router /rinner , and using the Hamaker model with various values for Hcut . The results of these simulations are then compared on the basis of predicted defluidization velocity and average particle displacement. Fig. 8 illustrates the configuration used to determine the effect of oblique collisions on the escape velocity. Fig. 9a reports the escape velocity determined via two-particle simulations using the square-well model with two different well depths and three different values of router /rinner . At larger router /rinner , the escape velocity varies somewhat, but overall the approach angle has little effect on the escape velocity except at very high approach angles. Similar results are obtained for the Hamaker model (with varied Hcut ) as displayed in Fig. 9b. As stated earlier, varying Hcut will affect the maximum Hamaker force and therefore affect the head-on escape velocity. To account for this variation, a different value of the Hamaker constant was used with each value of Hcut in order to achieve the same escape velocity for the head-on case (zero approach angle). The sets of parameters used are summarized in Table 4. While the
4522
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
0.01
router / rinner = 1.1
0.008
D=
1x10 -10 J
1.18
0.006 1.13 0.004 D = 1x10 -11 J 0.002 0
0
20
(a)
40 60 approach angle (degrees)
80
100
0.98 0.0001
Hcut = 4 x 10-11 m
0.04
Hcut = 4 x
1.08
1.03
0.05
vij,esc (m /s)
D=1x10-4 erg (vij,esc= 0.179 cm /s)
router / rinner = 1.001
U*mf, df
vij,esc (m /s)
D =1x10-3 erg (vij,esc= 0.566 cm/s)
router / rinner = 1.01
10-10 m
0.001
vij,esc= 0.00566 m/s
1.08
vij,esc= 0.00400 m/s vij,esc= 0.00179 m/s
0
(b)
20
40
60
80
100
approach angle (degrees)
U *mf, df
0
1
1.10 0.02 0.01
0.1
(router / rinner )-1
(a)
Hcut = 4 x 10-9 m
0.03
0.01
Fig. 9. Effect of approach angle on the escape velocity with (a) square-well model and (b) Hamaker model.
Table 4 Hamaker constants used to assess effect of cutoff distance at varying approach angle Hcut (m)
A (J)
vij ,esc (m/s)
4.0 × 10−9
1.04 × 10−16
0.0056 0.0056 0.0056 0.00179 0.00179 0.00179
4.0 × 10−10 4.0 × 10−11 4.0 × 10−9 4.0 × 10−10 4.0 × 10−11
1.06 × 10−18 1.06 × 10−20 3.15 × 10−19 3.34 × 10−19 3.36 × 10−21
escape velocity varies as the approach angle is increased (similar to the square-well model), the cutoff distance (Hcut ) has little very effect at each approach angle. For both cohesion models, the escape velocity deviates significantly from the “head-on” case only at very high collision angles. This deviation is partly due to the increased effects of friction as the particles develop a larger component of tangential relative velocity at higher approach angles. Two-particle simulations run without friction (results not shown) displayed a slight decrease in the escape velocity at very large approach angles. In addition to two-particle systems, the effects of length scale were also assessed on bulk parameters measured in many particle systems. The effects of router /rinner on the defluidization velocity are examined in Fig. 10a. At both well depths, the defluidization velocity changes very little as router /rinner is varied (if error bars are taken into account). Fig. 10b displays the effect of Hcut on the predicted defluidization velocity in simulations that utilize the Hamaker model.
1.06 1.04 1.02 1.00 1.E-11
(b)
1.E-10
1.E-09
1.E-08
1.E-07
Hcut (m)
Fig. 10. Effect of length scale on defluidization velocity: (a) effect of router /rinner for square-well model and (b) effect of Hcut for the Hamaker model.
In order to illustrate the effects of varying Hcut relative to varying the well depth (D), the results of simulations with higher and lower levels of cohesion are also presented. These results indicate the Hcut has little effect on the predicted defluidization velocity. 5. Summary Cohesive forces are known to be a continuous function of particle separation distance, and may result in enduring particle contacts (e.g., the formation of an agglomerate). Hence, the nature of cohesive forces conflicts with the assumption of binary, instantaneous interactions inherent in kinetic-theory (continuum) descriptions of particulate flows. In this work, a simplified model for cohesive forces, which has assumptions consistent with those of kinetic-theory treatments, has been assessed via comparison with more complex descriptions for cohesion. Specifically, a square-well model which implements cohesion as impulsive interactions between particles (thereby making it a viable method to incorporate cohesion into the kinetic theory) has been investigated. To gauge the effectiveness of the square-well model, cohesive forces have been
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
implemented into a discrete-particle fluidized-bed simulation using both the square-well model and the Hamaker model for van der Waals forces. Recall that the more complex Hamaker model is not suitable for incorporation into a continuum framework. Unlike the more elaborate Hamaker model, the square-well model is significantly limited by the lack of data characterizing real materials in terms of square-well parameters. Hence, a method to convert tabulated Hamaker constants into equivalent square-well parameters is required for the square-well model to be used as a predictive tool for particle flows. While there is no direct way to convert Hamaker constants to square-well depths, matching the escape velocity in a two-particle, headon collision was found to provide a good mapping between these two cohesion models. The effectiveness of this mapping technique was assessed by examining defluidization behavior, mixing behavior, fluidization velocities and average particle displacement. Wherever possible, simulation results from both models were also compared with experimental trends. Results of these simulations show that the square-well model makes predictions on par with those of the more elaborate Hamaker model. The incorporation of the square-well potential into a discrete-particle simulation is an important step in the possible implementation of cohesion into a kinetic-theory framework. Because the square-well model implements cohesion as impulsive interactions, it does not conflict with the assumption of binary, instantaneous interactions inherent in kinetic theory. Furthermore, the development of the escape-velocity mapping extends the usefulness of this model by providing a means to determine square-well parameters that correspond to material properties.
4523
where g is the void fraction, g is the gas density, vg is the gas velocity, S g is the gas-phase stress tensor, g is the acceleration due to gravity, and Igs is the rate of momentum transfer between the gas and solid phase per unit volume: vs − vg ), Igs = −s ∇Pg − Fgs (
(A.3)
where Pg is the gas-phase pressure, s is the solids volume fraction, and vs is the (locally averaged) solids velocity. The drag coefficient, Fgs , is determined using the expression developed by Symlal et al. (1994): Fgs =
3s g g 4Vt2 Dp
CDs | vs − vg |,
(A.4)
where Dp is the particle diameter and CDs is the single-sphere drag coefficient given by Dalla Valle (1948): 2 Vt CDs = 0.63 + 4.8 . (A.5) Re The terminal velocity, Vt , is given by the following correlation developed by Garside and Al-Dibouni (1977): Vt = 0.5 (A0 − 0.06Re 2 2 + (0.06Re) + 0.12Re(2B − A) + A ,
(A.6)
where A0 = 4.14 g
0.81.28 , 0.85, B = 2.65g g , g > 0.85
(A.7) (A.8)
and the Reynolds number, Re, is defined as Acknowledgments
Re =
The authors are grateful to the National Science Foundation for providing funding support for this work under Grant no. CTS 0226010. The authors would also like to thank Dr. Tom O’Brien and Dr. Dhanunjay S. Boyalakuntla for helpful discussions during the course of this work. Appendix A. Computational approach for (continuum) gas phase The multiphase flow with interphase exchanges (MFIX) program developed at the Department of Energy National Energy Technology Laboratory (DOE NETL) was used as a framework to model the gas–fluid system in this work (www.mfix.org). The fluid phase is modeled by solving locally averaged mass and momentum balances (Syamlal et al., 1994): j (g g ) + ∇ · (g g vg ) = 0, jt
(A.1)
j (g g vg ) + ∇ · (g g vg vg ) = ∇ · S g + g g g − Igs , (A.2) jt
Dp | vs − vg |g g
,
(A.9)
where g is the gas viscosity. For the purposes of calculating the void fraction, the twodimensional simulation domain is assumed to have a depth of one particle diameter (Hoomans et al., 1996; Kawaguchi et al., 1998; Tsuji et al., 1993; Tsuji, 2000, van Wachem et al., 2001; Xu and Yu, 1997). In each computational cell, the solids volume fraction is calculated as follows: nmn s,mn = i=1 , (A.10) Vmn where is the volume of the ith particle, s,mn is the solids volume fraction, Vmn is the volume of the grid cell and nmn is the total number of particles in the m, n grid cell. The average solids velocity is calculated by averaging over all the particles in a given grid cell: nmn vi vs,mn = i=1 , (A.11) nmn where vs,mn is the average solids velocity in grid cell m, n, and vi is the velocity of particle i.
4524
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
It should be noted that the method of calculating the void fraction (Xu and Yu, 1997) presented in Eq. (A.10) introduces some error since the simulation under consideration is twodimensional whereas drag correlations (which depend on the void fraction; see Eq. (A.4)) are based on three-dimensional systems. Equal numbers of particles will pack differently in two- and three-dimensional domains and therefore produce different void fractions. A discussion of the issues related to this mismatch, including possible methods to convert a twodimensional void fraction into an equivalent three-dimensional void fraction, is provided by van Wachem et al. (2001). In particular, they found that different strategies for converting the two-dimensional void fraction to a three-dimensional void fraction lead to differences in pressure fluctuations (bubble dynamics) and bed expansion. Such methods are not employed in this work, since the various methods are found to have little effect on the predicted minimum fluidization velocities (as was the focus in this work), which match closely with those obtained from empirical relations (Weber, 2004). The system of equations governing the gas-phase are solved numerically using a modified version of the SIMPLE method developed by Patankar (1980). For more details on the solution of the gas-phase equations, the reader is referred to the MFIX Numerical Theory Guide (Syamlal, 1998). Appendix B. Computational approach for (discrete) particle phase The three forces acting on a particle which are not described in the main text include the fluid-induced forces (Ff,i ), gravitational forces (Fg,i ), and forces arising from (repulsive) collisional contacts with other particles and/or the walls (Fcoll,i ). Fluid-particle forces (Ff,i ) are calculated using Eq. (A.4). The forces induced by the fluid are equal and opposite to the forces exerted by the particles on the fluid. For each computational cell, the total force imposed by the particles on the fluid is divided by the number of particles to determine the force imposed on each particle by the fluid (Boyalakuntla, 2003). Gravitational forces (Fg,i ) are calculated based on the volume of a spherical particle and the acceleration due to gravity: 4 3 Fg,i = rinner s g, 3
(B.1)
where g is the acceleration due to gravity. The repulsive, collisional force acting on particle i is composed of a normal and tangential component Fcoll,i = Fn,ij kij + Ft,ij kt,ij , (B.2) where Fn,ij and Ft,ij are the magnitudes of the normal and tangential components of the repulsive force when the ith and j th particles collide. The directions of these forces are along the line of the normal unit vector pointing from the center of the ith particle to the center of the j th particle (kij ) and the tangential unit vector (kt,ij ), which is perpendicular to the normal vector with the cross product (kij × kn,ij ) pointing out of the page.
The summations represent all collisions experienced by the ith particle during the given time step. The magnitude of the forces is calculated using a linear-spring-and-dashpot model (Tsuji et al., 1993) characterized by normal (Kn ) and tangential (Kt ) spring constants, as well as dissipation constants (n , t .). The basic equation for calculating the normal force between two particles is given by Fn,ij = − n,ij Kn − n
j n,ij , jt
(B.3)
where n,ij is the normal component of the overlap between particles i and j . The normal overlap can be calculated as the difference between the sum of the particle radii (2rinner ) and the distance between particle centers: n,ij = 2rinner − (xi − xj )2 + (yi − yj )2 , (B.4) where xi and yi are the x- and y- components of the position of the ith particle, and xj and yj are the x- and y- components of the position of the j th particle. The tangential component of the contact force is calculated as follows: Ft,ij = − t,ij Kt − t vt,ij ,
(B.5)
where Kt is the tangential spring constant and t is the tangential damping coefficient. In this expression, vt,ij is the tangential slip velocity between the two particles: vt,ij = vij · kt,ij − rinner (i + j ),
(B.6)
where i and j are the rotation rates of the ith and j th particles (positive is counter-clockwise), respectively, and rinner is the particle radius. Also, t,ij is the tangential overlap, defined as follows: t,ij = ts · vt,ij .
(B.7)
Tangential interactions are calculated as either sticking or sliding contacts based on the following inequality: |Ft,ij | |Fn,ij |,
(B.8)
where is the coefficient of friction. If this inequality is not satisfied, then the particles undergo a sliding contact and the tangential component of the collision force is calculated as Ft,ij = Fn,ij .
(B.9)
To reduce computational times, the spring constant is made somewhat softer than realistic materials (as is often done; see, for example, Tsuji et al., 1993; Gera et al., 1998; Kawaguchi et al., 1998). The damping coefficient is still set, though, to achieve a realistic coefficient of restitution (e = 0.94). Several researchers (Gera et al., 1998; Kawaguchi et al., 1998; Kuwagi et al., 2000) have indicated that the use of an artificially small spring constant does not have a significant effect on the nature of the flow in fluidized-bed simulations. Nonetheless, the results of the softened particles were spot-checked against simulations using spring constants that led to more realistic collision times to ensure that physical accuracy was not compromised. The spring constant does have a slight effect on the calculated
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
minimum fluidization velocity, but the pressure drop at fluidization is not affected. The observed effect on the minimum fluidization velocity is similar for both cohesion models (squarewell and Hamaker), and thus does not impact the comparison between the two models. Since an initial goal of the development phase was to reproduce the results of Rhodes et al. (2001a), the default parameters in Table 1 closely mimic the conditions in that effort.
Fcoll,i Ff,i Fg,i Fgs Fn,ij Fsw,ij Ft,ij FvdW ,ij g H1 H2 Hcut Hij Hiw I Igs Jcoh,∗ Jcoh,app,ij
Jcoh,esc,ij
kij kt,ij kij ,1 kij ,2 kij ,3
Notation ai A Ab Ajet D Dp e Fad Fi Fcoh,i
Jcoh,cap,ij
acceleration of ith particle, m/s2 Hamaker constant, J fluidized-bed cross-sectional area, m2 inlet jet cross-sectional area, m2 depth of square well, J particle diameter, m coefficient of restitution adhesion force, N sum of forces on ith particle, N cohesive interpaticle forces imposed on ith particle, N sum of forces from collisional contacts on ith particle, N forces induced by fluid on ith particle, N gravitational forces on ith particle, N drag coefficient, kg/m3 /s normal component of collisional force between ith and j th particle, N magnitude of force imposed on ith and j th particles from square-well interaction, N tangent component of collisional force between ith and j th particle, N cohesive van der Waals forces between the ith and j th particles, N accelerationdue to gravity, m/s2 average height of particles in family one, m average height of particles in family two, m cutoff distance for Hamaker van der Waals model, m surface-to-surface separation distance between ith and j th particles, m surface-to-surface distance between ith particle and wall, m mixing index rate of momentum transfer between gas and solid phase, N/m3 momentum impulse transferred between particles during square-well interaction, kg m/s momentum impulse transferred during approaching-cohesive interaction between ith and j th particles, kg m/s momentum impulse transferred during escaping-cohesive interaction between ith and j th particles, kg m/s
Kn Kt lavg li m mmax n nmax nmn nvdW P P ∗
Pg rinner router, Sg t ts tf Us Us∗ Umf ,df ∗ Umf ,dr Umf ,in ∗ Umf ,in Umf ,th ∗ Umf ,in vg vi vij ,a
4525
momentum impulse transferred during capturecohesive interaction between ith and j th particles, kg m/s normal unit vector pointing from center of the ith particle to center of the j th particle tangent unit vector normal unit vector at instance of approaching cohesive interaction normal unit vector at instance of physical contact normal unit vector at instance of departing cohesiveinteraction normal spring constant for particle–particle contacts (stiffness), N/m tangential spring constant for particle–particle contacts, N/m average particle movement, m distance moved by ith particle during current time step, m particle mass, kg number of computational cells in x-direction total number of particles in simulation number of computational cells in y-direction number of particles in grid cell m, n number of van der Waals contacts experienced by a particle at a given time step pressure drop across fluidized bed, Pa pressure drop across fluidized-bed normalized by weight of particles above lowest computational node, Pa gas-phase pressure, Pa particle radius, m width of square well, m gas-phase stress tensor, N/m2 time, s solid-phase time step, s gas-phase time step, s superficial velocity in fluidized bed, m/s normalized superficial gas velocity (Us /(Umf ,th,non-cohesive )) defluidization velocity, m/s defluidization velocity normalized by defluidization velocity for non-cohesive case incipient fluidization velocity,m/s incipient fluidization velocity normalized by incipient fluidization velocity for non-cohesive case theoretical minimum fluidization velocity, m/s theoretical minimum fluidization velocity normalized by theoretical minimum fluidization velocity for non-cohesive case gas velocity, m/s velocity of ith particle, m/s volume of ith particle, m3 relative velocity before approaching cohesive interaction, m/s
4526
vij ,c vij ,esc vs vt,ij vi Vjet Vt Vmn W W xi yi
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
relative velocity after inelastic collision, m/s minimum normal, relative velocity needed to escape cohesive forces, m/s locally averaged solid velocity, m/s slip velocity between ith and j th particles, m/s velocity change of the ith particle in a given time step, m/s specified as velocity at inlet jet, m/s terminal velocity, m/s volume of grid cell m, n, m3 total weight of particles in fluidized bed, kg total weight of particles above lowestcomputational node, kg x-position of ith particle, m y-position of ith particle, m
Greek letters n,ij t,ij g n t g g s i
surface energy, J normal component of overlap between ith and j th particles, m tangential component of overlap between ith and j th particles, m void fraction normal damping coefficient for particle–particle contacts, N s/m tangential damping coefficient for particle–particle contacts, N s/m gas viscosity, kgs1 m1 friction coefficient gas density, kg/m3 particle density, kg/m3 rotation rate of ith particle, 1/s
References Allen, M.P., Tildesley, D.J., 1987. Computer Simulation of Liquids. Oxford University Press, Oxford. Boyalakuntla, D., 2003. Simulation of granular and gas–solid flows using discrete element method. Ph.D. Thesis, Carnegie Mellon University, www.mfix.org. Dalla Valle, J.M., 1948. Micrometrics. Pitman, London. Forsyth, A.J., Hutton, S., Rhodes, M.J., 2002. Effect of cohesive interparticle force on the flow characteristics of granular material. Powder Technology 126, 150–154. Garside, J., Al-Dibouni, M.R., 1977. Velocity-voidage relationships for fluidization and sedimentation in solid–liquid systems. Industrial and Engineering Chemistry Process Design and Development 16, 206–214. Geankopolis, C.J., 1983. Transport Processes and Unit Operations. Allyn and Bacon, Boston. Geldart, D., 1973. Types of gas fluidization. Powder Technology 7, 285–292. Gera, D., Gautam, M., Tsuji, Y., Kawaguchi, T., Tanaka, Y., 1998. Computer simulation of bubbles in large-particle fluidized beds. Powder Technology 98, 38–47. Grace, J.R., Taghipour, F., 2004. Verification and validation of CFD models and dynamic similarity for fluidized beds. Powder Technology 139, 99–110. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach. Chemical Engineering Science 51, 99–118.
Israelachvili, J., 1991. Intermolecular and Surface Forces. Academic Press London, London. Kawaguchi, T., Tanaka, T., Tsuji, T., 1998. Numerical simulation of twodimensional fluidized beds using the discrete element method (comparison between two- and three-dimensional models). Powder Technology 96, 129–138. Kim, H., Arastoopour, H., 2002. Extension of kinetic theory to cohesive particle flow. Powder Technology 122, 83–94. Kobayashi, T., Mukai, T., Kawaguchi, T., Tanaka, T., Tsuji, Y., 2003. DEM analysis of flow patterns of Geldart’s group A particles in a fluidized bed. World Congress on Particle Technology, vol. 4. Sydney, Australia. Kuwagi, K., Mikami, T., Horio, M., 2000. Numerical simulation of metallic solid bridging particles in a fluidized bed at high temperature. Powder Technology 109, 27–40. McKeen, T., Pugsley, T., 2003. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Technology 129, 139–152. McLaughlin, L.J., Rhodes, M.J., 2001. Prediction of fluidized bed behaviour in the presence of liquid bridges. Powder Technology 114, 213–223. Mikami, T., Kamiya, H., Horio, M., 1998. Numerical simulation of cohesive powder behavior in a fluidized bed. Chemical Engineering Science 53, 1927–1940. Molerus, O., 1982. Interpretation of Geldart’s type A, B, C, and D powders by taking into account interparticle cohesion forces. Powder Technology 33, 81–87. Mutsers, S.M.P., Rietema, K., 1977. The effect of interparticle forces on the expansion of a homogeneous gas-fluidized bed. Powder Technology 18, 239–248. Nase, S.T., Vargas, W.L., Abatan, A.A., McCarthy, J.J., 2001. Discrete characterization tools for cohesive granular material. Powder Technology 116, 214–223. Passos, M.L., Mujumdar, A.S., 2000. Effect of cohesive forces on fluidized and spouted beds of wet particles. Powder Technology 110, 222–238. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere, New York. Rhodes, M., 1998. Introduction to Particle Technology. Wiley, West Sussex, England. Rhodes, M.J., 2001. Fluidization of particles by fluids. www.erpt.org/ 012Q/rhod-04.htm. Rhodes, M.J., Wang, X.S., Nguyen, M., Stewart, P., Liffman, K., 2001a. Use of discrete element method simulation in studying fluidization characteristics: influence of interpaticle force. Chemical Engineering Science 56, 69–76. Rhodes, M.J., Wang, X.S., Nguyen, M., Stewart, P., Liffman, K., 2001b. Onset of cohesive behavior in gas fluidized beds: a numerical study using DEM simulation. Chemical Engineering Science 56, 4433–4438. Seville, J.P.K., 1997. Processing of Particulate Solids. Blackie Academic and Professional, London. Seville, J.P.K., Clift, R., 1984. The effect of thin liquid layers on fluidization characteristics. Powder Technology 37, 117–129. Syamlal, M., 1998. MFIX Documentation Numerical Technique. US Department of Energy, Morgantown, West Virginia, www.mfix.org. Syamlal, M., O’Brien, T.J., 2003. Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed. A.I.Ch.E. Journal 49, 2793–2801. Syamlal, M., Rogers, W., O’Brien, T.J., 1994. MFIX Documentation Theory Guide. US Department of Energy, Morgantown, West Virginia, www.mfix.org. Tatemoto, Y., Mawatari, Y., Noda, K., 2005. Numerical simulation of cohesive particle motion in vibrated fluidized bed. Chemical Engineering Science 60, 5010–5021. Tsinontides, S.C., Jackson, R., 1993. The mechanics of gas fluidized beds with an interval of stable fluidization. Journal of Fluid Mechanics 255, 237–274. Tsuji, Y., 2000. Activities in discrete particle simulation in Japan. Powder Technology 113, 278–286. Tsuji, Y., Kawaguchi, T., Tanaka, T., 1993. Discrete particle simulation of two-dimensional fluidized bed. Powder Technology 77, 79–87. van Wachem, B.G.M., van der Schaaf, J., Schouten, J.C., Krishna, R., van den Bleek, C.M., 2001. Experimental validation of Lagrangian–Eulerian simulations of fluidized beds. Powder Technology 116, 155–165.
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Wang, X.S., Rhodes, M.J., 2004. Mechanistic study of defluidization by numerical simulation. Chemical Engineering Science 59, 215–222. Weber, M.W., 2004. Simulation of cohesive particle flows in granular and gas–solid systems. Ph.D. Thesis, University of Colorado, www.mfix.org. Weber, M.W., Hoffman, D.K., Hrenya, C.M., 2004. Discrete-particle simulations of cohesive granular flow using a square-well potential. Granular Matter 6, 239–254. Wright, P.C., Raper, J.A., 1998a. Role of liquid bridge forces in cohesive fluidization. Transactions of the Institutions of Chemical Engineers 76, 753–760. Wright, P.C., Raper, J.A., 1998b. Examination of dispersed liquid-phase threephase fluidized beds: part 1 Non-porous, uniform particle systems. Powder Technology 97, 208–226.
4527
Wu, S.Y., Baeyens, J., 1998. Segregation by size difference in gas fluidized beds. Powder Technology 98, 139–150. Xie, N., Battaglia, F., Fox, R.O., 2004. Simulations of multiphase reactive flows in fluidized beds using in situ adaptive tabulation. Combustion Theory and Modeling 8, 195–209. Xu, B.H., Yu, A.B., 1997. Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science 52, 2785–2809. Xu, B.H., Zhou, Y.C., Yu, A.B., Zulli, P., 2002. Force structures in gas fluidized beds of fine powders. World congress on particle technology, vol. 4. Sydney, Australia. Ye, M., van der Hoef, M.A, Kuipers, J.A.M., 2004. A numerical study of fluidization behavior of Geldart A particles using a discrete particle model. Powder Technology 139, 129–139.
Square-well model for cohesion in fluidized beds Michael W. Weber, Christine M. Hrenya ∗ Department of Chemical and Biological Engineering, University of Colorado, Boulder, CO 80309, USA Received 23 September 2005; received in revised form 10 December 2005; accepted 5 February 2006 Available online 8 May 2006
Abstract Cohesive forces are implemented into a discrete-particle, fluidized-bed simulation using a square-well potential. The square-well description treats cohesive interactions as instantaneous, binary events, thereby making it a viable option for the incorporation of cohesion into a kinetictheory-based continuum model. Cohesive forces are also incorporated into the simulation using the more elaborate Hamaker description of van der Waals forces in order to provide a basis for assessing the square-well model. Both cohesion models are implemented in the discrete-particle framework of the MFIX software package. A mapping method is also developed to convert material-specific Hamaker constants into equivalent square-well parameters. The corresponding results from the two models are compared both qualitatively and quantitatively. The predictions of the square-well model are on par with the Hamaker model with respect to mixing level, particle mobility and minimum fluidization velocity. Subtle differences are observed between the two models in cases that involved such high levels of cohesion that the particle bed could not fully fluidize. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Fluidization; Mathematical modeling; Multiphase flow; Particle; Cohesion
1. Introduction An important consideration in fluidized-bed operation is the potential impact of cohesive forces between particles. The presence of cohesion can lead to inefficient operation due to reduced gas–solid contacting and can also result in undesired behaviors such as plugging and channeling. Cohesive forces may arise from a variety of sources including liquid bridges, van der Waals forces, electrostatic forces and magnetic forces. Mathematical models (Mikami et al., 1998; Rhodes et al., 2001a,b; Tatemoto et al., 2005) for cohesive systems provide an important complement to experimental investigations since experimentally controlling the level of cohesion, independent of changing other system conditions, is non-trivial. For example, particle cohesion has been controlled experimentally by adding a liquid to the system such as water, oil or surfactant to control the volume of the liquid bridges formed between particles (Forsyth et al., 2002; McLaughlin and Rhodes, 2001; Nase et al., 2001). However, such changes may have secondary ∗ Corresponding author. Tel.: +1 303 492 7689; fax: +1 303 492 4341.
E-mail address: [email protected] (C.M. Hrenya). 0009-2509/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.02.008
effects on other properties. For example, humidification of the air changes the density and viscosity of the fluidizing gas, as well as the mass of the particles due to adsorbed water (Tsinontides and Jackson, 1993). Computational models, on the other hand, allow for independent control of cohesive forces. In models for fluidized-bed systems, the fluid is treated as a continuum and thus a single mass and momentum balance is used for this phase. With regard to the solid phase, both discrete and continuum treatments have been utilized. The discrete approach (Lagrangian) involves the solution of a force balance for every particle present in the system, whereas the continuum approach (Eulerian) involves the solution of a single mass and momentum balance. For systems of industrial size, the Eulerian approach is more computationally efficient than the Lagrangian approach, which is currently limited to approximately 100,000 particles for desktop machines. However, the incorporation of complex particle physics (e.g., cohesion) is a more difficult task with Eulerian models. In particular, the averaging technique used to derive the continuum balances gives rise to terms, like the solid-phase stress tensor, which require a constitutive relation. The incorporation of cohesion into such constitutive quantities is considerably less straightforward than its incorporation into a discrete-particle model.
4512
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
To date, Eulerian models that include the effects of cohesion are extremely limited. A kinetic-theory-based description requires that particle interactions are binary and instantaneous in nature, whereas cohesive forces are known to vary as a function of separation distance and often result in enduring contacts. Kim and Arastoopour (2002) used the kinetic theory framework to account for cohesion via a contact bonding energy. Their model accounts for the formation of spherical particle agglomerates but does not allow for deagglomeration. A less rigorous approach was taken by McKeen and Pugsley (2003) who developed an Eulerian model for the fluidization of cohesive, fluid-catalytic-cracking particles in which cohesion is accounted for by increasing the effective particle diameter used in the drag-force relation. Efforts to include cohesion in Lagrangian models are more numerous since, as mentioned earlier, the incorporation of such physics is more straightforward. Mikami et al. (1998) incorporated a liquid-bridge model into a fluidized-bed simulation. This effort produced pressure data illustrating an increase in the minimum fluidization velocity for wet particles relative to a dry system. Kuwagi et al. (2000) applied a sintering model and studied its effects on the formation of dead zones in the fluidized bed. Rhodes et al. (2001a,b) used the same fluidized-bed simulation developed by Mikami et al., but incorporated a simpler cohesion model. Instead of the liquid-bridge model, they devised a cohesion model in which cohesive forces are only implemented when the particles are in physical contact. The magnitude of the cohesive force is set to a constant and equal to an integer multiple of the particle weight. More recently, Wang and Rhodes (2004) developed new parameters to quantify the flowability of the fluidized bed. Using parameters such as the particle mobility and the average particle speed, they defined quantitative criteria for the transition from Geldart Type A to Type C particles. Xu et al. (2002) implemented a model for van der Waals forces into a fluidized-bed simulation. The results of this effort show the formation of static force chains within the fluidized bed. Ye et al. (2004) incorporated van der Waals forces into a soft-sphere discrete-particle fluidized-bed model using the Hamaker theory. They observed the formation of channels reminiscent of Geldart Type C behavior for high levels of cohesion. Tatemoto et al. (2005) studied the effects of cohesion in a discrete-element simulation of vibrated fluidized beds. They examined the effect of vibration frequency and strength on the pressure losses in the fluidized bed. To summarize, numerous Lagrangian models for cohesive particles in fluidized beds have been shown to compare well with experimental observations. Eulerian models for cohesion are relatively scarce, and their development has been hampered by the challenge of incorporating cohesion into a framework that requires binary and instantaneous interactions between particles, whereas cohesive interactions are often enduring in nature. For example, the Hamaker description of van der Waals is characterized by enduring contacts. With these ideas in mind, a primary goal of this effort is to investigate an alternative description of cohesive forces, namely the squarewell potential. The square-well description treats cohesion as arising from impulsive events between two particles, and thus
represents a possible means of incorporating cohesion into continuum (kinetic-theory) models. This model has previously been applied to granular flows (no interstitial fluid) where it was shown to produce qualitatively correct predictions (Weber et al., 2004). The performance of the square-well approximation for cohesive forces in fluidized beds is investigated here via discrete-particle simulations. For purposes of comparison, cohesive forces are also implemented in discrete-particle simulations using the more complex Hamaker description of van der Waals forces (recall that such enduring contacts cannot be incorporated directly into kinetic theory models). A mapping method is developed to determine the value of square-well inputs given the physical parameters required by the Hamaker description. The resulting discrete-particle simulations indicate the effectiveness of the simpler, square-well model relative to the more elaborate Hamaker model. 2. Computational method The simulation domain consists of a two-dimensional rectangular system with the origin of the x- and y- axis in the bottom left corner. Spheres are restricted to motion in a twodimensional plane (i.e., thickness of fluidized bed is a single particle diameter). The bed typically used is 0.08 m wide by 0.2 m tall. For this length and width, 27 computational cells in the x-direction and 30 computational cells in the y-direction are used. In the following text, any changes from this bed size are noted. The governing equations for the complete fluidizedbed simulation are composed of fluid-phase and particle-phase balances. The fluid-phase continuum equations are handled in a standard manner, the details of which are given in Appendix A. The boundary conditions for the gas phase consist of no-slip, impermeable walls on the vertical sides of the bed. For the outflow boundary condition at the top of the bed, zero gauge pressure is specified across the entire width. Thirteen separate inlet jets are used to specify the inlet boundary conditions. The boundaries of each inlet jet are required to coincide with the boundary of a computational cell, limiting the range of sizes that can be used for the inlet jets. In most cases, each jet consists of one computational cell (2.96 × 10−4 m in width). At each inlet jet, the void fraction is set to unity and the gas velocity is set to achieve the desired superficial velocity in the bed. This inlet jet configuration provides the most uniform fluidization with the least amount of dead space (unfluidized particles) in the bed. Between each jet, the distributor plate is represented by impermeable, no-slip walls. The particle-phase equations are described below. For further detail on the computational method, the reader is referred to Weber (2004). Table 1 includes a listing of all parameters used in the simulation, unless otherwise noted. 2.1. Particle advancement The discrete-particle simulation operates using a timestepped algorithm (Allen and Tildesley, 1987). The particles are modeled as soft spheres confined to movement in the x- and
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Table 1 Input parameters used in fluidized bed simulation Parameter
Value
Number of particles (n) Fluidized-bed width (xmax ) (m) Fluidized-bed height (ymax ) (m) Computational cells in x-direction (mmax ) Computational cells in y-direction (nmax ) Solid-phase time step (ts ) (s) Gas-phase time step (tf ) (s) Particle diameter (DP ) (m) Particle density (s ) (kg/m3 ) Normal spring constant (Kn ) (kg/s2 ) Tangential spring constant (Kt ) (kg/s2 ) Normal damping coefficient (n ) (kg/s) Tangential damping coefficient (n ) (kg/s) Friction coefficient () Coefficient of restitution (e) Collision duration (tc ) (s) Normal wall spring constant (Kn-wall ) (kg/s2 ) Tangential wall spring constant (Kt -wall ) (kg/s2 ) Normal wall damping coefficient (n-wall ) (kg/s) Tangential wall damping coefficient(n-wall ) (kg/s) Wall friction coefficient (wall ) Gas-phase viscosity (g ) (kg m s) Gas-phase density (g ) (kg/m3 )
4000 0.08 0.2 27 30 2.25 × 10−5 4.50 × 10−4 0.001 2650 800 800 0.79 × 10−3 0.79 × 10−3 0.3 0.94 9.25 × 10−5 800 800 0.79 × 10−3 0.79 × 10−3 0.3 1.75 × 10−5 1.205
y- directions only (i.e., two-dimensional motion). The velocity change of the ith particle ( vi ) within a given solid-phase time step (ts ) is determined from the acceleration ai as follows: vi = ai ts .
(1)
The acceleration is determined using the Newtonian equation of motion ai =
Fi , m
(2)
where m is the particle mass and the net force of the ith particle (Fi ) is a sum of the forces induced by the fluid (Ff,i ); gravitational forces (Fg,i ); (repulsive) collisional contacts with other particles and/or the walls (Fcoll,i ); and (attractive) cohesive contacts with other particles and/or the walls (Fcoh,i ): Fi = Ff,i + Fg,i + Fcoll,i + Fcoh,i .
(3)
The first three of the forces appearing on the right-hand side of Eq. (3) are treated in a typical manner, as described in Appendix B. The following section describes the implementation of cohesive forces. 2.2. Cohesive forces (Fcoh,i ) Two different cohesion models are applied to the fluidizedbed simulation in this work, namely a “simple” square-well model and a more physical, though more complex, Hamaker model for van der Waals forces. The following sections outline the mathematical implementation of these models in this fluidized-bed simulation.
4513
2.2.1. Hamaker model for van der Waals forces The Hamaker theory is used to estimate the van der Waals force between two equal-sized, spherical particles according to the following equation (Israelachvili, 1991): FvdW ,ij =
Ar inner , 12Hij2
(4)
where rinner is the particle radius, A is the Hamaker constant, which is specific to a given material and has typical values on the order of 10−20 J, and Hij is surface-to-surface separation distance between particles i and j: Hij = (xi − xj )2 + (yi − yj )2 − 2rinner , (5) where xi and yi are the x- and y- positions of the ith particle. The corresponding cohesive force between a spherical particle and a flat wall is (Israelachvili, 1991) FvdW ,iw =
Ar inner , 2 6Hiw
(6)
where Hiw is the shortest surface-to-surface distance between the particle and the wall. For both of these expressions, the cohesive force approaches infinity as the separation distance approaches zero. This unrealistic singularity incurred at particle contact is avoided by introducing a “cutoff” distance, Hcut . For this work, a minimum cutoff distance of 4 × 10−10 m is used, which is based on the average interatomic distance for many common solids (Seville, 1997). For separation distances below this cutoff distance, the interparticle cohesive force is given by a surface adhesion force (Fad ) model (Seville, 1997): Fad = 2rinner ,
(7)
where is a constant surface energy per unit area, the value of which is calculated such that the van der Waals force is continuous at the cutoff distance. Thus, the cohesive force is maintained at a constant value for any separation distances (based on Eq. (5)) below the minimum cutoff distance, which includes any “negative” separation distances that occur during actual particle contact. The direction of the Hamaker force is along the normal unit vector connecting the particle centers. The total cohesive force is the sum of all attractive contacts experienced by the given particle in the current time step. This leads to the following expression for the total cohesive force on the ith particle if the Hamaker model is active: Fcoh,i =
n vdW
FvdW ,ij kij ,
(8)
j =1
where nvdw is the number of Hamaker interactions experienced by the ith particle and kij is the normal unit vector pointing from the center of the ith particle to the center of the j th particle. This treatment is depicted graphically in Fig. 1 along with the square-well model, which is described below. 2.2.2. Square-well model Similar to the Hamaker model for van der Waals forces, the square-well model is characterized by two input parameters,
4514
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
Hamaker model
Hamaker model
Force Forc
Square-well model
H 2router− 2rinner
Hcut (a)
H
Potential Energy
(to ∞)
(b)
D
Square-well model Hcut
2router− 2rinner
Fig. 1. (a) Cohesive force and (b) potential energy for Hamaker and square-well models as a function of surface-to-surface separation distance (H).
>2router
No Interaction
Agglomerate
2router
Formation 2rinner
2router
Capture-Cohesive Interaction
2router ApproachingCohesive Interaction
Inelastic Collision
Escaping-Cohesive Interactions
Fig. 2. Square-well interaction sequence.
namely the length scale router and the well depth D. Based on this treatment, cohesion is incorporated into the simulation as instantaneous forces whenever the separation distance is equal to 2router (where router > rinner ), as displayed in Fig. 1. A typical interaction sequence is shown in Fig. 2. The approachingcohesive interaction occurs when distant particles become separated by 2router . At this point, particles experience an instantaneous, attractive force (which corresponds to an increase in the normal component of the relative velocity between particles). Particles then move in a straight-line trajectory which leads to a collisional (repulsive) interaction if the particle separation distance reduces to 2rinner . Following the collisional interaction, particles travel in a straight path until the departingcohesive interaction, which also occurs at a separation distance of 2router . If the magnitude of the relative, normal velocity between particles is not large enough, an internal reflection of the particles about the well occurs, which is referred to as a capture-cohesive interaction. In this case, the trajectory is symmetric about the outer wall of the square well, which may lead to another dissipative collision. Presuming no collisions with a third particle occur, this process continues to repeat.
Alternatively, if a departing interaction between particles is characterized by a normal, relative velocity which is great enough to overcome the cohesive well, then an escapingcohesive interaction occurs at a separation distance of 2router . At this instant, the particles experience a reduction in their normal, relative velocity and continue to separate. Although the physical picture of the interaction sequence shown in Fig. 2 is limited to two particles (for purposes of simplicity), the square-well description also allows for agglomeration and de-agglomeration of more than two particles. Furthermore, the agglomerates formed via the square-well approach are dynamic in nature. Animations of the simulated system indicate that agglomerates continuously form, rearrange and break apart. It is worthwhile to note that the dynamics of agglomerates formed using the square-well potential (i.e., a “bouncing” between physical contact and the outer, cohesive well) are different than those of more elaborate models (like a Hamaker model for van der Waals forces), which are characterized by cohesive interactions in which particles “stick” for extended periods of time. However, the maximum separation distance
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
between particles in a square-well agglomerate is kept small (∼ 0.1% of particle radius) via the chosen router value. Furthermore, the normal component of relative velocity is reduced after each inelastic collision, and hence the motion between the agglomerated particles becomes more tangential in nature. As a result, the dynamic appearance of the square-well agglomerates closely mimics that of more elaborate models. Namely, the particles remain extremely close and appear to rotate around a common axis. A mathematical description of these cohesive interactions is given below. Changes in particle velocity due to square-well interactions are determined via solution of the momentum and kinetic energy balance equations for the given two particles. For cohesive interactions, kinetic energy can be gained, lost or stay the same depending on the type of cohesive interaction experienced. In particular, (i) potential energy is converted to kinetic energy (resulting in kinetic energy gain) for approaching-cohesive interactions, (ii) kinetic energy is converted to potential energy (resulting in kinetic energy loss) for escaping-cohesive interactions and (iii) kinetic energy remains the same for capturecohesive interactions. The corresponding momentum and energy balances are detailed elsewhere (Weber et al., 2004); only a summary of the final results is given below for the sake of brevity. For each type of interaction, solution of the relevant equations leads to an expression for the momentum impulse imparted to each particle during the interaction (Jcoh ). This momentum impulse is converted to a force (Fsw ) by dividing by the time step used for particle advancement (ts ): Fsw,ij = Jcoh,∗,ij /ts .
(9)
This attractive force acts along the direction of the unit vector (kij ) pointing from the center of the ith particle to the center of the j th particle leading to the following expression for the total cohesive force on the ith particle: Fcoh,i = Fsw,ij kij .
(10)
Because the kinetic energy change is different for each type of interaction, different expressions are obtained for the momentum impulse for each type of interaction. For an approachingcohesive interaction, this treatment leads to the relation: 4D m 2 Jcoh,app,ij = kij ,1 · vij ,a − + (kij ,1 · vij ,a ) , (11) 2 m
4515
where Jcoh,esc,ij is the magnitude of the momentum impulse transferred to each particle during the escaping cohesive interaction, kij ,3 is the unit vector pointing from the center of the ith particle to the center of the j th particle when the departing particles are separated by 2rinner and vij ,c is the relative velocity after the inelastic collision and prior to the cohesive interaction. However, if the value inside the square root in Eq. (12) is negative, particle escape from the attractive well is not possible. Instead, internal reflection from the wall of the square well occurs and the particles form a dynamic agglomerate. In this case, the normal, relative velocity between the two particles does not change in magnitude, but simply changes its direction. The magnitude of the momentum impulse imposed (Jcoh,cap ) on both particles during a capture cohesive interaction is given by Jcoh,cap,ij = m(kij ,3 · vij ,c ).
(13)
Particle–wall contacts are modeled like particle–particle contacts except that the wall is assigned an infinite mass (which comes into effect in Eqs. (11)–(13)). Furthermore, to account for the additional attractive force resulting from the interaction of a plane with a sphere rather than between two spheres, the well depth is doubled in all particle–wall interactions. The doubled cohesive-well depth at the wall is consistent with the Hamaker model of van der Waals forces in which the cohesive force between a sphere and a plane (Eq. (6)) is twice that of the force between two spheres (Eq. (4)). To identify square-well (instantaneous) interactions in the context of the soft-sphere algorithm (time-stepped), particles are checked for square-well width (router ) overlap at each time step. Once a square-well overlap is detected, an approaching cohesive interaction is implemented at that time step. For each particle, a list of overlapped particles is kept so that once a particle leaves the square well, an escaping-cohesive or capturecohesive interaction can be implemented in the form of a cohesive force at that time step. The exact separation distance at which interactions are implemented is determined by the end of the time step. However, the time step is set small enough that the extra distance traveled by the particles is typically an order of magnitude smaller than the distance between the particle surface and the outer edge of the square well (i.e., router − rinner ). 2.3. Model verification and validation
where Jcoh,app,ij is the magnitude of the momentum impulse transferred to each particle, m is the particle mass, vij ,a is the relative velocity of the particles ( vi − vj ) before the cohesive interaction and kij ,1 is the unit vector pointing from the center of the ith particle to the center of the j th particle when the approaching particles are separated by 2router . For the departing interaction, two possibilities exist. First, the particles can continue to separate. In this case, the following relation is obtained: m 4D 2 Jcoh,esc,ij = kij ,3 · vij ,c + (kij ,3 · vij ,c ) − , (12) 2 m
Various test cases have been examined for the purposes of model verification and validation. Verification refers to the process of ensuring numerical accuracy of the model solution, where as validation refers to the ability of the model to accurately predict system behavior (i.e, verification deals with numerics while validation deals with physics) (Grace and Taghipour, 2004). In this work, verification and validation efforts focused on the discrete-element portion of the simulation, since the fluid-phase portion of MFIX has been thoroughly tested and documented (McKeen and Pugsley, 2003; Syamlal and O’Brien, 2003, Xie et al., 2004).
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
2.3.1. Verification methods To ensure that the fluidized-bed results are relatively insensitive to numerical parameters, the time step associated with the dynamic solution of both phases and the size of the spatial grid were systematically varied. The parameters varied include the gas-phase time step (tf ), the solid-phase time step (ts ), the number of computational cells in the x-direction (mmax ), and the number of computational cells in the y-direction (nmax ). For all parameters investigated, the simulation results do not vary significantly as the settings were changed from the default values (Table 1); see Weber (2004) for details. The discrete-particle portion of the fluidized-bed model was verified by simulating two- particle collisions and comparing the post-collisional velocities with calculations performed in a spreadsheet. This procedure was carried out for non-cohesive particles, as well as for simulations using the square-well and Hamaker model. For the cohesive cases, three-particle simulations were also compared to spreadsheet calculations to ensure that agglomeration and breakup of two-particle agglomerates by a third particle are mathematically correct. 2.3.2. Model validation Both cohesive and non-cohesive simulations exhibited the pressure drop hysteresis that is frequently observed in experimental fluidized beds at the point of fluidization (see Fig. 3, which is detailed later, for a sketch of hysteretic behavior). Starting with the bed in a bubbling state (Us = 1.20 m/s), the gas velocity was slowly decreased in 37 equal steps to a superficial velocity of Us = 0.56 m/s over the course of 3.7 s (i.e., each gas velocity was maintained constant for 0.1 s). The superficial gas velocity (Us ) is calculated by summing the overall volumetric flow rate through the inlet jets and dividing this quantity by the cross-sectional area of the bed (Ab ) as follows: njets Us =
i=1 Ajet Vjet
Ab
,
(14)
where the volumetric flow rate through each jet is the specified gas velocity at that jet (Vjet ) multiplied by the cross-sectional area (Ajet ) of that jet. Three additional steps were used to take the superficial velocity to zero in 0.3 s. The inlet gas velocity was then gradually increased in stepped increments, using the reverse of the defluidization procedure. At each step, the gas velocity was maintained long enough to allow the bed to come to a statistical steady state (0.1 s). The steady-state pressure drop was then measured at each inlet velocity to compile a plot of pressure drop versus superficial velocity. The pressure drop across the bed was estimated by the pressure calculated at the lowest computational node of the fluidized-bed (i.e., the first computational node above the distributor plate) and averaged over time. The resulting pressure drop versus superficial gas velocity curve displayed a hysteresis over a wide range of conditions (see Weber, 2004 for details), thereby mimicking the experimentally observed behavior (Rhodes, 1998) and providing a qualitative benchmark for the simulation. An important quantitative benchmark is the minimum fluidization velocity. The minimum fluidization generally refers to
B
A
C
W Ab
∆P
4516
Us Fig. 3. Points identified as minimum fluidization velocity in previous studies.
the velocity above which further increases in the superficial velocity do not cause a further increase in the bed pressure drop. In previous efforts, however, some ambiguity exists in the exact procedure used to determine the minimum fluidization velocity. Fig. 3 depicts the typical defluidization and fluidization branches of the pressure drop (P ) versus superficial velocity (Us ) plot. A dotted line in the plot represents the weight of the particle bed (W ) divided by the bed cross-sectional area (Ab ). As the superficial velocity is increased through the packed bed regime, a point will be reached at which the pressure drop is equal to the weight of the particles (point A in Fig. 3). Some efforts have identified this velocity as the minimum fluidization velocity (Rhodes et al., 2001b; Xu et al., 2002) even though the particles are not fluidized at this point. As the superficial velocity is increased further, a point is reached at which the pressure overshoot is sufficient to free the particles from their packed configuration (point B in Fig. 3). This point has been identified as the minimum fluidization velocity by some authors (Rhodes, 1998), although the pressure drop will decrease slightly beyond this point as the particle bed expands. As the superficial velocity is decreased from the bubbling state, a velocity is reached where the pressure drop starts to decrease (point C in Fig. 3). At this point, the fluidized particles are gradually transitioning to the packed state. Some references have identified this point as the minimum fluidization velocity (Gera et al., 1998; Kawaguchi et al., 1998; Mikami et al., 1998; Wright and Raper, 1998ba). While all three of these velocities loosely satisfy the vague definition of the minimum fluidization velocity given earlier, the associated numerical values are different. For the purposes of this project, point A is referred to as the theoretical minimum fluidization velocity (Umf ,th ), point B is referred to as the incipient fluidization velocity (Umf ,in ), and point C is referred to as the defluidization velocity (Umf ,df ). The procedure for calculating these critical velocities involves graphical analysis of the pressure drop versus superficial velocity curve for increasing superficial velocity (Umf ,im and Umf ,th )
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Table 2 Non-cohesive critical velocities (m/s) Umf ,th
Umf ,in
Umf ,df
0.719
0.743
0.73
or decreasing (Umf ,df ) superficial velocity. The details of these procedures are described elsewhere (Weber, 2004). Correlations for the minimum fluidization velocity are welldocumented (Geankopolis, 1983; Rhodes, 1998). These correlations are developed by equating the weight of the particles in the bed (divided by the bed cross-sectional area) with empirical relations for the pressure drop in a packed bed such as the Ergun equation (Rhodes, 1998). This intersection corresponds to the theoretical minimum fluidization velocity (Umf ,th ) described above. As a means of validating the simulation used in this effort, predictions based on the Ergun equation were compared to values of the theoretical minimum fluidization velocity. These results, which are not presented here (Weber, 2004) show good agreement between the simulation results and the correlation predictions. From this point on, pressure drop results will be presented as a normalized pressure drop (P ∗ ), which is the pressure drop divided by the weight of particles above the lowest computational node normalized by the bed area (W /Ab ). In order to be consistent with the location of the pressure measurement only the weight of particles above the lowest computational node (W ) are considered. Furthermore, critical fluidization veloci∗ ) which is the ties will be presented in a normalized form (Umf given critical velocity divided by the corresponding velocity for the non-cohesive case. The values of the non-cohesive critical fluidization velocities used for normalization are presented in Table 2. 3. Parameter mapping from Hamaker model to square-well model As stated earlier, the impulsive nature of the square-well model makes it a viable option to incorporate cohesive forces into the kinetic theory, thereby affording it an advantage over more elaborate cohesion models (with non-instantaneous interactions). However, the square-well model has a major disadvantage in that no experimental data exists to characterize cohesive forces in terms of square-well parameters. In order to make the square-well model a viable, predictive tool, a method is needed to relate square-well parameters to tabulated cohesion properties (e.g., Hamaker constant A). A basic gauge of the level of particle cohesion is the magnitude of the cohesive force experienced by both particles. Unfortunately, using this force as a common ground to determine equivalent parameters for the square-well and Hamaker models is not plausible. As displayed in Fig. 1a, the square-well model exhibits an infinite force, whereas the Hamaker model is finite. Another possible mapping method is to set the parameters in each model such that the overall potential energy change in a given particle–particle interaction is the same for each model.
4517
When using the square-well model, the overall potential energy change is equal to the depth of the square well. However, the overall potential energy change induced by the Hamaker model is more complicated (see Fig. 1b). An analytical examination (Weber, 2004) of the collision dynamics associated with the Hamaker model shows that the overall change in potential energy is a function of the relative, normal velocity, whereas no such dependency exists for the square-well model. Hence, a mapping between the two sets of parameters would require the choice of a single, representative velocity, whereas the fluidized bed will exhibit a range of relative velocities between particles. A mapping technique which overcomes the aforementioned limitations is one based on the critical velocity needed for agglomerate formation. With both models, for given particle trajectories, a unique (incoming) relative velocity exists that demarcates the onset of escaping-cohesive interactions (no agglomeration) from capture-cohesive interactions (agglomeration). More specifically, for head-on (normal) collisions, the parameter that controls the magnitude of the cohesive force (A in the Hamaker model and D in the square-well model) dictates this critical “escape” velocity. As expected, higher levels of cohesion (higher A or D) result in higher escape velocities. The influence of length scale (rinner and Hcut ) on oblique collisions will be detailed in the results section. For the square-well model, the escape velocity can be calculated analytically. During an interaction between the ith and j th particles, the criterion for agglomerate formation is defined by the definition of the square-well combined with the kinetic energy and momentum balance for the two-particle system, resulting in the following inequality: 4D (kij ,3 · vij ,c )2 − 0, (15) m where kij ,3 is the unit vector between particle centers at the instant of the departing cohesive interaction, vij ,c is the relative velocity after the inelastic collision and m is the particle mass. This expression defines the minimum departing velocity necessary to avoid agglomeration. By combining Eq. (15) with the energy balance employed during an approaching cohesive interaction (see Eq. (11)), the minimum approach velocity for escape (vij ,esc ) in a head-on collision can be obtained as follows (Weber, 2004): 4D 1 vij ,esc = − 1 , (16) m e2 where e is the coefficient of restitution:
−0 e = exp q
(17)
and the intermediate variables are defined in terms of the normal spring constant Kn and damping coefficient n as follows: 0 = 2Kn /m, (18) n , 2mK n q = 0 1 − 2 .
= √
(19) (20)
4518
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
Table 3 Escape velocities obtained using normal collisions for both the square-well and Hamaker models D (J)
vij ,esc (m/s)
A (J)
vij ,esc (m/s)
1.00 × 10−11 5.00 × 10−11 1.00 × 10−10
1.79 × 10−3 4.00 × 10−3 5.66 × 10−3
3.34 × 10−19 7.49 × 10−19 1.06 × 10−18
1.79 × 10−3 4.00 × 10−3 5.66 × 10−3
This set of Eqs. (16)–(20) can be used to analytically calculate the escape velocity for a given square-well depth. For the Hamaker model, an analytical determination of the escape velocity is not possible. Nonetheless, the escape velocity can be obtained numerically via a series of two-particle simulations. In these simulations, two particles are separated by a finite distance in the x-direction with the same y-coordinates. Both particles are given a zero velocity in the y-direction with equal and opposite velocities in the x-direction. The simulation is allowed to proceed in time until a determination is made as to whether the particles escape or agglomerate. To determine the escape velocity, a series of these simulations are run and the x-component of velocity of each particle is gradually increased (or decreased) in each simulation until the transition from agglomeration to escaping (or vice versa) is achieved. For the purposes of this work, Eq. (16) is used to develop a table of escape velocities for various cohesive well depths. Twoparticle simulations are used to develop a table of escape velocities associated with different levels of cohesion with Hamaker model. A comparison between both tables is then used to match equivalent cohesion parameters from each model (D and A). As detailed in the next section, many-particle simulations are then run with both models. The mapped parameter values used in this assessment process are presented in Table 3. 4. Results and discussion The “escape-velocity” mapping is derived from the behavior of two-particle interactions. The mapping method is now tested in many-particle systems to determine if the square-well model and the Hamaker model produce the similar results on a bulk scale when using equivalent parameters. These comparisons serve not only to assess the effectiveness of the escape-velocity mapping, but also to determine if the square-well model is capable of predicting similar hydrodynamics as the more elaborate Hamaker model. The range of cohesive forces investigated is chosen to obtain the maximum cohesive forces possible, without exhibiting plugging and/or channeling. The useful range of parameters for each model that produce these conditions provides an initial means of testing the escape-velocity mapping method. For these simulations, router /rinner is kept constant at 1.01 and Hcut is kept constant at 4 × 10−10 m. Again, the latter choice is based on the known interatomic distance for common solids, whereas the impact of the former choice is detailed below.
4.1. Qualitative comparison: defluidization From a practical standpoint, the most undesirable consequence of cohesion in fluidized beds is the inability to fluidize in a conventional manner (homogenous or bubbling) due to the formation of plugs or channels. In terms of the Geldart classification (Geldart, 1973), this behavior is exhibited by Type C (“cohesive”, or relatively small) particles. The simulation developed here is capable of predicting this type of behavior with both cohesion models. With the square-well model, plugging and channeling occur at roughly D > 10−10 J, while simulations using the Hamaker model exhibit plugging and channeling at approximately A > 10−18 J. These parameters map to escape velocities that are roughly same (see Table 3). It is worthwhile to note, though, that the transition from Type A (fluidizable) to Type C (not fluidizable) behavior is not sufficiently abrupt to be used as a quantitative benchmark between the two cohesion models. While defluidization occurs at relatively similar levels of cohesion (i.e., mapped to the same escape velocity) for both the models investigated here, the nature of the defluidization is different. As revealed in Fig. 4a, the square-well model forms plugs that extend across the bed diameter as a single mass. Plugs formed with the Hamaker model are not as strong and tend to break apart into stable channels as illustrated in Fig. 4b. For higher levels of cohesion, simulations using the Hamaker model exhibit the formation of a completely static bed structure. This phenomenon is not observed with simulations using the square-well potential. Video images of fluidized-bed experiments (Rhodes, 2001) with Type C particles indicate that the solid-plugs predicted by the square-well model more closely approximate the experimental system. However, some experiments have shown the formation of stable, expanded bed structures due to cohesive effects (Mutsers and Rietema, 1977). So, although it is clear that the square-well model is capable of predicting defluidization, an assessment of its ability to predict the more subtle aspects of defluidization is inconclusive due to the discrepancies in previous experimental observations. As an aside, the level of cohesion considered here (A = 10−18 J) in conjunction with the particle diameter (Dp = 0.001 m) corresponds to a cohesive-force-to-particleweight ratio of approximately 20. Such high levels of van der Waals force are not typical with such a large particle diameter. Data provided by Molerus (1982) indicate that most particles which exhibit plugging (Geldart Type C) have diameters between 4.0 × 10−6 and 4.5 × 10−5 m with a cohesive force to particle weight ratio of approximately 1900. To use such a small diameter would require more particles than is computationally tractable or the use of a reduced system size. The use of unnaturally large levels of cohesion with large particles in fluidized-bed simulations is a common means of examining the effects of cohesion in a computationally efficient manner. Rhodes et al. (2001b) observed Type C behavior at cohesiveforce-to-particle-weight ratios greater than 46 (Dp = 0.001 m) with a simplified surface cohesion model. Ye et al. (2004) observed this transition at a ratio of 56 (Dp = 0.0001 m) with a van der Waals model. Therefore, while the cohesive force to
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
(a)
(b)
4519
t = 0.4 sec
t = 0.5 sec
t = 0.6 sec
t = 0.7 sec
U*s =
U*s =
U*s =
1.10
U*s = 1.18
0.94
1.02
t = 0.4 sec
t = 0.5 sec
t = 0.6 sec
t = 0.7 sec
U*s = 0.94
U*s = 1.02
U*s = 1.10
U*s = 1.18
Fig. 4. Simulation snapshots displaying plugging using the (a) square-well model (D = 1 × 10−10 J) and the (b) Hamaker model (A = 1 × 10−18 J).
4.2. Quantitative comparison: mixing
4
vij,esc= 0.00566 m /s (A =1.06x10-18 J)
3.5 3 Mixing Index
particle weight ratios used in this study are made greater than that expected for particles of the size considered here, the system mimics the behavior of fluidized beds composed of (smaller) particles in which cohesive forces play a nonnegligible role.
2.5
vij,esc= 0.00179 m/s (A =3.34x10-19 J)
2 1.5 1
I=
h1 , h2
0
0
10
20
30
Time (s)
(a) 4 3.5
vij,esc= 0.00566 m/s (D =1x10-10 J)
3 2.5
vij,esc= 0.00179 m/s (D =1x10-11 J)
2 1.5 1
No cohesion
0.5 0
(21)
0 (b)
n1 1 h1 = yi , n1
(22)
n2 1 yi , n2
(23)
i=1
h2 =
No cohesion
0.5
Mixing Index
In systems that do not exhibit plugging and channeling, mild levels of cohesion may have a significant effect on the mixing behavior. In order to further investigate the nature of mixing in both models, a mixing index is used to quantify the extent of mixing (Wu and Baeyens, 1998). The mixing index used here is defined by dividing all the particles in the bed into two “families”. Family 1 includes all the particles that have a height equal to or greater than the median particle height for the bed in its initial, packed state (i.e., top half of bed). Family 2 includes all the particles with a height less than the median particle height (i.e., bottom half of bed). The mixing index, I , is defined as the ratio of average height of all the particles in family 1, h1 , to the average height of all the particles in family 2, h2 , as defined below:
i=1
where n1 is the number particles in family 1, n2 is the number of particles in family 2, and yi is the height of the ith
10
20
30
Time (s)
Fig. 5. Evolution of mixing index in time for simulations using the (a) square-well cohesion model and (b) Hamaker model.
particle. For a system that remains segregated, the mixing index has a value of 3, while a perfectly mixed system exhibits a mixing index of 1. In Fig. 5a, the evolution of the mixing index over time is shown for several different levels of
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
cohesion using the square-well model. Fig. 5b displays similar results for simulations using the Hamaker model, with input parameters mapped according to the escape-velocity criterion. This figure also shows the escape velocity (vij ,esc ) associated with each level of cohesion. A comparison between the two figures indicates that the mapping leads to similar mixing characteristics in the fluidized bed.
1.12 1.10 square well Hamaker
1.08 U*mf,df
4520
1.06 1.04
4.3. Quantitative comparison: minimum fluidization velocity
1.02
Experimental results have shown that cohesive forces can lead to an increase in the minimum fluidization velocity (Passos and Mujumdar, 2000; Seville and Clift, 1984, Wright and Raper, 1998a,b). In the current section, the effects of cohesion are examined in terms of its influence on the three velocities associated with minimum fluidization defined previously, namely the theoretical minimum fluidization velocity (Umf ,th ), the incipient fluidization velocity (Umf ,in ) and the defluidization velocity (Umf ,df ). Both the square-well model and the Hamaker model are used to simulate defluidization and fluidization curves, with parameters chosen such that both models map to the same escape velocity. In these simulations, the fluidized-bed starts in the bubbling state and the defluidization cycle is implemented first, before the fluidization cycle. The effects of cohesion on the defluidization velocity, theoretical minimum fluidization velocity, and incipient fluidization velocity are displayed in Fig. 6 for simulations using both the Hamaker model and the square-well model. The simulation parameters outlined in Table 1 were used in all of these simulations. For each model, three levels of cohesion were implemented using the A and D values given in Table 3. In order to gauge the effectiveness of the mapping, results from both models are presented on the same plots as a function of the escape velocity. Three replicates were run for each parameter set and the standard deviation of the replicates is represented by error bars at each data point. Additional replicates did not lead to a significant decrease in this standard deviation. The parameter mapping works well for the defluidization velocity, as is displayed in Fig. 6a, with results from both cohesion models showing reasonable agreement at each escape velocity. Furthermore, experimental efforts by Wright and Raper (1998a,b) and simulation efforts by Mikami et al. (1998) and Kobayashi et al. (2003) all examined the behavior of the defluidization velocity (as defined previously). In each of these efforts an increase in the defluidization velocity was observed as the magnitude of cohesion was increased, which is consistent with the simulation results presented here (since an increase in escape velocity corresponds to an increase in the level of cohesion, as prescribed by A and D). Results from the increasing-velocity branch of the fluidization curve are presented in Figs. 6b and c. These plots illustrate the increase in both the theoretical and incipient fluidization velocities as the level of cohesion is increased in both models. However, the mapping exhibits some limitations at the higher levels of cohesion. The small discrepancies between the cohesion models for fluidization curve (Figs. 6b and c) are not surprising since the initial state for this data involves particle
1.00 0
0.002
0.004
0.006
0.004
0.006
vij,esc (m/s)
(a) 1.08 square well Hamaker
U*mf,th
1.06
1.04
1.02
1.00 0
0.002 vij,esc (m/s)
(b) 1.12 square well Hamaker
U*mf,in
1.08
1.04
1.00
0.96 0 (c)
0.002
0.004
0.006
vij,esc (m /s)
Fig. 6. (a) Defluidization velocity, (b) theoretical fluidization velocity, and (c) incipient fluidization velocity as a function of cohesive forces for both the Hamaker and square-well models.
contacts that are multi-particle and enduring in nature. Such enduring contacts are inherently accounted for by the Hamaker model but not by the square-well model (in which the contacts are treated as instantaneous and binary). On the other hand, for mixing, defluidization velocity, and average particle displacement (see next section), the system is fluidized at its initial state and the square-well model provides results very close to the Hamaker model at similar values of escape velocity. A comparison of these results to previous experimental results, unfortunately, is not straightforward. A review of both experimental and simulation studies reveals some ambiguity in
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
0.0005
0.0003
square well
0.0002
Hamaker no cohesion
avg
(m)
0.0004
∆l
the observed behavior of the critical fluidization/defluidization velocities. The general increase in the minimum fluidization velocity with an increase in the magnitude of cohesion is consistent with many experimental (Passos and Mujumdar, 2000; Seville and Clift, 1984, Wright and Raper, 1998a,b) and simulation (Kobayashi et al., 2003; Mikami et al., 1998) studies. These results agree with the results shown here. Unfortunately, very few efforts have provided enough detail to determine which critical velocity has been measured (theoretical or incipient). Therefore, it is difficult to comment on the validity of the results presented here for the theoretical minimum fluidization (Fig. 6b) and incipient fluidization (Fig. 6c).
0.0001 0
0
0.002
0.004 v
4.4. Quantitative comparison: average particle displacement To further evaluate the effectiveness of the escape-velocity mapping, particle displacement over a set time interval is examined, similar to that used by Wang and Rhodes (2004). An extreme example of small displacement is the plugging and channeling exhibited by Type C particles. On the other end of the spectrum is a vigorous, bubbling fluidization, in which net particle displacement is relatively high. The level of particle displacement is calculated by recording the distance moved by each particle (ti ) over the time increment t1 − t2 : li = (xt1 − xt2 )2 + (yt1 − yt2 )2 , (24) where xt1 and yt1 are the x- and y-positions of a given particle at time t1 , and xt2 and yt2 are the x- and y-positions of a given particle at time t2 . This quantity is then averaged over every particle in the simulation to obtain the average particle displacement (lavg ) in the given time increment: 1 li , lavg = n n
(25)
i=1
where n is the number of particles in the simulation. The values reported in Fig. 7 represent a cumulative average of 30 temporal data collection intervals that each lasts 0.3 s. These results were taken using the simulation parameters presented in Table 1 and the cohesion parameters presented in Table 3 with a constant superficial velocity of 0.77 m/s. Furthermore, the system was allowed to come to a steady state for 0.65 s before data was collected and three replicates were run for each case. Fig. 7 illustrates that the square-well model predicts average particle displacement similar to that of the Hamaker at similar escape velocities. 4.5. Effect of length scale of cohesive force One drawback of the escape-velocity mapping is that it does not account for the length scale of the square-well (router ) and Hamaker (Hcut ) models. For the square-well model, a change in the well-width-to-particle-radius ratio (router /rinner ) does not affect the escape velocity for head-on (normal) collisions (see Eq. (15)), but does change the escape velocity for oblique impacts. For the Hamaker model, changes in the cutoff distance
4521
ij,esc
0.006
(m /s)
Fig. 7. Average particle displacement with both cohesion models compared to escape velocity.
2
Approach angle
1
1
2
Fig. 8. Layout of two-particle simulations used to determine effect of square-well width on escape velocity.
(Hcut ) will affect the head-on escape velocity. Two approaches are taken to assess the effect of the characteristic length scale in each model. First, two-particle simulations are used to directly predict the escape velocity for various angles of approach (i.e., oblique collisions). Second, many-particle (4000) simulations are run using the square-well model with various values for router /rinner , and using the Hamaker model with various values for Hcut . The results of these simulations are then compared on the basis of predicted defluidization velocity and average particle displacement. Fig. 8 illustrates the configuration used to determine the effect of oblique collisions on the escape velocity. Fig. 9a reports the escape velocity determined via two-particle simulations using the square-well model with two different well depths and three different values of router /rinner . At larger router /rinner , the escape velocity varies somewhat, but overall the approach angle has little effect on the escape velocity except at very high approach angles. Similar results are obtained for the Hamaker model (with varied Hcut ) as displayed in Fig. 9b. As stated earlier, varying Hcut will affect the maximum Hamaker force and therefore affect the head-on escape velocity. To account for this variation, a different value of the Hamaker constant was used with each value of Hcut in order to achieve the same escape velocity for the head-on case (zero approach angle). The sets of parameters used are summarized in Table 4. While the
4522
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
0.01
router / rinner = 1.1
0.008
D=
1x10 -10 J
1.18
0.006 1.13 0.004 D = 1x10 -11 J 0.002 0
0
20
(a)
40 60 approach angle (degrees)
80
100
0.98 0.0001
Hcut = 4 x 10-11 m
0.04
Hcut = 4 x
1.08
1.03
0.05
vij,esc (m /s)
D=1x10-4 erg (vij,esc= 0.179 cm /s)
router / rinner = 1.001
U*mf, df
vij,esc (m /s)
D =1x10-3 erg (vij,esc= 0.566 cm/s)
router / rinner = 1.01
10-10 m
0.001
vij,esc= 0.00566 m/s
1.08
vij,esc= 0.00400 m/s vij,esc= 0.00179 m/s
0
(b)
20
40
60
80
100
approach angle (degrees)
U *mf, df
0
1
1.10 0.02 0.01
0.1
(router / rinner )-1
(a)
Hcut = 4 x 10-9 m
0.03
0.01
Fig. 9. Effect of approach angle on the escape velocity with (a) square-well model and (b) Hamaker model.
Table 4 Hamaker constants used to assess effect of cutoff distance at varying approach angle Hcut (m)
A (J)
vij ,esc (m/s)
4.0 × 10−9
1.04 × 10−16
0.0056 0.0056 0.0056 0.00179 0.00179 0.00179
4.0 × 10−10 4.0 × 10−11 4.0 × 10−9 4.0 × 10−10 4.0 × 10−11
1.06 × 10−18 1.06 × 10−20 3.15 × 10−19 3.34 × 10−19 3.36 × 10−21
escape velocity varies as the approach angle is increased (similar to the square-well model), the cutoff distance (Hcut ) has little very effect at each approach angle. For both cohesion models, the escape velocity deviates significantly from the “head-on” case only at very high collision angles. This deviation is partly due to the increased effects of friction as the particles develop a larger component of tangential relative velocity at higher approach angles. Two-particle simulations run without friction (results not shown) displayed a slight decrease in the escape velocity at very large approach angles. In addition to two-particle systems, the effects of length scale were also assessed on bulk parameters measured in many particle systems. The effects of router /rinner on the defluidization velocity are examined in Fig. 10a. At both well depths, the defluidization velocity changes very little as router /rinner is varied (if error bars are taken into account). Fig. 10b displays the effect of Hcut on the predicted defluidization velocity in simulations that utilize the Hamaker model.
1.06 1.04 1.02 1.00 1.E-11
(b)
1.E-10
1.E-09
1.E-08
1.E-07
Hcut (m)
Fig. 10. Effect of length scale on defluidization velocity: (a) effect of router /rinner for square-well model and (b) effect of Hcut for the Hamaker model.
In order to illustrate the effects of varying Hcut relative to varying the well depth (D), the results of simulations with higher and lower levels of cohesion are also presented. These results indicate the Hcut has little effect on the predicted defluidization velocity. 5. Summary Cohesive forces are known to be a continuous function of particle separation distance, and may result in enduring particle contacts (e.g., the formation of an agglomerate). Hence, the nature of cohesive forces conflicts with the assumption of binary, instantaneous interactions inherent in kinetic-theory (continuum) descriptions of particulate flows. In this work, a simplified model for cohesive forces, which has assumptions consistent with those of kinetic-theory treatments, has been assessed via comparison with more complex descriptions for cohesion. Specifically, a square-well model which implements cohesion as impulsive interactions between particles (thereby making it a viable method to incorporate cohesion into the kinetic theory) has been investigated. To gauge the effectiveness of the square-well model, cohesive forces have been
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
implemented into a discrete-particle fluidized-bed simulation using both the square-well model and the Hamaker model for van der Waals forces. Recall that the more complex Hamaker model is not suitable for incorporation into a continuum framework. Unlike the more elaborate Hamaker model, the square-well model is significantly limited by the lack of data characterizing real materials in terms of square-well parameters. Hence, a method to convert tabulated Hamaker constants into equivalent square-well parameters is required for the square-well model to be used as a predictive tool for particle flows. While there is no direct way to convert Hamaker constants to square-well depths, matching the escape velocity in a two-particle, headon collision was found to provide a good mapping between these two cohesion models. The effectiveness of this mapping technique was assessed by examining defluidization behavior, mixing behavior, fluidization velocities and average particle displacement. Wherever possible, simulation results from both models were also compared with experimental trends. Results of these simulations show that the square-well model makes predictions on par with those of the more elaborate Hamaker model. The incorporation of the square-well potential into a discrete-particle simulation is an important step in the possible implementation of cohesion into a kinetic-theory framework. Because the square-well model implements cohesion as impulsive interactions, it does not conflict with the assumption of binary, instantaneous interactions inherent in kinetic theory. Furthermore, the development of the escape-velocity mapping extends the usefulness of this model by providing a means to determine square-well parameters that correspond to material properties.
4523
where g is the void fraction, g is the gas density, vg is the gas velocity, S g is the gas-phase stress tensor, g is the acceleration due to gravity, and Igs is the rate of momentum transfer between the gas and solid phase per unit volume: vs − vg ), Igs = −s ∇Pg − Fgs (
(A.3)
where Pg is the gas-phase pressure, s is the solids volume fraction, and vs is the (locally averaged) solids velocity. The drag coefficient, Fgs , is determined using the expression developed by Symlal et al. (1994): Fgs =
3s g g 4Vt2 Dp
CDs | vs − vg |,
(A.4)
where Dp is the particle diameter and CDs is the single-sphere drag coefficient given by Dalla Valle (1948): 2 Vt CDs = 0.63 + 4.8 . (A.5) Re The terminal velocity, Vt , is given by the following correlation developed by Garside and Al-Dibouni (1977): Vt = 0.5 (A0 − 0.06Re 2 2 + (0.06Re) + 0.12Re(2B − A) + A ,
(A.6)
where A0 = 4.14 g
0.81.28 , 0.85, B = 2.65g g , g > 0.85
(A.7) (A.8)
and the Reynolds number, Re, is defined as Acknowledgments
Re =
The authors are grateful to the National Science Foundation for providing funding support for this work under Grant no. CTS 0226010. The authors would also like to thank Dr. Tom O’Brien and Dr. Dhanunjay S. Boyalakuntla for helpful discussions during the course of this work. Appendix A. Computational approach for (continuum) gas phase The multiphase flow with interphase exchanges (MFIX) program developed at the Department of Energy National Energy Technology Laboratory (DOE NETL) was used as a framework to model the gas–fluid system in this work (www.mfix.org). The fluid phase is modeled by solving locally averaged mass and momentum balances (Syamlal et al., 1994): j (g g ) + ∇ · (g g vg ) = 0, jt
(A.1)
j (g g vg ) + ∇ · (g g vg vg ) = ∇ · S g + g g g − Igs , (A.2) jt
Dp | vs − vg |g g
,
(A.9)
where g is the gas viscosity. For the purposes of calculating the void fraction, the twodimensional simulation domain is assumed to have a depth of one particle diameter (Hoomans et al., 1996; Kawaguchi et al., 1998; Tsuji et al., 1993; Tsuji, 2000, van Wachem et al., 2001; Xu and Yu, 1997). In each computational cell, the solids volume fraction is calculated as follows: nmn s,mn = i=1 , (A.10) Vmn where is the volume of the ith particle, s,mn is the solids volume fraction, Vmn is the volume of the grid cell and nmn is the total number of particles in the m, n grid cell. The average solids velocity is calculated by averaging over all the particles in a given grid cell: nmn vi vs,mn = i=1 , (A.11) nmn where vs,mn is the average solids velocity in grid cell m, n, and vi is the velocity of particle i.
4524
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
It should be noted that the method of calculating the void fraction (Xu and Yu, 1997) presented in Eq. (A.10) introduces some error since the simulation under consideration is twodimensional whereas drag correlations (which depend on the void fraction; see Eq. (A.4)) are based on three-dimensional systems. Equal numbers of particles will pack differently in two- and three-dimensional domains and therefore produce different void fractions. A discussion of the issues related to this mismatch, including possible methods to convert a twodimensional void fraction into an equivalent three-dimensional void fraction, is provided by van Wachem et al. (2001). In particular, they found that different strategies for converting the two-dimensional void fraction to a three-dimensional void fraction lead to differences in pressure fluctuations (bubble dynamics) and bed expansion. Such methods are not employed in this work, since the various methods are found to have little effect on the predicted minimum fluidization velocities (as was the focus in this work), which match closely with those obtained from empirical relations (Weber, 2004). The system of equations governing the gas-phase are solved numerically using a modified version of the SIMPLE method developed by Patankar (1980). For more details on the solution of the gas-phase equations, the reader is referred to the MFIX Numerical Theory Guide (Syamlal, 1998). Appendix B. Computational approach for (discrete) particle phase The three forces acting on a particle which are not described in the main text include the fluid-induced forces (Ff,i ), gravitational forces (Fg,i ), and forces arising from (repulsive) collisional contacts with other particles and/or the walls (Fcoll,i ). Fluid-particle forces (Ff,i ) are calculated using Eq. (A.4). The forces induced by the fluid are equal and opposite to the forces exerted by the particles on the fluid. For each computational cell, the total force imposed by the particles on the fluid is divided by the number of particles to determine the force imposed on each particle by the fluid (Boyalakuntla, 2003). Gravitational forces (Fg,i ) are calculated based on the volume of a spherical particle and the acceleration due to gravity: 4 3 Fg,i = rinner s g, 3
(B.1)
where g is the acceleration due to gravity. The repulsive, collisional force acting on particle i is composed of a normal and tangential component Fcoll,i = Fn,ij kij + Ft,ij kt,ij , (B.2) where Fn,ij and Ft,ij are the magnitudes of the normal and tangential components of the repulsive force when the ith and j th particles collide. The directions of these forces are along the line of the normal unit vector pointing from the center of the ith particle to the center of the j th particle (kij ) and the tangential unit vector (kt,ij ), which is perpendicular to the normal vector with the cross product (kij × kn,ij ) pointing out of the page.
The summations represent all collisions experienced by the ith particle during the given time step. The magnitude of the forces is calculated using a linear-spring-and-dashpot model (Tsuji et al., 1993) characterized by normal (Kn ) and tangential (Kt ) spring constants, as well as dissipation constants (n , t .). The basic equation for calculating the normal force between two particles is given by Fn,ij = − n,ij Kn − n
j n,ij , jt
(B.3)
where n,ij is the normal component of the overlap between particles i and j . The normal overlap can be calculated as the difference between the sum of the particle radii (2rinner ) and the distance between particle centers: n,ij = 2rinner − (xi − xj )2 + (yi − yj )2 , (B.4) where xi and yi are the x- and y- components of the position of the ith particle, and xj and yj are the x- and y- components of the position of the j th particle. The tangential component of the contact force is calculated as follows: Ft,ij = − t,ij Kt − t vt,ij ,
(B.5)
where Kt is the tangential spring constant and t is the tangential damping coefficient. In this expression, vt,ij is the tangential slip velocity between the two particles: vt,ij = vij · kt,ij − rinner (i + j ),
(B.6)
where i and j are the rotation rates of the ith and j th particles (positive is counter-clockwise), respectively, and rinner is the particle radius. Also, t,ij is the tangential overlap, defined as follows: t,ij = ts · vt,ij .
(B.7)
Tangential interactions are calculated as either sticking or sliding contacts based on the following inequality: |Ft,ij | |Fn,ij |,
(B.8)
where is the coefficient of friction. If this inequality is not satisfied, then the particles undergo a sliding contact and the tangential component of the collision force is calculated as Ft,ij = Fn,ij .
(B.9)
To reduce computational times, the spring constant is made somewhat softer than realistic materials (as is often done; see, for example, Tsuji et al., 1993; Gera et al., 1998; Kawaguchi et al., 1998). The damping coefficient is still set, though, to achieve a realistic coefficient of restitution (e = 0.94). Several researchers (Gera et al., 1998; Kawaguchi et al., 1998; Kuwagi et al., 2000) have indicated that the use of an artificially small spring constant does not have a significant effect on the nature of the flow in fluidized-bed simulations. Nonetheless, the results of the softened particles were spot-checked against simulations using spring constants that led to more realistic collision times to ensure that physical accuracy was not compromised. The spring constant does have a slight effect on the calculated
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
minimum fluidization velocity, but the pressure drop at fluidization is not affected. The observed effect on the minimum fluidization velocity is similar for both cohesion models (squarewell and Hamaker), and thus does not impact the comparison between the two models. Since an initial goal of the development phase was to reproduce the results of Rhodes et al. (2001a), the default parameters in Table 1 closely mimic the conditions in that effort.
Fcoll,i Ff,i Fg,i Fgs Fn,ij Fsw,ij Ft,ij FvdW ,ij g H1 H2 Hcut Hij Hiw I Igs Jcoh,∗ Jcoh,app,ij
Jcoh,esc,ij
kij kt,ij kij ,1 kij ,2 kij ,3
Notation ai A Ab Ajet D Dp e Fad Fi Fcoh,i
Jcoh,cap,ij
acceleration of ith particle, m/s2 Hamaker constant, J fluidized-bed cross-sectional area, m2 inlet jet cross-sectional area, m2 depth of square well, J particle diameter, m coefficient of restitution adhesion force, N sum of forces on ith particle, N cohesive interpaticle forces imposed on ith particle, N sum of forces from collisional contacts on ith particle, N forces induced by fluid on ith particle, N gravitational forces on ith particle, N drag coefficient, kg/m3 /s normal component of collisional force between ith and j th particle, N magnitude of force imposed on ith and j th particles from square-well interaction, N tangent component of collisional force between ith and j th particle, N cohesive van der Waals forces between the ith and j th particles, N accelerationdue to gravity, m/s2 average height of particles in family one, m average height of particles in family two, m cutoff distance for Hamaker van der Waals model, m surface-to-surface separation distance between ith and j th particles, m surface-to-surface distance between ith particle and wall, m mixing index rate of momentum transfer between gas and solid phase, N/m3 momentum impulse transferred between particles during square-well interaction, kg m/s momentum impulse transferred during approaching-cohesive interaction between ith and j th particles, kg m/s momentum impulse transferred during escaping-cohesive interaction between ith and j th particles, kg m/s
Kn Kt lavg li m mmax n nmax nmn nvdW P P ∗
Pg rinner router, Sg t ts tf Us Us∗ Umf ,df ∗ Umf ,dr Umf ,in ∗ Umf ,in Umf ,th ∗ Umf ,in vg vi vij ,a
4525
momentum impulse transferred during capturecohesive interaction between ith and j th particles, kg m/s normal unit vector pointing from center of the ith particle to center of the j th particle tangent unit vector normal unit vector at instance of approaching cohesive interaction normal unit vector at instance of physical contact normal unit vector at instance of departing cohesiveinteraction normal spring constant for particle–particle contacts (stiffness), N/m tangential spring constant for particle–particle contacts, N/m average particle movement, m distance moved by ith particle during current time step, m particle mass, kg number of computational cells in x-direction total number of particles in simulation number of computational cells in y-direction number of particles in grid cell m, n number of van der Waals contacts experienced by a particle at a given time step pressure drop across fluidized bed, Pa pressure drop across fluidized-bed normalized by weight of particles above lowest computational node, Pa gas-phase pressure, Pa particle radius, m width of square well, m gas-phase stress tensor, N/m2 time, s solid-phase time step, s gas-phase time step, s superficial velocity in fluidized bed, m/s normalized superficial gas velocity (Us /(Umf ,th,non-cohesive )) defluidization velocity, m/s defluidization velocity normalized by defluidization velocity for non-cohesive case incipient fluidization velocity,m/s incipient fluidization velocity normalized by incipient fluidization velocity for non-cohesive case theoretical minimum fluidization velocity, m/s theoretical minimum fluidization velocity normalized by theoretical minimum fluidization velocity for non-cohesive case gas velocity, m/s velocity of ith particle, m/s volume of ith particle, m3 relative velocity before approaching cohesive interaction, m/s
4526
vij ,c vij ,esc vs vt,ij vi Vjet Vt Vmn W W xi yi
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527
relative velocity after inelastic collision, m/s minimum normal, relative velocity needed to escape cohesive forces, m/s locally averaged solid velocity, m/s slip velocity between ith and j th particles, m/s velocity change of the ith particle in a given time step, m/s specified as velocity at inlet jet, m/s terminal velocity, m/s volume of grid cell m, n, m3 total weight of particles in fluidized bed, kg total weight of particles above lowestcomputational node, kg x-position of ith particle, m y-position of ith particle, m
Greek letters n,ij t,ij g n t g g s i
surface energy, J normal component of overlap between ith and j th particles, m tangential component of overlap between ith and j th particles, m void fraction normal damping coefficient for particle–particle contacts, N s/m tangential damping coefficient for particle–particle contacts, N s/m gas viscosity, kgs1 m1 friction coefficient gas density, kg/m3 particle density, kg/m3 rotation rate of ith particle, 1/s
References Allen, M.P., Tildesley, D.J., 1987. Computer Simulation of Liquids. Oxford University Press, Oxford. Boyalakuntla, D., 2003. Simulation of granular and gas–solid flows using discrete element method. Ph.D. Thesis, Carnegie Mellon University, www.mfix.org. Dalla Valle, J.M., 1948. Micrometrics. Pitman, London. Forsyth, A.J., Hutton, S., Rhodes, M.J., 2002. Effect of cohesive interparticle force on the flow characteristics of granular material. Powder Technology 126, 150–154. Garside, J., Al-Dibouni, M.R., 1977. Velocity-voidage relationships for fluidization and sedimentation in solid–liquid systems. Industrial and Engineering Chemistry Process Design and Development 16, 206–214. Geankopolis, C.J., 1983. Transport Processes and Unit Operations. Allyn and Bacon, Boston. Geldart, D., 1973. Types of gas fluidization. Powder Technology 7, 285–292. Gera, D., Gautam, M., Tsuji, Y., Kawaguchi, T., Tanaka, Y., 1998. Computer simulation of bubbles in large-particle fluidized beds. Powder Technology 98, 38–47. Grace, J.R., Taghipour, F., 2004. Verification and validation of CFD models and dynamic similarity for fluidized beds. Powder Technology 139, 99–110. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidised bed: a hard-sphere approach. Chemical Engineering Science 51, 99–118.
Israelachvili, J., 1991. Intermolecular and Surface Forces. Academic Press London, London. Kawaguchi, T., Tanaka, T., Tsuji, T., 1998. Numerical simulation of twodimensional fluidized beds using the discrete element method (comparison between two- and three-dimensional models). Powder Technology 96, 129–138. Kim, H., Arastoopour, H., 2002. Extension of kinetic theory to cohesive particle flow. Powder Technology 122, 83–94. Kobayashi, T., Mukai, T., Kawaguchi, T., Tanaka, T., Tsuji, Y., 2003. DEM analysis of flow patterns of Geldart’s group A particles in a fluidized bed. World Congress on Particle Technology, vol. 4. Sydney, Australia. Kuwagi, K., Mikami, T., Horio, M., 2000. Numerical simulation of metallic solid bridging particles in a fluidized bed at high temperature. Powder Technology 109, 27–40. McKeen, T., Pugsley, T., 2003. Simulation and experimental validation of a freely bubbling bed of FCC catalyst. Powder Technology 129, 139–152. McLaughlin, L.J., Rhodes, M.J., 2001. Prediction of fluidized bed behaviour in the presence of liquid bridges. Powder Technology 114, 213–223. Mikami, T., Kamiya, H., Horio, M., 1998. Numerical simulation of cohesive powder behavior in a fluidized bed. Chemical Engineering Science 53, 1927–1940. Molerus, O., 1982. Interpretation of Geldart’s type A, B, C, and D powders by taking into account interparticle cohesion forces. Powder Technology 33, 81–87. Mutsers, S.M.P., Rietema, K., 1977. The effect of interparticle forces on the expansion of a homogeneous gas-fluidized bed. Powder Technology 18, 239–248. Nase, S.T., Vargas, W.L., Abatan, A.A., McCarthy, J.J., 2001. Discrete characterization tools for cohesive granular material. Powder Technology 116, 214–223. Passos, M.L., Mujumdar, A.S., 2000. Effect of cohesive forces on fluidized and spouted beds of wet particles. Powder Technology 110, 222–238. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere, New York. Rhodes, M., 1998. Introduction to Particle Technology. Wiley, West Sussex, England. Rhodes, M.J., 2001. Fluidization of particles by fluids. www.erpt.org/ 012Q/rhod-04.htm. Rhodes, M.J., Wang, X.S., Nguyen, M., Stewart, P., Liffman, K., 2001a. Use of discrete element method simulation in studying fluidization characteristics: influence of interpaticle force. Chemical Engineering Science 56, 69–76. Rhodes, M.J., Wang, X.S., Nguyen, M., Stewart, P., Liffman, K., 2001b. Onset of cohesive behavior in gas fluidized beds: a numerical study using DEM simulation. Chemical Engineering Science 56, 4433–4438. Seville, J.P.K., 1997. Processing of Particulate Solids. Blackie Academic and Professional, London. Seville, J.P.K., Clift, R., 1984. The effect of thin liquid layers on fluidization characteristics. Powder Technology 37, 117–129. Syamlal, M., 1998. MFIX Documentation Numerical Technique. US Department of Energy, Morgantown, West Virginia, www.mfix.org. Syamlal, M., O’Brien, T.J., 2003. Fluid dynamic simulation of O3 decomposition in a bubbling fluidized bed. A.I.Ch.E. Journal 49, 2793–2801. Syamlal, M., Rogers, W., O’Brien, T.J., 1994. MFIX Documentation Theory Guide. US Department of Energy, Morgantown, West Virginia, www.mfix.org. Tatemoto, Y., Mawatari, Y., Noda, K., 2005. Numerical simulation of cohesive particle motion in vibrated fluidized bed. Chemical Engineering Science 60, 5010–5021. Tsinontides, S.C., Jackson, R., 1993. The mechanics of gas fluidized beds with an interval of stable fluidization. Journal of Fluid Mechanics 255, 237–274. Tsuji, Y., 2000. Activities in discrete particle simulation in Japan. Powder Technology 113, 278–286. Tsuji, Y., Kawaguchi, T., Tanaka, T., 1993. Discrete particle simulation of two-dimensional fluidized bed. Powder Technology 77, 79–87. van Wachem, B.G.M., van der Schaaf, J., Schouten, J.C., Krishna, R., van den Bleek, C.M., 2001. Experimental validation of Lagrangian–Eulerian simulations of fluidized beds. Powder Technology 116, 155–165.
M.W. Weber, C.M. Hrenya / Chemical Engineering Science 61 (2006) 4511 – 4527 Wang, X.S., Rhodes, M.J., 2004. Mechanistic study of defluidization by numerical simulation. Chemical Engineering Science 59, 215–222. Weber, M.W., 2004. Simulation of cohesive particle flows in granular and gas–solid systems. Ph.D. Thesis, University of Colorado, www.mfix.org. Weber, M.W., Hoffman, D.K., Hrenya, C.M., 2004. Discrete-particle simulations of cohesive granular flow using a square-well potential. Granular Matter 6, 239–254. Wright, P.C., Raper, J.A., 1998a. Role of liquid bridge forces in cohesive fluidization. Transactions of the Institutions of Chemical Engineers 76, 753–760. Wright, P.C., Raper, J.A., 1998b. Examination of dispersed liquid-phase threephase fluidized beds: part 1 Non-porous, uniform particle systems. Powder Technology 97, 208–226.
4527
Wu, S.Y., Baeyens, J., 1998. Segregation by size difference in gas fluidized beds. Powder Technology 98, 139–150. Xie, N., Battaglia, F., Fox, R.O., 2004. Simulations of multiphase reactive flows in fluidized beds using in situ adaptive tabulation. Combustion Theory and Modeling 8, 195–209. Xu, B.H., Yu, A.B., 1997. Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science 52, 2785–2809. Xu, B.H., Zhou, Y.C., Yu, A.B., Zulli, P., 2002. Force structures in gas fluidized beds of fine powders. World congress on particle technology, vol. 4. Sydney, Australia. Ye, M., van der Hoef, M.A, Kuipers, J.A.M., 2004. A numerical study of fluidization behavior of Geldart A particles using a discrete particle model. Powder Technology 139, 129–139.