Modeling of the Minimum Fluidization Velocity and the Incipient Fluidization Pressure Drop in a Conical Fluidized Bed with Negative Pressure

Abstract

The modeling of the minimum fluidization velocity (U_0mf) and the incipient fluidization pressure drop (ΔP_mf) is a valuable research topic in the fluidization field. In this paper, first, a series of experiments are carried out by changing the particle size and material mass to explore their effects on U_0mf and ΔP_mf. Then, an Ergun equation modifying method and the dimensional analysis method are used to obtain the modeling correlations of U_0mf and ΔP_mf by fitting the experimental data, and the advantages and disadvantages of the two methods are discussed. The experimental results show that U_0mf increases significantly with increasing particle size but has little relationship with the material mass; ΔP_mf increases significantly with increasing material mass but has little relationship with the particle size. Experiments with small particles show a significant increase at large superficial gas velocity; we propose a conjecture that the particles’ collision with the fluidization chamber’s top surface causes this phenomenon. The fitting accuracy of the modified Ergun equation is lower than that of the dimensionless model. When using the Ergun equation modifying method, it is deduced that the gas drag force is approximately 0.8995 times the material total weight at the incipient fluidized state.

1. Introduction

In the pharmaceutical industry, conical fluidized beds are more widely used than cylindrical fluidized beds because the former have certain advantages [1,2,3]. First, the particle size distribution can be broad in a granulation/coating process in the pharmaceutical industry. Due to the axial velocity gradient of a conical cylinder, larger particles can be fluidized at the lower layer, while smaller particles can stay at the upper layer without being carried out. Second, the fluidized state of a conical fluidized bed is smoother than that of a cylindrical fluidized bed because the upper layer of the conical cylinder is more extensive and can hold more gas to avoid the slug phenomenon [4]. Except for the unusual shape of the fluidization chamber, most pharmaceutical fluidized beds are negative pressure fluidized beds since the gas-driven equipment (usually a high-pressure fan) is placed at the outlet of the gas path, which is different from the fluidized beds used in other fields. The reason for this design is that if there is leakage, a negative pressure vessel is safer than a positive pressure vessel, and the material does not scatter everywhere, which would cause contamination.

The gas pressure is one of the essential monitored parameters in a fluidized bed [5]. Through analysis of the gas pressure signal, workers can know whether the fluidized bed is in a normal working state. The gas pressure drop (ΔP) between the surface and the bottom of the material is always used to determine the flow regimes. When the superficial gas velocity (U₀) is low (the subscript “0” in this paper refers to the bottom section of the conical cylinder), the bed is in a static state, and ΔP increases with an increase in U₀; when U₀ reaches a certain velocity, ΔP stops increasing and maintains a relatively stable value. This stable value is called the incipient fluidization pressure drop (ΔP_mf) [6], and the corresponding U₀ at this point is called the minimum fluidization velocity (U_0mf). U_0mf and ΔP_mf are of great importance in fluidized bed operation, and many factors affect these two values, among which particle size and material mass are fundamental and essential factors.

Because the experimental methods to obtain U_0mf and ΔP_mf under specified operating conditions are time-consuming and laborious, many scholars have tried to obtain modeling correlations suitable for general or specific working conditions [7]. In recent years, the proposed modeling correlations can be divided into two categories [8]: one is based on modifying the Ergun equation, and the other is based on dimensional analysis. The Ergun equation was proposed by Sabri Ergun in 1952 [9] to calculate the pressure drop when a gas flows through a packed material [10], and the specific equation is as follows:

$$\frac{dP}{dh}=C_1\frac{{(1-\epsilon )}^2}{{\epsilon }^3}\frac{\mu U_{\mathrm{h}}}{{\left(\varphi d_{\mathrm{P}}\right)}^2}+C_2\frac{(1-\epsilon )}{{\epsilon }^3}\frac{{\rho }_{\mathrm{g}}{U}_{\mathrm{h}}^2}{\varphi d_{\mathrm{P}}}$$

(1)

The Ergun equation has high accuracy due to the introduction of parameters such as particle size, particle sphericity, and voidage, so it has been widely accepted and verified. In this equation, C₁ is the coefficient of the viscosity term, and C₂ is the coefficient of the inertial term. For the original Ergun equation, C₁ = 150 and C₂ = 1.75. Both C₁ and C₂ are empirical values, and many factors, such as equipment size/shape, the friction coefficient, and particle size distribution, may lead to inaccurate prediction results [11]. In 1966, Wen and Yu [12] proposed an empirical model based on the Ergun equation and fitting the experimental data, the specific correlation of which was:

$$\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=\sqrt{{33.7}^2+0.0408\mathrm{Ar}}-33.7$$

(2)

In 1982, Grace J.R. [13] argued that replacing the number 33.7 with 27.2 in Equation (2) is most suitable for the gas-solid system. Gibson I.A. [14] et al. also did an in-depth discussion of Wen and Yu’s equation.

For a conical fluidized bed, the original Ergun equation cannot be directly used because of the gas velocity gradient in the axial direction, which makes the pressure drop rate vary in the axial direction. To solve this problem, Peng and Fan [15] developed a model to calculate ΔP of a conical liquid–solid fluidized bed by integrating the pressure drop rates of different layers. Jing et al. [16,17] proposed a correlation based on the Ergun equation to calculate U_0mf and ΔP_mf of a conical fluidized bed using a mathematical processing method similar to Peng and Fan’s model. Rachadaporn Kaewklum et al. [18] used Peng and Fan’s model to predict U_0mf and ΔP_mf for air–sand conical beds, and the effects of the sand particle size, cone angle, and static bed height on the flow behavior were discussed. Andrei Koekemoer et al. [19] put forward a method of modifying the coefficients C₁/C₂ in Equation (1) to increase the prediction accuracy, which is adopted and improved in this paper. In this study, we use regression analysis to try to obtain a better combination of C₁ and C₂ that is more suitable for the working conditions of this study.

