Abstract

A coupled numerical code of the Euler-Euler model and the population balance model (PBM) of the liquid-liquid dispersions in a spray fluidized bed extractor (SFBE) has been performed to investigate the hydrodynamic behavior. A classes method (CM) and two representatively numerical moment-based methods, namely, a quadrature method of moments (QMOM) and a direct quadrature method of moments (DQMOM), are used to solve the PBE for evaluating the effect of the numerical method. The purpose of this article is to compare the results achieved by three methods for solving population balance during liquid-liquid two-phase mixing in a SFBE. The predicted results reveal that the CM has the advantage of computing the droplet size distribution (DSD) directly, but it is computationally expensive if a large number of intervals are needed. The MOMs (QMOM and DQMOM) are preferable to coupling the PBE solution with CFD codes for liquid-liquid dispersions simulations due to their easy application, reasonable accuracy, and high reliability. Comparative results demonstrated the suitability of the DQMOM for modeling the spray fluidized bed extractor with simultaneous droplet breakage and aggregation. This work increases the understanding of the chemical engineering characteristics of multiphase systems and provides a theoretical basis for the quantitative design, scale-up, and optimization of multiphase devices.

1. Introduction

Liquid-liquid dispersions are widely used in the chemical industry for conducting several operations with significant implications in the agrochemical, clinical medicine, and biopharmaceutical fields [1]. One of such processes is spray fluidized bed extractor (SFBE) which is widespread in chemical industry because the disperse phase in a SFBE has several special characteristics such as stable flow behavior and excellent capability of exchanging mass (energy or momentum) with the continuous phase. In SFBE, liquid (as continuous phase) is mixed by pumping oil (as disperse phase) from the bottom of the extractor. Droplets of two phases move in an apparently random manner in the bed and encounter collisions with each other [2]. These characteristics influence the mass (energy or momentum) transfer between two phases in the fluidized bed extractor. The distribution of disperse phase volume fraction controls the flow pattern of two phases and affects the intermixing degree. However, some knowledge regarding disperse phase distribution and continuous phase flow patterns is very difficult to achieve with experiments. It can be obtained by numerical simulations.

This information can be obtained with computational fluid dynamics (CFD) coupled with a population balance model (PBM) [3]. Nowadays, PBM is widely used to describe some procedures such as interphase mass (energy or momentum) transfer [4]. In multiphase flows involving a droplet size distribution, the population balance equation (PBE) is required to describe the variations in the droplet population, in addition to mass, energy, and momentum balances [5]. The numerical methods for solving PBE can be classified into the method of classes (CM) [6], sectional method (SM) [7], the Monte Carlo method (MCM) [8], method of moments (MOMs) [9], the method of characteristics [10], the method of weighted residual or orthogonal collocation [11], and the high-resolution finite volume methods [12].

In MOMs, the PBE is transformed into a series of transport equations for distribution moments. It is generally sufficient to solve a few moment equations. However, the MOMs encounter a technical barrier (closure problem), which limits its application in the academic field [13]. In order to solve this problem, a set of new improved MOMs have appeared such as the standard method of moments (SMOM) [14], the direct quadrature method of moments (DQMOM) [15], the fixed pivot quadrature method of moments (FPQMOM) [16], and the extended quadrature method of moment (EQMOM) [14]. The conservation equations about CM and MOMs are presented in Table 1.

The extensive research has been displayed in recent years to the application of these numerical methods to solve PBE in fluidized bed. Luo et al. simulated the gas-solid flows in polymerization FBR using a CFD-PBM and later extended the CFD-PBM to simulate the gas-solid flows in the multizone circulating polymerization reactor [17]. Professor Bart of Kaiserslautern in Germany used CM and QMOM to simulate two-phase flow behavior in a rotating disc extraction column which proved the population balance model is feasible in simulation of droplet size distribution in extraction columns [18]. Marchisio demonstrated that the QMOM and DQMOM are identical in simulated results for monovariate cases, whereas the DQMOM has an obvious advantage for multiphase problems [19]. Si–Si Liu used QMOM to study the effects of operating conditions on the growth rate in detail with the observation of fluidization phenomenon during the evolution of particle growth [20]. Drumm et al. used both the CM and the QMOM to solve the coupled CFD-PBM for an RDC. The results showed that the computational time of the QMOM was much lower [21]. Attarakih solved the breakage and coalescence of liquid droplets is proposed, which is based on SQMOM [22].

