Analysis of statistical quantities in simulation of fluidized beds

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

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

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)]. In our simulation of monolayer particles, the transition of fluidization is a continuous transition. Two types of fluidized phases, the channeling phase and the bubbling phase are observed. Our simulation suggests that the flow rate plays the role of the effective temperature and the process generating free volume is important for time-averaged statistical quantities. The flow rate dependence of the diffusion constant suggests the existence of a kind of the fluctuation-dissipation relation.

05.60.+w, 03.20.+i, 82.20.Mj, 83.70.Hq

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, where fluid is injected from a distributor 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. Consequently, the system is in the fluidized phase, which contains some 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 mod-

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

els and particle-dynamics models. In the two-fluid models, particles are treated as a fluid [5–10]. These models have benefit on analytical treatments and generalization of their discussion to other systems [11]. However, their basis, such as constitution equations for the particlephase 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 such as the distinct element method [12–14]. These models, however, cannot be the basis for the twofluid models. One of the main problems 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.

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 simulation for monolayer particles in fluidized bed and discuss the 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. In Sec. IV, we discuss the interpretations of statistical quantities observed in our simulation such as the bed height and the diffusion constant. In Sec. V, we conclude our results. In Appendix, we discuss the difficulty of the introduction of particle inertia.

II. SIMULATION METHOD

In this section, we briefly explain our model and how to simulate the dynamics of granular particles in fluid flows. The detailed explanation of our model can be seen in Refs. [15,16]. For simplicity, we only consider the cases of monodispersed spherical particles. We also neglect rotational motions of particles. We define the model by the following equation of motion for particles

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

where the bold-face letters without superscripts represents vectors in 3N-dimension in the N particles’ system. For example, the velocity of particles U in Eq.(2.1) has the following components

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

where the bold-face letters with superscripts represent vectors in 3-dimension. F_{c} in Eq.(2.1) represents hardcore elastic collision among particles which is calculated by the momentum exchange for contacting particles. St in Eq.(2.1) is the relaxation time of U to the terminal velocity V determined by

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

where u^{∞} is the flow rate of induced fluid which is equal to the superficial velocity conventionally used for the fluidized beds, and −E_{z} is the unit vector directed to the gravity. R^{↔} is the resistance matrix representing the hydrodynamic interaction among particles calculated by the method of the Stokesian dynamics [17], where periodic boundary condition [18,19,15]. We also introduce some fixed particles to express the effect of the bottom of chamber. 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 particles and g̃ = g(ρ_{p} − ρ_{f})/ρ_{p} with the gravitational acceleration g, and the densities of the particle ρ_{p} and the fluid ρ_{f}. Notice that our relaxation model (2.1) is not derived from basic equations, because a naive extension of Stokesian dynamics fails to describe collisions among particles (Appendix).

In our simulation, instead of integrating Eq.(2.1), we integrate the following ordinary differential equation

St\frac{du}{dt} = exp (+\frac{t}{St}) V,  (2.4)

where u = U exp(t/St), by the 4th-order Runge-Kutta scheme with fixed time-step Δt = 0.1. In each Runge-Kutta loop, we construct the resistance matrix R^{↔} from the positions, calculate the terminal velocity V, and then obtain new positions and velocities of particles. After each step Δt, we pick up the overlapping particles and make them collide by the momentum exchange. To avoid the divergence in lubrication function for overlapping particles when we calculate R^{↔}, we introduce minimum separation r_{c} = 10^{−17}. We recorded positions and velocities of particles every 10 steps and calculate the statistical quantities which will be discussed later.

Equation (2.1) contains two control parameters, u^{∞} and St. 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 unit cell in periodic boundary condition (L_{x}, L_{y}, L_{z}) = (34, 2, 100). In this situation, particles are confined in a vertical monolayer, 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 situations come from the limitation of our 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 are observed in 3-dimensional simulations.

III. SIMULATION RESULTS

In this section, we present the results of our simulations in details. We perform simulations at the points plotted in Fig. 1 in the parameter space within the range of 0.05 ≤ u^{∞} ≤ 0.8 and 0.1 ≤ St ≤ 100.

In the fixed phase at low flow rate, particles are at rest at the bottom. When the flow rate u^{∞} exceeds the critical flow rate u_{c}, the particles begin to be fluidized. The transition between the fixed phase and the fluidized phase seems to be 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 time evolution 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.

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 a 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, where the period and amplitude of oscillation become smaller. 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 the regular behavior in E(t). From Fig. 5, we observe a continuous transition at u^{∞} = u_{c}, where the date can be well fitted by

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

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

