Study of the Influence of the Mean Particle Diameter Choice and the Fractions Number on the Quality of Fluidized Bed Numerical Simulation

Abstract

We investigate the choosing of the fractions number for numerical simulation of a polydisperse bubbling fluidized bed using the Sauter mean diameter. The results were verified using experiments from a glass tube with a diameter of 2.2 cm and a height of 50 cm. As a fluidizing agent, air with a velocity of 0.0716 m/s to 0.1213 m/s was used. Polydispersed aluminum oxide particles with a diameter size of 20–140 µm were used as a solid phase. We propose a simple method for choosing the fractions number for the polydispersed granular phase in order to improve the quality of the numerical simulation results. In this study, we consider the Sauter mean diameter D₃₂ for each selected group of particles for the solid phase. By increasing the number of solid phase fractions, it is possible to obtain a mean boundary of the bubbling fluidized bed close to the observed experimental results. In our study, the division of polydispersed powder into four distinct solid-phase fractions enabled us to attain satisfactory agreement with experiments regarding the average value of the bed boundary.

1. Introduction

Fluidized bed reactors have many practical applications. These reactors have gained popularity in the following industries: oil refining [1,2], gasification [3,4], and petrochemistry [5]. The phase behavior of fluidized bed reactors can be studied both experimentally and mathematically.

Computational fluid dynamics (CFD) [6] is a modern tool for studying processes using numerical modeling. Currently, CFD is increasingly becoming the basis for understanding the physics of fluidized bed processes, replacing costly experimental studies. At the same time, modern CFD technologies are actively beginning to use artificial intelligence to increase accuracy and reduce computational costs [7,8].

A fluidized bed is a process that occurs in the presence of gas and solid particles. The selection and adjustment of the correct model for describing the behavior of a fluidized bed seems to be the primary task of a research engineer. One of the main models for describing the hydrodynamics of the behavior of gas (liquid) and solid particles in a fluidized bed is the Euler–Euler model (TFM—two-fluid model).

In the Euler–Euler model, each phase is regarded as a continuum, and the equations of mass and momentum conservation are formulated for each phase. The phases are interconnected through the inclusion of interphase forces in the equations, which comprise the lifting force, the drag force, and the virtual mass force [9]. The main disadvantage of the Euler–Euler model is that it is impossible to consider the behavior of particles separately when using this approach, since the model itself assumes the interaction of a continuum of interpenetrating phases [10,11,12,13,14]. Particles can move not only in a gas, but also in a liquid, and the same models of motion of multiphase media are applicable for calculations [15]. It is common to use the Euler–Euler multiphase model. This model requires less computation time but is able to describe the process with sufficient accuracy. Numerical modeling techniques can be used not only for the computation of vertical reactors, which typically take the form of a cylindrical column, but also for other issues. For instance, the authors of [16] conducted a simulation of a rotating drum using TFM.

Due to its superior efficacy, TFM has been extensively employed to simulate the thermal conversion of solid fuel in fluidized bed reactors [17]. The authors of [18] conducted an investigation into the fluid dynamics and heat transfer in a biomass fluidized bed gasifier by employing the two-dimensional Euler–Euler model. The authors of [19] conducted a simulation of the pyrolysis of combustible solid waste in a fluidized bed reactor by employing the Euler–Euler model. The numerical studies were conducted to examine the hydrodynamics, temperature, and emissions of harmful substances during the combustion of oxygen fuel in a combustion chamber equipped with a circulating fluidized bed [20]. The authors of [21] examined the process of biomass gasification in a bubble pseudo-fluidized bed using a combination of numerical modeling and experiments. In [22], the authors describe the results of an experimental study and numerical modeling of the process of isobutane dehydrogenation to isobutylene in a laboratory reactor with a fluidized bed and an alumina-chromium catalyst. In [23,24], the authors model an industrial reactor for the conversion of isobutane to isobutylene. In [25], the effect of internal grids on the expansion of the fluidized bed and particle separation in the reactor is considered. Fluidized bed processes have found wide applications not only in petrochemical applications, but also in energy applications, such as coal gasification processes. In [26], experiments are considered on different types of installations using nitrogen at low and high pressures.

Solid-phase particles may have non-uniform sizes in industrial reactors. Many catalysts use something called aluminum oxide, which can wear away over time. This will result in a modification of the particle size distribution during the operation of the reactor. Hence, it is imperative to consider particle polydispersity and the mechanisms of particle size distribution change when modeling fluidized bed reactors. In article [27], the author presents an overview and discusses the issues of attrition and heat transfer, their interaction, and the implications for the design of a fluidized bed system. The author also examine the effects of attrition on bed hydrodynamics and its impact on heat transfer and heat transfer correlations.

The utilization of diverse fraction groups with distinct particle sizes can be used in practice to alter the bed fluidization mode and enhance the efficiency of reactors. For example, studies [28,29,30,31] show bed expansion and an increase in the voidage of the dense phase due to the addition of small particles with a diameter of less than 45 µm. In other instances, the addition of these particles results in a decrease in the voidage of the solid phase [32,33]. The authors of [34] found that the addition of fine particles only slightly reduced the effective viscosity. In [35], experimental results showing a decrease in bubble size using fast X-ray tomography were obtained. The impact of fine particles on fluidization was examined in [36] through numerical simulations. It has been demonstrated that the addition of a significant number of fine particles suppresses bubbles, thereby enabling the apparatus to attain a stable solid-like operating mode. In studies [37,38], a bi-dispersed mixture of particles was employed in the isobutane dehydrogenation model, which had a positive impact on the expansion of the fluidized bed, heat transfer, and the yield of the reaction product. The value of 45 µm used in the study [28] cannot always be considered as a fine fraction. Such particles will not be fine for a very low gas velocity.

