A resolved Eulerian–Lagrangian numerical approach is used to study the heat transfer of 1204 heated spheres fluidized in a slit bed. This approach uses a direct numerical simulation combined with the immersed boundary method (DNS-IB). Pan et al. (2002, “Fluidization of 1204 Spheres: Simulation and Experiment,” J. Fluid Mech., **451**, pp. 169–192) studied the fluidization of 1204 spheres by a uniform flow without heat transfer using a fictitious domain-based DNS. The focus of this study is placed on the heat transfer between the heated spheres and fluid and also the fluidization by a jet flow. In the DNS-IB method, fluid velocity and temperature fields are obtained by solving the modified momentum and heat transfer equations, which result from the presence of heated spheres in the fluid. Particles are tracked individually and their velocities and positions are solved based on the equations of linear and angular motions; particle temperature is assumed to be a constant. The momentum and heat exchange between a particle and the surrounding fluid at its surface are resolved using the IB method with the direct forcing scheme. By exploring the rich data generated from the DNS-IB simulations, we are able to obtain statistically averaged fluid and particle velocity as well as particle heat transfer rate in a fluidized bed. Our results on the pressure drop and bed height are compared to the results of Pan et al. (2002, “Fluidization of 1204 Spheres: Simulation and Experiment,” J. Fluid Mech., **451**, pp. 169–192), which show good agreements. The case of the fluidization of 1204 spheres by a jet flow has also been studied and compared against the case of the fluidization by a uniform flow. The flow structures, drag, and heat transfer rate of two spheres placed along flow directions have been studied to understand the influence of a neighboring sphere. Results show that the trailing sphere has an insignificant effect on the leading sphere when it comes to the drag and heat transfer rate. On the contrary, the leading sphere can reduce the drag and heat transfer rate of the trailing sphere by more than 20% even when the two spheres are separated by six diameters. This demonstrates the need of a fully resolved DNS in accurately modeling dense particulate flows where a particle could be surrounded by multiple neighboring particles.

## Introduction

Because of enhanced solid–fluid heat and mass transfer and better mixing, fluidization of solid particles by fluid (gas or liquid) has been commonly used in many industrial applications, such as coal combustion, gasification, and fluid catalytic cracking. In designing these processes, numerical simulations can provide an efficient and accurate technique which is able to predict key operating parameters. For example, pressure drop, minimum fluidization velocity, solid fraction, and heat transfer rate in a fluidized bed could be obtained through computational fluid dynamics simulations. Simulation results can also be used to help better understand the flow dynamics and heat transfer of solid particles in fluidized beds. Compared to experimental studies, numerical simulations in general provide a low-cost and fast way to obtain many design parameters.

The existing numerical methods for modeling particulate flows can be divided primarily into three types: the resolved Eulerian–Lagrangian method or DNS, the unresolved Eulerian–Lagrangian method or discrete element method (DEM), and the Eulerian–Eulerian method or two-fluid model (TFM). Both the DNS and the DEM take the Lagrangian viewpoint by tracking particle dynamics while the carrying fluid is solved with Eulerian method. However, there is a significant difference between the DNS and the DEM in computing particle–fluid momentum and heat exchanges. For the DNS method, the flow fields, such as velocity, pressure, and temperature, are resolved around each particle based on the differential conservation of momentum and energy equations on a very fine grid; this allows the particle drag, lift, torque, and heat transfer rate to be computed directly from the resolved flow fields. The DEM, on the other hand, employs a much coarser grid, so the flow fields are not resolved; particle drag, lift, torque, and heat transfer rate have to be determined based on some empirical correlations. Compared to the DNS and the DEM, the TFM employs the Eulerian viewpoint for both solid particles and carrying fluid; it treats the large number of dispersed particles as another continuum “fluid” whose motion and temperature could be described by a set of partial differential equations derived from the laws of conservation.

Simulations of particulate flows require a multiscale approach. Large-scale industrial applications need the use of the TFM. For example, Zaman and Bergstrom recently used the TFM to study dilute gas–solid flow in pipes with rough walls [5]. However, the DEM results are needed to obtain constitution equations and transport properties of the particle phase, while the DNS results are needed to provide correlations for particle drag, lift, torque, and heat transfer rate. Furthermore, the DNS can also be used as a tool for validation and verification of the less detailed DEM and TFM.

