**Computational Fluid Dynamics Simulation of Gas–Solid Hydrodynamics in a Bubbling Fluidized-Bed Reactor: E**ff**ects of Air Distributor, Viscous and Drag Models**

**Ramin Khezri 1, Wan Azlina Wan Ab Karim Ghani 1,2,\*, Salman Masoudi Soltani 3,\*, Dayang Radiah Awang Biak 1, Robiah Yunus 1, Kiman Silas 1, Muhammad Shahbaz <sup>4</sup> and Shiva Rezaei Motlagh <sup>1</sup>**


Received: 8 July 2019; Accepted: 7 August 2019; Published: 8 August 2019

**Abstract:** In this work, we employed a computational fluid dynamics (CFD)-based model with a Eulerian multiphase approach to simulate the fluidization hydrodynamics in biomass gasification processes. Air was used as the gasifying/fluidizing agent and entered the gasifier at the bottom which subsequently fluidized the solid particles inside the reactor column. The momentum exchange related to the gas-phase was simulated by considering various viscous models (i.e., laminar and turbulence models of the re-normalisation group (RNG), k-ε and k-ω). The pressure drop gradient obtained by employing each viscous model was plotted for different superficial velocities and compared with the experimental data for validation. The turbulent model of RNG k-ε was found to best represent the actual process. We also studied the effect of air distributor plates with different pore diameters (2, 3 and 5 mm) on the momentum of the fluidizing fluid. The plate with 3-mm pores showed larger turbulent viscosities above the surface. The effects of drag models (Syamlal–O'Brien, Gidaspow and energy minimum multi-scale method (EMMS) on the bed's pressure drop as well as on the volume fractions of the solid particles were investigated. The Syamlal–O'Brien model was found to forecast bed pressure drops most consistently, with the pressure drops recorded throughout the experimental process. The formation of bubbles and their motion along the gasifier height in the presence of the turbulent flow was seen to follow a different pattern from with the laminar flow.

**Keywords:** gasification; fluidized bed; CFD; hydrodynamics; multiphase flow

#### **1. Introduction**

Owing to the renewable, sustainable and environmentally-friendly nature of biomass-based energy generation, over the past few years, it has attracted a great deal of interest within both the scientific community and the industrial sector. Biomass fuel can dampen the dependency on the energy production from fossil fuels, and hence, mitigate the associated adverse environmental emissions [1–4]. Biomass gasification is still facing various technical challenges, repelling potential commercialization. Such challenges comprise the complex chemical reactions propagation inside the fluidizing zone and the restrictions imposed through process control and difficulty in the prediction of momentum patterns caused by the turbulent behaviour of the gas flow and moving particles [5]. The complicated processes—inherently associated with embedded limitations—would make it almost economically and environmentally unviable to experimentally investigate all the factors in play, and therefore, necessitate the development of reliable mathematical models.

Reliable models can facilitate the investigation of thermochemical processes associated with biomass gasification; likewise, they help to study and evaluate the effects of the operational parameters on the gasifier performance and product quality [6]. Numerical models and simulation studies for biomass gasification using either Eulerian [7–11] or Lagrangian (discrete element) [12–14] models have been the focus of many researchers in recent years. A large number of numerical modelling and simulation studies have been reported in literature on the complex behaviour of gas–solid mixing, as well as on the fluidization hydrodynamics, in order to formulate the problems and to devise solutions. Furthermore, using computational fluid dynamics (CFD), model-based turbulent flow calculations can be done more efficiently and in a quicker manner. Standard (STD) k-ε, re-normalisation group (RNG) k-ε and k-ω are among the most commonly-used models used in systems involving fluid flows with large Reynolds numbers; e.g., flows in many fluidized beds. The RNG k-ε model is a modified STD k-ε model and in many cases, produces improved results.

In a previous study, the hydrodynamics of gas–solid mixing in fluidized bed reactors was evaluated through numerical simulation and experimental techniques in order to investigate the effect of drag models on the hydrodynamics of fluidization [15]. In another study, a CFD model was developed to determine the fluid-particle interaction during fast pyrolysis inside a fluidized bed reactor with a 150 gh−<sup>1</sup> capacity [16]. The model was capable of calculating the heat, mass and momentum transport between the fluidization agent and the bed material (i.e., silica sand). However, the characterization of momentum in both studies was merely limited by a laminar model. In another study a three dimensional CFD model was developed to simulate a high-flux circulating, fluidized bed riser [17]. The effects of different drag models on the hydrodynamics of the gas–solid fluidization were investigated. The results showed that the drag models of Gidaspow and Syamlal–O'Brien were able to accurately predict the solid volume fraction values, at almost every region inside the fluidized bed except for the zone near the bottom surface with the high concentration of the solids. In another study, a CFD simulation of flow through the gas distributor in a multi-stage biomass gasifier was conducted by adopting the Eulerian-multiphase model approach, where the effects of different bed materials, and a distributor perforated area, on the hydrodynamics of a solid–gas fluidized bed were examined [18]. Although the most critical aspects affecting the fluidized bed hydrodynamics were carefully addressed in the study, it lacked a study to reveal the impacts of different drag models.

This study aimed to simulate the hydrodynamics of gas–solid fluidization. The simulation results were experimentally validated on a bubbling fluidized bed gasifier. The effects of different momentum characteristics (i.e., laminar and turbulent models) on the hydrodynamics of gas–solid mixing, as well as on the fluidization of the dispersed solid phase, were studied. In addition, the pressure drop profile obtained from different drag models, including laminar, RNG k-ε and k-ω turbulent were compared with the experimentally-derived data. Furthermore, the effect of air distributors with different pore sizes was studied numerically to determine the turbulent viscosity of the fluid passing through the pores of the distributor plate. Ultimately, the motion patterns of the particles and bubble formation linked to the hydrodynamics of granular solid fluidization were investigated. The developed hydrodynamic model, upon successful experimental validation, is able to effectively pinpoint the optimum operating conditions for the investigated parameters, including the superficial velocity, bed expansion, pressure drop gradient throughout the reactor interior, and the momentum exchanged between the particles and the media. The above parameters are not easy to measure experimentally, yet are significant enough to be considered for accurate reactor optimization.

#### **2. Modelling and Experimental Methods**

#### *2.1. Governing Equations*

The basic equations used to in the development of the numerical models will be introduced and explained in this section. The numerical formulae are the governing equations used by Ansys Fluent Solver to calculate the corresponding parameters.

#### 2.1.1. Gas–Particle Flow Equations

An unsteady-state Eulerian–Eulerian multiphase flow model was used to compute the transient nature of gas–solid in the fluidization of granular solids. The different phases in the Eulerian model were treated mathematically as interpenetrating continua. The phasic volume fraction was considered to differentiate the volume occupied by each phase as a function of time and space, implying that although the involved phases were interacting with each other, the volume of each phase could not be occupied by another phase [19]. Conservation equations were derived for all phases and closed through constitutive relations defined via the application of kinetic theory. The gas phase was modelled for laminar, RNG k-ε turbulent and k-ω turbulent characteristics, while the particle phase was modelled using the kinetic theory of granular flow. The governing conservation equations are presented in Table 1.

**Table 1.** Governing equations of momentum and continuity for gas–solid phases [20–25].


The volume fraction, density and instantaneous velocity of gas/solid phase in the governing equations in Table 1 are represented by δ, ρ and u, respectively. Subscripts g and s denote gas and solid phases, respectively.

#### 2.1.2. Gas–Solid Interaction

The momentum exchange coefficient of solid–gas, *Ksg*, is defined as [18]:

$$\mathcal{K}\_{\mathfrak{s}\mathfrak{g}} = \frac{\mathfrak{e}\_{\mathfrak{s}} \mathfrak{p}\_{\mathfrak{s}} f}{\mathfrak{r}\_{\mathfrak{s}}} \tag{1}$$

where the particulate relaxation time, τ*s*, is defined as:

$$
\tau\_{\mathfrak{s}} = \frac{\rho\_{\mathfrak{s}} d\_{\mathfrak{s}}^2}{18\mu\_{\mathfrak{s}}} \left( d\_{\mathfrak{s}} \text{ is the particle size in the solid phase.} \right) \tag{2}
$$

#### 2.1.3. Drag Models

One of the significant terms in describing the hydrodynamics of fluidized bed in a Eulerian–Eulerian approach is the drag force. The significance of drag force is due to its effect on the characteristics of the interaction between the gas and solid phases. The drag models employed in this study to calculate the gas–solid momentum exchange coefficient are the Syamlal–O'brien, Gidaspow and energy minimum multi-scale method (EMMS). The Gidaspow model is a combination of the Wen and Yu model, and the Ergun equation [13,24], and is well-suited for densely packed fluidized bed applications:

$$K\_{\mathcal{S}\mathcal{F}} = \begin{cases} \frac{3}{4} \mathcal{C}\_{D} \Big( \frac{\delta\_{s} \delta\_{\mathcal{S}} \rho\_{\mathcal{S}} \left| \overleftrightarrow{v}\_{s} - \overrightarrow{v}\_{s} \right|}{d\_{s}} \Big) \delta\_{\mathcal{S}}^{-2.65}, & \delta\_{\mathcal{S}} > 0.8\\ 150 \Big( \frac{\delta\_{s} \left( 1 - \delta\_{\mathcal{S}} \right) \mu\_{\mathcal{S}}}{\delta\_{\mathcal{S}} d\_{s}^{2}} \Big) + 1.75 \Big( \frac{\delta\_{s} \rho\_{\mathcal{S}} \left| \overleftrightarrow{v}\_{s} - \overrightarrow{v}\_{s} \right|}{d\_{s}} \Big) & \delta\_{\mathcal{S}} < 0.8 \end{cases} \tag{3}$$

where, *CD*, is the drag coefficient and defined for a smooth particle as [26]:

$$\mathbf{C}\_{D} = \begin{cases} \frac{24}{\delta\_{\ $} Re\_{s}} \left[ 1 + 0.15 \left( \delta\_{\$ } Re\_{s} \right)^{0.687} \right], & Re < 1000\\ 0.44, & Re \ge 1000 \end{cases} \tag{4}$$

where, *Res*, is a function of slip velocity and is described as:

$$Re\_s = \frac{\rho\_{\mathcal{S}} d\_s \left| \overrightarrow{\upsilon\_s} - \overrightarrow{\upsilon\_{\mathcal{S}}} \right|}{\mu\_{\mathcal{S}}} \tag{5}$$

The Syamlal–O'Brien model is described as [15]:

$$\mathcal{K}\_{\mathcal{S}^\kappa} = \frac{3\alpha\_s \alpha\_{\mathcal{S}} \rho\_{\mathcal{S}}}{4\upsilon\_{r,s}^2 d\_s} \mathbb{C}\_d \Big(\frac{\mathcal{R}\mathcal{e}\_s}{\upsilon\_{r,s}}\Big| \overrightarrow{V}\_s - \overrightarrow{V}\_{\mathcal{S}}\Big|\tag{6}$$

where, *vr,s*, is the terminal velocity correlation for the solid phase as described below:

$$w\_{l,s} = 0.5 \left( A - 0.06 Re\_s + \sqrt{(0.06 Re\_s^2) + 0.12 Re\_5 (2B - A) + A^2} \right) \tag{7}$$

$$A = \;a\_{\mathcal{g}}^{4.14}, B = \;0.8a\_{\mathcal{g}}^{1.28} \text{ for } \alpha\_{\mathcal{g}} \le 0.85 \text{ and } B = \;0.8a\_{\mathcal{g}}^{2.65} \text{ for } \alpha\_{\mathcal{g}} > 0.85 \tag{8}$$

The EMMS model is described as [27]:

$$K\_{\mathbb{S}^8} = \frac{3}{4} c\_D \frac{\alpha\_{\mathbb{S}} \alpha\_{\mathbb{S}} \rho\_{\mathbb{S}} |\overrightarrow{V}\_{\mathbb{s}} - \overrightarrow{V}\_{\mathbb{S}}|}{d\_{\mathbb{S}}} \alpha\_{\mathbb{S}}^{-2.65} H\_d \text{ for } \alpha\_{\mathbb{S}} \ge 0.75 \tag{9}$$

$$K\_{\mathcal{S}^\text{s}} = 150 \frac{\alpha\_s^2 \mu\_\mathcal{g}}{\alpha\_\mathcal{g} d\_\mathcal{s}^2} + 1.75 \frac{\alpha\_s \rho\_\mathcal{g} |\overrightarrow{V}\_s - \overrightarrow{V}\_\mathcal{g}|}{d\_\mathcal{s}} \text{ for } \alpha\_\mathcal{g} < 0.75 \tag{10}$$

The heterogeneity index, *Hd*, is calculated based on a proposed method elsewhere [28].

#### 2.1.4. Turbulent Models

The transport equations for the RNG k-ε and k-ω models are written as below, where the buoyancy effect is neglected [29]:

$$\frac{\partial}{\partial t}(\rho k) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho k u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_l}{\sigma\_k}) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + P\_k + \rho \varepsilon \tag{11}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \frac{\partial}{\partial \mathbf{x}\_i}(\rho \varepsilon u\_i) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_t}{\sigma\_\varepsilon}) \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \right] + \mathbf{C}\_{\varepsilon 1} \frac{\varepsilon}{k} P\_k - \mathbf{C}\_{\varepsilon 2}^\* \rho \frac{\varepsilon^2}{k} \tag{12}$$

