Effect of Drag Models on Residence Time Distributions of Particles in a Wurster Fluidized Bed: a DEM-CFD Study †

Fluidized bed coating has been used to coat pellets or tablets with functional substances for a number of purposes. In this coating process, particle wetting, drying and film formation are coupled to particle motion. It is therefore of interest to study particle motion in such fluidized beds and to use the results to develop a model for predicting the quality of the final product. In this paper, we present results from DEM-CFD simulations, i.e. discrete element method and computational fluid dynamics simulations of particle motion in a laboratory-scale Wurster fluidized bed that was also employed in positron emission particle tracking (PEPT) experiments. As the drag force is the dominant interaction between the gas flow and the particle motion in this type of fluidized bed, the effect of drag models on the particle motion is investigated. More specifically, the particle velocity and residence time distributions of particles in different regions calculated from five different drag models are presented. It is found that the Gidaspow and Tang drag models predict both particle cycle and residence times well. The HKL and Beetstra drag models somewhat overestimate the particle velocity in the Wurster tube and therefore predict a reduced number of recirculations and a significantly shorter cycle time.


Introduction
The Wurster coating process (Wurster D.E., 1959), which takes place in a fluidized bed, has been utilized for pellet coating for several decades. In this process, particles are circulated in the fluidized bed and are coated in a spray zone for different functional needs (Teunou E. and Poncelet D., 2002;Ström D. et al., 2005). At the bottom of the fluidized bed, a fluidization air flow is supplied through a distributor plate. The distributor plate usually consists of a number of orifices that have been distributed so as to provide a specific flow distribution. An atomization air flow and a liquid solution in the form of droplets are introduced by one or more two-fluid nozzles at the bottom of the fluidized bed.
This type of fluidized bed can be divided into different regions (Christensen F.N. and Bertelsen P., 1997;Karlsson S. et al., 2009), as illustrated in Fig. 1. In an ideal coating cycle, a particle begins to receive coating in the spray zone, is then dried in the Wurster tube and the fountain region, and moves back towards the horizontal transport region through the downbed region. After travelling through these different regions, the particle is ready to begin the next coating cycle when it moves into the spray zone again. The sequence of coating and drying that occurs as the particles travel through these different regions is usually regarded as a particle coating cycle (Li L. et al., 2015a). It has been pointed out that the cycle time distribution (CTD) and the residence time distribution (RTD) in different regions are critical factors in determining the quality of the final product such as the film thickness and the film thickness variability (Cheng X.X. and Turton R., 2000;Shelukar S. et al., 2000;Li L. et al., 2015a). For example, the particles must spend sufficient time in the regions for drying before returning to the horizontal transport region; otherwise, the risk of agglomeration increases due to collisions between wet particles.
The CTD and RTD of particles in similar fluidized beds have been studied using both experimental (Mann V. and Crosby E., 1975;Shelukar S. et al., 2000;San José M.J. et al., 2013;Li L. et al., 2015a) and modeling (Fries L. et al., 2011;Yang S. et al., 2014;Li L. et al., 2015b) methods. It was found that particles spend most of the cycle time in the horizontal transport region and that different residence times in this region can cause a large variability in the cycle time, and therefore a large variability in the film thickness. Due to limitations of experimental techniques, it is difficult to measure how particles behave in detail in the spray zone. For example, in the positron emission particle tracking (PEPT) measurements (Parker D.J. et al., 1993), the γ-rays tend to be absorbed by the metal in the lower region of the fluidized bed, which results in poor data quality in the spray region (Li L. et al., 2015a). However, the particle behavior in the spray zone is important for coating. Thus advanced modeling methods, such as the discrete element method (DEM), in combination with computational fluid dynamics (CFD), offer an attractive way of addressing this kind of issue that cannot be handled experimentally. The discrete element method was initially proposed for granular assemblies by Cundall P.A. and Strack O.D. (1979), and was introduced into particle flow in fluidized beds by Tsuji Y. et al. (1993). In this approach, every particle is modeled as an individual element and the gas phase is treated as a continuum. Particle-particle and particle-wall interactions can be accounted for in a straightforward fashion and it is also simple to include particles of different sizes. Owing to this benefit, a number of reviews and applications of this method have been presented in the literature (Kafui K.D. et al., 2002;Deen N.G. et al., 2007;Van der Hoef M.A. et al., 2008;Link J.M. et al., 2008;Zhu H.P. et al., 2008;Ketterhagen W.R. et al., 2009;Turton R., 2010;Darabi P. et al., 2011;Fries L. et al., 2013).
On the other hand, since every particle is modeled, the discrete element method is not practical for systems involving a considerable number of particles. For many particle coating processes in pharmaceutical development, the systems often have a batch size ranging from a few grams to a few hundred grams, which corresponds to from tens of thousands to a few million particles. Nowadays, simulating tens of thousands particles can be done at a reasonable computational cost, while modeling systems with a few million particles is challenging, though still feasible (Jajcevic D. et al., 2013). Hence, DEM-CFD modeling offers a great opportunity to assist process and product development.
In DEM-CFD simulations, particles and gas flow are coupled through an interphase momentum transfer term via a suitable drag model. It is well known that in fluidized beds, particles interact strongly with the gas flow, and consequently the drag model implemented affects the predicted particle motion, as well as the CTD and RTD. In the assessment of drag models by Du W. et al. (2006) based on the two-fluid model (TFM), the Gidaspow drag model (Gidaspow D., 1994) gave the best agreement with experimental observations in modeling of spouted beds. With the recent increase in computational power, it has become feasible to derive the drag model by fully solving the flow between particles (Kriebitzsch S.H.L. et al., 2013). New drag models have been proposed and developed in recent years based upon simulations using the Lattice-Boltzmann method (LBM) or immersed boundary method (IBM) (Hill R.J. et al., 2001a;2001b;Van der Hoef M.A. et al., 2005;Benyahia S. et al., 2006;Beetstra R. et al., 2007;Tang Y. et al., 2015). In a combined experimental and simulation study, Link J.M. et al. (2005) evaluated the Gidaspow and Hill drag models and found that they could reproduce several important regimes in spoutfluid beds. At this moment, however, it is still unclear how the drag model affects the residence time distributions of particles in different regions of fluidized beds.
In the present study, the objective is to investigate the effect of drag models on the residence time distributions of particles in different regions of a Wurster fluidized bed. The CTD and RTD are obtained using DEM-CFD simulations with different drag models and are compared to recent PEPT experimental data (Li L. et al., 2015a). The performance of the selected drag models for fluidized bed coating applications is discussed.