Pan et al. [6] used a fictitious domain-based DNS method to study the fluidization of 1204 spheres in a slit bed by a uniform flow. Utilizing the large quantity data generated by the DNS, they were able to obtain a power-law correlation between the particle fluidization velocity and the fluid void fraction, as predicted by the well-known Richardson and Zaki [7] equation. The IB-based DNS or DNS-IB has become very popular in recent years for studying particulate flows [8,9]. This method represents the solid–fluid interface with a distributed forcing function [10]. The strength of the forcing function is chosen to enforce the no-slip velocity boundary condition; it makes the treatment of solid–fluid momentum exchange much easier. Feng and Michaelides extended DNS-IB to study particulate flows with heat and mass transfer by introducing a thermal density function to enforce no-temperature-jump conditions at the interface of particles and fluid [8,11–13]. Using the DNS-IB, Deen et al. [14] studied the heat transfer of particles in fixed and fluidized beds. These studies have provided an extensive validation of the DNS-IB method for solving particulate flows with heat transfer. The extension of fictitious domain-based DNS for particulate flows with heat transfer can be found in a recent article by Wachs [15], who studied particle motions due to strong free convection.

In the present paper, we extend the study by Feng and Musong [16] on heat and mass transfer in a fluidized bed. However, the focus of this study is as follows: to compare the hydrodynamic results from the present DNS-IB method with those obtained by fictitious domain-based DNS method [6] when a uniform flow is used for fluidization; to study the heat transfer between the heated spheres and fluid in a fluidized bed; to study how the bed properties and particle heat transfer rate change when a jet flow is used for fluidization; to study the drag and heat transfer rate of a pair of spheres when they are lined along flow direction but separated by different distances at Reynolds numbers 10, 100, and 400, respectively; and to study the influence on the drag and heat transfer rate of a sphere when a neighboring sphere is placed at various separation distances.

## Numerical Method

The DNS-IB has been implemented to solve particle motions and fluid velocity fields of particulate flows by Feng and Michaelides [8,9]. It has also been extended to solve heat transfer of moving particles with various applications [11–13]. Below, we provide a brief description of the method.

In the DNS-IB for particulate flows, two types of grids are used to solve the particle and fluid interactions: a fixed regular background grid or Eulerian grid for the entire computational domain occupied by both fluid and particles, and a moving grid or Lagrangian grid for outlining the immersed solid particles. A virtual boundary with a distribution force density function or an energy density function is used for the momentum or heat exchanges between the fluid and particle. The force density and energy density function are chosen in a way to enforce the no-slip velocity boundary condition and no-jump temperature boundary condition on the particle surface. The direct forcing scheme is used to determine the force and energy density [9,11,12]. Because the force or energy density function is defined as force or heat per unit volume, the actual boundary point should be placed beneath the surface at the centroid of the volume that the discrete force or energy density represents. This treatment has been implemented in the present study and it is able to yield more accurate results. For particles that have different temperatures than the fluid, the presence of heated particles will create a temperature gradient in the fluid that modifies the fluid's properties. For relatively small temperature differences between the particles and fluid, the Boussinesq approximation has been used successfully for coupling energy and momentum equations.

*D*suspended in a viscous flow of density $\rho f$ and viscosity $\mu f$. We chose the imposed fluidization velocity

*V*as the reference velocity, the sphere diameter

*D*as the reference length, and

*D*/

*V*as the reference time. The dimensionless temperature is defined as

For light particles considered in the present study, an explicit scheme [17] is used in updating particle motion. The particle temperature at its surface is assumed to be uniform at any time. The implication of this assumption requires the particle to have a very high thermal conductivity compared to that of the fluid, or a small Biot number, so that heat resistance inside the particle is less than the resistance between the particle and its surrounding fluid. A soft-sphere collision scheme [18] is used for particle–particle and particle–wall collisions to prevent unphysical particle overlaps.

## Drag and Heat Transfer Rate of a Pair of Heated Spheres

When the DEM is used in modeling particulate flows, the flow velocity and temperature fields are not resolved around a particle; instead the drag and heat transfer rate on a particle are modeled based on the correlations of a single particle immersed in an infinite flow. The effect of neighboring particles is taken into account by multiplying an influencing factor that is typically a function of solid fraction. In this section, we show that the drag and heat transfer rate of a spherical particle could be significantly changed in the presence of a neighboring particle.

We consider a uniform flow over a pair of spheres aligned along the flow direction and separated by a distance *δ*, ranging from 1 to 6 sphere diameters, at three Reynolds numbers of 10, 100, and 400, respectively. Results of drag coefficients and Nusselt numbers are then compared. The diameter of spheres is chosen to be *D* = 0.635 cm; fluid density is *ρ _{f}* = 1 g/cm

