# CFD-PBM Coupled Simulation of Liquid-Liquid Dispersions in Spray Fluidized Bed Extractor: Comparison of Three Numerical Methods

Authors: Dan Zheng; Wei Zou; Chuanfeng Peng; Yuhang Fu; Jie Yan; Fengzhen Zhang
Journal: International Journal of Chemical Engineering
 (2019)

---

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


---

## Body

## 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.Table 1
The numerical methods for PBM equations.
Method
Basic equation
Result
Advantage
Disadvantage
CM
∂
n
i
/
∂
t
+
∂
n
i
u
α
/
∂
x
α
−
∂
/
∂
x
α
Γ
t
∂
n
i
/
∂
x
α
=
n
i
−
1
G
L
i
−
1
/
L
i
−
L
i
−
1
−
n
i
G
L
i
/
L
i
+
1
−
L
i
+
S
i
¯
Accurate size distribution
Easy application
High computational cost
MOM
∂
m
0
/
∂
t
+
∂
m
0
u
α
/
∂
x
α
−
∂
/
∂
x
α
D
∂
m
0
/
∂
x
α
=
S
0
¯
t
k
=
0
Moments and mean size
Low computational cost
Nonclosed equation set
∂
m
k
/
∂
t
+
∂
m
k
u
α
/
∂
x
α
−
∂
/
∂
x
α
D
∂
m
k
/
∂
x
α
=
k
∫
0
∞
L
k
−
1
G
L
n
L
;
x
→
,
t
d
L
+
S
k
¯
t
k
>
0
SMM
∂
ρ
m
k
/
∂
t
+
∂
ρ
m
k
u
α
/
∂
x
α
=
∫
0
∞
L
k
n
L
;
x
→
,
t
d
L
Moments and mean size
Low computational cost
Nonclosed equation set
QMOM
∂
m
k
/
∂
t
+
∂
m
k
u
α
/
∂
x
α
−
∂
/
∂
x
α
D
∂
m
k
/
∂
x
α
=
k
∑
i
=
1
N
L
i
k
−
1
G
L
i
ω
i
+
S
k
¯
t
Moments and rough size distribution
Closed equation set and low computational cost
–
DQMOM
∂
ω
i
/
∂
t
+
∂
ω
i
u
α
/
∂
x
α
−
∂
/
∂
x
α
D
∂
ω
i
/
∂
x
α
=
a
i
Moments and rough size distribution
Closed equation set and low computational cost
Complex source term
∂
ξ
i
/
∂
t
+
∂
ξ
i
u
α
/
∂
x
α
−
∂
/
∂
x
α
D
∂
ξ
i
/
∂
x
α
=
b
i
FPQMOM
∂
m
k
/
∂
t
+
∂
m
k
u
α
/
∂
x
α
=
k
∑
i
=
1
N
V
i
k
−
1
G
V
i
ω
i
+
1
/
2
∑
i
=
1
N
ω
i
∑
j
=
1
N
ω
i
V
i
+
V
j
k
β
i
,
j
Moments and rough size distribution
Low computational cost and relatively accurate distribution
Complex source termThe 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](1)∂ρd∂t+∂ρdudα∂xα=0.The average velocities of the disperse phases are calculated by solving the corresponding momentum balance equations [25]:(2)∂ρdudα∂t+∂ρdudαudβ∂xα=−∂p∂xα+∂μd∂udα/∂xβ+∂udβ/∂xα∂xβ+ρdgα+Fα,where ud is the disperse phase velocity; ρd is the disperse phase density; Fα represents the interphase interaction; and αandβ represent the αandβ directions. Equations (3) and (4) give the transport equation for the turbulent kinetic energy term (k) and the transport equation for the turbulent dissipation term (ε), respectively [2]:(3)∂ρdk∂t+∂ρdkudα∂xα=∂∂xβρCμk2εσk∂k∂xβ+Gk−ρε+Gb−YM+Sk,(4)∂ρdε∂t+∂ρdεudα∂xα=∂∂xαρdCμk2εσε∂ε∂xα+εkC1εGk−C2ερdε+Sε.According to Launder and Spalding [26] recommended values, the constants in above models are equal to the numerical values below [27]:(5)C1ε=1.44,C2ε=1.92,Cμ=0.09,σk=1.0,σε=1.3.The energy conservation equation for the disperse phase [17] can be expressed as(6)32∂εdρdΘ∂t+∂εdρdudαΘ∂xα=−p+μd∂udα∂xαI¯¯+2μdτd¯¯:∇udα+∇⋅κd∇Θ−γΘ+φdd+φcd,where Θ is the granular temperature; εd is the disperse phase volume fraction; p is the disperse phase pressure; γΘ is the collisional rate of energy dissipation; κd is the droplet fluctuating energy; φdd denotes a dissipation term due to the interphase interaction; and φcd 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:(7)p=εdρdΘ21+egoεd+1,κd=εdρdΘ1+ego651+egoεd+12+2εd2ρddd1+egoΘπ1/2,γΘ=31−e2εd2ρdgoΘ4ddΘπ1/2−∂udα∂xα.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]:(8)Θ=12CxCx+CyCy,CxCx=1MN∑αM∑βNud,αβx→,t−udx→ud,αβx→,t−udx→,udx→=1MN∑αM∑βNud,αβx→,t.The interaction between two phases is considered with the interphase interaction (Fα), which includes lift force, turbulent diffusion force, drag force, and virtual mass force (Fα=Mlα+MTα+Mdα+MVα). The lift force (Mlα) 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] (aα represent the accelerations of two phases):(9)MVα=−CVMρcεdaαd−aαc.(2)
The drag force coefficientCD has a variety of models. The following Schiller-Naumann-Pb equation is employed for the drag force estimation [31]:(10)Mdα=34CDεdεcρddducα−udαεc−2.65,CD=24Re1+0.15Reb0.687+0.4131+16.3Re−1.09,Re=ρcucα−udαddμc.(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]:(11)MTα=CDμcPr∇εcεc−∇εdεd.