It is important to accurately model the behavior of particles in a reactor by correctly selecting the mean particle diameter. In general, the mean diameter is selected, which is called the Sauter diameter [39]. A comparison of the data presented in the literature for the Sauter diameter D₃₂ utilized indicates that, in certain instances, this choice may not provide a sufficient approximation of the results of calculating the concentration of particles in the layer to the results obtained in experimental studies.

The authors of [40] investigated the influence of particle size distribution on the results of computational particle fluid dynamics simulations. The authors found that an increase in average particle velocity leads to a decrease in pressure. They looked at how much the pressure drop varied depending on the average particle size (D₅₀), the geometric and arithmetic mean diameter, and the mean diameter according to Sauter (D₃₂). When the diameter was chosen according to Sauter, the relative standard deviation value was 5.8%, which is two times less than the relative standard deviation value obtained with D₅₀.

In a real fluidized bed study, the solid granular phase can consist of several different sizes. This can be a consequence of the destruction of the original particles or the intentional addition of fine particles to increase the efficiency of the apparatus operation. The correct choice of the particle mean diameter will allow for a highly accurate numerical simulation of the process under study. However, this is actually a true statement for a well-mixed fluidized bed, when we observe the same distribution of particles in each small volume. Obviously, fine particles are easier to lift upward than larger particles. Therefore, in some fluidization modes, particle segregation by size can be observed along the height of the apparatus. In [41], the authors developed a model for predicting the segregation directions of binary mixtures of particles with differences in density and size in bubbling fluidized beds. In [42], a comprehensive theoretical basis for regulating the mixing and segregation of particles in a fluidized bed was developed. Thus, in some zones of the studied area there will be more fine particles, and in other zones there will be more large particles. Then, the use of one mean diameter to describe the entire fluidized bed will be incorrect. In this case, it is possible to divide the studied particles into several fractions and consider each of them separately in numerical modeling. In previous studies [43,44], we identified fine particles using the terminal velocity parameter. Then, the average diameter D₃₂ of large and small particles was determined. The use of small particles allowed us to bring the bed boundary closer to the experimental values in numerical calculations. However, the results obtained have not yet allowed us to achieve good agreement with the experiment. This can be caused by several reasons: incorrect choice of the number of fractions, incorrect choice of the method of separating fractions by particle size, and incorrect choice of the particle mean diameter.

Thus, for numerical simulation of a fluidized bed with polydispersed solid particles, the researcher must solve additional problems: determine the number of particle fractions, determine the method of fraction choice, and determine the model of the mean particle diameter. When choosing the mean particle diameter, standard models D₃₂ or D₅₀ can be used. For a perfect determination of the fraction number and the method of their choice, the researcher must know the particle material, changes in the particle distribution function, the fluidization mode, and the geometric features of the apparatus. In this case, a positive result will most likely show a choice of two or three main fractions. In this study, we propose a simpler method when the researcher only needs to know the particle size distribution function. We will separate the fractions into equal parts by particle size. Thus, a positive modeling result can be achieved with a larger number of fractions. But this method is acceptable when there is no complete information about the particles or the fluidization mode.

The present study was conducted to determine the particle number fractions for numerical simulation of a bubbling fluidized bed with polydispersed powder. We are interested in determining when numerical simulation with the Sauter mean diameter D₃₂ is appropriate. To verify the findings, we utilize the outcomes of experimental investigations conducted on a glass tube with a diameter of 2.2 cm and a height of 50 cm. Air having a velocity ranging from 0.0716 m/s to 0.1213 m/s is employed as a fluidizing agent. Polydispersed aluminum oxide particles with a particle size of 20–140 µm are used as a solid phase. We intend to divide the particle size distribution function into various groups in order to enhance the outcomes of numerical modeling. Obviously, a large number of groups of different sizes can help us more closely approach real experimental data. However, an increase in the number of particle fractions will result in a significant increase in the cost of computing resources.

2. Materials and Methods

2.1. Experiment Setup

In this study of fluidized bed hydrodynamics, we used a laboratory reactor, which is a vertical glass tube with a height of 50 cm and a diameter of 2.2 cm (Figure 1a,b) [37,38,43,44]. The thickness of the glass tube is 1.2 mm. We selected the glass tube for the experiment owing to the excellent transparency of the wall material. In our study, we obtain video data from the experiment to determine the boundary of the particle bed at a moment in time. The air flow is supplied to the bottom of the glass cylinder through a multi-layer porous material, which ensures uniform flow distribution across the cross-section of the tube and prevents solid particles from entering the gas supply system. Aluminum oxide particles with a polydispersed fractional composition are placed in the reactor. Initially, a fixed bed of particles with a total mass of 30 g is poured. The density of the solid particle material is 2400 kg/m³. The height observed for the fixed bed is 0.079 m. Subsequently, air is introduced to the bottom of the tube, thereby initiating the fluidization process. The process of fluidization was recorded with a high-speed camera, and the bed height was estimated from the camera frames. The resulting frames were used to determine the average and boundary values. We used our own code, which is written in Python and based on intelligent computer vision technologies. The code used automatically detects the bed boundaries during the fluidization process.