Dimensional analysis is a widely used modeling method that only pays attention to the mathematical relationships between independent variables and dependent variables while ignoring the physical laws between physical parameters. A dimensionless model has a unified and straightforward equation form that is composed of a series of dimensionless groups. The general equation form is as follows:

$${\pi }_0={\pi }_1^{p_1}\cdot {\pi }_2^{p_2}\cdot {\pi }_3^{p_3}\cdots {\pi }_{\mathrm{n}}^{{\mathrm{p}}_{\mathrm{n}}}$$

(3)

where π₀, π₁, …, π_n are all dimensionless groups and P₁, P₂, …, P_n are the exponents of these dimensionless groups. Since the dimension of a dimensionless group is 1, no matter what the exponents are, the equation can satisfy the principle of dimensional harmony. The main difficulty of this method is determining the dimensionless groups, which requires knowing which factors affect the dependent variables. The value of each component can be determined by linear regression or nonlinear regression of the experimental data. Based on selecting appropriate dimensionless groups, a dimensionless model generally achieves higher accuracy under specific working conditions.

Sau et al. [20,21,22] used the dimensional analysis method to develop modeling correlations for U_0mf, ΔP_mf, and the bed expansion ratio of conical beds for different kinds of fluidization materials. Khani [23] proposed a dimensionless model using the term cos(α) instead of sin(α) compared with Sau’s model [20], considering that sin(α) approaches zero for mini-tapered beds, which may lead to an incorrect prediction result. Mojtaba Rasteh et al. [24] thoroughly considered the effect of PSD (particle size distribution) on U_0mf and proposed different correlations for different PSDs (Gaussian/flat/binary). Because the dimensional analysis method has strong applicability, it has been used by many scholars to derive modeling correlations in the field of fluidized bed [25,26,27,28,29,30,31]. Some of the correlations proposed by other scholars are listed in

Table 1. Correlations proposed by other scholars for predicting U_0mf and ΔP_mf.

Scholars	Correlations
Sau, D.C. et al. [20]	$\mathrm{Fr}=0.2714{(\mathrm{Ar})}^{0.3197}{(\sin \alpha )}^{0.6092}{\left(\frac{{\mathsf{\epsilon }}_0}{\varphi }\right)}^{-0.6108}$ $\Delta P_{\mathrm{mf}}=7.457{\left(\frac{D_1}{D_0}\right)}^{0.038}{\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.222}{\left(\frac{H_0}{D_0}\right)}^{0.642}{\left(\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{g}}}\right)}^{0.723}$
Khani, M.H. [23]	$\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=7.16{\left(\mathrm{Ar}\right)}^{0.393}{\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.987}{\left(\frac{{\epsilon }_0}{\varphi }\right)}^{-0.833}{(\cos \alpha )}^{-275.486},0\le \alpha \le 4.5°$ $\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=10.396{\left(\mathrm{Ar}\right)}^{0.367}{\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.889}{\left(\frac{{\epsilon }_0}{\varphi }\right)}^{-0.731}{(\cos \alpha )}^{-10.437},\alpha >4.5°$ $\frac{\Delta P_{\mathrm{mf}}}{{\rho }_{\mathrm{s}}g{H}_{\mathrm{s}}}=106.729(\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{g}}}{)}^{-0.522}{\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.309}{\left(\frac{H_0}{D_0}\right)}^{-0.379}{(\cos \alpha )}^{-10.858},0\le \alpha \le 4.5°$ $\frac{\Delta P_{\mathrm{mf}}}{{\rho }_{\mathrm{s}}g{H}_{\mathrm{s}}}=163.419(\frac{{\rho }_{\mathrm{s}}}{{\rho }_{\mathrm{g}}}{)}^{-0.524}{\left(\frac{d_{\mathrm{p}}}{{\mathrm{D}}_0}\right)}^{0.269}{\left(\frac{H_0}{D_0}\right)}^{-0.976}{(\cos \alpha )}^{-3.277},\alpha >4.5°$
Rasteh Mojtaba et al. [24]	$\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=0.254{\left(\mathrm{Ar}\right)}^{0.585}{\left(\frac{{\epsilon }_0}{\varphi }\right)}^{3.12}{\left(\frac{H_0}{D_0}\right)}^{0.264}{\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.0016}{(\cos \alpha )}^{-13.2}$

.

However, in the literature related to the modeling of U_0mf and ΔP_mf, most of the research equipment that was used was positive pressure cylindrical fluidized beds. Due to the differences in the gas path structure and the fluidization chamber shape, the previous research results cannot be directly applied to negative pressure conical fluidized beds. In addition, both the method based on modifying the Ergun equation and the method based on dimensional modeling have characteristics, advantages, and disadvantages, while few studies have compared these two methods. Based on these problems, we have carried out this study. The innovations of this study lie in the following: (a) the flow behaviors of the negative pressure conical fluidized bed are studied, (b) the Ergun equation modifying method and the dimensional analysis method are both used and compared, and (c) this study proposes the inference that the total gas drag force is 0.8995 times the particles’ gravity in the initial stage of fluidization. In this study, first, a series of experiments are carried out, and U_0mf and ΔP_mf are obtained as experimental empirical data; the effects of material mass and particle size on U_0mf and ΔP_mf are discussed. Second, we separately use the above two methods to fit the empirical data and derive models of U_0mf and ΔP_mf; by using these modeling correlations, workers can quickly determine U_0mf and ΔP_mf based on some variables, and then determine the production parameters. Finally, the prediction ability of the proposed correlations and the characteristics of the two modeling methods are discussed, which can be used as a reference for selecting an appropriate method in practical application.

2. Experimental Section

2.1. Experimental Equipment