### 2.2. Population Balance Model

#### 2.2.1. Balance Model

The PBM is based on the conservation equation of number density functionNDF [32]. A general form of PBM can be expressed as(12)∂nL;x→,t∂t+∂nL;x→,tuα∂xα=−∂∂LGLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where(13)BagL;x→,t=12∫0∞nλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,DagL;x→,t=∫0∞Lknλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,BbrL;x→,t=∫0∞Lk∫0∞aλbLλnλ;x→,tdLdλ,DbrL;x→,t=∫0∞Lkaλnλ;x→,tdL.The moments of method are defined as follows [33]:(14)mkx→,t=∫0∞LknL;x→,tdL.Substituting Equations (14) into (12)(15)∂mk∂t+∂mkuα∂xα=∫0∞kLk−1GLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where nL;x→,t is the number density function with droplet diameter L; GLnL;x→,t is the droplet flux due to growth rate; x→ is the droplet velocity vector; BagL;x→,t and DagL;x→,t are the birth and death rates of the droplets for aggregation; and BbrL;x→,t and DbrL;x→,t 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 frequencygV′ and the daughter PDF βV,V′, where(16)gV′=C2ε1/3erfcC3σρcε2/3d5/3+C4μcρcρdε2/3d5/3,βV,V′=30V′VV′21−VV′2.The constantsC2 = 2.52, C3 = 0.04, and C4 = 0.01. V and V′ 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 Vi and Vj:(17)ΩagVi,Vj=ωagVi,VjPagVi,Vj,ωagVi,Vj=π4di+dj2ninjuij¯,Pag=exp−cI0.751+xij21+xij31/2ρd/ρc+0.51/21+xij3Weij1/2,xij=di/dj.


### 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. Figure1 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].Figure 1
CFD-PBM coupled model.


### 2.5. Numerical Methods

#### 2.5.1. CM