*Processes* **2019**, *7*, 524

where *C*∗ <sup>ε</sup><sup>2</sup> <sup>=</sup> *<sup>C</sup>*ε<sup>2</sup> <sup>+</sup> *<sup>C</sup>*μη<sup>3</sup> <sup>1</sup><sup>−</sup> <sup>η</sup> η0 <sup>1</sup>+βη<sup>3</sup> , <sup>η</sup> <sup>=</sup> *Sk* <sup>ε</sup> , *<sup>S</sup>* <sup>=</sup> <sup>2</sup>*SijSij*1/2 and Pk is the production rate of turbulence. Transport equations for the k-ω model suggested by Wilcox (1991) [30] are described as:

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\rho l I\_j k) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\mu + \frac{\mu\_l}{\sigma\_k}) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + P\_k - \not{\rho} k \omega \nu + P\_{k\flat} \quad \text{``K-equation''} \tag{13}$$

$$\frac{\partial(\rho\omega)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j}(\rho \mathbf{U}\_j \omega) = \frac{\partial}{\partial \mathbf{x}\_j} \Big[ (\mu + \frac{\mu\_l}{\sigma\_\omega}) \frac{\partial \omega}{\partial \mathbf{x}\_j} \Big] + a \frac{\omega}{k} P\_k - \beta \rho \omega^2 + P\_{\alpha b} \quad \text{`` $\alpha$ -equation''} \tag{14}$$

