Article pubs.acs.org/IECR
Key Parameter Prediction and Validation for a Pilot-Scale RotatingDisk Contactor by CFD−PBM Simulation Hang Chen, Ze Sun,* Xingfu Song, and Jianguo Yu* National Engineering Research Center for Integrated Utilization of Salt Lake Resource, East China University of Science and Technology, Shanghai 200237, China ABSTRACT: Compared with traditional methods, computational fluid dynamics−population balance model (CFD−PBM) simulations are much more efficient for the design of extractors. In this study, the CFD−PBM framework was established for a pilot-scale rotating-disk contactor. The Euler−Euler approach incorporated with the realizable k−ε model was used for the twophase simulation. For droplet coalescence and breakage, predictive closures based on the model of Luo and Svendsen were employed for further verification. The species transport model was also solved to estimate axial mixing. The model was validated first by comparing the flow field with referenced particle image velocimetry data. Then the predicted key parameters of droplet diameter, dispersed-phase holdup, and Peclet number were compared with empirical correlations for a comprehensive validation. The results indicate that the predictive CFD−PBM is suitable for the design of extraction columns. Its average deviations on the three parameters are 15.7%, 12.9%, and 15.2%, respectively. The model also successfully predicts the regular variations of the droplet diameter and holdup with the rising agitation rate, which contribute to the validation.
1. INTRODUCTION Extraction is an important separation process, and it is widely used in chemical, biochemical, pharmaceutical, and environmental industries.1 Since the first liquid−liquid extraction column was designed, extraction equipment has experienced rapid development.2 Thousands of extraction columns have been constructed for different kinds of separation assignments. The layout of extractors is mainly based on the simplified concept of height of a transfer unit and number of transfer units,3 where expensive experiments are unavoidable for the flow nonidealities. The back and forward mixing models4 also face the same situation. Hence, it is necessary to develop a robust and highly efficient method. Computational fluid dynamics (CFD) has been applied in chemical engineering with great success. It is also used to investigate the hydrodynamics of extraction equipment.5 The investigations showed that the CFD model had the ability to reasonably predict the single-phase flow fields of extractors.6−9 For instance, the single-phase flow of a rotating-disk contactor (RDC) simulated by Fei et al.10 agreed well with laser doppler velocity measurements, and the internals were improved based on the simulation. With regard to multiphase flows, You et al.11 simulated the two-phase flow field of a Kühni column. Yadav et al.12 developed a two-phase CFD model for sieve and pulsedsieve plate extraction columns. The model was validated with holdup empirical correlations. Drumm et al.13,14 successfully developed the two-phase model for RDCs, and the model was validated with particle image velocimetry (PIV) measurements. Onink et al.15 simulated the two-phase flow field of a RDC with high-viscosity extraction systems. In general, the CFD is a useful tool for the hydrodynamic investigations of extraction columns. However, the assumption of a Sauter mean diameter is the main drawback of the two-phase CFD model. It is well-known that extraction column hydrodynamics involves complex droplet coalescence and breakage, which result in droplet-size distributions.16 Hence, the population © 2014 American Chemical Society
balance model (PBM) simulation of extractors has gotten attention. The PBM was introduced into chemical engineering for the description of physical processes involving dispersed particles, droplets, or bubbles.17 In extraction columns, it mainly focuses on the complicated droplet coalescence and breakage. In the past decades, a series of models were proposed for coalescence and breakage. They were reviewed particularly by Liao and Lucas.18,19 Special coalescence and breakage closures for RDCs were also developed by Simon et al.20,21 and Vikhansky et al.22 According to the breakage and coalescence theories, hydrodynamic physical quantities are essential for closure of the PBM. So, Kronberger et al.23 and Attarakih et al.24−29 coupled the experimental correlation-estimated hydrodynamics with the PBM. Although this way has reasonable simulation results,30,31 the establishment of hydrodynamic correlations requires numerous experiments. Hence, Bart32,33 suggested the coupling of CFD and PBM. Both of them can benefit from each other: the PBM helps to overcome the Sauter mean diameter assumption in the CFD model, and the CFD offers the necessary hydrodynamic parameters inversely. Up to now, only a few related researches were reported. GhaniyariBenis et al.34 established the coupled CFD−PBM for a RDC. Only the droplet breakage was considered because of the low dispersed-phase volume fraction. Both Vikhansky et al.35 and Drumm et al.36 simulated a RDC with the coupled CFD−PBM, and the droplet-size distributions were validated. Hlawitschka et al.37,38 employed the CFD−PBM to simulate a Kühni column, and the holdup was compared with that of laser-induced fluorescence measurements. Received: Revised: Accepted: Published: 20013
August 4, 2014 November 14, 2014 December 4, 2014 December 4, 2014 dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
is considered with the interphase forces F, which primarily includes the drag force, lift force, etc. Only the drag force was taken into account because the others are of minor importance.40 The following Schiller−Naumann model was employed for the drag-force estimation. ⎯u − → ⎯u |(→ ⎯ → ⎯ 3ρ αcαdC D|→ ⎯→ ⎯ d c ud − uc) Fc,d = c (4) 4d
Although the coupled CFD−PBM is a powerful tool for the multiphase simulation of extractors, laboratory-scale experiments are still inevitable because current coalescence and breakage models always contain several adjustable parameters. Drumm et al.36 and Chen et al.39 once presented a coalescence and breakage model without adjustable parameters. Further verifications are still necessary to confirm this model’s validity and predictability. Moreover, the existing validations of the CFD−PBM mainly concentrated on droplet-size distributions. Nevertheless, according to the two-film theory, both the interfacial area and concentration difference are significant to extraction columns. The interfacial area is determined by the dispersed-phase holdup and droplet mean diameter. The concentration difference depends on the axial mixing, which can be quantified with the dimensionless Peclet number. Hence, not only the droplet size but also the dispersed-phase holdup and Peclet number should be included for a comprehensive validation of the CFD−PBM. RDC is a typical countercurrent agitated extraction column. Therefore, on the basis of the above two points, the coupled CFD−PBM was established for a RDC. Predictive coalescence and breakage closures36,39 were used. For a comprehensive validation, the simulated flow field was compared with referenced PIV data first. Then the key hydrodynamic parameters of holdup, droplet diameter, and Peclet number were validated with empirical correlations. By the way, the results of a simplex CFD model were also included for comparison.
0.687 ⎧ ⎪ 24(1 + 0.15Re )/Re , Re ≤ 1000 CD = ⎨ ⎪ 0.44, Re > 1000 ⎩
where d is the droplet diameter and CD is the drag coefficient. In the CFD model, the experimental Sauter mean diameter of 2.5 mm was specified.13 For the turbulent simulation, the Reynolds stresses in the momentum equation must be modeled. The realizable k−ε mixture turbulence model was applied. The transport equations for k and ε are briefly described as follows:
2. MATHEMATICAL MODELS 2.1. Two-Phase CFD Model. For the flow-field validation with the PIV data of Drumm et al.,13 the RDC geometry and extraction system are kept the same. The important geometrical parameters of the pilot-scale RDC are listed in Table 1. A 2D axis-symmetric model was established.
size (mm)
column diameter D stator opening DS diameter rotor disk diameter DR compartment height H
150 105 90 30
(1)
∂ ⎯u ) + ∇·(α ρ → ⎯ → ⎯ (αkρk→ k k k uk uk ) ∂t ⎯τ + F ⃗ = −αk∇p + αkρk g ⃗ + ∇·→ k
(2)
where ρ is the fluid density, u is the fluid velocity, p is the pressure, g is the gravitational acceleration, and τ is the stress− strain tensor. The volume fraction αk represents the occurrence of the kth phase, and the sum of the phases has to be 1.
∑ αk = 1
(6)
∂ (ρ ε) + ∇·(ρm ⎯→ u m ε) ∂t m ⎛ μt,m ⎞ ε = ∇·⎜ ∇ε⎟ + ρm C1Sε − ρm C2 k vε σ + ⎝ ε ⎠ ε + C1ε C3εGb k
(7)
k2 ε
(8)
where k is the turbulence kinetic energy, ε is the turbulence dissipation rate, and μt is the turbulent viscosity. As shown in Figure 1, the column geometry was meshed with quadrilateral grids. The optimized mesh size was determined to be 0.5 mm according to the mesh independence results presented later. In view of the numerical stability, the boundary conditions of one pressure outlet and one velocity inlet were applied at the top and bottom of the column, respectively. For the simulation of agitation, the rotors were treated as moving walls, while the stators were regarded as stationary walls. All of the walls were defined with the no-slip wall boundary condition. The SIMPLE algorithm was used for the pressure−velocity coupling. All of the governing equations were discretized with the QUICK scheme. The convergence criteria was set as 10−3 (10−4 for the dispersed-phase volume fraction), and the time step of 0.001 s was chosen. 2.2. PBM. The PBM simulation is based on the conservation equation of a number density function n(x,V,t), which denotes the droplet spatial position x as an external coordinate and the droplet volume V as an internal coordinate. In extraction columns, the growth of droplets is negligible because of the liquid incompressibility. Then the governing equation is expressed as
The Euler−Euler approach was used for the two-phase simulation, where the transport of fluid is described by the continuity and momentum conservation equations as ∂ ⎯u ) = 0 (αkρk ) + ∇·(αkρk→ k ∂t
∂ (ρ k) + ∇·(ρm ⎯→ u mk ) ∂t m ⎛ μt,m ⎞ = ∇·⎜ ∇k ⎟ + Gk ,m + Gb ,m − ρm ε ⎝ σk ⎠
μt,m = ρm Cμ
Table 1. Geometric Sizes of the RDC geometry
(5)
(3)
∂ [n(V , t )] + ∇·[u ⃗ ·n(V , t )] = S(V , t ) ∂t
In the Eulerian framework, the conservation equations are solved for each phase. The interaction between different phases 20014
(9)
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
daughter particle-size distribution was calculated directly from the breakage rate.41 The breakage kernels are presented as g (d ) =
⎛ ε ⎞1/3 1 0.923(1 − α)⎜ 2 ⎟ ⎝d ⎠ 2 ⎛ ⎞ 1 1 (1 + ζ )2 12c f σ ⎜− ⎟ dζ df exp ⎜ v 2/3 5/3 11/3 ⎟ 0 ζmin ζ 11/3 ⎝ β0ρc ε d ζ ⎠
∫ ∫
(15) 1 (1 + ζ )2
2∫ β (f v v , v ) =
ζmin
1
ζ 11/3
1 (1 + ζ )2
v∫ ∫ 0 ζ
min
ζ 11/3
⎛ ⎞ 12c σ exp⎜ − 2/3 f5/3 11/3 ⎟ dζ ⎝ β0ρc ε d ζ ⎠ ⎛ ⎞ 12c σ exp⎜ − 2/3 f5/3 11/3 ⎟ dζ dfv ⎝ β0ρc ε d ζ ⎠ (16)
where σ is the interfacial tension, cf is the interface area coefficient, ζ is the droplet size ratio, and β0 is a numerical constant. It has been mentioned that the coalescence rate is the product of collision frequency h and collision efficiency λ. Usually, the droplet collision mechanism is analogous to collisions between molecules, as in the kinetic theory of gases.43 On the basis of the kinetic gas theory, the collision frequency is interpreted as the effective volume swept by the moving droplet per unit time. The relative velocity between two droplets is proportional to the mean square root of two equivalent eddy velocities, which can be calculated with the classical turbulent theories.19 Then the widely used expression for calculation of the droplet coalescence frequency is presented as follows: π h(d1 , d 2) = (d12 + d 2 2)β0ε1/3(d12/3 + d 2 2/3)1/2 (17) 4
Figure 1. Geometry and grid scheme of the RDC.
where S(V,t) is the source term for droplet coalescence and breakage. It contains the birth and death functions BC(V,t), DC(V,t), BB(V,t), and DB(V,t), which stand for the droplet birth and death rates caused by coalescence and breakage. S(V , t ) = BC(V , t ) − DC(V , t ) + BB(V , t ) − DB(V , t ) (10)
BC(V , t ) =
1 2
∫0
V
(11)
dV ′ DC(V , t ) =
∫0
B B (V , t ) =
∫V
On the other hand, there are several theories proposed to estimate the coalescence efficiency.19,43,44 The most famous and widely used one is the film drainage theory. This indicates that droplet coalescence only occurs for a droplet collision provided that their contact time exceeds the coalescence time required for drainage of the liquid film between them to a critical rupture thickness. With the parallel-film concept, Luo42 proposed his coalescence efficiency model with a constant of order unity. However, Drumm et al.36 and Chen et al.39 found that it overestimates the coalescence rate. Both of them got satisfactory results through reduction of the coalescence rate or enhancement of the breakage rate by a factor of 10. Hence, the modified Luo and Svendsen’s model was thought to be a potential predictive model for PBM simulations, and this scheme was also employed here for further verification. The modified coalescence efficiency is shown as
a(V − V ′, V ′) n(V − V ′, t ) n(V ′, t )
∞
a(V , V ′) n(V , t ) n(V ′, t ) dV ′
(12)
∞
pg (V ′) β(V |V ′) n(V ′, t ) dV ′
DB(V , t ) = g (V ) n(V , t )
(13) (14)
In the coalescence terms, the aggregation kernel a is used to account for the effective collision probability among droplets. It is the product of collision frequency h and collision efficiency λ. On the other hand, the droplet breakage rate is calculated based on two variables. One is the breakage frequency g, which represents the fraction of droplets breaking per unit time. The other is the probability density function β, which describes the size distribution of the daughter droplets resulting from breakage. 2.2.1. Coalescence and Breakage Closures. To solve the population balance equations, suitable coalescence and breakage closures are necessary. As mentioned above, current coalescence and breakage models always have several unknown parameters. The well-known model of Luo and Svendsen41,42 also has one adjustable parameter in the coalescence terms. It was developed based on the theories of isotropy and probability. The breakage frequency was constructed based on the arrival of turbulent eddies, whose kinetic energy exceeds the increase in the surface energy required for breakage. The
⎛ [0.75(1 + ζ 2)(1 + ζ 3)]1/2 ⎞ 1/2⎟ λ(d1 , d 2) = 0.1 exp⎜⎜ − We ⎟ (ρd /ρc + γ )1/2 (1 + ζ )3 ⎝ ⎠ (18)
where γ is the virtual mass coefficient and We is the Weber number. 2.2.2. Quadrature Method of Moments (QMOM). The QMOM was applied to numerically solve the population balance equations. It has the advantages of few tracked scalars and a dynamic calculation of the size bins. In order to apply the QMOM, the kth moment mk is defined, and it is expressed with N weights wi and N abscissas Li through the Gaussian quadrature approximation of order N. 20015
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
∫0
mk (x , t ) =
∞
Article
N
n(L , x , t ) Lk dL ≈
k ∑ wL i i i=1
the two-way-coupled CFD−PBM was established, and there is no requirement for experimental data or correlations. 2.3. Species Transport Model and Axial-Mixing Calculation. The species transport model was also established for the evaluation of axial mixing. The governing equation without mass transfer is shown as follows:
(19)
Then the number density conservation equation can be written as ∂ [mk ] + ∇·[u ⃗ ·mk ] = Sk ∂t Sk =
BkC
BkC =
DkC
− N
1 2
BkB =
−
DkB
(21)
N
∑ wi∑ wj(Li 3 + Lj3)k/3 a(Li , Lj) i=1
N
∑ Likwi∑ waj (Li , Lj) i=1
j=1
N
∞
∑ wi∫ i=1
0
(23)
Lk g (Li) β(L|Li) dL
(24)
N
DkB =
k ∑ wL i i g (Li ) i=1
(25)
The weights wi and abscissas Li can be calculated through the product−difference (PD) algorithm. Starting from the first 2N moments, it is possible to determine the N weights wi and N abscissas Li that minimize the error by deriving a tridiagonal matrix of rank N and by finding its eigenvalues and eigenvectors. Theoretically, the abscissas are equal to the eigenvalues, while the squares of the first component of the eigenvectors are the weights.45 Once the weights and abscissas are known, the source terms due to coalescence and breakage can be calculated. Then the transport equations for the moments can be solved. The droplet-size distribution can be reconstructed by the solved moments. First, the number density function n(L) is expressed as N−1
n(L) = exp( ∑ Ai Li) i=0
(26)
Then the equation for the kth moment is rewritten as mk =
∫0
∞
i=0
(28)
⎛ μ ⎞ Ji ⃗ = ⎜ρDm + t ⎟ ·∇Yi Sct ⎠ ⎝
(29)
∂Yi ∂ 2Y ∂Y = Dea 2i − u i ∂t ∂l ∂l
N−1
Lk exp( ∑ Ai Li) dL
⎯ ∂ ⎯u Y ) = ∇·→ (ραkYi ) + ∇·(ραk→ Ji k i ∂t
where Yi is the mass fraction of species i, J is the mass diffusion, and Dm is the molecular diffusion coefficient. In turbulent flows, the mass diffusion includes two parts: molecular diffusivity and turbulent diffusivity. The turbulent diffusivity is calculated with the turbulent Schmidt number Sct and turbulent viscosity μt. The default Schmidt number Sct of 0.7 was used. The hydrodynamic parameters in the species transport model come from the CFD simulation directly. For a pulsed tracer injection, the species mass fraction was specified as one in a small area between two stators after the flow field reached a steady state. Then the tracer concentration was monitored continuously at a downstream position to record the residence time distribution. Theoretically, the species concentration distribution has a direct effect on the fluid physical properties, which include the density, viscosity, and interfacial tension. Further, variations of these physical properties also influence the CFD and PBM calculations, more or less as shown in their governing equations. Particularly, a change in the interfacial tension has a significant effect on the droplet coalescence and breakage.19,42,46 In this study, the average species concentration in the whole column was kept at a low level of around 0.2%. Hence, its effect on the physical properties and CFD and PBM calculations was neglected. Tamhane et al.47 simulated the residence time distribution for a single-phase annular centrifugal extractor with the same assumption and demonstrated that it is reasonable. The QUICK discretization scheme was used to solve the species transport equation. For calculation of the Peclet number from the simulated residence time distribution, the following traditional axial diffusion model was employed, which is similar to the species transport equation.
(22)
j=1
N
DkC =
+
BkB
(20)
(27)
Pe =
Given N moments, the coefficients Ai can be found by a globally convergent Newton−Raphson method to obtain the droplet-size distribution. Previous simulations indicated that four moments were enough to accurately describe the droplet-size distribution. The discretization scheme of QUICK and convergence criteria of 10−3 were used for the moments. 2.2.3. Coupling of PBM and CFD. The hydrodynamic parameters of the fluid velocity u, turbulent dissipation rate ε, and dispersed-phase volume fraction α calculated by the CFD simulation were used to solve the population balance equations, which means a one-way coupling between the CFD and PBM. For a two-way coupling, the local Sauter mean diameter calculated with the simulated moments in each cell was returned to the CFD model for the drag-force estimation. So,
ul Dea
(30)
(31)
where l is the characteristic length and Dea is the axial diffusion coefficient. With an ideal unsteady-state tracer injection, the axial diffusion model has the following solution. 2 ⎡ Peτ ⎛ ⎛ Peτ 3 ⎞0.5 t ⎞⎟ ⎤ ⎜1 − ⎢ ⎥ exp − Yi = ⎜ ⎟ τ⎠ ⎦ ⎝ 4πt 3 ⎠ ⎣ 4t ⎝
(32)
Hence, the parameters of Peclet number Pe and average residence time τ can be obtained by regressing the residence time distribution with eq 32. For the above models, the solver employs Fluent 14.0, and all simulations were conducted on an Intel Xeon CPU W3520@ 2.67 GHz workstation running on Windows 7 Professional 64 20016
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
Figure 2. Effect of the mesh size on the axial and radial velocities as simulated by CFD.
Figure 3. Simulated velocity contours and vectors of the aqueous phase: (left) CFD; (right) CFD-PBM.
Figure 4. Simulated holdup contours and velocity vectors of the dispersed phase: (left) CFD; (right) CFD-PBM. 20017
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
Figure 5. Comparison of predicted aqueous-phase velocities with experimental PIV data: (left) at half-compartment height; (right) at the stirrer level.
velocities. The reason may be that the CFD model employed the experimental Sauter mean diameter of 2.5 mm, while the CFD−PBM used the droplet diameter returned by the PBM with a certain engineering error presented later. Although simulation of the droplet-size distribution is the main advantage of the CFD−PBM, it has evident effects on the holdup and axial mixing because of the droplet coalescence and forwardmixing effect, as shown later, but few effects on the flow-field structure and fluid velocity magnitude. Hence, the relatively accurate droplet mean diameter in the CFD model is the main reason why a little better prediction results on the fluid velocities. Nevertheless, the CFD−PBM-simulated result without any experiment is also satisfactory. Further comparisons and discussions are necessary. 3.3. Droplet-Diameter Validation. Estimation of the droplet-size distribution is the most important characteristic of the coupled CFD−PBM framework. In extraction columns, prediction of the Sauter mean diameter is much more important than the droplet-size distribution for calculation of the interface area. Hence, the validation presented here mainly focuses on the Sauter mean diameter d32 calculated from the simulated droplet-size distribution. The empirical correlation proposed by Kumar and Hartland48 was used for comparison. The following empirical correlation includes the effects of the column geometry, operating conditions, and physical properties.
bit with 15.9 GB of RAM. To get enough data for validation with empirical correlations, the agitation rate, flow rates, density, and viscosity were changed reasonably.
3. RESULTS AND DISCUSSION 3.1. Mesh Optimization. To determine the suitable mesh size, its effects on the simulated results were investigated first. The whole geometry was meshed with the sizes of 2.0, 1.0, 0.5, and 0.25 mm for the two-phase CFD simulations, which result in the total cells of 3273, 12976, 51132, and 204528, respectively. As shown in Figure 2, the effects of the mesh size on the CFD-simulated axial and radial velocities along line 1 (see Figure 1) are presented. It can be found that both the axial and radial velocities have apparent variations when the mesh size is reduced from 2 to 1 mm. Although the simulated velocities are still changed when the mesh size is decreased further, the variations and differences are small and insignificant. In contrast, the computational time increases significantly when the mesh size is reduced to 0.25 mm. Hence, the mesh size of 0.5 mm was selected considering its reasonable simulation results and moderate computational costs. 3.2. Flow-Field Validation. For the flow-field validation, the volume flow rates are 200 L/h aqueous phase and 100 L/h organic phase, while the agitation rate is 150 rpm. Figure 3 shows the simulated aqueous-phase velocity contours and vectors. The vortex structure between the disks and stators is predicted successfully by the CFD−PBM simulation. Figure 4 presents the dispersed-phase velocity vectors and holdup contours. It is obvious that the dispersed-phase fluid has a zigzag route, which results from both the rotational rotors and buoyancy force. It also shows that the ascending flow of the dispersed phase is an important factor, which causes the aqueous-phase vortices between internals. For a quantified validation, the aqueous-phase fluid velocities along the lines (as shown in Figure 1) were extracted and compared with the referenced PIV data.13 The results in Figure 5 indicate that the simulated aqueous-phase velocities agree well with those of the PIV measurements, which illustrates that the CFD−PBM has the ability to predict the two-phase flow field reasonably. An unexpected phenomenon is that the CFDsimulated result is a little better than the CFD−PBM-simulated
⎛ ND 2ρ ⎞−1.12 ⎛ μc d32 R c ⎜ ⎟⎟ = C3⎜⎜ ⎜ DR ⎝ μc ⎠ ⎝ σρc DR ⎛ D 2ρ g ⎞−0.05⎛ H ⎞0.42 ⎜⎜ R c ⎟⎟ ⎜ ⎟ , ⎝ DR ⎠ ⎝ σ ⎠
⎞−1.38⎛ Δρ ⎞−0.24 ⎟ ⎜⎜ ⎟⎟ ⎟ ⎝ ρc ⎠ ⎠
Re < 50000 (33)
⎛ ND 2ρ ⎞−0.55 ⎛ 0.23N 2D ⎞ d32 R c R ⎟⎟ ⎟ = C4⎜⎜ exp⎜ − DR g ⎠ ⎝ ⎝ μc ⎠ ⎛ μc ⎜ ⎜ σρ D ⎝ c R
⎞−1.30⎛ ρ ⎞0.75⎛ D 2ρ g ⎞−0.30⎛ ⎞0.28 H ⎟ ⎜⎜ d ⎟⎟ ⎜⎜ R c ⎟⎟ ⎜ ⎟ , ⎟ ρ σ D ⎝ ⎝ ⎠ ⎝ ⎠ R⎠ ⎠ c
Re ≥ 50000 20018
(34)
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
where DR is the column diameter, DS is the stator inner diameter, DR is the rotor diameter, H is the compartment height, N is the agitation rate, μ is the fluid viscosity, and Δρ is the density difference. The coefficients C3 and C4 are 0.18 and 7.01 × 10−3, respectively, for the cases without mass transfer. The correlation has the average deviation of 13%. The comparison presented in Figure 6 illustrates that the CFD−PBM model is appropriate to estimate the droplet
Figure 7. Cumulative droplet-size distributions at different agitation rates as predicted by the CFD−PBM.
where φ is the dispersed-phase holdup. The coefficient C5 in the correlation is equal to 1.37. The whole empirical correlation has an average deviation of 22.4%. Figure 8 presents the results. Although the predicted dispersed-phase holdup data locate in the reasonable error Figure 6. Comparison of the CFD-PBM-predicted Sauter mean diameter with the Kumar and Hartland correlation.
diameter. All of the data locate in the reasonable range with an average deviation of 15.7%. Thus, it can be concluded that the employed coalescence and breakage closures are suitable for the predictive simulation of extractors. Because the accurate droplet mean diameter should be specified to the CFD simulation, this predictive CFD−PBM approach is also thought to be an effective way of reducing experimental costs. The agitation rate is one of the most important operational parameters for a RDC. Figure 7 displays the droplet-size distributions at different agitation rates. It can be seen that the droplet diameter decreases with increasing agitation rate. The reason is that the growth of the turbulent dissipation rate results from the rising agitation rate makes droplet breakage dominant gradually. Hence, increasing the agitation rate is always beneficial to obtaining a large interface area. Of course, an excessive agitation also leads to serious axial mixing, even emulsification. 3.4. Dispersed-Phase-Holdup Validation. As an important physical variable, dispersed-phase holdup has direct effects on the performance of extraction columns. The timeaveraged holdup data predicted by the CFD−PBM simulation were compared with the empirical correlation of Kumar and Hartland.49
Figure 8. Comparison of the predicted dispersed-phase holdup data with the Kumar and Hartland correlation.
range of 20%, a slight overestimation clearly exists. This is mainly attributed to the experimental deviation. The empirical correlation was developed based on the experimental measurements with the shutdown procedure. As Brodkorb and Slater expounded,50 this procedure always gives low holdup results because some droplets can be trapped under stators, which is known as the static holdup. Usually, the static holdup is very small. The simulated data in this paper are the total holdup. Hence, the slight overestimation appearing here is acceptable, and the predicted holdup result is considered to be satisfactory. An interesting phenomenon was found according to the simulated results. With increasing agitation rate, the dispersed phase accumulates more and more under the rotors, while
0.26 ⎛ 0.59N 2D ⎞⎛ u 2 ⎞0.42 ⎛ uc ⎞ d R ⎟ ⎜1 + ⎟ ⎟⎜ φ = C5 exp⎜ g ud ⎠ ⎠⎝ D R g ⎠ ⎝ ⎝
⎛ μc ⎜ ⎜ σρ D ⎝ c R
⎞0.12 ⎛ Δρ ⎞−0.56 ⎛ D 2ρ g ⎞0.15⎛ ⎞−1.18 H ⎟ ⎜⎜ ⎜⎜ R c ⎟⎟ ⎜ ⎟ ⎟⎟ ⎟ ⎝ σ ⎠ ⎝ DR ⎠ ⎠ ⎝ ρc ⎠ (35) 20019
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
Figure 9. Time-averaged dispersed-phase holdup contours at different agitation rates, as predicted by the CFD−PBM.
Figure 4 that the CFD model has a low prediction on the holdup under the internals because of the absence of droplet coalescence. Second, the CFD−PBM-simulated dispersedphase vectors in Figure 4 indicate that some small droplets are easily trapped in the vortices, which are also included in the holdup. The CFD model cannot distinguish the motions of different-sized droplets. Therefore, the holdup is underestimated. 3.5. Axial-Mixing Validation. Especially in large-scale extraction columns, the axial mixing has significant effects on the mass-transfer efficiency. For the quantified estimation of aqueous-phase axial mixing, the residence time distribution was simulated by solving the species transport model and the Peclet number was estimated by regressing with the axial diffusion model. Figure 10 shows the residence time distribution curves at different agitation rates. The fitting result agrees well with the simulation, which indicates that this method is available to calculate the Peclet number. To validate the simulated aqueous-phase Peclet number, the next empirical correlation proposed by Strand et al.51 was used.
variation of the dispersed-phase volume under the stators is insignificant. The detailed time-averaged dispersed-phase holdup contours at different agitation rates in Figure 9 expressly show this difference. Drumm et al.14 also found the same observation in an industrial column hydrodynamic experiment. From this perspective, this simulated phenomenon also contributes to validation of the CFD−PBM. The intensified centrifugal force is considered to be the main reason because in the centrifugal field the light phase tends to move to the central region and the heavy phase tends toward the column wall. Along with the obstruction of rotors, more and more dispersed phase accumulated under the rotors because of the rising centrifugal effect. In addition, the increased holdup again illustrates that increasing agitation rate is advantagous to enlarging the interface area. The holdup data simulated by the CFD model are also shown in Figure 8. Apparently, the holdup result is underestimated. Some data even have great deviations. This is mainly attributed to the assumption of the Sauter mean diameter in the CFD model. First, the flow under the internals is usually weak, which means a high droplet coalescence rate. Then the droplets easily accumulate under the internals. So, it can be seen in 20020
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
PBM simulation. The result in Figure 12 shows that the Peclet number is almost constant when the molecular diffusion
Figure 10. CFD−PBM-simulated and axial-diffusion-model-fitted residence time distributions.
⎛ NDR ⎞⎛ DR ⎞2 ⎡⎛ DS ⎞2 ⎛ DR ⎞2 ⎤ Dea ⎟ ⎢⎜ ⎟ ⎥ ⎟⎜ ⎟ − ⎜ = 0.5 + 0.09⎜ ⎝D⎠⎦ ⎝ u ⎠⎝ D ⎠ ⎣⎝ D ⎠ uH
Figure 12. Effect of the molecular diffusion coefficient on the aqueousphase Peclet number.
(36)
The result in Figure 11 reveals that the CFD−PBM approach has a good prediction on the Peclet number. The CFD−PBM
coefficient is varied from the order of 10−10 to 10−5. So, the molecular diffusion coefficient has few effects on the axial mixing. The reason is that the flow in the extractor is turbulent. Hence, the axial mixing is mainly controlled by turbulent diffusion and molecular diffusion is negligible.
4. CONCLUSION In this study, a pilot-scale RDC was simulated by the CFD− PBM. Predictive droplet coalescence and breakage closures without any adjustable parameters were employed. A comprehensive validation was carried out for the model. By comparison with the referenced PIV data, the simulated flow field was proven to be reasonable. The predicted key parameters of droplet diameter, dispersed-phase holdup, and aqueous-phase Peclet number also have good agreement with the empirical correlations. The average deviations are 15.7%, 12.9%, and 15.2%, respectively. These physical quantities are the main parameters that determine the mass-transfer rate. Hence, their reasonable prediction is meaningful for further mass-transfer simulation. In addition, the model successfully predicted the decreased droplet diameter and increased holdup with rising agitation rate, which conform to physical facts. So, it also helps to demonstrate the model’s validity. The CFD simulation has large deviations in the prediction of the holdup and Peclet number. The assumption of the Sauter mean diameter is the main reason. Hence, coupling of the PBM is essential, and the predictive CFD−PBM approach is a highly efficient method for the design of extractors. In the future, the mass-transfer process is suggested to be coupled further to investigate the mass-transfer coefficient, together with the effects of mass transfer on droplet coalescence and breakage.
Figure 11. Comparison of the predicted aqueous-phase Peclet number with the Strand correlation.
approach has an average deviation of around 15%, which is acceptable. By contrast, the average prediction deviation of the CFD model exceeds 40%. The good prediction of the CFD− PBM is primarily ascribed to simulation of the droplet-size distribution. As is known, different-sized droplets have various terminal velocities, which intensify the axial mixing in extraction columns. The CFD model employed the Sauter mean diameter. Hence, it always underestimates the axial mixing and overestimates the Peclet number. The effect of molecular diffusion coefficient Dm on the aqueous-phase axial mixing was also investigated by the CFD−
■
AUTHOR INFORMATION
Corresponding Authors
*Tel: +86-21-64252826. Fax: +86-21-64252826. E-mail: zsun@ ecust.edu.cn. 20021
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
*Tel: +86-21-64252826. Fax: +86-21-64252826. E-mail: jgyu@ ecust.edu.cn.
(3) Pratt, H. Computation of stagewise and differential contactors. In Handbook of solvent extraction; Lo, T. C., Baird, M. H., Hanson, C., Eds.; Wiley: New York, 1983. (4) Neto, A. P.; Mansur, M. B. Transient modeling of zinc extraction with D2EHPA in a Kühni column. Chem. Eng. Res. Des. 2013, 91, 2323. (5) Schultes, M. Research on mass transfer columns: passé? Chem. Eng. Technol. 2013, 36, 1539. (6) Rieger, R.; Weiss, C.; Wigley, G.; Bart, H.-J.; Marr, R. Investigating the process of liquid−liquid extraction by means of computational fluid dynamics. Comput. Chem. Eng. 1996, 20, 1467. (7) Deshmukh, S. S.; Vedantam, S.; Joshi, J. B. Computational flow modeling and visualization in the annular region of annular centrifugal extractor. Ind. Eng. Chem. Res. 2007, 46, 8343. (8) Deshmukh, S. S.; Sathe, M. J.; Joshi, J. B. Residence time distribution and flow patterns in the single-phase annular region of annular centrifugal extractor. Ind. Eng. Chem. Res. 2009, 48, 37. (9) Tamhane, T. V.; Joshi, J. B.; Mudali, U. K.; Natarajan, R.; Patil, R. N. Axial mixing in annular centrifugal extractors. Chem. Eng. J. 2012, 207−208, 462. (10) Fei, W. Y.; Wang, Y. D.; Wan, Y. K. Physical modelling and numerical simulation of velocity fields in rotating disc contactor via CFD simulation and LDV measurement. Chem. Eng. J. 2000, 78, 131. (11) You, X. Y.; Xiao, X. Simulation of the three-dimensional twophase flow in stirred extraction columns by Lagrangian−Eulerian method. Chem. Biochem. Eng. Q. 2005, 19, 1. (12) Yadav, R. L.; Patwardhan, A. W. CFD modeling of sieve and pulsed-sieve plate extraction columns. Chem. Eng. Res. Des. 2009, 87, 25. (13) Drumm, C.; Bart, H.-J. Hydrodynamics in a RDC extractor: single- and two-phase PIV measurements and CFD simulations. Chem. Eng. Technol. 2006, 29, 1297. (14) Drumm, C.; Hlawitschka, M. W.; Bart, H.-J. CFD simulations and particle image velocimetry measurements in an industrial scale rotating disc contactor. AIChE J. 2011, 57, 10. (15) Onink, F.; Drumm, C.; Meindersma, G. W.; Bart, H.-J.; de Haan, A. B. Hydrodynamic behavior analysis of a rotating disc contactor for aromatics extraction with 4-methyl-butyl-pyridinium·BF4 by CFD. Chem. Eng. J. 2010, 160, 511. (16) Steinmetz, T.; Bart, H.-J. Scale-up of mini-plant extractors on energy dissipation-based population balances. Chem. Eng. Technol. 2009, 32, 693. (17) Valentas, K. J.; Amundson, N. R. Breakage and coalescence in dispersed phase systems. Ind. Eng. Chem. Fundam. 1966, 5, 533. (18) Liao, Y.; Lucas, D. A literature review of theoretical models for drop and bubble breakup in turbulent dispersions. Chem. Eng. Sci. 2009, 64, 3389. (19) Liao, Y.; Lucas, D. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem. Eng. Sci. 2010, 65, 2851. (20) Simon, M.; Bart, H.-J. Experimental studies of coalescence in liquid−liquid systems. Chem. Eng. Technol. 2002, 25, 481. (21) Simon, M.; Schmidt, S. A.; Bart, H.-J. The droplet population balance modelestimation of breakage and coalescence. Chem. Eng. Technol. 2003, 26, 745. (22) Vikhansky, A.; Kraft, M.; Simon, M.; Schmidt, S.; Bart, H.-J. Droplets population balance in a rotating disc contactor: an inverse problem approach. AIChE J. 2006, 52, 1441. (23) Kronberger, T.; Ortner, A.; Zulehner, W.; Bart, H.-J. Numerical simulation of extraction columns using a drop population model. Comput. Chem. Eng. 1995, 19, 639. (24) Attarakih, M. M.; Bart, H.-J.; Faqir, N. M. Numerical solution of the spatially distributed population balance equation describing the hydrodynamics of interacting liquid−liquid dispersions. Chem. Eng. Sci. 2004, 59, 2567. (25) Attarakih, M. M.; Bart, H.-J.; Faqir, N. M. Solution of the droplet breakage equation for interacting liquid−liquid dispersions: a conservative discretization approach. Chem. Eng. Sci. 2004, 59, 2547. (26) Attarakih, M. M.; Bart, H.-J.; Faqir, N. M. Numerical solution of the bivariate population balance equation for the interacting
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We acknowledge financial support from National Natural Science Foundation of China under Grant 21206038, a Specialized Research Fund for the Doctoral Program of High Education for New Lectures under Grant 20120074120014, and the Fundamental Research Funds for the Central Universities.
■
NOMENCLATURE cf = coefficient of the interface area C1, C2, C1ε, C3ε = constants CD = drag coefficient Dm = molecular diffusion coefficient, m2/s Dea = axial diffusion coefficient f v = breakup fraction g = gravity, m/s2, or breakage frequency G = generation of the turbulence kinetic energy, m2/s2 k = turbulence kinetic energy, m2/s2 l = characteristic length, m L = droplet size, m N = agitation rate, 1/s p = pressure, N/m2, or number of child droplets t = time, s u = fluid velocity, m/s, or superficial velocity V = droplet volume, m3 Y = species mass fraction
Greek Letters
α = volume fraction β0 = numerical constant γ = virtual mass coefficient ε = turbulence energy dissipation rate, m2/s3 μ = viscosity, Pa·s μt = turbulent viscosity, Pa·s ζ = size ratio, eddy/droplet, or droplet/droplet ρ = density, kg/m3 σ = interfacial tension, N/m, or turbulent Prandtl numbers for k and ε τ = shear stress, N/m2, or mean residence time, s φ = dispersed-phase holdup Dimensionless Numbers
Pe = Peclet number, ul/Dea Re = Reynolds number, ρul/μ or ρN2DR/μ We = Weber number, ρu2l/σ Subscripts
c = continuous phase d = dispersed phase m = mixture i = ith species
■
REFERENCES
(1) Treybal, R. Liquid extraction. Ind. Eng. Chem. 1961, 53, 161. (2) Bart, H.-J. Reactive extraction: principles and apparatus concepts. In Integrated chemical processes: synthesis, operation, analysis, and control; Sundmacher, K., Kienle, A., Morgenstern, A., Eds.; Wiley: New York, 2005. 20022
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023
Industrial & Engineering Chemistry Research
Article
hydrodynamics and mass transfer in liquid−liquid extraction columns. Chem. Eng. Sci. 2006, 61, 113. (27) Schmidt, S. A.; Simon, M.; Attarakih, M. M.; Lagar, L. G.; Bart, H.-J. Droplet population balance modellinghydrodynamics and mass transfer. Chem. Eng. Sci. 2006, 61, 246. (28) Attarakih, M. M.; Bart, H.-J.; Lagar, L. G.; Faqir, N. M. LLECMOD: a window-based program for hydrodynamics simulation of liquid liquid extraction columns. Chem. Eng. Process. 2006, 45, 113. (29) Attarakih, M. M.; Bart, H.-J.; Steinmetz, T.; Dietzen, M.; Faqir, N. M. LLECMOD: a bivariate population balance simulation tool for liquid−liquid extraction columns. Open Chem. Eng. J. 2008, 2, 10. (30) Jildeh, H. B.; Attarakih, M. M.; Bart, H.-J. Droplet coalescence model optimization using a detailed population balance model for RDC extraction column. Chem. Eng. Res. Des. 2013, 91, 1317. (31) Jaradat, M.; Attarakih, M. M.; Bart, H.-J. Population balance modeling of pulsed (Packed and Sieve-Plate) extraction columns: coupled hydrodynamic and mass transfer. Ind. Eng. Chem. Res. 2011, 50, 14121. (32) Bart, H.-J. Reactive extraction in stirred columnsa review. Chem. Eng. Technol. 2003, 26, 723. (33) Bart, H.-J. Extraction columns in hydrometallurgy. Hydrometallurgy 2005, 78, 21. (34) Ghaniyari-Benis, S.; Hedayat, N.; Ziyari, A.; Kazemzadeh, M.; Shafiee, M. Three-dimensional simulation of hydrodynamics in a rotating disc contactor using computational fluid dynamics. Chem. Eng. Technol. 2009, 32, 93. (35) Vikhansky, A.; Kraft, M. Modelling of a RDC using a combined CFD−population balance approach. Chem. Eng. Sci. 2004, 59, 2597. (36) Drumm, C.; Attarakih, M. M.; Bart, H.-J. Coupling of CFD with DPBM for an RDC extractor. Chem. Eng. Sci. 2009, 64, 721. (37) Hlawitschka, M. W.; Chen, F.; Bart, H.-J.; Hamann, B. CFD Simulation of Liquid−Liquid Extraction Columns and Visualization of Eulerian Datasets. Proc. IRTG 1131 2011, 27, 59. (38) Hlawitschka, M. W.; Jaradat, M.; Chen, F.; Attarakih, M. M.; Kuhnert, J.; Bart, H.-J. A CFD−Population Balance Model for the Simulation of Kühni Extraction Column. Comput.-Aided Chem. Eng. 2011, 29, 66. (39) Chen, P.; Sanyal, J.; Duduković, M. P. Numerical simulation of bubble columns flows: effect of different breakup and coalescence closures. Chem. Eng. Sci. 2005, 60, 1085. (40) Chen, P.; Sanyal, J.; Duduković, M. P. CFD modeling of bubble columns flows: implementation of population balance. Chem. Eng. Sci. 2004, 59, 5201. (41) Luo, H.; Svendsen, H. F. Theoretical model for drop and bubble breakup in turbulent dispersions. AIChE J. 1996, 42, 1225. (42) Luo, H. Coalescence, breakup and liquid circulation in bubble column reactors. Ph.D. Thesis, The University of Trondheim, Trondheim, Norway, 1993. (43) Coulaloglou, C. A.; Tavlarides, L. L. Description of interaction processes in agitated liquid−liquid dispersions. Chem. Eng. Sci. 1977, 32, 1289. (44) Klaseboer, E.; Chevaillier, J.-P.; Gourdon, C.; Masbernat, O. Film drainage between colliding drops at constant approach velocity: experiments and modeling. J. Colloid Interface Sci. 2000, 229, 274. (45) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation breakage processes. J. Colloid Interface Sci. 2003, 258, 322. (46) Chevaillier, J.-P.; Klaseboer, E.; Masbernat, O.; Gourdon, C. Effect of mass transfer on the film drainage between colliding drops. J. Colloid Interface Sci. 2006, 299, 472. (47) Tamhane, T. V.; Joshi, J. B.; Patil, R. N. Performance of annular centrifugal extractors: CFD simulation of flow pattern, axial mixing and extraction with chemical reaction. Chem. Eng. Sci. 2014, 110, 134. (48) Kumar, A.; Hartland, S. Prediction of drop size in rotating disc extractors. Can. J. Chem. Eng. 1986, 64, 915. (49) Mao, Z. Rotating disc contactor. In Handbook of solvent extraction; Wang, J. D., Chen, J. Y., Eds.; CIP: Beijing, 2001.
(50) Brodkorb, M. J.; Slater, M. J. Multicomponent and contamination effects on mass transfer in a liquid-liquid extraction rotating disc contactor. Chem. Eng. Res. Des. 2001, 79A, 335. (51) Strand, C. P.; Olney, R. B.; Ackerman, G. H. Fundamental aspects of rotating disc contactor performance. AIChE J. 1962, 8, 252.
20023
dx.doi.org/10.1021/ie503115w | Ind. Eng. Chem. Res. 2014, 53, 20013−20023