The CM divides the complete region of droplet size into several classes, where each classi corresponds to the size interval defined as [dp,dpi+1] [6]. The density function within class i and numbers of droplets Ni can be represented as(18)ndp,t=∑inidp,t=∑iNiδdp−dpi.The Sauter mean diameter (d32) is usually used as the mean droplet size. Sauter average size in extractor can be expressed as follows:(19)d32=∑idi3ni∑idi2ni.


#### 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:(20)mkx→,t=∫0∞fLnLdL=∫0∞nL;x→,tLkdL≈∑i=1NfLiωi.By applying QMOM, the transport equation for the moment ofkth is obtained as follows:(21)∂mk∂t+∂uαmk∂xα=k∑i=1NLik−1GiLiωi+12∑i=1Nωi∑j=1NωiLi3+Lj3k/3βLi,Lj−∑i=1NLikωi∑j=1NωiβLi,Lj+∑i=1Nωi∫0∞LkaLibLLidL−∑i=1NLikωiaLi.


#### 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 (ωi) and abscissas (ξi=Liωi) with respect to the two properties:(22)∂ωi∂t+∂ωiuα∂xα=ai,∂ξi∂t+∂ξiuα∂xα=bi.The above linear system can be written in the matrix form as(23)Aα=d,where the coefficient matrix A=A1A2 can be expressed as follows:(24)A1=1...10...0:...:−L12...−LN221−NL12N−1...21−NL12N−1,A2=1...10...02L1...LN:...:21−NL12N−2...21−NL12N−2.The2N vector of α can be represented as(25)α=a1...aNba...bNT=ab.The right side of Equation (23) is the source term involving droplet aggregation and breakage. The growth term is accounted directly in Equation (23):(26)d=S¯0N…S¯2N−1NT.The source term forkth the moment S¯kN (k = 0, ..., 2N − 1) can be expressed as(27)S¯kNx→,t=∫0∞LkSx→,tdL.


### 2.6. Reactor and Simulation Description

The geometry and computational meshes of the SFBE are shown in Figure2. The simulation is based on the axisymmetric model, with the mesh size of 1 mm. The following Tables 2–4 are the required boundary conditions, the model parameters, and the numerical schemes of simulation.Figure 2
(a) The equipment drawing of investigated spray fluidized bed extractor. (b) CFD grid.
(a)
(b)Table 2
Numerical schemes.
Variable
Schemes
Momentum
First-order upwind
Volume fraction
Quick
k
Second-order upwind
ε
Second-order upwind
Transient formulation
Second-order implicitTable 3
Boundary parameter.
Parameter
Value
Breakage kernel
Laakkonen
Aggregation kernel
Luo
Surface tension
0.05
Granular bulk viscosity
Lun
Granular viscosity
Syamlal O’Brien
Packing limit
0.63
Radial distribution
Lun
Solid pressure
Lun
Number of time step
1000
Time step
0.001 sTable 4
Model parameter.
Parameter
Value
Inlet condition
Velocity inlet
Outlet condition
Pressure outlet
Wall condition
No slip for oil
Free slip for liquid
Inlet turbulent intensity
0.0064 m2s−2
Inlet hydraulic diameter
0.01 m
Operating density
800 kgm−3
Continuous phase density
998 kgm−3
Disperse phase density
800 kgm−3
Granular temperature
0.0001 m2s−2
Inlet velocity
0.9 ms−1
Diameter
Sauter
Continuous phase viscosity
0.0012 kgm−1s−1
Disperse phase viscosity
0.0089 kgm−1s
Back volume fraction
1
Operating pressure
101325 pa
Inlet turbulent dissipation rate
0.009 m2s−2
Inlet turbulent kinetic energy
0.24 m2s−3