As can be seen from the above literatures, the QMOM and the DQMOM are the most promising MOMs to predict PSD or DSD in the fluidization system. Meanwhile, the CM has been one of the most popular numerical methods. However, there is so far no open article on assessment of two kinds of numerical methods in simulating liquid-liquid dispersions in SFBE. Moreover, the understanding of the effect of MOMs on the liquid-liquid flow behavior remains unclear, which needs further investigation. In this research, the method of classes (CM) and two representative MOMs (QMOM and DQMOM) are used to solve the PBE in the coupled mode in the ANSYS Fluent solution platform. The purpose of this study is to investigate on the liquid-liquid dispersion characteristics in a spray fluidized bed extractor.

2. Computational Model

2.1. Two-Fluid Model

The Eulerian-Eulerian approach is applied to calculate the flow field distribution of liquid-liquid dispersions in the SFBE [23]. In addition, it is required to make the following assumptions: (1) tow phases are treated as interpenetrating phases and (2) all fluids are assumed to be in an equilibrium state.

The continuity equation of the disperse phase is formulated as [24]

The average velocities of the disperse phases are calculated by solving the corresponding momentum balance equations [25]:where is the disperse phase velocity; is the disperse phase density; represents the interphase interaction; and represent the directions. Equations (3) and (4) give the transport equation for the turbulent kinetic energy term () and the transport equation for the turbulent dissipation term (), respectively [2]:

According to Launder and Spalding [26] recommended values, the constants in above models are equal to the numerical values below [27]:

The energy conservation equation for the disperse phase [17] can be expressed aswhere is the granular temperature; is the disperse phase volume fraction; is the disperse phase pressure; is the collisional rate of energy dissipation; is the droplet fluctuating energy; denotes a dissipation term due to the interphase interaction; and denotes a production term due to interactions with fluid turbulence. The first two terms on the right-hand side of Equation (6) represent the generation of granular energy due to production and diffusion. These disperse phase properties can be determined as a function of granular temperature according to the following expressions:

The granular temperature measures the random oscillations of droplets. Therefore, it can represent the fluctuation kinetic energy of droplets. The equations of granular temperature are given as follows [28]:

The interaction between two phases is considered with the interphase interaction (), which includes lift force, turbulent diffusion force, drag force, and virtual mass force (). The lift force () is so difficult to investigate. [29] It currently does not consider the lift force effect in order to facilitate the simulation. The following is an introduction for the interaction force:(1)When a droplet inside a fluid domain accelerates, part of the surrounding fluid volume is likely to be displaced during its motion. The virtual mass force represents the inertia force that is added to the system due to this situation. The virtual mass force for droplets can be expressed as follows [30] ( represent the accelerations of two phases):(2)The drag force coefficient has a variety of models. The following Schiller-Naumann-Pb equation is employed for the drag force estimation [31]:(3)The turbulent dispersion force results from the combined actions of the drag force and the surrounding turbulent eddies of the continuous phase. This effect of the turbulent dispersion force results in moving the droplets from high to low concentration areas. Lahey et al. launch the following equation to express turbulent dispersion force [31]:

2.2. Population Balance Model
2.2.1. Balance Model

The PBM is based on the conservation equation of number density function [32]. A general form of PBM can be expressed aswhere

The moments of method are defined as follows [33]:

Substituting Equations (14) into (12)where is the number density function with droplet diameter ; is the droplet flux due to growth rate; is the droplet velocity vector; and are the birth and death rates of the droplets for aggregation; and and are the birth and death rates of droplets for breakage, respectively.

2.3. Breakage and Aggregation Kernels

The Laakkonen breakage kernel is expressed as the product of the breakage frequency and the daughter PDF , where

The constants C2 = 2.52, C3 = 0.04, and C4 = 0.01. and are the daughter and parent droplets volumes, respectively. The majority of aggregation kernels are semiempirical formulas compared to the breakage kernel. One of the most popular is the Luo [34] model, the general aggregation kernel is defined as the rate of droplets volume formation as a result of binary collisions of droplets with volumes and :

2.4. CFD-PBM Coupled Model

In this study, CFD-PBM is used to describe the liquid-liquid dispersion characteristics between the two-phase flow fields and combine the advantages of two models. Figure 1 shows the algorithm for simulating the liquid-liquid mass (momentum or energy) transfer in the framework of the CFD-PBM-coupled model. The phase volume fraction, droplet velocity, and turbulent energy dissipation rate achieved by solving mass conservation, momentum conservation, and energy conservation equations are used to calculate the droplet breakage and aggregation kernels in PBE [35].

2.5. Numerical Methods
2.5.1. CM

The CM divides the complete region of droplet size into several classes, where each class corresponds to the size interval defined as [] [6]. The density function within class and numbers of droplets can be represented as

The Sauter mean diameter () is usually used as the mean droplet size. Sauter average size in extractor can be expressed as follows:

2.5.2. QMOM

QMOM is proposed by McGraw [36]. It replaces the exact closure needed by MOMs with an approximate closure. The QMOM are widely attractive for its special property and extensive application. This method is based on solving transport equations for the moments of the DSD:

By applying QMOM, the transport equation for the moment of is obtained as follows:

2.5.3. DQMOM

The DQMOM is the proper orthogonal decomposition of a number density distribution function by inserting the control equation, deducing the PBE by tracking orthogonal basis [37]. The moments of the population density are given by the quadrature approximation incorporating the weights () and abscissas () with respect to the two properties:

The above linear system can be written in the matrix form aswhere the coefficient matrix can be expressed as follows:

The vector of can be represented as

The right side of Equation (23) is the source term involving droplet aggregation and breakage. The growth term is accounted directly in Equation (23):

The source term for the moment (k = 0, ..., 2N − 1) can be expressed as

2.6. Reactor and Simulation Description

The geometry and computational meshes of the SFBE are shown in Figure 2. The simulation is based on the axisymmetric model, with the mesh size of 1 mm. The following Tables 24 are the required boundary conditions, the model parameters, and the numerical schemes of simulation.

3. Results and Discussion

3.1. Droplet Size Distribution

The mean diameter of the droplets is one of the most important hydrodynamic parameters for SFBE, and usually characterized by the Sauter mean diameter in the literature [38]. Figure 3(a) shows the droplet Sauter diameter distribution in the radial direction at the heights of 0.2 m using the three methods. As expected, the highest values of droplets mean diameter are deduced from DQMOM, the second largest values of droplets mean diameter are deduced from CM, and the lowest values of droplets mean diameter are found in QMOM. Figure 3(b) displays the predicted droplet size distribution (DSD) along the axial direction. As shown in Figure 3(b), using the QMOM and the DQMOM, a similar curve can be obtained. However, using the CM, at the bottom of the SFBE (0–0.05 m), the droplet Sauter diameter increases with the bed height. It indicates that the collision of droplets occurred mainly in the nozzle. Furthermore, in this region, the dynamic pressure, speed, granular temperature, and turbulent kinetic energy of droplets achieve maximum, and the collision probability of droplets is greatly increasing, which results in the coalescence frequency increases too. So the droplets easily break into many small pieces. However, the weakly influence of interdroplets makes the droplets generation rate equal to the death rate. Consequently, the droplet size does no longer change when droplets enter into the extractor. Compared with the axial direction of the nozzle, the droplet size of CM has little change in the radial direction, due to the smaller force between droplets in this direction probably. An appropriate method should be chosen in the study. When precise droplet size distribution is required, CM is applied.

