Dynamical simulation of fluidized beds – hydrodynamically interacting granular particles

Kengo Ichiki^{(1)} ^{∗} and Hisao Hayakawa^{(1),(2)} ^{†}

(1)Department of Physics, Tohoku University, Sendai 980-77, Japan

(2)Department of Physics, University of Illinois at Urbana-Champaign, 1110 West Green Street,

Urbana, IL61801-3080, USA ^{‡} (April 7, 1995)

Abstract

A numerical simulation of a gas-fluidized bed is performed without introduction of any empirical parameters. Realistic bubbles and slugs are observed in our simulation. It is found that the convective motion of particles is important for the bubbling phase and there is no convection in the slugging phase. From the simulation results, non-Gaussian distributions are found in the particle velocities and the relation between the deviation from Gaussian and the local density of particles is suggested. It is also shown that the power spectra of particle velocities obey power laws. A brief explanation on the relationship between the simulation results and the Kolmogorov scaling argument is discussed. 05.40.+j,47.27.-i,47.55.Kf

Typeset using REVTEX

∗E-mail address: ichiki@cmpt01.phys.tohoku.ac.jp †E-mail address: hisao@engels.physics.uiuc.edu ‡Present address

I. INTRODUCTION

Recently dynamics of granular systems has attracted much attention among physicists [1–5] as a typical object of non-equilibrium statistical physics. For example, in vibrating beds [6–10], collisions among particles produce convection and turbulence. However, the fluidization of granular particles immersed in a fluid stream, where the hydrodynamic interactions are relevant, exhibits richer phenomena, illustrative of dynamical phase transitions [11–13].

In an experiment on fluidized beds, we prepare a vessel containing granular particles and impose a gas flow from the bottom of the vessel. When the flow rate is small enough, particles do not move. This state is known as a fixed bed. Above a critical value of flow rate, the fixed bed is destabilized, and then is fluidized uniformly. At larger flow rate, the uniformly fluidized bed becomes unstable and bubbles appear. Increasing the flow rate further, the bubbles become larger and then become slugs which are horizontally spread bubbles. For further increase of flow rate the state becomes disordered, and finally reaches a dilute state of particles in which is recovered spatial homogeneity. These phenomena are similar to boiling of water. The phase transitions in fluidized beds, however, are not thermal phase transitions. Thus the mechanism of the phase transitions in fluidized beds must differ from that of the boiling of water.

These dynamical phase transitions have not been observed in the systems of smaller particles in flows, such as colloid particles which are much smaller than granular particle [14]. The reason such interesting phenomena are observed in granular systems may be absence of Brownian motion. Namely, granular systems cannot reach any equilibrium states, while colloid particles can reach. Thus, granular particles are one of interesting subjects in nonequilibrium statistical physics.

We do not have any established model which is suitable to describe fluidized beds. A modern and successful approach is simulation by the distinct element method (DEM) [15–17]. The DEM is also a powerful tool to describe vibrating beds [9,10]. In this approach, the interactions among particles are replaced by a mechanical model which consists of springs, dashpots and sliders. The fluid motion is assumed to obey a phenomenological model where the hydrodynamic interactions are replaced by the coarse-grained friction between the particles and the fluid. In this way, hundreds of thousands of particles have been successfully simulated and a realistic motion of particles is reproduced. However, we stress that the DEM is not independent of the experiments because experimental results are used to choose parameters.

Two-fluid models are often used to describe fluidized beds [12,13]. These models allow one to understand macroscopic pattern formation, using bifurcation analysis and hydrodynamic stability analysis. In addition, simulations based on the two-fluid models reproduce realistic motions [13,18,19], and some authors indicate that solitons play an important role near the onset of the instability of uniformly fluidized beds [20–24]. In spite of these successful results, the two-fluid model contains some difficulties. For example, it is difficult to choose a suitable two-fluid model [25]. The role of particle motion is not clear because particles are described as a fluid. Furthermore, the two-fluid model is supplemented by empirical laws for the choice of parameters. Although two-fluid models contain a continuous approximation, the simulation of two-fluid models is not easier than the DEM.

In this paper, we perform a simulation based on the model which does not contain any empirical parameters except for the particle radius, the mass densities of the fluid and the particles, and the shear viscosity of the fluid.

For this purpose, we neglect the complexities of granular systems, which are polydispersity, direct interactions among the particles such as the Coulomb interaction and intermolecular forces, and chemical reactions induced by mixing, etc. We only treat systems which contain monodisperse spheres with hard-core interactions among the particles and with hydrodynamic interactions described by the Stokes approximation [26].

In the next section we show how to construct the model and discuss its relevance in detail. In Sec. III, it is shown that bubbles and slugs are observed in the simulation. In Sec. IV we analyze the data obtained from the simulation. For example, we show the distribution functions and the power spectra of particle velocities. In Sec. V, we briefly explain how power laws in power spectra appear. In Sec. VI, we conclude and summarize our results. We will present the details of our calculation method for hydrodynamic interactions among particles in Appendix A, and the treatment of fixed particles in Appendix B.

II. SIMULATION METHOD

In this section, we summarize the algorithm of our simulation for granular particles immersed in a fluid stream. This section consists of two parts. The first part is devoted to the general perspective of the motion of particles in a fluid. In the next part, we discuss the validity of our approximations. The hydrodynamic interactions are calculated by the Stokesian dynamics method [27–29] which is briefly described in Appendix A.

In general, classical particles with mass m in a fluid obey the Langevin equation

m\frac{d}{dt}U = F_{f} + F_{g} + F_{i} + F_{b},  (1)

where the velocity U, the position x and the forces F_{l} (l = f, g, i, b) represent vectors containing N particle elements as

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

where the superscript (i) represents the index of the particle. In Eq. (1) there are four kinds of forces: F_{f} is the drag force from the fluid, F_{g} is the gravitational force exerted on the particles, F_{i} is the force due to direct interactions among particles, and F_{b} is the Brownian force coming from the thermal motion of the fluid. Although the size distribution and the shape of particles are important factors in technology, we restrict ourselves to the motion of monodisperse spherical suspensions.

We define the following dimensionless quantities,

P e = \frac{6πμa^{2}V}{k_{B}T}  (3)

