Behavior of fluidized beds similar to equilibrium states

Kengo Ichiki^{∗} & Hisao Hayakawa^{†}

Graduate School of Human and Environmental Studies, Kyoto University, Kyoto 606-01, Japan (April 25, 1997)

Abstract

Systematic simulations are carried out based on the model of fluidized beds proposed by the present authors [K.Ichiki and H.Hayakawa, Phys. Rev. E 52, 658 (1995)]. From our simulation, we confirm that fluidization is a continuous transition. We also confirm the existence of two types of fluidized phases, the channeling phase and the bubbling phase. We find the close relations between the averaged behaviors in fluidized beds and quasi equilibrium states in dense liquids. In fluidized beds, (i) the flow rate plays the role of the effective temperature, and (ii) the existence of a kind of the fluctuation-dissipation relation is suggested. 05.60.+w, 03.20.+i, 82.20.Mj, 83.70.Hq

Typeset using REVTEX

∗e-mail: ichiki@phys.h.kyoto-u.ac.jp †e-mail: hisao@phys.h.kyoto-u.ac.jp

I. INTRODUCTION

Recently granular materials have been studied extensively from both experimental and theoretical point of views in the context of the nonequilibrium statistical physics [1–3]. Since the granular materials are dissipative, energy injections are necessary to preserve steady states. Many of recent studies for granular materials are focused on the behavior of systems excited by mechanical activations such as vibration or rotation of vessels. On the other hand, the researches on fluidized beds, where systems are excited by the fluid flow, are not relatively advanced in spite of their variety of dynamical behaviors [4,5].

Fluidized beds are widely used in chemical industries since early 19’s, and they have been studied from technological point of views. Fluidized beds consist of granular particles confined in a tall chamber with distributor for the fluid flow at the bottom. In experiments, energy injection to the system is controlled by the flow rate of fluid. At low flow rate, the system is in the fixed phase where particles rest on the bottom. When the flow rate exceeds the critical value, particles start moving. This state is called the fluidized phase, which contains sub phases, for instance, the homogeneous phase, the bubbling phase, the channeling phase, etc.

There are many models to describe fluidized beds, which can be classified in two categories, two-fluid models and particle-dynamics models. In the two-fluid models, particles are treated as a fluid [6–10]. These models have benefit on analytical treatments and to generalize their discussion to other systems [11]. However, their bases, such as constitution equations for the particle-phase pressure and the stress tensor, have not been established. On the other hand, the particle-dynamics models describe the direct motion of particles. There are various models, which are kinetic theories [5], the discrete element method [12–14], etc. However, these models cannot be the basis for the two-fluid models. The main problem of them is that hydrodynamic interactions among particles are over simplified. For instance, the boundary condition between particles and fluid is not satisfied in the scale of particles, and the fluid equation is calculated under inviscid limit.

Recently the present authors have proposed a numerical model of fluidized beds, where hydrodynamic interaction among particles is calculated with reliable accuracy [15,16]. In this paper, we will present the results of our systematic simulations and behaviors on statistical quantities obtained from the simulations.

The contents of this paper are as follows. In Sec. II, we review the method of our simulation. We show the results in Sec. III, where we observe the transition of fluidization and the existence of two fluidized phases. We also discuss statistical quantities, which are analogous to equilibrium correspondences. In Sec. IV we give an interpretation of the averaged quantities by the hole theory for simple liquids [17,18]. In Sec. V, we conclude our results. In Appendices, we summarize the method of our simulation and discuss theoretical difficulties in the modeling of fluidized beds.

II. SIMULATION METHOD

In this section, we briefly explain our model and how to simulate the dynamics of granular particles in fluid flows. The detail explanation of our model can be seen in Ref. [15,16] and Appendices. For simplicity, we only consider the cases of monodispersed spherical particles. We assume the following equation of motion for the particles

St\frac{dU}{dt} = −U + V + F_{c},  (2.1)

where St is the effective Stokes number, F_{c} represents hard-core collision among particles, U is the velocity of particles. V is the terminal velocity determined by

V − u^{∞} = −R^{↔}^{−1} · E_{z},  (2.2)

where u^{∞} is the flow rate of induced fluid which is equal to the superficial velocity conventionally used for the fluidized beds. R^{↔} is the resistance matrix representing the hydrodynamic interaction among particles calculated by the method of the Stokesian dynamics [19], where periodic boundary condition is adopted as the effect of chamber (see Appendix A 2).