From the experimental data, we can determine the initial volume fraction of the solid dispersed phase. We will use this information to adjust the numerical simulation model. We know the inner diameter of the tube D = 2.2 cm and observe the height of the fixed bed H = 7.9 cm. The calculated volume V = H × π × D²/4 = 30.02 cm³ contains solid particles and air. The density of the particle material is determined during production [45] and is ρ = 2400 kg/m³ = 2.4 g/cm³. The measured mass of the poured bed of particles is m = 30 g. Then, the volume of the solid phase is V_s = m/ρ = 12.5 cm³. Thus, we determine the initial solid volume fraction of the bed as α_s = V_s/V = 0.41639.

In laboratory experiments, we investigated a polydispersed powder of microspherical particles with a diameter of 20–140 µm. The particles employed are intended for the production of a catalyst for isobutane dehydrogenation [45,46,47]. The granules belong to Geldart group B. The cumulative distribution function, the probability distribution function and the discrete probabilities of the particle size are shown in Figure 1b. Particle size distribution was determined using a Mastersizer 2000 analyzer and its software v5.60.

For a wide range of air flow rates, experimental studies were carried out to observe bubble fluidization with the formation of a height-limited particle bed. In this investigation, we examine four values of air velocity, namely 0.0716 m/s, 0.0892 m/s, 0.1088 m/s, and 0.1213 m/s.

As a result of the experimental study, we evaluated the height of the particle bed. During the fluidization process, bubbles of various sizes are generated, which are destroyed at the upper boundary of the bed. Therefore, the boundary of the particle bed is non-stationary. The minimum and maximum boundaries of the bed are observed, as well as the average value of the bed boundary during the observation period. The results are presented in

Table 1. Measurement result from experiments.

Air Velocity, m/s	Observed Minimum Height, m	Observed Average Height, m	Observed Maximum Height, m
0.0716	0.107	0.1386	0.18
0.0892	0.12	0.1433	0.19
0.1088	0.13	0.1683	0.23
0.1213	0.13	0.1721	0.24

.

2.2. Particle Mean Diameter

The Euler–Euler multiphase model necessitates the discretization of particle characteristics. The discretization of particle characteristics can be accomplished through a multigroup approach, wherein the probability distribution function transforms into a collection of scalars that correspond to a comprehensive group of particle sizes. For each group, it is necessary to solve the transport equations.

The simplest CFD approaches that employ phase interaction assume the existence of particles with the same diameters. It is common practice to identify the mean size of polydispersed particles using statistical or other models. When applying the multigroup approach, the mean particle size must be determined for each group.

First, consider the models used for selecting one mean particle diameter size. If the particle size distribution function by diameter sizes exhibits a normal form, then all models for determining the average diameter will exhibit the same value. In our case, the size distribution function (Figure 2) differs from the normal form. We calculate the median diameter D₅₀ by separating half of the particles by volume fraction. We also consider the statistical modal diameter D_mode. Moreover, it is imperative to take into consideration the Sauter mean diameter D₃₂.

When modeling a fluidized bed, the Sauter mean diameter D₃₂ is frequently employed, which is determined by utilizing the particle diameter distribution density function [39]:

$$D_{32}={\displaystyle \frac{{\displaystyle \underset{0}{\overset{\infty }{\int }}f(D)D^{3}dD}}{{\displaystyle \underset{0}{\overset{\infty }{\int }}f(D)D^{2}dD}}}.$$

(1)

The distribution density function satisfies the condition:

$${\displaystyle \underset{0}{\overset{\infty }{\int }}f(D)dD}=1.$$

(2)

For the distribution function shown in Figure 1, we obtained the following values: D₅₀ = 56.705 µm, D_mode = 56.3685 µm, and D₃₂ = 69.8087 µm. The results are also shown in Figure 2.

From Figure 2 and the obtained numerical values, it is evident that D₅₀ and D_mode have close values, and the D₃₂ value is approximately 13 µm greater. In [37], we studied the selection of one effective diameter for similar experimental parameters. It was shown that good agreement with experimental results can be achieved by choosing D = 60 µm. This value is incompatible with any analytical model.

We will then consider the use of a multi-group approach, with division into several fractions and determination of the mean particle size for each group. It is obvious that a large number of different groups of different sizes can allow us to more closely approach the real problem. At the same time, if there are more fractions, it will make numerical modeling more expensive. The particle distribution will be divided into several groups of fractions in this investigation. In this case, we will not use physical parameters to separate the groups, but will instead use uniform division. First, the number of particle groups will allow us to bring the numerical modeling closer to the experimental data. Second, such a simple division is based only on knowledge of the minimum and maximum size of the particles used and does not require additional knowledge of the gas and solid phases. The mean diameter values D₃₂ are presented in

Table 2. Values of mean diameters D₃₂ for different numbers of separated fractions.