As shown in Figure 1, in this study, a pharmaceutical granulation/drying fluidized bed (FBL-10, Xiaolun Pharmaceutical Machinery, Wenzhou, China) was used, which was a negative pressure conical fluidized bed. The conical cylinder was detachable with a bottom diameter of 180 mm and a top diameter of 300 mm, and the cross-sectional area was round; the conical angle α = 8.88°. There was a distribution plate made of stainless steel with a thickness of 3 mm for supporting materials and ventilation. The holes on the distribution plate were arranged in a square array mode, the diameter of each hole was 1.6 mm, and the total number of holes was 1617, so the overall opening ratio was 12.78%. A double-layer fine screen covered the distribution plate to prevent fine materials from falling through the holes and to force the airflow to be uniformly distributed. A high-pressure fan (2RB-730N-7AH26, Gzling Mechanical & Electrical Equipment, Wuxi, China) was placed at the outlet of the gas path, so the pressure in the fluidization chamber was negative with respect to the ambient pressure. The high-pressure fan had a rated power of 3 kW and could generate a maximum negative pressure of −22 kPa. The power of the fan was controlled by a frequency converter (ATV310HU30N4, Schneider Electric SA, Paris, France); thus, the gas velocity could be changed. The control range of the fan power was 1–100%, and the minimum control interval was 1%. The ambient pressure was standard atmospheric pressure, the air humidity was controlled by a humidity controller (MDH-40Y, Senjing Motor Manufacturing, Hangzhou, China) to 50–60%, and the ambient temperature was controlled via air conditioning to 25 °C ± 1 °C.

A vortex flowmeter (DN100, Meikong Automation Technology, Hangzhou, China) was connected in series at the exhaust pipe to measure the flowrate with an accuracy of 1%. Two air pressure measurement ports were set in the center of the plenum and the front of the filter bag, named P1 and P2, respectively. Two identical high-precision air pressure sensors were connected to the measurement ports, and the measurement range was −10 kPa–10 kPa with an accuracy of 0.5%. The pressure sensors’ signals were sampled by a high-speed data acquisition card (USB-6009, National Instruments, Austin, TX, USA) with a sampling frequency of 40 Hz. The pressure transfer hose (not shown in Figure 1) had a length of 100 mm and an inner diameter of 4 mm, so the transmission time delay of the pressure fluctuations in the hose could be ignored [32].

Micro glass beads were used as fluidized materials in this study, which are nonporous particles. The true density of these particles was 2409 ± 8.8 kg/m³, measured by the drainage method. Using different screen combinations, we screened out a total of 5 kinds of materials with different particle sizes, which were numbered #1, #2, #3, #4, and #5. The average diameter (randomly measuring the diameters of 50 particles and taking the average) and sphericity (replaced by the circularity of the projection photo) of each kind of particle were determined by a digital microscope (GP-660V, Gaopin Precision Instruments, Shanghai, China), and the natural bulk voidage of each kind of particle was measured by a bulk density meter (XF-16913, LICHEN-BX Instrument Technology, Shanghai, China), as shown in

Table 2. Physical properties of materials and Geldart particle type.

Material Number	#1	#2	#3	#4	#5
Screen size range (mm)	0.653–0.855	0.527–0.653	0.377–0.527	0.177–0.265	0.140–0.177
Particle diameter (mm)	0.744 ± 0.061	0.632 ± 0.046	0.454 ± 0.046	0.235 ± 0.025	0.164 ± 0.015
Bulk density (kg/m3)	1579 ± 2.1	1498 ± 2.7	1500 ± 2.3	1420 ± 0.5	1471 ± 3.3
Static voidage (-)	0.344 ± 1.7 × 10−3	0.378 ± 2.0 × 10−3	0.377 ± 1.9 × 10−3	0.410 ± 1.6 × 10−3	0.389 ± 2.2 × 10−3
Sphericity (-)	0.979 ± 0.031	0.966 ± 0.055	0.970 ± 0.058	0.903 ± 0.162	0.992 ± 0.037
Angle of repose (°)	24.1 ± 0.23	23.3 ± 0.15	26.2 ± 0.11	24.3 ± 0.30	23.8 ± 0.12
Hausner ratio (-)	1.011 ± 0.8 × 10−3	1.039 ± 1.1 × 10−3	1.035 ± 0.3 × 10−3	1.077 ± 0.8 × 10−3	1.058 ± 0.3 × 10−3
Material mass (kg)	2/2.5/3/3.5/4	2/2.5/3/3.5/4	2/2.5/3/3.5/4	2/2.5/3/3.5/4	2/2.5/3/3.5/4
Geldart particle type (-)	D	B	B	B	B

. Since the screen size range of each kind of particles was relatively narrow, the particle distributions of all kinds of particles were considered a typical narrow distribution. The microscopic photograph of the experimental particles is shown in Figure 2, and the bright spots on the particles in the photograph are caused by the microscope’s illumination light.

According to the Geldart classification method [33], #1 was class D particles, and #2–#5 were class B particles, where #2 was very close to the BD boundary. According to common knowledge in the field of fluidization, the main characteristic of class B and D particles is that there is no obvious expansion before the fluidization occurs. Additionally, the minimum bubbling velocity and the minimum fluidization velocity are the same; in other words, there is almost no particulate fluidization process in the fluidized bed.

2.2. Experimental Design and Procedure

Since port P1 is below the distribution plate, the pressure difference between ports P1 and P2 includes the pressure drop caused by the distribution plate. However, ΔP refers to the pressure drop between the top and bottom sections of the material, so it is necessary to obtain the pressure drop caused by the distributor first. For this reason, we first carried out a separate experiment in which there was no material in the fluidized bed, and the number of this experiment was 00.