3.2. Spray Process

Figure 4 shows the spray processes from 0.1s to 1s using all three numerical methods. A comparison between the simulation results and previous research performed by Bright A [39] is also shown. It displays the complete process of formation, development, fragmentation, and the jet flow behavior of oil bubble (disperse phase) and shows the conformation of fountain visually. It can be observed from Figure 4, see that the distribution of disperse phase has somewhat difference which is caused by different droplet sizes. Firstly, a thread of oil comes out of the nozzle to be injected into the fluidized bed. Then, the oil bubble appears and moves along the flow direction and develops under the pressure of continuous phase to become thinner at the downstream. Subsequently, the oil bubble breaks under the increasing of shear force to give birth of first droplet. When the disperse phase is entering into the bed completely, oil bubble bursts occur, and oil fully mixes with the liquid. The oil bubble formation is drastically affected by the absence and presence of droplets. Inertial force that impedes the bubble detachment acting on forming bubbles increases because of the presence of droplets. As can be seen from Figure 4(b), the movement represented by DQMOM is more slowly due to the larger droplet size, and the volume fraction distribution is more extensive when compared with the case of QMOM. The difference of droplet size overcomes the interfacial tension easily. In this condition, the disperse phase thread flows slowly. The simulated results are in good agreement with Bright A’s experiment in Figure 4(d). Meanwhile, the process is different in the methods of classes. The droplets acquire a massive amount of kinetic energy and inertia force in the nozzle. The reason is that force caused by the collisions between equal sizes of the droplets is no more than it from different sizes of the droplets, leading to its jet distance that sprays into the continuous phase is more far away compared to MOMs. The oil bubble breaks up to give birth of the first droplet, and the latter is followed by the droplet with small size; the effect of continuous phase viscosity could be observed on the droplet size seem to be a little bigger. The decrease in the viscous stress of continuous phase allows us to approximately reduce the shear forces and results in slower rupture of the disperse phase; consequently some droplets are generated after the mother droplet with bigger size breaks. Based on these results, it can be concluded that the droplet size influences on the complete spray process. In MOMs, it has a certain resistance force caused by the wall effect, which leads to reducing the diffusion rate of the disperse phase and inhibits the increase of oil bubble.

3.3. Evolution of Mixture Density

The time evolutions for the mixture density predicted with three methods are compared in Figure 5. It is obvious to figure out that the density of the mixture arrived at an equilibrium state after fluidizing for 6s. As it is seen for both MOMs, the density of mixture is quite similar. Besides this, the disperse phase is quite homogeneously distributed. This implies that mixing in the SFBE is quite fast and capable of distributing the droplets around.

3.4. Velocity Distribution

Figure 6 illustrates the comparison of the disperse phase velocity profile in the radial direction at various heights. The existence of high velocity at the center and low velocity near the walls is clearly predicted. This means droplets are carried up by the liquid in the center region and droplets fall down along the wall in the annular region. In comparison, Figure 6 reveals that the velocity profile predicted by the QMOM is flatter throughout the radial position than other numerical methods. A high inlet liquid velocity will also give a high droplet velocity due to more energy input into the bed. The oil droplets induced turbulence is strong and the flow structure close to isotropic.

3.5. Numerical Simulation of the Nozzle Area
3.5.1. Plunging Pressure Distribution

The plunging pressure is of high importance in the operation of an extractor. The simulation results of the plunging pressure are shown in Figure 7. The pressure radial distribution of QMOM and DQMOM is simulated, which indicates that MOMs have the high degree of flattening and the preponderance of plunging pressure from jet center is no longer existing under the flow rate about 1.8 m/s obviously. However, the pressure values of CM are heterogeneous, that is to say the CM based on different droplet sizes is not suitable for the investigation, and the simulation results of the MOMs are more accurate.