Fraction Numbers	D32, µm, Group 1	αs Group 1	D32, µm, Group 2	αs Group 2	D32, µm, Group 3	αs Group 3	D32, µm, Group 4	αs Group 4	D32, µm, Group 5	αs Group 5
1	69.8087	0.41639	-	-	-	-	-	-	-	-
2	48.6256	0.160518	76.4034	0.255872	-	-	-	-	-	-
3	34.0683	0.00941	59.6621	0.316956	81.3161	0.090024	-	-	-	-
4	26.3030	0.00201	48.6609	0.158508	69.7646	0.219063	88.0555	0.036809	-	-
5	23.0010	0.00083	38.7122	0.033186	55.2646	0.214233	74.0327	0.13133	93.0555	0.036809

and Figure 3.

The particle mean diameter for each group was determined using the Sauter model. For this, we used Formula (1), but the integration was performed within the limits of the size variation for each group. The dashed lines in Figure 3 show the boundaries of the selected groups.

In Figure 3d, we see that the selected finest fraction is almost not represented in the volume fraction. We can obtain the same information from Table 2. In fact, the boundary regions are of special interest for our problem. In the present study, we took the cumulative function, which shows the presence of particles in the range from 0 to 26 μm. When we use the Sauter mean diameter model for the entire particle distribution (the case of one group), we integrate from zero to infinity. When dividing into several groups, we also started the integration from zero for the first group boundary. When there were few groups (for example, 2 or 3), this did not play a significant role. However, when dividing into five groups, we obtained a fine fraction with a volume fraction of 0.00083. Obviously, this fraction will not make a significant contribution to the model. With further division into groups, we obtain the finest phase with the near-zero fraction. This will also be a disadvantage of our model. The model should be improved in further studies to avoid similar situations.

The values of the solid-phase volume fraction presented in Table 2 are calculated for the total volume, which initially contains 0.58361 air, and the initial solid-phase fraction is α_0s = 0.41639. Let us consider the normalized solid-phase fractions for each group relative to the initial total solid-phase fraction α_0s. The results are presented in

Table 3. Normalized solid fractions for each group in present model.

Fraction Numbers	αs/α0s Group 1	αs/α0s Group 2	αs/α0s Group 3	αs/α0s Group 4	αs/α0s Group 5
1	1	-	-	-	-
2	0.3854992	0.6145008	-	-	-
3	0.0225990	0.7611998	0.2162012	-	-
4	0.0048272	0.3806720	0.5261005	0.0884003	-
5	0.0019933	0.0796993	0.5145008	0.3154014	0.0884003

.

From Table 3, we see that for the case of three groups, the normalized fraction of group 1 will be about 2%. This group will make a minor contribution to the overall numerical model of the fluidized bed. For the case of four groups, the normalized fraction of group 1 will be less than 0.5%. For the case of five groups, the normalized fraction of group 1 will be less than 0.2%. This also shows the limitation of the present model when groups with very small volume fractions can be considered in the calculation.

2.3. Numerical Simulation Model

In order to perform numerical simulation, we must define a geometry model. Since the experimental tube has a simple shape, we choose to use a cylinder as the computational domain. The scheme is shown in Figure 4a. Within the lower surface of the cylinder, we define the condition of supplying the gas velocity value; for the granular solid phase, the surface will be impermeable. Within the upper surface of the cylinder, we define a simple exit condition. Since heat transfer is not considered in this problem, the walls are set to be impermeable without additional thermal conditions. As the initial conditions for numerical modeling, we established a fixed bed of solid particles with a specified volume fraction within the lower zone of the cylinder. From the experiment, the initial bed height of 0.079 m and the initial volume fraction of the solid substance of 0.41639 were determined. The density of the solid phase corresponds to the material value from the experiment of 2400 kg/m³. The standard for air density and viscosity is determined at a temperature of 20 °C.

After the computation start, upon the introduction of gas from below, the particle bed will commence to expand. After a certain period, the expansion of the bed will cease, and we shall consider it to have reached the fluidization mode at a given gas velocity. Next, we observe the parameters of the fluidized bed for some time and average the results.

2.4. Fluidized Bed Model

In the numerical simulation of the fluidized bed, a two-fluid continuous Euler–Euler multiphase model was used, supplemented by the kinetic theory of gases to take into account the collisions of solid particles. It is assumed that all considered phases in total occupy the entire volume of the computational domain. Let us denote the volume fraction of the i-th phase by ${\alpha }_i$, respectively; the sum of the volume fractions of all phases is equal to 1. For each of the phases, the equations of conservation of mass, momentum and energy are satisfied [48,49]. In the study of this problem, the following equations were solved.

The mass conservation equation for the gas phase:

$${\displaystyle \frac{\partial {\alpha }_g{\rho }_g}{\partial t}}+\nabla \cdot \left({\alpha }_g{\rho }_g{\overrightarrow{v}}_g\right)=0,$$

(3)

where ${\rho }_g$ is the real density of the gas phase, and ${\overrightarrow{v}}_g$ is the velocity of the gas phase.

The mass conservation equation for the solid phase:

$${\displaystyle \frac{\partial {\alpha }_s{\rho }_s}{\partial t}}+\nabla \cdot \left({\alpha }_s{\rho }_s{\overrightarrow{v}}_s\right)=0,$$

(4)

where ${\rho }_s$ is the real density of the solid phase, and ${\overrightarrow{v}}_s$ is the velocity of the solid phase.

The momentum conservation equation for the gas phase:

$${\displaystyle \frac{\partial {\alpha }_g{\rho }_g{\overrightarrow{v}}_g}{\partial t}}+\nabla \cdot \left({\alpha }_g{\rho }_g{\overrightarrow{v}}_g{\overrightarrow{v}}_g\right)=-{\alpha }_g\nabla p+\nabla \cdot {\overline{\overline{\tau }}}_g+{\alpha }_g{\rho }_g\overrightarrow{g}+{\displaystyle \sum _{j=1}^N{K}_{g\hspace{0.17em}s}\left({\overrightarrow{v}}_g-{\overrightarrow{v}}_{s,j}\right)},$$

(5)

where $p$ is the pressure, ${\overline{\overline{\tau }}}_g$ is the stress tensor in gas phase, and $K_{gs}$ is the coefficient of interaction between the gas and discrete solid phases. Here, the subscript s,j means the interaction of the gas with the j-th solid phase for cases where more than one solid phase is present in the model. In Equation (5), the stress tensor is calculated as follows:

$${\overline{\overline{\tau }}}_g={\alpha }_g{\mu }_g\left(\nabla {\overrightarrow{v}}_g+\nabla {\overrightarrow{v}}_g^T\right)+{\alpha }_g{\displaystyle \frac{2}{3}}{\mu }_g\nabla \cdot {\overrightarrow{v}}_g\overline{\overline{I}},$$

(6)

where ${\mu }_g$ is the shear viscosity, and $\overline{\overline{I}}$ is the unit tensor.

The momentum conservation equation for the solid phase is

$${\displaystyle \frac{\partial {\alpha }_s{\rho }_s{\overrightarrow{v}}_s}{\partial t}}+\nabla \cdot \left({\alpha }_s{\rho }_s{\overrightarrow{v}}_s{\overrightarrow{v}}_s\right)=-{\alpha }_s\nabla p-\nabla p_s+\nabla \cdot {\overline{\overline{\tau }}}_s+{\alpha }_s{\rho }_s\overrightarrow{g}+K_{s\hspace{0.17em}g}\left({\overrightarrow{v}}_s-{\overrightarrow{v}}_g\right)+{\displaystyle \sum _{j=1}^N{K}_{s\hspace{0.17em}s,j}\left({\overrightarrow{v}}_s-{\overrightarrow{v}}_{s,j}\right)},$$

(7)

where $p_s$ is the pressure of solid-phase granules, ${\overline{\overline{\tau }}}_s$ is the stress tensor in the solid phase, $K_{sg}$ is the coefficient of interaction between discrete solid phases and gas, and $K_{ss.j}$ is the coefficient of interaction between operating solid phases and another j-th solid phase. For the interaction of gas and solid phases, we have $K_{gs}=K_{sg}$. In Equation (7), the stress tensor is calculated as follows:

$${\overline{\overline{\tau }}}_s={\alpha }_s{\mu }_s\left(\nabla {\overrightarrow{v}}_s+\nabla {\overrightarrow{v}}_s^T\right)+{\alpha }_s\left({\lambda }_s-{\displaystyle \frac{2}{3}}{\mu }_s\right)\nabla \cdot {\overrightarrow{v}}_s\overline{\overline{I}},$$

(8)

where ${\mu }_s,\hspace{0.17em}\hspace{0.17em}{\lambda }_s$ are shear and bulk viscosity of the solid phase.

The conservation of solid particle phase fluctuating energy [50] is calculated as follows:

$${\displaystyle \frac{3}{2}}\left[{\displaystyle \frac{\partial }{\partial t}}\left({\alpha }_s{\rho }_s{\Theta }_s\right)+\nabla \cdot \left({\alpha }_s{\rho }_s{\overrightarrow{v}}_s{\Theta }_s\right)\right]=\left(-p_s\overline{\overline{I}}+{\overline{\overline{\tau }}}_s\right):\nabla {\overrightarrow{v}}_s+\nabla \cdot \left(k_{{\Theta }_s}\nabla {\Theta }_s\right)-{\gamma }_{{\Theta }_s}+{\varphi }_{g\hspace{0.17em}s}+{\displaystyle \sum _{j=1}^N{\varphi }_{s\hspace{0.17em}s,j}},$$

(9)

where $k_{{\Theta }_s}$ is the granule energy diffusion coefficient, ${\Theta }_s$ is the solid granule temperature, ${\gamma }_{{\Theta }_s}$ is the particle collisions energy dissipation, ${\varphi }_{gs}$ is the energy exchange between solid and gas phases, and ${\varphi }_{ss,j}$ is the energy exchange between a pair of solid phases.

Further, we consider the constitutional relations for the closure of the equations system presented above. In the case of the interaction of two phases of the gaseous–solid granular type, the model as described in [49] is used, as follows:

$$K_{sg}=\left\{\begin{array}{l}{\displaystyle \frac{3}{4}}{C}_D{\displaystyle \frac{{\alpha }_s{\alpha }_g{\rho }_g\left|{\overrightarrow{v}}_s-{\overrightarrow{v}}_g\right|}{d_s}}{\alpha }_g^{-2/65},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\alpha }_g>0.8,\\ 150{\displaystyle \frac{{\alpha }_s\left(1-{\alpha }_g\right){\mu }_g}{{\alpha }_g{\mu }_s^2}}+1.75{\displaystyle \frac{{\rho }_g{\alpha }_s\left|{\overrightarrow{v}}_s-{\overrightarrow{v}}_g\right|}{d_s}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\alpha }_g\le 0.8,\end{array}\right.$$