In this study, a total of 5 kinds of material masses were set to 2/2.5/3/3.5/4 kg. The five kinds of material masses were cross-combined with the five kinds of particle sizes, and therefore, a total of 25 experiments were developed. The experimental numbers are given in

Table 3. Experimental number correspondence table.

Material Number	Experimental Number
	Material Mass 2 kg	Material Mass 2.5 kg	Material Mass 3 kg	Material Mass 3.5 kg	Material Mass 4 kg
#1	01	02	03	04	05
#2	06	07	08	09	10
#3	11	12	13	14	15
#4	16	17	18	19	20
#5	21	22	23	24	25

. For example, in experiment 14, #3 particles with a mass of 3.5 kg were used as the experimental material.

The experimental steps of experiments 01–25 were as follows: starting from a fan power opening (called fan opening hereafter) of 70%, decrease the fan opening by 1% each time. After the fluidized state is stable, record the signals (${\mathrm{x}}_{\mathrm{P}1}$, ${\mathrm{x}}_{\mathrm{P}2}$, ${\mathrm{x}}_{\mathrm{f}}$) for 30 s. Repeat the above operations until the fan power reaches 11%. Each experiment was repeated three times, and the most credible set of data was selected as the final data during data screening. To ensure that all experiments were carried out under the same conditions, before each experiment, the distribution plate, double-layer screen, and filter bag were cleaned, and the equipment was prefluidized at 50% fan power for 15 min to keep the material in a dry and loose state. For experiment 00, the fan power stopped at 0%; the rest of the procedure was the same as that in experiments 01–25.

3. Signal Processing Methods

This paper ignores the influence of gas compression and considers the gas flow rate at any cross section of the conical cylinder to be the same. The superficial gas velocity at any section is calculated using the following equations:

$$U_{\mathrm{h}}={\overline{x}}_{\mathrm{f}}/S_{\mathrm{h}}$$

(4)

$$U_0={\overline{x}}_{\mathrm{f}}/S_0$$

(5)

where ${\overline{x}}_{\mathrm{f}}$ is the average value of the flow rate time-series signal ($x_{\mathrm{f}}$) and S_h is the sectional area of the conical cylinder at the height of h.

The number of the empty bed experiment is 00, and the relationships among the fan opening, U₀, and ΔP_p are shown in Figure 3. As shown in Figure 3, with an increase in the fan opening, U₀ and ΔP_p continue to increase.

As shown in Equation (6), by subtracting the average values of $x_{\mathrm{P}1}$ and $x_{\mathrm{P}2}$ in the empty bed experiment (material mass W = 0), the pressure drop generated by the gas distribution plate (ΔP_p) can be calculated. For each fan opening, ΔP_p and U₀ were recorded so that the mapping relationship between a series of ΔP_p and a series of U₀ can be established. Therefore, for a given U₀, the value of ΔP_p can be determined by this mapping relationship. As shown in Equation (7), for experiments 01–25, by subtracting the average values of $x_{\mathrm{P}1}$ and $x_{\mathrm{P}2}$ and then subtracting ΔP_p (under the current U₀), ΔP can be obtained.

$$\mathsf{\Delta }{P}_{\mathrm{P}}=({\overline{x}}_{\mathrm{P}1}-{\overline{x}}_{\mathrm{P}2})\left|{}_{W=0}\right.(\mathrm{empty}\mathrm{bed}\mathrm{experiment},\mathrm{number}00)$$

(6)

$$\mathsf{\Delta }P={\overline{x}}_{\mathrm{P}1}-{\overline{x}}_{\mathrm{P}2}-\mathsf{\Delta }{P}_{\mathrm{p}}\left|{}_{U_0}\right.(\mathrm{experiment}01–25)$$

(7)

The intersection method is used to determine U_0mf and ΔP_mf. The specific determination process is as follows. first determine an inflection point, and from the second point before the inflection point, take as many points as possible near a straight line to fit the first straight line; then, from the second point after the inflection point, take 4–6 continuous points, then plot a horizontal straight line as the second line, whose height is the average of the taken points’ ordinate values; the intersection point of the two straight lines is considered to be the point (U_0mf, ΔP_mf). Take experiment 06 shown in Figure 4 as an example. The two black lines in Figure 4 are the two straight lines, and U_0mf and ΔP_mf are the abscissa and ordinate of the intersection point, respectively.

In this paper, the root mean square error (RMSE) and R² are used to describe the fitting accuracy quantitatively. RMSE means the root mean square error; the smaller the value, the closer the fitting value is to the actual value. R² is the determination coefficient, and its value is between 0 and 1. The higher the value is, the higher the fitting accuracy.

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{(x_{\mathrm{i}}-{\widehat{x}}_i)}^2}}$$

(8)

$${\mathrm{R}}^2=1-\frac{{\displaystyle \sum _{i=1}{(x_i-{\widehat{x}}_i)}^2}}{{\displaystyle \sum _{i=1}{(x_i-\overline{x})}^2}}$$

(9)

4. Results and Discussion

4.1. Effects of Material Mass and Particle Size on Flow Behavior

As shown in Figure 5, the 2nd row and 2nd column of the data in Table 3 are taken to display the effect of the material mass and the particle size on the flow behavior, and the experimental numbers are 06/07/08/09/10 and 02/07/12/17/22, respectively. A series of vertical dotted lines on the left side of Figure 5 represent U_0mf of each experiment, and their colors and markers are the same as those of the curves of the corresponding experiments.

