Hydrodynamic and Heat Transfer Study of a Fluidized Bed by Discrete Particle Simulations

: A numerical simulation study was carried out to study the combined thermal behavior and hydrodynamics of a pseudo-2D fluidized bed using a computational fluid dynamics–discrete element method (CFD-DEM). To mimic the effect of heterogeneous exothermic reactions, a constant heat source was implemented in the particle energy equation. The effects of superficial gas velocity, bed height and heat source distribution were analyzed with the aid of averaged volume fraction and temperature distributions and velocity profiles. It was found that both the gas superficial velocity and the bed aspect ratio have a profound influence on fluidization behavior and temperature distributions.


Introduction
Fluidized beds find widespread application in various industrial processes because of their favorable heat transfer characteristics. Many of the important applications of fluidized beds involve reactions with large heat effects. This may give rise to high rates of heat removal or production. Examples can be found in fluidized bed coal combustion, fluid catalytic cracking and polymerization processes for the production of polyethylene. For a proper understanding of the role of heat production and transport, detailed knowledge is required. This includes considering studies of the fluidized bed hydrodynamics and heat transfer (see, e.g., the work of Kunii and Levenspiel [1]; Basu and Nag [2]; Zhou et al. [3,4]). Most of the previous work was based on temperature probes placed inside, mounted or placed on the fluidized bed walls (see, e.g., Borodulya et al. [5]; Valenzuela and Glicksman [6]). Limtrakul et al. [7], Zhou et al. [8] and Patil et al. [9] previously used a computational fluid dynamics-discrete element method (CFD-DEM) to study the hydrodynamics combined with heat transfer in fluidized beds. However, only recently has CFD-DEM with heat transport been used to investigate bubble-to-emulsion heat transfer (see Patil et al. [9]).
In this paper, we use a CFD-DEM approach to study heat transport in fluidized beds. This approach combines a computational fluid dynamics (CFD) model to describe the flow of the gas phase with a discrete element method (DEM) to describe the particle motion (Deen et al. [10]; van der Hoef et al. [11]). We present CFD-DEM simulations of a pseudo-2D fluidized bed with spherical particles, through which we studied the influence of particle size and superficial gas velocity on fluidized bed thermal behavior. To mimic processes with heat production inside the particles due to a heterogeneous exothermic chemical reaction, a source term was added to the thermal energy equation of the particles. The simulations were carried out with an open-source package, OpenFOAM-CFDEM-LIGGGHTS. The objective of this work is to study the influence of superficial gas velocity, bed aspect ratio and heat production distribution on the hydrodynamics and heat transport in the bed.