^{3}, its dynamic viscosity is

*μ*= 0.01 g/cm s, and Prandtl number Pr = 7. Periodic boundary conditions are applied at the two directions other than the flow direction. Other simulation parameters used for each case are listed in Table 1.

_{f}Re | 10 | 100 | 400 |
---|---|---|---|

Domain size (cm) | 4 × 4 × 12 | 3 × 3 × 12 | 3.5 × 3.5 × 14 |

Grid size | 80 × 80 × 240 | 100 × 100 × 300 | 150 × 150 × 600 |

Time step (s) | 0.001 | 0.001 | 0.0005 |

D/Δx | 12.7 | 21.2 | 27.2 |

Inlet velocity, V (cm/s) | 0.15748 | 1.5748 | 6.2992 |

Re | 10 | 100 | 400 |
---|---|---|---|

Domain size (cm) | 4 × 4 × 12 | 3 × 3 × 12 | 3.5 × 3.5 × 14 |

Grid size | 80 × 80 × 240 | 100 × 100 × 300 | 150 × 150 × 600 |

Time step (s) | 0.001 | 0.001 | 0.0005 |

D/Δx | 12.7 | 21.2 | 27.2 |

Inlet velocity, V (cm/s) | 0.15748 | 1.5748 | 6.2992 |

Table 2 shows the computed drag coefficients of a pair of spheres at different Reynolds numbers and different separation distances. The drag coefficient of an isolated sphere in an infinite uniform flow has been studied extensively both experimentally and numerically. There are many correlations between the drag coefficient and Reynolds number that can be found in the literature. Table 2 also lists the results from the numerical simulations by Feng and Michaelides [19] and from the correlation by Clift et al. [20]

Re = 10 | Re = 100 | Re = 400 | |||||
---|---|---|---|---|---|---|---|

Single sphere | Clift et al. [20] | 4.15 | 1.092 | 0.612 | |||

Feng and Michaelides [19] | 4.30 | 1.103 | 0.579 | ||||

Present | 4.50 | 1.124 | 0.711 | ||||

A pair of spheres at δ/D | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | |

1 | 3.81 | 2.23 | 1.061 | 0.247 | 0.734 | −0.064 | |

2 | 4.19 | 2.94 | 1.054 | 0.430 | 0.710 | 0.028 | |

3 | 4.28 | 3.34 | 1.081 | 0.548 | 0.690 | 0.190 | |

4 | 4.45 | 3.51 | 1.105 | 0.620 | 0.679 | 0.331 | |

5 | 4.47 | 3.71 | 1.117 | 0.671 | 0.702 | 0.382 | |

6 | 4.49 | 4.00 | 1.121 | 0.709 | 0.710 | 0.421 |

Re = 10 | Re = 100 | Re = 400 | |||||
---|---|---|---|---|---|---|---|

Single sphere | Clift et al. [20] | 4.15 | 1.092 | 0.612 | |||

Feng and Michaelides [19] | 4.30 | 1.103 | 0.579 | ||||

Present | 4.50 | 1.124 | 0.711 | ||||

A pair of spheres at δ/D | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | |

1 | 3.81 | 2.23 | 1.061 | 0.247 | 0.734 | −0.064 | |

2 | 4.19 | 2.94 | 1.054 | 0.430 | 0.710 | 0.028 | |

3 | 4.28 | 3.34 | 1.081 | 0.548 | 0.690 | 0.190 | |

4 | 4.45 | 3.51 | 1.105 | 0.620 | 0.679 | 0.331 | |

5 | 4.47 | 3.71 | 1.117 | 0.671 | 0.702 | 0.382 | |

6 | 4.49 | 4.00 | 1.121 | 0.709 | 0.710 | 0.421 |

It is seen that the present simulation results for a single sphere case agree well with those found in the literature at Re = 10 and 100, with an overall difference of a few percentages. However, more than 15% difference is observed at Re = 400. Two factors may contribute to the deviation. The first factor is that the current simulation domain is much smaller compared to those used in the literature; as a matter of fact the configuration of the current simulation domain yields a higher passing flow velocity and thus higher drag on the sphere. The second factor is that because of the limitation of available computational power, the grid used for Re = 400 is not sufficiently fine to fully resolve the laminar boundary layer near the surface of the sphere [16].