Equations of particle motion
For every particle, the particle motion is described using Newton's second law where m i is the particle mass, i v  is the particle velocity, g u  is the gas velocity, g  is the gravitational acceleration, V i is the particle volume, ∆ P is the gas pressure gradient, ε s is the particle volume fraction, c,i F  is the contact force during particle-particle and particle-wall collisions and β is the interphase momentum transfer coefficient. It is noted that other forces such as a cohesion force (Thornton C. and Yin K.K., 1991) can be incorporated into Eqn. (1) in a straightforward manner.
The angular momentum of each particle is described using where I i is the moment of inertia of the particle, i ω  is the rotational velocity and i T  is the total torque acting on the particle.

Equations of gas flow
As discussed by Li L. et al. (2015b), turbulence can be neglected in this study. Hence, the continuity and momentum equations are written as follows where ε g = 1 -ε s is the gas volume fraction, ρ g is the density of air and μ g is the dynamic viscosity of air, respectively.

Soft sphere model
In accounting for interparticle forces, i.e. c,i F  (see Eqn. (1)), there are two common models available (Crowe C.T. et al., 2012), the hard sphere model and the soft sphere model. In the current study, the latter is selected since it is more appropriate for describing the behavior in a dense fluidized bed where particles can remain in contact for a long time. In this model, particle-particle collisions are described via a spring-dashpot system (Tsuji Y., 1993), as sketched in Fig. 2. In the present work, the contact forces in the normal and tangential directions are calculated using the Hertz-Mindlin theory (Hertz H., 1882).