Re = \frac{ρ_{f}a^{3}V^{2}}{μa^{2}V} = \frac{ρ_{f}aV}{μ}  (4)
St = \frac{mV^{2}}{6πμa^{2}V} = \frac{2}{9} \frac{ρ_{p}aV}{μ}  (5)

where P e,Re and St are respectively the Péclet number, the Reynolds number and the Stokes number. In Eqs.(3)-(5), a, ρ_{p}, ρ_{f} and μ are the particle radius, the mass densities of the particle and the fluid, and the shear viscosity of the fluid respectively, k_{B} is the Boltzmann’s constant. V is the characteristic velocity, where we choose the sedimentation rate of one particle

U_{0} = \frac{mg̃}{6πμa} = \frac{2a^{2}ρ_{p}g̃}{9μ}  (6)

with the effective gravitational acceleration g̃ = g(ρ_{p} −ρ_{f})/ρ_{p}. These dimensionless numbers actually have definite meanings. The Péclet number is the ratio of work done by the drag force over the size of a particle (6πμa^{2}V) to thermal energy (k_{B}T). The Reynolds number is well known as the ratio of inertia to drag. The Stokes number represents the relative importance of the kinetic energy of particles (mV^{2}) to work done by the drag force. It is also recognized as the ratio of the time scales,

St = \frac{T_{r}}{T_{p}} ,  (7)

where T_{r} = m/6πμa is the relaxation time of particle velocity due to the drag force, and T_{p} = a/V is the passing time of the particle scale a with the velocity V.

The Froude number, which has been widely used on this problem, is given by

F r = \frac{V^{2}}{g̃L},  (8)

with L being the linear size of system. The meaning of the Froude number is the ratio of the kinetic energy to the gravitational potential, and F r is proportional to St if we adopt V = U_{0}. For later convenience, we also introduce the effective Reynolds number of particles as

Re_{(p)} ≡ \frac{VL}{ν_{p}} = \frac{9L}{2a}St,  (9)

where ν_{p} = μ/ρ_{p}. The Stokes number is related to the particle size, while the Froude number and the effective Reynolds number Re_{(p)} are related to the system size. The role of Re_{(p)} will be discussed in Sec. V. Thus P e,Re and St are fundamental and independent parameters to determine the motion of particles.

Let us consider the explicit form of forces in (1) and discuss their relevance. The gravitational force F_{g} acting in the z direction is simply given by

F_{g} = −mg̃E_{z},  (10)

where E_{z} is the generalization of unit vector e_{z} as in (2). The direct interaction among particles F_{i} is assumed to be due to hard-core interactions. Therefore it can be treated as exchanging velocities during elastic collisions among particles. We neglect the random force F_{b}, because in this paper we treat the system where the Péclet number is large enough. The validity of P e ≫ 1 will be discussed later.

The most relevant and complicated force is F_{f}, which is determined by the Navier-Stokes equation

ρ_{f} {∂_{t}u + (u · ∇)u} = μ∇^{2}u − ∇p,  (11)

with the incompressible condition

∇ · u = 0.  (12)

The incompressibility is valid even if the fluid is air, when the particle motion is much slower than the sound velocity of fluid. When Re ≪ 1, the Navier-Stokes equation is reduced to the Stokes equation

−μ∇^{2}u + ∇p = 0.  (13)

In this case the particle velocity U is connected with the force induced by the particles on the fluid F_{p} in the following linear relation [30],

F_{p} = R (x) · (U − u^{∞}) ,  (14)

where u^{∞} is the fluid velocity in the absence of particles and R is the 3N × 3N resistance matrix which depends only on the particle configuration x. The details of the construction of R are described in Appendix A. The force F_{p} induced by the particle is related to F_{f} by

F_{p} = −F_{f}.  (15)

In later discussions, we assume Re ≪ 1. As we will show, this approximation is valid for the motion of relatively small particles.

In the case of Re ≪ 1 and P e ≫ 1 and with hard-core interactions, the Langevin equation (1) can be reduced to

St\frac{d}{dt̂}Û(t̂) = −R̂(x(t̂)) · (Û(t̂) − û^{∞}) − E_{z}.  (16)

Here we scale the velocities by V, the length by the radius a and the resistance matrix by the drag factor 6πμa. We denote the dimensionless value of U as Û. From (16) it is obvious that the Stokes number St is an important parameter for our model. The role of St is also pointed out for dilute monodisperse suspensions with low Reynolds number [31,32]. For vibrating beds in which the particle inertia is dominant and the hydrodynamic interaction is negligible, St should be large while for liquid-fluidized beds St is small.

In our simulation, we need to integrate Eq.(16) for a small time interval numerically. We divide Û(t̂) into two parts,

Û(t̂) = Û_{0} + Û_{1}(t̂).  (17)

Û_{0} is determined by

0 = −R̂(x) · (Û_{0} − û^{∞}) − E_{z}.  (18)

It is obvious that U_{0} is the same solution as that for St ≪ 1. Although Û_{1}(t̂) is determined by

St\frac{d}{dt̂}Û_{1}(t̂) = −R̂(x) · Û_{1}(t̂),  (19)

we adopt the simplest form for computational efficiency,

Û_{1}(t̂) = Û_{1}(0) exp (−\frac{t̂}{St}) ,  (20)

which is the solution of (19) on the assumption of R̂ ≃ I, where I is the unit tensor. From the initial condition Û(0) = Û_{0} + Û_{1}(0), we thus obtain

Û(t̂) = Û_{0} + (Û(0) − Û_{0}) exp (−\frac{t̂}{St}) .  (21)

This solution shows us that the particle velocities are damped from the initial value Û(0) to the terminal velocity Û_{0}. The simplification in Eq.(20) is crucial. We expect, however, the error from this simplification is small when we choose a small enough time interval for the numerical integration, because the resistance matrix is calculated as a function of particle configuration at each step.