As shown in Figure 5a, when U₀ < U_0mf, with an increase in U₀, ΔP continues to increase linearly; when U₀ > U_0mf, ΔP stops growing and enters a stable stage. In addition, the higher the material mass, the greater ΔP_mf is. The material mass has little effect on U_0mf; with an increase in the material mass, U_0mf increases slightly. As shown in Figure 5b, when U₀ < U_0mf, the particles with a larger particle size have a lower growth rate. When U₀ > U_0mf, the experiments have similar ΔP_mf values. When U₀ reaches a certain velocity, experiments 17/22 show significant ΔP growth, while the other experiments show little or no increase. We suppose the possible reason is as follows. When the gas velocity is high enough, the fine particles move to the upper part of the fluidized chamber and collide with the top surface. Assuming that all the particles are a whole body, this whole body is now subjected to the downward gravity force and the downward support force exerted by the top surface, so the pressure drop continues to increase because the required gas flow drag force becomes larger. For experiments with larger particle sizes (experiment 02/07/12), the particles cannot be brought high enough to collide with the top surface. Particle size has a significant influence on U_0mf; as the particle size increases, U_0mf increases significantly.

U_0mf and ΔP_mf obtained in all experiments are shown in

Table 4. U_0mf and ΔP_mf of each experiment.

Material Number	U0mf (m/s) ΔPmf (kPa)
	Material Mass 2 kg	Material Mass 2.5 kg	Material Mass 3 kg	Material Mass 3.5 kg	Material Mass 4 kg
#1	0.403 0.650	0.408 0.810	0.413 0.964	0.419 1.064	0.433 1.213
#2	0.302 0.640	0.328 0.793	0.344 0.933	0.361 1.093	0.375 1.213
#3	0.232 0.615	0.245 0.756	0.249 0.903	0.255 1.043	0.258 1.170
#4	0.106 0.643	0.125 0.771	0.136 0.899	0.149 1.020	0.157 1.116
#5	0.113 0.665	0.116 0.829	0.125 0.966	0.134 1.117	0.143 1.189

. 3D views of the particle size and material mass versus U_0mf and ΔP_mf are plotted to facilitate the observation of the variation regularity in Figure 6.

Table 4 and Figure 6 show that for U_0mf, the main influencing factor is the particle size; the larger the particle size, the larger U_0mf is. The influence of the material mass is small, and U_0mf increases slightly with an increase in the material mass. For ΔP_mf, the main influencing factor is the material mass, and ΔP_mf is approximately proportional to the material mass. The influence of the particle size is small and complex; for the same material mass, ΔP_mf of #1–#4 particles is almost the same, and ΔP_mf of #5 particles is 2–13% larger than the former.

4.2. Mathematical Modeling

4.2.1. Ergun Equation Modifying Method

The fluidized state could be observed through the glass window embedded in the sidewall of the conical cylinder, and the expansion of the bed surface could be clearly observed with the help of a thin film ruler that was pasted on the inner wall of the conical cylinder. When U₀ < U_0mf, the bed was in a fixed state. When U₀ approached U_0mf, slight expansion could be observed in the bed, and the bed height expansion ratio was approximately 2%. After U₀ reached U_0mf, bubbling was observed in the central part of the bed surface. As U₀ continued to increase, significant bubbling occurred throughout the entire cross section, and the fluidization height also increased as more particles were brought into the freeboard.

To simplify the model, the following assumptions are made for the point at which U₀ = U_0mf: the material expands uniformly, the bed height expansion ratio in the height direction is 2%, and the surface of the bed is flat.

From Equations (4) and (5), it can be further deduced that U_h is the function value of r_h, and the deduced equation is as follows:

$$U_{\mathrm{h}}=U_0\frac{S_0}{S_{\mathrm{h}}}=U_0\frac{r_0^2}{r_{\mathrm{h}}^2}$$

(10)

As Figure 7 shows, according to the trigonometric function equation:

$$r_{\mathrm{h}}=r_0+h\tan (\alpha )$$

(11)

$$r_1=r_0+H_0\tan (\alpha )$$

(12)

To obtain a precise H₀ value for each experiment in Table 3, the following equation is drawn based on the relationship of volume times density is equal to mass:

$$\frac{1}{3}{H}_0\pi (r_0^2+r_0{r}_1+r_1^2\left)\right(1-{\epsilon }_0){\rho }_{\mathrm{s}}=W$$

(13)

By solving the combined equations of Equations (12) and (13), the H₀ value can be obtained. In this study, an exhaustive method is used, and the programming and calculation are conducted in MATLAB. The specific method is to assign all the values of every 0.0001 in 0–0.1 to H₀ (unit: m) and substitute H₀ into Equations (12) and (13); then, calculate the equation error of Equation (13); finally, the best-matched H₀ value can be obtained by finding the minimum absolute equation error. The calculated H₀ values of all experiments are listed in Appendix A.

Since the bed expansion ratio is 2%, H_mf, r_mf, and ε_mf can be determined by the following equations:

$$H_{\mathrm{mf}}=1.02H_0$$

(14)

$$r_{\mathrm{mf}}=r_0+H_{\mathrm{mf}}\tan (\alpha )$$

(15)

$${\epsilon }_{\mathrm{mf}}=1-\frac{W}{{\rho }_{\mathrm{s}}\pi H_{\mathrm{mf}}(r_0^2+r_0{r}_{\mathrm{mf}}+r_{\mathrm{mf}}^2)/3}$$

(16)

Under the condition of U₀ = U_0mf, the original Ergun equation is as follows:

$$\frac{dP}{dh}=150\frac{{(1-{\epsilon }_{\mathrm{mf}})}^2}{{\epsilon }_{\mathrm{mf}}^3}\frac{\mu U_{\mathrm{h}}}{{\left(\varphi d_{\mathrm{P}}\right)}^2}+1.75\frac{(1-{\epsilon }_{\mathrm{mf}})}{{\epsilon }_{\mathrm{mf}}^3}\frac{{\rho }_{\mathrm{g}}{U}_{\mathrm{h}}^2}{\varphi d_{\mathrm{P}}}$$

(17)