## 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](1)∂ρd∂t+∂ρdudα∂xα=0.The average velocities of the disperse phases are calculated by solving the corresponding momentum balance equations [25]:(2)∂ρdudα∂t+∂ρdudαudβ∂xα=−∂p∂xα+∂μd∂udα/∂xβ+∂udβ/∂xα∂xβ+ρdgα+Fα,where ud is the disperse phase velocity; ρd is the disperse phase density; Fα represents the interphase interaction; and αandβ represent the αandβ directions. Equations (3) and (4) give the transport equation for the turbulent kinetic energy term (k) and the transport equation for the turbulent dissipation term (ε), respectively [2]:(3)∂ρdk∂t+∂ρdkudα∂xα=∂∂xβρCμk2εσk∂k∂xβ+Gk−ρε+Gb−YM+Sk,(4)∂ρdε∂t+∂ρdεudα∂xα=∂∂xαρdCμk2εσε∂ε∂xα+εkC1εGk−C2ερdε+Sε.According to Launder and Spalding [26] recommended values, the constants in above models are equal to the numerical values below [27]:(5)C1ε=1.44,C2ε=1.92,Cμ=0.09,σk=1.0,σε=1.3.The energy conservation equation for the disperse phase [17] can be expressed as(6)32∂εdρdΘ∂t+∂εdρdudαΘ∂xα=−p+μd∂udα∂xαI¯¯+2μdτd¯¯:∇udα+∇⋅κd∇Θ−γΘ+φdd+φcd,where Θ is the granular temperature; εd is the disperse phase volume fraction; p is the disperse phase pressure; γΘ is the collisional rate of energy dissipation; κd is the droplet fluctuating energy; φdd denotes a dissipation term due to the interphase interaction; and φcd 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:(7)p=εdρdΘ21+egoεd+1,κd=εdρdΘ1+ego651+egoεd+12+2εd2ρddd1+egoΘπ1/2,γΘ=31−e2εd2ρdgoΘ4ddΘπ1/2−∂udα∂xα.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]:(8)Θ=12CxCx+CyCy,CxCx=1MN∑αM∑βNud,αβx→,t−udx→ud,αβx→,t−udx→,udx→=1MN∑αM∑βNud,αβx→,t.The interaction between two phases is considered with the interphase interaction (Fα), which includes lift force, turbulent diffusion force, drag force, and virtual mass force (Fα=Mlα+MTα+Mdα+MVα). The lift force (Mlα) 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] (aα represent the accelerations of two phases):(9)MVα=−CVMρcεdaαd−aαc.(2)
The drag force coefficientCD has a variety of models. The following Schiller-Naumann-Pb equation is employed for the drag force estimation [31]:(10)Mdα=34CDεdεcρddducα−udαεc−2.65,CD=24Re1+0.15Reb0.687+0.4131+16.3Re−1.09,Re=ρcucα−udαddμc.(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]:(11)MTα=CDμcPr∇εcεc−∇εdεd.


## 2.2. Population Balance Model

### 2.2.1. Balance Model

The PBM is based on the conservation equation of number density functionNDF [32]. A general form of PBM can be expressed as(12)∂nL;x→,t∂t+∂nL;x→,tuα∂xα=−∂∂LGLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where(13)BagL;x→,t=12∫0∞nλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,DagL;x→,t=∫0∞Lknλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,BbrL;x→,t=∫0∞Lk∫0∞aλbLλnλ;x→,tdLdλ,DbrL;x→,t=∫0∞Lkaλnλ;x→,tdL.The moments of method are defined as follows [33]:(14)mkx→,t=∫0∞LknL;x→,tdL.Substituting Equations (14) into (12)(15)∂mk∂t+∂mkuα∂xα=∫0∞kLk−1GLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where nL;x→,t is the number density function with droplet diameter L; GLnL;x→,t is the droplet flux due to growth rate; x→ is the droplet velocity vector; BagL;x→,t and DagL;x→,t are the birth and death rates of the droplets for aggregation; and BbrL;x→,t and DbrL;x→,t are the birth and death rates of droplets for breakage, respectively.


## 2.2.1. Balance Model