Hard-core collisions are assumed to be elastic and calculated by the momentum exchange for contacting particles in simulations (see Appendix A 1). The bold-face letters without superscripts represents vectors in 3N-dimension, where N is the number of particles in the unit cell of periodic boundary condition. For example, the velocity U has the following components

U = [ U^{(1)} ; U^{(2)} ; ⋮ ; U^{(N)} ] ,  (2.3)

where the bold-face letters with superscripts represent vectors in 3-dimension. In this paper, we use dimensionless quantities with the aid of the particle radius a and the sedimentation velocity of a single particle in a viscous fluid U_{0} = mg̃/6πμa, where μ is the viscosity of the fluid, m is the mass of the particles and g̃ = g(ρ_{p} − ρ_{f})/ρ_{p} with the gravitational acceleration g, and the densities of the particle ρ_{p} and the fluid ρ_{f}. Equation (2.1) represents the relaxation process of U to V with the time-scale St.

In Eq.(2.1), there are two control parameters, the effective Stokes number St and the flow rate u^{∞}. For the parameters related to the system size, we choose the number of mobile particles N_{M} = 256, that of fixed particles N_{F} = 10 and the size of the unit cell in periodic boundary condition (L_{x}, L_{y}, L_{z}) = (34, 2, 100). In this situation, particles are confined in the vertical plane, while hydrodynamic interactions are considered in 3-dimensional space. We adopt the fixed phase as initial conditions of our simulations, which is constructed from simulations with u^{∞} = 0. The choice of the system size and these artifact situations come from the limitation of computer resources. We have checked that statistical quantities seems to be insensitive to the choice of L_{z} within the range of 50 ≤ L_{z} ≤ 100 and the choice of the initial conditions is not relevant from the comparison of results with other initial conditions. We have also confirmed that qualitatively similar behaviors to this situation are observed in 3-dimensional simulations and in the case of N_{M} = 133.

III. SIMULATION RESULTS

A. General behavior

In this section, we present the results of our simulations in details. We perform simulations at the points in the parameter space (Fig. 1) within the range of 0.05 ≤ u^{∞} ≤ 0.8 and 0.1 ≤ St ≤ 100, where we observe fixed, bubbling and channeling phases. Transitions among these phases will be discussed below.

In the fixed phase at low flow rate, particles are rest at the bottom. At the critical flow rate u^{∞} = u_{c}, the particles begin to be fluidized. It seems that the transition between the fixed phase and the fluidized phase is independent of St. We observe two fluidized phases. One is the channeling phase observed for small St where we can see a channel or a path of fluid flow. Another is the bubbling phase observed for large St where bubbles raise through the particle beds. Typical snapshots of them are shown in Figs. 2 and 3. We show the area of channel-bubble transition observed in our simulations as the transition area in Fig. 1, where we will see that statistical quantities qualitatively change their behaviors.

To characterize the transition of fluidization quantitatively, we calculate the kinetic energy per particle E(t) defined by

E(t) = \frac{1}{N_{M}} ∑_{α=1}^{N_{M}} |U^{α}(t)|^{2}.  (3.1)

A typical behavior of E(t) in the bubbling phase is shown in Fig. 4. We observe regular behavior after the minimum point following the first peak of E(t). In the channeling phase, we also observe qualitatively similar behavior of E(t) to Fig. 4, though the period of peaks is smaller and each peak is not distinguishable. We now introduce the average of E(t) defined by

Ē = \frac{1}{ΔT} ∫_{ΔT} dt E(t),  (3.2)

where ΔT is the period of regular behavior in E(t). From Fig. 5, we observe a continuous transition at u^{∞} = u_{c} and the linear behavior of Ē(u^{∞}) in the fluidized phase (u^{∞} > u_{c}). This behavior can be well fitted by