Our results suggest that u_{c} is independent of St and Ē(u^{∞}) is linear function of u^{∞} for any St. 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 how the system is fluidized 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 in the frame of u^{∞} = 0. If the flow rate u^{∞} is smaller than U_{fall}, the mobile particles are supported by the fixed particles moving downward with u^{∞}. Therefore the system is not fluidized. 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 may grow into a bubble and propagate upward through the particles, or may construct a channel. Thus, the system is fluidized. Our picture means that the critical flow rate u_{c} is determined by the falling velocity U_{fall} which is independent of St, because U_{fall} is evaluated by configuration of particles in fixed phase as the sedimentation rate of suspensions [20] and the configuration is the same for any St. We numerically evaluate U_{fall} = 0.106 as the average of V in Eq. (2.3) for the configuration of fixed phase, which is nearly identical to u_{c}.

The consideration of the Galilei invariant properties also gives us some insight on the linearity of Ē with u^{∞}. In the coordinate with U_{F} = −u^{∞}E_{z}, it is obvious that the gain of particle energy per unit time from the gravitational potential is proportional to u^{∞}, because the center of mass also fall with the velocity u^{∞}.

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 channelbubble transitions are observed in the area shown in Fig. 1.

We also observe some qualitative changes in the statistical quantities corresponding to the channel-bubble transitions. Figure 7 shows Ē(St) for u^{∞} = 0.15 and 0.3 scaled by the value at St = 10 for each case. From this figure, 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 means that 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 may be understood as follows. From Figs. 6 and 7 in the previous paper [15] or in Fig. 4, we have seen that the kinetic energy is generated by bubbles with convective motion of particles. This is because the local difference of volume fraction of particles causes the relative motion of particles through hydrodynamic interaction. As St increases, trajectories of particles become ballistic and the difference of volume fraction tends to be leveled, since St represents particle inertia. Thus, the convective motion caused by bubbles becomes small when St becomes large. In fact, we see a clear 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. Taking into account the increase of Ē(St) in the channeling phase as a function of St, Ē(St) has a peak around the transition point.

IV. DISCUSSIONS

In this section, we present two characteristic results of our simulation, which are the height of center of mass H̄(St) and the diffusion constant D(u^{∞}), and then discuss the relation between them.

Figure 10 show H̄(St) for u^{∞} = 0.15, 0.2, 0.3 and 0.5. We see the qualitative change in the transition area in Fig. 1. The clear logarithmical increase of H̄ in the bubbling phase is remarkable. We fit the data in Fig. 10 by

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

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

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

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. 11. From this figure, C_{H} may be proportional to (u^{∞} − u_{c}). Since St is the characteristic time τ of the system, we may rewrite Eq.(4.2) as

τ ∝ exp (F \frac{ΔV}{u^{∞} − u_{c}} ) ,  (4.3)

where F is a constant. Equation (4.3) suggests that the activation process exists, where τ is the time to generate the volume expansion ΔV. The importance of activation process similar to the hole theory for simple liquids [21,22] has been indicated in compaction process of vibrating beds [23–25]. For dense particle system, it is natural that the viscosity is proportional to the time τ which needs to generate a certain free volume. Therefore, we expect that the effective viscosity μ_{e} may be written as

μ_{e}(u^{∞}) ∝ e^{ε′/(u^{∞}−u_{c})},  (4.4)

where ε′ is a constant. Equation (4.4) is similar to the Vogel-Fulcher law [26,27] in glass-forming liquids. In the experiment of fluidized beds [28], the shear viscosity measured by the modified Stormer viscometer also obeys similar behavior, although it is proposed simply the Arrhenius function exp(ε′/u^{∞}) rather than Eq.(4.4). The result of Ē(u^{∞}) in Fig. 5 suggests that the flow rate u^{∞} behaves as the effective temperature of the environment as that of the heat bath in equilibrium systems. Therefore, the critical flow rate u_{c} may correspond to the glass transition temperature.

In Fig. 12, we show D(u^{∞}) in the simulation calculated by

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

where Ũ^{α}(ω) is the Fourier transform of U^{α}(t). We see that D = 0 for u^{∞} < u_{c} and it increases monotonously for u^{∞} > u_{c}. We also try to fit the data by

D(u^{∞}) ∝ u^{∞}e^{−ε/(u^{∞}−u_{c})},  (4.6)

where ε is a fitting parameter. From Fig. 12, Eq. (4.6) is a good fitting function of the simulation results.

Let us discuss the relation between H̄(St) and D(u^{∞}). Equations (4.4) and (4.6) suggest that the diffusion constant D may be related to the viscosity μ_{e} as

D ∝ \frac{u^{∞}}{μ_{e}} .  (4.7)

Equation (4.7) suggests that the Einstein relation or the fluctuation-dissipation relation in fluidized beds may exist with the replacement of temperature by the effective temperature 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.