The PBM is based on the conservation equation of number density functionNDF [32]. A general form of PBM can be expressed as(12)∂nL;x→,t∂t+∂nL;x→,tuα∂xα=−∂∂LGLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where(13)BagL;x→,t=12∫0∞nλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,DagL;x→,t=∫0∞Lknλ;x→,t∫0∞βλ,Lnλ;x→,tdLdλ,BbrL;x→,t=∫0∞Lk∫0∞aλbLλnλ;x→,tdLdλ,DbrL;x→,t=∫0∞Lkaλnλ;x→,tdL.The moments of method are defined as follows [33]:(14)mkx→,t=∫0∞LknL;x→,tdL.Substituting Equations (14) into (12)(15)∂mk∂t+∂mkuα∂xα=∫0∞kLk−1GLnL;x→,t+BagL;x→,t−DagL;x→,t+BbrL;x→,t−DbrL;x→,t,where nL;x→,t is the number density function with droplet diameter L; GLnL;x→,t is the droplet flux due to growth rate; x→ is the droplet velocity vector; BagL;x→,t and DagL;x→,t are the birth and death rates of the droplets for aggregation; and BbrL;x→,t and DbrL;x→,t 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 frequencygV′ and the daughter PDF βV,V′, where(16)gV′=C2ε1/3erfcC3σρcε2/3d5/3+C4μcρcρdε2/3d5/3,βV,V′=30V′VV′21−VV′2.The constantsC2 = 2.52, C3 = 0.04, and C4 = 0.01. V and V′ 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 Vi and Vj:(17)ΩagVi,Vj=ωagVi,VjPagVi,Vj,ωagVi,Vj=π4di+dj2ninjuij¯,Pag=exp−cI0.751+xij21+xij31/2ρd/ρc+0.51/21+xij3Weij1/2,xij=di/dj.


## 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. Figure1 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].Figure 1
CFD-PBM coupled model.


## 2.5. Numerical Methods

### 2.5.1. CM

The CM divides the complete region of droplet size into several classes, where each classi corresponds to the size interval defined as [dp,dpi+1] [6]. The density function within class i and numbers of droplets Ni can be represented as(18)ndp,t=∑inidp,t=∑iNiδdp−dpi.The Sauter mean diameter (d32) is usually used as the mean droplet size. Sauter average size in extractor can be expressed as follows:(19)d32=∑idi3ni∑idi2ni.


### 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:(20)mkx→,t=∫0∞fLnLdL=∫0∞nL;x→,tLkdL≈∑i=1NfLiωi.By applying QMOM, the transport equation for the moment ofkth is obtained as follows:(21)∂mk∂t+∂uαmk∂xα=k∑i=1NLik−1GiLiωi+12∑i=1Nωi∑j=1NωiLi3+Lj3k/3βLi,Lj−∑i=1NLikωi∑j=1NωiβLi,Lj+∑i=1Nωi∫0∞LkaLibLLidL−∑i=1NLikωiaLi.


### 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 (ωi) and abscissas (ξi=Liωi) with respect to the two properties:(22)∂ωi∂t+∂ωiuα∂xα=ai,∂ξi∂t+∂ξiuα∂xα=bi.The above linear system can be written in the matrix form as(23)Aα=d,where the coefficient matrix A=A1A2 can be expressed as follows:(24)A1=1...10...0:...:−L12...−LN221−NL12N−1...21−NL12N−1,A2=1...10...02L1...LN:...:21−NL12N−2...21−NL12N−2.The2N vector of α can be represented as(25)α=a1...aNba...bNT=ab.The right side of Equation (23) is the source term involving droplet aggregation and breakage. The growth term is accounted directly in Equation (23):(26)d=S¯0N…S¯2N−1NT.The source term forkth the moment S¯kN (k = 0, ..., 2N − 1) can be expressed as(27)S¯kNx→,t=∫0∞LkSx→,tdL.


## 2.5.1. CM

The CM divides the complete region of droplet size into several classes, where each classi corresponds to the size interval defined as [dp,dpi+1] [6]. The density function within class i and numbers of droplets Ni can be represented as(18)ndp,t=∑inidp,t=∑iNiδdp−dpi.The Sauter mean diameter (d32) is usually used as the mean droplet size. Sauter average size in extractor can be expressed as follows:(19)d32=∑idi3ni∑idi2ni.