Ē(u^{∞}) = \{ \begin{matrix}0 & (u^{∞} < u_{c}) \\ A_{E}(u^{∞} − u_{c}) & (u^{∞} > u_{c}) \end{matrix} ,  (3.3)

where A_{E} and u_{c} are the fitting parameters which depend on St. Equation (3.3) defines u_{c} shown in Fig. 1.

Now we show that the transition of fluidization can be understood as the process generating the free volume around the particles in the fixed phase. It is useful to remember that our model is Galilei invariant, that is, the system with the fixed particles of U_{F} = 0 under the flow rate u^{∞} is equivalent to that with the fixed particles of U_{F} = −u^{∞}E_{z} under the flow rate 0. Let us consider the process under the latter situation. First we define U_{fall} which is the falling velocity of the mobile particles in the fixed phase without the support of fixed particles. If the flow rate u^{∞} is smaller than U_{fall}, the mobile particles cannot pass over the fixed particles moving downward with −u^{∞}. Therefore the mobile particles hold on the fixed particles and the gap between the mobile particles and the fixed particles is not generated. While the flow rate u^{∞} is larger than U_{fall}, the mobile particles apart from the fixed particles and then the gap between them is generated. The gap causes the Rayleigh-Taylor instability observed when a heavy fluid exists above a light fluid [20]. Then the gap may grow into a bubble and propagate upward through the particles, or may construct a channel. From this discussion, the critical flow rate u_{c} is determined by the falling velocity U_{fall}. This suggests that u_{c} is independent of St because the falling velocity U_{fall} can be evaluated as the sedimentation rate of suspensions [21].

Next we discuss the channel-bubble transition. In view of Figs. 2 and 3, it is hard to distinguish the channeling phase from the bubbling phase. At first, we show the variance V_{H} defined by

V_{H} = \frac{1}{ΔT} ∫_{ΔT} dt (((H(t) − H̄)))^{2},  (3.4)

where H and H̄ are the height of the center of mass and its average defined by

H(t) = \frac{1}{N_{M}} ∑_{α=1}^{N_{M}} z^{α}(t)  (3.5)

and

H̄ = \frac{1}{ΔT} ∫_{ΔT} dt H(t).  (3.6)

We expect that V_{H} is small in the channeling phase and large in the bubbling phase. Figure 6 shows the corresponding behavior for u^{∞} = 0.3 and the transition is observed around St = 5. For other cases, the channel-bubble transitions are observed in the area shown in Fig. 1.

We also observe some qualitative changes in the physical quantities corresponding to the channel-bubble transitions. Here we discuss the behavior of Ē(St). From Fig. 7, we see that Ē increases with St in the channeling phase, and Ē decreases with St in the bubbling phase. The behavior in the channeling phase can be understood by the following reason. At St = 0, we observe a steady channel, no relative motions among particles. When St increases, the particles on the channel can be fluidized. This is because St is the relaxation time of the particle velocity U to their terminal velocity V, which is determined under the case of St = 0, and the non-relaxed particles cause the collapse of the channel. Therefore, Ē(St) increases with St in the channeling phase. On the other hand, we observe that Ē decreases with St in the bubbling phase. This can be understood as follows. From the previous paper [15] or in Fig. 4, we have seen that the relative motion of particles or the kinetic energy is generated by bubbles which have the difference of the local volume fraction in the system. In fact, bubbles are always accompanied by the convective motion of particles. Because St is the response time of the velocity for the change of configurations, bubbles become obscure as St increases. In fact, we see a sharp bubble and the definite convective motion of particles for St = 10 (Fig. 8), while we see a relatively obscure bubble and the weaker convection for St = 100 (Fig. 9). Therefore Ē(St) in the bubbling phase decreases with St. As a result, Ē(St) has a peak around the transition point. This reflects on the St dependence of various properties such as μ_{e}(St) shown in Sec. III B and H̄(St) discussed in Sec. IV.

B. The analogy to equilibrium systems

In this section, we demonstrate the existence of surprising correspondences in statistical quantities between fluidized beds which is in highly nonequilibrium states and quasi equilibrium systems. The result of Ē(u^{∞}) in Fig. 5 suggests that the flow rate u^{∞} behaves as the effective temperature of the environment such as that of the heat bath for equilibrium systems. Therefore, the critical flow rate u_{c} may correspond to the critical temperature and Ē may be the order parameter of the fluidization. In the following, we will interpret the results of our simulations using this effective temperature u^{∞}.

In nearly equilibrium systems at temperature T, the resistance ζ of a tracer particle is given by

ζ = \frac{kT}{D} ,  (3.7)

where D is the diffusion constant and k is the Boltzmann constant [22]. This is the Einstein relation and is the simplest form of the fluctuation-dissipation theorem which relates the correlation functions in the equilibrium state and the transport coefficients.

Similarly, we introduce for our simulations the effective viscosity μ_{e} as

μ_{e} = \frac{u^{∞}}{D_{p}} ,  (3.8)

where D_{p} is the diffusion constant of our simulations defined by

D_{p} = \frac{1}{2} \frac{1}{N_{M}} ∑_{α=1}^{N_{M}} |Ũ^{α}(ω = 0)|^{2},  (3.9)

where Ũ^{α}(ω) is the Fourier transform of U^{α}(t) calculated by the standard FFT. Equation (3.8) can be regarded as an extension of Eq. (3.7) with the aid of the effective temperature u^{∞}, because the viscosity is usually proportional to the resistance. The observed viscosity μ_{e} defined by (3.7) is shown in Figs. 10 and 11.

If the definition of the viscosity (3.8) is self-consistent, the Einstein relation or the fluctuation-dissipation relation in fluidized beds may be valid with the replacement of kT by u^{∞}. This statement is interesting because the system is in a highly nonequilibrium state and there is no reason of the existence of the fluctuation-dissipation relation in the sense of linear nonequilibrium statistical mechanics.

We find that the flow-rate dependence of the viscosity μ_{e}(u^{∞}) obtained by Eq. (3.8) obeys the Arrhenius function

μ_{e}(u^{∞}) ∝ e^{ε/u^{∞}},  (3.10)

where ε is a fitting parameter. The fitting by Eq. (3.10) with ε = 0.113±0.017 is also shown in Fig. 10. We compare our result of μ_{e} with experimental result in fluidized beds [23]. In the experiment of fluidized beds, the shear viscosity measured by the modified Stormer viscometer also obeys the Arrhenius function of Eq.(3.10) [23]. Therefore, our result from Eq. (3.8) is consistent with the experiment. This behavior which can be understood by a dense liquid theory in part will be explained in Sec. IV.

The connection between the non-Gaussian property and the dissipation in the system have been discussed for granular materials [24,25,15]. Here we check the non-Gaussian property in the velocity distribution functions P (U_{x}) obtained from our simulations. To characterize the non-Gaussian property of the velocity distribution, we calculate the 4th cumulant C_{4} defined by

C_{4}(U_{x}) = ⟨U_{x}^{4}⟩ − 3⟨U_{x}^{2}⟩^{2} − 4⟨U_{x}⟩⟨U_{x}^{3}⟩
+12⟨U_{x}⟩^{2}⟨U_{x}^{2}⟩ − 6⟨U_{x}⟩^{4}.  (3.11)

The resultant behaviors of C_{4} are shown in Figs. 12 and 13, where they are scaled by the square of variance (or 2nd cumulant) defined by

C_{2}(U_{x}) = ⟨U_{x}^{2}⟩ − ⟨U_{x}⟩^{2}.  (3.12)

This non-Gaussian parameter C_{4}/(C_{2})^{2} is zero for the Gaussian distribution and 3 for the exponential distribution. These behaviors of C_{4}/(C_{2})^{2} in Figs. 12 and 13 are similar to those of the effective viscosity μ_{e} in Figs. 10 and 11. In fact, we can also fit C_{4}/(C_{2})^{2} by the Arrhenius function

C_{4}/(C_{2})^{2} ∝ exp(ε′/u^{∞}),  (3.13)

as shown in Fig. 12 with ε′ = 0.175 ± 0.045. It will be an interesting subject that we will study a quantitative relation between (3.10) and (3.13).

IV. DISCUSSIONS

In this section, we demonstrate that the qualitative understanding of the results of our simulation is possible with the aid of the hole theory applied to simple liquids.

First we discuss the flow rate dependence of the viscosity μ_{e}(u^{∞}). The hole theory, which is used for the behavior of simple liquids, is based on the following picture. A molecule in a liquid can move when a free volume or a hole is generated around it. From this picture, the empirical relation of the viscosity μ_{l} [17] is derived as

μ_{l}(T) ∝ exp (\frac{ε_{l}}{kT} ) ,  (4.1)

where ε_{l} is the activation energy to create the free volume. If we regard u^{∞} as the temperature, Eq. (3.10) is identical to Eq. (4.1). This suggests that the averaged behavior in the fluidized beds may be understood by the hole theory in simple liquids [26–28].

Next we discuss St dependence of the height of center of mass H̄(St) defined by (3.6). In Fig. 14 we show a typical behavior of H̄(St) at u^{∞} = 0.3. We see the qualitative change of behavior in the transition area in Fig. 1, where H̄ is almost constant in the channeling phase, and H̄ increases logarithmically in the bubbling phase. Since the change of behavior with St in the channeling phase is only how the channel collapses, H̄(St) is expected to be independent of St. While the behavior in the bubbling phase is interesting. We can fit the data in Fig. 14 as

H̄(St) = C_{H} log(St) + D_{H},  (4.2)

where the fitting parameters C_{H} and D_{H} depend on u^{∞}. Equation (4.2) can be rewritten as

St = exp (\frac{H̄ − D_{H}}{C_{H}} ) .  (4.3)

Here H̄ − D_{H} can be understood as the volume expansion ΔV, since D_{H} is the height at St = 1. We show u^{∞} dependence of C_{H} in Fig. 14. From this figure, C_{H} is approximately proportional to the effective temperature u^{∞}. Since St is the characteristic time τ of the system, we rewrite Eq.(4.3) as

τ ∝ exp (F \frac{ΔV}{u^{∞}} ) ,  (4.4)

where F is a constant. Equation (4.4) is consistent with Eq.(4.1) under the reasonable assumption where the viscosity is characterized by the time τ to generate the free volume, and the activation energy to generate the free volume or expansion is proportional to ΔV.

Before closing this section, we give some remarks. Although the transition of fluidization in experiments seems to be the discontinuous phase transition [29], our simulations suggests a continuous phase transition. Also it is an open problem that at present we cannot reproduce homogeneous phase in our simulation. For these problems, we need to examine carefully the difference between the experiments and the simulations, and we must investigate the behavior near the critical flow rate u_{c} in detail, because the discontinuous phase transition and the homogeneous phase are observed there in experiments.

V. CONCLUSIONS

In this paper, we have carried out systematic simulations with the change of two control parameters, the flow rate of the fluid u^{∞} and the effective Stokes number St. When the flow rate u^{∞} is small, particles rest in the fixed phase. Above the critical flow rate u_{c}, particles are fluidized. The critical value u_{c} is independent of St. We have found two fluidized phases, the channeling phase and the bubbling phase, where the former changes to the latter as St increases.

We have found that the flow rate u^{∞} plays the role of the effective temperature. In terms of the effective temperature u^{∞}, we have defined the effective viscosity μ_{e} with the aid of the Einstein relation. The flow-rate dependence of the viscosity μ_{e} is similar to that in the experiments in real fluidized beds. We also find that the viscosity μ_{e}(u^{∞}, St) can be an index of the non-Gaussian property in velocity distribution of particles. This property is consistent with the behavior on granular materials or the system of inelastic particles. Qualitative behavior of fluidized beds such as μ_{e}(u^{∞}) and H̄(St) can be understood by means of the hole theory which has been used for simple liquids.

ACKNOWLEDGMENTS

We would like to thank S. Sasa, Y-h. Taguchi, T. Ooshida for fruitful discussions on the subject of this work. K.I. also thanks M. Doi and J. F. Brady for encouragements and useful comments on this work. This work has been partially supported by JSPS Research Fellowships for Young Scientists and Grant in Aid of Ministry of Education of Science and Culture in Japan. Numerical simulations was carried out by the facilities of the Supercomputer Center at the Institute for Solid State of Physics at the University of Tokyo.

APPENDIX A: THE IDEAL MODEL OF FLUIDIZED BEDS

Here we review our model of fluidized beds, presented in the previous paper [15], where only essences are extracted from real systems and all other irrelevant mechanisms are neglected.

1. Equation of motion

We construct the model by only four mechanisms, which are the inertial effect of the particles, hydrodynamic interaction through the fluid, the gravitational force and the contact force. Therefore the equation of motion can be written as

St_{0} \frac{dU}{dt} = F_{f} + F_{g} + F_{c},  (A1)

where F_{f}, F_{g}, F_{c} are the force from the fluid, the gravitational force and the contact force respectively. St_{0} is the bare Stokes number defined by

St_{0} = \frac{mU_{0}}{6πμa^{2}} .  (A2)

Particles are assumed to be monodisperse and hard-core spheres and rotational motions and torques acting on particles are neglected.

The gravitational force F_{g} can be written as

F_{g} = −E_{z},  (A3)

where E_{z} is the unit vector directed to the z axis by the notation of Eq. (2.3). It is assumed that the direction of the gravity is −z.

It is assumed that F_{c} is the impulse by prefect elastic collisions. Because collisions are inelastic in real systems, this assumption of elastic collision means our standpoint of modeling that the essential mechanism of fluidized beds is not the inelasticity in collisions but the hydrodynamic interaction. Even in this case, the model is dissipative, because the hydrodynamics interaction is nothing less than the friction. In our simulation, we represent the collision by the momentum exchange at contact in stead of the contact force F_{c}. Therefore we do not write F_{c} explicitly in the equations in the following discussion.

2. Hydrodynamic interaction

In our model, hydrodynamic interaction among particles through the fluid is considered under the Stokes approximation where the viscous effect of the fluid dominates the inertia of the fluid. The reason to adopt the Stokes approximation is as follows: The hydrodynamic interaction is the friction between the particles and the fluid and the friction is originated by the viscosity of the fluid.

In the Stokes approximation, the force acting on the particles from the fluid F_{f} and the velocity of the particles U are related by the resistance matrix R^{↔} as

− F_{f} = R^{↔} · (U − u^{∞}),  (A4)

where u^{∞} is the velocity of the fluid without particles. R^{↔} contains all information about the interaction and depend only on the configuration of particles.

To calculate R^{↔}, we adopt the method in the Stokesian dynamics, which is developed by J.F.Brady and his collaborators for dense colloidal particles [30,31,19]. In the Stokesian dynamics, R^{↔} is constructed from the two contributions in limited cases, which are the mobility matrix in dilute limit M^{↔}^{∞} and the exact resistance matrix in two-body problem R^{↔}_{2B} [32], as follows,

R^{↔} = (M^{↔}^{∞})^{−1} + R^{↔}^{lub},  (A5)

where R^{↔}^{lub} is constructed by the pairwise-additive manner from the two-body lubrication matrix R^{↔}_{2B}^{lub} which is defined by

R^{↔}_{2B}^{lub} = R^{↔}_{2B} − (M^{↔} _{2B}^{∞})^{−1}.  (A6)

M^{↔} _{2B}^{∞} is the two-body mobility matrix in dilute limit. In general M^{↔}^{∞} is formulated by the multipole expansion [30]. In the model of fluidized bed, however, we need to introduce the effect of the chamber. The chamber in the real fluidized beds has two contributions which are to bound the fluid by the vertical wall and to support the particles by the bottom.

We introduce the contribution of bounding the fluid in terms of the periodic boundary condition. Because M^{↔}^{∞} has the long-range interaction, we use the Ewald summation technique [33]. We can construct the resistance matrix under the periodic boundary condition also by Eq. (A5) only replacing M^{↔}^{∞} to the Ewald summed tensor for N particles in the unit cell [31].

On the other hand, the contribution of the bottom supporting the particles is introduced by the particles fixed in space. In this case, we get the force from the fluid to the mobile particles F_{f} in Eq. (A1) as follows,

−F_{f} = R^{↔}_{MM} · (U_{M} − u^{∞}) − R^{↔}_{MF} · u^{∞},  (A7)

where the subscripts “M” and “F” represent mobile and fixed particles respectively. The complete form of the resistance relation is

− [ \begin{matrix}F_{f} \\ F_{F} \end{matrix} ] = [ \begin{matrix}R^{↔}_{MM} & R^{↔}_{MF} \\ R^{↔}_{FM} & R^{↔}_{FF} \end{matrix} ] · [ \begin{matrix}U_{M} − u^{∞} \\ 0 − u^{∞} \end{matrix} ] ,  (A8)

where F_{F} is the force acting on the fixed particles.

For simplicity we only discuss the case without fixed particles in the following. However we can get the correct forms only the replacement of F_{f} of (A4) by (A7).

3. Effective inertia

The inertia of particles causes the relaxation process on the velocity from the initial value to the optimal value where the inertia corresponds to the relaxation time. The optimal value is usually called the terminal velocity determined on the steady state. For the fluidized beds, the terminal velocity V is determined by Eq. (A1) with dU/dt = 0 and get Eq.(2.2). This is the case that all forces acting on particles balance. With this terminal velocity we can write the relaxation process as Eq.(2.1).

Equation (2.1) is the same as used in the previous paper [15]. Although we had argued that this equation might be justified in some approximation, we would state here that (2.1) contains all of essential processes in fluidized beds. On this point, we will discuss in Appendix B.

APPENDIX B: THE INERTIAL EFFECT OF PARTICLES

Here we discuss the difficulties arising when we use Eq. (A1) with F_{f} in the Stokes approximation (A4).

If we write Eq. (A1) with F_{f}, we get

St_{0}R^{↔}^{−1} · \frac{d}{dt}U = −U + V.  (B1)

From the simulation of Eq.(B1), we observe no collision between particles even in the case with large St_{0}. Particles form a cluster and relative motions among them almost disappear.

This situation may be understood by the following model,

St_{0} \frac{dU}{dt} = −\frac{1}{r − 2a}U.  (B2)

Here we extract the radial component of the motion between two particles. r denotes the separation between the centers of the pair and U denotes the relative velocity. The resistance 1/(r − 2a) reflects the lubrication effect which diverge at contact (r → 2a). In this model, we need initially the infinite energy to approach the contact point (r = 2a) even in the case of large St_{0}. Thus the model (B2) cannot contains any collision.

Our model (2.1) can be understood as the renormalization of the singularity in the lubrication, because we get Eq. (2.1) by multiplying the singularity R^{↔} on the inertial term of (B1). (Here we note that R^{↔} is dimensionless because it scaled by 6πμa.) Our model behave reasonably like the real fluidized beds where collisions are occurred so frequently. This suggests that the singularity of the lubrication must be prevented from some mechanisms in the real systems.

We can imagine several possible mechanisms preventing the singularity. For example, if there are some dimples on the surface of the particles, they collide before the mean surfaces contact. From another point of view we can also say that the continuous description of the fluid in the gap between the particles breakdowns when particles approach closely and the gap becomes comparable to the mean-free path of the molecules of the fluid [34].

Recently a model in this context has been presented [35]. They introduce a cut-off length, which may correspond to the height of the dimples on the surface or the mean-free path of the fluid molecules. If particles approach with each other within the cut-off length, the gap between them is assumed to be the cut-off in calculation of hydrodynamic interaction. The result of their simulation is suggestive even though the situation, which is the behavior in the shearing flow without the gravity, is different. Their results are characterized by the parameter St_{s}/⟨R_{c}⟩, where St_{s} = mγ/6πμa is the Stokes number in shear flow with the shear rate γ and ⟨R_{c}⟩ is the averaged resistance of a tracer with the same volume fraction and the cut-off length.

From the above discussion, we get the meaning of the effective Stokes number St as

St = St_{0}/⟨R_{c}⟩,  (B3)

which depends on the cut-off length of the real systems. We need more delicate investigations for the dependence of the cut-off length or the mechanism preventing the singularity in the lubrication.

REFERENCES

[1] H. M. Jaeger, S. Nagel, and R. P. Behringer, Rev. Mod. Phys. 68, 1259 (1996).

[2] Granular Matter – An Interdisciplinary Approach, edited by A. Mehta (Springer-Verlag, New York, 1993).

[3] H. Hayakawa, H. Nishimori, S. Sasa, and Y-h. Taguchi, Jpn. J. Appl. Phys. 34, 397 (1995).

[4] Fluidization, 2nd ed., edited by J. F. Davidson, R. Clift, and D. Harrison (Academic Press, London, 1985).

[5] D. Gidaspow, Multiphase Flow and Fluidization (Academic Press, London, 1994).

[6] D. A. Drew, Annu. Rev. Fluid Mech. 15, 261 (1983).

[7] G. K. Batchelor, J. Fluid Mech. 193, 75 (1988).

[8] S. Sasa and H. Hayakawa, Europhys. Lett. 17, 685 (1992).

[9] M. F. Göz, J. Fluid Mech. 240, 379 (1992).

[10] T. S. Komatsu and H. Hayakawa, Phys. Lett. A 183, 56 (1993).

[11] T. S. Komatsu, Phys. Rev. E 52, 5574 (1995).

[12] Y.Tsuji, T.Tanaka, and T.Ishida, Powder Technol. 71, 239 (1992).

[13] Y. Tsuji, T. Kawaguchi, and T. Tanaka, Powder Technol. 77, 79 (1993).

[14] T. Tanaka, T. Kawaguchi, and Y. Tsuji, Int. J. Mod. Phys. B 7, 1889 (1993).

[15] K. Ichiki and H. Hayakawa, Phys. Rev. E 52, 658 (1995).

[16] K. Ichiki, Ph.D. thesis, Tohoku University, 1996.

[17] J. Frenkel, Kinetic Theory of Liquids (Oxford University Press, London, 1946).

[18] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic Press, London, 1986).