3.5.2. Granular Temperature Distribution

Figure 8 compares the numerical granular temperature of droplets and the experimental data of He et al [40] at various nozzle heights. The granular temperature reveals the intensity of droplets fluctuation in a fluidized bed. The droplets fluctuation is depending on the interdroplet collision, convection, diffusion, shear, and interaction with the continuous phase. Figure 8 indicates the granular temperature is decreased with an increase of nozzle height. It shows that the granular temperature is high in the bed center and decreases toward the walls up to 0 values, which means the velocity fluctuation in this region is much stronger than in other regions. At impact, the granular temperature increases rapidly as the droplet suddenly applies forces. The granular temperature decreases as the disperse phase penetrates in to the fluidized bed and finally comes to rest at a finite penetration depth. At the height of 0.045 m, the granular temperature is close to zero because droplets move hardly near the bottom of the nozzle. Generally, all the numerical methods can qualitative predict the droplets granular temperature distribution. Due to the difference of the droplet size calculated by three methods, the granular temperature predicted is also distinctive. As shown in Figure 8, the droplet granular temperature simulated by CM is overall lower than that simulated by MOMs. In addition, the droplet granular temperature simulated by DQMOM is slightly higher than that simulated by QMOM in the bed center. The reason is that the strong fluctuation of the small droplets arises from oil bubble breakage near the bed surface. The figure also shows that the droplet granular temperatures obtained using DQMOM agree well with the measured results from He et al.

3.6. Numerical Simulation of the Bed
3.6.1. Droplet Dynamic Pressure Distribution

Comparisons of droplet dynamic pressure profiles between the height direction and the radical distance are presented in Figure 9. An increase in extractor axial height induces a decrease in the droplet dynamic pressure. The relationships between droplet dynamic pressure and axial height can be expressed as , where is the bed voidage. Since the disperse phase holdup by the three methods has only a little difference in axial height in Figure 4, the distribution of the void fraction is also similar. When using the nozzle with 0.01 m inside diameter, increasing violent motions of droplets which generates an increase in dynamic forces of droplets that causes the fluctuations of simulated pressure values. The DMOM can predict more severe fluctuation than the other two methods.

3.6.2. Mixture Solution Turbulent Properties

In Figure 10, the mixture phase turbulent performance along the radical distance is compared to the experiment data of Chen et al [41]. The turbulent kinetic energy changing trend is roughly the same as the turbulent dissipation rate in the radial direction. In general, the turbulent kinetic energy is proportional to the droplet velocity. From the above analysis, the radical velocity of droplets is firstly increasing and then decreasing. The change of turbulent kinetic energy is approximately the same. The droplets fall down rapidly in the near wall region due to the wall effect. However, the turbulent dissipation rate is increasing in the area close to the wall since the larger values of the dissipation rate are predicted by the CM near the wall. In general, the dramatic increasing of the turbulent dissipation rate in MOMs has two main reasons: (1) the speed fluctuation caused by the droplets wake flow and vortex is violent and (2) the instability of buoyancy caused by selective aggregation phenomenon of the density gradient is also increasing. The turbulent dissipation rate has direct effects on the calculations of aggregation and breakage rates. The larger the turbulent dissipation rate is, the higher the droplet breakage rate becomes. Hence, the calculation results estimated by QMOM are smaller than the results estimated by DQMOM in Figure 10(b).

3.6.3. Disperse Phase Volume Fraction

Figure 11 shows the simulated disperse phase volume fraction profile in the radial direction at different heights of the extractor using the three numerical methods. It concluded that the radial distribution of disperse phase volume is heterogeneous, lower at the central region, and higher near the wall region due to the wall attachment effect. Meanwhile, as there is an increase in the axial height, the disperse phase volume decreased slightly and the peak radial disperse phase volume disappeared, so better homogeneous radial distribution of the disperse phase volume is observed. Furthermore, the phenomenon that MOMs can predict the bigger disperse phase volume fraction than the CM when the extractor height is relatively low while obtaining the smaller one when the height is high can be seen in Figure 11.

