CFD-population balance modeling and simulation of coupled hydrodynamics and mass transfer in liquid extraction columns

Accepted Manuscript CFD-Population Balance Modelling and Simulation of Coupled Hydrodynamics and Mass Transfer in Liquid Extraction Columns Menwer Attarakih, Mark Hlawitschka, Mazen Abu-Khader, Samer Al-Zyod, Hans-Jörg Bart PII: DOI: Reference:

S0307-904X(15)00232-2 http://dx.doi.org/10.1016/j.apm.2015.04.006 APM 10525

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

19 October 2013 7 March 2015 1 April 2015

Please cite this article as: M. Attarakih, M. Hlawitschka, M. Abu-Khader, S. Al-Zyod, H-J. Bart, CFD-Population Balance Modelling and Simulation of Coupled Hydrodynamics and Mass Transfer in Liquid Extraction Columns, Appl. Math. Modelling (2015), doi: http://dx.doi.org/10.1016/j.apm.2015.04.006

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Ver 5: 06.03.2015

CFD-Population Balance Modelling and Simulation of Coupled Hydrodynamics and Mass Transfer in Liquid Extraction Columns Menwer Attarakih1,2*, Mark Hlawitschka2,3 , Mazen Abu-Khader4, Samer Al-Zyod1 and HansJörg Bart2,3 1 The University of Jordan, Faculty of Engineering & Technology, Department of Chemical Engineering, 11942 Amman, Jordan 2 Chair of Separation Science and Technology, TU Kaiserslautern 3 Centre of Mathematical and Computational Modelling, TU Kaiserslautern, P.O. Box 3049 - 67653 Kaiserslautern, Germany 4

Al-Balqa Applied University, Faculty of Eng. Tech. Chem. Eng. Dept., POB 15008, 11134 Amman, Jordan *Corresponding author

Menwer Attarakih The University of Jordan, Faculty of Engineering and Technology Chemical Engineering Department, Amman 11942, Jordan. Tel.: 00962 65158587; fax: 00962 64894294. E-mail addresses: [email protected], [email protected]

Ver 5: 06.03.2015

CFD-Population Balance Modelling and Simulation of Coupled Hydrodynamics and Mass Transfer in Liquid Extraction Columns Menwer Attarakih1,2*, Mark Hlawitschka2,3 , Mazen Abu-Khader4, Samer Al-Zyod1 and HansJörg Bart2,3 1 The University of Jordan, Faculty of Engineering & Technology, Department of Chemical Engineering, 11942 Amman, Jordan 2 Chair of Separation Science and Technology, TU Kaiserslautern 3 Centre of Mathematical and Computational Modelling, TU Kaiserslautern, P.O. Box 3049 - 67653 Kaiserslautern, Germany 4

Al-Balqa Applied University, Faculty of Eng. Tech. Chem. Eng. Dept., POB 15008, 11134 Amman, Jordan *Corresponding author

Menwer Attarakih The University of Jordan, Faculty of Engineering and Technology Chemical Engineering Department, Amman 11942, Jordan. Tel.: 00962 65158587; fax: 00962 64894294. E-mail addresses: [email protected], [email protected] Abstract A hierarchical approach for modelling and simulation of coupled hydrodynamics and mass transfer in liquid extraction columns using detailed and reduced bivariate population balance models is presented. The hierarchical concept utilizes a one-dimensional CFD model with detailed bivariate population balances. This population balance model is implemented in the PPBLAB software, which is used to optimize the column hydrodynamics. The optimized droplet model parameters (droplet breakage and coalescence) are then used by a two-dimensional CFD reduced population balance model. As a reduced bivariate population balance model, OPOSPM (One Primary and One Secondary Particle Method) is implemented in the commercial FLUENT software to predict the coupled hydrodynamics and mass transfer of an RDC extraction column with 88 compartments. The simulation results show that the coupled two-dimensional-OPOSPM model produces results that are very close to that of the one-dimensional PPBLAB detailed population balance model. The advantages of PPBLAB are the ease of model setup, implementation and the reduced simulation time (order of minutes), when compared to the computational time (order of weeks) and computational resources using FLUENT software. The advantages of the two-dimensional CFD model is the direct estimation of the turbulent energy dissipation using the k-ε model and the local resolution of continuous phase back mixing. Keywords: Liquid extraction; OPOSPM, RDC; CFD; population balance; PPBLAB; FLUENT .

Ver 5: 06.03.2015