[19] J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech. 20, 111 (1988).

[20] S. Chandrasekhar, Hydrodynamics and Hydromagnetic Stability (Oxford University Press, London, 1961).

[21] H. Hayakawa and K. Ichiki, Phys. Rev. E 51, R3815 (1995).

[22] A. Einstein, Ann. Phys. (Leipzig) 17, 549 (1905).

[23] J. Furukawa and T. Ohmae, Ind. Eng. Chem. 50, 821 (1958).

[24] Y-h. Taguchi and H. Takayasu, Europhys. Lett. 30, 499 (1995).

[25] J. J. Brey, F. Moreno, and J. W. Dufty, Phys. Rev. E 54, 445 (1996).

[26] H. Hayakawa and D. C. Hong, Phys. Rev. Lett. 78, 2764 (1997).

[27] S. Warr and J. P. Hansen, Phys. Rev. Lett. (1997), in press.

[28] G. H. Ristow, G. Strassburger, and I. Rehberger, Phys. Rev. Lett. (1997), in press.

[29] A. Lippert, Chem. Ind. Tech. 38, 350 (1966); I. Muchi, S. Mori and M. Horio, Reaction-Technology in Fluidized Beds (Baihu-kan, Tokyo, 1984) (in Japanese).

[30] L. J. Durlofsky, J. F. Brady, and G. Bossis, J. Fluid Mech. 180, 21 (1987).