Set $A=150\frac{{\left(1-{\epsilon }_{\mathrm{mf}}\right)}^2}{{\epsilon }_{\mathrm{mf}}^3}\frac{\mu }{{\left(\varphi d_{\mathrm{P}}\right)}^2}$ and $B=1.75\frac{(1-{\epsilon }_{\mathrm{mf}})}{{\epsilon }_{\mathrm{mf}}^3}\frac{{\rho }_{\mathrm{g}}}{\varphi d_{\mathrm{P}}}$ so that Equation (17) can be simplified to:

$$\frac{dP}{dh}=A{U}_{\mathrm{h}}+B{U}_{\mathrm{h}}^2$$

(18)

The drag force of the gas on the particles is:

$$dF=dP\cdot S_{\mathrm{h}}$$

(19)

Integrate the left and right sides of Equation (19) to h at the same time with an integral interval of 0-H_mf, and simplify it as follows:

$$\begin{array}{ll}{F}_{\mathrm{mf}}& ={\displaystyle {\int }_0^{H_{\mathrm{mf}}}\left(A{U}_{0\mathrm{mf}}\frac{r_0^2}{{\left(r_0+h\tan \alpha \right)}^2}+B{U}_{0\mathrm{mf}}^2\frac{r_0^4}{{\left(r_0+h\tan \alpha \right)}^4}\right)\pi {\left(r_0+h\tan \alpha \right)}^2\mathrm{dh}}\\ & =\pi r_0^2{H}_{\mathrm{mf}}\left(A{U}_{0\mathrm{mf}}+B{U}_{0\mathrm{mf}}^2\frac{r_0}{r_{\mathrm{mf}}}\right)\end{array}$$

(20)

Since the inner wall of the conical cylinder is relatively smooth, the friction force of the wall against the particles can be ignored. According to previous studies [15,16], when U₀ = U_0mf, the drag force exerted on the particles is equal to the gravity force of the particles; that is:

$$F_{\mathrm{mf}}=G$$

(21)

Substitute Equation (20) into Equation (21); that is:

$$\pi r_0^2{H}_{\mathrm{mf}}(A{U}_{0\mathrm{mf}}+B{U}_{0\mathrm{mf}}^2\frac{r_0}{r_{\mathrm{mf}}})=Wg$$

(22)

By solving Equation (22), the calculated U_0mf can be obtained. (H_mf, r_mf, and ε_mf can be obtained in advance by Equations (14)–(16), respectively)

Integrate the left and right sides of Equation (18) to h at the same time with an integral interval of 0-H_mf, and simplify it as follows:

$$\begin{array}{ll}\Delta P_{\mathrm{mf}}& ={\displaystyle {\int }_0^{H_{\mathrm{mf}}}\left(A{U}_{0\mathrm{mf}}\frac{r_0^2}{{\left(r_0+h\tan \alpha \right)}^2}+B{U}_{0\mathrm{mf}}^2\frac{r_0^4}{{\left(r_0+h\tan \alpha \right)}^4}\right)dh}\\ & =A{U}_{0\mathrm{mf}}{H}_{\mathrm{mf}}\frac{r_0}{r_{\mathrm{mf}}}+B{U}_{0\mathrm{mf}}^2{H}_{\mathrm{mf}}\frac{r_0(r_0^2+r_0{r}_{\mathrm{mf}}+r_{\mathrm{mf}}^2)}{3r_{\mathrm{mf}}^3}\end{array}$$

(23)

The calculated ΔP_mf can be obtained by substituting the calculated U_0mf into Equation (23).

For every experiment in Table 3, conduct the above calculation steps to obtain the calculated U_0mf and calculated ΔP_mf values. To facilitate observation, the experimentally measured values of U_0mf and ΔP_mf are taken as the abscissa and the calculated values as the ordinate, as Figure 8 shows.

Figure 8 shows that the prediction accuracy of U_0mf is relatively high (absolute error < 20%) for the data of experiments 06–10 and very poor for experiments 16–25 in which small particles were used, and all the prediction values are smaller than the experimental values. For the prediction of ΔP_mf, the prediction accuracy is high, with similar errors between 0–+20%. Substituting all prediction and experimental values into Equations (8) and (9), the overall RMSE and R² can be obtained, as shown in

Table 5. Coefficients and prediction accuracies of the original and modified Ergun equations.

	C1	C2	C3	U0mf		ΔPmf
				RMSE (m/s)	R2 (-)	RMSE (kPa)	R2 (-)
Original	150	1.75	1	0.0872	0.4209	0.1076	0.6908
Modified	46.61	3.25	0.8995	0.0402	0.8771	0.0291	0.9774

. There are many possible reasons for the low prediction accuracy: (a) The measured values of some variables, especially the sphericity, are difficult to accurately determine. (b) There may be differences in some environmental parameters compared with the experiments for deducing the Ergun equation, such as gas humidity, which greatly affects the gas viscosity and the drag force. (c) The conical wall may have an upward support force for the material in the initial fluidization stage.

To increase the prediction accuracy of U_0mf and ΔP_mf, the original coefficients at two sites are changed. The coefficients C₁/C₂ in Equation (1) are changed at the first site, whose original values are 150/1.75, respectively. At the second site, considering that the conical wall may support the material, a multiplicative correction coefficient C₃ is added to Equation (21) as follows:

$$F_{\mathrm{mf}}=C_3\cdot G$$

(24)

Therefore, Equation (22) is modified as follows:

$$\pi r_0^2{H}_{\mathrm{mf}}(A{U}_{0\mathrm{mf}}+B{U}_{0\mathrm{mf}}^2\frac{r_0}{r_{\mathrm{mf}}})=C_3\cdot Wg$$

(25)