(10)

where drag coefficient is calculated as follows:

$$C_D={\displaystyle \frac{24}{{\alpha }_g{\mathrm{Re}}_s}}\left[1+0.15{\left({\alpha }_g{\mathrm{Re}}_s\right)}^{0.687}\right],$$

(11)

The Reynolds number is calculated as follows:

$${\mathrm{Re}}_s={\displaystyle \frac{{\alpha }_g{\rho }_g{d}_s\left|{\overrightarrow{v}}_s-{\overrightarrow{v}}_g\right|}{{\mu }_g}}.$$

(12)

The application of the interphase interaction model [49] has shown good agreement with the experimental and industrial results in our previous studies [24,25,26].

For the solid phase, the values of shear and bulk viscosity, granule pressure, diffusion and energy dissipation are determined using standard models [49,50,51,52].

The coefficient of interphase interaction between a pair of solid phases is determined according to the model [53]:

$$K_{ss,j}={\displaystyle \frac{3\left(1+e_{ss,j}\right)\frac{\pi }{2}{\alpha }_s{\rho }_s{\alpha }_{s,j}{\rho }_{s,j}{\left(d_s+d_{s,j}\right)}^2{\vartheta }_{0ss,j}}{2\pi \left({\rho }_s{d}_s^3+{\rho }_{s,j}{d}_{s,j}^3\right)}}\left|{\overrightarrow{v}}_{s,j}-{\overrightarrow{v}}_s\right|$$

(13)

where $e_{ss,j}$ is the coefficient of recovery of solid-phase particles after a collision, and ${\vartheta }_{0ss,j}$ is the radial distribution coefficient. Typically, the default value is $e_{ss,j}$ = 0.9, but this value can be adjusted to suit the particle type. The radial distribution coefficient ${\vartheta }_{0ss,j}$ is determined using model [49].

Considered in the glass tube, the fluidized bed component movement is turbulent. In the solution model used for calculations, a dispersed $k-\epsilon$ turbulence model was used, in which the motions of the “secondary” solid granular phases are generated against the turbulent motion of the “primary” gas phase. The ANSYS Fluent 19.2 software has been used to solve the written system of equations. Many of the component property parameters are taken from the database of the software.

2.5. Terminal Velocity

Let us consider the formula for determining the value of the terminal velocity. This study will help us determine which particles will be carried away from the apparatus at a given gas velocity. We consider the terminal velocity in the form [5]:

$$v_t={\left[{\displaystyle \frac{4gD({\rho }_s-{\rho }_g)}{3{\rho }_g{C}_D}}\right]}^{0.5},$$

(14)

where drag coefficient can be used from Formula (11). By solving Equation (14) for a given value of $v_t$, we can find the particle diameter D_t. All particles whose size is smaller than D_t will be carried out of the apparatus by the gas flow. The results of the calculated diameters for the studied velocity values are presented in

Table 4. Fine particle diameter for adjusted air velocity.

Air Velocity, m/s	Fine Particle Diameter, µm
0.0716	31.2
0.0892	34.5
0.1088	37.3
0.1213	39.1

.

Comparing the values in Table 4 and Table 2, we see that some fraction of the fine particles must be carried out of the fluidized bed apparatus for an adjusted gas velocity. For the two-fraction division, all of the selected particles are large enough to form a bounded fluidized bed. For the three-fraction division, we see that the fine fraction must be carried out of the apparatus for gas velocities of 0.1088 m/s and 0.1213 m/s. For the four-fraction division, the finest fraction will be carried out for all gas velocities studied. For the five-fraction division, the finest fraction will also be carried out for all gas velocities studied. In this case, the second particle group will most likely be located throughout the entire height of the apparatus and will be slowly carried out of it.

2.6. Mesh Size

We need to cover the computational domain with a mesh for using computational fluid dynamics methods. Since the geometric model of the computational domain has small dimensions and simple shapes, we used a uniform mesh with a cell size of $\Delta x$ = 0.001 m (Figure 4b). We conducted a mesh study in order to select the optimal size for solving our problem. Let us consider the cell sizes $\Delta x$ = 0.005 m, 0.0025 m, 0.001 m, and 0.0005 m. Mesh examples are presented in Figure 5. We will examine the volume fraction of the solid phase in the tube and the pressure drop between the inlet and outlet surface.

The numerical simulation represents 15 s of real time. In order to determine the average volume fraction of the solid phase, we record and average the results from the 5th to the 15th second of the calculation, with a step of 0.025 s. The results of the mesh study for a case with a gas velocity of 0.0892 m/s are shown in Figure 6. The results of testing the mesh size $\Delta x$ for the time-step value $\Delta t$ = 0.0005 s are shown in Figure 6a. It is observed that reducing the mesh size almost does not alter the volume fraction profile of the solid phase after the value of $\Delta x$ = 0.001.