[31] J. F. Brady, R. J. Phillips, J. C. Lester, and G. Bossis, J. Fluid Mech. 195, 257 (1988).

[32] D. J. Jeffrey and Y. Onishi, J. Fluid Mech. 139, 261 (1984).

[33] C. W. J. Beenakker, J. Chem. Phys. 85, 1581 (1986).

[34] R. R. Sundararajakumar and D. L. Koch, J. Fluid Mech. 313, 283 (1996).

[35] A. S. Sangani, G. Mo, H.-K. Tsao, and D. L. Koch, J. Fluid Mech. 313, 309 (1996).

FIGURES

FIG. 1. Simulations performed in the parameter space (u^{∞}, St) with N_{M} = 256, N_{F} = 10 and (L_{x}, L_{y}, L_{z}) = (34, 2, 100). In this figure we show the observed phases, fixed phase (•), channeling phase (□), bubbling phase (∘) and transition phase (△). We also show the transition line of fluidization u_{c} with solid line and the area of the channel-bubble transition with dotted line, which are discussed in the text.

FIG. 2. Typical snapshots of channeling phase with St = 0.5, u^{∞} = 0.15. The time proceeds from left to right with the interval 20 dimensionless time.

FIG. 3. Typical snapshots of bubbling phase with St = 10.0, u^{∞} = 0.15. The time proceeds from left to right with the interval 20 dimensionless time.