where $A=C_1\frac{{\left(1-{\epsilon }_{\mathrm{mf}}\right)}^2}{{\epsilon }_{\mathrm{mf}}^3}\frac{\mu }{{\left(\varphi d_{\mathrm{P}}\right)}^2}$ and $B=C_2\frac{\left(1-{\epsilon }_{\mathrm{mf}}\right)}{{\epsilon }_{\mathrm{mf}}^3}\frac{{\rho }_{\mathrm{g}}}{\varphi d_{\mathrm{P}}}$. Equation (25) is the equation used in the following fitting operations.

All the calculations are performed by MATLAB software, a method of exhaustive optimization is used, and the optimization goal is to minimize the sum of the RMSE values of U_0mf and ΔP_mf (RMSE_sum).

The specific calculation processes are shown in Figure 9. First, a preset value of 0.8 is given for C₃, then C₁/C₂ are obtained by fitting the experimental data using a linear regression method (“regress” function in MATLAB). Then, this C₁/C₂/C₃ combination is used to calculate the value of RMSE_sum, and the C₁/C₂/C₃ combination and RMSE_sum are recorded. Then, add 0.0001 to the previous C₃ preset value, and perform the same calculation and recording for each C₃ preset value. When the C₃ preset value reaches 1.2, the exhaustive optimization calculation is ended. Finally, the C₁/C₂/C₃ combination corresponding to the minimum value of RMSE_sum is found. The final results are C₁ = 46.61, C₂ = 3.25, and C₃ = 0.8995. The coefficients and modeling accuracies of the original and modified Ergun equations are listed in Table 5, and the modeling results of the modified Ergun equation are shown in Figure 10.

Substitute C₃ = 0.8995 into Equation (24), i.e.,

$$F_{\mathrm{mf}}=0.8995\cdot G$$

(26)

We speculate that the reason for Equation (26) is that when fluidization occurs, the conical wall still has support for the particles, which can decompose into an upward component force; therefore, when the drag force of the gas is less than the gravity force of particles (0.8995 < 1), the particles can still be fluidized, and the coefficient C₃ should be related to the cone angle. Some literature has proposed separation models of conical fluidized beds, which may help analyze this problem [34,35]. In future research, we will try to measure the possible supporting force and study the value trend of the coefficient C₃.

To sum up, the specific calculation process of the Ergun equation modifying method is that U_0mf can be calculated by solving the following equation:

$$\pi r_0^2{H}_{\mathrm{mf}}(A{U}_{0\mathrm{mf}}+B{U}_{0\mathrm{mf}}^2\frac{r_0}{r_{\mathrm{mf}}})=0.8995\cdot Wg$$

(27)

where $A=46.6\frac{{(1-{\epsilon }_{\mathrm{mf}})}^2}{{\epsilon }_{\mathrm{mf}}^3}\frac{\mu }{{\left(\varphi d_{\mathrm{P}}\right)}^2}$ and $B=3.25\frac{(1-{\epsilon }_{\mathrm{mf}})}{{\epsilon }_{\mathrm{mf}}^3}\frac{{\rho }_{\mathrm{g}}}{\varphi d_{\mathrm{P}}}$ (H_mf, r_mf, and ε_mf can be obtained by Equations (14)–(16), respectively). After U_0mf is calculated, ΔP_mf can be calculated by substituting the calculated U_0mf into Equation (23).

4.2.2. Dimensional Analysis

To obtain a dimensionless correlation, we first determine the factors that influence the hydrodynamic behavior of the fluidized bed, which are the particle size (d_p), initial bed height (H₀), bottom diameter (D₀), initial bed voidage (ε₀), particle sphericity (φ), particle true density (ρ_s), gas density (ρ_g), gas viscosity (μ), angle of repose (θ), Hausner ratio (R). These factors are recast into several dimensionless groups containing the Reynolds number ($\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=\frac{{\rho }_{\mathrm{g}}{U}_{0\mathrm{mf}}{d}_{\mathrm{p}}}{\mu }$), Archimedes number ($\mathrm{Ar}=\frac{g{d}_{\mathrm{p}}^3{\rho }_{\mathrm{g}}({\rho }_{\mathrm{s}}-{\rho }_{\mathrm{g}})}{{\mu }^2}$), $\frac{\Delta P_{\mathrm{mf}}}{{\rho }_{\mathrm{s}}g{H}_0}$, $\frac{{\epsilon }_0}{\varphi }$, $\frac{H_0}{D_0}$, $\frac{d_{\mathrm{p}}}{D_0}$, $\frac{\theta }{R}$, and the correlations among these above dimensionless groups with undetermined coefficients are as follows:

$$\mathrm{R}{\mathrm{e}}_{\mathrm{mf}}=a_1\cdot {\left(\mathrm{Ar}\right)}^{b_1}\cdot {\left(\frac{{\epsilon }_0}{\varphi }\right)}^{c_1}\cdot {\left(\frac{H_0}{D_0}\right)}^{d_1}\cdot {\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{e_1}\cdot {\left(\frac{\theta }{R}\right)}^{f_1}$$

(28)

$$\frac{\Delta P_{\mathrm{mf}}}{{\rho }_{\mathrm{s}}g{H}_0}=a_2\cdot {\left(\mathrm{Ar}\right)}^{b_2}\cdot {\left(\frac{{\epsilon }_0}{\varphi }\right)}^{c_2}\cdot {\left(\frac{H_0}{D_0}\right)}^{d_2}\cdot {\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{e_2}\cdot {\left(\frac{\theta }{R}\right)}^{f_2}$$

(29)

To obtain the undetermined coefficients of Equations (28) and (29), both linear regression and nonlinear regression methods can be used. For the linear regression method, first, take the logarithm of both sides of Equation (28) as follows:

$$\lg (\mathrm{R}{\mathrm{e}}_{\mathrm{mf}})=\lg (a_1)+b_1\cdot \lg (\mathrm{Ar})+c_1\cdot \lg \left(\frac{{\epsilon }_0}{\varphi }\right)+d_1\cdot \lg \left(\frac{H_0}{D_0}\right)+e_1\cdot \lg \left(\frac{d_{\mathrm{p}}}{D_0}\right)+f_1\cdot \lg \left(\frac{\theta }{R}\right)$$

(30)

Then, the “regress” function in MATLAB can be used to fit the undetermined coefficients, i.e., a₁, b₁, c₁, d₁, e_1, and f₁

For the nonlinear regression method, the “nlinfit” function in MATLAB is used to fit the undetermined coefficients of Equations (28) and (29) directly, and the results of the linear regression are set as the initial values of iteration. The results of the linear and nonlinear regression methods are shown in

Table 6. Results of the linear and nonlinear regression methods of U_0mf for the dimensionless model.

	a1	b1	c1	d1	e1	f1	RMSE (m/s)	R2 (-)
Linear regression	1	0.504	−1.355	0.326	0.236	−0.635	0.011	0.991
Nonlinear regression	1.033	0.509	−1.142	0.186	0.232	−0.645	0.009	0.994

and

Table 7. Results of the linear and nonlinear regression methods of ΔP_mf for the dimensionless model.

	a2	b2	c2	d2	e2	f2	RMSE (kPa)	R2 (-)
Linear regression	1	−0.003	−0.642	0.020	0.010	−0.354	0.015	0.994
Nonlinear regression	0.997	−0.003	−0.648	0.017	0.010	−0.352	0.015	0.994

and Figure 11. Based on the values of R² and RMSE, there is little difference between the fitting accuracy of the linear and nonlinear regression methods, especially for ΔP_mf, where the models are almost the same.

Since the linear regression has better convergence ability than the nonlinear regression and does not need to assign initial values, it is more practical. Therefore, we choose the results of the linear regression method, and the final dimensionless correlation to predict U_0mf and ΔP_mf is as follows:

$$U_{0\mathrm{mf}}={\left(\mathrm{Ar}\right)}^{0.504}\cdot {\left(\frac{{\epsilon }_0}{\varphi }\right)}^{-1.355}\cdot {\left(\frac{H_0}{D_0}\right)}^{0.326}\cdot {\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.236}\cdot {\left(\frac{\theta }{R}\right)}^{-0.635}\cdot \frac{\mu }{d_{\mathrm{p}}{\rho }_{\mathrm{g}}}$$

(31)

$$\mathsf{\Delta }{P}_{\mathrm{mf}}={\left(\mathrm{Ar}\right)}^{-0.003}\cdot {\left(\frac{{\epsilon }_0}{\varphi }\right)}^{-0.642}\cdot {\left(\frac{H_0}{D_0}\right)}^{0.020}\cdot {\left(\frac{d_{\mathrm{p}}}{D_0}\right)}^{0.010}\cdot {\left(\frac{\theta }{R}\right)}^{-0.354}\cdot {\rho }_{\mathrm{s}}g{H}_0$$

(32)

4.2.3. Discussion of the Two Modeling Methods

Based on the values of R² and RMSE of the two modeling methods (Table 5, Table 6 and Table 7), the dimensional analysis method is better than the Ergun equation modifying method. This result mainly occurs because the dimensionless correlation contains six undetermined coefficients, which is more than the three undetermined coefficients in the modified Ergun equation. More undetermined coefficients can effectively improve the fitting accuracy but also increase the complexity of the fitting calculation.

The physical laws between physical quantities are considered in the Ergun equation modifying method, so its correlation has relatively strong logicality. When using this method, by analyzing and deducing the formulas, the current fluidized state can be analyzed more deeply, and even some other inferences or conclusions can be obtained. For example, the inference of F_mf = 0.8995G is obtained in this study, which offers some new information on analyzing the fluidized state. The dimensional analysis only fits the relationship between physical quantities from the numerical point of view. Under the condition of the same equipment and operation parameters, it has high fitting accuracy. However, the adaptability of this method is inadequate. For example, Equations (30) and (31) have higher prediction accuracy for the negative pressure conical fluidized bed with a glass bead material, but they may not be able to predict accurately for other types of fluidized beds or materials. These two methods have their advantages and disadvantages, and the appropriate method can be chosen according to the needs of actual production.

5. Conclusions

The research equipment in this paper is a pharmaceutical conical fluidized bed with negative pressure. First, the relationship between particle size/material mass and U_0mf/ΔP_mf is explored. Then, an Ergun equation modifying method and the dimensional analysis method are used to establish prediction models for U_0mf/ΔP_mf, and linear and nonlinear regression analysis methods are both used to determine the coefficients in the models. Finally, the advantages and disadvantages and applicable scope of the two modeling methods are compared and evaluated.

The conclusions of this study can be summarized as follows:

For U_0mf, the main influencing factor is the particle size, and the larger the particle size, the larger the U_0mf is; the influence of the material mass is small. For ΔP_mf, the main influencing factor is the material mass, and ΔP_mf is approximately proportional to the material mass; the influence of the particle size is relatively small. Experiments with small particles show a significant increase at large superficial gas velocity; we propose a conjecture that the particles’ collision with the fluidization chamber’s top surface causes this phenomenon.

The Ergun equation modifying method takes into account the physical laws between physical quantities. The original Ergun equation has low prediction accuracy for U_0mf and ΔP_mf for this study. Modifying the coefficients of the Ergun equation can significantly increase the prediction accuracy. By analyzing the data, we draw the inference of F_mf = 0.8995G, which may be due to the upward supporting effect of the conical cylinder wall on the material and needs to be further explored and verified.

The dimensional analysis method considers only the mathematical relationships between independent variables and dependent variables. Compared with the Ergun equation modifying method, the dimensional analysis method has a higher fitting accuracy. There is little difference in the fitting accuracy between the linear regression method and the nonlinear regression method.