Before closing this section, we give some remarks. Although the transition of fluidization in some experiments seems to be the discontinuous phase transition [29,30], our simulations suggests a continuous phase transition. To remove such a contradiction, we indicate that recent experiment [31] shows that the continuous transitions are observed in the system with small particles, conventionally called A particles of Geldart’s classification [32], where the Stokes approximation and our model would be justified. It is also an open problem that at present we cannot reproduce homogeneous phase in our simulation. The reason of the lack of the phase might be the lack of the contact force among particles [30].

V. CONCLUSIONS

In this paper, we have carried out systematic simulations of fluidized bed with the change of two control parameters, u^{∞} and St. When the flow rate u^{∞} is small, particles are at rest. Above the critical flow rate u_{c}, particles are fluidized. The critical value u_{c} is independent of St. We have observed 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. The height of center of mass H̄(St) in the bubbling phase increase logarithmically with St, which suggests that the effective viscosity obey the Vogel-Fulcher law (4.4). The diffusion constant D(u^{∞}) is well fitted by Eq. (4.6) and this suggests the existence of a kind of fluctuation-dissipation relation in fluidized beds.

ACKNOWLEDGMENTS

We would like to thank S. Sasa, Y-h. Taguchi and Ooshida T. 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: THE INERTIAL EFFECT OF PARTICLES

Here we discuss the difficulties of the introduction of particle inertia into the Stokes approximation.

Let us consider directly the equation of motion of particles with only three mechanisms, which are the inertial effect of the particles, hydrodynamic interaction through the fluid under the Stokes approximation, the gravitational force. Then the equation of motion can be written as

St_{0} \frac{dU}{dt} = −R^{↔} · (U − u^{∞}) − E_{z},  (A1)

where St_{0} is the bare Stokes number defined by

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

Multiplying R^{↔}^{−1} to Eq. (A1) from left-hand side, we get

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

where V is defined by Eq.(2.3). From the simulation of Eq.(A3), 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}{2(r − 2a)}U.  (A4)

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/2(r − 2a) reflects the lubrication effect which diverge at contact (r → 2a). In this case, we need initially the infinite energy to approach the contact point (r = 2a) even in the case of large St_{0}. Thus the model (A4) cannot contains any collision.

Our model (2.1) may 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 (A3). (Here we note that R^{↔} is dimensionless because it scaled by 6πμa.) Our model behave reasonably like real fluidized beds where collisions are occurred so frequently. 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 [33].

Recently a model in this context has been presented [34]. To remove the singularity, they introduce a cutoff length, which may correspond to the height of the dimples on the surface or the mean-free path of the fluid molecules. Their simulation shows that the results are characterized by a scaled Stokes number which strongly depends on the cut-off length rather than the bare Stokes number.

We need to give the existence of renormalization procedure on particle inertia problem in Stokes approximation, which we accept a priori in Eq.(2.1).

[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, Phys. Rev. E 52, 3697 (1995).

[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, Sendai, 1996 (in Japanese).

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

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

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

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

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

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

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

[24] E. Caglioti, V. Loreto, H. J. Herrmann, and M. Nicodemi, Phys. Rev. Lett. 79, 1575 (1997).

[25] T. Boutreux and P. G. de Gennes, Physica A 244, 59 (1997).

[26] H. Vogel, Physik. Z. 22, 645 (1921).

[27] G. S. Fulcher, J. Am. Ceram. Soc. 8, 339 (1925).

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

[29] A. Lippert, Chem. Ind. Tech. 38, 350 (1966).

[30] S. C. Tsinontides and R. Jackson, J. Fluid Mech. 255, 237 (1993).

[31] G. D. Cody, D. J. Goldfarb, J. G. V. Storch, and A. N. Norris, Powder Technol. 87, 211 (1996).

[32] D. Geldart, Powder Technol. 7, 285 (1973).

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

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

FIG. 1. Plots of points in the parameter space (u^{∞}, St) where we perform simulation 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 time evolution of pattern in channeling phase at St = 0.5, u^{∞} = 0.15. The time proceeds from left to right with the interval 20 dimensionless time.

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

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

FIG. 5. The flow rate dependence of the averaged energy Ē(u^{∞}) at 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} at 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 Ē scaled by Ē(St = 10) for u^{∞} = 0.15(□) and 0.3 (∘). For each case, we can observe a peak around the channel-bubble transition.

FIG. 8. The convective motion of particles in the bubble at 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 at 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. St dependence of H̄ at u^{∞} = 0.15 (∘), 0.2 (□), 0.3 (•) and 0.5 (△). The fittings by Eq. 4.1 are also shown.

FIG. 11. u^{∞} dependence of the fitting parameter C_{H}. The solid line shows the fitting by C_{H} ∝ (u^{∞} − u_{c}).

FIG. 12. The flow rate dependence of diffusion constant D at St = 10.0. The fitting by Eq.(4.6) with ε = 0.113 ± 0.017 is also shown as the solid line.