Hydrodynamic Behavior Analysis of Agitated-Pulsed Column by CFD-PBM

ABSTRACT Recently, a newly designed type of energy-input extraction column, agitated-pulsed column (APC), has been achieved high mass-transfer efficiency due to the small-size droplets and high dispersed-phase holdup. In view of the lack of feasible population-balance-model (PBM) kernel functions in APC, parameter optimization is conducted by a simplified PBM method. Then the optimized PBM kernel functions are implemented in the computational fluid-dynamics (CFD) code to investigate local two-phase flow behaviors in a 25 mm APC. The results show that the CFD-PBM successfully predicts the drop-size distribution measured in the experiments. CFD-PBM also gives good prediction of Sauter mean diameter and dispersed-phase holdup. The local flow behaviors are illustrated to understand the effects of operating conditions on the hydrodynamic performance. This work demonstrates a possibility for prediction of drop-size distribution by the combination of simplified PBM and CFD simulation.


Introduction
Solvent extraction is an important separation technique for a range of industries, in which one liquid phase is dispersed as droplets in another phase by gravity and energy input. [1] Various types of solvent extraction columns have been widely used in nuclear, hydrometallurgy, and chemical sectors. [2] Because of the complexity of two-phase interaction, the mass-transfer performance inside is significantly affected by the energy input such as stirring, [3] rotation, [4] oscillation, [5] and pulsation. [6,7] Recently, a newly designed type of energy-input extraction column, agitated-pulsed column (APC), [8] has been successfully used in separation of chiral molecules [9,10] and 5-hydroxymethylfurfural, [11] which exhibits excellent mass-transfer characteristics. It is measured to be 17 and 25 theoretical stages per meter for the solute transfer from direction dispersed-to-continuous phases and continuous-to-dispersed phases, respectively. The high masstransfer efficiency is attributed to that the reduction of droplet diameter leads to a high dispersed-phase holdup and thus increases the interfacial area of contact between the phases.
Knowledge of the hydrodynamic parameters including drop size and holdup is crucial for design of liquid-liquid extraction columns as they are related to the interfacial area for mass transfer and the allowable throughputs. [12,13] Recently, computational fluid dynamics (CFD) has been considered as a powerful mathematical tool to performance modeling and help in scale-up of extractor. [14][15][16][17][18][19][20][21][22] Moreover, the coupling of CFD and populationbalance model (PBM) provides an opportunity to reproduce the two-phase flow characteristics more realistically. [23] The main challenge of the application of PBM is to determine its kernel functions including the drop-breakage frequency, daughter droplet-size distribution, and coalescence rates.
In the past decades, different kernel functions have been derived theoretically to express the behavior of droplet breakup and coalescence. [24][25][26][27] Although several kernel functions have been applied to the CFD-PBM simulation of two-phase flow in extraction columns, [28,29] validation was done only for Sauter mean diameter, which may be insufficient to identify which one is feasible in CFD-PBM simulation. [30] Considering there is a lack of feasible PBM kernel functions in the PDDC, Zhou et al. [31] developed a direct measurement method for the breakage of droplets by directly counting the breakup frequency in the real contactor. Empirical kernel functions were established for PBM. However, different parameters of kernel functions may be required for different internal structure and materials where expensive experiments are inevitable. [32] Thus, development of an inverse populationbalance model that can optimize the parameters of the existing models for extraction column is necessary. [33,34] In order to solve the inverse population-balance problem, a simplified PBM method was developed in previous studies. [35] Thus, the main objectives of this study are to combine this method with CFD-PBM simulation and get a deeper understanding of the flow characteristics in the APC, which can be helpful for optimization of operation as well as design of this type of extractor.
In this study, CFD-PBM simulations are conducted to investigate the local hydrodynamics and drop dynamics in APC. The PBM kernel functions obtained from literatures are improvement by a simplified PBM method. The drag law is modified to account for the effect of turbulence. Simulation results of drop-size distribution are presented and compared to the experimental measurements of our previous work. The results of simulated dispersed-phase holdup, x d , and Sauter mean drop diameter, d 32 are also used to compare with experimental data for verification. The local variations of liquid-liquid flow characteristics are discussed with simulation results as well as profiles of experimental dispersed-phase holdup.