Next, we present the results of the study of the pressure drop for gas passing through a fluidized bed of particles (Figure 6b). We observe significant deviations for the value of $\Delta x$ = 0.005 m. With a further decrease in the mesh size, the pressure drop values reach a constant value. There are small differences for the cases of $\Delta x$ = 0.0025 m and $\Delta x$ = 0.001 m. For the cases of $\Delta x$ = 0.001 m and $\Delta x$ = 0.0005 m, the differences are less than 0.5 Pa.

We also looked at the y+ parameter. This value is important for problems with the formation of a boundary layer of gas when moving along a wall. The study of y+ shows the need to use additional wall functions with insufficient grid resolution. In our problem, there is no full-fledged continuous gas motion along the walls. The presence of solid particles in the fluidization mode can form gas flows in any direction near the wall. However, the study of y+ can additionally evaluate the quality of the selected cell size. In our problem, we found that the y+ value varies from 1.3 to 3.65 for different values of $\Delta x$. For the selected mesh size $\Delta x$ = 0.001 m, the average value of y+ = 1.56. These values are obtained from the ANSYS Fluent software during numerical simulation. The selected mesh size does not require the use of additional wall functions in the turbulence model used.

The optimal mesh size for calculations in our model is $\Delta x$ = 0.001 m.

2.7. Time Step

Also, for all sizes $\Delta x$, we used four time-step values: $\Delta t$ = 0.0025 s, 0.001 s, 0.0005 s, and 0.00025 s. The results of the test on the step $\Delta t$ for the mesh size of $\Delta x$ = 0.001 m are shown in Figure 7a. It is seen that changing the time step almost does not affect the volume fraction of the solid phase.

Figure 7b shows the results of the pressure drop calculation for all the values of $\Delta x$ and $\Delta t$ used. A change in the value of $\Delta t$ also does not greatly affect the pressure drop value.

To select the time step $\Delta t$, we estimated the Courant number [54]:

$$N_C={\displaystyle \frac{\left|\overrightarrow{v}\right|\Delta t}{\Delta x}}.$$

(15)

The values $\Delta t$ were chosen considering the fulfillment of the criteria for the Courant number N_C < 1. However, numerical modeling of a fluidized bed may have its own characteristics, different from modeling the motion of a homogeneous gas. A maximum Courant number of 0.3 was suggested by the authors of [55]. The authors of [56] proposed establishing a limit on the Courant number for monodispersed liquid–fluidized suspensions, ranging from 0.03 to 0.3. The gas velocity ranges from 0.0716 to 0.1213 m/s in our studies. The time step should be between 0.000247 and 0.00419 s for such $\Delta x$ = 0.001 m and velocity values. The optimal time step for calculations in our model is $\Delta t$ = 0.0005 s. In all calculations, sufficient convergence is observed for all parameters for the selected values $\Delta x$ and $\Delta t$.

3. Results and Discussion

The results of a numerical simulation with one mean diameter D₃₂ are presented first. Figure 8 shows the volume fraction fields for the solid phase for a gas velocity of 0.0892 m/s at different times. The scale for the volume fraction is selected from 0 to 0.5. It is evident that the bubbling fluidized bed is already formed at the fourth second of the numerical simulation calculation. Similar results were obtained for other gas velocity values. We recorded and averaged the results from the 5th to the 15th second of the calculation with a step of 0.025 s, to analyze the average volume fraction profile. The results obtained for the average volume fraction of the solid phase are depicted in Figure 9 for four different gas velocity values.

Figure 9 shows that the calculation results are within the minimum and maximum values of the bed boundary that were found in the experiments. Note that these boundaries may be different for each experiment because they only cover one extreme case. At the same time, our numerical simulation results do not match the average experimental values. Numerical calculation results are lower than experimental values in all cases. All particles of the solid phase are located in the lower part of the apparatus and are not carried away by the gas flow. The chosen diameter is too large, and the bed is not sufficiently expanded as a result of the numerical simulation. In the results of numerical simulations, the average height of the bed boundary varies from 12.26 cm for a gas velocity of 0.0716 m/s to 13.75 cm for a gas velocity of 0.1213 m/s.

Next, we will consider numerical modeling options for cases where several solid-phase fractions are selected. It should be noted that we are currently considering the cases of two, three, four, and five fractions. The volume fraction fields for the solid phase for a gas velocity of 0.0892 m/s at a time of 14 s are shown in Figure 10 for the cases of two and five groups. The volume fraction scale is selected from 0 to 0.2. It is observed that, for two solid-phase groups, particle fractions are well mixed and fill the entire height of the fluidized bed. For the case of five groups, the observed height of the fluidized bed is higher. In this case, the phase of the smallest particles is almost absent. The gas flow from the tube removes particles of the smallest size. According to the theoretical analysis of the terminal velocity, we can conclude that the finest fraction is carried out of the apparatus. Therefore, we do not visually see the fraction of group 1/5 in Figure 10.

Let us examine the results of the average volume fraction profiles for all the studied cases. Figure 11 shows the results of modeling using two solid groups. Next, Figure 12 illustrates the numerical simulation results using three solid-phase groups. Results of the numerical simulation using four solid-phase groups are presented in Figure 13. Figure 14 shows the results of numerical modeling using five groups of the solid phase.

A good approximation can be obtained for the gas velocity cases of 0.0716 and 0.0892 m/s by using two solid particle fractions. However, for the other two gas velocities (0.1088 m/s and 0.1213 m/s), the fluidized bed height remains insufficient to reach the average experimental line. For the gas velocity of 0.1213 m/s, particle separation can be observed in the upper part of the bed. Most of the small particles are located over a bed of larger particles. Despite the partial separation of different fractions, all particles remain in a bounded bed in the lower part of the apparatus.