## 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:(20)mkx→,t=∫0∞fLnLdL=∫0∞nL;x→,tLkdL≈∑i=1NfLiωi.By applying QMOM, the transport equation for the moment ofkth is obtained as follows:(21)∂mk∂t+∂uαmk∂xα=k∑i=1NLik−1GiLiωi+12∑i=1Nωi∑j=1NωiLi3+Lj3k/3βLi,Lj−∑i=1NLikωi∑j=1NωiβLi,Lj+∑i=1Nωi∫0∞LkaLibLLidL−∑i=1NLikωiaLi.


## 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 (ωi) and abscissas (ξi=Liωi) with respect to the two properties:(22)∂ωi∂t+∂ωiuα∂xα=ai,∂ξi∂t+∂ξiuα∂xα=bi.The above linear system can be written in the matrix form as(23)Aα=d,where the coefficient matrix A=A1A2 can be expressed as follows:(24)A1=1...10...0:...:−L12...−LN221−NL12N−1...21−NL12N−1,A2=1...10...02L1...LN:...:21−NL12N−2...21−NL12N−2.The2N vector of α can be represented as(25)α=a1...aNba...bNT=ab.The right side of Equation (23) is the source term involving droplet aggregation and breakage. The growth term is accounted directly in Equation (23):(26)d=S¯0N…S¯2N−1NT.The source term forkth the moment S¯kN (k = 0, ..., 2N − 1) can be expressed as(27)S¯kNx→,t=∫0∞LkSx→,tdL.


## 2.6. Reactor and Simulation Description

The geometry and computational meshes of the SFBE are shown in Figure2. The simulation is based on the axisymmetric model, with the mesh size of 1 mm. The following Tables 2–4 are the required boundary conditions, the model parameters, and the numerical schemes of simulation.Figure 2
(a) The equipment drawing of investigated spray fluidized bed extractor. (b) CFD grid.
(a)
(b)Table 2
Numerical schemes.
Variable
Schemes
Momentum
First-order upwind
Volume fraction
Quick
k
Second-order upwind
ε
Second-order upwind
Transient formulation
Second-order implicitTable 3
Boundary parameter.
Parameter
Value
Breakage kernel
Laakkonen
Aggregation kernel
Luo
Surface tension
0.05
Granular bulk viscosity
Lun
Granular viscosity
Syamlal O’Brien
Packing limit
0.63
Radial distribution
Lun
Solid pressure
Lun
Number of time step
1000
Time step
0.001 sTable 4
Model parameter.
Parameter
Value
Inlet condition
Velocity inlet
Outlet condition
Pressure outlet
Wall condition
No slip for oil
Free slip for liquid
Inlet turbulent intensity
0.0064 m2s−2
Inlet hydraulic diameter
0.01 m
Operating density
800 kgm−3
Continuous phase density
998 kgm−3
Disperse phase density
800 kgm−3
Granular temperature
0.0001 m2s−2
Inlet velocity
0.9 ms−1
Diameter
Sauter
Continuous phase viscosity
0.0012 kgm−1s−1
Disperse phase viscosity
0.0089 kgm−1s
Back volume fraction
1
Operating pressure
101325 pa
Inlet turbulent dissipation rate
0.009 m2s−2
Inlet turbulent kinetic energy
0.24 m2s−3


## 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.Figure 3
Droplet Sauter diameter distribution along (a) the radical distance (h = 0.2 m) and (b) the height in the extractor.
(a)
(b)


### 3.2. Spray Process

Figure4 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.Figure 4
Comparison of the spray process from 0.1s to 1s using (a) CM. (b) DQMOM and (c) QMOM with the (d) the simulation results predicted by Bright A.


### 3.3. Evolution of Mixture Density

The time evolutions for the mixture density predicted with three methods are compared in Figure5. 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.Figure 5
Comparison of the time evolution of the density of mixture using three methods.