Let us estimate the values of the dimensionless parameters. We adopt ρ_{p} = 2.5[g cm^{−3}] which is a typical value of glass beads. Fluids are assumed to be air (ρ_{f} = 1.2×10^{−3}[g cm^{−3}] and μ = 1.82 × 10^{−4}[g/cm sec]) and water (ρ_{f} = 1.0[g cm^{−3}] and μ = 1.0 × 10^{−2}[g/cm sec]) at room temperature. Substituting k_{B} = 1.38 × 10^{−16}[erg K^{−1}], T = 293[K] = 20[^{∘}C] and g = 981[cm sec^{−2}], we obtain Table I and Fig. 1. We also show the characteristic times, T_{r} and T_{p}, in (7).

The lines in Fig. 1 connect the points where dimensionless parameters are equal to unity for air and water. Therefore, there is no meaning in the vertical axis. In the left of the line of P e = 1 in Fig. 1, the random force is dominant and the particles can be regarded as typical colloidal suspensions. On the other hand, the particles to the right of P e = 1 can be regarded as typical granular particles where the hydrodynamic interaction among particles dominates the Brownian force. In the left region of Re = 1, the fluid motion can be described by the Stokes equation. On the other hand, in the right of Re = 1, the inertia term in the Navier-Stokes equation is important. In the left region of St = 1 where T_{r} < T_{p}, we can neglect the inertia effect of particles. Thus, the particles move with their terminal velocities. In the right of St = 1 where T_{r} > T_{p}, the effect of particle inertia exceeds the drag from the fluid. Therefore, collisions among the particles are important.

We choose the particle radius as 10μm because we are interested in the case of P e ≫ 1 and Re ≪ 1. In this case the dimensionless parameters have the values presented in Table I for air and water. As anticipated in Fig. 1, the difference between air and water appears in the value of the Stokes number. We now interested in the phenomena driven by the drag force, which is the behavior in the time scale T_{p}. Therefore, for gas-fluidized beds, the collisions among particles are not negligible, while for liquid-fluidized beds, the collisions are not important. We focus on gas-fluidized beds in this paper. We will discuss the properties of liquid-fluidized beds elsewhere.

We notice that the particle radius a ≃ 10^{−3}[cm] belongs to the group C in the classification by Geldart [33,13] in which all cohesive powders are difficult to fluidize. This difficulty in experiments arises because direct interparticle forces dominate the hydrodynamic drag force. If the interparticle forces, except for the hard-core interactions can be removed, a collection of small particles must be fluidized as will be shown.

Finally we summarize the adopted assumptions.

1. the diameter of all particles is identical.

2. the interaction among particles is described by hard-core interactions.

3. the thermal motion is negligible.

4. the inertia of the fluid is negligible.

Assumptions 3 and 4 are valid when the particle radius is between 1[μm] and 10[μm] in both air and water.

III. SIMULATION RESULTS

In this section, we collect three typical results of simulations. In many of our simulations, the configuration of particles is restricted to lie in a plane parallel to the direction of gravity for efficiency of calculation, although periodic boundary conditions are also applied in the direction perpendicular to this plane. Thus, in such cases, particles are influenced by threedimensional hydrodynamic interactions.

At first we show the result of our simulation using rectangular cells with periodic boundary conditions in which the ratio of the height to the width is large compared with unity. We observe stable slugs which move upward in this simulation (Fig. 2). This result is also observed in the three-dimensional (not monolayer) simulation with slender cells as well (Fig. 3). From these results, particles in the dilute region relatively fall down, so that in the denser region the sedimentation velocity is smaller. This tendency agrees with the standard theory [34–36] and the experiments of sedimentation [14,37]. In our simulation for slugs the configuration of particles is nearly closed packed and particles are immobile in concentrated regions. On the other hand in the dilute region, particles monotonically fall down. We cannot observe any convection of particles in this simulation. Thus slugs can be produced by a pseudo one-dimensional motion of particles, after the uniform state has become unstable.

Next we show the result of our simulation using square unit cells with periodic boundary conditions. Figure 4 is a typical snapshot. In contrast to the previous result, we see that particles in the dilute region float up, while particles in the dense region fall down. Thus particle convection exists around the dilute region, which may be regarded as a bubble. This kind of convection inside bubbles has been observed in experiments [38]. Thus, we infer that convection is important to create bubbles. This bubble, however, is not stable and it will disappear soon afterwards. We can also see periodic birth and death processes of bubbles.

We have also performed simulations introducing fixed particles in our system. The treatment of fixed particles is described in Appendix B. In real systems, particles are settled in a vessel. For this purpose, we introduce particles fixed at the bottom of the unit cell. In our simulation, the fixed particles are placed horizontally, spaced as far as 3.5 radius apart. We perform this simulation as follows. At first we randomly position free particles and allow them to drift to the bottom under the influence of gravity. Thus the free particles fall down and produce a fixed bed. Next we inject the flow upward with the velocity u^{∞}. The result of our simulation after the injection of flow is shown in Fig. 5. We notice that the unit cell drawn in Fig. 5 is connected with its mirrors in all directions as in the previous figures.

From Fig. 5, we observe that at first the fixed bed floats like a single cluster, but the cluster becomes unstable. Then the fixed bed becomes a fluidized bed. In this fluidized bed, we observe the formation of bubbles periodically at the bottom, which float up, and travel through the bed. We also observe particle convection around bubbles. This kind of bubble formation from fixed beds is similar to that in the real experiments [38] and large size simulation [11,15–17].

Figure 6 shows the standard deviation of the particle velocities and Figure 7 shows the number of particles in the region at 8/25 of the cell height above the bottom. These figures show that bubble formation occurs periodically at the peaks of the standard deviation.

At the end of this section, we summarize important characteristics in fluidized beds. When we compare the results in Figs. 2 and 4 or 5, we observe that the particle motion has different characteristics. In Fig. 4, particles float up in dilute regions and fall down in concentrated regions. This tendency can be understood as follows. The concentrated regions can be regarded as clusters, where flow cannot penetrate inside clusters. Thus clusters may have fast sedimentation velocity as in the definition of V ∝ a^{2} with the radius a. On the other hand, slugging motion is understood from the standard theory of sedimentation of homogeneous suspensions, where dilute regions have larger sedimentation rate than dense regions. In Fig. 2 we see no isolated clusters, and particle distribution for horizontal direction is almost uniform. If particles are dispersed uniformly, fluid flow feels larger drags from many particles than that from small particles. In other word, the difference between the bubbling phase and the slugging phase is whether or not the convection of particles exists. To produce convections, we need both characteritics of clusterings and sedimentation. Thus, the particle motions in real systems are determined by complex combinations of two different characteristics.