Table 2 clearly shows that the drag of the trailing sphere is significantly reduced compared to the leading sphere. In general, the influence on the drag of the leading sphere due to the presence of a trailing sphere is not significant, especially when the separation distance is more than 3*D*. However, the influence of the leading sphere onto the drag of the trailing sphere can be large, which also increases as flow Reynolds number increases even when the separation distance remains the same. At a separation distance of 2*D*, the drag on the trailing sphere is 70% at Re = 10, 41% at Re = 100, and 4% at Re = 400 compared to the drag on the leading sphere; at a separation distance of 5*D*, the drag on the trailing sphere is 83% at Re = 10, 60% at Re = 100, and 54% at Re = 400. At Re = 400, the drag to the leading sphere is largest when it is in contact with the trailing sphere; it gradually decreases as the separation distance increases; however, after the separation distance exceeds 4*D*, the drag starts to increase again. This nonmonotonic behavior is also observed in the case of two cylinders placed in tandem at moderate Reynolds numbers between 1 and 40 in a recent study by Vakil and Green [21].

Figure 1 shows the temporal behavior of drag coefficient at Re = 400 for a single sphere and a pair of spheres with a separation distance of 1*D*, 3*D*, and 6*D*. The accompanying vortex structures are shown in Fig. 2 by using the *λ*_{2}-method of vortex identification. There are a few interesting observations in these figures. In the case of a single sphere, periodical vortex loops or hairpins are constantly shed from the surface of the sphere at a period of 0.87 s, resulting in a Strouhal number St = 0.12. This agrees well with Wu and Faeth, who found Strouhal number around 0.13 for $335\u2264Re\u2264620$ [22]. Other researchers reported higher Strouhal numbers. In the case of two spheres placed along the flow direction, a slight increase in the drag coefficient on the leading sphere and a negative drag coefficient on the trailing sphere are observed when they are in contact (separation distance *δ* = *D*). The negative drag coefficient on the trailing sphere implies a “push” effect to the leading sphere which explains the higher drag coefficient of the leading sphere. It is also found that the trailing sphere is able to stabilize the flow when placed closely with a gap of one or two diameters; no hairpin vortex shedding is seen. However, as the gap is increased to three diameters or more, vortex shedding occurs again. At *δ* = 3*D*, the wake of the leading sphere has vortex hairpins beginning to form, but they are interrupted by the presence of the trailing sphere; the shedding of vortex hairpins is largely due to the leading sphere. However at *δ* = 6*D*, the wake of the leading sphere is able to be fully developed before it is interrupted by the trailing sphere; in addition to the Strouhal number of 0.12 due to the leading sphere vortex shedding, an additional higher Strouhal number of 0.5 is also observed, which is likely due to the flow instability of the trailing sphere.

Table 3 shows the Nusselt numbers for a single sphere and a pair of spheres at different separation distances. In the case of forced convection of an isolated sphere in a large flow domain, Ranz and Marshall developed the following correlation experimentally [23]:

Re = 10 | Re = 100 | Re = 400 | |||||
---|---|---|---|---|---|---|---|

Single sphere | Ranz and Marshall [23] | 5.61 | 13.40 | 24.8 | |||

Feng and Michaelides [24] | 5.98 | 14.54 | 27.4 | ||||

Present | 5.43 | 12.70 | 29.8 | ||||

A pair of spheres at δ/D | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | |

1 | 4.97 | 3.00 | 12.04 | 4.98 | 29.3 | 5.2 | |

2 | 5.31 | 3.52 | 12.34 | 6.26 | 29.3 | 9.9 | |

3 | 5.39 | 3.83 | 12.56 | 6.85 | 29.3 | 13.8 | |

4 | 5.41 | 4.04 | 12.63 | 7.22 | 29.4 | 17.4 | |

5 | 5.42 | 4.16 | 12.67 | 7.51 | 29.7 | 18.1 | |

6 | 5.42 | 4.30 | 12.69 | 7.75 | 29.8 | 18.6 |

Re = 10 | Re = 100 | Re = 400 | |||||
---|---|---|---|---|---|---|---|

Single sphere | Ranz and Marshall [23] | 5.61 | 13.40 | 24.8 | |||

Feng and Michaelides [24] | 5.98 | 14.54 | 27.4 | ||||

Present | 5.43 | 12.70 | 29.8 | ||||

A pair of spheres at δ/D | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | Sphere 1 | Sphere 2 | |

1 | 4.97 | 3.00 | 12.04 | 4.98 | 29.3 | 5.2 | |