The numerical values of the coefficients for the RNG k-ε and k-ω models are available in the literature [31].

#### *2.2. Models' Descriptions and Simulation Methods*

#### 2.2.1. Initial Conditions, Boundary Conditions, Geometry and Grids

The geometry, boundary conditions and grids of the reactor at the initial state are defined and shown in Figure 1a,b, in order to study the fluidization of the granular biomass particles.

**Figure 1.** (**a**) Reactor geometry, (**b**) boundary zones and grid design at initial state (static bed).

Air (i.e., the fluidizing agent) enters the reactor from the bottom and passes through a distributor plate prior to its injection into the chamber. The producer gas exits from the top of the reactor where the pressure is set to atmospheric with no backflow pressure. The walls are thermally insulated and defined with a specularity coefficient of φ = 0.5 for the solid phase, as recommended for multiphase granular flow on bubbling fluidized beds and a no-slip condition for air [32]. Up to 20% of the reactor volume was initially packed with bed material (i.e., silica sand). The volume fraction of solid particles within the patched zone was adjusted to 60%. The ambient air temperature and pressure were set at the reference values.

In order to avoid instability and poor convergence in the simulation of the multiphase flows, a small time step-size of 0.001 s with 100 iterations per time-step was used until the relative convergence criterion of 1 <sup>×</sup> 10−<sup>3</sup> was satisfied. The physical specifications of the bed material are summarised in Table 2.