IV. ANALYSIS

In this section, we analyze the data obtained from our simulation. At first we discuss the velocity distribution functions (VDFs), where not only Gaussian-like but non-Gaussian distributions are observed and close relationship between the deviation from Gaussian and the local density of particles is indecated. Next we discuss the power spectra obtained from the data, which suggest the existence of tails obeying power laws. These power law tails may correspond to those reported by Taguchi [39], and taken as an evidence for powder turbulence [40].

We display VDFs in our simulation for systems with the square periodic cells in Fig. 8. From this figure, it is obvious that the VDFs are near Gaussian but anisotropic. In this case, the vertical VDF has two branches, each obeying a different Gaussian distribution. This transmission in the VDF has also been observed in the simulation of the two-fluid model [19]. The anistropic poperty of VDF indicates that the system does not have any local equilibrium.

For the systems with fixed particles, the form of VDFs is similar to an exponential one (Fig. 9). It is interesting that the introduction of fixed particles produces an exponential VDF. In this system as shown in Fig. 6, active states which has large standard deviation of particle velocities and inactive states of motion of particles emerge in turn. From Figs. 6 and 7, we indicate that the active state corresponds to the appearance of a bubble and in the inactive state contains is no bubble. To see the quantitative difference between two states, we investigate VDF in each state separately. Here we define the active state as the time region where the standard deviation of particle velocities is more than 0.4U_{0}, and the inactive state as the period where the standard deviation is less than 0.3U_{0}. Figure 10 is the VDFs for the active and the inactive state. From this figure, VDF in the inactive state is close to an exponential distribution, while that in the active state is a Gaussian distribution.

For the systems with slugs, VDFs are far from Gaussian futher and also different from simple exponential distribution (Fig. 11). We attempt to fit them to t-distribution,

f(U) ∼ (1 + aU^{2})^{−b} ,  (22)

where a and b are parameters. This t-distribution is observed by Taguchi and Takayasu [41] for VDFs in vibrating beds and by Sinai and Yakhot [42] for passive scalars in turbulence. From this figures, although the horizontail VDF can be fitted by the t-distribution, the vertical VDF is far from any t-distribution. However the tail for negative U_{y} can be understood by the following simple picture. In the slugging state, most of particles are included in the dense region and only a few particles are falling in the slug. Here we assume that falling particles start from zero relative velocity to the dense region. The particles are accelerated by the gravity then they collide with the top of the dense region and join in it. If this assumptions are valid, the VDF of falling particles is uniform in the range between the initial and the terminal falling velocities, which may correspond to the tail for U_{y} observed in Fig. 11 . To investigate the VDF in the dense region, we substract the uniform distribution from the original VDF (Fig. 11 (b)) in the range of −0.4 < Û < −0.115. Here the probability of the uniform distribution is approximated by the original VDF in the range of −0.7 < Û < −0.4, and Û = −0.115 is the location of the peak of the original VDF. The resultant VDF is shown in Fig. 12 in the range of −0.4 < Û < 0.3. We can fit this well by (22) for each side of the peak, although there is still anisotopy in the figure. Thus, t− distribution seems to be applicable to our system.

From our results of VDFs, we may indicate that non-Gaussian properties are closely related to the local density of particles. Because in relatively dilute case such as the system with square unit cells and the active state in the system with fixed particles, VDFs are close to a Gaussian, while in the dense case such as the inactive state in the system with fixed particles and the system which has slugs, VDFs are far from Gaussian. These two characteristics may be understood as follows. In active states, the density of particles are relatively low. As a result, particles can move almost freely without influence of lubrication effects. Then collisions among particles occur at random, and inelastic effects from viscous terms which is mainly from the lubrication force may be suppressed. On the other hand, in inactive states, the density of particles is extremely high, and the lubrication effect also becomes important. Therefore, the inelastic effects dominate the randomness to produce Gaussian distributions.

It is interesting that we observe non-Gaussian, exponential-like, VDFs in our systems. Non-Gaussian probability distributions are observed in various systems. In fluid turbulence non-Gaussian distributions are found in the probability distribution functions of velocity differences [43] and passive scalars [44,45]. In simulations of vibrating beds, VDFs of particles are found to be described by a non-Gaussian VDFs [41]. For the simple models consisting of hard-spheres, VDFs also show non-Gaussian distribution [46,47]. In astronomy non-Gaussian VDFs are also reported [48].

We also indicate that in the region obeying non-Gaussian VDFs, the kinetic temperature which is defined through the deviation of the Gaussian VDF cannot be used. Although it is possible to introduce the granular temperature from the different context, we need to be aware of diffrence between this granular temperature and the usual kinetic temperature.

Next, we investigate the power spectra in frequency E(ω) defined through

E(ω) = \frac{1}{N} ∑_{α}^{N} < Ũ^{α}(ω) · Ũ^{α∗}(ω) >,  (23)

where

Ũ(ω) = ∫_{−∞}^{∞} dt e^{−iωt}U(t).  (24)

The results shown in Figs. 13, 14 and 15 are obtained by a standard fast Fourier transform routine with the Parzen window [49]. All of these figures indicate that there are three regions; the spectra seems to be white in low ω, there are some peaks in ther middle, and the spectrum obeys a power law in high ω. To understand the mechanism to make three regions, we consider the following three characteristic frequencies,

ω_{L} ≡ 2π \frac{Ū}{L} ,  (25)
ω_{p} ≡ \frac{2π}{T_{p}} ,  (26)
ω_{r} ≡ \frac{2π}{T_{r}} ,  (27)

where we adopt the average particle velocity relative to the fluid for the systems without fixed particles and the induced flow rate for the system with fixed particles as Ū. T_{p} and T_{r} introduced in (7) are the passing time of the particle scale and the relaxation time respectively. These frequencies are also shown in the figures. From these figures, ω_{L} and ω_{r} seem to correspond to boundaries of three regions; the white spectrum region, the region with some peaks, and the region obeying power law. The peak near ω/ω_{0} ≃ 13 in Fig. 15 corresponds to the frequency of bubble formations, where ω_{0} = 2π/2048 × 10^{−4}[sec]is the smallest frequency produced from our entire simulation. In the higher frequency range, all of these figures show the existence of a power law E(ω) ∼ ω^{β} between ω_{r} and ω_{p} whose exponent is fluctuated between −1.49 and −1.63. This is not far from the Kolmogorov spectrum E(ω) ∼ ω^{−5/3} in fluid turbulence. About the power laws, we will discuss futher in the next section.