2 | 5.31 | 3.52 | 12.34 | 6.26 | 29.3 | 9.9 | |

3 | 5.39 | 3.83 | 12.56 | 6.85 | 29.3 | 13.8 | |

4 | 5.41 | 4.04 | 12.63 | 7.22 | 29.4 | 17.4 | |

5 | 5.42 | 4.16 | 12.67 | 7.51 | 29.7 | 18.1 | |

6 | 5.42 | 4.30 | 12.69 | 7.75 | 29.8 | 18.6 |

Results in Table 3 show that the wake created by the leading sphere also has a significant effect on the heat transfer rate of the sphere behind. Even at *δ* = 6*D* away from the leading sphere, the Nusselt number only recovers 60% compared to the leading sphere. However, the trailing sphere has an insignificant effect on the leading sphere for heat transfer rate. It must be pointed out that compared to the placement of boundary nodes on the surface of spheres as implemented in our previous study [16], the use of the centroid of a control volume associated to the body force as the boundary node, typically near half of the grid spacing, is able to improve simulation results.

The temperature contours for the case of a single sphere and a pair of spheres with a separation distance of 1*D*, 2*D*, 3*D*, and 6*D* are plotted in Fig. 3. The vortex shedding for a single sphere and a pair of spheres at *δ* = 3*D* and 6*D* causes nonsymmetric temperature contours behind the spheres. However, at a separation distance of 1*D* and 2*D* flows remain stable and symmetric even after a long period of time.

In summary, the results presented in this section show that the drag and heat transfer rate of a sphere could be significantly affected by its neighboring spheres. A simple influencing factor that depends only on solid fraction may not be sufficient to account for such effect. In order to obtain accurate representations of drag and heat transfer rate in a complex flow filled with many neighboring particles, the fully resolved discrete method or DNS needs to be applied. In Sec. 4, we employ the DNS to study the fluidization of 1204 spheres in a slit fluidized bed.

## Fluidization of 1204 Spheres by a Uniform Flow

The DNS-IB method described in Sec. 2 is used in this section to study the particulate flows with heat transfer in a slit bed. The bed size is [20.3 cm, 0.6858 cm, 70.22 cm] and it contains 1204 spheres. The diameter of the spheres is 0.635 cm. This is identical to the setup used by Pan et al. [6]. We used a uniform grid of $357\xd713\xd71229$, which results in a total number of 5,703,789 nodes. The grid spacing is thus $\Delta x=0.0571\u2009cm$, and one diameter of the sphere is outlined by 11.11 grid elements. Pan et al. used a finite element-based fictitious domain approach in their study with a grid of $297\xd711\xd71025$ or a total of 3,348,675 nodes. It took 1660 hrs (69 days) to simulate 26 s of physical time on four R12000 processors in a SGI Origin 2000 computer cluster [6]. The current computation is performed on a cluster composed of Intel Xeon 5160 Quad Core 3 GHz processors without use of parallelization; it takes about 290 hrs (12 days) to simulate 26 s of physical time.

We placed 1024 spheres randomly in the lower portion of the bed. At the beginning, the flow velocity and particle velocities are zero. The physical and simulation parameters used in this study are listed in Table 4. Two cases with fluidization velocities 4 cm/s and 4.5 cm/s are considered. However, most discussions below are based on the results of *V* = 4 cm/s. A soft-sphere collision scheme [18] with a spring constant *k* = 2000 dyn/cm and damping coefficient $\eta =1\u2009dyn\u2009s/cm$ is employed in the simulations. Because of a relatively low solid fraction and fluidization velocity, the use of a different collision scheme and the selection of collision input constants are expected to have insignificant effects on the simulation results.

Number of particles, N_{p} | 1204 |

Particle diameter, d | 0.635 cm |

Particle density, $\rho p$ | 1.14 g/cm^{3} |

Fluid density, $\rho f$ | 1.00 g/cm^{3} |

Fluid dynamic viscosity, $\mu $ | 0.01 g/cm s |

Domain size, l_{x} × l_{y} × l_{z} | 20.3 cm × 0.6858 cm × 70.22 cm |

Grid size, N_{x} × N_{y} × N_{z} | 297 × 19 × 1025 |

Grid spacing, Δx | 0.0571 cm |

Time step, Δt | 0.001 s |

Prandtl number, Pr | 7.0 |

Initial bed height, H_{0} | 40 cm |

Fluidization velocity, V | 4 cm/s and 4.5 cm/s |

Number of particles, N_{p} | 1204 |

Particle diameter, d | 0.635 cm |