Using three solid particle fractions shows a similar result for the bed height compared to using two fractions for the gas velocity cases of 0.0716 m/s and 0.0892 m/s. However, for the gas velocity cases of 0.1088 m/s and 0.1213 m/s, the average particle bed height decreased. Here, we cannot see how much of the finest fraction is there. The separation of the fine particle group plays a negative role in this case. According to the theoretical analysis of the terminal velocity, we can conclude that the finest particles will be completely removed from the apparatus for the cases of gas velocity of 0.1088 m/s and 0.1213 m/s. That is, we will lose particles that in our model occupy about 2.5% of the total solid volume fraction. This is not much. But the remaining fractions have a large particle size. At the same time, the presence of fine particles for the cases of gas velocity of 0.0716 m/s and 0.0892 m/s contributes to a greater expansion of the total particle bed.

The use of four groups of solid particles shows good agreement with the experimental results for the average height of the fluidized bed for all cases of gas velocity. The smallest fraction is almost completely carried away by the gas flow from the tube, but its volume fraction is significantly less than in the case of using three groups. It is already possible to observe separation of particles in the upper part of the bed for two gas velocity values of 0.1088 m/s and 0.1213 m/s.

The use of five groups of solid particles shows good agreement with the experimental results for all cases of gas velocity. But for a gas velocity of 0.1213 m/s, the boundary of the fluidized bed is lower than when using four groups of the solid phase. Particle separation in the upper part of the bed can be observed already for three gas velocities of 0.0892 m/s, 0.1088 m/s and 0.1213 m/s. For the cases of gas velocity of 0.0716 m/s and 0.0892 m/s, we see in the graphs a sufficient volume fraction of the four groups of the largest particle sizes. They form a limited fluidized bed. The smallest fraction is almost invisible. We saw the same result in Figure 10. For the case of gas velocity of 0.1088 m/s, we see that the three large fractions form the boundary of the bed. Particles of the smallest fraction are almost completely absent. The next fraction is presented mainly in the upper part of the bed and above the bed of large particles. For the case of a gas velocity of 0.1213 m/s, we see that two small groups of solid phase are carried out from the bed. As a result, the remaining fractions form a fluidized bed, the height of which is lower than for the case of using four groups of solids.

The removal of the finest group of particles from the apparatus for all values of gas velocity shows the disadvantage of the proposed model. One fraction, identified in the model, does not participate in the formation of the fluidized bed. However, such a situation may correspond to the real behavior of particles. We observed the removal of particles from the apparatus during experiments. But for numerical modeling, this is a negative point, since conservation equations will be solved for the fine fraction, which will increase computational costs.

By increasing the number of solid-phase fractions in a numerical simulation, we can obtain an average height of the fluidized bed close to the observed average experimental results. We can choose four solid-phase groups as a good result for our task. A further increase in the number of solid-phase groups could contribute to improving the accuracy of calculations, but it would require large computational resources. And the removal of some solid fractions from the apparatus leads to ineffective selection of some groups of particles. Further research is required to determine the size and number of small particle groups. It has been revealed that when selecting the mean diameter and volume fraction, it is imperative to consider the factors of gas and solid particle interaction, such as particle carryover and particle size separation by bed height.

4. Conclusions

This paper presents the investigation into the methodology employed to determine the mean particle size and the number of fractions for numerical simulation of a fluidized bed. We performed experimental studies of a bubbling fluidized bed with polydispersed powder on a glass tube. By selecting a different number of solid-phase groups, we can achieve a satisfactory agreement with the experimental results. We can more closely approach the real results by having a large number of groups of different sizes. A significant increase in computational costs will be caused by an increase in the number of fractions.

We decided to focus on the Sauter mean diameter model and selected the number of solid-phase groups as a control parameter. In this paper, we considered cases consisting of one fraction, as well as from two to five groups of solid-phase fractions. We have separated groups using only the particle size distribution function. This simple division does not require additional knowledge of the properties of the gas or solid phase. We found that the choice of four solid-phase groups is suitable for numerical simulation of a bubbling fluidized bed. Note that we have considered a powder sample to make an industrial catalyst. Our powder has a good unimodal particle size distribution with a pronounced high content of medium particles. Furthermore, each particle size distribution for a polydispersed powder may have its own characteristics that will contribute to the fluidization regime.

The proposed method for uniformly dividing the particle size distribution function is simple to use for numerical modeling. However, some particles may be carried away from the bed and not participate in the fluidization of the bed. This situation is shown in results when we divide the solid-phase fractions from three to five groups. In this situation, it is better to create methods that link the size of solid particles to the properties of the gas and solid phases. Similar problems will be considered in future studies in order to clarify the method considered in this study. Also, some particle ranges may have a very small volume fraction. Such groups make a minor contribution to the simulation result but take up additional computational resources. In this situation, it is better to create methods that relate the size of solid particles, the volume fraction of particle groups, and the properties of the gas and solid phases. Thus, to improve the efficiency, the proposed method of separating the distribution function requires further development. Such issues will be addressed in future studies to improve the method considered in this study.