V. DISCUSSION

Now we discuss the origin of the power laws observed in the power spectra, which are similar to Kolmogorov scaling [50]. Taguchi [39] also observed the Kolmogorov-like scaling in vibrating beds. We briefly explain the relationship between our observed spectra and the Kolmogorov scaling [51]. In this sense, we need to clarify whether the power law observed in our simulations is the same as that of Taguchi.

Our system contains an energy source on the scale of the lowest wave number, and the energy is dissipated in the fluid on the scale of the highest wave number. Vibrating beds have also common feature. Therefore, we might expect to have a cascade process, as assumed by Kolmogorov and proposed by Richardson, where the energy dissipation rate ϵ determines the statistical properties of small particle motion. For this purpose, we assume that the motion of particles is determined by only the energy dissipation rate ϵ and the viscousity μ. As discussed in Sec. II, the particle motion is described by only one relevant dimensionless parameter, the Stokes number, St = 2ρ_{p}aV/9μ in the case of P e ≫ 1 and Re ≪ 1. Although the Stokes number may be regarded as the effective Reynolds number for the particle fluid, it contains only the energy dissipation length scale. Instead of St it is convenient to use the effective Reynolds number for particle fluid Re_{(p)} ≡ VL/ν_{p} = 9LSt/2a defined in Eq.(8). From dimensional analysis, we obtain the following scaling for the energy spectrum

E(k) = ϵ^{1/4}ν_{p}^{5/4}F (\frac{k}{k_{d}} , Re_{(p)}) ,  (28)

where F (x, y) is a dimensionless function, k_{d} is the dissipation scale which is the maximum of a^{−1} and the Kolmogorov scale ϵ^{1/4}ν_{p}^{−3/4}. In the region of high wave numbers we might guess that there is a balance between energy injection and cascade process, where the dissipation is not important. Therefore, we obtain

E(k) = ϵ^{2/3}k^{−5/3}.  (29)

This result is equivalent to that by Kolmogorov [50].

In order to ensure the scaling ansatz, we need to impose the condition

Lk_{d} ∼ Re_{(p)}^{3/4} ≫ 1.  (30)

Equation (30) means that the inertia spectrum can be observed in our system (Table II). Although there is an ambiguity in choosing L, it is interesting that the values of Re_{(p)}^{3/4} correspond to the strength of the non-Gaussain VDFs. We also expect the cutoff a^{−1} is larger than the dissipation scale ϵ^{1/4}ν_{p}^{−3/4}, because we do not observed any dissipation range (see also [39]). We also comment on the reason to observe turbulent charactersitics in relatively small Re_{(p)}, where in usual fluid, the value of Re_{(p)} in Table II is not enough large to show turbulent charcateristics. We guess that inelatic scatterings of particles in our system enhance chaotic characteristics, while the molecules in pure fluid only have elastic scatterings. We can use this argument to vibrating beds, where μ is replaced by the effective viscosity of particle flow.

We only discuss the spectrum for E(k) and not for our observed E(ω). Although the general relationship between E(k) and E(ω) is unclear, Viecelli [52] demonstrates

E(ω) ∼ ω^{−5/3}.  (31)

for the inertia range in both Euclerian and Lagrangian coordinate frames can be derived from (29) with the assuption of space-time symmetry. His argument is independent of the basic equations and assumes that the system can be charcaterized only through the kinetic viscoisty and ϵ as in our argument. Therefore, we may expect that his argument is applicable to our case.

We also comment on the reason that we observe the Kolmogorov like scaling even in two-dimensional systems. In our case the enstrophy is not a conserved quantity in two dimensional case. Therefore it is not surprizing that we observe the Kolmogorov-like scaling in two-dimensional systems, because the derivation of (29) does not contain any information about the spatial dimensionality.

VI. CONCLUSION

In this paper, we have succeeded in simulating granular systems with fluid flow without any empirical parameters. Our simulation reproduces realistic slugs and bubbles. According to the results of our simulation, we confirm that particle convection is important for bubble formation. From our simulation we observe several velocity distribution functions, Gaussian and non-Gaussian or exponential distribution. We also observe power spectra obeying power laws. We briefly explain the relationship between our observed spectra and Kolmogorov scaling. In this sense, our paper has confirmed the universality of the powder turbulence proposed by Taguchi [40].

Before closing this paper, we comment on our approximation in using the Stokes equation for fluid flow. In real experiments and simulations by engineers, the particle radius, typically a ≃ 10^{−1} [cm], is often much larger than that we assume (a ≃ 10^{−3} [cm]). This is because the direct interaction among particles except for hard-core interactions dominates all other forces, and it is difficult to observe fluidized beds in such small particles [13,33]. For the fluidized beds with such large particles, we need to consider the advection term in the fluid motion and the drag which should contain a term in proportion to the square of the difference of velocities between fluid and particles [13]. To simulate this system without introduction of empirical laws is almost impossible. We, however, have shown that these kinds of complexities are irrelevant for the formation of bubbles and slugs. The essence of the physics in fluidization can be understood when we drop irrelevant complexities such as direct interparticle forces and the advection term in fluid flows.

ACKNOWLEDGEMENTS

The authors thank fruitful discussions with T.S. Komatsu, Y-h. Taguchi, S. Sasa, H. Takayasu and H. Nishimori. We also appreciate the useful comments by H.J. Herrmann and A.S. Sangani, and are grateful to N. Goldenfeld for his critical reading and comments. One of the authors (K.I.) thanks J.F. Brady for explaining the details of his calculations to the author. This work was partially supported by the Foundation of Promotion for Industrial Science and the U.S. National Science Foundation through grant number NSF-DMR-93-14938.

APPENDIX A: HYDRODYNAMIC INTERACTION

