Three-Dimensional VOF-DEM Simulation Study of Particle Fluidization Induced by Bubbling Flow

Abstract

The bubbling flow plays a key role in gas–liquid–solid fluidized beds. To understand the intrinsic fluidization behaviors at the discrete bubble and particle scale, coupled simulations with the volume of fluid model and the discrete element method are performed to investigate the effects of the gas inlet velocity, particle properties and two-orifice bubbling flow on particle fluidization. Three-dimensional simulations are carried out to accurately capture the dynamic changes in the bubble shape and trajectory. A bubbling flow with a closely packed bed is simulated to study the onset of particle fluidization. The obvious phenomena of particle fluidization are presented by both the experiment and simulation. Although an increasing gas inlet velocity promotes particle fluidization, the good fluidization of particles cannot be achieved solely by increasing the gas inlet velocity. When the channel is packed with more particles, the bubbles take a longer time to pass through the higher particle bed, and the bubbles grow larger in the bed. The increase in particle density also extends the time needed for the bubbles to escape from the bed, and it is more difficult to fluidize the particles with a larger density. Even if more particles are added into the channel, the percentage of suspended particles is not significantly changed. The percentage of suspended particles is not increased with a decrease in the particle diameter. The particle suspension is not significantly improved by the bubbling flow with two orifices, while the particle velocity is increased due to the more frequent bubble–particle collisions. The findings from this study will be beneficial in guiding the enhancement of particle fluidization in multiphase reactors.

1. Introduction

Gas–liquid–solid fluidized beds have found wide application over the years for many industrial reactions. They provide excellent heat and mass transfer performance owing to the bubble and particle dynamics. It has been recognized that the bubble dynamics plays a key role in the transport phenomena and fluidization behaviors [1,2,3,4]. Although an overwhelming amount of research has been dedicated to describing the complexity of the gas–liquid–solid hydrodynamics during the past few decades, the challenges in real industries call for new first-principles-based alternative solutions to more deeply understand the intrinsic fluidization behaviors. Due to the limited freedom in experimental measurements, computational fluid dynamics (CFD) has proven to be useful in describing the hydrodynamics in gas–liquid–solid fluidized beds [5].

Significant effort has been made to carry out CFD simulation studies of gas–liquid–solid fluidized beds. For engineering purposes, the Eulerian multi-fluid models are very popular due to the low computational effort required [6,7,8,9,10]. However, these models fall short of providing bubble-scale or particle-scale information. Alternatively, the discrete particle method (DPM) with the volume of fluid (VOF) model has been developed to resolve gas–liquid–solid flows. Fan’s group has conducted pioneering research on gas–liquid–solid fluidization systems. A combined CFD-VOF-DPM method was developed to simulate single bubble or bubble wake behaviors [11,12,13,14]. The topological change in the bubble surface was directly calculated with the VOF model, while the motion of particles was described by the DPM. Their simulations focused on bubble or bubble wake behaviors in suspended particle systems, while particle fluidization was not thoroughly discussed. Liu and Luo [15] used a combined CFD–VOF–DPM method to simulate the bubble rising and coalescence phenomenon in low-hold-up particle suspension systems. The gas and particle phases were simulated by the VOF and DPM methods, respectively. Only a few particles were distributed in the two-dimensional simulation domain. Particle entrainment by a single bubble or two bubbles was investigated in particle suspension systems. However, the bubble formation process was simplified by only patching a circular region in the simulation domain. Chen’s group also used the combined CFD–VOF–DPM method to investigate the wake structure and particle entrainment behavior of a single bubble in gas–liquid–solid systems [16,17]. Due to the very low solid hold-up, the particle–particle interactions were ignored in their studies. A small number of particles were sparsely distributed in the simulation domain, while the particles were closely packed in the experiment. This discrepancy can result in the inappropriate prediction of particle fluidization. Bo et al. [18] simulated the particle motion and entrainment process under the action of single and double bubbles using the computational fluid dynamics–deterministic particle track model–volume of fluid (CFD–DPTM–VOF) method. The influence of bubble motion on particle entrainment was described qualitatively and quantitatively. The particle–particle collisions were calculated with the direct Monte Carlo method. However, the particle collision model still needs to be verified in three-dimensional simulations. Although the success of the previous models was achieved with good agreement between the simulations and experiments, these two-dimensional simulations were still oversimplified for real multiphase flows with a three-dimensional nature. Furthermore, the sparse distribution of particles is not truly representative of the industrial situation in gas–liquid–solid fluidized beds. Liao et al. [19] coupled the volume of fluid (VOF) method with the multiphase particle in cell (MP-PIC) method to simulate bubble pair coalescence in a slurry. However, the particle’s penetration into the gas bubble was represented by the VOF method, which is not physically reasonable. Van Sint Annaland et al. [20] presented a hybrid model for the simulation of three-dimensional gas–liquid–solid flows using a combined front tracking (FT) and discrete particle (DP) approach. The hard sphere DP model was extended for the particles, and the FT model was combined for complex free surface flows. Although three-dimensional computational results were provided, revealing the particle entrainment behaviors, the coupled FT–DP simulations were mainly used for demonstration purposes and a qualitative discussion was given to illustrate the bubble-induced particle fluidization.