4. Conclusions

In this work, the numerical simulations are conducted to investigate the effect of hydrodynamic behavior of the liquid-liquid dispersions in the SFBE. This work aims at choosing the most appropriate method in CFD-PBM model. Three representative methods, CM, QMOM, and DQMOM, are implemented in our coupled model for comparison. The simulated results of the disperse phase spray process and some relevant parameters are discussed to illustrate the differences. The study acquires the fundamental procedure and advantages and disadvantages of the three methods. Our simulation results show that all the three methods can obtain the reasonable mixture solution turbulent properties, time-averaged flow field, and granular temperature distribution at steady state. Without considering the precise DSD, both the QMOM and DQMOM are more suitable and comparing with CM. DQMOM is more dominant due to the small amount of calculation. Generally, the application of DQMOM for solving FBE in CFD-PBM is reasonable. The more reliable results can be obtained relative to CM and QMOM. Further researches on the improvement of the CFD-PBM in SFBE are in progress in our group.

Nomenclature

:Source term of the weight control equation of the direct quadrature method of moments
:Source term of the source orthogonal basis weight control equation of the direct quadrature method of moments
:The fluctuating velocity of droplets
:Coefficients in the turbulence model
:Drag coefficient
:Coefficients in virtual mass
:Coefficients in turbulence diffusion
:Droplets diameter (m)
:Reactor diameter (m)
:Restitution coefficient
:Gravitational constant (m·s−1)
:Radial distribution function
:Breakage frequency, the fraction of particles of volume V breaking per unit time
:The generation of turbulence kinetic energy resulting as the effect of buoyancy
:The production rate of turbulent kinetic energy
:Specified number of moments
:The 0th moment of number density function
:The kth moment of the number density function
:The total number of data over a given time period
:Drag force (N⋅m−3)
:Turbulent dispersion force (N·m−3)
:Virtual mass force (N·m−3)
:The number of droplets per unit volume
:Particulate phase pressure (Pa)
:Breakup probability for a droplet of volume Vi breaking when hit by a volume Vj
:Turbulent Prandtl number
:Reynolds number
:The source term
:Flow time (s)
:The characteristic velocity of collision of two droplets with diameters di and dj (m⋅s−1)
:The instantaneous droplet velocity
:th components of the space position vector
:Diameter ratio of di and dj
:The contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate
Greek letters
:Interphase momentum transfer coefficient (kg⋅m−3⋅s−1)
:Probability density function (PDF) of particles breaking from volume V to a droplets of volume V
:Turbulence dissipation rate (m2⋅s−2)
:Viscosity of disperse phase (pa⋅s)
:Continuous phase density (kg⋅m−3)
:Disperse phase density (kg⋅m−3)
:Turbulent Prandtl numbers for k
:Turbulent Prandtl numbers for ε
:Diffusion coefficient of the number density distribution function (m2⋅s−1)
:Frequency of collision
:Discrete weight of number density distribution function
Subscripts
:Continuous phase
:Disperse phase.

Data Availability

All data generated or analyzed during this study are included in this published article. The data included in this study are available upon request from the corresponding author.

Disclosure

An earlier version of this study was presented as a poster in 2016 AIChE Annual Meeting proceedings.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the Science and Technology Program of Sichuan, China (no. 2018CZ0025), the Foundation for High-level Talents in Higher Education of Sichuan University of Science and Engineering (no. 2017RCL68), the Foundation of Sichuan Educational committee of China (no. 17ZA0280), the Sichuan Province Universities and Colleges Key Laboratory Fund of Process Equipment and Control Engineering (no. GK201709).