Particle density, $\rho p$ | 1.14 g/cm^{3} |

Fluid density, $\rho f$ | 1.00 g/cm^{3} |

Fluid dynamic viscosity, $\mu $ | 0.01 g/cm s |

Domain size, l_{x} × l_{y} × l_{z} | 20.3 cm × 0.6858 cm × 70.22 cm |

Grid size, N_{x} × N_{y} × N_{z} | 297 × 19 × 1025 |

Grid spacing, Δx | 0.0571 cm |

Time step, Δt | 0.001 s |

Prandtl number, Pr | 7.0 |

Initial bed height, H_{0} | 40 cm |

Fluidization velocity, V | 4 cm/s and 4.5 cm/s |

### Hydrodynamics of a Fluidized Bed.

The bed properties, such as the pressure drop and bed height, can be determined based on the resolved flow velocity fields and particle positions. As the simulation begins and fluid enters the bed from the bottom, the particles start to fall downward because of gravity. However, the flow exchanges momentum with the particles and exerts an upward force onto them. The particles are being fluidized and able to move within the flow. Figure 4 shows the snapshots of particle distributions at several different time instances. For better visualizations, particles are differentiated based on their initial positions. It is seen that the bed rises little in the first 5 s when compared to the next 5 s. This is because the drag on particles is less than the apparent weight of particles, causing particles to settle down at the beginning. After about 15 s, the bed height becomes stable; however, the particles continue to mix. The velocity contours of fluid and particle mixtures corresponding to these time instances are plotted in Fig. 5. It shows that fluid velocity is higher within the bed than the region above, while particle velocity is lower compared to fluid velocity inside the bed.

We further decompose the bed into a number of $Lz/d$ layers along the vertical direction. By counting the number of the spheres in each layer, we are able to determine the solid volume fraction along the vertical direction and track the bed height. Figure 6 shows the evolution of bed height from *t* = 0 to *t* = 25 s at *V* = 4 cm/s. It is seen that the solid fraction at front decreases as the bed expands. After 15 s, the bed height is steady. The solid fraction in the bed direction fluctuates around a constant mean value; this shows that the particles are nearly uniformly distributed along the bed direction.

We also computed the average *x* and *z* velocity components of the particles, *V _{px}* and

*V*. The results are plotted in Fig. 7. For the case of

_{pz}*V*= 4 cm/s, the vertical velocity component

*V*becomes very small after 15 s and eventually it approaches zero after 20 s. Furthermore, it takes longer for

_{pz}*V*to reach steady-state as the fluidization velocity increases from 4 cm/s to 4.5 cm/s.

_{pz}Results of the DNS are able to provide the total pressure drop over the bed or the total pressure gradient term $\u2212(dp/dz)$ on the left side of Eq. (10); the first term on the right side of Eq. (10) can also be calculated straightforwardly. Hence, the wall pressure gradient can be determined. We compute the total pressure gradient and the wall pressure gradient for the two cases of different fluidization velocities as shown in Table 5. The final bed height $H$ and the solid volume fraction $\varphi $ are also listed in the table. As expected, an increase in fluidization velocity causes an increase in the bed height and a decrease in the solid volume fraction. This causes fewer particle–wall collisions, which explains the decrease of the wall pressure gradient. Results from Pan et al. [6], who used a fictitious domain method and a repulsive sphere–sphere collision scheme, are also listed in the table. A good agreement is found in the case of *V* = 4.5 cm/s, with less than 5% difference for all the variables listed. In the case of *V* = 4 cm/s, the bed height we predicted is about 7% higher than Pan et al.'s predicted value; the wall pressure gradient is nearly 10% lower than what Pan et al. computed. We also noticed that our simulation results show a much more uniformly distributed bed compared to the one obtained by Pan et al. as seen in Fig. 8.

Figure 9 shows the solid volume fraction within the fluidized bed along the *x*-direction at *V* = 4 cm/s after the bed becomes steady. Using the technique described in our previous study [25], we further compute the time- and space-averaged particle velocity in the *z*-direction and fluid–particle composite velocity inside the bed. These results are plotted in Fig. 10. The averaged particle velocity fluctuates around zero, indicating a stable fluidized bed. In addition, particle velocities near the walls are positive, implying a upward particle slip velocity along the wall. The composite fluid velocity again shows a large upward velocity near the wall boundary. This is consistent to the findings by Pan et al. [6].

### Heat Transfer in a Fluidized Bed.