In this section, we explain how to treat hydrodynamic interactions. When we adopt the Stokes approximation, the calculation of the hydrodynamic interactions is equivalent to constructing the resistance matrix. Our aim is to calculate the resistance matrix with acceptable accuracy taking into account computational efficiency. We thus adopt the Stokesian dynamics developed by Brady and his coworkers [27–29] to describe colloidal dispersions. They distinguish hydrodynamic long-range interactions from the lubrication force which represents short range hydrodynamic repulsive interactions. For the long-range part, we use the multipole expansions, while for the lubrication part, we use the exact solution of two-body problem [53,54] as

[ \begin{matrix}F^{(1)} \\ F^{(2)} \end{matrix} ] = R_{2B} · [ \begin{matrix}U^{(1)} − u^{∞} \\ U^{(2)} − u^{∞} \end{matrix} ] ,  (A1)

where u^{∞} is the fluid velocity without particles. The lubrication part is separated from the exact solution as

R_{2B}^{lub} = R_{2B} − (M_{2B}^{∞})^{−1} ,  (A2)

where M_{2B}^{∞} is the Rotne-Prager tensor [55], which represents the long-range interaction. Let R^{lub} be the linear combination of all possible pairs of R_{2B}^{lub}. Then, an approximate resistance matrix is represented by

R = (M^{∞})^{−1} + R^{lub}.  (A3)

The validity of this method has been confirmed in various numerical experiments for colloid systems. This idea is also valid in the theoretical calculation of the sedimentation rate [36].

For the long-range interaction, we apply periodic boundary conditions to describe the system containing infinite number of particles. For simplicity, we consider only the contribution of force. Note that higher order corrections from torque have been discussed by Brady et al. [29]. Beenakker [56] has shown that the mobility matrix can be written Ewald’s summation as

6πμa(U_{i}^{α} − u_{i}^{∞}) = F_{i}^{α} + ∑_{γ} \begin{matrix}N \\ ∑ & ′ \\ β=1 \end{matrix} M_{ij}^{(1)}(x^{α} − x^{β} + r_{γ})F_{j}^{β}
+\frac{1}{V} ∑_{λ≠0} ∑_{β=1}^{N} cos (k_{λ} · (x^{α} − x^{β})) M_{ij}^{(2)}(k_{λ})F_{j}^{β} − M_{ij}^{(2)}(r = 0)F_{j}^{α}, (A4)

where N is the number of particles in the unit cell and M_{ij}^{(1)}(r), M_{ij}^{(2)}(k) and M_{ij}^{(2)}(r = 0) are respectively given by

M_{ij}^{(1)}(r) = erfc(ξr) \{δ_{ij} (\frac{3}{4} \frac{a}{r} + \frac{1}{2} \frac{a^{3}}{r^{3}} ) + r̂_{i}r̂_{j} (\frac{3}{4} \frac{a}{r} − \frac{3}{2} \frac{a^{3}}{r^{3}} )\}

+\frac{e^{−ξ^{2}r^{2}}}{\sqrt{π}} [δ_{ij} (4a^{3}ξ^{7}r^{4} + 3aξ^{3}r^{2} − 20a^{3}ξ^{5}r^{2} − \frac{9}{2}aξ + 14a^{3}ξ^{3} + a^{3}ξ \frac{1}{r^{2}} )

+r̂_{i}r̂_{j} (−4a^{3}ξ^{7}r^{4} − 3aξ^{3}r^{2} + 16a^{3}ξ^{5}r^{2} \frac{3}{2}aξ − 2a^{3}ξ^{3} − 3a^{3}ξ \frac{1}{r^{2}} )] ,  (A5)
M_{ij}^{(2)}(k) = 6πa (δ_{ij} − k̂_{i}k̂_{j}) \frac{1}{k^{2}} (1 − \frac{1}{3}a^{2}k^{2}) (1 + \frac{k^{2}}{4ξ^{2}} + \frac{k^{4}}{8ξ^{4}} ) e^{− \frac{k^{2}}{4ξ^{2}}} ,  (A6)

and

M_{ij}^{(2)}(r = 0) = \frac{δ_{ij}}{\sqrt{π}} (6aξ − \frac{40}{3}a^{3}ξ^{3}) .  (A7)

Here erfc(x) is the complimentary error function given by

erfc(x) = \frac{2}{\sqrt{π}} ∫_{x}^{∞} exp(−z^{2})dz,  (A8)

and ξ with units of inverse length is an arbitrary parameter which we choose to minimize the number of lattice-sums in both real space and k-space. In the simulation, we use ξ = \sqrt{π}/L̄ where L̄ is the average of the length of the periodic cell. The suffix γ = (n_{1}, n_{2}, n_{3}) represents a periodic cell in the real space and r_{γ} is the lattice vector given by

r_{γ} = (n_{1}L_{1}, n_{2}L_{2}, n_{3}L_{3}) ,  (A9)

where L_{1}, L_{2}, L_{3} denotes the length of unit cell in each direction. The wave vector k_{λ} in the reciprocal cell λ = (m_{1}, m_{2}, m_{3}) is given by

k_{λ} = ( \frac{2πm_{1}}{L_{1}} , \frac{2πm_{2}}{L_{2}} , \frac{2πm_{3}}{L_{3}} ) .  (A10)

The summation ^{∑′} M_{ij}^{(1)} in Eq.(A4) means that the contribution of the sum in which α = β at γ = 0 is eliminated. The contribution of k = 0 in the lattice sum in reciprocal space is canceled by that of the average force acting on the fluid [29].

We have checked the validity of our program to examine several sedimentation velocities under regular configurations like simple cubic, body-center cubic and face-center cubic. We compared our results with the Stokesian dynamics by Brady et al. [29] in which we choose the corresponding model, neglect higher order moment of hydrodynamic interactions, and with the exact solution by Zick and Homsy [57]. From this test, we recover the corresponding result by Brady et al. [29], that is, our program seems to work correctly. Furthermore, our result is very close to the exact solution by Zick and Homsy [57] except for extremely concentrated region. Thus the simplification of neglecting higher order moments does not cause serious problems.

APPENDIX B: FIXED PARTICLES

In this section, we show how to calculate the terminal velocities when we introduce particles fixed in the space. If we obtain the terminal velocities, the motions of free particles are determined by (21).