Drag models
In Eqn. (1), the interphase momentum transfer coefficient can be specified using different drag models. In this study, five drag models are examined: the Gidaspow, Min, HKL, Beetstra and Tang drag models. Their respective expressions are presented in Table 1 and a brief summary is given below.
As mentioned in the introduction, Ergun S. (1952) obtained an expression for the average fluid particle force based on a classic study of pressure drop through packed beds. Wen C. and Yu Y. (1966) also conducted a series of fluidization experiments and presented their correlation. Later, Gidaspow D. (1994) suggested a model based on the Ergun equation for the dense regime and the Wen-Yu correlation for the dilute regime. The modification made by Gidaspow D. was found to be the best in modeling of spouted beds (Du W. et al., 2006) and has been widely used in engineering practice (Deen N.G. et al., 2007).
In other cases, the minimum value given by the Ergun equation and the Wen-Yu correlation is adopted (the Min drag model), for example, in the work of Link J.M. et al. (2005).
Using accurate numerical data from LBM simulations, Hill R.J. et al. (2001b) derived a new drag model for moderate Reynolds-number flows in ordered and random arrays of spheres. In this model, a dimensionless drag force, F(ε s ,Re p ), is defined and solved via a force balance. That is, the total force exerted on a particle by the gas flow, g s F   , is the sum of the drag force, d F  , and the buoyancy-type force due to the gas pressure gradient, b F  . Van der Hoef M. A. et al. (2005) pointed out that the drag force and the total force differ by a factor of ε g , i.e. d g g s ε F F     , and that the dimensionless drag force can be converted to the interphase momentum transfer coefficient in Eqn. (1) via the following expression: Later on, this drag model was extended by Benyahia S. et al. (2006) to cover a wider range of Reynolds numbers and particle volume fraction. The model proposed is complicated and in Table 1, a simplification suggested by Benyahia S. et al. (2006) is provided 1 . Beetstra R. et al. (2007) used the same method to derive a drag model where a more complex functional form was proposed instead of a linear scaling with Re to describe the inertial effects. Moreover, their extensive LBM simulations covered a wider range of parameters and were performed for Reynolds numbers up to 1000.
Quite recently, Tang Y. et al. (2015) proposed a new drag model based on simulation results obtained using fully resolved IBM simulations. In order to obtain highly accurate (essentially grid-independent) results, they resolved the flow past fixed assemblies of monodisperse spheres. It was stated that the grid resolution effects are consistently better resolved compared to the previous LBM simulations (Hill R.J. et al., 2001b;Beetstra R. et al., 2007).
For these different drag models, the interphase momentum transfer coefficient, β, is plotted in Fig. 3 as a function of the particle volume fraction, ε s . Since the characteristic slip velocity and the particle volume fraction are to a large extent different in different regions of the present fluidized bed, β is calculated for a slip velocity of 5 m/s and a particle volume fraction up to 0.2 for the dilute region, e.g. the Wurster tube. For the dense region, e.g. the horizontal transport region, β is calculated for a slip velocity of 1 m/s and a particle volume fraction ranging from 0.2 to 0.6. It can be seen that in the dilute region, β decreases in the following order: the HKL, Beetstra, Tang, Gidaspow and Min drag models. More specifically, β for ε s < 0.2 due to the HKL and Beetstra drag models is nearly twice that for the other three. For the moderately dense region with 0.2 ≤ ε s < 0.4, there is only a slight difference in the interphase momentum transfer coefficient between the different drag models. For the very dense region with ε s ≥ 0.4, the Beetstra, Tang and HKL drag models give a significantly larger interphase momentum transfer coefficient than the Gidaspow and Min drag models.