We further investigate the heat transfer between particles and fluid in the fluidized bed. The temperature of fluid and particles is made dimensionless following Eq. (2). Fluid initial temperature is zero and particle temperature is one. However, these heated particles retain a constant temperature, though this assumption can be relaxed to a time-dependent temperature [12]. At the surfaces other than the inlet and the outlet, we use a zero temperature gradient boundary condition.

The resolved DNS generates a huge amount of data, including the flow pressure, velocity, and temperature at any grid point and time step as well as the particle positions, velocities, and forces. Such detailed information could be used to extract statistical values of particle phase, such as the time-averaged or spatial-averaged particle velocity and particle heat transfer rate. Figure 11 shows the flow temperature and pressure contours in several *yz*-planes at *t* = 1 s for the fluidization velocity of 4 cm/s. Heat transfer rate from a particle to the fluid is directly proportional to temperature gradient at the interface of the particle, hence the particles near the inlet of the bed where fluid temperature is low in general have larger heat transfer rates compared to those particles above. The flow temperature fields at the vertical central plane of the bed recorded at several different time instants are shown in Fig. 12. Because of strong heat convection between cold fluid and heated particles, fluid temperature within the bed rises quickly as it moves past heated particles. This lowers the temperature difference between fluid and particles and decreases the heat transfer rate of particles as fluid moves toward the outlet. Particles near the inlet continue to have strong heat transfer with fluid. As demonstrated in Sec. 2, the Nusselt number of each sphere can be significantly different when two spheres are placed along the flow direction; the trailing sphere has a considerable drop in heat transfer rate compared to the leading sphere. An analogy may be applied to the case of particulate flows with heat transfer in a fluidized bed; the heat transfer rate of particles far away from the inlet of a fluidized bed is expected to be affected significantly by those particles near the inlet. However, these particles are in motion so the phenomena in a fluidized bed are much more complex.

*i*. The particle-averaged Nusselt number is a time-dependent function. Figure 13 shows the particle-averaged Nusselt number as a function of time for both

*V*= 4 cm/s and 4.5 cm/s. The average $Nu\xaf$ starts at a large number because of the large temperature gradients at the beginning. Then, it drops off quickly and starts to fluctuate around a steady mean value. It is also noted that $Nu\xaf$ becomes steady at $t\u22483s$ for

*V*= 4 cm/s, which is much shorter compared to the time needed ($t\u224820s$) for overall particle velocity to stabilize as shown in Fig. 7. This is due to the higher particle Peclet numbers compared to the particle Reynolds numbers since the fluid Prandtl number is 7. It is also found that as the fluidization velocity increases from 4 cm/s to 4.5 cm/s (a 12.5% increase), the average Nusselt number increases from 10.7 to 12.1 (a 13.1% increase).

## Fluidization of 1204 Spheres by a Jet Flow

Other than uniform flows, jet flows are also commonly used in particle fluidization process. In this section, we study how the flow characteristics and properties of the bed change when the uniform flow is replaced by a jet flow in the fluidization of 1204 spheres. The high velocity jet is placed at the center of the inlet; the jet orifice has a width of 2.03 cm and a depth of 0.6858 cm; the jet velocity is 40 cm/s, resulting in the same mass flow rate as that of fluidization by a uniform flow with *V* = 4 cm/s. A smaller time step $\Delta t$ = 0.0002 s is used in the simulation because of much higher maximum flow velocity in a fluidized bed with a jet flow. All the other physical properties and simulation parameters are chosen the same as those used in Sec. 4.

Figure 14 shows the distributions of 1204 spherical particles at five different time instants: 1 s, 3 s, 5 s, 7 s, and 10 s. For better visualizations, particles are colored with four different colors based on their initial locations. From these snapshots, it is clearly seen that the particles are well mixed during the fluidization; the bottom region where the particles are dense or so-called “dead zone” accounts for roughly a quarter of the entire bed. However, the particles inside the dead zone are not stationary; some of the particles originally trapped in the dead zone are being fluidized out while “fresh” particles are being entrained into the dead zone.

The time and space (*z*-direction) averaged particle velocity and fluid–particle composite velocity for a steady bed are computed and their results are shown in Fig. 15. The space-averaging in the *z*-direction is done by averaging all the vertical velocity components of each node or each particle within a vertical column [25]. As expected, the particles in the middle of the bed are being entrained into the jet stream and move upward; however, as the jet losses momentum the entrained particles start to fall downward near the left and right walls of the bed. Furthermore, the averaged particle velocity approaches zero while the averaged mixture velocity approaches 4 cm/s as the bed becomes stable.