If we know the resistance matrix which contains both free particles and fixed particles, we obtain the following equation,

[ \begin{matrix}F_{m} \\ F_{f} \end{matrix} ] = [ \begin{matrix}R_{mm} & R_{mf} \\ R_{fm} & R_{ff} \end{matrix} ] · [ \begin{matrix}U_{m} − u_{m}^{∞} \\ U_{f} − u_{f}^{∞} \end{matrix} ] ,  (B1)

where the suffix m represents free particles and the suffix f represents fixed particles. We already know the force acting on free particles F_{m}, the velocity of fixed particles U_{f} = 0 and the velocity of induced fluid u^{∞}. (B1) can be solved for U_{m} as

U_{m} − u_{m}^{∞} = R_{mm}^{−1} · (F_{m} + R_{mf} · u_{f}^{∞}) .  (B2)

Therefore the terminal velocities of free particles can ben represented by known variables. This procedure is applicable to systems with periodic boundary conditions. In fact, from Eq.(A4) the mobility matrix can be composed of particles in the unit cell.

REFERENCES

[1] H. M. Jaeger and S. Nagel, Science 255, 1523 (1992).

[2] Y-h. Taguchi, H. Hayakawa, S. Sasa, and H. Nishimori eds., Int. J. Mod. Phys. B 7, 1757 (1993).

[3] Disorder and Granular Media, edited by D. Bideau and A. Hansen (North-Holand, Amsterdam, 1993).

[4] Powders and Grains 93, edited by C. Thornton (Balkema, Rotterdam, 1993).

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

[6] P. Evesque, E.Sznatula, and J.-P. Denis, Europhys. Lett. 12, 623 (1990).

[7] P. Evesque and J. Rajchenbach, Phys. Rev. Lett. 61, 44 (1989).

[8] C. Laroche and S. Douady, J. Phys. (France) 50, 699 (1989).

[9] Y-h. Taguchi, Phys. Rev. Lett. 69, 1367 (1992).

[10] J. A. C. Gallas, H. J. Hermann, and S. Sokołowski, Phys. Rev. Lett. 69, 1371 (1992).

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

[12] D. Gidaspow, Appl. Mech. Rev. 39, 1 (1986).

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

[14] W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions (Cambridge University Press, Cambridge, 1989).

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

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

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

[18] J. Ding and D. Gidaspow, AIChE J. 36, 523 (1990).

[19] T. S. Komatsu, in preparation.

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

[21] H. Hayakawa and S. Sasa, in Complex Fluids, edited by L. Garrido (Springer, Berlin, 1993), p. 319.

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

[23] H. Hayakawa, T. S. Komatsu, and T. Tsuzuki, Physica A 204, 277 (1994).

[24] S. Harris and D. Crighton, J. Fluid Mech. 266, 243 (1994).

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

[26] Preliminary modeling based on Stokesian dynamics has been proposed by K.Ichiki and H.Hayakawa, Int. J. Mod. Phys. B 7, 1899, (1993). However, the treatment by them is not complete, because they did not consider periodic effects as explained in Appendix A.

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

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

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

[30] J. Happel and H. Brenner, Low Reynolds number hydrodynamics (Martunus Nihhoff, Dordrecht, 1973).

[31] D. L. Koch, Phys. Fluids A 2, 1711 (1990).

[32] D. L. Koch, Phys. Fluids A 4, 1337 (1992).

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

[34] G. K. Batchelor, J. Fluid Mech. 52, 245 (1972).

[35] J. F. Brady and L. J. Durlofsky, Phys. Fluids 31, 717 (1988).

[36] H. Hayakawa and K. Ichiki, (1995), to be published in Phys. Rev. E. RC. (adaporg/9502003).

[37] R. H. Davis and A. Acrivos, Ann. Rev. Fluid Mech. 17, 91 (1985).

[38] S. Mori, private communications.

[39] Y-h. Taguchi, Europhys. Lett. 24, 203 (1993).

[40] Y-h. Taguchi, Physica D 80, 61 (1995).

[41] Y-h. Taguchi and H. Takayasu, (1995), submitted to Europhys. Lett. (chaodyn/9406012).

[42] Y. G. Sinai and V. Yakhot, Phys. Rev. Lett. 63, 1962 (1989).

[43] C. W. Van Atta and W. Y. Chen, J. Fluid Mech. 44, 145 (1970).

[44] M. Sano, X. Z. Wu, and A. Libchaber, Phys. Rev. A 40, 6421 (1989).

[45] B. Castaing et al., J. Fluid Mech 204, 1 (1989).

[46] Y-h. Taguchi and H. Takayasu, (1995), submitted to Phys. Rev. Lett. (adaporg/9501003).

[47] C. Jarzynski and W. J. Swiatecki, Nucl. Phys. A552, 1 (1993).

[48] M. S. Miesch and J. M. Scalo, (1995), submitted to Astrophys. J. (astro-ph/9412042).

[49] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (Cambridge University Press, Cambridge, 1988).

[50] A. N. Kolmogorov, C. R. Acad. Sci, USSR 30, 301 (1941).

[51] Needless to say, the turbulent like characteristics is observed in motion of granular particles. Therefore, there is no wonder that turbulent like characteristics can produced by motion of particles obeying the Stokes equation. Recall the inertia of particles is important and the effective Reynolds number for particle fluid is large.

[52] J. A. Viecelli, Phys. Fluid A 3, 2698 (1991).

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

[54] S. Kim and R. T. Mifflin, Phys. Fluids 28, 2033 (1985).

[55] J. Rotne and S. Prager, J. Chem. Phys. 50, 4831 (1969).

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

[57] A. A. Zick and G. M. Homsy, J. Fluid Mech. 115, 13 (1982).

FIGURES

FIG. 1. The relationship between particle radius and dimensionless parameters. The straight lines connects the points in which dimensionless parameters are 1.

FIG. 2. Snapshots of the slug. The number of particles in the unit cell is 72 and the volume fraction is 0.348. The ratio of the height to the width is 6. Time interval is 2.0 × 10^{−2} sec.