1. Introduction Population balance modelling (PBM) has recently gained an increased attention from the academic and industrial communities [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. This is because the future needs for equipment design and operation require the introduction of more fundamental transport laws for next-generation of particle processing equipment models [5, 11]. This increased the demand on coupling the population balance equation and detailed Computational Fluid Dynamics (CFD) to simulate two-phase flow problems. In the last decade many researches were concerned in coupling CFD and PBM on one hand, and validating experimentally the coupled model prediction on the other hand [e.g. 3, 12 , 13, 14, 15, 16, 17, 18]. Concerning agitated liquid extraction columns, an early application of CFD simulation was performed to predict the nonideal dispersed phase flow in stirred extraction columns by Modes and Bart [19]. Later on, Drumm and Bart [20] and Drum et al. [21] used CFD to simulate single and two-phase flow field in a short segment of a rotating disc contactor (RDC). The predicted velocity fields were validated by particle image velocimetry (PIV) measurements. These authors used the commercial software FLUENT and the Finite Pointset method (FPM), which is essentially a Lagrangian mesh-free CFD software to simulate the hydrodynamics of RDC column. Hlawitschka and Bart [16] used the particle image velocimetry (PIV) to measure the velocity field and hence the energy dissipation using mini-plant Kühni extraction column. The local averaged phase fraction was measured with a laser induced fluorescence (LIF) system. A further development is achieved by a two-way coupling of the PBM and CFD, where detailed flow fields and turbulent energy dissipation are predicted in the presence of the particulate phase [3, 13, 14]. Tiwari et al. [3] coupled the SQMOM [22] to simulate the hydrodynamics of a fivecompartment RDC extraction column. In a pioneering step, Drumm et al. [14] used a one-group PBM, which is a special case of the SQMOM. The authors were able (for the first time) to simulate the hydrodynamics of a full scale pilot plant RDC column with fifty compartments. This one-group reduced population balance model is essentially a one primary and one secondary particle method (OPOSPM). This provides an efficient CFD-PBM model from implementation and computational point of views [14, 17]. Despite these efforts in coupling CFD and PBM, it is still far from reality to simulate full scale liquid-liquid extraction columns using full discrete population balance models. This is especially when the mass transfer process is included. This is due to the increase in the number of transport equations dictated by the

Ver 5: 06.03.2015

discrete PBM, and the long simulation time to get steady state mass transfer solute concentration profiles. In order to exploit the full potential of the PBM approach, one dimensional PBM-CFD models were developed and proved to be efficient and had good scale up and prediction capabilities [4, 5,6, 7, 8]. In this regard, Attarakih et al. [8] presented the PPBLAB software: Particulate Population Balance LABoratory as a new comprehensive windows-based MATLAB environment. PPBLAB uses the bivariate PBM with respect to particle internal properties: Concentration and size. To resolve the PBM, an extended fixed-pivot technique and the QMOM (Quadrature Method Of Moments) were coupled in a one-dimensional PBM-CFD model [4,, 22]. In the new PPBLAB software, recent advances in column design and performance using new population balance models and solvers are utilized. These include the One Primary and One Secondary Particle Model (OPOSPM) and the bivariate fixed-pivot technique [14, 6]. In this article we will present the modelling and simulation of an extraction column of RDC type, which can process chemical systems with low to moderate viscosity and moderate to high surface tension systems [23, 24]. The coupled PBM and CFD models are implemented using the one-dimensional CFD model of PPBLAB [8] and the commercial software FLUENT 13.0. In a hierarchical step, PPBLAB with a detailed bivariate population balances is used to optimize the column hydrodynamics by solving an inverse population balance problem. The optimized droplet model parameters (droplet breakage and coalescence) are then used by the two-dimensional CFD model to predict steady state solute concentration profile along the column height. 2.Multivariate Population Balance Modelling Let f be the distribution function of droplets in a two-phase liquid dispersion, where f(r, t, v, cy)drdvdcy is the number of droplets at time t in a volume dr around r with size v in the neighbourhood of dv and solute concentration cy in the neighbourhood of dcy. To simplify the algebraic manipulations, the droplets are considered to be spherical with volume v, surface area sp and a mean velocity uy (v,cy), which is conditioned on the particle internal properties v and cy. This means for any given droplet of size v and solute concentration cy at a given point (r,t), there is only one mean velocity uy with a zero dispersion around uy. This defines the droplet distribution function as f(r,t,v,cy)δ(u-uy(v,cy)) with the population and momentum balances are given as

Ver 5: 06.03.2015

∂  f    + ∇e ∂t  uy f 

 uy f   + ∇T ⋅  i T   uy uy f 

 ζ
f   S  ⋅   =    0   Sm 

(1)

This equation describes the evolution of the particulate phase and its interaction with the surrounding continuous phase. The first term on the left hand side of Eq.(1) represents the rate of change of number of droplets and its corresponding momentum per unit volume of the dispersion. The second term represents the movement of droplet clouds due to the mean droplet velocity uy and its corresponding momentum through faces of a unit cube which is parallel to the spatial vector r = [x y z]T. The third term is the convection of the particulate phase along particle internal coordinate: ζ
= ( c
y

v
) , where the over dot represents differentiation with respect to T

time. The rate of change of droplet volume and solute concentration take into account droplet growth or contraction and its internal contents due to interphase solute transfer. The source term S describes droplet-droplet interactions, which lead to droplet breakage and coalescence, while Sm takes into account the net force acting on a cloud of particles by the continuous phase (interfacial forces) and body forces (gravity, buoyancy, lift, virtual mass and Basset forces) [25]. Assuming that the liquid droplet is well mixed (internal concentration gradients are destroyed by breakage and coalescence or internal circulation) [4, 26], the rate of change of solute concentration and droplet volume of a given individual droplet are found by balancing solute mas flux at the interface:

 1 − y    c
y    = K s  v  c* − c oy p  1 y) ( y v
       ρ p 

(2)

Where c* is the droplet equilibrium solute concentration and the overall mass transfer coefficient (Koy) is estimated using the two film resistance theory. This is expressed in terms of the individual mass transfer coefficients in both phases and the solute partition coefficient (m = dcy/dcx). The solute mass fraction y is defined as cy/ρp. Since for liquid droplets ρp >>1, the droplet growth rate v
is small when compared to that of c
y .

Ver 5: 06.03.2015

The bivariate source term (S) in Eq.(1) takes into account droplet breakage and coalescence due to birth and death of droplets of internal states d (droplet diameter) and cy. The details of these interactions are given in [4]. 2.1Momentum balances The momentum conservation of the dispersed phase droplets is given as the second component of the vector equation (Eq.(1)) in an Eulerian frame of reference. To avoid the limitations of the mixture model, it is desirable to track particle trajectories in Lagrangian frame of reference [21, 3]. To set up the Lagrangian model, the source term (Sm) is assumed to consist of interfacial interaction force (Drag force) and the gravity and buoyancy as body forces. Other forces such as the virtual mass force, the lift force and the Basset force can be neglected for liquid–liquid interactions as shown by Wang and Mao [25], where they simulated a two-phase liquid dispersion in a stirred tank. By writing the second component of Eq.(1) in Lagrangian form, one gets the following Equation of motion for single droplets:

τp

Duy Dt

+ u y = k s h (α y )u 0 − u x

dr = uy dt

(3) (4)

where ux is the continuous phase velocity. By expressing the single droplet slip velocity (u0) with respect to the continuous phase in terms of the drag coefficient (CD) and Reynolds number (Re), u0 can be expressed as: g (ρ x − ρy )d 2  24   u0 =   C D Re  18 µx Re =

ρx u y − u x d µx

(5)

(6)

The drag coefficient (CD) can be estimated using the Schiller and Nauman [21] correlation, which is implemented in most CFD simulation packages [3, 21]. This correlation can describe the terminal droplet velocity in liquid-liquid dispersions [22]. The function h(αy) is used to account for the hindering effect due to the presence of other droplets in the dispersion. The simplest form of this function is: h(αy) = (1-αy)b where the

Ver 5: 06.03.2015

exponent (b) is function of the droplet Reynolds number and is introduced to compensate for the increase of the drag coefficient [2, 27]. Brennen [27] reported that the exponent n was in the order of magnitude of 2, while Garthe [2] modified the drag coefficient to take into account the effect of droplet swarm on the relative droplet velocity. For the purpose of one-dimensional liquid extraction column simulation, we introduced the function ks, which takes into account the slowing of droplet motion due to neglecting the effect of the internal column geometry on the continuous and dispersed velocity fields. This is a function of droplet diameter and column internal geometry, which has a mean value of 0.5 for and RDC extraction column [2]. In Eq.(3), the droplet relaxation time (τp) is given by: ρyd 2  24  ρy u 0  = τp = 18µx  C D Re  g ∆ρ

(7)

It is clear from the above equation that the relaxation time is only explicit (independent of droplet velocity) in the viscous region where Stokes’s law applies. In [28], it was shown that the Wesselingh and Bollen[29] approximation (the viscous and inertial regions of the particle motion are covered by one equation) of Eq.(7) was in good agreement with droplet terminal velocity as predicted by Eqs.(3) and (4). Moreover, Eq.(3) is reduced to the algebraic slip velocity model at steady state [4]. The particle relaxation time is found to be in the order of 1 second for the chemical system water-toluene, where toluene is the dispersed phase. Therefore, steady state approximation of Eq.(3) is used in PPBLAB. Other single droplet velocity correlations for liquid droplets are also available in PPBLAB based on the algorithm developed by Godfrey and Slater [23]. Using the Lagrangian frame of reference, the momentum balance for the continuous phase is given by [3, 21].

Dux = −αx ∇e p + ∇e ⋅ τ + αx ρx g + FD Dt 1   τ = αx µx  ∇e ux + (∇e ux )T − (∇e ux )I  3  

αx ρx

(8)

Ver 5: 06.03.2015

Where τ is the stress tensor, µx is the dynamic viscosity of the continuous phase, I is the identity matrix and the drag force per unit volume (FD) is given by:

FD =

C 3 αy ρx D uy − ux (uy − ux ) 4 d

(9)

The continuity equation of the continuous phase holdup (αx) is coupled to the first component of Eq.(1) through the source term: ∞ Dαx
(r, t, d )∂d + αx ( ∇e ⋅ ux ) = −∫ vf 0 Dt

(10)

On the other hand, the continuity equation of the dispersed phase holdup is obtained by multiplying the first component of Eq.(1) by v(d) and integrating both sides from 0 to ∞ w.r.t. v and cy and by using the definition of the dispersed phase holdup αy = ∫∫vf(r,t,v,cy)dvdcy one gets:

D αy + αy ( ∇e ⋅ uy ) = Dt

∞

∫0


(r, t, d )∂d vf

(11)

Note that the source term in Eqs.(10) and (11) takes into account total phase growth or contraction due to solute transfer from the continuous to the population of moving droplets in the continuous phase.

Moreover, Eqs.(10) and (11) are constrained by the sum of the phase

fractions: αx + αy = 1. 2.2 Droplet number balance The single droplet number balance n(d) = f(r, t, d) can be obtained by integrating the first component of Eq.(1) with respect to cy from 0 to c *y :

∂n ∂
) = Sv (n) + ∇e ⋅ ( uyn ) + ( vn ∂t ∂v

(12)

The above equation presents the number continuity equation, which describes the evolution of the droplet phase along spatial domain because of convection (second term on the left hand side) and along the droplet size due to particle growth or contraction (third term on the left hand side) and due to the net generation of particles by breakage and coalescence (right hand side).

Ver 5: 06.03.2015

2.3 Chemical species balances The transport equation, which describes the solute concentration in the moving droplets, is obtained by characterizing each droplet by a mean solute concentration, which is conditioned on the particle size (d). This results in expressing the droplet distribution function as f(r,t,d)δ(uuy(d,cy))δ(c-cy(d)), which means that each moving droplet of size d has its own velocity uy(d,cy)) and mean solute concentration cy(d). This is because the concentration variation around its mean was shown to be small in the case of mass transfer in liquid-liquid dispersions [4, 30]. This assumption was validated by Ribeiro [30] who simulated coupled hydrodynamics and mass transfer in a spatially homogeneous continuous vessel and a non-homogeneous single compartment of a Kuhni liquid extraction column respectively. They found negligible differences between the full bivariate model (w.r.t. d and c) and the reduced model, which is expressed by f(r,t,d)δ(u-uy(d,cy)) δ(c-cy(d)). Attarakih et al. [4] solved the bivariate population balance equation, which describes the coupled hydrodynamics and mass transfer in an agitated RDC extraction column with five compartments. The authors used the bivariate QMOM, where the covariance between droplet size and solute concentration was simulated. It was found to be constant and of small values for the case of simultaneous droplet breakage and coalescence. Using this reduced projected density function (f(r,t,d)δ(u-uy(d,cy)) δ(c-cy(d))), one can integrate the first component of Eq.(1) from 0 to cy,max with respect to cy and from 0 to ∞ with respect to uy after multiplying both sides by the solute mass (v(d)cy). The result is a transport equation which is expressed in terms of the conserved quantity q(d) = ∫vcy f(r,t,d) δ(c-cy(d))dcy = vf(r,t,d) cy :