Experimental
Recent PEPT experimental data presented by Li L. et al. (2015a) are used in this study to evaluate different drag models by comparing the calculated particle velocity, cycle and residence time distributions with the measured values.
In the PEPT experiments, a tracer particle labeled with Fluorine-18 is incorporated into the fluidized bed. As the tracer particle moves around in the vessel, the radioisotope decays and releases back-to-back γ-rays. These γ-rays are detected via two large position-sensitive detectors located on both sides of the fluidized bed. Each pair of γ-rays can be used to calculate the location of the tracer particle by means of three-dimensional triangulation. More details about the technique and the algorithms for data processing can be found in the literature (Parker D.J. et al., 1993;Forster R. et al., 2000).
The fluidized bed used in the PEPT experiments was based on the STREA-1 TM laboratory fluidized bed from Aeromatic-Fielder. The fluidized bed is 380 mm high, with top and bottom diameters of 250 mm and 114 mm, respectively; detailed dimensions are provided elsewhere (Li L. et al., 2015a). A nozzle with a diameter of 5 mm was employed in order to avoid the numerical difficulties in simulating sonic or even supersonic flow while retaining the feature of high jet velocity due to the atomization air flow. The nozzle was only a big circular orifice and no liquid solution was introduced. The fluidization air flow was supplied through a bowl-shaped distributor, consist- 1 It must be noted that this extension was developed following the original derivation by Hill R. J. et al, where the total force exerted on the particle was referred to as the drag force. Thus another factor of the gas volume fraction is required when converting the dimensionless drag force to the interphase momentum transfer coefficient using Eqn. (5).
ing of a central base and an outer annular region. The central base of the distributor is fully open, while the outer annular region of the distributor has a number of orifices that correspond to approximately 4 % of the area of the outer annulus region. A wire mesh screen was put over the distributor to prevent particles from falling down below the distributor.
Microcrystalline cellulose (MCC) pellets, a common material in the pharmaceutical industry, were employed in the PEPT experiments. The particle size was measured using QICPIC Particle Size Analysis (Sympatec GmbH), and the particle density was determined using mercury intrusion porosimetry (Micromeritics AutoPore III 9410). The material parameters are given in Table 2; they were selected to match those of MCC. A photograph of a sample of the MCC pellets is provided in Fig. 4.
The shape of particles may affect the drag force acting on the particles. In this case, the so-called sphericity factor can be used to account for the effect of the shape of particles (see, e.g. (Zastawny M. et al., 2012)). In order to evaluate the drag force more accurately, correlations are required according to, e.g. Crowe C.T. (2006), using the sphericity of particles provided in Table 2. In addition, the shape of particles may affect the porosity in the dense downbed region and therefore the air flow passing through this region. However, this is beyond the scope of the current study.

General
A fully coupled multi-phase flow solver, MultiFlow (www.multiflow.org), was employed in this study to solve the equations of particle motion and gas flow. The fluid mesh, which is refined near the nozzle and in the Wurster tube, contains approximately 65000 hexahedral cells. Simulations of single-phase gas flow were performed to ensure grid independence.
In order to specify the boundary conditions in the DEM-CFD simulations, the flow rates of the atomization and fluidization air flow, as well as the flow distribution at the distributor are required. The global flow rates were measured using a rotameter in the PEPT experiments (Li L. et al., 2015a). However, it was not possible to measure the flow distribution in the experiments. Instead, a singlephase CFD model including the air supplying chamber, the distributor, the wire mesh screen, and the fluidized bed was developed in Ansys Fluent to determine the flow distribution at the distributor. The results (not shown) obtained using this single-phase CFD model show that the air flow passing through the central base of the distributor is between 53 % and 64 % of the total fluidization air flow. Since the particle velocity profile is an indication of the flow distribution, DEM-CFD simulations with different flow distributions at the distributor were performed to compare the particle velocity at different heights above the distributor with the corresponding experimental data. These additional simulations were used to verify that a flow distribution for which approximately 55 % of the fluidization flow enters normal to the boundary in the central base of the distributor and 45 % enters normal to the boundary in the outer annular region (Li L. et al., 2015b) gave adequate agreement with the experimental data.
The atomization nozzle is modeled as a flow normal to the inlet boundary, while the outlet of the domain is considered to be a pressure outlet. The walls of the fluidized bed and the Wurster tube are set to be no-slip for the gas flow.
In order to obtain a reasonable compromise between Table 2 The material parameters used in the numerical simulations (Roberts R. et al., 1994;Bharadwaj R. et al., 2010;Darelius A. et al., 2007;Li L. et al., 2015b) Variable Unit Value Air density kg/m 3 1.204 Air viscosity Pa·s 1.837 × 10 -5

Particle diameter μm 1749
Particle density kg/m 3 1420 data quality and computational cost, the CTD and RTD in different regions were evaluated for different simulation times. It was found that the CTD becomes fairly steady after a simulation time of 25 s. It was also found that a simulation time of 10 s was appropriate to determine the RTD in different regions, except in the horizontal transport region.

