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

: The modeling of the minimum ﬂuidization velocity ( U 0mf ) and the incipient ﬂuidization pressure drop ( ∆ P mf ) is a valuable research topic in the ﬂuidization ﬁeld. In this paper, ﬁrst, a series of experiments are carried out by changing the particle size and material mass to explore their e ﬀ ects 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 ﬁtting the experimental data, and the advantages and disadvantages of the two methods are discussed. The experimental results show that U 0mf increases signiﬁcantly with increasing particle size but has little relationship with the material mass; ∆ P mf increases signiﬁcantly with increasing material mass but has little relationship with the particle size. Experiments with small particles show a signiﬁcant increase at large superﬁcial gas velocity; we propose a conjecture that the particles’ collision with the ﬂuidization chamber’s top surface causes this phenomenon. The ﬁtting accuracy of the modiﬁed 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 ﬂuidized state. method and the nonlinear regression method.


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 0 ) 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 0 ; when U 0 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 0 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: dP 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 1 is the coefficient of the viscosity term, and C 2 is the coefficient of the inertial term. For the original Ergun equation, C 1 = 150 and C 2 = 1.75. Both C 1 and C 2 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: Re mf = 33.7 2 + 0.0408Ar − 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 1 /C 2 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 1 and C 2 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 Appl. Sci. 2020, 10, 8764 3 of 19 equation form that is composed of a series of dimensionless groups. The general equation form is as follows: where π 0 , π 1 , . . . , π n are all dimensionless groups and P 1 , P 2 , . . . , 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. 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.

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 3 , 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. 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.

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 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.

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. 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 (x P1 , x P2 , x 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.

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: where x f is the average value of the flow rate time-series signal (x 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 0 , and ∆P p are shown in Figure 3. As shown in Figure 3, with an increase in the fan opening, U 0 and ∆P p continue to increase.
where x ̅ f is the average value of the flow rate time-series signal (x f ) and Sh 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, U0, and ΔPp are shown in Figure 3. As shown in Figure 3, with an increase in the fan opening, U0 and ΔPp continue to increase. As shown in Equation (6), by subtracting the average values of x P1 and x P2 in the empty bed experiment (material mass W = 0), the pressure drop generated by the gas distribution plate (ΔPp) can be calculated. For each fan opening, ΔPp and U0 were recorded so that the mapping relationship between a series of ΔPp and a series of U0 can be established. Therefore, for a given U0, the value of ΔPp can be determined by this mapping relationship. As shown in Equation (7), for experiments 01-25, by subtracting the average values of x P1 and x P2 and then subtracting ΔPp (under the current U0), ΔP can be obtained.
The intersection method is used to determine U0mf and ΔPmf. 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 (U0mf, ΔPmf). Take experiment 06 shown in Figure  4 as an example. The two black lines in Figure 4 are the two straight lines, and U0mf and ΔPmf are the abscissa and ordinate of the intersection point, respectively. As shown in Equation (6), by subtracting the average values of x P1 and x P2 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 0 were recorded so that the mapping relationship between a series of ∆P p and a series of U 0 can be established. Therefore, for a given U 0 , 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 P1 and x P2 and then subtracting ∆P p (under the current U 0 ), ∆P can be obtained.
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 2 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 2 is the determination coefficient, and its value is between 0 and 1. The higher the value is, the higher the fitting accuracy. In this paper, the root mean square error (RMSE) and R 2 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 2 is the determination coefficient, and its value is between 0 and 1. The higher the value is, the higher the fitting accuracy.

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.  Table 3.
U0mf and ΔPmf obtained in all experiments are shown in Table 4. 3D views of the particle size and material mass versus U0mf and ΔPmf are plotted to facilitate the observation of the variation regularity in Figure 6.   Table 3.
As shown in Figure 5a, when U 0 < U 0mf , with an increase in U 0 , ∆P continues to increase linearly; when U 0 > 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 0 < U 0mf , the particles with a larger particle size have a lower growth rate. When U 0 > U 0mf , the experiments have similar ∆P mf values. When U 0 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. 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. U0mf and ΔPmf obtained in all experiments are shown in Table 4. 3D views of the particle size and material mass versus U0mf and ΔPmf are plotted to facilitate the observation of the variation regularity in Figure 6.    Figure 6 show that for U0mf, the main influencing factor is the particle size; the larger the particle size, the larger U0mf is. The influence of the material mass is small, and U0mf increases slightly with an increase in the material mass. For ΔPmf, the main influencing factor is the material mass, and ΔPmf is approximately proportional to the material mass. The influence of the particle size   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.

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 0 < U 0mf , the bed was in a fixed state. When U 0 approached U 0mf , slight expansion could be observed in the bed, and the bed height expansion ratio was approximately 2%. After U 0 reached U 0mf , bubbling was observed in the central part of the bed surface. As U 0 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 0 = 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. (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:

From Equations
As Figure 7 shows, according to the trigonometric function equation: Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 19 Under the condition of U0 = U0mf, the original Ergun equation is as follows: so that Equation (17) can be simplified to: The drag force of the gas on the particles is: Integrate the left and right sides of Equation (19) to h at the same time with an integral interval of 0-Hmf, and simplify it as follows: 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 U0 = U0mf, the drag force exerted on the particles is equal to the gravity force of the particles; that is: Substitute Equation (20) into Equation (21); that is: To obtain a precise H 0 value for each experiment in Table 3, the following equation is drawn based on the relationship of volume times density is equal to mass: By solving the combined equations of Equations (12) and (13), the H 0 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 0 (unit: m) and substitute H 0 into Equations (12)  Since the bed expansion ratio is 2%, H mf , r mf , and ε mf can be determined by the following equations: Under the condition of U 0 = U 0mf , the original Ergun equation is as follows: so that Equation (17) can be simplified to: The drag force of the gas on the particles is: 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: 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 0 = U 0mf , the drag force exerted on the particles is equal to the gravity force of the particles; that is: Substitute Equation (20) into Equation (21); that is: 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: 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][17][18][19][20][21][22][23][24][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 2 can be obtained, as shown in Table 5. 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. r r r + r r + r = AU H + BU H r r The calculated ΔPmf can be obtained by substituting the calculated U0mf into Equation (23). For every experiment in Table 3, conduct the above calculation steps to obtain the calculated U0mf and calculated ΔPmf values. To facilitate observation, the experimentally measured values of U0mf and ΔPmf are taken as the abscissa and the calculated values as the ordinate, as Figure 8 shows.  Figure 8 shows that the prediction accuracy of U0mf is relatively high (absolute error < 20%) for the data of experiments 06-10 and very poor for experiments [16][17][18][19][20][21][22][23][24][25] in which small particles were used, and all the prediction values are smaller than the experimental values. For the prediction of ΔPmf, 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 2 can be obtained, as shown in Table 5. 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 U0mf and ΔPmf, the original coefficients at two sites are changed. The coefficients C1/C2 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 C3 is added to Equation (21) as follows: Therefore, Equation (22) is modified as follows:  To increase the prediction accuracy of U 0mf and ∆P mf , the original coefficients at two sites are changed. The coefficients C 1 /C 2 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 3 is added to Equation (21) as follows: Therefore, Equation (22) is modified as follows: where 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 3 , then C 1 /C 2 are obtained by fitting the experimental data using a linear regression method ("regress" function in MATLAB). Then, this C 1 /C 2 /C 3 combination is used to calculate the value of RMSE sum , and the C 1 /C 2 /C 3 combination and RMSE sum are recorded. Then, add 0.0001 to the previous C 3 preset value, and perform the same calculation and recording for each C 3 preset value. When the C 3 preset value reaches 1.2, the exhaustive optimization calculation is ended. Finally, the C 1 /C 2 /C 3 combination corresponding to the minimum value of RMSE sum is found. The final results are C 1 = 46.61, C 2 = 3.25, and C 3 = 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.
C3, then C1/C2 are obtained by fitting the experimental data using a linear regression method ("regress" function in MATLAB). Then, this C1/C2/C3 combination is used to calculate the value of RMSEsum, and the C1/C2/C3 combination and RMSEsum are recorded. Then, add 0.0001 to the previous C3 preset value, and perform the same calculation and recording for each C3 preset value. When the C3 preset value reaches 1.2, the exhaustive optimization calculation is ended. Finally, the C1/C2/C3 combination corresponding to the minimum value of RMSEsum is found. The final results are C1 = 46.61, C2 = 3.25, and C3 = 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.  Figure 9. Specific calculation processes to determine the coefficients C 1 /C 2 /C 3.
("regress" function in MATLAB). Then, this C1/C2/C3 combination is used to calculate the value of RMSEsum, and the C1/C2/C3 combination and RMSEsum are recorded. Then, add 0.0001 to the previous C3 preset value, and perform the same calculation and recording for each C3 preset value. When the C3 preset value reaches 1.2, the exhaustive optimization calculation is ended. Finally, the C1/C2/C3 combination corresponding to the minimum value of RMSEsum is found. The final results are C1 = 46.61, C2 = 3.25, and C3 = 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.
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 3 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 3 .
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: where A =46.6 (1 −ε mf ) 2 ε 3 mf µ (φd P ) 2 and B =3.25 (1 −ε mf ) ε 3 mf ρ g φd 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).

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 0 ), bottom diameter (D 0 ), initial bed voidage (ε 0 ), 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 (Re mf = ρ g U 0mf d p µ ), Archimedes number , and the correlations among these above dimensionless groups with undetermined coefficients are as follows: 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: Then, the "regress" function in MATLAB can be used to fit the undetermined coefficients, i.e., a 1 , b 1 , c 1 , d 1 , e 1, and f 1 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 Tables 6 and 7 and Figure 11. Based on the values of R 2 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 U0mf and ΔPmf is as follows: 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.

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.