### 3.4. Velocity Distribution

Figure6 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.Figure 6
The velocity of disperse phase along the radical distance (a)h = 0.04 m and (b) h = 0.2 m in the extractor.
(a)
(b)


### 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 Figure7. 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.Figure 7
The plunging pressure along the radical distanceh = 0.04 m on the nozzle outlet section.


#### 3.5.2. Granular Temperature Distribution

Figure8 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.Figure 8
The granular temperature in the axial direction along the nozzle zone.


### 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 Figure9. 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 ΔPd=ρd−ρc1−ηgH, 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.Figure 9
The droplets dynamic pressure along (a) the height and (b) the radical distanceh = 0.2 m in the extractor.
(a)
(b)


#### 3.6.2. Mixture Solution Turbulent Properties

In Figure10, 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).Figure 10
The mixture phase turbulent performance along the radical distanceh = 0.2 m.
(a)
(b)


#### 3.6.3. Disperse Phase Volume Fraction

Figure11 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.Figure 11
The disperse phase volume along the radical direction at (a)h = 0.09 m and (b) h = 0.2 m in the extractor.
(a)
(b)


## 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.Figure 3
Droplet Sauter diameter distribution along (a) the radical distance (h = 0.2 m) and (b) the height in the extractor.
(a)
(b)


## 3.2. Spray Process

Figure4 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.Figure 4
Comparison of the spray process from 0.1s to 1s using (a) CM. (b) DQMOM and (c) QMOM with the (d) the simulation results predicted by Bright A.


## 3.3. Evolution of Mixture Density

The time evolutions for the mixture density predicted with three methods are compared in Figure5. 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.Figure 5
Comparison of the time evolution of the density of mixture using three methods.


## 3.4. Velocity Distribution

Figure6 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.Figure 6
The velocity of disperse phase along the radical distance (a)h = 0.04 m and (b) h = 0.2 m in the extractor.
(a)
(b)


## 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 Figure7. 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.Figure 7
The plunging pressure along the radical distanceh = 0.04 m on the nozzle outlet section.


### 3.5.2. Granular Temperature Distribution

Figure8 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.Figure 8
The granular temperature in the axial direction along the nozzle zone.


## 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 Figure7. 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.Figure 7
The plunging pressure along the radical distanceh = 0.04 m on the nozzle outlet section.


## 3.5.2. Granular Temperature Distribution

Figure8 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.Figure 8
The granular temperature in the axial direction along the nozzle zone.


## 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 Figure9. 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 ΔPd=ρd−ρc1−ηgH, 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.Figure 9
The droplets dynamic pressure along (a) the height and (b) the radical distanceh = 0.2 m in the extractor.
(a)
(b)


### 3.6.2. Mixture Solution Turbulent Properties

In Figure10, 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).Figure 10
The mixture phase turbulent performance along the radical distanceh = 0.2 m.
(a)
(b)


### 3.6.3. Disperse Phase Volume Fraction

Figure11 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.Figure 11
The disperse phase volume along the radical direction at (a)h = 0.09 m and (b) h = 0.2 m in the extractor.
(a)
(b)


## 3.6.1. Droplet Dynamic Pressure Distribution

Comparisons of droplet dynamic pressure profiles between the height direction and the radical distance are presented in Figure9. 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 ΔPd=ρd−ρc1−ηgH, 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.Figure 9
The droplets dynamic pressure along (a) the height and (b) the radical distanceh = 0.2 m in the extractor.
(a)
(b)


## 3.6.2. Mixture Solution Turbulent Properties

In Figure10, 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).Figure 10
The mixture phase turbulent performance along the radical distanceh = 0.2 m.
(a)
(b)


## 3.6.3. Disperse Phase Volume Fraction

Figure11 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.Figure 11
The disperse phase volume along the radical direction at (a)h = 0.09 m and (b) h = 0.2 m in the extractor.
(a)
(b)


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



---

*Source: 4836213-2019-03-03.xml*