Dq ∂
)= ( vq + q ∇e ⋅ uy + Dt ∂v

∞

∫0

c
y vf (r, t, d, cy )∂cy + Sq

(13)

Where Sq = ∫vcy f(r,t,d) δ(c-cy(d))Sdcy is an integral source term, which takes in into account the effect of mixing on the interphase mass transfer due to droplet breakage and coalescence. The details of this source term is given in [6]. The velocities along particle property space in Eq.(13) were given by Eq.(2). The solute concentration in the continuous phase (cx) is transported by making a solute balance in the same phase at constant density. This is because the solute quantity, which is transferred from the continuous phase to the dispersed phase is relatively small. The solute balance on the continuous phase is given by:

Ver 5: 06.03.2015

∞ cy max

∂ηx − ∇ ⋅ ( ux ηx + Dax ∇ηx ) = −∫ ∂t 0 ηx = αxcx

∫

υc
y f (r, t, d, cy )∂d ∂cy

0

(14)

In the above equation, Dax is the well-known axial dispersion coefficient, which was introduced in one-dimensional CFD models to count for the departure from the ideal plug flow behavior. The deviation is caused by circulatory flow, small eddies, channeling, and velocity profile [23]. In detailed CFD modeling, the axial dispersion is taken into account by resolving the detailed flow field, which is geometry dependent. Using CFD modelling, it is possible to use local information about fluid flow structure to correlate the axial dispersion coefficient using particle tracing methods. In his recent work, Hlawitschka [45] predicted the axial dispersion coefficient of the continuous phase in a Kuhni extraction column using OpenFOAM software. 2.4 Droplet breakage and Coalescence 2.4.1 Droplet breakage Recent application of PPBLAB software to simulate agitated liquid extraction columns showed that existing droplet breakage and coalescence models can be used successfully to simulate the coupled hydrodynamics and mass transfer in a Kuhni

extraction column [8].

In agitated

columns, droplet breakage due to the interaction with agitated column internals and the countercurrent continuous phase is governed mainly by two factors: The first one is due to the shearing effect of the agitator blade (in the case of RDC column), and the latter one is attributed to the impact of the turbulent eddies due the energy dissipation in the continuous phase [2, 23, 25]. The effect of column internals (agitator and stators) on droplet breakage was observed experimentally and found to take place above a certain critical agitator speed which was dependent on column geometry and chemical system physical properties [2, 23, 25]. The droplet breakage frequency in an RDC extraction column is described by the product of breakage probability and the inverse of droplet residence time in a single column compartment [2]. Since at low dispersed phase holdup droplet breakage is dominant, the droplet breakage probability, which is function of modified Weber number, was successfully correlated based on single droplet experiments in small scale devices [2, 23].

Detailed description of the breakage

probability correlations for extraction columns are found in Garthe [2] and Schmidt et al. [26].

Ver 5: 06.03.2015

The daughter droplet distribution and the mean number of daughter droplets were also experimentally correlated [25]. Another way of modelling droplet breakage frequency is through the use of Coulaloglou and Tavlarides [31] model, which is based on droplet-eddy collisions in an isotropic turbulent field. 2.4.2 Droplet coalescence In agitated liquid extraction columns, droplet coalescence is influenced by many factors which include droplet-droplet, droplet-continuous phase interactions, interfacial forces, and the presence of surfactants [23]. The combined effect of the continuous phase turbulent eddies and external flow field could result in successful collision of two or more liquid droplets [31, 32]. Driven by the complexity of this phenomenon, many theories were developed to predict the droplet-droplet coalescence rate and efficiency [23, 31, 32]. Among these theories was the film drainage theory of Shinnar and Church [31, 32]. In this theory, the liquid droplets may encounter successful coalescence event if the contact time exceeds that required for the drainage of the thin continuous film entrapped between these droplets [33]. One of the fundamental and widely accepted models for binary droplet coalescence is due to Coulaloglou and Tavlarides [31]: 4    C2µx ρx ε  (dd ')   ε1/3  2 2/3 2/3 1/2   ω(d,d ', αy ) = C1 (d + d ') (d + d ' )  exp  −     σ2(1 + αy )3  (d + d ')    1 + αy 

(15)

The above model was developed based on the kinetic theory of gases which consists of the product of two terms: The coalescence rate and efficiency. The coalescence rate describes the rate of collision of two liquid droplets due to the turbulent continuous phase eddies in a Kolmogoroff’s local isotropic turbulent field. The coalescence efficiency was developed using the concepts from the film drainage theory which depends on continuous phase physical properties (viscosity µ, density ρ) and the system interfacial tension (σ). The damping of the continuous phase turbulence at high holdup of the dispersed phase (αy) was taken into account by terms such as (1/(1 + αy)) in Eq.(15). Equipped with the two fitting parameters (C1 and C2), the above model is flexible enough to reproduce experimental data in agitated liquid extraction columns [8, 26]. For more details, the reader can consult reference [33] mechanisms and models.

for a comprehensive and critical review of droplet coalescence

Ver 5: 06.03.2015

2.5 Interphase mass transfer In liquid-liquid extraction columns, the mechanism for the rate of interphase mass transfer is complicated by the internal state of droplet, the accumulation of impurities at the droplet interface and its interaction with the surrounding continuous phase. In general, these mechanisms include drop formation at the tip of distributer nozzles, droplet rise or fall in the continuous phase and droplet breakage and coalescence. The relative importance of these mechanisms depends on the details of the equipment internals. For example, mass transfer during droplet formation and coalescence (with short free path length) is dominant in columns with internals, while in spray columns droplet rise or fall is the dominant interphase mass transfer mechanism[23,34]. In [41] it was was discussed that the droplet side mass transfer coefficient depends on the behaviour of individual droplets, which can be stagnant, circulating or oscillating [23, 35, 36]. This may be distinguished based on the droplet’s Reynolds number as suggested by [37], where this can be used to select individual film mass transfer models. In addition to this, the intensity of energy dissipation in stirred extraction columns was reported to increase the mass transfer coefficient below a certain critical agitator speed [23].

It is believed that most of the

enhancement of the mass transfer rate occurs in regions near the rotor edges where droplet motion is governed by the turbulent continuous phase [23]. Following these lines, Young and Korchinsky [38] and Wegener et al.[34] developed a droplet side mass transfer model to predict the mass transfer coefficient by modifying the molecular diffusivity inside spherical droplets. Both models have the disadvantage of using an experimental fitting parameter, which needs experimental data from single droplet or pilot plant experiments. Concerning an RDC extraction column, we found in our previous work [4] that the model of Handlos and Baron [39] and the proposed correlation of Kumar and Hartland [35] were adequate to predict similar mass transfer profiles in a short segment of an RDC extraction column (0.15 m diameter). In contrast to this work, a bivariate Quadrature Method of Moments and an extended fixed-pivot technique were used in [4] to predict the coupled hydrodynamics and mass transfer dynamic performance of a short segment RDC extraction column. On the other hand, the film mass transfer coefficient for the continuous phase is subjected to the same selection procedure and limitations as described above.

Ver 5: 06.03.2015

Based on the aforementioned brief discussion we used in this work the model of Handlos and Baron [39] and that of Kumar and Hartland [35] to predict the film mass transfer coefficient for the continuous phase. For the dispersed phase, the model of Kronig and Brink [40] (for diffusional and internal circulation states) and the Kumar and Hartland correlation [35] (for circulating and oscillating states) were used. With these individual mass transfer coefficients at hand, the overall mass transfer coefficient (Koy) is then estimated based on the two-resistance film theory [23]. The thermodynamic equilibrium solute concentration ( cy* = (dcy / dcx )cx ) can be estimated using simple correlations or even more complicated thermodynamic models

such as the

UNIQUAC model [10]. 3. Numerical resolution of the one-dimensional PBM The mathematical model of the liquid-extraction column, which is given by equations 3 through 14 along with the consecutive equations for droplet breakage and coalescence, is implemented in PPBLAB software. The PPBLAB software offers a bivariate population balance solver, which is based on the extended-fixed pivot technique [4]. The underlying idea behind this solver is to accommodate the source term, which results from solute material balance for droplets in the size range: d ± ∆d. The balance should guarantee the conservation of total number and solute concentration (mass). Since at the discrete droplet level the newly born droplets by breakage and/ or coalescence events could not be represented exactly by any available grid point. Direct numerical discretization of droplet size can only conserve one total property [4, 6]. To overcome this problem, these newly born particles are distributed between the nearest two representative (pivots) contiguous grid points to conserve any two desired properties. For detailed numerical presentation of this population balance solver, the reader may refer to the references [4, 6]. On the other hand, the one-dimensional spatial model is coupled through the convective and source terms and is dominated by convection transport with discontinuous breakage and coalescence events. The time scale of these events is considered too small when compared to that of convection and diffusion transport mechanisms. Therefore, moving front capturing methods such as the upwind and non-oscillatory central difference schemes are required. In this paper, we are restricted to the steady state simulation of agitated extraction columns. Therefore, sharp profiles (dispersed phase holdup and solute concentration) are not expected to appear along the column height [10, 41]. This allows the use of first-order upwind scheme with flux vector

Ver 5: 06.03.2015

splitting to distinguish upward and downward moving droplets (small droplets are entrained against the positive direction of the large ones), which is implemented in PPBLAB as a default steady state one dimensional space solver. The hierarchy of the PPBLAB numerical solution algorithm is shown in Fig.(1). The two main loops for the population balance equation (PBE) solver and the one dimensional space solver are shown within the dashed line loops. The result after internal (droplet diameter) and external (column height) discretization is a system of ODEs as shown in the center circle of Fig.(1). These three systems of ODEs are: One system for dispersed phase holdup (dϕ ϕ /dt), one system for solute concentration in the dispersed phase (dq/dt) and another one equation is for solute concentration in the continuous phase (dX/dt)). These are further discretized with respect to time using explicit scheme which is further stabilized by lagging the nonlinear time dependent terms. The size of each system is Nc×Np, where Nc is the number of actual compartments or numerical cells along the column height and Np is the number of pivots (number of droplet diameter grids). 4.Numerical results and discussion

In this section the PPBLAB software is used to simulate a pilot plant RDC extraction column using the population balance framework. The underlying idea behind this is based on the decomposition of the different hydrodynamics mechanisms (droplet breakage and coalescence) and mass transport processes into isolated or interacted phenomena in small scale devices [2, 5, 25]. Therefore, the RDC column correlations for droplet breakage, coalescence, droplet slowing factor and mass transfer correlations as obtained from single droplet experiments were implemented using PPBLAB. The chemical system is the toluene-acetone-water, where its physical and transport properties at 20 °C, are summarized in Table (1). The pilot plant RDC column operating conditions are given in Table (2 ), which are necessarily the same as those used by Garthe [2]. The pilot plant internal and external geometry are shown in Table (3). In the present work, a hierarchical simulation approach is proposed as shown in Fig.(2). This is realized by fitting the RDC column hydrodynamics to the experimental data using a onedimensional reduced PBM-CFD model as a fast and cheap solver. This is because the coupled detailed PBM to the two-dimensional CFD model is a very expensive task in terms of simulation time and implementation. It is found that the computational time of the two-dimensional CFD model for the column hydrodynamics using a detailed population balance model is in the order

Ver 5: 06.03.2015

of one week. On the other hand, the computational time is reduced to around one day using OPOPM (One Primary and One Secondary Particle Model), which is essentially a one-group reduced population balance model [14]. Fig.(2) shows this hierarchical approach, where the prediction of the mass transfer profiles is the final objective in the simulation cycle. Figure (1): Hierarchy of PPBLAB numerical solution algorithm.

Table (1 ): Physical properties of the EFCE chemical test system at 20 °C [2] for the simulation of an RDC column using PBLAB and FLUENT software.

Table (2): RDC pilot plant operating conditions [2] using PBLAB and FLUENT software. Table (3): RDC internal and external column geometry [2] as used in PPBLAB and FLUENT software. Table (4): User defined numerical parameters and correlations for the simulation of an RDC column using PBLAB and FLUENT software.

4.1 RDC Simulation using PPBLAB

According to the hierarchical simulation approach of Fig.(2), PPBLAB is used to simulate the fast RDC column hydrodynamics in one-dimensional space, which approaches steady state faster than the diffusional and convective mass transfer. This provides a convenient way to optimize the droplet breakage and coalescence model parameters (e.g. C1 and C4, which appear in Eq.(15)). Using the experimental data of Garthe [2] for the dispersed phase holdup and the mean droplet diameter, a weighted least square method is used. This is formed by the weighted sum of residuals based on the difference between the measured and predicted hydrodynamics data [7]:

min eT We , C

s.t. e = (∆α1...∆αL ...∆d1...∆dL )T

W=

(16)

diag(w1α ...w Lα ..w1d ...w Ld )

In the above equation the operator ∆ denotes the difference between the measured and experimental data points along the RDC column height at steady state. Note that different set of weights are used in the matrix W due to the complication of optical droplet size distribution

Ver 5: 06.03.2015

when compared to the holdup measurement [42]. This results in an uncertainty in droplet mean diameter higher than that in the measured holdup. Fig.(3) shows the results of simultaneous fitting of the dispersed phase holdup and the mean droplet diameter profiles along the column height using the data of Tables (1) to (5). As reported by Attarakih et al. [8], the effect of number of pivots on the predicted hydrodynamics profiles is found minimal when using more than 10 pivots. The estimated values of the droplet coalescence parameters (C1 and C2) are shown in Table (4). As mentioned above the mean droplet diameter was given less weight due to relative high uncertainty in the measurements. The objective function given by Eq.(16) is found insensitive to changes in the parameter C2, which is in agreement with the published literature [43]. In previous works [8, 41] our own results showed that fitting of the column hydrodynamics resulted in a very good prediction of the solute concentration profiles along the column height. Figure (2): A hierarchical simulation approach for liquid extraction columns: From 1D experimental hydrodynamics fitting to the 2D prediction of mass transfer profiles.

Figure (3): Fitting the PPBLAB one dimensional CFD model for the RDC extraction column using the experimental hydrodynamic data of Garthe [2] with simulation conditions as given by Tables (1 to 5). The fitting parameters are only of Coulaloglou and Tavlarides [31] model. The upper panel is the dispersed phase holdup and the lower panel is the mean droplet diameter along the column height.

Figure (4): Volume fraction distribution along the column height as predicted by PPBLAB one dimensional CFD model for the RDC extraction column using the experimental hydrodynamic data of Garthe [2] with simulation conditions as given by Tables (1 to 5).

To get an insight into the local volume droplet distribution inside the column, Fig.(4) shows the mean volume fraction distribution as predicted by PPBLAB using the full population balance model. The upper panel of Fig.(4) shows the steady state evolution of the two-dimensional volume fraction distribution at steady state. The inlet mean diameter is increased from 2.5 mm to an approximately constant value of 2.9 mm over a short distance along the column height (≈ 1.5 m) due to droplet-droplet coalescence in the first few compartments (refer to Fig. (3-lower panel)

Ver 5: 06.03.2015

). The lower and upper settling zones are identified by a sudden discontinuity in space along the column height as shown in Fig.(4-upper panel). It is worthwhile to mention here that these local information are only the mean values averaged over the whole column cross-sectional area. The detailed local mean droplet diameter and the dispersed phase volume fraction are depicted in Fig.(6) using the detailed 2D-OPOSPM model. In this figure, the local droplet diameter exceeded the mean droplet diameter below the stator rings where the the droplet coalescence is dominant. Fig.(8) shows the one-dimensional PPBLAB prediction of the solute concentration profiles along the column height as compared to the experimental data of Garthe [2]. Despite the simplicity of this model, which is reflected by the required simulation time (order of 4 minutes for real simulation time of 10000 s running on a lap top of Intel (R) Core (TM) i5-2450M [email protected] GHz and 4GB RAM), the predicted results are in close agreement with the experimental data. Nevertheless, the concentration profile of the dispersed phase shows higher values in the upper part of the column. In spite of this, this shows again the efficiency of the hierarchical concept of Fig.(2). More RDC simulation cases using the one-dimensional CFD model as coupled with OPOSPM were presented in [41].

4.2 RDC Simulation using FLUENT

The pilot plant RDC extraction column, with operating conditions and geometry as given in Tables (2) and (3) respectively, is simulated using a two-dimensional CFD model coupled with OPOSPM as reduced population balance model [14]. The basic equations that are solved using FLUENT 13.0 are those to track the fluid flow, which include the conservation equation of mass, conservation equation of momentum and solute component mass balances. The turbulent energy dissipation (ε), which is required to close Eq.(15) is calculated using the realizable k-ε model. The mixture turbulent model is based on the average density and velocity of the phases in each cell. As numerical solvers, the first order upwind differencing schemes were used to discretize the convection terms initially and QUICK schemes are used for the final solution. PRESTO was applied as the discretization method for pressure. The first order implicit scheme was used for time advancing and the pressure-velocity coupling was done using the SIMPLE algorithm in two-phase flow. The viscous boundary sublayer is accounted by the enhanced wall treatment model, whereas the sublayer is fully resolved by refining the mesh such that y+ (the wall distance unit) is approximately equal to 1.

Ver 5: 06.03.2015

A two-way coupling is realized between the two phases using the mean droplet diameter (d30), which is calculated using the one-group model of OPOSPM [14]. On the other hand, the column hydrodynamics is coupled to the solute transfer through the mixture physical properties. These are: Mixture density, viscosity and the interfacial tension, which is calculated using the correlation of Misek et al. [44]. The coupled CFD-OPOSPM RDC column model is run on a single non optimized server PC with 2 AMD Opteron (4184) processors (each 6 cores) and a total of 48 GB 667MHz RAM. FLUENT 13.0 is running under UBUNTU 2.6.27-14-server version in parallel (10 cores) using the double precision solver. In the present calculations, the inflow and outflow zones are reduced (to reduce the size of the computational grid) and the column is described by a two dimensional rotational symmetrical mesh which contains at the start of the simulation 69363 cells. At the beginning, the computational time step is set to ∆t = 0.05 s, which is a compromise between accuracy and convergence time of the simulation. As shown in Table (5), the mesh is refined gradually during the simulation to avoid the dispersed phase flooding which may result from the accumulation of dispersed phase underneath the stators. The time step is reduced at a simulation time of 5000 s to obtain a better convergence of the simulation. After 5400s, the mesh is continuously refined to fulfill the near wall distance of y+ = 1. The final resulting mesh is shown for the whole column and in detail for a section in Fig. (5).

Table 5: Two-dimensional gradual grid refinement during the simulation of an RDC extraction column using FLUENT 13.0 software. Figure(5): Details of the CFD simulation grid for the simulation of an RDC extraction column using coupled OPOSPM-FLUENT . a- Half plane of the RDC column showing the rotors and stators. b- Enlarged sample of the computational grid. The predicted column hydrodynamics using CFD-OPOSPM model is shown in Fig.(6). In Fig.( 6-a), the disperse phase enters the compartment from below the stator, rises to the stirrer plate where it is pushed to the outside of the column and finally accumulates again under the stator plate. Due the high phase fraction and the reduced agitation intensity, the liquid droplets form a high phase holdup under stators. Therefore, the droplet size depends on the location inside the compartment. In the active height, the droplet size under the stators reaches a maximum value of 4.7 mm. The coalesced droplets tend to break up at the stirrer plate resulting in a droplet size of

Ver 5: 06.03.2015

3.2 mm in the lower half of the compartment and 2.2 mm in the upper part of the compartment. A high accumulation of dispersed phase can also be observed under the stirrer plate which even exceeds the accumulation under the stator plate resulting from the centrifugal force and the difference in density between the two phases. The average simulated dispersed phase holdup is found to be 5 percent higher than that of the experimental value [2]. Figure (6): CFD simulations of the dispersed phase holdup and the mean droplet size at the bottom of the RDC extraction column. a- Dispersed phase fraction (-). b- Mean droplet size (mm). The detailed velocity fields of the continuous and dispersed phases as predicted by the simulated CFD-OPOSPM model are shown in Fig.(7) for a single compartment. The continuous phase velocity forms a double vortex structure above and below the stirrer. The vortex structure of the continuous phase velocity is thereby influenced by the dispersed phase, which results in a fully developed circular structure in the upper part of the compartment. On the other hand, in the lower part of the compartment, the dispersed phase counteracts the continuous phase velocity, which results in a deformed vortex structure. These vertices and the different rise velocities of the dispersed phase droplets, result in what so called back mixing in both phases. This results in the reduction of the driving force available for solute transfer. It is clear from Fig.(7) that the highest velocity is at the tip of the rotor which helps, along with the shear forces, to break up the dispersed phase droplets. This results in a high dispersed phase volume fraction in these regions as can be seen in Fig.(6-a) due to the increased residence time of the small liquid droplets. On the other hand, the small velocity magnitudes of the continuous phase velocity (Fig.(7-a)) below the stator rings results in the accumulation of the small dispersed phase droplets. This enhances the droplet-droplet coalescence which in turn results in an increased droplet velocity away from the stagnant corners (dead zones) as can be seen from Fig.(7-b). The solute concentration profiles along the column height are shown in Fig. (8). It is clear that the CFD-OPOSPM model fits well the experimental data of Garthe [2]. This good agreement can be attributed to the detailed resolution of the velocity fields and the use of the turbulent energy dissipation in both droplet breakage and coalescence models of Coulaloglou and Tavlarides [31]. This eliminates the need for lumped correlations to predict the energy dissipation, which is the case in the PPBLAB one-dimensional CFD model.

Ver 5: 06.03.2015

However, the increased accuracy of the 1D PPBLAB software when compared to the 2DOPOSPM model in FLUENT is due to the optimized one-dimensional correlations based on actual experimental data. On the other hand, the 2D-OPOSPM model is still promising in the case of the lack of experimental data. This is the case when new column internal geometry is encountered or when the ranges of the column operating conditions lie outside the range of the experimental correlations. Unfortunately, both 1D and 2D modelling approaches still relies on experimentally correlated models when droplet interactions are to be included (breakage and coalescence).

On the other hand, the slight gain in accuracy (when using 2D-OPOSPM model) is at the expense of implementation complexity and the drastic increase in the simulation time. The total computational simulation time using CFD-OPOSPM with the mesh refinements is around 27 days, when compared to 4-6 minutes using PPBLAB. Nevertheless, a quiet good convergence of the concentration profile could be obtained without further mesh refinement after 5000 s of simulation time, which required a computational time of 4.63 days. Figure (7): CFD simulations of the continuous and dispersed velocity fields in compartment of an RDC extraction column using coupled OPOSPM-FLUENT. a-The continuous phase velocity field. b- Dispersed phase velocity field. The colors of the vectors indicate the velocity in m/s Figure (8): CFD simulations of the solute (acetone) mass fraction profiles in the continuous and dispersed phases of an RDC extraction column using the one dimensional PPBLAB and two-dimensional coupled OPOSPM-FLUENT models. The lower curve is the solute mass fraction in the dispersed phase, while the upper one is that in the continuous phase. 5. Summary and Conclusions

In this work the, the importance of modelling of liquid extraction columns (RDC column) highlighted and discussed using coupled PBM and CFD. A hierarchical approach is introduced both at the population balance and CFD levels. At the population balance level, OPOSPM (as a reduced population balance model) is used successfully to predict the coupled hydrodynamics and mass transfer in a pilot plant RDC column using two-dimensional CFD model. At the CFD level, a detailed bivariate population balance model is used to predict the one-dimensional RDC performance with results that are in close agreement with that using the two-dimensional CFDOPOSPM model. At the experimental level, the model ability to predict the experimental

Ver 5: 06.03.2015

performance is demonstrated to predict the dispersed phase holdup, mean droplet diameter and solute concentration profiles. For the one-dimensional CFD model, the only adjusted and sensitive parameter is found to be the pre-exponential parameter in the droplet coalescence model of Coulaloglou and Tavlarides [31]. However, for the two-dimensional CFD model droplet breakage and coalescence parameters are required. This is because droplet breakage probability based on single droplet experiments is only valid for one-dimensional CFD models. These parameters were optimized using the one-dimensional PPBLAB CFD model, which is extremely fast simulator when compared to the two-dimensional CFD models. The column mass transfer profiles were then predicted independently without any further parameter adjustment. By comparing the two-dimensional (FLUENT 13.0) and one-dimensional (PPBLAB) CFD models, it is found that both models are able to predict the RDC hydrodynamics and mass transfer. The advantage of PPBLAB is twofold: First, ease of setting up and implementation of the model, and Second the small computational time (order of 4-6 minutes) when compared to that of FLUENT 13.0 (order of weeks). However, PPBLAB requires experimental correlations to estimate the lumped turbulent energy dissipation and backmixing of the continuous phase. This information is available using the two-dimensional CFD model where local distribution of droplets velocity of both phases and the k-ε model were used to account for back mixing and turbulent energy dissipation. Unfortunately, both 1D and 2D modelling approaches still relies on experimentally correlated models when droplet interactions are to be included (breakage and coalescence).

Acknowledgments The authors would like to thank the Deutsche Forschungsgemeinschaft (DFG) for the financial support for research summer stays in 2012 and 2013 at the University of Kaiserslautern/ Chair of separation sciences. Nomenclature b c c
C1, C2 d d32, d30

droplet swarm exponent solute concentration, kg/m3 droplet concentration time derivative, kg/s empirical parameters in the droplet coalescence model droplet diameter, m volume to surface and volume number mean droplet diameters respectively, m

Ver 5: 06.03.2015

D f F FD g hc I k, Ko ks n Nc Np p pr q r Re s S t u U v v
W y z Greek symbols α ε ζ
η µ ρ σ ω τ ϕ Subscripts a D e i in p r

axial dispersion coefficient, m2/s Number concentration function, 1/m/m3 flux vector of conserved variables drag force vector, N/m3 acceleration of gravity, m/s2 compartment height, m identity matrix individual and overall mass transfer coefficients respectively, m/s droplet slowing factor due to column internal geometry droplet number concentration function, 1/m3 number of compartments in the RDC column number of pivots pressure, N/m2 droplet breakage probability solute mass distribution in the dispersed phase, kg/m3 spatial coordinate vector Reynolds number particle specific surface area, 1/m population or moment balance source term time, s velocity vector, m/s conservative variables vector, m/s particle volume, m3 Time derivative of droplet volum,m3/s weight matrix solute mass fraction in the light phase column spatial coordinate dispersed phase holdup turbulent energy dissipation, m2/s3 velocity vector along droplet space property, (property unit)/time solute mass distribution in the continuous phase, kg/m3 viscosity, Pa.s density, kg/m3 interfacial tension, N/m droplet coalescence frequency, m3/s relaxation time, s local droplet holdup vector axial dispersion drag coefficient external internal Inlet particle relative

Ver 5: 06.03.2015

x y 0

heavy phase light phase terminal or slip

Superscripts *

equilibrium

References: 1. A. Gerstlauer, C. Gahn, H. Zhou, M. Rauls, M. Schreiber, Application of population balances in the chemical industry—current status and future needs, Chem. Eng. Sci. 61 (2010) 205 - 217. 2. D. Garthe, Fluid Dynamics and Mass Transfer of Single Particles and Swarms of Particles in Extraction Columns, Dissertation, TU Munchen, Germany, 2006. 3. S. Tiwari, C. Drumm, M. M. Attarakih, J. Kuhnert, H.-J. Bart, Coupling of the CFD and the Droplet population balance equation with finite poinset method, in: M. Griebel; M.A. Schweitzer (Eds.), Lecture Notes in Computational Science and Engineering: Meshfree Methods for Partial Differential Equations IV, vol. 65, Springer Verlag, Heidelberg, 2008, pp. 315-334. 4. M. Attarakih , H.-J. Bart, N. Faqir, Numerical Solution of the Bivariate Population Balanced Equation for the Interacting Hydrodynamics and Mass Transfer in LiquidLiquid Extraction Columns, Chem. Eng. Sci. 61 (2006) 113-133. 5. M. Kalem, F. Buchbender, A. Pfennig, Simulation of hydrodynamics in RDC extraction columns using the simulation tool “ReDrop”, Chem. Eng. Res. Des. 89 (2011) 1-9. 6. M. Attarakih, H.-J. Bart, T. Steinmetz, M. Dietzen, N. Faqir, LLECMOD: A bivariate population Balance simulation tool for liquid-liquid extraction columns, Open Chem. Eng. J. 2 (2008)10-34. 7. M. Attarakih, H. B. Jildeh, M. Mickler, H.-J. Bart, The OPOSPM as a nonlinear autocorrelation population balance model for dynamic simulation of liquid extraction columns, Comp. Aided. Chem. Eng. 31 (2012) 1216-1220. 8. M. Attarakih, S. Al-Zyod, M. Abu-Khader, H.-J. Bart, PPBLAB: A New Multivariate Population Balance Environment for Particulate System Modelling and Simulation, Procedia Engineering 42 (2012) 1574-1591. 9. M. Attarakih, Integral formulation of the population balance equation: Application to particulate systems with particle growth, Comp. Chem. Eng. 48 (2013) 1-13.

Ver 5: 06.03.2015

10. M. Attarakih, T. Albaraghthi, M. Abu-Khader, Z. Al-Hamamre, H.-J. Bart, Mathematical Modeling of High- Pressure Oil-Splitting Reactor using a Reduced Population Balance Model, Chem. Eng. Sci. 84 (2013) 276–291. 11. S. Mohanty, Modelling of liquid-liquid extraction column: A review, Rev. Chem. Eng. 16 (2000) 199-248. 12. L. F. L. R. Silva, P. L. C. Lage, Development and implementation of a polydispersed multiphase flow model in OpenFOAM., Comp. Chem. Eng. 35 (2011) 2653-2666. 13. C. Drumm, M. M. Attarakih, H.-J. Bart, Coupling of CFD with DPBM for an RDC extractor, Chem. Eng. Sci. 64 (2009) 721-732. 14. C. Drumm, M. Attarakih, M. W. Hlawitschka, H.-J. Bart, One-Group Reduced Population Balance Model for CFD Simulation of a Pilot-Plant Extraction Column. Ind. Eng. Chem. Res. 49 (2010) 3442–3451. 15. M. W. Hlawitschka, M. Jaradat, F. Chen, M. M. Attarakih, J. Kuhnert, H.-J. Bart (2011), A CFD-Population balance model for the simulation of Kühni extraction column. Comp. Aided Chem. Eng. 29 (2011) 66-70. 16. M. W. Hlawitschka, H. J. Bart, Determination of local velocity, energy dissipation and phase fraction with LIF- and PIV-measurement in a Kühni miniplant extraction column, Chem. Eng. Sci. 69 (2013) 138-145. 17. M.W. Hlawitschka, H.-J. Bart, CFD-Mass transfer simulation of an RDC column, Comp. Aided Chem. Eng. 31 (2013) 920–924. 18. A. Buffo, M. Vanni, M., D. L. Marchisio, Multidimensional population balance model for the simulation of turbulent gas-liquid systems in stirred tank reactors,. Chem. Eng. Sci. 70 (2012) 31-44. 19. G. Modes, H.-J. Bart, CFD simulation of nonideal dispersed phase flow in stirred extraction columns. Chem. Eng. Tech. 24 (2001) 1342-1345. 20. C. Drumm, H.-J. Bart, Hydrodynamics in a RDC extractor: Single and two-phase PIV measurements and CFD simulations, Chem. Eng. Tech. 29 (2006) 1397-1302. 21. C. Drumm, S. Tiwari, J. Kuhnert, H.-J. Bart, Finite Pointset Method for simulation of the liquid-liquid flow field in an extractor. Comp. Chem. Eng., 32 (2008) 2946-2957. 22. M. M. Attarakih, C. Drumm, H.-J. Bart H.-J., Solution of the population balance equation using the sectional quadrature method of moments (SQMOM). Chem. Eng. Sci. 64 (2009) 742-752. 23. J. C. Godfrey, M. J. Slater, Liquid-Liquid Extraction Equipment, John Wiley & Sons, Chichester, 1994.

Ver 5: 06.03.2015

24. V.S. Kislik, Solvent Extraction: Classical and Novel Approaches, Elsevier, Amsterdam, 2011. 25. F. Wang, Z.-S. Mao, Numerical and experimental investigation of liquid-liquid twophase flow in stirred tanks, Ind. Eng.Chem. Res. 44 (2005) 5776-5787. 26. S. A. Schmidt, M. Simon, M. M. Attarakih, G. Lagar, H.-J. Bart, Droplet population balance modelling—hydrodynamics and mass transfer, Chem. Eng. Sci. 61 (2006) 246 256. 27. C. E. Brennen, Fundamentals of Multiphase Flows, Cambridge University Press, New York, 2005. 28. M. M. Attarakih, M. Jaradat, C. Drumm, H.-J. Bart, S. Tiwari, V. K. Sharma, J. Kuhnert, A. Klar, A multivariate population balance model for liquid extraction columns. Comp. Aided Chem. Eng. 26 (2009) 1339-1344. 29. J. A. Wesselingh, A. M. Bollen, Single particles, bubbles and drops: Their velocities and mass transfer coefficients. Trans. Ins.t Chem. Eng.77 (1999) 89-96. 30. L. M. Ribeiro, P. F. R. Regueiras, M. M. L. Guimaraes, C. M. N. Madureira, J. J. C. Cruz-Pinto, The dynamic behavior of liquid-liquid agitated dispersions II. Coupled hydrodynamics and mass transfer, Comp. Chem. Eng. 21 (1997) 543-558. 31. C. A. Coulaloglou, L. L. Tavlarides, Description of interaction processes in agitated liquid– liquid dispersions, Chem. Eng. Sci. 32 (1977) 1389–1397.

32. R. Shinnar, J. M. Church, Predicting particle size in agitated dispersions, Ind. Eng. Chem. 52 (1960) 253–256. 33. Y. Liao, D. Lucas, A literature review on mechanisms and models for the coalescence process of fluid particles, Chem. Eng. Sci. 65 (2010) 2851-2864. 34. M. Wegener, A. R. Paschedag, M. Kraume,Mass transfer enhancement through marangoni instabilities during single drop formation, Int. J. Heat and Mass Transfer 52(2009) 2673–2677. 35. A. Kumar, S. Hartland, Correlations for prediction of mass transfer coefficients in single drop systems and liquid-liquid extraction columns, Chem. Eng. Res. Des.77(1999) 372384. 36. M. Amanabadi, H. Bahmanyar, Z. Zarkeshan, M. A. Mousavian, Prediction of Effective Diffusion Coefficient in Rotating Disc Columns and Application in Design, Chinese J. Chem. Eng. 17 (2009) 366-372.

Ver 5: 06.03.2015

37. S. H. Zhang, S. C. Yu, Y. C. Zhou, Y. F. Su, Model for liquid-liquid extraction column performance - The influence of drop size distribution on extraction efficiency, Canad. J. of Chem. Eng. 63 (1985) 213-226. 38. C. H. Young, W. J. Korchinsky, Modelling drop-side mass transfer in agitated polydispersed liquid-liquid systems. Chem. Eng. Sci. 44 (1989) 2355-2361. 39. A. E. Handlos, T. Baron, Mass and heat transfer from drops in liquid-liquid extraction. AIChE J. 3 (1957) 127-136. 40. R. Kronig, J. Brink,

On the theory of extraction from falling droplets, Appl. Sci. Res. A2

(1950) 142-154.

41. M. Attarakih, M. Abu-Khader, H.-J. Bart, Modelling and dynamic analysis of a rotating disc contactor (RDC) extraction column using One Primary and One Secondary Particle Method (OPOSPM), Chem. Eng. Sci. 91 (2013) 180-196. 42. M. Mickler, H. B. Jildeh, M. Attarakih, H.-J. Bart, Online monitoring, simulation and prediction of multiphase flows, Can. J. Chem. Eng. In press (2013) doi: 10.1002/cjce.21893. 43. V. Alopaeus, J. Koskinen, K. I. Keskinen, J. Majander, Simulation of the population balances for liquid-liquid systems in a nonideal stirred tank: Part 2- parameter fitting and the use of multiblock model for dense dispersions, Chem. Eng. Sci. 57 (2002) 18151825. 44. T. Misek, R. Berger, J. Schröter, Standard test systems for liquid extraction, Second Ed., European Federation of Chemical Engineering, London, 1985. 45. M. Hlawitschka, Computational fluid dynamics aided design liquid-liquid extraction columns. , PhD Thesis, Technical University of Kaiserslautern, Kaiserslautern, 2013.

Ver 5: 06.03.2015

List of Figures captions:

Figure (1): Hierarchy of PPBLAB numerical solution algorithm. Figure (2): A hierarchical simulation approach for liquid extraction columns: From 1D experimental hydrodynamics fitting to the 2D prediction of mass transfer profiles. Figure (3): Fitting the one dimensional CFD model of the RDC extraction column using the hydrodynamic data of Garthe [2] using PPBLAB with simulation conditions as given Tables (2 to 5). The fitting parameters are only of Coulaloglou and Tavlarides [31] model. The upper panel is the dispersed phase holdup and the lower panel is the mean droplet diameter along the column height. Figure (4): Volume fraction distribution along the column height as predicted by PPBLAB one dimensional CFD model for the RDC extraction column using the experimental hydrodynamic data of Garthe [2] with simulation conditions as given by Tables (1 to 5). Figure(5): Details of the CFD simulation grid for the simulation of an RDC extraction column using coupled OPOSPM-FLUENT . a- Half plane of the RDC column showing the rotors and stators. b- Enlarged sample of the computational grid. Figure (6): CFD simulations of the dispersed phase holdup and the mean droplet size at the bottom of the RDC extraction column. a- Dispersed phase fraction (-). b- Mean droplet size (mm). Figure (7): CFD simulations of the continuous and dispersed velocity fields in compartment of an RDC extraction column using coupled OPOSPM-FLUENT. a-The continuous phase velocity field. b- Dispersed phase velocity field. The colors of the vectors indicate the velocity in m/s

Figure (8): CFD simulations of the solute (acetone) mass fraction profiles in the continuous and dispersed phases of an RDC extraction column using the one dimensional PPBLAB and twodimensional coupled OPOSPM-FLUENT models. The lower curve is the solute mass fraction in the dispersed phase, while the upper one is that in the continuous phase.

Ver 5: 06.03.2015

List of Tables:

Table (1): Physical properties of toluene-acetone-water chemical test system at 20 °C [2] for the simulation of an RDC column using PBLAB and FLUENT software. Property Density (Kg/m3) Viscosity (Pa.s) Distribution coefficient (concentration/concentration units) Interfacial tension (N/m) Binary diffusion coefficient (m2/s)

Water 992

Toluene 863

0.71 0.024 (calculated bsed on the correlation of Misek et al. [44] Dacetone-water = 1.152×10-9 Dacetone-toluene = 2.778×10-9

Table (2): RDC pilot plant operating conditions [2] using PBLAB and FLUENT software. Continuous Phase flow Rate(L/h) Dispersed Phase Flow Rate (L/h) Continuous Phase inlet Concentration (Kg/m3) Dispersed Phase inlet Concentration (Kg/m3) Rotation Speed (rpm)

61.3 74.5 60.61 0 400

Table (3): RDC internal and external column geometry [2] as used in PPBLAB and FLUENT software. Column internal geometry Stator Diameter(m) Rotor Diameter (m) Compartment Height (m)

0.050 0.045 0.050 Column external geometry

Column Diameter(m) Upper Settling Zone Diameter (m) Lower Settling Zone Diameter (m) Total Number of Compartment Continuous Phase Inlet Compartment Dispersed Phase Inlet Compartment

0.08 0.220 0.150 88 76 17

Ver 5: 06.03.2015

Table (4): User defined numerical parameters and correlations for the simulation of an RDC extraction column using PBLAB and FLUENT software. Number of pivots Population mean (mm) Population variance (mm) Minimum droplet diameter (mm) Maximum droplet diameter (mm) Mean slowing factor [2] Axial dispersion coefficient (continous phase) Terminal droplet velocity Breakage frequency correlation

Coalescence frequency correlation Optimized coalescence frequency constants C1, C2 Continuous phase mass transfer correlation Dispersed phase mass transfer model correlation

15 2.5 0.5 0.01 6 0.52 Kumar and Hartland [35] Vignes (1965) Schmidt et al. [26] for PPBLAB simulation and Coulalogou and Tavlarides [31] for FLUENT with parameters 0.62 (pre-exponential factor) and 0.078 for exponential factor Coulalogou and Tavlarides [31] 0.06, 1.33×1011m-2 (see Eq.(15)) Handlos and Baron [39] Kumar and Hartland [35]

Table (5): Two-dimensional gradual grid refinement during the simulation of an RDC extraction column using FLUENT 13.0 software. Simulation time (s)

Number of Cells

Time step

0-250

69363

0.05 s

250

79026

0.05 s

5000

79026

0.01 s

5400

120036

0.01 s

5800 s

152412

0.01 s

6000

182010

0.01 s

6100

413787

0.01 s

6500

182010

0.01 s

Total actual simulation time

26.85 days

Fig. 1 

Fig. 2 

Fig. 3  20

 (-)

15

10

Line: PPBLAB simulation  Symbols: Exp. Data (Garthe, 2006) 

5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

3

3.5

4

4.5

d32 (mm) 

3

2 Line: PPBLAB simulation  Symbols: Exp. Data (Garthe, 2006) 

1

0

0

0.5

1

1.5

2 2.5 Column Height (m)

Volume Concentration

Fig. 4  0.06 0.04 0.02 0 6 6

4

4

2

2 0

Droplet Diameter (mm)

0

Column Height (m)

0.06 Volumn Concentration

Inet stage

0.05

Middle of the column Top stage

0.04 0.03 0.02 0.01 0

1

2

3 4 Droplet Diamtere (mm)

5

6

Fig. 5 

(a) 

(b) 

Fig. 6 

(a) 

(b) 

Fig. 7 

(a) 

(b) 

Fig. 8 

Conc entration (wt% )

8

PPBLAB CFD Exp.(Garthe 2006)

6

4

2

0

1

1.5

2

2.5 3 Column Height (m)

3.5

4

4.5

Graphical abstract