**Table 2.** Physical specifications of bed materials (silica sand particle) [33].

A fine grid was selected to capture the steep gradients between the neighbouring areas. The non-uniform grid spacing was generated using Ansys Desing Modeler (Figure 1). The second order upwind scheme was used to discretise the convective terms in the momentum and energy equations. The values and parameters used in the grid generation and numerical model development are shown in Table 3.

**Table 3.** Parameters used in grid generation and model development of fluidization hydrodynamics.


#### 2.2.2. The Effect of Viscous Models on Pressure Drop Profile

The pressure drop profiles against different superficial velocities were derived for each viscous model (i.e., laminar, RNG k-ε and k-ω) and were then compared against the experimental data. It must be noted that the occurrence of pressure drop was due to continual fission and fusion of the bubbles in the fluidization zone [35]. The pressure drop values were recorded over the hypothetical vertical line aligned in the centre of the reactor's geometry (Figure 1a).

The best viscous model was selected based on how accurately it could describe the pressure drop profile when compared with the experimentally-obtained data. That model was subsequently used for further analysis. In order to identify such a model, a statistical method was employed to evaluate a two-sample-T-test and to calculate the mean squared error (MSE) for each curve. In order to ensure that the bed has been effectively fluidized, the velocities <minimum fluidization velocity (Umf) were discarded in the statistical comparison.

#### 2.2.3. Experimental Validation of the Numerical Analyses

Controlled experiments were conducted using a fluidized bed reactor in order to validate the simulation results. A reactor with similar dimensions and physical characteristics assumed in the modelling work, was used in experimental data collection. Pressure drop profiles were generated from the cold test evaluation when the reactor was initially packed with bed material with similar particle diameter (i.e., 550 μm) as assumed in the modelling, to the same level (i.e., 20% of total reactor volume) as described earlier in this work. Further details on the reactor specifications and experimental setup can be found in our previous study [36]. Each incremental drop in pressure throughout the reactor was caused by the gradual increase in the superficial velocity of the air between the air distributor and the reactor outlet.

#### 2.2.4. The Effect of Air Distributor on Fluid Properties