Table 6 lists the bed properties and particle-averaged Nusselt numbers when 1204 spheres are fluidized by a uniform flow of 4 cm/s and a jet flow of 40 cm/s. There are a few interesting findings that are worth mentioning. Simulation results show that the bed is slightly higher when the uniform flow is used for fluidization. However, the wall pressure gradient which is an indication of wall friction is much smaller in the jet flow fluidized bed. This can be explained by realizing a much smaller velocity gradient of mixture in the jet flow fluidized bed as shown in Fig. 15 compared to a higher velocity gradient of mixture in the uniform flow fluidized bed as shown in Fig. 10. The pressure drop over the jet flow fluidized bed is also much lower than that of the uniform flow fluidized bed. The final observation is that the particle-averaged Nusselt number in the jet fluidized bed is significantly higher than that of the uniform flow fluidized bed. One plausible explanation is that the fluid warms up quickly as it moves upward in the case of uniform flow fluidized bed which reduces heat transfer rate of particles away from the inlet, as seen in Fig. 12; however, in the case of the jet flow fluidized bed, high jet stream is able to force cold fluid to penetrate deeper into the bed, resulting in higher heat transfer rate when it interacts with particles along its path, as shown in Fig. 14.

Uniform flow (V = 4 cm/s) | Jet flow (V = 40 cm/s) | |
---|---|---|

Bed height (cm) | 50.4 | 47.4 |

Solid fraction | 0.2303 | 0.2166 |

−dp/dz (dyn/cm^{2}) | 42.071 | 31.093 |

−dp/d_{w}z (dyn/cm^{2}) | 10.469 | 1.327 |

Particle-averaged Nusselt number | 10.7 | 16.8 |

Uniform flow (V = 4 cm/s) | Jet flow (V = 40 cm/s) | |
---|---|---|

Bed height (cm) | 50.4 | 47.4 |

Solid fraction | 0.2303 | 0.2166 |

−dp/dz (dyn/cm^{2}) | 42.071 | 31.093 |

−dp/d_{w}z (dyn/cm^{2}) | 10.469 | 1.327 |

Particle-averaged Nusselt number | 10.7 | 16.8 |

## Conclusions

The fluidization and heat transfer of 1204 heated spheres in a slit bed have been studied using a resolved Eulerian–Lagrangian method or DNS-IB method. We first investigated the effect on drag and heat transfer rate of a sphere when a neighboring sphere is placed along the flow direction. These two spheres are separated by a distance ranging from 1 diameter of the sphere to 6 diameter of the sphere and flow Reynolds numbers are chosen at 10, 100, and 400, respectively. Our results show that the trailing sphere has an insignificant influence on the drag and heat transfer rate of the leading sphere. On the contrary, the leading sphere can reduce the drag and heat transfer rate of the trailing sphere by more than 20% even when the two spheres are separated by 6 diameters. This demonstrates the need of a fully resolved DNS in accurately modeling particulate flows where a particle could be surrounded by multiple neighboring particles. We studied the fluidization of 1204 heated spheres in a slit bed at two different uniform fluidization velocities of 4 cm/s and 4.5 cm/s. The hydrodynamics properties of the bed, such as the bed height, the solid fraction, and the wall pressure gradient, agree reasonably well with the simulation results of Pan et al. [6]. The particle-averaged Nusselt number which is the sum of Nusselt numbers of all the particles divided by the number of the particles is computed to determine the total heat transfer rate between the particles and fluid. It is found that the particle-averaged Nusselt number increases by 13% as the fluidization velocity increases from 4 cm/s to 4.5 cm/s. The fluidization of 1204 heated spheres by a jet flow at a velocity of 40 cm/s has also been studied. Compared to the uniform flow fluidized bed, our results show a much smaller velocity gradient of mixture near a wall in the jet flow fluidized bed and hence a lower wall pressure gradient; our results also show a higher particle-averaged Nusselt number in the jet flow fluidized bed likely due to the high velocity jet being able to force cold fluid to penetrate deeper into the bed, resulting in more heat being transferred from particles to fluid.

## Acknowledgment

This research work was supported by grants from the U.S. Department of Energy (Award No. DE-FE0011453) and National Science Foundation of China (Award No. 31271002). This work also received computational support from the Computational Systems Biology Core funded by the National Institutes of Health National Institute on Minority Health and Health Disparities (G12MD007591), from the National Institutes of Health.