Governing Equations
The CFD-DEM model was first introduced by Hoomans et al. [12]. In this work, we use the same model equations as those reported by Patil et al. [9] and Li et al. [13]. In the CFD-DEM, the gas phase hydrodynamics is described by the Navier-Stokes equations, while the motion of each particle is described by Newton's second law. The volume-averaged conservation equations for the gas phase mass and momentum are given by   where represents the source term for momentum originating from the particulate phase. In this work, the Beetstra [14] drag force correlation is used to calculate the interphase momentum transfer coefficient: The solid phase is considered to be discrete, and consequently, the modeling of this phase is based on the tracking of individual particles, accounting for particle-particle and/or particle-wall encounters. The motion of each spherical particle with mass and moment of inertia is described by Newton's equations for spheres: where is the particle position. The terms on the right-hand side of Equation (4) represent forces due to the pressure gradient, drag, gravity and inter-particle contact. For the contact force, a linear spring dashpot model (i.e., Hertz model) is used. is the torque, and is the angular velocity.
The thermal energy equation of the gas phase is given by where represents the source term, accounting for the gas-particle heat transfer, whereas is the effective thermal conductivity of the gas phase. The latter is related to the microscopic fluid thermal conductivity ( ) as follows: This equation was originally proposed by Syamal and Gidaspow [15]. The source term due to the gas-particle heat transfer is obtained by summing and distributing the contributions of all particles using a discrete Dirac delta-function D in a similar manner as is done for the drag force: where is the particle temperature, is the temperature of the gas phase at the Eulerian position , and is the fluid-particle heat transfer coefficient obtained from the empirical correlation by Gunn [ , Re , Pr where and are the particle Reynolds number and Prandtl number.
The thermal energy equation for each particle phase is given by where is the heat production rate due to exothermic intra-particle reaction.

Model Verification
In this section, two test cases-without and with heat source term, respectively-were simulated to verify the heat transfer implementation. The test system was taken as a fixed bed with initially cold particles (273 K) that is heated up by a hot gas flow (373 K), which passes through the bed uniformly from the bottom of the bed. The particles were uniformly arranged in a block structure and maintained stationary, which is a similar geometry to that used by Van Sint Annaland et al. [17], Patil et al. [9] and Li et al. [13]. The physical properties of the particles and the gas in the CFD-DEM verification simulations are listed in Table 1. The geometrical characteristics of the bed are summarized in Table 2. Table 1. Particle and gas phase properties in the fixed bed verification.

Property (Unit), Symbol
Value Gas density (kg/m 3 ), ρg 1000 Inlet gas temperature (K), Tg,0 373 Gas viscosity (Pa s), µg 1.0 × 10 −3 Gas heat capacity (J/kg/K), Cpg 4187 Gas thermal conductivity (W/m/K), kg 0.5 Particle diameter (m), dp 3.95 × 10 −3 Particle density (kg/m 3 ), ρp 8400 Initial temperature of gas and particles (K) 273 Particle heat capacity (J/kg/K), Cpp 385 The simulation results were verified by comparing them with results from a semi-analytical model. The latter is based on a one-dimensional convection equation with interphase heat transfer between the gas and solid phase. For the 1D model, the gas and particle phase energy balances are given by where as is the specific interfacial area, given by as = 6(1 − εg)/dp. The details of the analytical solution of Equations (12) and (13) can be found in Bird et al. [18].
Dimensionless axial temperature profiles of both phases without heat production are shown in Figure 1. The simulation results are compared with the semi-analytical data from the 1D model. Figure 2 shows simulation results in the case with heat production.
It was found that the results of both models (both without and with heat production) agree very well. Hence, we conclude that the heat transfer model is correctly implemented in the 3D CFD-DEM model and can be used for heat transfer studies of fluidized beds.

Simulation Settings
A series of fluidization simulations were performed using the OpenFOAM-CFDEM-LIGGGHTS package (we refer to this package as CFDEM hereafter) of a pseudo-2D fluidized bed, which is schematically presented in Figure 3. Details of the geometry are shown in Table 3. The simulations were carried out with particles of 0.995 mm in diameter, which are classified as Geldart D type particles. Simulations were performed for initial bed heights of 0.04 m (aspect ratio 0.5), 0.08 m (aspect ratio 1) and 0.16 m (aspect ratio 2). The number of particles in the system amounted to 37,225, 74,450 and 148,900, respectively. Free-slip boundary conditions were applied to the front and back walls, and no-slip boundary conditions were used for the left and right walls. Table 4 lists the physical properties of the gas and particles that were used in the simulations. Simulations were performed for several gas superficial velocities, namely 0.4, 0.5 and 0.6 m/s (umf = 0.262 m/s). The particles were initialized at a temperature of 345 K and stacked on a cubic lattice at the bottom of the bed, after which a constant gas stream with Tg,0 = 330 K was supplied through the gas supply at the bottom of the bed. Each particle had a constant heat source of 0.0209 W/(K*m), which is the same as that used by Li et al. [13]. All simulations were performed for 10 s.   Figure 4 shows snapshots of the typical volume fraction distribution under conditions for each of the three simulated aspect ratios and superficial gas velocities. From this figure, it is apparent that bubbles grow in size as the gas superficial velocity is increased. From animations, it was observed that bubbles move more chaotically when the superficial gas velocity increases. This effect becomes more apparent with increasing H0/W.  To quantify the effect of the gas superficial velocity and initial bed height on bubble behavior, the temporally and spatially averaged gas porosity was obtained from the last 8 s of simulation data. It was found that 8 s of simulation time was sufficient to obtain well-converged results. In Figure 5, it is shown that the overall mean bed voidage increases more or less proportionally with increasing u0. When the superficial gas velocity increases, the differences of averaged voidage between different aspect ratios increase, which is in line with the observations from Figure 4. To quantify the effect of ug,0 and H0/W on the solids motion, the temporally and spatially averaged particle velocity was obtained from the last 8 s of simulation data. The time-averaged vertical particle velocities at z = H0, as a function of superficial gas velocity and aspect ratio, are depicted in Figure 6. It can be seen that particle velocity is affected by superficial gas velocity as follows: both the up-flow in the core and the down-flow in the annulus are enhanced. This can be attributed to enhanced bubble action. The bed hydrodynamics is also influenced by the aspect ratio, as indicated by the altered shape of the time-averaged vertical velocity. When the initial bed height increases, the up-flow in the core region becomes larger. This reveals that a higher aspect ratio promotes intensified upward solids motion.

Effect of Superficial Velocity and Aspect Ratio on Heat Transfer Behavior
In this section, we focus on the effect of the gas superficial velocity and bed aspect ratio on the gas-particle heat transfer behavior. Figure 7 shows snapshots of the particle temperature at different gas superficial velocities and aspect ratios. With increasing superficial velocity, it is obvious that the particle temperature is more homogeneous. This can be explained by the fact that the particles become better-mixed and the heat transfer rate between gas and particles is more intense. The whole bed temperature increases with increasing aspect ratio, which can be explained by the fact that with increasing H0/W, the bed mass increases, and hence, the associated overall heat source in the bed also increases. To characterize the bed thermal behavior, the temporal-and spatial-averaged particle temperature was calculated, along with the standard deviation of the particle temperature. Figure 8 shows the standard deviation (defined in Equation (14)) of the particle temperature as a function of the superficial velocity and bed aspect ratio. As can be seen, when the gas superficial velocity increases, the temperature standard deviation decreases in all cases. This can be explained by the fact that at higher u0, the bed is better-mixed and hence becomes more isothermal. Moreover, at higher aspect ratios, the relative influence of entrance effects is reduced and particles become better-mixed, which further reduces the standard deviation of the particle temperature. To quantify the effect of ug,0 and H0/W on the thermal behavior, the probability density function (PDF) of the temperature was calculated. Prior to this, all particle temperatures were nondimensionalized by subtracting the inlet gas temperature (Tp-Tg,0) and dividing by the driving force, i.e., the particle melting point (Tp,m = 380 K) minus the inlet gas temperature (Tp,m-Tg,0). Dimensionless temperatures of 0 and 1 correspond to a particle temperature equal to the inlet gas temperature and the melting point of the solid material, respectively. Subsequently, the dimensionless temperatures were used as an input to create a probability density function (PDF). This was done by dividing the dimensionless temperature range between 0 and 1 into 25 equal-sized bins. During a simulation period of 0.25-2 s, the number of particles in each of the temperature bins was counted (Li et al., [13]). When ug,0 was increased, the particle temperature PDF changed considerably (see Figure 9). The mean value of the temperature dropped as ug,0 was increased, implying that more heat is removed as more gas is supplied. When the aspect ratio increased, the mean value of dimensionless particle temperature increased, and the dimensionless temperature exceeded 1 at H0/W = 2.0. This can be explained by the fact that at higher initial bed height, there are more particles and, consequently, a stronger heat source in the bed, which will further increase the equilibrium temperature and dimensionless particle temperature. Figure 9. Dimensionless particle temperature probability density function (PDF) for varying aspect ratios: 0.5, 1 and 2 (top to bottom) and different gas superficial velocities.

Effect of the Heat Source Term on Heat Transfer Behavior
In the previous sections, it is assumed that all particles have uniform heat production. Under realistic conditions, however, heat production would depend on the exact amount of catalyst material in each particle, the degree of deactivation and differences in mass transfer limitations. To investigate the effect of differences in heat production, we defined a case with an extreme distribution in heat production properties: 98% "inactive" particles were set to have no heat production, whereas the remaining 2% "active" particles were set to have a heat production 50 times higher than a particle in a bed with uniform heat production. In this way, both the uniform and non-uniform bed will have the same total heat production.
As can be seen from the snapshots in Figure 10, all particles mix well in the bed. Because of the very high heat transfer coefficient, the heat transfer from the active particles to the gas and from the gas to the inactive particles is very efficient.
The mean temperature reaches an equilibrium state after 10 s. The steady-state bed-averaged gas temperature can be derived analytically, assuming ideal mixing of both gas and particles. The thus obtained equilibrium temperatures of the gas phase and the inactive (1) and active (2) particles are given by where X is the percentage of "inactive" particles. The terms on the right hand represent the gas inlet temperature, the adiabatic temperature rise and the temperature difference between gas and particles, respectively. Because of the very effective heat transfer, the relatively small temperature difference between "active" and "inactive" particles is an important factor because, in reality, many types of particles with different activities would be present in the bed simultaneously. The equilibrium temperature of the gas and the two types of particles with different heat production under different conditions are listed in Table 5. Table 5. Equilibrium temperature of the gas phase (Tg) and the active/inactive particles (Tp1/Tp2). From the snapshots of gas and particle temperature in Figure 10, the gas temperature and the distribution of particles with and without the heat source term, respectively, are shown. Beyond the entrance region, the gas temperature shown in Figure 10 is practically homogeneous and approximately equal to the temperature of the inactive particles, as shown in Figure 11. This aligns well with the equilibrium temperatures given in Table 5. Furthermore, Figure 11 shows the active particles distributed through the bed uniformly. As expected from the equilibrium temperatures, there is an obvious temperature difference between inactive and active particles. Figure 12 shows the probability density function of the dimensionless particle temperatures, which was created during a period of 10-20 s (in a similar way to that described by Li et al. [13]). The equilibrium temperatures of the (in)active particles are included for reference.

Non-Uniform Cases Uniform Cases
By comparing the temperature PDFs in Figure 9 (uniform heating) and the left column of Figure  12 (non-uniform heating), we can see that, in both cases, the temperature distribution is more uniform at lower superficial velocity. In the case of the non-uniform heating, a small extra peak appears associated with the small number of (very) active particles.
To inspect the behavior of the active particles, separate PDFs were obtained for both the inactive and active particles only (see the middle and right column of Figure 12, respectively). With lower superficial velocity and lower H0/W, inlet effects are dominant, leading to long tails towards low temperatures for both particle classes. From the comparison between simulation results and equilibrium temperatures, we observe that there still are differences, and the differences become larger at lower superficial velocity. This is because the assumption of a thermally ideally mixed bed is not valid. This also illustrates the added value of fluidized bed CFD-DEM simulations for the inspection of hot and cold zones in the bed.

Conclusions
In this work, CFD-DEM simulations of a non-isothermal gas-fluidized bed were presented. A heat source was incorporated in the particle phase thermal energy equation to mimic the heat liberation due to chemical conversion.
The influence of the gas superficial velocity and bed aspect ratio on the hydrodynamics and thermal behavior of a pseudo-2D fluidized bed were analyzed. With increasing superficial velocity and initial bed height, the overall mean bed voidage increases more or less linearly. The standard deviation of particle temperature provides quantitative information on the particle temperature distribution. At higher aspect ratio and with increasing ug,0, the particles become better-mixed, reducing the standard deviation. Furthermore, the time-averaged particle velocity profiles strongly depend on ug,0 and H0/W. Enhanced particle circulation with increasing superficial velocity is observed, whereas more particles are transported upwards in the central region for increasing initial bed heights. The fluidized bed is more isothermal when the superficial gas velocity and aspect ratio increase. The diversity of particles with and without heat production also has an apparent effect on heat transfer behavior. Because of the very good heat transfer characteristics, the heat transfer from very active particles via the gas phase to the non-active particles is very effective. The gas and particle temperatures predicted by CFDEM simulations agree reasonably well with the equilibrium temperatures.