Compared with the CFD–VOF–DPM method, the VOF model coupled with the discrete element method (DEM) is more sophisticated and can reveal the flows in three-phase fluidized beds. The particle–particle collisions are described using the soft sphere model in DEM. Sun and Sakai [21] proposed a VOF–DEM coupled model to simulate gas–liquid–solid flows. This model was validated by the simple examples of water entry, a dam break, and a rotating tank. Pozzetti and Peters [22] developed a multiscale VOF–DEM model with a dual-grid approach to resolve the bulk scale from the fine fluid scale. Wu et al. [23] used a VOF–DEM model to study the turbulent free surface flow in a stirred mixing system. In the above VOF–DEM simulations, no bubble–particle interaction was taken into account, and only the free surface flow was simulated with particle physics. Ge et al. [24] investigated the fluid flow and particle dynamics in a Taylor flow microreactor, and the VOF–DEM model predicted the dilute particle flow well. Li et al. [25] simulated a single bubble’s behavior with a focus on bubble formation and the bubble rise velocity in gas–liquid–solid mini-fluidized beds through the VOF–DEM method. Moreover, the VOF–DEM method was used again to study the influence of the gas velocity and initial solid hold-up on the mini-fluidized beds with a bubbling flow [26]. In these mini-fluidized beds, small spherical bubbles rose along a rectilinear path in a very small channel, and the particle fluidization phenomenon was not obviously observed because of the small bubble rise velocity. Since bubble-induced particle fluidization has not been clearly elucidated, three-dimensional simulation studies are required to fully understand the underlying interaction between the particle fluidization and bubbling flow.

To reveal the influence of the bubble dynamics on particle fluidization, a continuous bubbling flow with particle motion is investigated with the three-dimensional VOF–DEM method, rather than the simple two-dimensional one. Instead of the sparsely distributed particles used in previous research, a closely packed particle bed is initially simulated to capture the onset of particle fluidization. To ensure consistency with the experiment, the dynamic process of bubble formation is also simulated. To evaluate the applicability of the VOF–DEM simulation method, an experiment is carried out to observe the bubble motion and particle fluidization. Deformed bubbles are generated with a single orifice or double orifices, and their motion trajectories in the channel are recorded by a high-speed camera. The effects of the gas inlet velocity and particle properties are numerically investigated with a focus on particle fluidization. The simulation also reveals the particle fluidization by the bubbling flow generated by two orifices.

2. Experimental Section

In order to directly observe the bubble and particle dynamics in a three-phase fluidized bed, an experimental facility is established, as shown in Figure 1. The experimental container is composed of a transparent acrylic channel with a cross-section of 8 mm × 8 mm and a height of 100 mm. The bubbles are generated by an orifice of 2.3 mm i.d. located at the bottom of the container. A gas flow is continuously supplied by a syringe pump, while liquid is added into the container at one time. By adjusting the motor speed of the pump, three gas flow rates of Q = 20, 40, and 80 mL/min are tested in the experiment. Initially, the height of the static liquid is maintained at 40 mm in each experiment. The motion of bubbles and particles is photographed and recorded via a high-speed camera (Revealer 5F04M, HF Agile Device Co., Ltd., Hefei, China) at a maximum frame rate of 1000 fps. An LED light source is used as the backlight. From the video, the average bubble rise velocity U_bubble is calculated by the bubble displacement at different time intervals:

$$U_{bubble}=\frac{z_2-z_1}{t_2-t_1}$$

(1)

where z is the position of the bubble center, and t is the time. The experimental error of the average bubble rise velocity is determined from the deviation from the mean of the five repetitions.

The particles are composed of polymethyl methacrylate. The particle density is 1200 kg/m³, and the particle diameter is almost uniform (d_p = 0.5 mm). Before the experiment, about 2500 particles are added to the channel and closely packed, corresponding to a bed height of 5 mm. After adding the particles and water in the channel, the experiment is started by supplying air through the orifice. Air (ρ_g = 1.2 kg/m³, μ_g = 1.8 × 10⁻⁵ Pa·s) and purified water (ρ_l = 1000 kg/m³, μL = 1 × 10⁻³ Pa·s, σ = 0.072 N/m) are used as the gas phase and liquid phase, respectively.

3. Mathematical Models

3.1. Volume of Fluid Model

The interface of the bubbles is captured based on the volume of fluid (VOF) method, which simulates a two-phase immiscible fluid by solving a set of continuity and momentum equations. Since the particles are initially closely packed, the effect of solid volume fraction α_p should not be ignored in the model. The continuity equation and momentum equation are given as follows:

$$\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f\right)+\nabla •\left({\alpha }_f{\rho }_f{\mathbf{u}}_f\right)=0$$

(2)

$$\frac{\partial }{\partial t}\left({\alpha }_f{\rho }_f{\mathbf{u}}_f\right)+\nabla •\left({\alpha }_f{\rho }_f{\mathbf{u}}_f{\mathbf{u}}_f\right)=-{\alpha }_f\nabla p+{\alpha }_f\nabla •\mathit{\tau }+{\alpha }_f\sigma \kappa \nabla {\alpha }_f+{\alpha }_f{\rho }_f\mathbf{g}$$

(3)

where ρ_f is the fluid density, kg/m³; u_f is the fluid velocity, m/s; p is the pressure, Pa; τ is the stress tensor, Pa; σ is the surface tension coefficient, N/m; κ is the mean curvature of the free interface, 1/m; α_f is the volume fraction of the fluid (gas or liquid), α_f = 1 − α_p. In the experiment, the bubbling flow is assumed to be a laminar flow due to the small gas inlet velocity, and τ is modeled by Newton’s law of viscosity:

$$\mathit{\tau }={\mu }_f\left[\nabla {\mathbf{u}}_f+{\left(\nabla {\mathbf{u}}_f\right)}^T\right]-\frac{2}{3}{\mu }_f\left(\nabla •{\mathbf{u}}_f\right)\mathbf{I}$$

(4)

where µ_f is the molecular viscosity, Pa·s; I is the unit tensor. In the continuum surface model [27], κ is calculated by the following equation:

$$\kappa =-\nabla •\left(\frac{\nabla {\alpha }_f}{\left|\nabla {\alpha }_f\right|}\right)$$

(5)

In the VOF model, different fluids are marked by the volume fraction, and the tracking of the phase interface is achieved by solving the continuous equation of each phase’s volume fraction in each computing unit:

$$\frac{\partial {\alpha }_f}{\partial t}+\nabla •\left({\alpha }_f{\mathbf{u}}_f\right)=0$$

(6)

where the value of α_f represents the fluid state in the computing unit:

$${\alpha }_f=\left\{\begin{array}{c}1\\ 0<{\alpha }_f<1\\ 0\end{array}\right.$$

(7)

The liquid is in the cells where α_f = 1, while the gas corresponds to α_f = 0. The interface is smeared over the cells where 0 < α_f < 1. The properties of the fluid, such as density ρ_f and viscosity µ_f, are calculated with the weighted average value of the volume fraction of two-phase fluids:

$${\rho }_f={\alpha }_f{\rho }_l+\left(1-{\alpha }_f\right){\rho }_g$$

(8)

$${\mu }_f={\alpha }_f{\mu }_l+\left(1-{\alpha }_f\right){\mu }_g$$

(9)

where l represents the liquid phase and g represents the gas phase.

3.2. Discrete Element Model

The translational and rotational motion of each individual particle are calculated with Newton’s second law [28]:

$$m_p\frac{d{\mathbf{u}}_i}{dt}=m_p\mathbf{g}+{\displaystyle \sum _j{\mathbf{F}}_{c,ij}}+{\mathbf{F}}_{pf}$$

(10)

$$I_p\frac{d{\mathbf{\omega }}_i}{dt}={\displaystyle \sum _j{\mathbf{T}}_{c,ij}}$$

(11)

where m_p is the particle mass; I_p is the rotational inertia; u_i is the translational velocity of the i-th particle; ω_i is the angular velocity; T_c;ij is the torque that is caused by the particle collision force; F_c,ij is the particle–particle or particle–wall contact forces, which are determined by the soft sphere model [29]; F_pf is the force generated by the fluid–particle interactions, which will be explained in the following section.

3.3. Fluid–Particle Interaction

In this study, the drag force, viscous force, pressure gradient force, and bubble interface force are taken into account in the fluid–particle interactions. F_pf is expressed as follows:

$${\mathbf{F}}_{pf}={\mathbf{F}}_d+{\mathbf{F}}_{\nabla p}+{\mathbf{F}}_{\nabla ·\tau }+{\mathbf{F}}_{interface}$$

(12)

The equation of drag force F_d is given by

$${\mathbf{F}}_d=-\frac{1}{8}\pi {\rho }_f{d}_p^2{C}_d\left|{\mathbf{u}}_p-{\mathbf{u}}_f\right|\left({\mathbf{u}}_p-{\mathbf{u}}_f\right)$$

(13)

where the drag coefficient C_d is calculated with the correlation proposed by Schiller and Naumann [30]:

$$C_d=\left\{\begin{array}{cc}\frac{24}{{\mathrm{Re}}_p}\left(1+0.15{\mathrm{Re}}_p^{0.687}\right)\hfill & {\mathrm{Re}}_p<1000\hfill \\ 0.44& {\mathrm{Re}}_p\ge 1000\hfill \end{array}\right.$$

(14)

$${\mathrm{Re}}_p=\frac{{\rho }_f\left|{\mathbf{u}}_f-{\mathbf{u}}_p\right|d_p}{{\mu }_f}$$

(15)

The forces from the pressure gradient and viscous stress are calculated as

$${\mathbf{F}}_{\nabla p}=-V_p\nabla p$$

(16)

$${\mathbf{F}}_{\nabla \cdot \tau }=-V_p\nabla ·\mathit{\tau }$$

(17)

where V_p is the particle volume.

The bubble surface region is marked with the local volume fraction α_f in a computational cell. To prevent particles from penetrating into the bubble, the bubble interface is modeled according to Sun and Sakai [21]:

$${\mathbf{F}}_{interface}=C{m}_p\nabla {\alpha }_f$$

(18)

where C is an empirical coefficient and it is set as 0.1 in this work.

4. Simulation Details

4.1. Geometry and Mesh

According to the geometry of the experimental container, the inside region of the channel is selected as the computational domain. Figure 2a shows the computational domain with a single orifice as a gas inlet. The cross-section of the domain with a single orifice is an 8 mm × 8 mm square. The diameter of the gas circular inlet is the same (2.3 mm i.d.) as that of the orifice in the experiment. Figure 2b shows the computational domain with two circular inlets, which is not experimentally used but is created for the simulation to investigate the effect of bubble formation from two orifices on particle fluidization. The unstructured meshes of these domains are generated with the blockMesh, surfaceFeatureExtract, and snappyHexMesh tools in OpenFOAM-5.x. First, a background mesh of hexahedral cells is created by the blockMesh tool. Second, the feature edges are extracted by the surfaceFeatureExtract tool for local mesh refinement. Finally, snappyHexMesh is executed to perform cell removal, cell splitting, and surface snapping. To sufficiently resolve the bubble surface boundary, the ratio of the bubble size to the mesh cell size is about 20. Figure 3 shows the meshes of the two computational domains, and the cell number is 245,760 and 491,583, respectively. The mesh size of the bottom surface is 0.125 mm, and the mesh inside the channel has a size of 0.25 mm.

4.2. Numerical Settings

The simulations in this study are carried out based on the VOF–DEM solver developed by Ge et al. [24] within the CFDEMcoupling framework [31], which combines the volume of fluid method of OpenFOAM with the discrete element method of LIGGGHTS. The interFoam solver in the open-source software OpenFOAM-5.x [32] is used to solve the governing equations of the gas–liquid flow, and the phase volume fraction transport equation is solved by the Multidimensional Universal Limiter with Explicit Solution (MULES) scheme. The temporal terms are discretized with the first-order Euler scheme. The convective term of the phase velocity is discretized with a second-order upwind scheme with limited gradients (linearUpwind) for numerical stability, while the convection term of the phase volume is discretized with the vanLeer scheme. The detailed numerical schemes are summarized in

Table 1. The numerical schemes used in the simulations.

Term	Numerical Scheme
$\frac{\partial \phi }{\partial t}$	Euler
$\nabla \phi$	Gauss linear
$\nabla •\left({\alpha }_f{\rho }_f{\mathbf{u}}_f{\mathbf{u}}_f\right)$	Gauss linear Upwind grad (U)
$\nabla •\left({\alpha }_f{\mathbf{u}}_f\right)$	Gauss vanLeer
$\nabla •\mathit{\tau }$	Gauss linear
${\nabla }^2\phi$	Gauss linear corrected

Note: φ is a variable; ∇² is the Laplacian term.

.

When the VOF solution is converged, the fluid–particle interactions are calculated and the open-source DEM solver LIGGGHTS [31] performs the computation of the particle motion with particle–particle or particle–wall collisions. For each VOF solution step, information exchange after 100 steps of DEM calculation is executed. The detailed settings of the VOF–DEM simulation are given in

Table 2. The simulation parameters.

Parameter	Value
Particle diameter (mm)	0.3, 0.5, 0.7
Particle density (kg/m3)	1200, 1500, 2000, 2500
Particle number	2000, 2500, 3000, 3500
Particle–particle normal stiffness (N/m)	100
Particle–particle tangential stiffness (N/m)	100
Restitution coefficient	0.9
Friction coefficient	0.3
DEM time step (s)	1 × 10−6
CFD time step (s)	5 × 10−5
Coupling interval	100

.

At the walls, the no-slip boundary condition is used for the fluid velocity. At the inlet, the fixed value boundary condition is specified for the gas inlet velocity. The other variables at the walls and inlet are specified with the zero gradient boundary condition.

5. Results and Discussion

5.1. Effect of Gas Inlet Velocity

Figure 4a presents the experimental and simulated bubble rising velocities at different gas inlet velocities (U_g = 0.08 m/s, 0.16 m/s, and 0.32 m/s). The error of the experimental results is about 5%. Compared with the experimental data, the simulated results of the bubble rising velocity are acceptable. As the gas inlet velocity increases from 0.08 m/s to 0.32 m/s, the bubble rising velocity also increases. The increase in bubble rising velocity is not particularly large. The rising bubble collides with the wall of the channel, which will reduce the rising velocity of the large bubbles at a high gas inlet velocity. For different gas inlet velocities, all bubbles rise in a zigzag trajectory, as shown in Figure 4b. During the rising process of the bubble, the bubble size also increases with the increase in the gas velocity, which can be seen from

Table 3. The simulated and experimental results of the bubble number and bubble diameter.

Ug (m/s)	Nbubble (exp.)	dbubble (exp.)	Nbubble (sim.)	dbubble (sim.)
0.08	3	4.19 mm	3	3.56 mm
0.16	5	4.95 mm	6	4.44 mm
0.32	6	5.71 mm	6	5.12 mm

. At the same time, the number of bubbles in the channel increases with the gas inlet velocity, which is due to the larger frequency of bubble formation at a higher gas inlet velocity. The number of bubbles is also influenced by the gas inlet velocity. When the gas inlet velocity is 0.08 m/s, three bubbles exist in the channel, while a gas inlet velocity of 0.32 m/s results in six bubbles. Moreover, at a large gas inlet velocity, the simulation reveals that bubble coalescence is pronounced, as shown in Figure 5. Based on the averaged bubble rising velocity and bubble diameter, the Re number of the bubbles is determined as 600, 750, and 878, and the Eo number is 1.7, 2.6, and 3.5, respectively, for the gas inlet velocities of 0.08 m/s, 0.16 m/s, and 0.32 m/s. It is clear that the bubbles near the orifice are ellipsoidal, and wobbling bubbles are also found due to the bubble–wall collisions. In the two-dimensional bubble simulations reported in the literature [12,16], it was found that the shapes of the bubbles were symmetrical, while asymmetrical shapes were found for the bubbles in the three-dimensional simulations of this study. Compared with the particles in the mini-fluidized beds studied by Li et al. [25], in this study, both the simulation and experiment showed the obvious fluidization of the particles caused by the bubbling flow.

In Figure 6a, it is obvious that an increasing gas inlet velocity will promote the fluidization of the particles. The particles near the orifice are mainly lifted by the elastic bubble surface. In the upper region of the channel, the bubble wake plays a role, as shown in Figure 7. Compared with the wake following a single bubble reported in the literature [16], the continuous bubbling flow generates significantly more wakes in the channel. For a larger gas inlet velocity, the bubble wake becomes stronger to entrain more particles with the bubbles. Although an increasing gas inlet velocity is conducive to particle fluidization, most of the particles still stay below 0.03 m, indicating that these particles are not sufficiently fluidized, as shown in Figure 5. The maximum fluidization height is about 30 mm. The particles that are not fluidized account for 83%, 81%, and 76%, respectively, for the gas inlet velocities of 0.08 m/s, 0.16 m/s, and 0.32 m/s. Figure 6b shows the distribution of the particle velocities influenced by the gas inlet velocity. Clearly, a larger gas inlet velocity causes more particles to be fluidized.

5.2. Effect of Particle Properties

The effect of the particle number on particle fluidization is investigated, as shown in Figure 8. For different particle numbers, the bubbles still rise in a zigzag trajectory, and the bubble size does not significantly change with the particle number.

Table 4. The effect of the particle number on the bubble number and diameter.

Particle Number	Nbubble	dbubble
2000	6	4.15
2500	6	4.12
3000	5	4.54
3500	5	4.51

illustrates the effect of the particle number on the bubble number and bubble diameter. When the particle number is increased to 3000, the number of bubbles existing in the channel is reduced, and the bubble size becomes slightly larger. This can be explained by the bed height of the particles that are not fluidized. As shown in Figure 9, as the number of particles is increased, the particle bed is thickened, which will result in a longer time required for the bubbles to pass through the particle bed. This equivalently extends the time for bubble formation. As a result, the bubbles should become much larger to escape from the particle bed.

Figure 10 shows the distributions of the particle positions and velocities for different numbers of particles. When more particles are added into the channel, the number of suspended particles is increased. Moreover, when the particle number is increased to 3500, the percentage of unsuspended particles is obviously reduced. This is attributed to the large bubbles in the particle bed. It is seen from Figure 10b that there is no significant difference in the particle velocity for different particle numbers. This can be explained by the more frequent collisions of suspended particles.

Figure 11 shows the distributions of the particle positions and velocities for different particle densities. When the particle density is larger, it is more difficult to fluidize the heavy particles. As can be seen from

Table 5. The effect of the particle density on the bubble number and diameter.

Particle Density (kg/m3)	Nbubble	dbubble
1200	6	3.52
1500	6	3.48
2000	5	3.55
2500	5	3.61

, with the increase in the particle density, the size of the bubble does not change significantly. As the particle density increases, the mass of a single particle increases, and the resistance force acting on the bubble also becomes stronger. Therefore, as shown in Figure 12, for the bed with particles having a large density, the bubbles experience a stronger resistance force from the particles and they take a longer time to pass through the particle bed. Similarly, the frequency of bubble formation becomes smaller, and therefore the number of bubbles existing in the channel is smaller. As a result, the bubble size becomes larger due to the extended time of bubble formation in the bed.

The distributions of the particle positions and velocities influenced by different particle densities are presented in Figure 13. As the particle density is increased, the percentage of suspended particles is decreased. Moreover, the particles are fluidized to reach a maximum height of about 30 mm. In Figure 13b, it is interestingly shown that when the particle density is increased, more particles have large velocities. The heavy particles can exert stronger gravitational forces on the elastic bubble surface. The stronger elastic deformation of the bubbles will result in a larger velocity in the particles.

Figure 14 shows the effects of different particle sizes (d_p = 0.3 mm, 0.5 mm, and 0.7 mm) on the bubble behavior and particle suspension. To keep the height of particle bed at about 5 mm, the number of particles in the beds is 13,888, 2500, and 1000, respectively, when adding the particles with the three diameters (d_p = 0.3 mm, 0.5 mm, and 0.7 mm). It is seen that the small particle of d_p = 0.3 mm can be easily fluidized by the bubbling flow. The fine particles are more susceptible to the action of bubbles or bubble wakes and are more likely to be suspended than the large particles. This finding is consistent with that in the study of Bo et al. [18]. However, there are still a large number of particles packed at the bed bottom. The percentage of suspended particles is not increased, which is also shown in Figure 15. It is found from

Table 6. The effect of the particle diameter on the bubble number and diameter.

Particle Diameter (mm)	Nbubble	dbubble
0.3	6	4.05
0.5	6	4.10
0.7	4	5.30

and Figure 14 that when the particle diameter is increased to 0.7 mm, there are very large bubbles in the channel. The large bubbles are formed due to the bubbles’ coalescence, which is indicated in Figure 16.

5.3. Bubbling Flow Generated by Two Orifices

As discussed above, the fluidization of particles is strongly dependent on the bubbling flow. Therefore, the bubbling flow generated by two orifices is simulated for a comparison with that generated by a single orifice. The cross-sectional area of the computational domain is two times larger than that with a single orifice, as shown in Figure 2. The particle beds have the same height as the beds in the above simulations. Figure 17 presents the comparison of these two scenarios for three different gas inlet velocities. From Figure 17, it is difficult to qualitatively judge the differences in particle fluidization. For the two-orifice bubbling flow, there are still many particles packed at the bottom. For the large channel, the increase in the bubble number does not significantly improve the particle fluidization. The quantitative results of the particle positions and velocities are provided in Figure 18. Although the improvement in particle fluidization is not significant, the use of two orifices indeed can move slightly more particles upward in the channel. At the gas inlet velocities of 0.16 m/s and 0.32 m/s, the percentage of particles with a large velocity is increased. This is because the increased bubble number results in a greater possibility of bubble–particle collision.

6. Conclusions

This study focuses on the three-dimensional simulation of particle fluidization induced by a continuous bubbling flow in a gas–liquid–solid fluidized bed. The simulations of the three-phase flow are carried out using the VOF–DEM coupling method. Instead of studying the bubble wakes, the bubble-induced particle fluidization is thoroughly discussed. By using the three-dimensional model, the dynamic changes in the bubble shape and trajectory are well captured. The bubbling flow with a closely packed particle bed is simulated to capture the onset of particle fluidization. By experimental validation, it is found that the three-dimensional VOF–DEM model can reliably simulate gas–liquid–solid flows. The effects of the gas inlet velocity, particle properties, and two-orifice bubbling flow on the bubble behavior and particle fluidization are investigated. The current study will be helpful to guide the process intensification of gas–liquid–solid fluidized bed reactors. The main conclusions are given as follows.

(1) The gas inlet velocity plays a role in controlling bubble formation and particle fluidization. By increasing the gas inlet velocity, the fluidization of the particles can be enhanced by the bubble–particle interaction and bubble wake. The continuous bubbling flow generates considerably more wakes in the channel. However, the method of solely increasing the gas inlet velocity cannot effectively realize the excellent fluidization of the particles.

(2) When more particles are added into the channel, the higher bed will result in a longer time required for the bubbles to pass through the bed. The increase in the particle density also extends the time required for the bubbles to leave the bed. It is more difficult to fluidize particles with a larger density. When adding more particles into the channel, the percentage of suspended particles is not significantly changed. For finer particles, the percentage of suspended particles is not increased.

(3) If the same channel region is used for fluidization, the bubbling flow generated by two orifices only causes slightly more particles to be suspended. There are still many particles packed at the bottom. The frequent bubble–particle collisions can contribute to the increase in particle velocity.

In a future study, the bubble–particle interaction should be further investigated with more sophisticated models to capture the attachment and detachment of the particles on the bubble surface. For example, the attachment of particles to bubbles plays a key role in mineral flotation processes.