FIG. 3. Snapshots of the slug in three-dimensional simulation. The number of particles is 50 in the unit cell and the volume fraction is 0.45. The ratio of the height to the width is 4, and time interval between them is 2.0 × 10^{−2} [sec].

FIG. 4. A snapshot of the simulation with square cells with the periodic boundary condition. It is shown with four periodic images. The number of particles in the unit cell is 90 and the volume fraction is 0.327.

FIG. 5. The sequence of the simulation with fixed particles. The number of free particles and fixed particles in the unit cell are 128 and 10, respectively. The velocity of induced fluid is u^{∞} = 0.3U_{0}, where U_{0} = mg̃/6πμa is the one particle sedimentation velocity. The time interval from the left-up to the right-down is 0.112 sec. We can see three bubbles flowing up through the bed (except for the initial instabilities).

FIG. 6. The time evolution of standard deviation of particle velocities in the simulation of Fig.

5. Velocities are normalized by one particle sedimentation velocity U_{0}. 1 step means 10^{−4} sec.

FIG. 7. The time evolution of the volume fraction of particles in the area that is the horizontal region at 8/25 of the cell height above the bottom in the simulation of Fig. 5. 1 step means 10^{−4} seconds.

FIG. 8. The velocity distribution functions (VDFs) of (a) the horizontal direction and (b) the vertical direction in the simulation of Fig. 4 (with square unit cells). The velocities are normalized by one particle sedimentation velocity U_{0}. The vertical component of velocity is plotted in the coordinate system in which the average velocity of the fluid is zero and the direction of gravity is negative direction. The horizontal VDF can be approximated by exp(−20.0Û^{2}), where Û ≡ U/U_{0}. The left hand side and the right hand side of the vertical VDF can be approximated by exp(−10.0(Û + 0.2)^{2}) and exp(−4.5(Û + 0.2)^{2}), respectively. The VDF f(Û) is normalized by ∫ dÛ f(Û) = 1.

FIG. 9. The velocity distribution functions (VDFs) in the simulation of Fig. 5(with fixed particles). The velocities are normalized by one particle sedimentation velocity U_{0} and we use the coordinate system in which the fixed particle’s velocity become zero. The VDFs in horizontal direction and vertical direction can be approximated by exp(−5.0|Û|) and exp(−3.6|Û|), respectively. Here Û ≡ U/U_{0}. The normalization is the same as that in Fig. 8.

FIG. 10. The velocity distribution function (VDF) in (a) the active state and (b) the inactive state of the simulation with fixed particles which is the same of that in Fig. 9. The velocities are normalized by one particle sedimentation velocity U_{0} and we use the coordinate system in which the fixed particle’s velocity become zero. The Gaussian-fitted function is also ploted.

FIG. 11. The velocity distribution function (VDF) of (a) the horizontal direction and (b) the vertical direction in the simulation of Fig. 2 (with rectangular unit cells). The velocities are normalized by one particle sedimentation velocity U_{0}. The VDF in the horizontal direction can be approximated by exp(−25.0|Û|) in the exponential form or (1.0 + 700.0Û^{2})^{−1.5} in the t-distribution form. The VDF in the vertical direction for the right-hand-side by can be approximated by exp(−50.0|Û + 0.115|), in the exponential form or by (1.0 + 800.0(Û + 0.115)^{2})^{−2.0} in the t-distribution form. Here Û ≡ U/U_{0}. The normalization is the same as that in Fig. 8.

FIG. 12. The substracted velocity distribution function (VDF) from the VDF in Fig. 11 (b). We make this figure by substracting the uniform distribution, representing the contribution of falling particles in the slug, from the original VDF in the range of −0.4 < Û < −0.115. This VDF can be fitted for the left-hand-side by (1.0 + 700.0(Û + 0.115)^{2})^{−1.5} and for the right-hand-side by (1.0 + 800.0(Û + 0.115)^{2})^{−2.0}. Here Û ≡ U/U_{0}.

FIG. 13. Power spectrum E(ω) of (a) the holizontal velocity and (b) the vertical velocity in the simulation of Fig. 4 (with square unit cells). The frequency is normalized by ω_{0} = 2π/2048 [step^{−1}] ≃ 30.7[sec^{−1}]. The straight lines are least-square fits between the frequency from ω_{r} to ω_{p} with −1.597 ± 0.009 and −1.626 ± 0.009 slope for U_{x} and U_{y} respectively.

FIG. 14. Power spectrum E(ω) of (a) the holizontal velocity and (b) the vertical velocity in the simulation of Fig. 2 (with rectangular unit cells). The frequency is normalized by ω_{0} = 2π/2048 [step^{−1}] ≃ 30.7[sec^{−1}]. The straight lines are least-square fits between the frequency from ω_{r} to ω_{p} with −1.494 ± 0.01 and −1.513 ± 0.009 slope for U_{x} and U_{y} respectively.

FIG. 15. Power spectrum E(ω) of (a) the holizontal velocity and (b) the vertical velocity in the simulation of Fig. 5(with fixed particles). The frequency is normalized by ω_{0} = 2π/2048 [step^{−1}] ≃ 30.7[sec^{−1}]. The straight lines are least-square fits between the frequency from ω_{r} to ω_{p} with −1.514 ± 0.007 and −1.490 ± 0.007 slope for U_{x} and U_{y} respectively.

TABLES

TABLE I. Values of the dimensionless parameters and the characteristic time T_{r}[sec] and T_{p}[sec] when the particle radius is 10^{−3}[cm].

Air  Water
St  9.1 × 10^{0}  1.8 × 10^{−3}
P e  2.5 × 10^{5}  1.5 × 10^{5}
Re  2.0 × 10^{−2}  3.3 × 10^{−3}
T_{r}  3.1 × 10^{−3}  5.6 × 10^{−5}
T_{p}  3.3 × 10^{−4}  3.1 × 10^{−2}

TABLE II. Values of Re_{(p)} in our simulations, where L is measured by the maximum length of the unit cell.

square cell(Fig. 4)  rectangular cell(Fig. 2)  with fixed particles(Fig. 5)
Re_{(p)}  9.9 × 10^{2}  2.1 × 10^{3}  2.1 × 10^{3}
Re_{(p)}^{3/4}  1.8 × 10^{2}  3.1 × 10^{2}  3.1 × 10^{2}