FIG. 4. A typical behavior of E(t) in the fluidized phase. Parameters used are St = 20, u^{∞} = 0.3. The oscillation corresponds to the generation of bubbles. 1 step means 1 dimensionless time.

FIG. 5. The flow rate dependence of the averaged energy Ē(u^{∞}) with St = 10. Error bars are their standard deviation. We can see that the transition of fluidization at u^{∞} = u_{c} and the linear behavior of Ē for the flow rate u^{∞} in the fluidized phase (u^{∞} > u_{c}).

FIG. 6. St dependence of the variance V_{H} with u^{∞} = 0.3. We can see the transition around St = 5. At the channeling phase the variance is small, while at the bubbling phase the variance become large.

FIG. 7. St dependence of the averaged energy Ē(St) with u^{∞} = 0.3. We can observe a peak around the channel-bubble transition shown in Fig. 1.

FIG. 8. The convective motion of particles in the bubble with u^{∞} = 0.3, St = 10. Here two periodic images are shown. We can observe the sharp edge of the bubble and the definite convection.

FIG. 9. The convective motion of particles in the bubble with u^{∞} = 0.3, St = 100. Here two periodic images are shown. Comparing to the Fig. 8, we can observe broader edge of the bubble and weaker convection, where the scale of the velocities are the same.

FIG. 10. Effective viscosity μ_{e}(u^{∞}). This result is calculated on the simulation with St = 10.0. The fitting by the Arrhenius type function (3.10) with ε = 0.113 ± 0.017 is also shown.

FIG. 11. Effective viscosity μ_{e}(St). This result is calculated on the simulation with u^{∞} = 0.3.

FIG. 12. The non-Gaussian parameter C_{4}/(C_{2})^{2} with St = 10 for the change of u^{∞}. We also show the fitting by exp(ε′/u^{∞}) with ε′ = 0.175 ± 0.045.

FIG. 13. The non-Gaussian parameter C_{4}/(C_{2})^{2} with u^{∞} = 0.3 for the change of St.

FIG. 14. St dependence of H̄ with u^{∞} = 0.3. For the bubbling phase, we fit by Eq.(4.2) with C_{H} = 3.087 ± 0.055, D_{H} = 16.51 ± 0.14.

FIG. 15. u^{∞} dependence of the fitting parameter C_{H}. We can see C_{H} is linear for u^{∞}.