Euler-Euler mode
In the present work, Eulerian-Eulerian model is used to describe the twophase flow, which assumes two phases as continuous medium by introducing the concept of local volume fraction, α. The conservation equations for the ith phase are given as: Continuity equation: Momentum conservation equation: where u * is phase velocity, ρ is density, τ represents the stress-strain tensor, p represents the pressure shared by all phases, g is gravitational acceleration, F * i is the external body force, F * i;lift is the lift force, F * i;vm the virtual mass force acting on the ith phase. R * ij is the interaction force between ith and jth phase. The subscript i represents phase i.
Because the two-phase condition is used in this study (the aqueous phase is set as the dispersed-phase d, and the organic phase is the continuous phase, c), volume fractions of the two phases should be satisfied as: Among the interaction force R * ij , drag force is the dominant force for liquidliquid interactions, while the virtual mass force and the lift force can be neglected. [20,36,37] The drag force term between the continuous and dispersed phase is defined in terms of the interphase exchange coefficient (F c,d ) as: in which the interphase exchange coefficient is calculated as: where d 32 is the Sauter mean diameter of the dispersed phase, which is calculated as: n i is the number density of drops of size d i . C D is the drag coefficient related to the relative Reynolds number Re: For the drag-coefficient calculating, the standard correlation of Schiller and Naumann is used [38 This basic drag correlation predicts droplets moving in a still liquid well but may not be suitable for droplets moving in turbulent liquid. [39] Thus, a modified drag law is used to take into account the effect of turbulence. [40] Re ¼ where C is an empirical introduced to represent the reduction of relative velocity due to turbulent flow. [14,39,40] The suggested value of 0.5 is set to C. [19,41] The variable μ t;m is the turbulent viscosity, which can be calculated from Equation (15).

Turbulence model
The multiphase turbulence in this study is described by the realizable k-ε mixture turbulence model, which includes solving the partial differential equations separately for turbulence kinetic energy k and its dissipation rate ε.
where S represents the strain rate tensor, G k accounts for the effect of mean velocity gradients in generation of turbulence kinetic energy, and G b represents the effect of buoyancy in generation of turbulence kinetic energy. The subscript m represents the mixture of the two phases. The mixture averaged values of the density and velocity are calculated from:

Population-balance model
The PBM describes the evolution of number-density of different drop sizes with source terms accounting for droplet breakage and coalescence and is given as: . u * and Γ t are the mean droplet velocity and diffusion coefficient of the number density, respectively. S d; t ð Þ is the source term accounting for the coalescence and breakup of the droplets and is generally expressed as: where D B d; t ð Þ and B B d; t ð Þ are, respectively, the death and birth rates due to breakage. D C d; t ð Þ and B C d; t ð Þ are, respectively, the death and birth rates due to coalescence. The general equations of the birth and death rates are given as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where g d ð Þ is the breakup rate and β d; d 0 ð Þ is the distribution of daughter droplets. λ d; d 0 ð Þ and h d; d 0 ð Þ are the coalescence efficiency and collision frequency, respectively, between drops of diameter d and d 0 : The product of two terms is defined as coalescence rate: In the literature, several kernels have been derived theoretically. One of the fundamental and widely accepted models for droplet breakup and coalescence is due to Coulaloglou and Tavlarides [25] :

Numerical procedure
The schematic representation of the experimental apparatus is shown in Figure 1(a) and the key dimensions of the pilot-plant column are shown in Table 1. The organic phase (continuous phase) used in this study is 30% (tributyl phosphate, supplied by Myer, Shanghai) in Shellsoll 2046. The aqueous phase (dispersed phase) used is tap water. Table 2 shows the physical properties of the experimental system used in this study.
A simulation domain of the laboratory-scale APC is designed as shown in Figure 1(b) to conduct the augmented CFD-PBM simulations. In this work, four agitated cells have been used to keep the computational time within reasonable limits. The 3D periodic grid is used because of the rotational symmetry of plates and stirred. A section part of 60 degrees is designed as shown in Figure 1(b). The two sides of the calculation domain are used as periodic boundary condition. The mesh is constructed with quadrilateral grids. A mesh size of 0.5 mm was used in the effective domain of column, and a mesh size of 1 mm was used at the inlet and outlet of the domain. The mesh independent test was conducted (Supplemental Materials, Figure S1)   and mesh sizes of 251,645 cells are found to be optimum and used in the final simulations. All internal walls are set with a no-slip condition, and the nearwall region is modelled with standard wall function. The boundary conditions are properly set as shown in Figure 1(b). An annular domain is created around the impeller to simulate the agitation of impeller as per the moving reference frame (MRF) model. Moving-wall boundary condition is given for the wall of shaft. A pressure outlet condition is defined at the top boundary for organic outlet. The aqueous inlet, aqueous outlet, and organic inlet boundaries are defined as velocity inlet boundaries. A sinusoidal normal velocity is superposed on the inlet of organic phase to represent pulsation input that is given as where V c represents the flow rate of continuous phase, V p represents the pulsation velocity calculated on its amplitude A and frequency f. Although a universal size range has no remarkable influence on steady simulation results, appropriate droplet-size range for each condition would give a better simulation result. [30] The advisable drop-size range for each operating condition is estimated from experimental results and is listed in Table 3. Ten bins are used to represent the drop-size range. The initial dropsize distribution of dispersed-phase inlet is set as 100% maximum diameter. The exponent, q, is used in the discretization of the growth-term diameter coordinate: Simulations are performed in double-precision mode using commercial CFD software ANSYS STUDENT FLUENT 2019R2. The pressure-based solver with implicit formulation is used to numerically solve the model equations. A coupled SIMPLE scheme is chosen to realize pressure-velocity coupling. The time step size used in the model is 0.002 s, while all governing equations are discretized a second-order upwind method. The simulation results are considered to be converged only when all the residuals of the equations are less than 10 -4 .

Modified of kernel functions
Considering there is a lack of feasible droplet-breakage and coalescencefrequency functions in the APC, a simplified PBM method is developed to optimize the parameters of PBM kernel functions. [35] In the regression of parameters, the empirical correlation of energy dissipation rate ε is important, as it directly affects the droplet-breakage and coalescence frequency. So, precise energy dissipation rate ε correlation must be determined before the optimized PBM kernel functions can be implemented in CFD simulation. The average energy dissipation rate under different operation conditions obtained by CFD simulations is presented in Figure 2. Due to the inherent transient nature, simulation results are obtained by the time-average over a periodic cycle. To predict energy dissipation rate using this type of column, correlation is proposed based on Kumar and Hartland [42] and Coulaloglou and Tavlarides [25] : The average absolute relative deviation (AARD) of the calculated values by applying Equation (29) to the simulated results is used to evaluate the accuracy of the prediction. The parameters suggested by Kumar and Hartland [42] and Coulaloglou and Tavlarides [25] result in an AARD of 59.1%, which may be owed to the difference of internal structure and materials between this work and the previous work. In order to obtain a better prediction, the parameters in Equation (29) are refitted to reduce the AARD to about 1.87%. The correlation is given in Equation (31), and the comparison of the simulated energy dissipation rate with the predicted results is shown in Figure 2.
Based on Equation (31), the parameters of Coulaloglou and Tavlarides [25] are refitted with the simplified PBM proposed in our previous studies. [35] The sensitivity analysis of Coulaloglou and Tavlarides model to adjusted parameters is presented in Figure S2 (Supplemental Materials). The sketch of the regression with simplified PBM model is shown in Figure 3, and the optimized parameters are listed in Table 4. It can be seen that the parameters vary with types of column and physical properties (Supplemental Materials, Figure S3). Thus, one needs to choose the most appropriate parameters for a given system.

Validation and comparison
For validation, the simulated cumulative volume droplet-size distributions (DSD) are compared with our previous experimental results measuring at the middle part of column. [35] It can be seen from Figure 4 that an acceptable agreement in the DSD between CFD-PBM simulation and experimental data is obtained, which indicates the reliability of regressed parameters to calculate the drop-size distribution in APC. DSD predicted by the simplified PBM are observed to be more precise than that of CFD-PBM simulation. The deviations between CFD-PBM simulation and experimental results might be attributed to the assumption of homogeneous distribution of turbulent energy-dissipation rate in the regression of parameters differs from the real conditions. It also can be seen from Figure 4(d) that the proportion of droplets with small diameters simulated from CFD-PBM is overestimated compared with experimental data. The reason is probably that the non-homogeneous distribution of turbulent energy-dissipation rate shown in Figure 5 leads to energy focusing and enhances the breakup of smaller droplets. The dispersed holdup and Sauter mean drop diameter, d 32 , are also obtained by the droplet-size distribution to compare with the experimental results obtained in our previous work [35,44] It can be seen from Figure 6, which indicates that the simulated dispersed holdup and Sauter mean diameter with CFD-PBM have the same trend with the experimental results. Figure 7 depicts a more systematic comparison between experimental and simulated results, in which the AARD is always within ±15% for dispersed holdup and ±20% for Sauter mean drop diameter, respectively. The average absolute relative error, considering all data points, is found to be about 10.84% for dispersed-phase holdup and 9.40% for Sauter mean drop diameter, which suggests that the established computational approach is reliable to simulate the two-phase flow behaviors in APC.

Local behaviors of two-phase flow
The local behaviors of two-phase flow can provide useful insights into twophase flow characteristics in the APC, which is important to understand the effect of operation conditions on hydrodynamic performance. Figure 8 depicts the velocity field of the dispersed phase during a pulsation period at agitation speed of 300, 600 and 900 rpm. With increase of agitation speed, the flow gradient field becomes more intensive and recirculation appears in the dispersed phase (Figure 9). The maximum velocity is clearly observed at the impeller tip.
At high agitation speed (900 rpm), two big vortexes one above and one underneath the impeller are observed, which extend the residence time of dispersed-phase droplets leads to an increase of dispersed-phase holdup as well as axial dispersion. These are in agreement with the PIV-measurement results presented by Hlawitschka and Bart [45] with a Kuhni column. At the negative peak moment of the pulsation cycle (t = 0.25 s), the vortex underneath the impeller is stronger than the vortex above impeller, which is the inverse of that at the positive peak moment (t = 0.75 s). It is noted that as the agitation speed increases, there are insignificant changes in vortex scale, while the intensity is significantly enhanced.
For APC, droplet coalescence and breakage rates are closely linked to energy-dissipation rate ε. Figure 10 depicts the contours of energydissipation rate during one pulsation period at agitation speed of 300, 600, and 900 rpm. It can be seen that, with the increase of agitation speed, the turbulent energy-dissipation rate rapidly increases. The turbulence dissipation inside the APC is not homogenously distributed, and the relatively higher energy-dissipation rate can be observed at the edge of the impellers as well as the apertures between impellers and column walls under higher agitation speed. This agrees well with that of the velocity vector of the dispersed phase shown in Figure 8, which is primarily due to the remarkable changes of velocity magnitude and direction leading to a high turbulence kinetic energy and hence increases the energy-dissipation rate. Figure 10(a) shows the distribution of Sauter mean diameter at agitation speed of 300, 600, and 900 rpm. The droplet size decreases with increasing agitation speed due to the increase of turbulence energy-dissipation rate shown in Figure 10. Furthermore, the droplet distribution becomes more uniform at higher agitation speed, which is in agreement with the experimental results shown in Figure 10(b).
At the negative peak moment of the pulsing cycle (t = 0.25 s), the mean droplet diameter at the top domain is always found larger than that at the lower domain due to the consecutive introducing of maximum size droplets from the dispersed-phase inlet. [30] At the positive peak moment of the pulsing cycle (t = 0.75 s), the droplet swarm is entrained into the vortex, which increases the residence time of dispersed phase and enhances the droplet breakage. So, the local mean diameter during positive pulsation period is smaller than that during negative pulsation period.  The distribution of dispersed-phase volume fractions at agitation speed of 300, 600, and 900 rpm are shown in Figure 11. It can be observed that the dispersed-phase holdup increases with the increase of agitation speed, which is in agreement with the experimental results shown in Figure 10(b). A more uniform distribution of dispersed holdup can be observed at higher agitation speed. This is attributed to the increase of drag force acting on the fine drops.
At the negative peak of the pulsing cycle (t = 0.25 s), the dispersed phase in the vicinity of the plates is carried away with the continuous phase into the following downward cell to prevent flooding of single cells and the entire column. At the positive peak pulsing cycle (t = 0.75 s), the upward flow of continuous phase and plate holes can prevent droplets from moving to the next plate, which causes an accumulation of the dispersed phase on the plate surfaces. Moreover, at high agitation speed (N = 900 rpm), the dispersed phase is involved into the vortex and mostly circles in a single cell of the column, leading to a longer residence time of dispersed phase and hence increases holdup.
The excellent mass-transfer performance in APC is attributed to the enhancement of droplet breakage and coalescence at high agitation speed. Figure 12 shows the droplet breakage and coalescence frequency distribution under high agitation speed (N = 900) during one pulsation period. The required droplet size is set to 0.74 mm, which is the experimental Sauter mean droplet diameter. It can be seen from Figure 12 that the droplet breakage frequency is most intensive at the impeller tip, which coincides with the distribution of the turbulent energy-dissipation rate shown in Figure 9. The maximum value of droplet coalescence frequency locates at the zones between the impellers and plates.

Conclusion
In this study, the optimized PBM kernel functions obtained by a simplified PBM method are implemented in the CFD to investigate the two-phase flow behaviors in the APC. The following conclusions can be drawn from this study: (a) The improvement of breakage and coalescence kernel parameters with simplified PBM is used in CFD-PBM simulation of APC. CFD-PBM results successfully predict the influence of varying pulsation intensity on droplet-size distribution with acceptable deviation, which validates the regressed parameters with the simplified PBM method.
(b) Favorable agreements between calculated results (with simplified PBM) and experimental data indicate the ability of simplified PBM method to predict dropsize distribution in APC, which is a result highly interesting for design purposes.
(c) CFD-PBM was found to give a good prediction for dispersed-phase holdup and Sauter mean diameter in APC with AARD of 10.8% and 9.4%, respectively.
(d) Comparison between the experimental and simulated local flow behaviors has shown the feasibility of the established Euler-Euler. Both local variations of hydrodynamics and velocity field provide valuable insights into the effect of operating conditions on hydrodynamic performance.
In the further work, the effect of physical properties on simulation should be considered, and the mass-transfer model needs to be implemented into CFD-PBM to investigate the mass-transfer performance in APC.