The nature of the turbulence caused by the passage of the fluid through the air distributor was investigated. The area below and above the air distributor plate located at the bottom of the reactor was modelled. The dynamics of the air passing through the distributor plates (with pore diameters of 2, 3 and 5 mm) was selected for our simulation. The properties of the fluid with a superficial velocity of 0.2 ms−<sup>1</sup> was described by the RNG K-ε turbulent model in this section. The selected superficial velocity was about five times greater than the minimum fluidization velocity to ensure complete fluidization [9,37].

A mesh-structured air distributor plate was placed at the bottom of the reactor (Figure 2a). It was necessary to provide a regulated momentum of air inside the reactor chamber to avoid the falling back of the solid particles into the distributor zone, which could consequently have blocked the air nozzles. Distributor plates with different pore sizes (Figure 2b) were used in this simulation to investigate the effect of each perforated plate on the turbulent viscosity of the air. Although the diameter of the plate's air holes was variable, the total area of the pores was intended to be constant by adjusting the number of the pores accordingly.

**Figure 2.** Geometry, boundary condition and grid design for the air distributor plate used in (**a**) experimentation and (**b**) numerical studies.

#### 2.2.5. The Effect of Drag Models on Fluidization Hydrodynamics

Drag force is an key dominant force between the solid particles and the gas phase in a fluidized bed reactor [38]. Typically, the drag law describing the inter-phase momentum exchange is described by experimentally-developed models. However, for a fluidized bed with particular characteristics, it was necessary to investigate whether such models are still applicable. Therefore, three drag models (i.e., Syamlal–O'Brien, Gidaspow and EMMS) were selected in this study to numerically investigate the effect of drag functions on the fluidization hydrodynamics. Initially, the pressure drop throughout the bed was recorded over time for each drag model. That helped to understand when the impact of pressure drop began to diminish and a steady-state was reached. That indicates complete fluidization where the bed has been expanded to its maximum target capacity. At that point, the hydrodynamics of gas–solid fluidization can be defined on a time-averaged basis. Subsequently, the effect of drag models on the volume fraction of the particles was investigated.

#### 2.2.6. Integrated Turbulent Model and Model Validation

A comprehensive model of the fluidization of the granular particles was be developed via the integration of the optimum parameters obtained for the most suitable viscosity, drag and pore size of the air distributor pore models for each stage. The contours of particles' volume fraction obtained from the turbulent model were compared with the standard laminar model in a similar setup. It must be noted that the term "solid" in this study refers to the particle materials in the bed (i.e., granular silica sands). It is of key importance to achieve a well-developed fluidized bed in terms of bed expansion and pressure gradient prior to injecting the biomass into the reactor. Moreover, in the actual process, the biomass particles degrade immediately as injected into the hot reactor and undergo a drastic loss in density while converting into the solid biochar within a very short period of time. The density of the produced char is about 40 times less than the bed material (i.e., 2636 kg m−1) (e.g., Napier grass biochar has the density of 71 kg m−<sup>1</sup> [36]). Therefore, it was not possible for the fluidizing

agent to overcome the drag force, and as a consequence, the particles would fly out of the reactor. Therefore, the characteristics of fluidization, such as the minimum fluidization condition, maximum bed expansion and the bubble formation are mostly governed by the motion pattern of the sand particles and are less effected by the biomass particles.

Once the final integrated model has been set up, it is capable of monitoring the solid volume fraction profile inside the reactor throughout the gasification process. It enabled us to also predict the parameters associated with the particles' motion and the physical aspects of the fluidization, including the kinetic parameters, pressure-related characteristics, formation and transportation of bubbles and the parameters associated with the fluidized bed expansion.

#### 2.2.7. Grid Resolution Study

Grid resolution study was conducted in order to check the accuracy of the results. This was done by pursuing the balance between coarseness of the mesh elements and the computation time required for a solution. Indeed, coarser grids require less computation time and vice versa. The procedure starts with selecting a coarse grid with the least number of elements and performing the computation until the results have been converged with lower relative errors than the pre-defined limit (i.e., <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>−</sup>4) for all models involved in this study. Similar grid sizes and numbers of elements were used in both studies of fluidization hydrodynamics and evaluation of temperature gradient. Therefore, a single grid study was valid for both models. The results of grid resolution tests are shown in Figures 3 and 4, respectively.

**Figure 4.** Grid resolution test for re-normalisation group (RNG) K-ε turbulent flow.

The output factor selected in this grid study was the maximum bed height and is related to the particles' superficial velocities along the y-axis. Maximum bed expansion occurs when a particle with the highest y value in all time steps has the velocity of zero. Zero y-velocity with highest possible elevation means that the particle is not travelling to any higher zones, and thus, merely bounces within the fluidization region. An accurate grid quality for a reasonable computation time ensures the highest bed expansion and better fluidization characteristics. For the laminar model, the computations were completed for six runs while each run was committed to a gradual increase in element numbers and eventually converged within a specific time window. Figure 3 shows that the bed height does not significantly increase with finer grid structure above 16,000 elements. Any further increase in element quantities above this point leads to an unnecessarily prolonged computing time. The same trend was observed in the turbulent model but with 58,000 elements (Figure 4). Those optimum grid qualities were employed in the evaluation of the models in this work.

#### **3. Result and Discussion**

#### *3.1. The E*ff*ect of Viscous Models on Pressure Drop Profiles*

The modelled pressure drop profiles against the superficial velocities corresponding to each viscous model are shown in Figure 5.

**Figure 5.** Pressure drop versus superficial velocity for different fluid regimes.

In the laminar regime, the uniform flow of the inlet air stream shows insignificant acceleration, and hence, appears to follow the Hagen–Poiseuille equation order [39]. According to Figure 5, although different pressure drop values are indicated by each model, the minimum fluidization velocities (Umf) are just slightly different. The values of Umf are 0.030, 0.029, 0.031 and 0.037 (ms<sup>−</sup>1) for experimental, laminar, RNG K-ε and K-ω models, respectively.

At lower superficial velocities, an increase in the pressure drop in the laminar model is seen with a constant slope until the maximum Δp is reached. This can be due to the uniform and vibrating motion of the static bed as well as the regular pressure drop gradient [40]. For the velocities below Umf, an increase in the superficial velocity leads to an increase in the pressure drop; however, the slope was observed to be steeper for the RNG K-ε and K-ω models compared to the laminar model. That can be explained by the different functions of each model in describing the particle interactions and parameters, such as friction coefficient, shear stress, formation of eddies and the associated transferred energy [41]. All models tend to converge and surpass at the onset of fluidization when the superficial velocity is over 0.03 ms−1. The Δp is then recorded at constant pressure of 1.016 kPa during the experimental trial: 0.980 kPa for laminar, 1.001 kPa for RNG K-ε and 1.032 kPa for K-ω models.

According to the statistical analyses, only the two laminar and RNG K-ε models with *p*-values of <0.05 were capable of adequately representing the experimental data. The RNG K-ε with a smaller MSE of 539 was selected as the most desirable model for this study to further describe the turbulent characteristics of the fluidization hydrodynamics.

#### *3.2. The E*ff*ect of Air Distributor on the Fluid Properties*

In this study, the effect of air distributor on the fluid properties was investigated and an effective turbulent viscosity was calculated to analyse the transfer of momentum caused by turbulent eddies. The corresponding contours of the turbulent viscosity across the distributor area are presented in Figure 6.

**Figure 6.** Effect of distributor plate with pore diameters of (**a**) 5 mm, (**b**) 3 mm and (**c**) 2 mm on the turbulent viscosity of the fluid.

The air distributor with larger pores (Figure 6a) was observed to be less effective in increasing the magnitude of the kinetic energy of the fluid. Moreover, the pressure gradient is not uniformly distributed above the surface of the plate. The energy that is introduced by eddies are maximized around the corners which may result in pushing the solid particles to travel along the wall zones and less towards the centre of the reactor. This may cause a sudden pressure drop across the centre which simultaneously reduces the bed expansion. Plates with pore diameters of <5 mm are, therefore, desirable. However, according to Figure 6c, very small pore diameters also tend to be non-ideal as they increase the turbulent viscosity, as well as the gas pressure below the distributor plate. Therefore, they create the red zone—an indication of the air failing to pass through the nozzles. This is due to the fact that the sudden pressure drop imposed by the very small pores leads to the overall pressure of the fluid below the plate to significantly increase [42]. Larger eddies are created across the plate and the chance of a back pressure developing increases. The velocity of the fluid above the distributor plate, on the other hand, is too small to vibrate the solid particles, hence is not able to expand the bed. A plate with a pore diameter of 3 mm (Figure 6b) was observed to have a better turbulent viscosity profile than the other two cases discussed. The distribution of energy above the plate is uniform and the viscosities are maximal (i.e., a 1.8 fold increase). Although the intensity of turbulent viscosity is higher around the wall area, it is relatively constant towards the centre of the plate.

A plot of the turbulent viscosity against the centre line of Y-coordinates is shown in Figure 6b. The turbulent viscosity of the air was slightly reduced from 1.79 <sup>×</sup> 10−<sup>4</sup> to 1.65 <sup>×</sup> 10−<sup>4</sup> kg(ms)−<sup>1</sup> from the heights of 0 to 0.18 m. However, as the fluid touches the distributor plate, the turbulent viscosity,

encountered by an intensive pulse, passes through the grate pores. At the beginning, it rapidly increases to 1.7 <sup>×</sup> 10−<sup>4</sup> kg(ms)−<sup>1</sup> at a height of 0.196 m and then dramatically drops to 0.58 <sup>×</sup> 10−<sup>4</sup> kg(ms)−<sup>1</sup> at a height of 0.2 m. It then increases to the maximum value of 2.69 <sup>×</sup> 10−<sup>4</sup> kg(ms)−<sup>1</sup> at point "A" at a height of 0.217 m and then falls again as it exits from the pores of the mesh structured plate, and flows along the height of reactor until it reaches 1.21 <sup>×</sup> 10−<sup>4</sup> kg(ms)−<sup>1</sup> at the height of 0.3 m. The resultant pulse was similar to the performance of an air nozzle, indicating a rapid rise in the kinetic energy by applying a sudden pressure drop to the fluid that is passing through the narrow channels of the nozzle [43]. For a better visualization of the phenomena, a closer view of the fluid momentum when passing through the nozzles is shown in Figure 7 with the assistance of velocity vectors.

**Figure 7.** Vector velocity coloured by turbulent viscosity on the pores of the distributor plate.

A sudden pressure drop is observed when the upwards-moving fluid reaches the small air nozzles. The air molecules are contracted as a result and their kinetic energy reduces accordingly while passing through the nozzles. Molecules with a high potential energy eject from the pores to above the distributor plate as their energy converts to kinetic energy in a very short time span, creating the eddies with the highest motions [42]. Since the fluid velocity is adequate and the pore diameters are defined properly, the turbulent viscosity profile becomes consistent from across the top of the plate and creates a steady pressure profile for fluidization. An air distributor plate with a pore diameter of 3 mm was, therefore, considered desirable in this study and will be used for further hydrodynamics investigations. Also, all the calculated quantities corresponding to the turbulent fluidizing fluid have been employed to derive the properties at point "A" in Figure 6b.

#### *3.3. The E*ff*ect of Drag Models on Fluidisation Hydrodynamics*

The results of pressure drop versus time for the three drag models used in this study are shown in Figure 8. It is observed that the pressure drop corresponding to all the three drag models have considerably reduced and hit a minimum at the beginning of the fluidization in approximately 1.5 s.

However, the pressure drop began to slightly increase until it reached a relatively steady-state condition with least fluctuation after almost 3 s.

Figure 9 indicates the variation of time-averaged solid volume fraction as the result of using the selected drag models along the height of the reactor at the central axis. It is seen that the models predicted by Syamlal–O'Brien and Gidaspow showed relatively lower solid volume fractions at any heights below 0.4 m as compared to the EMMS. However, the overall bed expansion for the model predicted by the EMMS was smaller (i.e., 0.54 m) than those forecast by Syamlal–O'Brien (i.e., 0.58) and Gidaspow (i.e., 0.6). Furthermore, the solid volume fraction at the surface of the fluidized bed gradually dropped to zero for the models predicted by Syamlal–O'Brien and Gidaspow, which seemed to be more physically realistic than those predicted by the EMMS.

**Figure 8.** Comparison of simulated bed pressure drop using Syamlal–O'Brien, Gidaspow and energy minimum multi-scale (EMMS) drag models (U <sup>=</sup> 0.2 ms<sup>−</sup>1; i.e., <sup>∼</sup>5 Umf).

**Figure 9.** Effect of drag models on the variation of time-averaged solid volume fraction along the reactor height at the central axis.

The EMMS model unrealistically predicts too sharp of a change in the solid volume fraction at similar zones. According to the pressure drop values corresponding to the three drag models (Figure 8), the formation of bubbles and the motion of particles are expected to follow a similar trend; however, the effects of each on the bed expansion and solid volume fraction are significant. The pressure drop corresponding to the minimum fluidization condition obtained experimentally (i.e., 1.166 kPa) was compared to that from the simulation of the drag models. The value obtained from the Syamlal–O'Brien model (i.e., 1.123 kPa) is observed to be more consistent with the experimental data, and therefore, it was considered as the reference drag model in the development of the hydrodynamic model.

#### *3.4. Development of the Fluidisation Hydrodynamics Turbulent Model*

The hydrodynamics of fluidization was simulated with the RNG K-ε turbulent model and the Syamlal–O'Brien drag model. An air distributor with a pore diameter of 3 mm was assumed to define the momentum exchange between the gas and solid phases, inter-particle interactions, bubble formation and motion pattern of the bed material in the state of fluidization. The initial air velocity and the turbulent viscosity were 0.2 ms−<sup>1</sup> and 2.69 <sup>×</sup> 10−<sup>4</sup> kg(ms)−1, respectively. The process was simulated for the initial 3 s. The contours for the solid volume fractions corresponding to different time steps-starting from the initial condition (i.e., t = 0) and continuing for three seconds of complete fluidization, are illustrated in Figure 10.

**Figure 10.** Solid volume fraction contours in a (**a**) laminar model and a (**b**) RNG k-ε turbulent model.

The formations of bubbles as well as the maximum bed expansion were noticeable at each time step. Both systems gradually reached a steady state and became fully-developed within three seconds once the fluidization was initiated. The particles in the laminar model (Figure 10a) are seen to have an axisymmetric distribution. The bubbles began to form at the surface of the air distributor and were enlarged as they moved upwards due to the surface tension gradient. Similar hydrodynamic behaviour has been reported previously [10,13]. However, as depicted in Figure 10b, the formation and development of bubbles as described by the RNG k-ε turbulent model seems to be non-uniform when compared to that described by the laminar model. Large quantities of small bubbles are firstly formed at the surface of the air distributor as opposed to the laminar model. These bubbles then gradually distributed towards the bed medium—more likely closer to the wall. The enlargement of the bubbles close to the wall was due to the friction caused by the wall roughness. The momentum imparted to the walls was introduced by the specularity coefficient for the solid phase at wall conditions. However, the bubbles in the middle section of the chamber were developed in size and quantity due to the interaction with the smaller nearby bubbles. All bubbles formed eventually tend to migrate to and collapse at the bed surfaces. Similar behaviour of bubbling fluidized beds is described in a previous study [44] as the bubbles started to form and expand in diameters prior to upwards movement through the bed medium towards the free board region. The movements of the bubbles, as they collide with

the stationary particles, transfers the kinetic energy and causes them to travel in random directions along the height of the gasifier and hence, establishing the fluidization phenomenon.

#### **4. Conclusions**

A numerical model with a Eulerian—Eulerian multiphase approach was employed to simulate biomass gasification in a bubbling fluidized bed reactor. Different turbulent models were investigated to study the hydrodynamics of the gasification process. The pressure drop was modelled along the gasifier's height for various superficial velocities, and the results were compared and validated against the experimentally-obtained data. The values corresponding to the RNG k-ε turbulent model were seen to have a better description of the experimental results. Moreover, the impact of the air distributor, with different pore diameters on the momentum characteristics of the air as it enters the gasifier, was investigated. The distributor plate with a pore diameter of 3 mm was selected for generating the largest turbulent viscosities above the distributor's plate (i.e., 2.69 <sup>×</sup> 10−<sup>4</sup> kg(ms)−1) with a uniform distribution profile. The drag models of Syamlal–O'brien, Gidaspow and EMMS were investigated in solid fluidization in terms of pressure drop profiles and the volume fraction. In comparison, the Syamlal–O'Brien model corresponding to a pressure drop of 1123 Pa was observed to be the most consistent with the experimental data (i.e., 1166 Pa). Upon the integration of the models' outputs, the contours for the solid volume fractions were generated and compared with the laminar model. The formation of the bubbles and their movement along the gasifier's height as well as the pressure drop profiles in the presence of the turbulent flow followed a different pattern compared to the laminar flow. The developed hydrodynamic model was able to predict the numerical terms associated with the particles' motions and the physical aspects of fluidization, including the kinetics, pressure-related characteristics, formation, and transportation of bubbles and fluidized bed expansion. The derived hydrodynamic model can also help to reduce the number of empirical trials required for process optimization for a fluidized bed reactor—providing accurate details on the physical aspects of fluidization which are challenging to experimentally measure.

**Author Contributions:** Conceptualization, R.K. and W.A.W.A.K.G.; Funding Acquisition, W.A.W.A.K.G.; Methodology, R.K.; Supervision, W.A.W.A.K.G., D.R.A.B., and R.Y.; Writing—Original Draft, R.K.; Writing—Review & Editing, W.A.W.A.K.G., S.M.S., M.S., S.R.M., and K.S.

**Funding:** This work has been funded and supported under the Ministry of Higher Education (MOHE) Malaysia via LRGS grant (Grant No. LRGS/2013/UKM/PT). Also, this research has been partially funded by the Engineering and Physical Sciences Research Council through the BEFEW project (Grant No. EP/P018165/1).

**Acknowledgments:** The authors would like to acknowledge the Chemical and Environmental Engineering Departments at University Putra Malaysia to support conducting this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