Simulation parameters
The numerical parameters and the operating conditions are presented in Tables 2 and 3, respectively. The partition gap is set to 15 mm and the Wurster tube length is 150 mm. The fluidization and atomization air flow rates are 73.3 m 3 /h and 3.50 m 3 /h, respectively. A base case with a batch size of 200 g was selected in this study to investigate the drag models. In addition to the base case, which uses the Gidaspow drag model (Run #1), DEM-CFD simulations with the same operating conditions are performed for the other drag models in Runs #2-5 and for larger batch sizes in Runs #6-9. A typical simulation was run in parallel using 32-40 CPU cores and took about 12 hours wall clock time per simulated second.

Post-processing
From the point of view of a typical coating cycle, we define the start of a particle coating cycle as the point when a particle enters the spray zone from the horizontal transport region or the lower region in the Wurster tube. The particle coating cycle then ends when the particle reappears in the spray zone after having subsequently traveled through the Wurster tube, the fountain, downbed, and horizontal transport regions. When one particle coating cycle ends, the next one begins. The particle cycle time is thus defined to be the time it takes for a particle to complete one cycle.
In the PEPT experiments, it was difficult to detect the tracer particle in the lower region of the fluidized bed. In the experiments, it was therefore not possible to identify a spray region and the particle cycle was defined so that it begins and ends in the Wurster tube. A spray zone is therefore not defined in the simulations either, and the particle cycle is defined in the same way as in the experiments.
It is noted that in the PEPT experiments, only one tracer particle was incorporated into each run and continuously followed for 1.5 hours. In the simulations, on the other hand, different particle trajectories are followed for 25 s. As a result, the cycle and residence time distributions were calculated from different cycles for the same tracer particle in the experiments, while they are obtained from many different particles in the DEM-CFD simulations. For the pseudo-steady particle dynamics in this fluidized bed, it is thus assumed that the information obtained (such as the particle velocity or the cycle time distribution) by following one particle for a long time is equivalent to the information obtained by following all particles for a short period of time (Li L. et al., 2015b).

Results and discussion
In Table 4, the mean cycle time and mean residence times of particles in different regions are summarized for the batch size of 200 g; the corresponding experimental data is also provided. Figs. 5 and 6 show the cycle time distributions and mean residence times of particles in different regions, respectively, including values both measured in the PEPT experiment and calculated from the DEM-CFD simulations for different drag models. For all selected drag models, the general feature of the particle residence times in different regions is captured. As can be seen in Table 4, the fraction of the cycle time that particles spend in different regions corresponds closely to PEPT measured values for the different drag models.
However, the cycle time and the residence time in the different regions, particularly in the Wurster tube and in the horizontal transport region, exhibit distinct differences for the different drag models. It can be seen that both the mean cycle time in Table 4 and the CTD in Fig. 5 calculated with the Gidaspow and Tang drag models agree quantitatively with the experimental data. For the Min drag model, there is a greater fraction of particles with a longer cycle time, which can be attributed to a somewhat slower acceleration in the Wurster tube and a longer residence time in the Wurster tube (see in Table 4). Interestingly, the CTDs predicted using the HKL and Beetstra drag models peak at a shorter cycle time and are narrower than the others. These two effects result in a much shorter mean cycle time for these two drag models.
In order to understand the cause of these differences, the vertical particle velocity at selected heights above the bottom of the fluidized bed, i.e. 45 mm, 75 mm and 105 mm, is plotted in Fig. 7. At 45 mm, which is the lowest location where reliable PEPT data could be obtained, the particle velocity is close to the PEPT measured velocity in the Wurster tube for all drag models studied, while the HKL drag model shows better agreement with the PEPT data in the downbed region. At 75 mm, all the drag models give very similar results. At 105 mm, only the particle velocity in the Wurster tube is shown, and it can be seen that all the drag models overpredict the particle velocity to a great extent in a central region corresponding to a radial location smaller than approximately 5 mm. This latter difference is probably due to the fact that in the experiments, there is turbulence in this central region induced by the atomization air flow, while the simulations employ a laminar model to predict the air flow. Nevertheless, since this central region contains only a few particles, it is not expected to affect the general particle motion in the fluidized bed and the predicted CTD and RTD. It can also be seen that the HKL and Beetstra drag models give a distinctly higher velocity than the other drag models in this central region. Outside of this central region, the particle velocity for all of these drag models shows better agreement with the PEPT measured velocity. But a larger overprediction of the particle velocity compared to the PEPT measurement-though smaller than that in the central region-can be found for the HKL drag model than for the other drag models. This is similar to the results presented by Link J.M. et al. (2005) for a pseudo-two-dimensional spout-fluid bed, where it was concluded that the HKL drag model has the best predictive capabilities. Fig. 8 shows the particle volume fraction examined at two heights above the bottom of the fluidized bed, corresponding to the inlet (15 mm) and outlet (165 mm) of the Wurster tube, respectively. It can be seen that the particle volume fraction is lower for the HKL and Beetstra drag models at the inlet of the tube and vice versa at the outlet of the tube. Since the Beetstra, Tang and HKL drag models give a larger interphase momentum exchange coefficient in the dense downbed region (cf. Fig. 3), there is a larger pressure drop in the downbed region for these models. As a result of the larger pressure drop in the Fig. 6 The residence times of particles in different regions for the base case (the secondary y-axis is used for the residence time in the horizontal transport region due to the higher values). downbed region, a certain amount of fluidization air flow is diverted from this region into the Wurster tube. This is supported by a lower air velocity at 15 mm in the downbed region for these models, as shown in Fig. 9.
Since the particles are accelerated by the air flow through the Wurster tube, the particle volume fraction at the outlet of the tube decreases compared to that at the inlet for each drag model. For the HKL and Beetstra drag models, the larger interphase momentum exchange coefficient in the dilute Wurster tube (cf. Fig. 3) results in a higher particle velocity in the tube compared to other drag models. Correspondingly, a larger mass flow rate through the Wurster tube is expected and a higher particle volume fraction at the outlet of the tube is obtained for the HKL and Beetstra models than for the other drag models.
The higher particle volume fraction at the inlet of the Wurster tube for the Tang, Gidaspow and Min drag models can be attributed to a lower air flow rate in the Wurster tube and a lower interphase momentum transfer coefficient. The resulting slower acceleration of particles in the Wurster tube leads to particles having a greater risk of recirculating in the Wurster tube and moving back towards the horizontal transport region from the lower edge of the Wurster tube. This type of behavior was also observed in the PEPT experiments (Li L. et al., 2015a) and is sketched in Fig. 10. Thus, it is of interest to examine the number of particle recirculations in the Wurster tube, as shown in Fig. 11. It is significant in Fig. 11 that for the HKL and Beetstra drag models, the fraction of ideal cycles, i.e. the fraction of particle cycles without any recirculation, is approximately 63 % and 43 %, which is much higher than   the PEPT measured value of 12 %. The reduced number of particle recirculations for the HKL and Beetstra drag models can thus decrease the residence time in the Wurster tube and prevent a very long residence time in the horizontal transport region. This again emphasizes how important the detailed particle motion in the Wurster tube is to the operation of this type of fluidized bed. In order to explore a scenario where the dense region is likely to have a greater effect, simulations of larger batch sizes of 400 g and 600 g were performed for the HKL and Tang drag models. The results obtained for the mean particle cycle times are presented in Fig. 12 for different batch sizes. In agreement with the recent PEPT experiments, all selected drag models predict a shorter cycle time for a larger batch size. Notably, the Tang drag model shows the closest agreement for the increased batch size of 400 g. The results also indicate that, as the batch size increases, the deviation from the experiments is reduced for the HKL drag model, since the effect due to particle recirculations in the Wurster tube is reduced.

Conclusions
In this study, DEM-CFD simulations of particle motion in a Wurster fluidized bed have been performed in order to study the effect of drag models on particle cycle and residence times. A comparison between the calculated particle cycle and residence time distributions from the simulations and recent PEPT experimental data was made. It was found that the Gidaspow, Min and Tang drag models predict the particle cycle time, the CTD and the RTDs qualitatively and quantitatively. The Tang drag model also demonstrates better agreement with experimental data for a larger batch size. The HKL and Beetstra drag models slightly overestimate the particle velocity in the Wurster tube. This overestimation of the particle velocity in the Wurster tube results in a much smaller number of recirculations, and therefore a much shorter residence time in the horizontal transport region and a much shorter cycle time. For the larger batch sizes, where the effect of particle recirculations is less pronounced and the dense region becomes more dominant, the agreement between the experimental data and the predictions of the HKL drag model improves. The present results thus suggest that the drag model must be selected carefully for the purpose of studying residence time distributions and, hence, coating processes in fluidized beds.    i   rotational velocity of the particle (rad/s) ε g gas volume fraction (-) ε s particle volume fraction (-) β interphase momentum transfer coefficient (kg/ (m 3 ·s)) μ g dynamic viscosity (Pa·s) ρ g air density (kg/m 3 )